微分方程数值解 第三章（2）

稳定性与收敛性

ODE - IVP 的稳定性回顾

Euler 法：误差 $e_n$ 满足 $\displaystyle |e_n|\leq|e_{n - 1}|+L\tau|e_{n - 1}|=(1 + L\tau)|e_{n - 1}|$ ，递推可得

$$|e_n|\leq(1 + L\tau)^n|e_0|\leq e^{LT}|e_0|$$

这里 $L$ 是 Lipschitz 常数，$\tau$ 是步长，$T$ 是总时间，说明在一定条件下 Euler 法的误差可控

多步法：设 $\displaystyle E_n=(e_{n + k - 1},\cdots,e_n)^T$ ，则 $E_{n + 1}=CE_n + B_n$ ，要求 ${C^n}$ 一致有界，这样多步法的误差增长才能得到控制，保证方法的稳定性

椭圆方程 BVP 差分法稳定性回顾

对于椭圆方程边值问题

$$\begin{cases}L_hv_i = f_i\\v_0 = v_N = 0\end{cases}$$

若满足 $|v_h|\leq M|f_h|_R$ ，则称该差分法按右端稳定，即解的范数能被右端项范数控制，保证了方法在求解椭圆方程时的稳定性

抛物方程初边值问题的稳定性类型

按初值稳定：关注初始条件误差对解的影响，若初始条件有小扰动，解不会产生过大偏差，则格式按初值稳定

按边值稳定：考虑边界条件误差对解的作用，当边界条件存在误差时，解仍能保持在合理范围，格式按边值稳定

抛物方程差分格式的统一形式

向前差分格式

$$u_j^{n + 1}=ru_{j - 1}^n+(1 - 2r)u_j^n+ru_{j + 1}^n+\tau f_j$$

可写成矩阵形式

$$\begin{align} U^{n + 1} &=\begin{bmatrix}1 - 2r&r& & \\r&1 - 2r&r& \\& \ddots&\ddots&\ddots \\& &r&1 - 2r\end{bmatrix}U^n+\tau F\\ &=BU^n+\tau F \end{align}$$

其中 $B=(1 - 2r)I + rS$ ，$S$ 是特定的 $(J - 1)\times(J - 1)$ 矩阵

向后差分格式

$$-ru_{j - 1}^{n + 1}+(1 + 2r)u_j^{n + 1}-ru_{j + 1}^{n + 1}=u_j^n+\tau f_j$$

写成矩阵形式为

$$\begin{bmatrix}1 + 2r&-r& & \\-r&1 + 2r&-r& \\& \ddots&\ddots&\ddots \\& &-r&1 + 2r\end{bmatrix}U^{n + 1}=U^n+\tau F$$

即

$$AU^{n + 1}=U^n+\tau F$$

Crank - Nicolson 格式

$$-\frac{r}{2}u_{j - 1}^{n + 1}+(1 + r)u_j^{n + 1}-\frac{r}{2}u_{j + 1}^{n + 1}=\frac{r}{2}u_{j - 1}^n+(1 - r)u_j^n+\frac{r}{2}u_{j + 1}^n+\tau f_j$$

矩阵形式为

$$\begin{bmatrix}1 + r&-\frac{r}{2}& & \\-\frac{r}{2}&1 + r&-\frac{r}{2}& \\& \ddots&\ddots&\ddots \\& &-\frac{r}{2}&1 + r\end{bmatrix}U^{n + 1}=\begin{bmatrix}1 - r&\frac{r}{2}& & \\\frac{r}{2}&1 - r&\frac{r}{2}& \\& \ddots&\ddots&\ddots \\& &\frac{r}{2}&1 - r\end{bmatrix}U^n+\tau F$$

即

$$AU^{n + 1}=BU^n+\tau F$$

Richardson 格式

$$u_j^{n + 1}=2r(u_{j - 1}^n - 2u_j^n + u_{j + 1}^n)+u_j^{n - 1}+2\tau f_j$$

令 $W^n=(U^{n + 1},U^n)^T$ ，写成矩阵形式为

$$\begin{bmatrix}U^{n + 1}\\U^n\end{bmatrix}=\begin{bmatrix}D&I\\ I&0\end{bmatrix}\begin{bmatrix}U^n\\U^{n - 1}\end{bmatrix}+2\tau\begin{bmatrix}F\\0\end{bmatrix}$$

即

$$W^{n + 1}=CW^n+\tau\begin{bmatrix}2F\\0\end{bmatrix}$$

统一形式

抛物方程差分格式统一为

$$AU^{n + 1}=BU^n+\tau F$$

变形得到

$$\begin{align} U^{n + 1} &=A^{-1}BU^n+\tau A^{-1}F \\ &= CU^n+\tau A^{-1}F \end{align}$$

其中 $C = A^{-1}B$ 为增长矩阵，其与空间步长 $h$ 和时间步长 $\tau$ 有关 。增长矩阵在稳定性分析中起关键作用，它刻画了每一步数值解的变化情况

按初值稳定性定义推导

为研究按初值稳定性，设 $F = 0$ （消除非齐次项影响，专注初始值对解的影响 ）。假设有两组解 $U^n$ 和 $V^n$ ，分别满足

$$\begin{align} U^{n + 1} &=CU^n+\tau A^{-1}F\\ V^{n + 1} &=CV^n+\tau A^{-1}F \end{align}$$

令 $E^n = U^n - V^n$ ，则

$$E^{n + 1}=CE^n$$

当 $F = 0$ 时

$$U^{n + 1}=CU^n$$

递推可得

$$U^{n + 1}=C^{n + 1}U^0$$

若存在 $\tau_0>0$ 以及与 $h$、$\tau$ 无关的常数 $k>0$ ，使得当 $0 < \tau\leq\tau_0$ 时，对于任意初始向量 $U^0\in\mathbb{R}^{J - 1}$ （$J$ 为空间节点数 ），在 $0 < n\tau\leq T$ （$T$ 为总时间 ）的范围内，都有

$$\|U^{n + 1}\|=\|C^{n + 1}U^0\|\leq k\|U^0\|$$

就称该差分格式按初值稳定 。这等价于存在 $\tau_0>0$ 和 $k>0$ ，使得对于 $0 < n\tau\leq T$ 且 $0 < \tau\leq\tau_0$ ，有

$$\|C^n\|\leq k$$

这里

$$\|C^n\|=\sup_{\|v\| = 1}\|C^nv\|$$

直观来说，就是初始值的误差在计算过程中不会无限制放大，保证了数值解的稳定性

按右端稳定性定义推导

对于差分格式 $U^{n + 1}=CU^n+\tau A^{-1}F$ ，为研究按右端稳定性，可设 $U^0 = 0$ （消除初始值影响，专注右端项 $F$ 对解的作用 ）。假设有两组解 $U^n$ 和 $V^n$ ，分别满足

$$\begin{align} U^{n + 1} &=CU^n+\tau A^{-1}F_1\\ V^{n + 1} &=CV^n+\tau A^{-1}F_2 \end{align}$$

令

$$E^n = U^n - V^n$$

则

$$E^{n + 1}=CE^n+\tau A^{-1}(F_1 - F_2)$$

且 $E^0 = 0$ ，所以可专注于 $U^n$ 受 $F$ 影响的情况

若存在 $\tau_0>0$ 以及与 $h$、$\tau$ 无关的常数 $k>0$ ，使得当 $0 < \tau\leq\tau_0$ 时，对于任意右端项向量 $F\in\mathbb{R}^{J - 1}$ （$J$ 为空间节点数 ），在 $0 < n\tau\leq T$ （$T$ 为总时间 ）的范围内，都有

$$\|U^{n + 1}\|\leq k\|F\|$$

则称该格式按右端稳定 。这表明右端项的变化不会导致数值解无限制地增大，保证了解对右端项的稳定性

按初值稳定与按右端稳定的关系证明

前提条件：已知 $|A^{-1}|\leq k’$

推导过程：设 $U^0 = 0$ ，则

$$\begin{align} U^1 &=\tau A^{-1}F\\ U^2 &= CU^1+\tau A^{-1}F=(C + I)\tau A^{-1}F\\ U^3 &= CU^2+\tau A^{-1}F=(C^2 + C + I)\tau A^{-1}F\\ & \quad \vdots\\ U^{n + 1} &=(C^n + C^{n - 1}+\cdots + C + I)\tau A^{-1}F \end{align}$$

若差分格式按初值稳定，即存在常数 $k’’$ ，使得

$$\|C^n\|\leq k''$$

那么

$$\|U^{n + 1}\|\leq(n + 1)k''\cdot\tau\cdot k'\cdot\|F\|$$

当 $n\tau\leq T$ 时，可找到合适的 $k$ ，使得

$$\|U^{n + 1}\|\leq k\|F\|$$

即差分格式按右端稳定 。所以在给定条件下，按初值稳定可推出按右端稳定，也解释了为何稳定性常指按初值稳定

稳定性分析总体思路

要检验矩阵族 $\displaystyle {C^n: 0 < \tau\leq\tau_0, 0 < n\leq\frac{T}{\tau}}$ 的一致有界性，常用方法有矩阵法、Fourier 方法、代数准则法

矩阵法相关命题

命题 2.1（必要条件）：设增长矩阵 $C$ 的谱半径为 $\rho(C)$ ，若差分格式稳定，则存在 $M > 0$ ，使得 $\rho(C)\leq1 + M\tau$

证明思路：因为格式稳定，存在 $\tau_0$ 和 $k\geq1$ ，当 $0 < \tau\leq\tau_0$ ，$\displaystyle 0 < n\leq\frac{T}{\tau}$ 时，$|C^n|\leq k$ 。又因为

$$\rho(C)=\rho(C^n)^{\frac{1}{n}}\leq\|C^n\|^{\frac{1}{n}}\leq k^{\frac{1}{n}}$$

对 $k^{\frac{1}{n}}$ 进行放缩，如

$$k^{\frac{1}{n}}\leq k^{\frac{1}{\frac{T}{\tau}+1}}\leq k^{\frac{\tau}{T+\tau}} = e^{\frac{\tau}{T + \tau}\ln k}$$

进一步放缩可得

$$e^{\frac{\tau}{T + \tau}\ln k}\leq e^{\frac{\ln k}{T - \tau_0}\cdot\tau}\leq1 + M\tau$$

命题 2.2（充分条件）：设 $C^HC = CC^H$ （$C$ 为正规矩阵），若 $\rho(C)\leq1 + M\tau$ ，则差分格式稳定

证明思路：此时 $|C|_2 = \rho(C)$ ，则

$$\|C^n\|\leq\|C\|^n=\rho(C)^n\leq(1 + M\tau)^n\leq(e^{M\tau})^n = e^{Mn\tau}$$

当 $n\tau\leq T$ 时

$$e^{Mn\tau}\leq e^{TM}=k$$

说明矩阵族有界，格式稳定

四种格式稳定性分析

向前差分格式

向前差分格式中

$$C=(1 - 2r)I + rS$$

$S$ 为特定矩阵且 $C^HC = CC^H$ ，$S$ 的特征值

$$\lambda_j^S = 2\cos j\pi h \quad (j = 1,\cdots,J - 1,h=\frac{L}{J} )$$

特征向量

$$V_j^S = (\sin j\pi h,\sin 2j\pi h,\cdots,\sin (J - 1)j\pi h)^T$$

由此可得 $C$ 的特征值

$$\begin{align} \lambda_j^C &= 1 - 2r + r\cdot2\cos j\pi h \\ &= 1 - 4r\sin^2\frac{j\pi h}{2} \end{align}$$

要使 $|\lambda_j^C|\leq1 + M\tau$ ，分析

$$|1 - 4r\sin^2\frac{j\pi h}{2}|\leq1 + M\tau$$

右边 $1 + M\tau\geq1$ 恒成立，只需考虑左边。当 $4r\leq2$ ，即 $r\leq\frac{1}{2}$ 时

$$|1 - 4r\sin^2\frac{j\pi h}{2}|\leq1$$

能保证，此时向前差分格式稳定；当 $r > \frac{1}{2}$ 时，$|1 - 4r\sin^2\frac{j\pi h}{2}|$ 可能大于 $1$ ，格式不稳定

向后差分格式

向后差分格式增长矩阵

$$C = [(1 + 2r)I - rS]^{-1}=A^{-1}$$

其特征值

$$\begin{align} \lambda_j^C &=(1 + 2r - r\lambda_j^S)^{-1}\\ &=(1 + 2r - r\cdot2\cos j\pi h)^{-1}\\ &=(1 + 4r\sin^2\frac{j\pi h}{2})^{-1} \end{align}$$

因为 $(1 + 4r\sin^2\frac{j\pi h}{2})^{-1}\leq1$ 恒成立，所以

$$\rho(C)\leq1\leq1 + M\tau$$

即向后差分格式恒稳定

Crank - Nicolson 格式

增长矩阵

$$\begin{align} C &= A^{-1}B\\ &=[(1 + r)I - \frac{r}{2}S]^{-1}[(1 - r)I + \frac{r}{2}S] \end{align}$$

其中

$$(1 + r)I - \frac{r}{2}S = V^T\begin{bmatrix}\alpha_1& & \\ & \ddots& \\ & &\alpha_{J - 1}\end{bmatrix}V$$ $$(1 - r)I + \frac{r}{2}S = V^T\begin{bmatrix}\beta_1& & \\ & \ddots& \\ & &\beta_{J - 1}\end{bmatrix}V$$ $$\alpha_j = 1 + r - \frac{r}{2}\cdot2\cos j\pi h = 1 + 2r\sin^2\frac{j\pi h}{2}$$ $$\beta_j = 1 - r + \frac{r}{2}\cdot2\cos j\pi h = 1 - 2r\sin^2\frac{j\pi h}{2}$$

则

$$\begin{align} C &= V^T\begin{bmatrix}\alpha_1^{-1}& & \\ & \ddots& \\ & &\alpha_{J - 1}^{-1}\end{bmatrix}V\cdot V^T\begin{bmatrix}\beta_1& & \\ & \ddots& \\ & &\beta_{J - 1}\end{bmatrix}V \\ &= V^T\begin{bmatrix}\alpha_1^{-1}\beta_1& & \\ & \ddots& \\ & &\alpha_{J - 1}^{-1}\beta_{J - 1}\end{bmatrix}V \end{align}$$

特征值

$$|\lambda_j^C| = \left|\frac{\beta_j}{\alpha_j}\right|=\left|\frac{1 - 2r\sin^2\frac{j\pi h}{2}}{1 + 2r\sin^2\frac{j\pi h}{2}}\right|\leq1$$

故 Crank - Nicolson 格式恒稳定

Richardson 格式

Richardson 格式的增长矩阵

$$C=\begin{bmatrix}2r(S - 2I)&I\\I&0\end{bmatrix}$$

对于特征值问题

$$C\begin{bmatrix}w_1\\w_2\end{bmatrix}=\lambda\begin{bmatrix}w_1\\w_2\end{bmatrix}$$

得到方程组

$$\begin{cases}2r(S - 2I)w_1+w_2=\lambda w_1\\w_1=\lambda w_2\end{cases}$$

将 $w_1=\lambda w_2$ 代入第一个方程并整理，可得

$$S w_2=(2+\frac{\lambda}{2r}-\frac{1}{2\lambda r})w_2$$

已知 $S$ 的特征值

$$\mu = 2\cos j\pi h$$

则

$$\mu = 2+\frac{\lambda}{2r}-\frac{1}{2\lambda r}$$

进一步变形为关于 $\lambda$ 的二次方程

$$\lambda^2 + 8r\sin^2\frac{j\pi h}{2}-1 = 0$$

解上述二次方程，得到

$$\lambda_{1,2}=-4r\sin^2\frac{j\pi h}{2}\pm\sqrt{16r^2\sin^2\frac{j\pi h}{2}+1}$$

谱半径

$$\rho(C)=\max_j\{|\lambda_1^j|,|\lambda_2^j|\}$$

通过分析可知

$$\rho(C)>r+\sqrt{1 + r^2}>r + 1 = 1+\frac{a\tau}{h^2}$$

根据矩阵法中差分格式稳定的必要条件（若差分格式稳定，则 $\rho(C)\leq1 + M\tau$ ），这里不存在 $M$ 使得 $\rho(C)\leq1 + M\tau$ 成立

故 Richardson 格式恒不稳定

收敛性与敛速估计

考虑抛物方程初边值问题

$$\begin{cases}Lu = u_t - au_{xx}=f(x)\\u(x,0)=\varphi(x)\\u(0,t)=u(l,t)=0\end{cases}$$

对应的差分格式为

$$\begin{cases}L_hu_j^n = f_j\\u_j^0=\varphi_j\\u_0^n = u_J^n = 0\end{cases}$$

其中截断误差定义为

$$R_j^n(u)=L_hu(x_j,t_n)-[Lu]_j^n$$

当 $h,\tau\rightarrow0$ 时，若 $R_j^n(u)\rightarrow0$ ，则称差分格式具有相容性

定理 2.1

定理内容：对于线性抛物方程的差分格式，相容性 + 稳定性 $\Rightarrow$ 收敛性

证明思路：由相容性可知

$$\lim\limits_{h,\tau\rightarrow0}\|R^n(u)\| = 0$$

令

$$e_j^n = u(x_j,t_n)-u_j^n$$

为误差函数

根据截断误差的定义

$$\begin{align} R_j^n(u) &=L_hu(x_j,t_n)-[Lu](x_j,t_n)=L_hu(x_j,t_n)-f_j\\ &= L_hu(x_j,t_n)-L_hu_j^n \\ &= L_he_j^n \end{align}$$

可得

$$\begin{cases}L_he_j^n = R_j^n(u)\\e_j^0 = 0, \quad j = 1,\cdots,J - 1\end{cases}$$

写成矩阵向量形式为

$$AE^{n + 1}=BE^n+\tau R^n$$

进一步得到

$$\begin{align} E^{n + 1} &=A^{-1}BE^n+\tau A^{-1}R^n \\ &= CE^n+\tau A^{-1}R^n \end{align}$$

其中

$$E^n=(e_1^n,\cdots,e_{J - 1}^n)^T$$

因为按初值稳定可推出按右端稳定，根据按右端稳定性质有

$$\|E^n\|\leq k\|R^n\|$$

又因为相容性保证

$$\lim\limits_{h,\tau\rightarrow0}\|R^n\| = 0$$

所以

$$\lim\limits_{h,\tau\rightarrow0}\|E^n\| = 0$$

即证明了差分格式的收敛性

推论 2.2 （收敛阶）

向前格式：当 $r\leq\frac{1}{2}$ 时，向前格式稳定，截断误差为 $O(\tau + h^2)$ ，所以收敛阶为 $O(\tau + h^2)$ 。这意味着随着时间步长 $\tau$ 和空间步长 $h$ 趋于 0 ，数值解与精确解的误差以 $\tau + h^2$ 的速度趋近于 0

向后格式：向后格式恒稳定，截断误差是 $O(\tau + h^2)$ ，其收敛阶为 $O(\tau + h^2)$

Crank - Nicolson 格式：该格式恒稳定，截断误差为 $O(\tau^2 + h^2)$ ，收敛阶为 $O(\tau^2 + h^2)$ ，表明其误差减小速度比前两者在时间方向上更快，精度更高

Richardson 格式：由于恒不稳定，即便截断误差有一定阶数，但误差会随计算过程无限增长，不存在收敛性