第 10 章有限差分与稳定性入门

第 10 章有限差分与稳定性入门 | 自在学

第 8 章里，我们已经按“上一层的三个温度算下一层的一个温度”推进过热方程。那个公式很短，真正需要花时间的却是公式之外的判断：空间网格变细，为什么计算反而可能爆掉？隐式方法不爆掉，为什么图像仍可能离真实解很远？

这一章把这些判断接起来。我们从导数近似的误差出发，写出实际需要求解的方程组，再看每一步的误差怎样传播。读完以后，你应该能够同时检查网格、边界、时间步和误差，而不只认出一个差分模板。

显式热格式的三个输入来自旧时间层，端点另按边界条件更新。

网格值与精确取样分开记

在 $0\le x\le L$ 上取

x_j=jh,\quad j=0,1,\ldots,N,\qquad h=L/N,

时间节点为 $t^n=n\Delta t$ 。本章用大写 $U_j^n$ 表示计算机得到的数值，用

v_j^n=u(x_j,t^n)

表示精确解在同一节点上的取样。两者通常不同；它们的差才是我们想控制的离散误差。

$j=0,N$ 是边界点， $1\le j\le N-1$ 是内部点。若给的是 Dirichlet 条件，端点值已经由题目规定，通常只把内部值列为未知数。如果边界随时间变化，每一层都要使用相应时刻的端点值，不能只在初始时刻设置一次。

对称取样为什么能消掉低阶误差

先离开时间变量，只看一个光滑函数 $v(x)$ 。在 $x$ 附近展开 $v(x+h)$ 与 $v(x-h)$ 。为明确下面余项的阶数，暂假设 $v$ 在邻域内有连续且有界的六阶导数。于是

v(x\pm h)=v(x)\pm hv'(x)+\frac{h^2}{2}v''(x) \pm\frac{h^3}{6}v'''(x)+\frac{h^4}{24}v^{(4)}(x) \pm\frac{h^5}{120}v^{(5)}(x)+O(h^6).

相减时，偶数次项消失。除以 $2h$ 后得到

D_xv(x):=\frac{v(x+h)-v(x-h)}{2h} =v'(x)+\frac{h^2}{6}v'''(x)+O(h^4).

相加后再减去 $2v(x)$ ，消掉的则是奇数次项：

D_{xx}v(x):=\frac{v(x+h)-2v(x)+v(x-h)}{h^2} =v''(x)+\frac{h^2}{12}v^{(4)}(x)+O(h^4).

两个近似都是二阶精度，但一个近似一阶导数，一个近似二阶导数。“二阶精度”说的是 $h$ 缩小时误差的阶，并不是所求导数的阶数。

对 v(x)=x⁴，在 x=1/2 处取 h=0.2，中心二阶差分为 3.08，精确导数为 3。

例如 $v(x)=x^4$ ，在 $x=1/2$ 处有 $v''=3$ ，而

D_{xx}v=3+2h^2.

$h=0.2$ 时得到 $3.08$ ，误差 $0.08$ ； $h=0.1$ 时得到 $3.02$ ，误差 $0.02$ 。步长减半，误差成为四分之一。这次结论还是精确等式，因为更高次导数已经为零。

换成 $v(x)=|x|$ ，在原点却有

D_xv(0)=0,\qquad D_{xx}v(0)=2/h.

第一个数不是原点的导数，因为原点根本不可导；第二个数还会随 $h\downarrow0$ 发散。对称公式照样可以算出数字，是否近似了一个存在的导数，需要另外检查。

实验：误差缩小的规律什么时候失效

从四次函数预设开始，把 $h$ 连续减半，核对误差比是否接近 $4$ 。切换到尖角，比较“差分算出了值”和“导数存在”这两句话。指数函数的极小步长预设展示另一种困难：相近浮点数相减会丢失有效数字，再除以 $h^2$ ，舍入误差可能被放大。因此减小 $h$ 不能无限改善结果。

离散 Laplacian 与边界方程组

一维中，二阶差分还可以写成

\Delta_hU_j=\frac{2}{h^2}\left(\frac{U_{j-1}+U_{j+1}}2-U_j\right).

它比较的是当前点与邻居平均值。当前点高于平均值时，离散 Laplacian 为负；低于平均值时为正。这与热方程中局部峰顶下降、谷底上升的方向一致。

二维正方形网格的五点模板来自两个方向相加：

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

四个邻点在上下左右，对角点并不参与这个模板。若两个方向的步长不同，则应分别除以各自的平方：

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

正方形网格的五点 Laplacian：中心权重为 −4，四个轴向邻点权重为 1，再除以 h²。

例题：先把三个内部未知量算出来

求 $-u''=1$ 、 $u(0)=u(1)=0$ 的差分近似，取 $h=1/4$ 。内部未知量为。

方程是负的二阶导数，因此内部行应写成

-U_{j-1}+2U_j-U_{j+1}=h^2.

若端点改为 $u(0)=a_0$ 、 $u(1)=b_0$ ，第一行出现 $-a_0$ ，移到右侧变成；最后一行同理加上。对一般源项，右端为

\begin{pmatrix}h^2f(x_1)+a_0\\h^2f(x_2)\\h^2f(x_3)+b_0\end{pmatrix}.

例如 $a_0=1,b_0=3,f=1$ 时，右端为 $(17,1,49)^{\mathsf T}/16$ ，解为

(U_1,U_2,U_3)=\left(\frac{51}{32},\frac{17}{8},\frac{83}{32}\right).

它们也正好等于 $1+2x+x(1-x)/2$ 的取样值。端点不是被丢掉了，而是进入右端，影响所有内部未知量。

已知端点的负项移到右端，变成正的边界贡献；内部矩阵列只对应未知内点。

这个线性系统为什么有唯一解

记 $T$ 为对角线 $2$ 、上下副对角线 $-1$ 的矩阵。对内部向量 $z=(z_1,\ldots,z_{N-1})$ ，补上，整理有限求和可得

z^{\mathsf T}Tz=\sum_{j=0}^{N-1}(z_{j+1}-z_j)^2.

若右侧为零，相邻值全部相同，再由两端为零得到 $z=0$ 。所以非零向量都有 $z^{\mathsf T}Tz>0$ ， $T$ 正定，从而可逆。这与上一章用梯度平方积分证明唯一性的思路相近，只是积分变成了有限求和。

实验：残差为零，为什么还有误差

在二次解预设中核对三行方程，再改动端点，观察右端改变了哪几行。实验会实际求解三对角系统。

换成四次解，离散方程残差仍接近零，但数值节点开始偏离精确曲线。选择“求解 N 与 2N，比较误差”，比较同一源项、同一边界下的两组结果。这里有两种不同的检查：线性系统残差检验方程有没有解对，离散误差检验所解的离散方程与连续问题相差多少。

显式热格式：全部读取旧时间层

对 $u_t=\kappa u_{xx}$ ，前向时间差分和中心空间差分给出

\frac{U_j^{n+1}-U_j^n}{\Delta t} =\kappa\frac{U_{j+1}^n-2U_j^n+U_{j-1}^n}{h^2}.

令 $r=\kappa\Delta t/h^2$ ，则

\boxed{U_j^{n+1}=rU_{j-1}^n+(1-2r)U_j^n+rU_{j+1}^n.}

所有右端值都来自第 $n$ 层。实际写程序时，用一个新数组保存下一层；不要算完 $U_1^{n+1}$ 后，马上拿它去代替第二行所需的 $U_1^n$ 。那样就改变了算法。

当 $0\le r\le1/2$ 时，三个权重非负且和为 $1$ ，新值落在三个旧值的最小值与最大值之间。若边界值也在固定区间 $[m,M]$ 内，逐层归纳可知所有数值都留在 $[m,M]$ 。这就是该格式的离散最大值性质。

对具有相同边界的两次计算，差向量端点为零。同一个凸组合还给出

\|e^{n+1}\|_\infty\le\|e^n\|_\infty,

其中 $\|e\|_\infty=\max_j|e_j|$ 。它说明初始误差不会因为时间推进而被放大，是本章将使用的一种稳定性。

例题：算出比值，还要把权重代对

设 $\kappa=0.2$ 、 $h=0.1$ 。采用与网格无关的充分条件 $r\le1/2$ ，时间步应满足

r=20\Delta t\le\frac12,\qquad \Delta t\le0.025.

取 $\Delta t=0.02$ ，有 $r=0.4$ ，所以更新式为

\boxed{U_j^{n+1}=0.4U_{j-1}^n+0.2U_j^n+0.4U_{j+1}^n.}

若三个旧值依次为 $1,4,2$ ，新值是 $0.4+0.8+0.8=2$ 。用 $0.2,0.6,0.2$ 虽然也得到一个加权平均，却对应 $r=0.2$ ，已经不是这道题给出的步长。

若把 $\Delta t$ 改成 $0.04$ ，则 $r=0.8$ ，中心权重变成负值，不能再使用上述最大值论证。网格减半而时间步不变，也会让 $r$ 增大到原来的四倍。要维持同一个 $r$ ， $h$ 减半时应把 $\Delta t$ 缩为四分之一。

用一个误差模态检查放大因子

在无限均匀网格或周期网格上，考虑误差模态

e_j^n=G^n e^{\mathrm i j\theta}.

这里 $\theta$ 是相邻节点之间的相位增量。代入显式格式并消去公共因子，得到

G_{\rm FE}(\theta) =1-2r+r(e^{-\mathrm i\theta}+e^{\mathrm i\theta}) =1-4r\sin^2(\theta/2).

每推进一步，该模态就乘一次 $G$ 。在这里的常系数、正交模态分析中，要求所有 $\theta$ 都有 $|G|\le1$ ，便得到不放大误差的条件：

-1\le1-4r,\qquad 0\le r\le\frac12.

最危险的是 $\theta=\pi$ 的交错模态：相邻点正负交替。 $r>1/2$ 时，其放大因子小于 $-1$ ，于是它既翻转符号，又逐步增长。

有限区间的端点会改变可用波数

零 Dirichlet 端点下，不能把处处交错的 $(−1)^j$ 当作满足边界的误差向量。真正的离散模态是

z_j^{(m)}=\sin\frac{m\pi j}{N},\qquad m=1,\ldots,N-1.

利用正弦加法公式，

Tz^{(m)}=4\sin^2\frac{m\pi}{2N}\,z^{(m)}.

这些向量是实对称矩阵 $T$ 的一组正交基；不同的特征值保证正交， $N-1$ 个非零独立向量正好覆盖全部内部空间。因此显式推进矩阵 $I-rT$ 对每个模态的因子为

G_m=1-4r\sin^2\frac{m\pi}{2N}.

固定 $N$ 时，所有模态均不放大的准确条件是

0\le r\le\frac1{2\cos^2(\pi/(2N))}.

右端略大于 $1/2$ ，但随 $N\to\infty$ 趋向 $1/2$ 。因此 $r\le1/2$ 是适用于所有这些网格的充分条件，而且还保证非负权重与最大值性质。不要把“某个固定网格上的平方范数不增长”和“更新是凸组合”当成完全相同的要求。

$r=1/2$ 也有值得留意的现象：周期网格的交错模态满足 $G=-1$ ，幅度完全不衰减；连续热方程中的高频却应快速衰减。稳定只保证不增长，没有保证把衰减速度算准。有限 Dirichlet 网格没有恰好这个 $\theta=\pi$ 模态，但最高模态在细网格上同样可能缓慢翻转衰减。

无限网格或偶数节点周期网格上，同一个交错模态在显式与隐式算法中的一步因子。

隐式 Euler：新时间层要一起求

把空间差分放在第 $n+1$ 层，得到

\frac{U_j^{n+1}-U_j^n}{\Delta t} =\kappa\frac{U_{j+1}^{n+1}-2U_j^{n+1}+U_{j-1}^{n+1}}{h^2},

即

-rU_{j-1}^{n+1}+(1+2r)U_j^{n+1}-rU_{j+1}^{n+1}=U_j^n.

相邻新值同时未知，不能像显式格式那样直接扫一遍代入。每一步要解

(I+rT)U^{n+1}=b^n.

对零端点， $b^n$ 就是旧的内部向量。如果两端在新时刻的值为 $a(t^{n+1})$ 和 $b(t^{n+1})$ ，右端首行要加，末行加。这次使用的是的边界，因为被移项的邻居来自新层。

矩阵可逆也有直接理由：

z^{\mathsf T}(I+rT)z =\sum_{j=1}^{N-1}z_j^2+r\sum_{j=0}^{N-1}(z_{j+1}-z_j)^2>0

对所有非零 $z$ 、 $r\ge0$ 成立。因此每一步的内部值唯一确定。

三对角系统怎样解

这里不必每次对一个稠密大矩阵做通用消元。记三对角方程为

\ell_jX_{j-1}+d_jX_j+c_jX_{j+1}=b_j.

第一行没有左邻未知量。由上一行消去当前行的 $X_{j-1}$ 时，乘数是

m_j=\ell_j/\widetilde d_{j-1},

并更新

\widetilde d_j=d_j-m_jc_{j-1},\qquad \widetilde b_j=b_j-m_j\widetilde b_{j-1}.

消到最后一行后，从 $X_{N-1}=\widetilde b_{N-1}/\widetilde d_{N-1}$ 开始向前回代：

X_j=\frac{\widetilde b_j-c_jX_{j+1}}{\widetilde d_j}.

这就是三对角追赶法，每一步只需与内部节点数成正比的运算。在当前隐式热系统中， $d_j=1+2r$ 、 $\ell_j=c_j=-r$ ，消元主元满足

\widetilde d_1=1+2r,\qquad \widetilde d_j=1+2r-\frac{r^2}{\widetilde d_{j-1}}.

$r>0$ 时可归纳得到 $\widetilde d_j>1+r>0$ ：若上一主元大于，减掉的就小于。因此这里不会遇到零主元。则是单位矩阵，直接保留旧值。

无条件稳定的两个检查

代入同一个 Fourier 模态，空间新层多出一个 $G$ ，所以

G=1-4rG\sin^2(\theta/2),\qquad \boxed{G_{\rm BE}=\frac1{1+4r\sin^2(\theta/2)}.}

只要 $r\ge0$ ，就有 $0<G_{\rm BE}\le1$ 。在有限区间上，只需换成前面允许的离散波数，结论仍成立。

最大范数也能直接检查。若 $(I+rT)z=b$ 、端点补零，取 $|z_j|$ 最大的那一行：

(1+2r)|z_j| \le|b_j|+r|z_{j-1}|+r|z_{j+1}| \le\|b\|_\infty+2r\|z\|_\infty.

移项得 $\|z\|_\infty\le\|b\|_\infty$ 。因此这个隐式求解操作不会放大差向量的最大范数，不需要 $r\le1/2$ 。

但稳定不代表准确。对一个连续模态 $\sin(kx)$ ，真实的一步衰减因子是 $e^{-\kappa k^2\Delta t}$ ；数值因子使用的却是离散波数和 Euler 时间近似。时间步很大时，两者可能相差明显。隐式 Euler 还会有自己的耗散误差，不能只凭曲线没有爆炸就认定它可信。

实验：同一时刻比较显式、隐式和精确解

实验中两种算法使用同一个网格和同一个时间步。隐式曲线通过真正求解三对角方程组得到；精确曲线只用来验算，没有拿来替代算法。

从稳定步长开始推进，比较三条曲线。然后使用不稳定预设，留意很小的最高模态如何在显式结果中增长。大步长预设则用来检查另一面：隐式结果保持有界，却可能明显偏离精确解。改为非零端点时，检查边界是否在每一层都保持正确，以及隐式方程组残差是否仍然很小。

从截断误差到固定时刻的总误差

现在把精确取样 $v_j^n=u(x_j,t^n)$ 代入显式格式，通常会留下一个残差：

\tau_j^n= \frac{v_j^{n+1}-v_j^n}{\Delta t}-\kappa D_{xx}v_j^n.

在固定时空范围内，若所需的时间、空间导数连续有界，Taylor 展开给出

\tau_j^n= \frac{\Delta t}{2}u_{tt}(x_j,t^n) -\frac{\kappa h^2}{12}u_{xxxx}(x_j,t^n) +O(\Delta t^2+h^4).

因此存在与小步长无关的常数 $C$ ，使 $\|\tau^n\|_\infty\le C(\Delta t+h^2)$ 。这里把除过的方程残差称为截断误差；。把这两个量混淆，会在累积误差时多算或少算一个时间步。

把相同边界产生的已知向量记为 $q^n$ 。精确取样满足

v^{n+1}=Av^n+q^n+\Delta t\,\tau^n,

数值解满足 $U^{n+1}=AU^n+q^n$ ，其中 $A=I-rT$ 。向量在首行含、末行含，零边界时为零；相同边界的这些已知项在相减后消失。于是满足

e^{n+1}=Ae^n-\Delta t\,\tau^n.

当 $r\le1/2$ 时，前面的最大范数收缩性给出

\|e^{n+1}\|_\infty \le\|e^n\|_\infty+\Delta t\,C(\Delta t+h^2).

连续使用 $n$ 次，并限制在同一个物理终点 $n\Delta t\le T$ ，得到

\boxed{\|e^n\|_\infty \le\|e^0\|_\infty+TC(\Delta t+h^2).}

初值精确取样时 $e^0=0$ 。因此在光滑性、正确边界与稳定步长等条件下，随着 $\Delta t,h\to0$ ，数值解确实趋向精确解，典型全局误差是 $O(\Delta t+h^2)$ 。这个结论是由一致性和稳定性共同得到的，不能只从 Taylor 展开的一行局部公式直接宣布。

稳定的误差递推把每步缺陷累积到同一个物理终点；一步缺陷含额外的时间步因子。

对隐式格式，也可以在新时刻定义

\tau_j^{n+1}=\frac{v_j^{n+1}-v_j^n}{\Delta t}-\kappa D_{xx}v_j^{n+1}.

展开后时间首误差是 $-\Delta t\,u_{tt}(x_j,t^{n+1})/2$ ，空间项仍为 $- κ h^{2} u_{x}$ ，所以同样是。差方程为

(I+rT)e^{n+1}=e^n-\Delta t\,\tau^{n+1}.

使用刚证明的逆矩阵最大范数界，又得到同样形式的累积估计。这次不受 $r\le1/2$ 限制，但要让误差趋零，时间步仍必须趋零。

网格加密必须比较同一道题

比较 $h$ 与 $h/2$ 时，应固定方程、初值、边界和观察终点 $T$ 。若显式格式保持 $r$ 不变，时间步同时缩为四分之一，步数则约增为四倍，时间和空间的典型误差都会降到原来的四分之一。若只保持步数不变，较细网格实际走过的物理时间变短，就不是同一时刻的收敛比较。

误差比例也有适用边界。二次多项式的稳态例子在节点上本来就精确，剩下的浮点误差可能让两个极小数的比值毫无意义。更合适的检验是 $-u''=12x^2$ 、零端点，此时精确解为

u(x)=x-x^4.

因为 $D_{xx}(x^4)=12x^2+2h^2$ ，代入可得

-D_{xx}u=12x^2+2h^2.

离散解 $U$ 满足 $-D_{xx}U=12x^2$ ，所以误差 $e=U-u$ 满足、零端点。二次函数的差分精确，因而

\boxed{e_j=-h^2x_j(1-x_j).}

偶数 $N$ 的网格包含中点，最大节点误差恰为 $h^2/4$ 。这就给出一个可以独立验算的二阶收敛实验，同时说明残差接近零与离散误差不为零完全可以并存。

下面哪项最准确地描述稳定性与准确性的关系？

三对角方程组的残差接近零，就说明数值节点与连续方程精确解的节点值完全相同。

练习

1. 步长限制

对 $u_t=0.5u_{xx}$ ，取 $h=0.2$ 。采用 $r\le1/2$ 的显式步长条件时，中哪些满足要求？

2. 稳定而不衰减的高频

在含偶数个节点的周期网格上， $r=1/2$ ，初始误差为 $e_j^0=(-1)^j$ 。写出 $e_j^n$ 。这说明“稳定”和“正确模拟扩散”有什么区别？

该模态 $\theta=\pi$ ，所以 $G=1-4r=-1$ ，得到 $e_j^n=(-1)^{n+j}$ 。误差幅度不增长，因此满足这里的不放大要求，却也不衰减；连续热方程会迅速削弱高频。故稳定并不保证衰减过程算得准确。周期与偶数节点条件使这个交错模态确实符合边界，不能直接把它当成有限零 Dirichlet 网格的模态。

3. 二维的步长条件

对二维热方程，用五点 Laplacian 与显式 Euler。若两个方向的网格步长分别为 $h_x,h_y$ ，令 $r_x=\kappa\Delta t/h_x^2$ 、。求适用于所有 Fourier 波数的不放大条件，并给出时的结果。

4. 改用隐式后还需要检查什么

空间网格很细，显式稳定性要求的时间步小到无法接受。提出一个本章学过的替代办法，并说明它解决了什么问题、还留下什么问题。

5. 手算差分误差

令 $v(x)=x^4$ 。在 $x=1$ 处求 $D_xv$ 、与精确导数的差。若把减半，两项误差分别怎样变化？

6. 换源项与非零端点

求 $-u''=2$ 、 $u(0)=1,u(1)=3$ 的 $h=1/4$ 差分系统，并解出三个内部值。

7. 一步显式与一步隐式

零端点下， $N=4,r=1/2$ ，旧内部向量为 $(0,1,0)^{\mathsf T}$ 。分别计算一步显式和隐式结果。

8. 比较同一终点

取 $\kappa=1,r=0.4$ ，希望模拟到 $T=0.1$ 。 $h=0.1$ 时需要多少步？改为 $h=0.05$ 并保持不变后，需要多少步？

9. 用四次解验证收敛阶

对 $-u''=12x^2$ 、零端点，计算 $N=8$ 与 $N=16$ 时的最大节点误差。为什么不能改用正文的常数源二次解，再把几乎为零的误差相除来估计阶数？

10. 新时刻的边界进入哪一层

隐式热格式取 $N=4,\kappa=1,\Delta t=1/32$ ，初始内部值全为零，左边界为零，右边界为 $u(1,t)=t$ 。写出第一步的矩阵与右端向量。

$h=1/4$ ，所以 $r=(1/32)/(1/16)=1/2$ 。矩阵对角为 $2$ 、副对角为 $-1/2$ 。第一步的右边界是，因此右端为；不能仍代入旧时刻的零边界。这个题的初始角点满足值的兼容，但初始二阶空间导数与边界时间导数不兼容，不能直接用覆盖初始角点的高阶光滑性来宣称前述误差界。

11. 负的中心权重会发生什么

显式更新取 $r=0.6$ ，某内部点的左邻点、当前点、右邻点的旧值依次为 $0,1,0$ 。求新值，并说明它破坏了哪条性质。

x

x

x

/

12

-\kappa h^2u_{xxxx}/12

G=1-4r_x\sin^2(\theta_x/2)-4r_y\sin^2(\theta_y/2).

D_{x x}

v

D_{xx}v