偏微分方程有限差分法

有限差分法把偏微分方程中的导数，替换成邻近网格点函数值的线性组合。一个点只与少数邻点发生联系；把整张网格上的这些局部关系连起来，就得到一个稀疏代数系统，或一个随时间演化的大型常微分方程组。

真正的 PDE 差分格式由内部离散、边界闭合和未知量组织共同组成。写出一个漂亮的内部差分公式，还没有完成一个可求解的边值问题。

形式陈述

从导数近似得到局部模板

设网格间距为 $h$ 。在 $x$ 两边做 Taylor 展开：

\begin{aligned} u(x+h)&=u(x)+hu'(x)+\frac{h^2}{2}u''(x) +\frac{h^3}{6}u'''(x)+\frac{h^4}{24}u^{(4)}(x)+\frac{h^5}{120}u^{(5)}(x)+O(h^6),\\ u(x-h)&=u(x)-hu'(x)+\frac{h^2}{2}u''(x) -\frac{h^3}{6}u'''(x)+\frac{h^4}{24}u^{(4)}(x)-\frac{h^5}{120}u^{(5)}(x)+O(h^6). \end{aligned}

相加消去奇数阶项，得到中心二阶差分

\frac{u(x-h)-2u(x)+u(x+h)}{h^2} =u''(x)+\frac{h^2}{12}u^{(4)}(x)+O(h^4),

这里假定六阶导数在所讨论的邻域内有界；五阶项相消后，除以 $h^2$ 的余项为 $O(h^4)$ 。这说明该公式在光滑函数上具有二阶截断精度；它还没有说明求解整个离散系统后，数值解的误差也是二阶。后者需要稳定性估计。[1]

考虑矩形区域内的 Poisson 方程

-\Delta u=f,\qquad u|_{\partial\Omega}=g.

在横纵间距都为 $h$ 的规则网格上，分别近似 $u_{xx}$ 和 $u_{yy}$ ，得到

\frac{4U_{i,j}-U_{i-1,j}-U_{i+1,j}-U_{i,j-1}-U_{i,j+1}}{h^2} =f_{i,j}.

中心点和东、西、南、北四个邻点构成五点模板。记忆模板时，最好同时记住它来自 $-\Delta$ ：中心系数为正，四个邻点系数为负。若写的是 $\Delta$ ，整组符号就要反转。

横纵步长不同的时候，不能再直接使用同一个 $1/h^2$ ：

\frac{2U_{i,j}-U_{i-1,j}-U_{i+1,j}}{h_x^2} + \frac{2U_{i,j}-U_{i,j-1}-U_{i,j+1}}{h_y^2} =f_{i,j}.

模板是离散算子的局部表达，不是脱离网格几何的一幅固定图案。

负 Laplace 的五点模板：中心系数为 4，四邻点系数为 −1，统一除以 h²。编号把局部连接转换成稀疏矩阵行。

直觉

边界值怎样进入方程

若某个邻点位于 Dirichlet 边界，其值已经由 $g$ 给出，就应从未知量中去掉。例如 $U_{0,j}=g_{0,j}$ 时，紧邻左边界的一行变为

\frac{4U_{1,j}-U_{2,j}-U_{1,j-1}-U_{1,j+1}}{h^2} =f_{1,j}+\frac{g_{0,j}}{h^2}.

右端出现正号，是把左端的 $-g_{0,j}/h^2$ 移过去的结果。边界不会因为没有进入未知向量就消失；它通过右端继续影响内部解。

给内部节点选定全局编号，例如先沿一行从左到右，再逐行向上，便可写成

A_hU=b_h.

编号改变会置换矩阵行列，却不改变离散问题本身。每个内部节点最多连接四个内部邻点，因此对 $N$ 个未知量，五点格式的非零元数不超过 $5N$ 。采用稀疏矩阵存储，才与这个结构相匹配；构造一个 $N\times N$ 的稠密数组会丢掉差分模板带来的主要计算优势。

一个可以手算的二维系统

在单位正方形上解 $-\Delta u=1$ ，四边取零。令 $h=1/3$ ，内部只有四个点，按

U=(U_{1,1},U_{2,1},U_{1,2},U_{2,2})^T

编号。方程为

9\begin{pmatrix} 4&-1&-1&0\\ -1&4&0&-1\\ -1&0&4&-1\\ 0&-1&-1&4 \end{pmatrix}U = \begin{pmatrix}1\\1\\1\\1\end{pmatrix}.

对称性给出四个分量相等。设它们都为 $c$ ，每行便成为 $9(4c-c-c)=1$ ，所以 $c=1/18$ 。这里的“两个邻点”指两个内部邻点；另外两个邻点的零边界值已经被消去了。

这个例子虽然很小，却包含了完整实现中的全部骨架：模板、编号、边界贡献、稀疏结构和解方程。加密网格只是增加同样的局部关系，而不是换一种基本原理。

为什么矩阵可解，为什么误差不会随意放大

对上述零 Dirichlet Poisson 离散， $A_h$ 是对称正定矩阵。将网格函数在边界补零后，二次型可以写成相邻节点差值平方之和：

U^TA_hU =\frac1{h^2}\sum_{\{p,q\}}(U_p-U_q)^2,

其中每条内部网格边或连接内部与边界的网格边只计一次。右端非负；若为零，网格函数在连接区域上为常数，再由零边界推出这个常数只能为零。因此系统唯一可解。

这一结构也解释了数值解的“平滑倾向”：相邻点差异越大，离散能量越大。不过，正定性只说明每个固定网格可逆；要证明随网格加密仍然稳定，还需在选定的离散范数中控制 $A_h^{-1}$ ，不能只检查一次矩阵求逆成功。

设精确解在网格上的取值为 $u_h^*$ ，离散残差为

r_h=A_hu_h^*-b_h.

若数值解满足 $A_hU=b_h$ ，定义 $e_h=U-u_h^*$ ，则

A_he_h=-r_h,\qquad \|e_h\|_h\le\|A_h^{-1}\|_h\,\|r_h\|_h.

这一步把局部截断误差与全局解误差接上了。对规则区域上的稳定二阶格式，若解和边界足够光滑，残差的 $O(h^2)$ 才能转化为相应范数下的 $O(h^2)$ 解误差。完整的逻辑见一致性、稳定性与收敛性。

用一个精确模态核对二阶精度

在单位正方形上取

u(x,y)=\sin(\pi x)\sin(\pi y),\qquad f=2\pi^2u,

边界仍为零。它在等距网格上的取值向量 $v$ 恰好也是五点矩阵的特征向量：

A_hv=\lambda_hv,\qquad \lambda_h=\frac8{h^2}\sin^2\frac{\pi h}{2} =2\pi^2-\frac{\pi^4}{6}h^2+O(h^4).

于是数值解可以直接写成

U=\frac{2\pi^2}{\lambda_h}v =\left(1+\frac{\pi^2}{12}h^2+O(h^4)\right)v.

这里误差主要体现为振幅稍大，且相对误差按 $h^2$ 缩小。网格步长减半后，进入渐近区间时误差应约为原来的四分之一。用这个已知解检查程序，比仅看热力图是否“足够平滑”更能发现符号、尺度或边界装配错误。

补充示意图

网格五点 stencil 的局部耦合映射为全局稀疏矩阵的一行 — 五点差分模板到稀疏矩阵行

例子与边界

演化问题：空间模板还要配时间推进

对热方程 $u_t=\kappa u_{xx}$ ，中心差分先给出

\dot U_i(t)=\kappa\frac{U_{i-1}(t)-2U_i(t)+U_{i+1}(t)}{h^2}.

这仍是连续时间的方程组。再采用显式 Euler，才得到完全离散格式

U_i^{n+1}=(1-2r)U_i^n+rU_{i-1}^n+rU_{i+1}^n, \qquad r=\frac{\kappa\tau}{h^2}.

当 $0\le r\le1/2$ 时，新值是旧值的非负加权平均。这既提供了离散最大值原理的直觉，也与对应常系数格式的 Fourier 稳定性阈值一致。若把时间方法改成后向 Euler，同一个空间模板则产生一个每步都要解的线性系统，而不再有这一显式时间步限制。[2]

因此，“二阶中心差分”只说明空间部分。讨论整个 PDE 算法时，应同时说明时间方法、步长关系、边界更新和误差范数。线法专门组织这两个离散阶段；von Neumann 分析则解释不同频率的误差如何被完全离散格式放大。

边界闭合与规则网格之外

Neumann 条件指定的是导数，通常不能像 Dirichlet 值那样直接删除未知量。例如在左端可使用二阶单边公式

\frac{-3U_0+4U_1-U_2}{2h}=g_0,

也可以引入幽灵点并结合边界导数消元。两种做法都必须与内部方程一起检查稳定性和精度。内部是二阶，边界处却随意换成一阶近似时，整体误差阶不一定还保持二阶。

对于纯 Neumann 问题，连续模型的常数自由度通常会保留为离散矩阵的零空间。需要检查离散兼容条件并固定平均值，而不是任意删掉一条方程后便认为问题已被完整处理。

非均匀网格上，左右步长不同，原来的中心模板需要重推。曲边区域上，阶梯状边界可能让几何误差主导总误差；高阶内部公式不能自动弥补低精度几何表示。映射网格、紧致差分或嵌入边界方法会改变局部结构和分析条件，不能只替换一个系数继续沿用原来的结论。

推论与应用

变系数问题：保留通量结构

对

-(\kappa(x)u'(x))'=f(x),\qquad \kappa(x)>0,

直接写 $-\kappa_i(U_{i-1}-2U_i+U_{i+1})/h^2$ 一般不对，因为原算子还包含系数变化的影响。更自然的做法是先离散界面通量，再计算通量差：

-\frac{ \kappa_{i+1/2}(U_{i+1}-U_i) -\kappa_{i-1/2}(U_i-U_{i-1}) }{h^2}=f_i.

同一个界面通量，在左侧单元的方程中流出，在右侧单元的方程中流入，符号相反。这样组织方程会保留局部守恒结构，也使正扩散系数对应的矩阵具有合适的对称性。

若界面两侧分别是系数 $\kappa_L,\kappa_R$ 的等长半段，稳态通量通过两个半段时应相同。两段“阻力”相加，得到有效系数

\kappa_{i+1/2} =\frac{2\kappa_L\kappa_R}{\kappa_L+\kappa_R}.

这就是调和平均的来源。它并不是任意场景下都必须使用的经验公式，而是对应这项界面位置与分段常系数模型的正确等效系数。若两段长度不同，应改用按长度加权的调和平均。系数跃变时，通量连续往往比点值导数连续更接近原问题的物理要求。

如何判断一次计算是否可信

一条有用的核对路线是：先用已知光滑解制造右端和边界，检查残差与观测收敛阶；再检查完整离散算子的稳定性；最后把线性或非线性求解误差控制到不掩盖离散误差的量级。三者针对的是不同环节。

例如把线性系统解到极小残差，只说明已经非常接近这张网格的离散解。它不能排除网格太粗、边界处理不当或模型本身写错。反过来，若迭代残差过大，即使格式理论上二阶，也可能在加密实验中看不到预期的四倍误差下降。

有限元方法从弱形式组织类似的局部耦合，有限差分从导数或通量近似组织耦合。对简单规则网格，它们的矩阵可能非常接近；选择方法时真正要比较的是几何表示、守恒或能量结构、边界处理，以及是否有与之匹配的稳定性和误差分析。

模板中的导数近似来自有限差分公式，但局部公式的精度还不能代替完整边值问题的稳定性分析。

当解足够光滑且几何规则时，可以用全局插值代替局部差分模板。谱配点法在节点处满足微分方程，再以边界方程补齐系统；其稠密耦合、逼近精度和稳定性需要与局部网格方法分别比较。

参考资料

[1] Randall J. LeVeque, Finite Difference Methods for Ordinary and Partial Differential Equations, SIAM, 2007，差分公式、边值问题与演化方程的统一处理。

[2] Tobin A. Driscoll and Richard J. Braun, Fundamentals of Numerical Computation, Method of lines，空间差分与时间积分的具体连接。

[3] John C. Strikwerda, Finite Difference Schemes and Partial Differential Equations, 2nd ed., SIAM, 2004，稳定性、边界条件及收敛分析。

[4] Lloyd N. Trefethen, Finite Difference and Spectral Methods for Ordinary and Partial Differential Equations, 1996，差分算子与 Fourier 模态分析的作者教材。

形式陈述 ​

从导数近似得到局部模板 ​

边界值怎样进入方程 ​

一个可以手算的二维系统 ​

为什么矩阵可解，为什么误差不会随意放大 ​

用一个精确模态核对二阶精度 ​

演化问题：空间模板还要配时间推进 ​

边界闭合与规则网格之外 ​

变系数问题：保留通量结构 ​

如何判断一次计算是否可信 ​

反向链接