抛物型方程有限差分法

方程类型与定解问题分类

方程类型：抛物型方程

$\frac{\partial u}{\partial t}=a\Delta u + f(x) \quad(a>0)$

其中时间变量是一阶导数，空间变量是二阶导数。当 $t\rightarrow\infty$ 时，$\displaystyle \frac{\partial u}{\partial t}\rightarrow0$ ，方程变为 $0 = a\overline{u}+f(x)$ ，即椭圆方程，说明椭圆方程是抛物方程的定态情形，抛物方程关注中间过程

定解问题分类（一维情形）

初值问题（Cauchy）：给定初始条件

$u(x,0)=\varphi(x),\quad -\infty<x<+\infty$

初边值问题：初始条件

$u(x,0)=\varphi(x),\quad 0\leq x\leq l$

边界条件有多种形式，如

$\begin{align} &① u(0,t)=u(l,t),\quad 0\leq t\leq T\\ &② u(0,t)=\alpha,\quad u_x(l,t)=\beta\\ &③ u(0,t)+\alpha_1u_x(0,t)=\beta_1,\quad u(l,t)+\alpha_2u_x(l,t)=\beta_2 \end{align}$

最简差分格式

考虑方程

$\begin{cases}\displaystyle \frac{\partial u}{\partial t}=a\frac{\partial^2u}{\partial x^2}+f(x),&0<t\leq T,x\in(0,l)\\u(x,0)=\varphi(x),&x\in[0,l]\\u(0,t)=0,u(l,t)=0\end{cases}$

区域离散化

对于区域 $\overline{G}={(x,t):0\leq x\leq l,0\leq t\leq T}$ ，在 $x$ 方向将区间 $[0, l]$ 进行 $J$ 等分，步长 $\displaystyle h = \frac{l}{J}$ ，得到节点 $x_j = jh$（$j = 0, 1, \cdots, J$ ）；在 $t$ 方向将区间 $[0, T]$ 进行 $N$ 等分，步长 $\displaystyle \tau = \frac{T}{N}$ ，得到节点 $t_n = n\tau$（$n = 0, 1, \cdots, N$ ）。用 $\displaystyle u_j^n$ 近似表示 $u(x_j, t_n)$

向前差分格式

对抛物型方程 $\displaystyle \frac{\partial u}{\partial t}=a\frac{\partial^2u}{\partial x^2}+f(x)$ ，时间方向用向前差商 $\displaystyle \frac{u_j^{n + 1}-u_j^n}{\tau}$ 近似 $\displaystyle \frac{\partial u}{\partial t}$ ，空间方向用中心差商 $\displaystyle \frac{u_{j - 1}^n - 2u_j^n + u_{j + 1}^n}{h^2}$ 近似 $\displaystyle \frac{\partial^2u}{\partial x^2}$ ，得到向前差分格式：

$\frac{u_j^{n + 1}-u_j^n}{\tau}=a\frac{u_{j - 1}^n - 2u_j^n + u_{j + 1}^n}{h^2}+f_j$

整理后为

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

令 $\displaystyle r = \frac{a\tau}{h^2}$（网比），则

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

该格式也称为显格式，因为可以根据已知的 $n$ 时刻的值直接计算出 $n + 1$ 时刻的值

初始条件

$\begin{align} u_j^0 &= \varphi_j = \varphi(x_j)\\ &f_j = f(x_j) \end{align}$

结合边界条件可依次递推求解 $u_j^n$

截断误差分析

定义微分算子

$Lu=\frac{\partial u}{\partial t}-a\frac{\partial^2u}{\partial x^2}$

差分算子

$L_h^{(1)}u_j^n=\frac{u_j^{n + 1}-u_j^n}{\tau}-a\frac{u_{j - 1}^n - 2u_j^n + u_{j + 1}^n}{h^2}$

截断误差

$R_j^n(u)=L_h^{(1)}u(x_j, t_n)-[Lu]_j^n$

通过泰勒展开：

$\frac{u(x_j, t_{n + 1})-u(x_j, t_n)}{\tau}=\frac{\partial u}{\partial t}(x_j, t_n)+\frac{\tau}{2}\frac{\partial^2u}{\partial t^2}(x_j, t_n)+O(\tau^2)$ $\frac{u(x_{j - 1}, t_n)-2u(x_j, t_n)+u(x_{j + 1}, t_n)}{h^2}=\frac{\partial^2u}{\partial x^2}(x_j, t_n)+\frac{h^2}{12}\frac{\partial^4u}{\partial x^4}(x_j, t_n)+O(h^3)$

则

$L_h^{(1)}u(x_j,t_n)=[Lu]_j^n+\frac{\tau}{2}(\frac{\partial^2u}{\partial t^2})_j^n-\frac{ah^2}{12}(\frac{\partial^4u}{\partial x^4})_j^n+O(\tau^2)+O(h^3)$

利用方程关系化简

已知抛物型方程 $u_t = au_{xx}+f(x)$ ，对其求导可得：

$\begin{align} u_{tx}&=au_x^{(4)}+f''\\ u_{tt}&=au_{xxt} \end{align}$

进而推出

$\begin{align} au_x^{(4)} &= u_{tx}-f''\\ &=\frac{1}{a}u_{tt}-f'' \end{align}$

将 $\displaystyle au_x^{(4)} = \frac{1}{a}u_{tt}-f’’$ 代入截断误差表达式 $R_j^n(u)$ ：

$R_j^n(u)=\frac{\tau}{2}(\frac{\partial^2u}{\partial t^2})_j^n-\frac{h^2}{12}[\frac{1}{a}(\frac{\partial^2u}{\partial t^2})_j^n - f_j'']+O(\tau^2)+O(h^3)$

整理得

$R_j^n(u)=\tau(\frac{\partial^2u}{\partial t^2})_j^n(\frac{1}{2}-\frac{1}{12r})+\frac{h^2}{12}f_j''+O(\tau)+O(h^3)$

其中 $\displaystyle r = \frac{a\tau}{h^2}$

当 $\tau$ 和 $h$ 趋于 0 时

$R_j^n(u)=O(\tau + h^2)$

表明向前差分格式对时间是一阶精度，对空间是二阶精度

向后差分格式

对于抛物型方程 $\displaystyle \frac{\partial u}{\partial t}=a\frac{\partial^2u}{\partial x^2}+f(x)$ ，时间方向用向后差商 $\displaystyle \frac{u_j^{n + 1}-u_j^n}{\tau}$ 近似 $\displaystyle \frac{\partial u}{\partial t}$ ，空间方向用中心差商 $\displaystyle \frac{u_{j - 1}^{n + 1}-2u_j^{n + 1}+u_{j + 1}^{n + 1}}{h^2}$ 近似 $\displaystyle \frac{\partial^2u}{\partial x^2}$ ，得到向后差分格式：

$\frac{u_j^{n + 1}-u_j^n}{\tau}=a\frac{u_{j - 1}^{n + 1}-2u_j^{n + 1}+u_{j + 1}^{n + 1}}{h^2}+f_j$

整理可得

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

令 $\displaystyle r = \frac{a\tau}{h^2}$ ，进一步化为

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

该格式是隐格式，因为在求解 $n + 1$ 时刻的 $u_j^{n + 1}$ 时，需要同时考虑 $j - 1$、$j$、$j + 1$ 节点在 $n + 1$ 时刻的值，即需求解联立线性方程组。已知 $u_j^n$ ，结合边界条件（如 $u_0^{n + 1}=0$ ，$u_J^{n + 1}=0$ ）来求解 $u_j^{n + 1}$

截断误差

在点 $(x_j,t_{n + 1})$ 处对相关函数进行泰勒展开。代入 $L_h^{(2)}u(x_j,t_{n + 1})$ 并与 $[Lu]_j^n$ 作差：

$L_h^{(2)}u(x_j,t_{n + 1})=\frac{u_j^{n + 1}-u_j^n}{\tau}-a\frac{u_{j - 1}^{n + 1}-2u_j^{n + 1}+u_{j + 1}^{n + 1}}{h^2}$ $\begin{align} R_j^n(u) &= L_h^{(2)}u(x_j,t_{n + 1})-[Lu]_j^n \\ &=-\tau(\frac{\partial^2u}{\partial t^2})_j^n(\frac{1}{2}+\frac{1}{12r})+O(\tau^2 + h^2)\\ &=O(\tau + h^2) \end{align}$

Crank - Nicolson 格式

Crank - Nicolson 格式是通过将向前差分格式和向后差分格式进行算术平均得到的。对于抛物型方程 $\displaystyle \frac{\partial u}{\partial t}=a\frac{\partial^2u}{\partial x^2}+f(x)$ ，其差分格式表达式为

$\frac{u_j^{n + 1}-u_j^n}{\tau}=a[\frac{1}{2}\frac{u_{j - 1}^{n + 1}-2u_j^{n + 1}+u_{j + 1}^{n + 1}}{h^2}+\frac{1}{2}\frac{u_{j - 1}^n - 2u_j^n + u_{j + 1}^n}{h^2}]+f_j$

进一步整理为

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

其中 $\displaystyle r = \frac{a\tau}{h^2}$ ，再化为

$-\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$

这是一种隐格式，也叫对称六点差分格式

边界条件：已知边界条件如 $u_0^{n + 1}=0$ ，$u_J^{n + 1}=0$ ，求解时需解线性方程组

截断误差

在点 $\displaystyle (x_j,t_{n+\frac{1}{2}})$ 处展开分析截断误差 $R_j^n(u)$ 。根据泰勒展开，将差分算子 $L_h^{(3)}$ 作用于 $u(x_j,t_{n+\frac{1}{2}})$ 与微分算子 $Lu$ 在该点的值作差

$\begin{align} R_j^n(u) &= L_h^{(3)}u(x_j,t_{n + 1})-[Lu]_j^{n+\frac{1}{2}} \\ &=(\frac{\partial u}{\partial t})_j^{n+\frac{1}{2}}-a(\frac{\partial^2u}{\partial t^2})_j^{n+\frac{1}{2}} +O(\tau^2)+O(h^2)\\ &=O(\tau^2 + h^2) \end{align}$

这表明该格式在时间和空间上均具有二阶精度，相比向前差分格式（$O(\tau + h^2)$ ）和向后差分格式（$O(\tau + h^2)$ ），精度更高

加权平均形式

该格式还可表示为加权平均形式

$\frac{u_j^{n + 1}-u_j^n}{\tau}=\frac{a}{h^2}[\theta\delta_x^2u_j^{n + 1}+(1 - \theta)\delta_x^2u_j^n]$

其中

$\delta_x^2u_j^n = u_{j - 1}^n - 2u_j^n + u_{j + 1}^n$

当 $\displaystyle \theta = \frac{1}{2}$ 时，就是 Crank - Nicolson 格式。通过调整 $\theta$ 值，可以研究不同差分格式的稳定性、精度等性质

Richardson 格式

Richardson 格式表达式为

$\frac{u_j^{n + 1}-u_j^{n - 1}}{2\tau}=a\frac{u_{j - 1}^n - 2u_j^n + u_{j + 1}^n}{h^2}+f_j$

整理可得

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

其中 $\displaystyle r = \frac{a\tau}{h^2}$

截断误差：截断误差为 $O(\tau^2 + h^2)$ ，在时间和空间上均具有二阶精度

差分格式考量因素

计算简单：格式的计算复杂度低，易于编程实现和计算

收敛性与敛速：格式的解是否收敛到原问题的解，以及收敛速度的快慢

稳定性：计算过程中误差是否会随时间或迭代次数无限增长，不稳定的格式即使有高精度也难以应用

稳定性分析

以Richardson 格式的稳定性为例

设 $e_j^n$ 表示 $u_j^n$ 的误差，假设 $f_j^n$ 精确，误差源于初始值，即 $e_j^{-1}=0$ ，$e_0^0 = \varepsilon$ ，$e_j^0 = 0$（$j\neq0$ ）。根据 Richardson 格式，误差满足

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

从相关误差传播分析可知，当 $n\rightarrow\infty$ 时，误差会无限增长，所以 Richardson 格式是不稳定的。这也体现了在构造差分格式时，稳定性是一个关键考量因素，不稳定的格式在实际计算中会产生无意义的结果