Skip to content

算法Algorithm

自适应有限元的估计、标记与加密

Adaptive finite element method · AFEM

先复算一维局部加密,再在固定二维Poisson模型中证明估计量缩减、能量正交及加权收缩,并解释误差不变的边界加密例子。

均匀加密把计算资源平均分给整个区域;自适应有限元则先观察当前解在哪里还不够准确,再把新自由度放到这些位置。它的关键不是让网格看起来不均匀,而是让“误差证据—单元选择—重新求解”形成可核查的闭环。

本条先完成一个一维任务:从四单元网格出发,计算每个残差指标,选出需要二分的单元,重解系统,并用已知精确解检查误差变化。最后用相同误差、相同自由度两种比较说明加密位置的作用。所有积分和小型线性系统都按精确算术处理。随后为二维协调三角网格证明可逐步复核的加权收缩定理,区分一次算例的收益、收敛保证与最优复杂度。

形式陈述 ​

估计量和标记规则 ​

沿用一维残差后验估计的零端点 Poisson 模型、连续一次元及指标

ηK2=hK2‖f‖K2+12∑xi∈∂K∩(0,1)min(hi,hi+1)Ji2,η2=∑KηK2.

一次循环由四项操作组成:在当前网格上求解离散方程;计算各 ηK2;选择标记集合 M;二分标记单元并在新空间中重新求解。每个新循环都重新计算指标,不能永久沿用第一轮的排序。

这里采用平方形式的 Dörfler 标记:给定 0<θ≤1,要求

∑K∈MηK2≥θη2.

为使任务完全确定,先按 ηK2 从大到小排序,相等时从左到右;取达到门槛的最短前缀。它具有最小的标记单元数,因为任意同样数量的单元都不会超过最大若干项的总和。若 η=0,直接停止。

这里 θ 表示要覆盖的平方指标总量比例。有些文献写成 η(M)≥ϑη,其中 η(M)2=∑K∈MηK2;两者对应 θ=ϑ2,不能把参数数值直接混用。[1]

先明确停止证书 ​

对精确离散解,已有 E≤η/π,所以可以用 η/π≤ε 作为本模型的充分停止条件。指标下降是计算进展,达到具有可靠性常数的门槛才构成误差保证。实际迭代求解的额外误差需要另行控制,不能用一个未经解释的矩阵残差阈值替代。

二维收缩定理需要的明确模型 ​

以下结论固定有界 Lipschitz 多边形区域上的 −Δu=f、齐次 Dirichlet 边界及同一个 f∈L2(Ω)。每轮 uk 都是精确积分、精确求解的协调一次元解。网格 Tk+1 细化 Tk,没有悬挂节点,所得网格族具有统一形状正则界。因此 Vk⊂Vk+1⊂H01(Ω)。

采用二维残差页已经证明可靠性的指标:

hK=|K|1/2,ηk(v,K)2=hK2‖f‖K2+hK∑e⊂∂K, e 内部‖Je(v)‖e2,ηk2=∑Kηk(uk,K)2.

一条内部边在两个单元中各计一次,全局权重为 hK+hL。仍用平方标记 ηk(uk,Mk)2≥θηk2。二维可用固定单元编号打破排序平局;收缩只需要标记不等式,不需要集合基数最小。

再固定整数 b≥1,要求每个标记单元 K 的每个最终后代 K′⊂K 都满足

|K′|≤2−b|K|,hK′≤ρhK,ρ=2−b/2<1.

这是每条后代路径至少经历 b 次等面积二分,不能解释成“整个父单元里任意做了 b 次切割”。为恢复协调性而额外细化邻居是允许的。比如二维最新顶点二分配合协调补全,可以实现这样的细化模块;证明直接使用上述可检查性质。[1]

选面积尺度有实质作用:等边三角形的一次二分仍可能让两个子单元各保留一条原长边,直径都没有下降。形状正则下直径和面积开方虽等价,却不能由此推断直径逐步严格收缩。

先固定旧解,观察估计量如何缩减 ​

暂取任意旧空间函数 v∈Vk,把它视为细网格上的同一个函数。一个旧三角形内的梯度恒定,故其内部新产生的边上 Je(v)=0;旧边上的跳跃只被分段,没有改变数值。

对标记单元 K,将它的所有细后代相加,体积项满足

∑K′⊂KhK′2‖f‖K′2≤ρ2hK2‖f‖K2.

沿旧边分段的 L2 积分可相加,属于 K 这一侧的权重均至多为 ρhK,所以同侧边项至多变为原来的 ρ。无需假设旧边本身全部被切短。对未标记父单元,后代尺度至多等于原尺度,因此相应项不增。因 ρ2≤ρ,

ηk+1(v)2≤ηk(v)2−(1−ρ)ηk(v,Mk)2.

这个不等式比较的是同一个函数在两张网格上的指标;重新求解后函数改变,还要补上下一步。

新旧离散解的变化要付多少代价 ​

在同一张细网格上比较 v+,v,令 d=v+−v,gK=∇d|K。因为二者均为一次元,体积强残差仍都是 f,只需估计跳跃差。每条内部边有

|Je(d)|2≤2(|gK|2+|gL|2).

乘边长及双侧权重并求和,得到

∑e 内部(hK+hL)‖Je(d)‖e2≤Cs‖∇d‖2,

其中可取任一对网格族统一的正上界

Cs≥maxK2|K|∑e⊂∂K, e 内部(hK+hL(e))|e|.

例如把右侧上界再与 1 取最大即可保证 Cs>0。形状正则使共享边两侧尺度可比、周长至多为常数乘 hK,故存在不随加密轮数增长的 Cs。

把各加权体积残差和跳跃组成一个乘积 L2 空间向量,它的范数就是估计量。由三角不等式及 Young 不等式,对任意 δ>0,

ηk+1(v+)2≤(1+δ)ηk+1(v)2+(1+δ−1)Cs‖∇(v+−v)‖2.

取 v=uk,v+=uk+1,结合固定旧解的缩减及标记,令

s=(1−ρ)θ∈(0,1),Dk=‖∇(uk+1−uk)‖,

就有

ηk+12≤(1+δ)(1−s)ηk2+(1+δ−1)CsDk2.

这也解释了为什么不能直接声称重新求解后的估计量每一步都按 1−s 缩减。

正交恒等式把扰动项抵消 ​

记 Ek=‖∇(u−uk)‖。嵌套性使 uk+1−uk∈Vk+1,细空间 Galerkin 正交给出

a(u−uk+1,uk+1−uk)=0,Ek+12=Ek2−Dk2.

这是展开 u−uk=(u−uk+1)+(uk+1−uk) 得到的勾股恒等式。现在明确选择

δ=s2(1−s),γ=s(2−s)Cs,Qk=Ek2+γηk2.

于是 (1+δ)(1−s)=1−s/2 且 γ(1+δ−1)Cs=1。将估计量不等式乘 γ 后与能量恒等式相加,Dk2 项恰好抵消:

Qk+1≤Qk−γs2ηk2.

可靠性 Ek2≤Crelηk2 进一步给出 Qk≤(Crel+γ)ηk2,从而

 Qk+1≤qQk,q=1−s22((2−s)CsCrel+s)∈(0,1). 

分母严格大于 2s>s2,所以参数确实合法。递推得到 Qk≤qkQ0,进而 Ek→0 且 ηk→0。例如 b=2,θ=1/2 时,

ρ=12,s=14,δ=16,γ=17Cs,q=1−18(7CsCrel+1).

这些常数是一个充分选择,不是最佳收缩率;没有另行认证具体网格族的 Crel 时,不能把它任意填成 1。证明依赖固定方程、对称能量、嵌套空间、精确求解、可靠性和规定的细化。它保留完整 f,所以不需另按数据振荡标记;这不删除局部效率定理中的振荡项。若每轮先投影荷载、允许非协调空间或未收敛的线性解,上述等式和指标都需重新分析。[1]

直觉

局部指标把连续方程尚未满足的部分分配到单元上;标记规则再把这些分散证据转成有限的加密选择。覆盖一半平方指标总量,并不要求标记一半单元:如果误差证据集中在少量位置,只增加少量节点就可能显著改善逼近。重新求解后,斜率和跳跃都会改变,因此下一轮应依据新指标重新作决定。

一维二分只需在被选区间中加入中点,相邻区间仍然形成协调划分。高维加密还涉及邻居配合、悬挂节点与形状控制,这些操作由上面的二维细化假设明确约束,不能从一维循环自动推得。

例子与边界

四单元网格上的求解与估计 ​

局部荷载让误差集中在左端 ​

取

f(x)={4,0<x<1/4,0,1/4<x<1,u(0)=u(1)=0.

精确解为

u(x)={7x/8−2x2,0≤x≤1/4,(1−x)/8,1/4≤x≤1.

两段在 x=1/4 的函数值都是 3/32,导数都是 −1/8,因此解及通量连续。左端是抛物线,右侧已经是一次函数。这意味着一次元需要改善的弯曲只存在于左端;右侧即使单元很长也能精确表达这一段解。

初始节点为 0,1/4,1/2,3/4,1。荷载断点与节点对齐,因此每个单元的数据振荡为零,局部真误差可以使用 EK2=cK2hK3/12 精确核验。

组装、求解并计算跳跃 ​

消去边界值后,系统为

A=4(2−10−12−10−12),b=(1/200),U=(3/321/161/32).

第一个荷载分量来自 4∫01/4ϕ1dx=1/2。代入矩阵可逐行验证 AU=b。从节点值计算的四个斜率和三个跳跃是

s=(3/8,−1/8,−1/8,−1/8),J=(−1/2,0,0).

第一个单元的体积项为

h12‖f‖K12=(1/4)2⋅16⋅(1/4)=1/4.

唯一非零跳跃的加权平方是 (1/4)(1/2)2=1/16,左右各收到 1/32。所以

(ηK12,…,ηK42)=(9/32,1/32,0,0),η2=5/16.

真实误差只在第一单元:

(EK12,…,EK42)=(1/48,0,0,0),E=1/48≈0.144338.

第二单元的真误差为零,但因共享左端跳跃而有正指标。后验指标为选择提供证据,并不逐单元复制真实误差。

标记、二分与重新求解 ​

第一轮只需选一个单元 ​

设 θ=1/2,需覆盖的平方指标总量是 5/32。最大项 9/32 单独就达到要求,因此 M={K1}。二分后节点变为

0,1/8,1/4,1/2,3/4,1.

新系统与解分别是

A+=(16−800−812−400−48−400−48),b+=(1/21/400),U+=(5/643/321/161/32).

这次有四个内部自由度。新增节点不是简单取旧折线的中点值;旧插值在那里取 3/64,重新求解后得到 5/64。空间增加之后,求解才能利用新增方向修正原来的误差。

新指标与新真误差 ​

新斜率与跳跃为

s+=(5/8,1/8,−1/8,−1/8,−1/8),J+=(−1/2,−1/4,0,0).

前两个单元的体积项各为 16(1/8)3=1/32。两个非零节点项分别为 1/32 和 1/128,仍各半分给左右单元,因此

(ηK2)+=(364,13256,1256,0,0),η+2=13128.

两个带荷载的单元各有 EK2=16(1/8)3/12=1/384,总计

E+2=1/192,E+=1/192≈0.072169.

能量误差恰好减半,而 η+/η=13/40≈0.5701。估计量与真误差并没有同比下降,可靠性和效率也没有要求它们同比。

如果再运行一轮,同样的标记门槛是 η+2/2=13/256。现在最大项是第二单元的 13/256,恰好单独达到门槛,因此下一步应二分 (1/8,1/4),不是机械地继续二分最左单元。二分后总误差平方为 1/384+1/1536=5/1536;这由单元误差的四分之一规律得到。

加密位置决定新增自由度是否能改善本例的误差。蓝色粗段标出荷载所在区域,新增中点以空心圆表示;E2 是能量误差平方,DOF 是内部自由度数。左端二分与均匀二分具有相同误差,但分别只需 4 与 7 个内部自由度。

相同误差与相同预算的比较 ​

将初始全部四个单元都二分,会得到八个单元。其前两个单元与局部加密网格相同;剩余区间的精确解本来就是直线,多加节点不会降低已经为零的误差。因此均匀二分也有 E2=1/192,并且本例中全局估计量也恰好相同。

反过来,若只二分最后一个单元 (3/4,1),同样使用五个单元、四个内部自由度,却没有改善左端的逼近。新增节点值为 u(7/8)=1/64,其相邻斜率仍为 −1/8,所有非零残差贡献仍在原位置。

网格选择 单元数 内部自由度 η2 E2
初始四等分 4 3 5/16 1/48
只二分左端单元 5 4 13/128 1/192
所有单元均匀二分 8 7 13/128 1/192
只二分右端单元 5 4 5/16 1/48

这张表的结论是具体的:对该荷载,第一次指标标记用较少自由度取得了均匀加密的同一误差;相同预算放错位置则没有改善。它不证明任意荷载、任意标记和任意网格族上的最优复杂度。一般自适应收敛理论还需要分析估计量缩减、网格性质、振荡及相邻离散解之间的关系。[1][2]

代数求解误差不能混进加密判断 ​

能量正交给出的精确分解 ​

设 U 是某张网格的精确离散系数,U^ 是近似求解结果,r=b−AU^。令 d=U−U^=A−1r,则对应有限元函数之间的能量距离是

Ealg2=dTAd=rTA−1r.

由于 uh−u^h∈Vh,Galerkin 正交消去交叉项,得到

‖u′−u^h′‖22=‖u′−uh′‖22+‖uh′−u^h′‖22=Edisc2+Ealg2.

因此两种工作解决不同部分:改善线性求解消除 Ealg,改变网格才可能改善 Edisc。用 ‖r‖2 代替 rTA−1r 会丢掉矩阵对残差方向和尺度的作用。

在初始网格上核验一次 ​

故意取 U^=U+(1/32,0,0)T,则

r=(−1/4,1/8,0)T,Ealg2=(1/32)2A11=1/128.

此时总能量误差平方是 1/48+1/128=11/384。把系统精确解完后,r 和代数误差都变为零,但离散误差仍为 1/48。这个计算也解释了为什么本条先明确“精确求解”,再使用残差指标的证明;对于近似离散解,要额外估计代数误差,不能直接沿用节点精确性。

可复算任务的完成标准 ​

可以独立实现这四张网格的组装与求解。每个单元累加 hK−1(1−1−11),常荷载 cK 的局部右端为 (cKhK/2)(1,1)T;消去两端自由度后求解。用所得节点值计算斜率、跳跃、局部指标,再与表中分数比较。

真误差必须在单元内部积分 ∫K(u′−sK)2,不能只比较节点值,因为本模型在节点处恰好精确。完成后应能解释四件事:为什么第一轮只标记左端;为什么下一轮改标第二单元;为什么均匀加密浪费了本例右侧的新增自由度;为什么矩阵残差为零仍有正的 PDE 误差。这些解释与分数计算共同构成一次完整的自适应任务。

二维边界加密:误差不变,估计量下降 ​

沿用二维残差算例的单位正方形和 f=1,初始网格由中心连向四角。唯一内部自由度为 u0(c)=1/12,且

η02=14+29.

现在标记全部四个单元,把每个三角形从中心连到其外边界边的中点。所得八个三角形协调、形状正则,每个旧单元的两个后代面积都减半,满足 b=1 的条件。

新增节点全部属于 Dirichlet 边界,内部自由度仍只有中心。在每个旧三角形上,旧中心帽函数沿外边界恒为零,因而也在新增边界中点取零;它限制到两个子单元后,正是细网格的中心帽函数。因此 V1=V0,并非仅维数碰巧相同,故

u1=u0,D0=0,E1=E0>0.

严格正误差是因为 u0 在各单元内的 Laplacian 为零,不能满足内部荷载为 1 的连续方程。新边两侧梯度相同,跳跃为零;四条旧内部边的跳跃及边长不变。新单元面积为 1/8,面积尺度为 1/(22),所以

η12=8(18)(18)+4(12)(1182)=18+19=1772.

因此

η12η02=1718+82≈0.579933<2−1/2.

全部标记满足任意 0<θ≤1 的门槛。这个例子并不反驳收缩定理:Q 中的能量部分保持不变,估计量部分严格下降,二者的加权和仍可收缩。它反驳的是不加条件的“每轮能量误差必须严格下降”。反复加密的长期结论还要求整个网格序列保持定理中的形状与细化条件。

推论与应用

一次循环实际需要多少计算 ​

设当前网格有 N 个单元。在各单元所需的荷载积分与 ‖f‖K2 已提供、每次标量运算与比较按常数成本计的模型下,组装、计算斜率与全部局部指标、一维二分及节点重编号各需 O(N) 操作。消去端点后,刚度矩阵是对称正定三对角矩阵,三对角 LDLT 分解与回代需 O(N) 标量运算,存储需 O(N);二分后单元数至多为 2N,重解仍具有同一线性量级。本条指定的完整排序需 O(Nlog⁡N) 比较,之后扫描前缀需 O(N),所以这一实现的单轮成本由排序主导,辅助存储可保持 O(N)。

这些是算术操作计数,不把精确有理数的位数增长算作常数。对任意 f,取得精确单元积分也不是免费步骤;若改用数值积分,其计算量及引入的误差须另外计算。经过多轮、各轮网格大小为 N0,…,NL 时,应累加各轮成本,而不能仅用最后一张网格的大小代替总工作量。前面的同误差比较只核验了这个算例的自由度收益,没有证明达到给定误差所需的总计算量最优。

二维定理给出按迭代轮数的几何收缩,不等于按最终自由度或总工作量达到最优速率。证明这种更强结论还需控制协调补全成本、标记集合大小及可达到的逼近类;本条的终点停在已证明的加权收缩。

参考资料

[1] Ricardo H. Nochetto, Kunibert G. Siebert and Andreas Veeser, Theory of Adaptive Finite Element Methods: An Introduction, 2009,作者讲义。第 93 页的 Mark 条件使用估计量范数形式的参数;§§8.1–8.3,第109–116页,引理18–21和定理17给出正交、可靠性、估计量缩减及加权收缩;第32页定义面积/体积开方尺度,第43页式(54)记录二分代数。本文展开单位扩散、一次元的证明,平方标记参数已明确换算,显式参数选择和一维/二维账本由正文独立计算。

[2] Willy Dörfler, “A Convergent Adaptive Algorithm for Poisson’s Equation,” SIAM Journal on Numerical Analysis 33(3), 1106–1124, 1996,论文页面。这是 Dörfler 标记的原始工作;本条不把该论文特定算法的收敛结论直接赋予上述一维实现。

关系图谱7 个相邻概念 · 2 类关系

拖动节点调整位置。

显示关系

显示:依赖

  1. 前置三跳
  2. 前置二跳
  3. 前置一跳
  4. 当前条目
  5. 后续一跳
  6. 后续二跳
  7. 后续三跳
文字版关系按与当前条目的最短距离分组
类型化关系

使用的工具