有限元方法

网格、形状函数与自由度 ​

把 $[0,1]$ 分成节点 $0=x_0<x_1<\cdots <x_N=1$。最简单的选择是：每个小区间上使用一次多项式，相邻区间在公共节点处连续。零边界条件下，只有内部节点值需要求解。

对内部节点 $x_i$ 定义“帐篷函数” ${\varphi }_i$：它在 $x_i$ 取 1，在其他节点取 0，在每个相邻区间上线性变化，离开这两个区间便为零。于是

这里 $U_i$ 是自由度，${\varphi }_i$ 是与自由度对应的形状函数。节点值是最直观的一类自由度，但并不是唯一的一类；高阶元可以使用边上的节点值，混合元也可能使用边或面上的通量积分。

更一般地，一个有限元由三元组 $(K,P_K,{\mathcal{N}}_K)$ 描述：$K$ 是单元，$P_K$ 是单元上的局部函数空间，${\mathcal{N}}_K$ 是一组线性自由度。关键要求是单值可解性：给定这些自由度，必须唯一确定 $P_K$ 中的函数。若自由度为 $N_1,\dots ,N_m$，相应局部基满足 $N_i({\varphi }_j)={\delta }_{ij}$。

局部函数还要按所需的连续性拼接。连续分片一次元拼成 $H^1$ 协调空间；需要法向通量连续的空间采用另一种拼接规则。有限元的共同结构是“局部空间与自由度”，而不是“所有场量都在所有单元交界处连续”。

帽函数只在邻近单元上非零；共享单元产生非对角耦合，共享节点的局部贡献在组装时相加。

从局部积分到全局矩阵 ​

把 $u_h=\sum _j{U}_j{\varphi }_j$ 代入弱问题，再依次取 $v_h={\varphi }_i$，得到

只有支撑相交的两个基函数才可能产生非零矩阵元。因此，在局部连接数有界的低阶网格上，每行通常只有少量非零项。稀疏性来自基函数的局部支撑，不是求解器额外施加的近似。

在长度为 $h_e$ 的一维单元 $K_e=[x_e,x_{e+1}]$ 上，两个局部形状函数是

对常系数 $-u^{″}$，局部刚度矩阵与质量矩阵分别为

刚度矩阵积分的是导数乘积，质量矩阵积分的是函数乘积。后者会出现在含 $cu$ 的反应项中，也会出现在演化问题的 $M\dot{U}+AU=b$ 中；它并不是每一个静态 Poisson 系统都要额外相加的矩阵。

组装就是把局部编号映射到全局编号，并把同一全局位置收到的贡献相加。共享节点的系数只存一份，但相邻两个单元都会向它对应的矩阵行和列贡献积分。对稀疏矩阵而言，这是一项“累加”操作，而不是后来的单元覆盖前面的单元。

二维三角形同样如此。例如参考三角形顶点为 $(0,0),(1,0),(0,1)$，线性形状函数为 $1-x-y,x,y$。其梯度在单元内为常向量，所以单位扩散系数下

这个小矩阵的每行和为零，因为常函数的梯度为零。局部矩阵因此可以是奇异的；施加适当边界条件后，组装出的全局内部自由度系统才成为正定系统。把局部奇异误判为整体不可解，会混淆这两个层次。

补充示意图

一个从网格到解的完整例子 ​

取 $f\equiv 1$，把 $[0,1]$ 等分成三段，$h=1/3$。内部节点是 $x_1=1/3,x_2=2/3$。每个单元的荷载向量为

组装并代入两端零边界值后，方程只有两个未知数：

由对称性 $U_1=U_2$，因而 $U_1=U_2=1/9$。数值解是依次连接

的折线。精确解为 $u(x)=x(1-x)/2$，在这两个内部节点上也恰好取 $1/9$。

这种节点吻合不代表整个函数精确：在中间单元上，$u_h$ 是常数 $1/9$，精确解却仍是一段抛物线。它也不是任意高维问题或任意变系数问题都具有的性质。这个例子展示的是，未知系数的准确程度与单元内函数的准确程度需要分别考察。

节点完全吻合时，函数还有多少误差 ​

将刚才的三单元计算推广到 $N$ 个等长单元，$h=1/N$。内部方程为

把 $U_i=x_i(1-x_i)/2$ 代入：二次多项式的二阶中心差恰好给出 $2U_i-U_{i-1}-U_{i+1}=h^2$，所以每个方程都成立。刚度矩阵正定，解唯一，故 $u_h=I_{h}u$，即数值解正好是精确解的节点插值。这个结论在这里由常系数一维模型和精确积分得到，并不能据此推断一般变系数或高维模型也逐节点精确。

在任一单元 $K=[a,a+h]$ 上令 $t=x-a$。精确解与弦线之差是一个两端为零、二阶导数为 $-1$ 的二次函数，因此

误差本身是一段非负小弧线，误差导数则从 $h/2$ 线性下降到 $-h/2$。由此直接积分：

共有 $N=1/h$ 个单元，平方误差先相加、再开平方，得到三个不同范数中的精确答案：

最后一个最大值在每个单元中点取得。对 $N=3$，三项分别为 $1/\sqrt{108}$、$1/\sqrt{9720}$、$1/72$。节点误差全为零，能量误差与函数误差却都严格为正：只打印节点值无法检验整个函数的精度。

零节点误差仍对应非零的函数误差。左图比较 $u=x(1-x)/2$ 与 $N=3$ 的连续分片一次解；右图将单元坐标归一化为 $s=t/h$，展示 $e/h^2=s(1-s)/2$ 和 $e^{\prime }/h=1/2-s$，两条纵轴分别保留各自的尺度。

加密表与可复算任务 ​

下表由上述精确公式计算，不依赖某个软件求解器的输出。把步长减半，能量误差减半，$L^2$ 误差变为四分之一。

可以用这个算例检查一次独立实现：逐单元累加 $A^{(e)}$ 和 $b^{(e)}$，消去两端自由度，分别求解 $N=3,6,12$ 的系统；再在单元内部积分 $(u-u_h{)}^2$ 与 $(u^{\prime }-u_h^{\prime }{)}^2$，比较相邻网格误差比以及 ${\log }_2(E_h/E_{h/2})$。前者应得到 $2$ 与 $4$，后者应得到 $1$ 与 $2$。这里函数误差平方是四次多项式，三点 Gauss 积分即可精确积分；导数误差平方是二次多项式，两点即可。若只用节点样本估计误差，会错误地报告零误差。

完成任务时还应解释：$AU-b=0$ 只表示已准确解出离散方程，它不会消除上表中的离散误差。由此才能把“组装正确”“代数求解准确”和“逼近足够精细”三个判断分别核对。

不知道精确解时，怎样决定在哪里加密 ​

上面的误差表使用了已知精确解，适合核验实现；实际问题通常无法直接计算 $u-u_h$。残差后验估计公理库有限元残差后验估计Finite element residual estimator · Finite element a posteriori residual estimator从单元残差与斜率跳跃构造可计算指标，在一维模型中给出显式常数，再以二维准插值和泡函数证明可靠性与含振荡的局部效率。改为读取当前折线的单元残差与节点斜率跳跃，并在同一个零端点 Poisson 模型中证明能量误差上界。它还区分可靠性与含数据振荡的效率，避免把可计算指标直接当作真误差。

有了局部指标，就可以把新增自由度放到最需要的位置。自适应有限元的完整任务公理库自适应有限元的估计、标记与加密Adaptive finite element method · AFEM先复算一维局部加密，再在固定二维Poisson模型中证明估计量缩减、能量正交及加权收缩，并解释误差不变的边界加密例子。使用只作用于左四分之一区间的荷载：从四个单元出发，只二分左端单元便取得均匀八单元的相同能量误差。这个任务承接本页的组装与求解，增加的是“估计、标记、重解和比较”的决策环节。

边界、几何与求解：三个不能省略的环节 ​

边界条件进入哪里 ​

Dirichlet 条件指定未知函数本身的值，因此要限制试探空间，或把边界自由度作为已知数。若把全局系统按内部 $I$ 与边界 $B$ 分块，并给定 $U_B=g_B$，则实际求解

这样既把非零边界贡献移入右端，也保留了 $A_{II}$ 的对称性。只把边界行改成单位行而保留对应列，虽然可能得到同一个解，却会破坏整体矩阵的对称外观，影响后续求解器选择。

Neumann 条件指定边界通量，它通常通过分部积分留下的边界积分进入右端。纯 Neumann Poisson 问题还需要源项与边界通量满足兼容条件，解只确定到一个常数；添加平均值约束是去除这个自由度，而不是“修复一个计算错误”。Robin 条件则可能同时贡献边界矩阵和荷载。

单元如何搬到真实几何上 ​

实际程序往往在参考单元上准备基函数和积分点，再通过映射 $F_K$ 搬到真实单元。对标量函数，梯度按 $J_K^{-T}$ 变换，体积元乘上 $|\det J_K|$。因此，几何映射不仅改变积分范围，也改变梯度方向和尺度。

畸形单元会使这些变换放大误差。形状正则性限制的正是这种退化；它不等于所有单元必须同样大小。局部加密可以保留良好形状，并把自由度集中到误差较大的位置。对 $H(\mathrm{div})$、$H(\mathrm{curl})$ 元，还需使用匹配通量或切向分量的 Piola 型映射，不能照搬标量节点元的变换。

离散误差与代数误差分别控制 ​

数学上的 $u_h$ 假设单元积分和线性系统已经准确处理。数值积分不足、曲边几何近似或迭代求解提前停止，会引入另外的误差。实际计算应让这些误差与目标离散误差相称，而不是把线性残差压到极小后，就认定 PDE 解已经足够准确。

在准均匀、形状正则的网格上，对典型二阶椭圆问题和常用节点基，未经预条件的刚度矩阵条件数常随 $h^{-2}$ 增长。网格加密因此同时提高逼近能力并增加求解难度。稀疏存储公理库稀疏矩阵表示与运算Sparse matrix只存非零项及其位置，以稀疏模式组织矩阵运算、图结构、重排序与填充成本。解决的是存储与矩阵运算组织；预条件、迭代法和多重网格解决的是更深一层的求解效率。[3]

高阶单元还可能有只在单元内部耦合的未知量。Schur 补与静态凝聚公理库Schur 补与静态凝聚Schur complement · Static condensation把内部变量的精确响应折入界面矩阵与右端，推导 Schur 补、恢复公式及最小能量性质，并逐项凝聚五节点链。先精确消去这些内部自由度，把它们的响应折入界面矩阵和右端，再回代恢复。这与给定 Dirichlet 边界值后直接移项不同：被凝聚的变量原先仍是未知数。

为什么它会逼近原来的函数 ​

协调空间 $V_h\subset V$ 保留连续问题的测试结构。若双线性型的连续常数为 $M$、强制常数为 $\alpha >0$，Galerkin 方法中的 Céa 证明公理库Galerkin 方法Galerkin method · Galerkin discretization在有限维试探与检验空间中离散连续变分问题，并用 Galerkin 正交、Céa 准最优性及稳定条件组织误差。给出

有限元负责构造具有局部支撑、又有良好逼近能力的 $V_h$。对称正定问题中，离散解就是能量内积下的最佳逼近；对本条一维 Poisson 模型，Galerkin 条目已经用单元斜率平均证明 $|u-I_{h}u{|}_{H^1}\le h\Vert u^{″}{\Vert }_2$，再用零端点对偶问题证明 $\Vert u-u_h{\Vert }_2\le h^2\Vert u^{″}{\Vert }_2$。本条的精确公式给出这两种阶的实际常数，把抽象逼近估计与一张可复算的网格接在了一起。

在 Poisson 例子里，同一个解还最小化势能

有限元只是在“全部允许位移”中选解，改成在“这张网格能表达的位移”中选解。这是矩阵正定性、能量最小化与误差正交性背后的同一个结构。[2]

对形状正则网格上的连续一次元，若解有 $H^2$ 正则性，通常可得能量误差 $O(h)$；若相应对偶问题也有足够正则性，则可进一步得到 $L^2$ 误差 $O(h^2)$。更高的 $L^2$ 阶不是仅凭“一次多项式”四个字就能推出的。区域的凹角、系数跳跃或奇异荷载都可能降低正则性，从而改变实际收敛阶。椭圆内部正则性公理库椭圆内部正则性Interior elliptic regularity · 椭圆方程内部H2正则性对平方可积右端的Laplace弱方程，用能量截断和Fourier逆乘子证明内部二阶正则性，并以凹角奇性区分内部估计与全局边界估计。中的 $3\pi /2$ 凹角算例给出零边界、$f\in L^2$ 而 $u\notin H^2$ 的明确计算；因此只有内部二阶估计，还不足以代入本段要求整个区域 $H^2$ 的能量误差界。

二维一次元：把参考元估计搬回网格 ​

现在固定有界 Lipschitz 多边形区域 $\mathrm{\Omega }\subset {\mathbb{R}}^2$、齐次 Dirichlet Poisson 问题，以及一族协调、形状正则的三角剖分。记 $d_K=\mathrm{diam}K$，$d_h=\underset{K}{max}{d}_K$。本节假设 $u\in H^2(\mathrm{\Omega })\cap H_0^1(\mathrm{\Omega })$；下面会看到，这项正则性用于节点插值，不能从弱解存在性自动获得。

对顶点 $z_0,z_1,z_2$，写

若 $\hat{v}=v\circ F_K$，链式法则和变量替换给出

形状正则性使 $\Vert B_K\Vert \le C{d}_K$、$\Vert B_K^{-1}\Vert \le C/d_K$ 且 $|\det B_K|\asymp d_K^2$，常数只依赖允许的形状。因而

这里第二导数半范数可取 Hessian 的 Frobenius 范数；等价范数只改变固定常数。二维中梯度的一次缩放恰好与面积开方抵消，二阶导数则多留下一个 $d_K$。[5]

还需证明参考元上确实有可搬运的估计。令 $\hat{I}$ 为三顶点线性插值。对 $\hat{u}\in H^2(\hat{K})$，取

$\hat{u}-p$ 及其两个一阶导数的平均值均为零。先对梯度、再对函数应用参考三角形上的 Poincaré 不等式，得到

固定二维三角形上的 $H^2\hookrightarrow C^0$ 使三个点值有界；三个固定基函数的 $H^1$ 范数也有限，所以 $|\hat{I}w{|}_{H^1}\le C\Vert w{\Vert }_{H^2}$。又因 $\hat{I}p=p$，

物理元与参考元的插值在三个顶点取同样的值，因此 $(I_{K}u)\circ F_K=\hat{I}\hat{u}$。代入缩放关系，逐元平方求和，即得

齐次边界使 $I_{h}u\in V_h$。对任意 $v_h\in V_h$，Galerkin 正交给出

取 $v_h=I_{h}u$ 并约去非零误差，得到完整的二维先验结论

证明只需形状正则，不要求所有单元大小可比。它也不说二维离散解等于节点插值：插值这里只提供一个可比较的候选函数。解仅在 $H_0^1$ 中时，二维点值一般不足以定义这种插值；二维残差估计公理库有限元残差后验估计Finite element residual estimator · Finite element a posteriori residual estimator从单元残差与斜率跳跃构造可计算指标，在一维模型中给出显式常数，再以二维准插值和泡函数证明可靠性与含振荡的局部效率。改用局部平均准插值，直接从当前离散解证明可靠性，不需要未知解具有全局 $H^2$ 正则性。

与其他离散方法的联系 ​

有限差分法公理库偏微分方程有限差分法Finite-difference method for PDE · Finite-difference PDE discretization · Difference stencil由 Taylor 展开构造 PDE 差分模板，完整处理编号与边界，并用小系统、离散能量和制造解检验计算。从点上的导数近似组织方程，有限元从弱问题与局部函数空间组织方程。在规则网格和简单系数下，两者有时会给出相同或成比例的矩阵；这种重合并不会消除两种推导方式在边界处理、几何和误差分析上的区别。

用于时间依赖问题时，有限元先产生 $M\dot{U}+AU=b(t)$，再由线法公理库线法：先空间离散再时间积分Method of lines · MOL · Semi-discretization in space先把 PDE 变为空间半离散 ODE，再分析时间推进、边界输入、刚性、质量矩阵和空间与时间误差的配合。选择时间积分器。用于自适应计算时，后验误差估计指导“在哪里加密”，而 Céa 与插值估计主要解释“这样的空间为何能够逼近”。这些都是同一条主线的延伸：构造合适的有限维空间，再控制在这个空间中求解所造成的误差。