4. 一维模型问题的先验估计实践

先验误差估计给出的是“误差上界随网格尺度 h 或多项式阶次 p 如何衰减”的规律，而数值验证要回答的是：实际算出来的误差是否真的按这个规律下降。本节以最简单的一维 Poisson 方程为例，把估计式、实现代码和结果判读串成一条可复现的验证链。

取模型问题：−u''(x) = f(x)，x ∈ (0,1)，边界条件 u(0)=u(1)=0。选择真解 u(x)=sin(πx)，则 f(x)=π²sin(πx)。这样做的目的是让离散误差可以精确计算，从而直接检验收敛阶，而不必依赖后验估计器。约定单元长度为 h，采用连续分片多项式空间 V_h^p（p 为多项式次数），在均匀网格上做 Galerkin 有限元。

先验估计给出的预期

对上述问题，能量范数（H¹ 半范数）误差满足 |u−u_h|_{1} ≤ C h^{p} |u|_{p+1}，L² 范数误差满足 ‖u−u_h‖_{0} ≤ C h^{p+1} |u|_{p+1}。也就是说，固定 p 加密网格（h 收敛）时，能量范数误差的收敛阶为 p，L² 误差的收敛阶为 p+1；固定网格放大 p（p 收敛）时，误差随 p 呈指数型衰减（对解析解而言）。数值验证的任务就是分别测出这两条斜率。

h 收敛的代码实现

下面用 Python 实现 p=1（线性元）与 p=2（二次元）两种情形。为便于教学，采用均匀网格、单元刚度矩阵解析积分，并用 5 点 Gauss 积分计算误差。代码只依赖 NumPy，可直接运行。

import numpy as np

def fem_1d(p, n_elem):
    # 均匀网格，节点数 = p*n_elem + 1
    n_node = p * n_elem + 1
    x = np.linspace(0.0, 1.0, n_node)
    K = np.zeros((n_node, n_node))
    F = np.zeros(n_node)
    # 单元局部形函数在参考单元 [-1,1] 上的 Gauss 积分
    gp, gw = np.polynomial.legendre.leggauss(p + 2)
    for e in range(n_elem):
        idx = np.arange(e * p, e * p + p + 1)
        xe = x[idx]
        he = xe[-1] - xe[0]
        Ke = np.zeros((p + 1, p + 1))
        Fe = np.zeros(p + 1)
        for xi, w in zip(gp, gw):
            # Lagrange 形函数及其导数（参考单元）
            N, dN = lagrange_basis(p, xi)
            J = he / 2.0
            Ke += (dN @ dN.T) / J * w
            xq = 0.5 * (1 - xi) * xe[0] + 0.5 * (1 + xi) * xe[-1]
            Fe += N * (np.pi**2 * np.sin(np.pi * xq)) * J * w
        K[np.ix_(idx, idx)] += Ke
        F[idx] += Fe
    # 齐次 Dirichlet：去掉首末自由度
    free = np.arange(1, n_node - 1)
    U = np.zeros(n_node)
    U[free] = np.linalg.solve(K[np.ix_(free, free)], F[free])
    return x, U

def lagrange_basis(p, xi):
    # 返回参考单元上 p+1 个形函数值及对 xi 的导数
    nodes = np.linspace(-1, 1, p + 1)
    N = np.ones(p + 1); dN = np.zeros(p + 1)
    for i in range(p + 1):
        for j in range(p + 1):
            if j != i:
                N[i] *= (xi - nodes[j]) / (nodes[i] - nodes[j])
        # 导数用乘积法则逐项求
        for k in range(p + 1):
            if k == i: continue
            term = 1.0 / (nodes[i] - nodes[k])
            for j in range(p + 1):
                if j != i and j != k:
                    term *= (xi - nodes[j]) / (nodes[i] - nodes[j])
            dN[i] += term
    return N, dN

误差计算时，把每个单元再细分为若干采样点，用形函数重构 u_h，与真解比较，分别累加 L² 与 H¹ 半范数误差的平方。对 h 序列（如 n_elem = 4, 8, 16, 32, 64）计算误差后，用相邻两点的对数比估计收敛阶：

rate = np.log(err[i] / err[i+1]) / np.log(h[i] / h[i+1])

若实现正确，p=1 时 H¹ 误差阶应接近 1、L² 误差阶接近 2；p=2 时分别接近 2 和 3。若测得的阶次明显偏低，常见原因包括：边界自由度未正确消去、单元刚度矩阵积分阶次不足、误差采样点过稀导致积分不准，或网格编号与形函数节点顺序不一致。

p 收敛的验证方式

p 收敛固定网格（例如 n_elem=4），只提升多项式次数 p=1,2,3,4,5。由于真解是解析函数，误差应随 p 快速下降，在半对数坐标下近似直线（指数衰减）。实现上只需复用同一套代码，把 p 作为参数传入；注意高阶时 Gauss 积分点数要足够（建议 p+2 以上），否则刚度矩阵积分误差会污染收敛阶。若 p 提升后误差不再下降甚至反弹，通常说明积分不足或条件数恶化，而不是理论失效。

结果解读与判据

验证是否通过，建议同时看三条：其一，实测收敛阶与理论阶的偏差在合理范围（一般 0.1～0.2 以内可接受，受网格较粗和常数影响）；其二，误差随 h 单调下降，没有异常跳变；其三，p 收敛曲线在半对数坐标下呈近似直线。需要强调的是，先验估计中的常数 C 与具体问题有关，因此绝对误差大小不能直接与理论界比较，能比较的是斜率。另外，能量范数误差对奇异性更敏感，若把真解换成 x^α 型奇异函数，h 收敛阶会退化到 min(p, α) 量级，这正是先验估计中正则性假设的体现。

通过本节这套“真解—离散—误差—斜率”的闭环，先验估计不再只是纸面上的不等式，而成为可以逐项核对、可定位实现错误的工程判据。后续引入后验估计器时，这套验证框架同样可以作为检验估计器有效性的基准。