Fourier 方法分析差分格式稳定性

问题设定与格式回顾

考虑抛物方程

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

以向前差分格式

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

为例

设精确解为 $u_j^n$ ，近似解为 $v_j^n$ ，误差 $e_j^n = u_j^n - v_j^n$ ，误差满足

$e_j^{n + 1}=re_{j - 1}^n+(1 - 2r)e_j^n+re_{j + 1}^n$

分离变量法求解原方程

对于不含源项的初边值问题

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

使用分离变量法设

$u(x,t)=X(x)T(t)$

代入方程可得

$\frac{T'(t)}{aT(t)}=\frac{X''(x)}{X(x)}=-\lambda$

$\lambda$ 为分离常数

求解 $X(x)$ ：结合边界条件 $X(0)=X(l)=0$ ，可得

$X_k(x)=b_k\sin\frac{k\pi}{l}x,\quad k = 1,2,\cdots$

求解 $T(t)$ ：

$T_k(t)=a_ke^{-(\frac{k\pi}{l})^2at}$

则原方程解

$\begin{align} u(x,t)&=\sum_{k = 1}^{\infty}a_kb_k e^{-(\frac{k\pi}{l})^2at}\sin(\frac{k\pi}{l})x\\ &=\sum_{k = 1}^{\infty}c_ke^{-(\frac{k\pi}{l})^2at}e^{i\frac{k\pi}{l}x} \end{align}$

其中 $c_k$ 由初始条件确定

$\varphi(x)=\sum_{k = 1}^{\infty}c_ke^{i\frac{k\pi}{l}x}$

已知

$\begin{align} u(x,t_{n + 1})&=\sum c_ke^{-a(\frac{k\pi}{l})^2(t_n+\tau)}e^{i\frac{k\pi}{l}x}\\ u(x,t_n)&=\sum c_ke^{-a(\frac{k\pi}{l})^2t_n}e^{i\frac{k\pi}{l}x} \end{align}$

对应级数每一项的 “振幅” 和增长因子为

$e^{\displaystyle -a(\frac{k\pi}{l})^2\tau}$

将差分格式解表示为 Fourier 级数

把 $u_j^n$ 看成是区间上的分片常数函数 $u^n(x)$ ，即

$u^n(x)=u_j^n \quad (x_{j - \frac{1}{2}}<x<x_{j + \frac{1}{2}})$

且

$u^n(x)=\sum_{k = -\infty}^{+\infty}v_k^ne^{i\frac{k\pi}{l}x}\quad (Fourier)$ $u^{n + 1}(x)=\sum_{k = -\infty}^{+\infty}v_k^{n + 1}e^{i\frac{k\pi}{l}x}$

由向前差分格式

$u^{n + 1}(x)=ru^n(x - h)+(1 - 2r)u^n(x)+ru^n(x + h)$

将 Fourier 级数代入可得：

$\begin{align} \sum_{k = -\infty}^{+\infty}v_k^{n + 1}e^{i\frac{k\pi}{l}x}=&r\sum_{k = -\infty}^{+\infty}v_k^ne^{i\frac{k\pi}{l}(x - h)}\\ &+(1 - 2r)\sum_{k = -\infty}^{+\infty}v_k^ne^{i\frac{k\pi}{l}x}\\ &+r\sum_{k = -\infty}^{+\infty}v_k^ne^{i\frac{k\pi}{l}(x + h)} \end{align}$

整理得到

$\begin{align} v_k^{n + 1}&=v_k^n(r e^{-i\frac{k\pi}{l}h}+(1 - 2r)+r e^{i\frac{k\pi}{l}h}) \\ &= G(\lambda h)v_k^n \end{align}$

其中 $\displaystyle \lambda=\frac{k\pi}{l}$ ，$G(\lambda h)$ 为增长因子。进而有

$v_k^n = G^n(\lambda h)v_k^0$

推导增长因子并分析稳定性

计算增长因子

$\begin{align} G(\lambda h) &=1 - 2r + 2r\cos\frac{k\pi}{l}x \\ &=1 - 4r\sin^2\frac{\lambda h}{2} \end{align}$

进一步近似

$\begin{align} G(\lambda h)&\approx1 - 4r(\frac{\lambda\pi}{2l}\cdot h)^2 \\ &=1 - a\tau(\frac{\lambda\pi}{h})^2\\ &\approx e^{\displaystyle-a\tau(\frac{\lambda\pi}{h})^2} \end{align}$

对一般的抛物方程差分格式，增长因子 $G$ 还可能与 $\tau$ 有关，记为 $G(\alpha h,\tau)$

例如带低阶项抛物方程差分格式

对于带低阶项的抛物方程

$u_t = au_{xx}+bu$

其差分格式为

$\frac{u_j^{n + 1}-u_j^n}{\tau}=a\cdot\frac{1}{h^2}\cdot\delta_x^2u_j^n + bu_j^n$

其中 $\delta_x^2$ 是二阶中心差分算子

推导出增长因子

$G = 1 - 4r\sin^2\frac{\lambda h}{2}+b\tau$

这里 $\displaystyle r=\frac{a\tau}{h^2}$ 。通过分析 $|G|\leq1$ 的条件来确定格式的稳定性，例如当 $b$ 和 $r$ 满足一定关系时，可保证增长因子的模不大于 1，从而使格式稳定

一般二层差分格式稳定性判定

对于线性常系数一维抛物方程（且 $f = 0$ ）的二层差分格式

$\sum_{m\in N}a_mu_{j + m}^{n + 1}=\sum_{m\in N_0}b_mu_{j + m}^n \quad (j = 0,1,\cdots)$

格式稳定的充要条件是存在 $\tau_0$ 和 $k$ ，使得当 $0 < n\leq\frac{T}{\tau}$ ，$0 < \tau\leq\tau_0$ 时

$\|U^n\|\leq k\|U^0\|$

这等价于在 $L^2$ 范数下

$\|u^n(x)\|_{L^2}^2\leq k^2\|u^0(x)\|_{L^2}^2$

将 $u^n(x)$ 用 Fourier 级数表示为

$u^n(x)=\sum_{k = -\infty}^{+\infty}v_k^ne^{i\frac{k\pi}{l}x}$

进而

$k_l\sum_{k = -\infty}^{+\infty}|v_k^n|^2\leq k^2k_l \sum_{k=-\infty}^{+\infty}|v_k^0|^2$ $\sum_{k = -\infty}^{+\infty}|G^n(\alpha h,\tau)v_k^0|^2\leq k^2 \sum_{k=-\infty}^{+\infty}|v_k^0|^2$

取特殊的 $v_k^0$ （只有一个 $v_k^0 = 1$ ，其他为 0 ），可得到稳定性等价于

$|G^n(\lambda h,\tau)|\leq k$

即增长因子的幂次一致有界。反之，若 $G^n(\lambda h,\tau)$ 一致有界，也能推出

$\|u^n(x)\|_{L^2}^2\leq k^2\|u^0(x)\|_{L^2}^2$

从而确定差分格式稳定

命题3.1

差分格式稳定等价于 $G^n(\lambda h,\tau)$ 一致有界，也等价于

$|G(\lambda h,\tau)|\leq1 + M\tau$

称为 Von Neumann 条件

增长因子的计算与稳定性判定

将

$u_j^n = v^n e^{i\lambda jh}$

代入向前差分格式

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

经过推导可得增长因子

$\begin{align} G(\lambda h) &=1 - 2r + 2r\cos\lambda h \\ &= 1 - 4r\sin^2\frac{\lambda h}{2} \end{align}$

当 $G(\lambda h)$ 满足 Von Neumann 条件，即 $|G(\lambda h)|\leq1$ 时，可推出 $r\leq\frac{1}{2}$

所以当 $r\leq\frac{1}{2}$ 时，向前差分格式稳定

多层差分格式增长因子的计算

二层差分格式：一般二层差分格式，同样用 $u_j^n = v^n e^{i\lambda jh}$ 代入格式，通过化简计算得到增长因子 $G$ ，进而根据 Von Neumann 条件判断稳定性

三层差分格式（以 Richardson 格式为例）：Richardson 格式

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

令

$w_j^{n + 1}=u_j^n$

改写成方程组

$\begin{cases}u_j^{n + 1}=2r(u_{j - 1}^n - 2u_j^n + u_{j + 1}^n)+w_j^n\\w_j^{n + 1}=u_j^n\end{cases}$

设

$\begin{align} u_j^n &= v_1^n e^{i\lambda jh}\\ w_j^n &= v_2^n e^{i\lambda jh} \end{align}$

代入方程组

化简得到

$\begin{cases}v_1^{n + 1}=4r(\cos\lambda h - 1)v_1^n + v_2^n\\v_2^{n + 1}=v_1^n\end{cases}$

写成矩阵形式

$\begin{bmatrix}v_1^{n + 1}\\v_2^{n + 1}\end{bmatrix}=\begin{bmatrix}4r(\cos\lambda h - 1)&1\\1&0\end{bmatrix}\begin{bmatrix}v_1^n\\v_2^n\end{bmatrix}$

此时增长矩阵

$G(\lambda h)=\begin{bmatrix}4r(\cos\lambda h - 1)&1\\1&0\end{bmatrix}$

通过求该矩阵的特征值来分析稳定性，其特征方程

$\lambda^2 - 4r(\cos\lambda h - 1)\lambda - 1 = 0$

特征值会使谱半径大于 $1 + O(\tau)$ ，不满足 Von Neumann 条件，所以 Richardson 格式不稳定

一般的 s + 1 层差分格式稳定性分析

对于一般的 $s + 1$ 层差分格式，其稳定性分析将归结于考虑如下差分方程组

$\sum_{m\in N_1}A_m\overline{U}_{j + m}^{n + 1}=\sum_{m\in N_0}B_m\overline{U}_{j + m}^n$

设

$\overline{U}^n = V^n e^{i\lambda jh}$

代入消去 $e^{i\lambda jh}$ ，得到 $s$ 阶矩阵 $G(\lambda h,\tau)$

命题 3.2：差分格式稳定等价于

$\{G^n(\lambda h,\tau):0 < n\leq\frac{T}{\tau},0 < \tau\leq\tau_0\}$

一致有界，这是从稳定性定义出发，强调增长矩阵的幂次在一定时间步长和步长范围内不会使误差无限增长

命题 3.3：${G^n(\lambda h,\tau)}$ 一致有界可推出 $\rho(G)\leq1 + O(\tau)$ （Von Neumann 条件）。当 $G(\lambda h,\tau)$ 是正规矩阵时，两者等价