后验误差估计

后验误差估计（a posteriori error estimation）利用已经计算出的数值解、方程数据和网格，估计它与真解之间的误差。它与事先给出的渐近阶数不同：后验估计应反映当前这一次计算的精度，并能指出哪些局部区域需要继续加密。一个可用的指标还应说明估计的范数、可靠性条件和可能的过估计。

当前网格上的真实误差

考虑 $-u^{″}=2$、$u(0)=u(1)=0$，真解为 $u=x(1-x)$。用两段等长线性有限元，节点为 $0,1/2,1$，离散解在中点为 $1/4$，所以两段斜率分别为 $1/2,-1/2$。节点值虽与真解相同，单元内部仍有误差。

在任一宽为 $h$、中点为 $m$ 的单元上，${u_h}^{\prime }=1-2m$，真导数 $u^{\prime }=1-2x$，所以误差导数为 $-2(x-m)$。积分得 $${\int }_K|u^{\prime }-{u_h}^{\prime }{|}^2\mathrm{d}x=4{\int }_{-h/2}^{h/2}{s}^2\mathrm{d}s=\frac{h^3}{3}.$$ 两单元总能量误差为 $\Vert u^{\prime }-{u_h}^{\prime }{\Vert }_{L^2}=1/\sqrt{12}\approx 0.288675135$。这次已知真解，能直接复算；一般问题中未知真解，估计器必须只用可得数据构造。

残差不是误差本身

弱形式为：求 $u\in H_0^1(0,1)$，使 $\int u^{\prime }{v}^{\prime }=\int fv$ 对所有 $v\in H_0^1$ 成立。对保形近似 $u_h$ 定义残差泛函 $$R(v)={\int }_0^{1}fv\mathrm{d}x-{\int }_0^1{u_h}^{\prime }{v}^{\prime }\mathrm{d}x.$$ 令 $e=u-u_h$，便有 $R(v)=\int e^{\prime }{v}^{\prime }$。Cauchy–Schwarz 给出 $|R(v)|\le \Vert e^{\prime }\Vert \Vert v^{\prime }\Vert$；取 $v=e$ 又达到等号。因此 $$\Vert R{\Vert }_{(H_0^1,\Vert v^{\prime }\Vert {)}^{\ast }}=\Vert e^{\prime }{\Vert }_{L^2}.$$ 这给出残差与能量误差的精确关系，但求这个对偶范数本身通常又是一个全局问题。常用局部残差估计器会将单元内部的 $f+{u_h}^{″}$ 和界面导数跳跃组合起来，并配以网格权重和可证明的常数。只检查线性代数系统的求解残差接近零，并不能证明离散误差小：上面的有限元方程可被精确解出，空间误差仍非零。

平衡通量给出的可靠上界

在一维设重构通量 $\sigma \in H^1(0,1)$ 满足 $-{\sigma }^{\prime }=f$，定义可计算量 $$\eta =\Vert \sigma -{u_h}^{\prime }{\Vert }_{L^2}.$$ 真通量是 $u^{\prime }$。因 $(\sigma -u^{\prime }{)}^{\prime }=0$，并且误差 $e$ 在端点为零，分部积分得 $${\int }_0^1(\sigma -u^{\prime })e^{\prime }\mathrm{d}x=[(\sigma -u^{\prime })e{]}_0^1-{\int }_0^1(\sigma -u^{\prime }{)}^{\prime }e\mathrm{d}x=0.$$ 将 $\sigma -{u_h}^{\prime }=(\sigma -u^{\prime })+e^{\prime }$ 展开平方，交叉项消失，得到 $${\eta }^2=\Vert e^{\prime }{\Vert }_{L^2}^2+\Vert \sigma -u^{\prime }{\Vert }_{L^2}^2\ge \Vert e^{\prime }{\Vert }_{L^2}^2.$$ 这是一条常数为一的可靠上界，前提是精确平衡、正确边界类型、保形近似及准确积分。多维对应的平衡通量需要合适的 $H(\mathrm{div})$ 空间；非平衡数据还需加入数据振荡等项。这一正交恒等式是 Prager–Synge 型估计的一维形式。Vohralík 的讲义系统讨论了能量对偶范数、通量重构及该等式；讲义取物理通量为负梯度，符号与本篇的正导数约定相反。

f=2 时可取 σ=1−2x。两段离散导数为正负0.5，差的平方积分分别为1/24，总估计为1/√12。

本例可取 $\sigma =1-2x$，恰好等于真导数，所以估计等于真实能量误差。若改取同样满足平衡的 $\sigma =2-2x$，则额外差为常数一，得到 ${\eta }^2=1/12+1=13/12$。上界仍可靠，却明显过大。这一有解变式解释了可靠性与效率的区别：平衡保证不会低估，良好重构还要减少无用的通量误差。

局部指标与自适应

将 ${\eta }^2=\sum _K{\eta }_K^2$，其中 ${\eta }_K^2={\int }_K|\sigma -{u_h}^{\prime }{|}^2$，便能按单元比较误差贡献。指标用于自适应网格加密，但它们之和才控制这里声明的全局范数；一个单元指标不一定等于该处真实误差。对尖角奇性、系数跳跃或边界层，需要与模型相匹配的重构和权重。

后验估计还应区分离散误差、迭代求解误差和数据近似误差。若迭代尚未充分收敛，使用假定 Galerkin 正交的估计公式可能失去所需条件；若近似右端代替原始右端，也应计算相应的数据项。

当前例子的精确积分代码

下列函数接受严格递增的节点、对应有限元节点值及平衡通量常数 $C$，针对 $f=2$ 使用 $\sigma =C-2x$，返回局部平方指标和总指标。它不需要真解，也不修改输入；每个单元积分二次多项式，时间和输出存储为 $O(M)$。为了避免不必要的相消，积分写成中点值与局部方差的和。

Python 3

import math


def equilibrated_estimator(nodes, values, C=1.0):
    x, u = list(map(float, nodes)), list(map(float, values))
    C = float(C)
    if len(x) < 2 or len(x) != len(u) or not all(math.isfinite(v) for v in x+u+[C]):
        raise ValueError("matching finite nodes and values required")
    if any(b <= a for a, b in zip(x, x[1:])) or x[0] != 0 or x[-1] != 1 or u[0] != 0 or u[-1] != 0:
        raise ValueError("increasing [0,1] nodes and zero endpoint values required")
    local = []
    for a, b, ua, ub in zip(x, x[1:], u, u[1:]):
        h = b-a; slope = (ub-ua)/h; midpoint = (a+b)/2
        squared = h*(C-2*midpoint-slope)**2+h**3/3
        if not math.isfinite(squared):
            raise OverflowError("indicator outside floating-point range")
        local.append(squared)
    return local, math.sqrt(sum(local))


print(equilibrated_estimator([0, 0.5, 1], [0, 0.25, 0]))
# ([1/24, 1/24], 0.28867513459481287)
print(equilibrated_estimator([0, 0.5, 1], [0, 0.25, 0], C=2)[1])
# sqrt(13/12) = 1.0408329997330663