有限元法

有限元法（finite element method，FEM）把连续方程的弱形式限制到由网格上局部基函数张成的有限维空间，求得分片函数近似。它既是一种局部装配方式，也是一种函数空间近似方法；网格节点上的数值最终应与基函数一起解释，不能仅当作一张离散数据表。

从拉紧细弦的平衡方程开始

考虑单位区间上的教学边值问题 $- u^{″} = 1$ 、 $u (0) = u (1) = 0$ 。真解是 $u = x (1 - x) / 2$ 。先只允许在 $x = 1 / 2$ 折一次的连续分段直线。它由一个未知高度 $U$ 决定，可写成 $u_{h} = U ϕ$ ，其中帽函数 $ϕ$ 在两端为零、中点为一。

把方程乘以两端为零的测试函数 $v$ ，分部积分后得到 $\int_{0}^{1} u^{'} v^{'} d x = \int_{0}^{1} v d x .$ 在一维空间中取 $v = ϕ$ ，左边为 $U \int (ϕ^{'})^{2} = 4 U$ ，右边为帽形面积 $1 / 2$ ，所以 $U = 1 / 8$ 。有限元解的中点恰等于真解，但单元内仍是一条折线，并不是原来的抛物线。

两段线性帽函数在中点取一，两端为零，左右导数为正负二 — 一个内部自由度对应一个帽函数。测试函数和试探函数取自同一有限维空间，弱方程成为单个代数方程。

弱形式与局部装配

对一般 $f \in L^{2} (0, 1)$ ，弱问题是在 $H_{0}^{1} (0, 1)$ 中求 $u$ ，使 $a (u, v) = ℓ (v)$ 对全部测试函数成立，其中 $a (u, v) = \int u^{'} v^{'}$ 、 $ℓ (v) = \int f v$ 。零端点条件与 Poincaré 不等式使导数内积具有所需的强制性，弱解存在唯一。

取节点 $0 = x_{0} < \dots < x_{N} = 1$ ，空间 $V_{h}$ 由连续分段线性且端点为零的函数构成，写 $u_{h} = \sum_{j = 1}^{N - 1} U_{j} ϕ_{j}$ 。对每个 $ϕ_{i}$ 测试，得到 $K U = F, K_{i j} = \int {ϕ_{j}}^{'} {ϕ_{i}}^{'}, F_{i} = \int f ϕ_{i} .$ 基函数只在相邻单元重叠，因此矩阵稀疏。宽度为 $h_{K}$ 的线性单元，其两个局部基函数导数为 $- 1 / h_{K}, 1 / h_{K}$ ，局部刚度为 $K^{(K)} = \frac{1}{h_{K}} (\begin{matrix} 1 & - 1 \\ - 1 & 1 \end{matrix}) .$ 若单元上 $f = c_{K}$ 为常数，局部载荷为 $c_{K} h_{K} (1, 1)^{T} / 2$ 。把共享节点的局部贡献相加，并固定边界自由度，便完成装配。MFEM 的官方 Poisson 示例明确展示网格、有限元空间、双线性型、线性型和边界条件的对应关系。^[1]

Galerkin 正交与最佳逼近的证明

真解与离散解都对 $v_{h} \in V_{h}$ 满足同一个弱方程，相减得到 $a (u - u_{h}, v_{h}) = 0.$ 这叫 Galerkin 正交。对任意 $w_{h} \in V_{h}$ ，写 $u - w_{h} = (u - u_{h}) + (u_{h} - w_{h})$ ，两部分在 $a$ 内积下正交。因此 $‖ u^{'} - {w_{h}}^{'} ‖_{L^{2}}^{2} = ‖ u^{'} - {u_{h}}^{'} ‖_{L^{2}}^{2} + ‖ {u_{h}}^{'} - {w_{h}}^{'} ‖_{L^{2}}^{2} \geq ‖ u^{'} - {u_{h}}^{'} ‖_{L^{2}}^{2} .$ 离散解是在所选空间中能量误差最小的函数。这不保证每个节点或每一点的误差分别最小；结论对应的是明确的导数二范数。若积分近似、非保形空间或迭代求解误差破坏了同一个弱方程，等式也需相应修正。

两单元和四单元的有限元折线逐渐接近二次真解 — 方程 −u″=1。均匀细化后节点数增加，分片直线能更准确表示曲率；本例节点精确并不代表单元内部没有误差。

主例在宽度为 $h$ 的每个单元上，导数误差等于 $- (x - m_{K})$ ，其中 $m_{K}$ 是中点。因此单元误差平方为 $h^{3} / 12$ ，均匀网格总误差为 $h / \sqrt{12}$ 。两单元时约 $0.144337567$ ；四单元时约 $0.072168784$ ，减少一半。对一般解，这种一阶能量收敛要求相应正则性；尖角奇性或粗糙系数会影响阶数和网格选择。

有解变式将右端改为二。线性方程右端翻倍，有限元节点值和真解都翻倍：两单元的中点为 $1 / 4$ ，能量误差为 $1 / \sqrt{12}$ 。这检验了装配中的单元长度因子，漏掉它会给出错误尺度。

非均匀一维网格的可运行装配

代码以每个单元上的常数右端为输入，使用精确局部积分和零端点，返回全部节点值，不修改输入。为清楚展示装配，使用 NumPy 稠密矩阵求解，时间约 $O (N^{3})$ 、存储 $O (N^{2})$ ；本例实际为三对角结构，可改用线性时间和空间的专门求解器。

Python 3

import numpy as np


def fem_poisson_1d(nodes, cell_loads):
    x = np.array(nodes, dtype=float, copy=True)
    loads = np.array(cell_loads, dtype=float, copy=True)
    if x.ndim != 1 or x.size < 2 or loads.shape != (x.size-1,):
        raise ValueError("one load per interval required")
    if not np.isfinite(x).all() or not np.isfinite(loads).all() or not (np.diff(x) > 0).all():
        raise ValueError("finite data and increasing nodes required")
    if x[0] != 0 or x[-1] != 1:
        raise ValueError("domain must be [0,1]")
    count = x.size
    K = np.zeros((count, count)); F = np.zeros(count)
    for i, (h, f) in enumerate(zip(np.diff(x), loads)):
        K[i:i+2, i:i+2] += np.array([[1, -1], [-1, 1]])/h
        F[i:i+2] += f*h/2
    result = np.zeros(count)
    if count > 2:
        result[1:-1] = np.linalg.solve(K[1:-1, 1:-1], F[1:-1])
    if not np.isfinite(result).all():
        raise ArithmeticError("nonfinite finite-element solution")
    return result.tolist()


print(fem_poisson_1d([0, 0.5, 1], [1, 1]))  # [0, 0.125, 0]
print(fem_poisson_1d([0, 0.25, 0.5, 0.75, 1], [1]*4))
# [0, 0.09375, 0.125, 0.09375, 0]

高维问题仍遵循局部弱形式、装配与边界条件，但单元映射、积分规则和函数空间会更复杂。弹性、流体或电磁问题还需要与物理量相匹配的空间；把所有变量都放进相同连续线性空间未必满足稳定性条件。误差估计与自适应网格加密则帮助判断当前网格应怎样改进。

参考资料

↑ MFEM authors. Example 1: Poisson Problem。

[mfem-1] MFEM authors. Example 1: Poisson Problem。

[1]