自适应网格加密

自适应网格加密（adaptive mesh refinement，AMR）根据当前解的局部误差或目标量贡献，选择部分网格进行细化。它把求解、估计、标记和加密组成一个反馈循环，力求把计算资源放到误差较大的区域。网格变细本身不是验收标准；需要用与方程和目标一致的指标重新检查精度。

只在半个区间需要细网格的例子

考虑 $-u^{″}=f$、$u(0)=u(1)=0$，其中 $f=8$ 于 $0<x<1/2$，其余处为零。要求解和通量在界面连续，得到 $$u(x)=\{\begin{array}{cc}-4x^2+3x,& 0\le x\le 1/2,\\ 1-x,& 1/2\le x\le 1.\end{array}$$ 可逐项核对：两侧在 $x=1/2$ 的值均为 $1/2$，导数均为 $-1$，端点为零，单元内部二阶导数符合右端。右半段本来就是直线，用线性单元已经能精确表示；左半段有曲率，需要更细网格。

初始网格 $\{0,1/2,1\}$ 的线性有限元解节点值为 $(0,1/2,0)$。一维常系数问题在精确积分下，这些节点值与真解一致，但单元内部仍有误差。左单元离散导数为一，真导数为 $3-8x$；右单元两者均为负一。因此局部能量误差平方为 $${\eta }_{left}^2={\int }_0^{1/2}(2-8x{)}^2\mathrm{d}x=\frac{2}{3},{\eta }_{right}^2=0.$$ 本例用已知精确通量展示指标，数值上也可按后验误差估计中的平衡通量方法构造。一般问题不能预先使用未知真解给网格打分。

右半段为直线，无需新增节点；左半段在 x=1/4 增加一个节点后，能量误差从 √(2/3) 降至 √(1/6)。

求解、估计、标记和再求解

每轮先在当前网格上求解离散方程，再构造可计算指标 ${\eta }_K$。标记选择需要细化的单元，加密后重新建立离散空间并求解。旧解插值可以作为新迭代初值，但不能代替新方程的求解。误差指标既要衡量总误差，也要指出误差集中在哪些单元，才能为局部加密提供依据；Vohralík 的讲义讨论了这种误差定位及其在自适应加密中的作用。

常见的 Dörfler 标记选择集合 $ℳ$，使 $$\sum _{K\in ℳ}{\eta }_K^2\ge \theta \sum _{K\in 𝒯}{\eta }_K^2,0<\theta \le 1.$$ 按指标平方从大到小取单元，直到达到阈值，是一种直接实现。它控制的是被标记部分占总指标平方的比例，不是固定取“误差最大的百分之几单元”。本文把参数直接记为平方和比例 $\theta$；有些软件用 ${\theta }^2$，调用时应核对约定。

本例只需标记左半单元。在 $x=1/4$ 加节点后，新网格有三个单元，节点值为 $(0,1/2,1/2,0)$。新的左两单元各贡献 $1/12$，总误差平方为 $1/6$。若把右侧也均匀二分，则需要四个单元才能得到相同能量误差；多出的右侧自由度在本例并未减少误差。

局部细化收益的推导

在常右端 $f=c$ 的单元上，真解为二次函数，线性插值的导数等于真导数在单元中点的值。记宽度 $h$，中点 $m$，导数差为 $-c(x-m)$，所以 $${\eta }_K^2=c^2{\int }_{-h/2}^{h/2}{s}^2\mathrm{d}s=\frac{c^2{h}^3}{12}.$$ 将这个单元对半分，两个子单元指标之和为 $2c^2(h/2{)}^3/12={\eta }_K^2/4$。这完整解释了左单元每次二分的局部收益；若 $c=0$，原指标就为零，二分也不改善该范数的误差。

有解变式从三单元网格继续，取 $\theta =1/2$。两个左单元指标相等，各为 $1/12$；任取其中一个便达到总量的一半。二分其中一个后，总指标平方变为 $1/12+1/48=5/48$，误差为 $\sqrt{5/48}$。这说明 Dörfler 标记可能有并列选择，确定性程序应规定相等时的次序，而不把不同合法网格视为算法错误。

可复算的标记与加密

下面函数对严格递增的一维节点和对应非负指标平方作标记，返回细化网格及被选单元下标。不修改输入。排序需 $O(M\log M)$ 时间，其余为 $O(M)$；输出存储为 $O(M)$。零指标总和时不标记任何单元，避免“全都为零”仍无限加密。

Python 3

import math


def mark_and_bisect(nodes, indicator_squared, theta=0.5):
    x, eta2 = list(map(float, nodes)), list(map(float, indicator_squared))
    theta = float(theta)
    if len(x) < 2 or len(eta2) != len(x)-1:
        raise ValueError("one indicator per interval required")
    if not all(math.isfinite(v) for v in x+eta2+[theta]) or not 0 < theta <= 1:
        raise ValueError("finite data and 0 < theta <= 1 required")
    if any(b <= a for a, b in zip(x, x[1:])) or any(v < 0 for v in eta2):
        raise ValueError("increasing nodes and nonnegative squared indicators required")
    total = math.fsum(eta2)
    if not math.isfinite(total):
        raise OverflowError("indicator sum outside floating-point range")
    if total == 0:
        return x, []
    order = sorted(range(len(eta2)), key=lambda i: (-eta2[i], i))
    marked = []; accumulated = 0.0
    for i in order:
        marked.append(i); accumulated += eta2[i]
        if accumulated >= theta*total:
            break
    selected = set(marked); refined = [x[0]]
    for i, (a, b) in enumerate(zip(x, x[1:])):
        if i in selected:
            midpoint = a+(b-a)/2
            if not a < midpoint < b:
                raise ArithmeticError("interval too small for floating-point bisection")
            refined.append(midpoint)
        refined.append(b)
    return refined, sorted(marked)


print(mark_and_bisect([0, 0.5, 1], [2/3, 0]))
# ([0.0, 0.25, 0.5, 1.0], [0])
print(mark_and_bisect([0, 0.25, 0.5, 1], [1/12, 1/12, 0]))
# ([0.0, 0.125, 0.25, 0.5, 1.0], [0])

这段代码只执行一维标记和二分，不包含有限元装配、重新求解或一般后验估计，不能单独称为完整 PDE 自适应求解器。二维三角网格还要保证相邻单元的相容细化；悬挂节点、单元形状质量和多处理器负载均衡也是实现的一部分。

精度目标与停止条件

能量范数小不一定保证某一点的峰值、边界通量或目标积分已经足够准确。面向目标量的自适应可结合对偶问题，使用与目标误差相关的权重。时间依赖问题还应区分空间误差、时间误差和移动结构，可能需要粗化及守恒转移。指标可靠性、网格正则性和求解精度等条件不足时，仅重复四阶段流程并不能自动得到最优收敛定理。