有限差分法

有限差分法（finite difference method）用网格点上函数值的线性组合近似导数，从而把微分方程转成代数方程或时间递推。方法由差分模板、边界处理和求解方式共同构成；一个导数公式精度很高，并不自动保证整个边值或初值问题稳定收敛。

相邻三个点近似曲率

在均匀节点 $x_i=ih$ 上，Taylor 展开给出 $$u(x_i\pm h)=u(x_i)\pm h{u}^{\prime }(x_i)+\frac{h^2}{2}{u}^{″}(x_i)\pm \frac{h^3}{6}{u}^{‴}(x_i)+\frac{h^4}{24}{u}^{(4)}({\xi }_{\pm }).$$ 两式相加，奇数项抵消，因此 $$D_h^2{u}_i=\frac{u_{i-1}-2u_i+u_{i+1}}{h^2}=u^{″}(x_i)+O(h^2).$$ 若四阶导数在相邻区间有界为 $M$，差值不超过 $M{h}^2/12$。它说明二阶中心差分测量局部曲率；对线性函数结果为零，对二次函数恰好精确。

例如 $u=x^2$，三点公式的分子为 $(x-h{)}^2-2x^2+(x+h{)}^2=2h^2$，除以 $h^2$ 得二。若误将分母写为 $h$，即使曲线图看起来平滑也会算错量纲和结果。

一个正弦边值问题

求 $-u^{″}={\pi }^2\sin \pi x$、$u(0)=u(1)=0$，真解为 $u=\sin \pi x$。取 $h=1/4$，三个内部未知量满足 $$16\left(\begin{array}{ccc}2& -1& 0\\ -1& 2& -1\\ 0& -1& 2\end{array}\right)\left(\begin{array}{c}{U}_1\\ U_2\\ U_3\end{array}\right)={\pi }^2\left(\begin{array}{c}\sqrt{2}/2\\ 1\\ \sqrt{2}/2\end{array}\right).$$ 离散正弦仍是矩阵的特征向量，因为 $$2\sin (\pi x_i)-\sin (\pi (x_i-h))-\sin (\pi (x_i+h))=4{\sin }^2(\pi h/2)\sin (\pi x_i).$$ 因此不用直接求三阶逆矩阵，就能写出 $$U_i=\frac{{\pi }^2{h}^2}{4{\sin }^2(\pi h/2)}\sin (\pi x_i).$$ 中点近似为 ${\pi }^2/[64{\sin }^2(\pi /8)]\approx 1.05302929$，比真值一略大。这个误差来自微分算子的近似，而非线性系统未解准。

正弦问题取4和8个等长区间。连接节点的线只帮助阅读，有限差分未知量本身首先是网格点值。

从截断误差到全局误差

局部公式的一致性还要配合稳定性。记离散算子 $A_h=-D_h^2$，边界为零，真解采样为 $u_i$。若 $\tau =A_{h}u-f$，数值误差 $e=U-u$ 满足 $A_{h}e=-\tau$。

这个一维算子有离散比较原理：若端点 $v\ge 0$、内部 $A_{h}v\ge 0$，则全部 $v\ge 0$。否则取内部负的最小值，$2v_i-v_{i-1}-v_{i+1}\le 0$，与非负要求合起来迫使两个邻居也等于该最小值。重复向边界传播，就与端点非负矛盾。

取屏障函数 $B_i=x_i(1-x_i)/2$，二次函数的差分精确性给出 $A_{h}B=1$，且 $0\le B\le 1/8$。设 $T=\Vert \tau {\Vert }_{\mathrm{\infty }}$，则 $A_h(TB\pm e)=T\mp \tau \ge 0$，边界为零。比较原理推出 $$|e_i|\le T{B}_i\le \frac{1}{8}\Vert \tau {\Vert }_{\mathrm{\infty }}\le \frac{M{h}^2}{96}.$$ 这完成了从局部二阶缺陷到全局二阶网格误差的证明，并指出使用了零 Dirichlet 边界、均匀一维网格、比较原理及四阶导数控制。Netlib 的 Poisson 模型说明了这类差分矩阵与迭代求解器的联系。

有解变式为 $-u^{″}=2$，仍取零端点。真解 $x(1-x)$ 是二次函数，截断误差恰为零，所以任意均匀网格的内部节点解都精确等于它。只有在网格点上精确，不意味着任意节点之间的插值也等于二次真解。

可运行的边值解算器

以下代码对任意有限采样右端使用二阶中心差分和零边界，返回节点及节点值，不修改输入。为突出方程结构，采用 NumPy 稠密求解，成本约为 $O(N^3)$ 时间和 $O(N^2)$ 空间；三对角专门解法则只需 $O(N)$。

Python 3

import numpy as np


def fd_poisson_1d(forcing, intervals):
    if not isinstance(intervals, int) or intervals < 2:
        raise ValueError("at least two integer intervals required")
    x = np.linspace(0, 1, intervals+1); h = 1/intervals
    rhs = np.array([forcing(float(t)) for t in x[1:-1]], dtype=float)
    if rhs.shape != (intervals-1,) or not np.isfinite(rhs).all():
        raise ValueError("forcing must return finite scalars")
    m = intervals-1
    A = 2*np.eye(m)-np.eye(m, k=1)-np.eye(m, k=-1)
    values = np.zeros(intervals+1)
    values[1:-1] = np.linalg.solve(A, h*h*rhs)
    if not np.isfinite(values).all():
        raise ArithmeticError("nonfinite difference solution")
    return x, values


x, u = fd_poisson_1d(lambda x: np.pi**2*np.sin(np.pi*x), 4)
print(u[2])  # 1.0530292875...
x, u = fd_poisson_1d(lambda x: 2, 4)
print(u)     # [0, 0.1875, 0.25, 0.1875, 0]

在时间依赖问题中，空间模板还要与时间推进匹配：热方程显式更新受 $\Delta t=O(h^2)$ 约束，输运差分需考虑传播方向和 CFL 条件。非均匀网格不能原样使用同一个三点系数，曲线坐标还会产生几何系数。边界附近单边公式的阶数与稳定性也应独立核验。