四种有限差分方法的数学推导

简介

本文默认读者对有限差分方法有一定了解，因此不做过多的前置知识介绍。考虑如下初值问题：

{\begin{cases} 𝑢^{'} = 𝑓 (𝑡, 𝑢), 𝑡 \in [0, 𝑇] \\ 𝑢 (0) = 𝑢^{(0)} \end{cases}

将 $[0, 𝑇]$ 划分为 $𝑀$ 个等距区间， $ℎ = \frac{𝑇}{𝑀}$ ，我们有以下几种有限差分格式用于求其数值解：

向前欧拉格式：
$𝑢^{(𝑛 + 1)} = 𝑢^{(𝑛)} + ℎ 𝑓 (𝑡^{(𝑛)}, 𝑢^{(𝑛)}) .$
向后欧拉格式：
$𝑢^{(𝑛 + 1)} = 𝑢^{(𝑛)} + ℎ 𝑓 (𝑡^{(𝑛 + 1)}, 𝑢^{(𝑛 + 1)}) .$
跃点（Leap frog）格式：
$𝑢^{(𝑛 + 1)} = 𝑢^{(𝑛 - 1)} + 2 ℎ 𝑓 (𝑡^{(𝑛)}, 𝑢^{(𝑛)}) .$
梯形格式：
$𝑢^{(𝑛 + 1)} = 𝑢^{(𝑛)} + ℎ \frac{𝑓 (𝑡^{(𝑛)}, 𝑢^{(𝑛)}) + 𝑓 (𝑡^{(𝑛 + 1)}, 𝑢^{(𝑛 + 1)})}{2} .$

接下来，对以上四种有限差分格式，我们将进行其准确性（Accuracy）与稳定性（Stability）的数学推导。

准确性推导

首先定义局部截断误差（Local Truncation Error） $𝑅^{(𝑛)}$ ，其定义为：其他点处的精确解在 $𝑡^{(𝑛)}$ 处按照差分格式计算得到的结果与 $𝑡^{(𝑛)}$ 处的精确解的差值。

警告：上述局部截断误差与通常定义的局部截断误差相差一个 $ℎ$ ，参考时请注意。

向前欧拉格式

向前欧拉格式的局部截断误差 $𝑅^{(𝑛)}$ 可表示为如下等式

𝑅^{(𝑛)} = \frac{𝑢 (𝑡^{(𝑛 + 1)}) - 𝑢 (𝑡^{(𝑛)})}{ℎ} - 𝑢^{'} (𝑡^{(𝑛)}) .

对 $𝑢 (𝑡^{(𝑛 + 1)})$ 进行 Taylor 展开有

𝑢 (𝑡^{(𝑛 + 1)}) = 𝑢 (𝑡^{(𝑛)}) + ℎ 𝑢^{'} (𝑡^{(𝑛)}) + \frac{ℎ^{2}}{2} 𝑢^{″} (𝜉^{(𝑛)}), 𝜉^{(𝑛)} \in [𝑡^{(𝑛)}, 𝑡^{(𝑛 + 1)}] .

于是有

𝑅^{(𝑛)} = \frac{ℎ}{2} 𝑢^{″} (𝜉^{(𝑛)}) = 𝑂 (ℎ), 𝜉^{(𝑛)} \in [𝑡^{(𝑛)}, 𝑡^{(𝑛 + 1)}] .

误差 $𝐸^{(𝑛)}$ 可表示为如下形式：

\frac{𝐸^{(𝑛 + 1)} - 𝐸^{(𝑛)}}{ℎ} + 𝑅^{(𝑛)} = 𝑓 (𝑡^{(𝑛)}, 𝑢 (𝑡^{(𝑛)})) - 𝑓 (𝑡^{(𝑛)}, 𝑢^{(𝑛)}) .

令 $𝑅 = {max}_{𝑛} | 𝑅^{(𝑛)} |$ ，对误差 $𝐸^{(𝑛)}$ 有如下估计：

\begin{aligned} | 𝐸^{(𝑛 + 1)} | & \leq | 𝐸^{(𝑛)} | + ℎ 𝑅 + ℎ 𝐿 | 𝑢 (𝑡^{(𝑛)}) - 𝑢^{(𝑛)} | = (1 + ℎ 𝐿) | 𝐸^{(𝑛)} | + ℎ 𝑅 \\ \leq {(1 + ℎ 𝐿)}^{𝑛 + 1} | 𝐸^{(0)} | + \frac{𝑅}{𝐿} ({(1 + ℎ 𝐿)}^{𝑛 + 1} - 1) \\ \leq 𝑒^{𝐿 𝑇} | 𝐸^{(0)} | + \frac{𝑅}{𝐿} (𝑒^{𝐿 𝑇} - 1) \leq 𝑂 (ℎ), 𝑛 = 0, 1, \dots, 𝑀 - 1 . \end{aligned}

向后欧拉格式

向后欧拉格式的局部截断误差 $𝑅^{(𝑛)}$ 可表示为如下等式

𝑅^{(𝑛)} = \frac{𝑢 (𝑡^{(𝑛)}) - 𝑢 (𝑡^{(𝑛 - 1)})}{ℎ} - 𝑢^{'} (𝑡^{(𝑛)}) .

对 $𝑢 (𝑡^{(𝑛 - 1)})$ 进行 Taylor 展开有

𝑢 (𝑡^{(𝑛 - 1)}) = 𝑢 (𝑡^{(𝑛)}) - ℎ 𝑢^{'} (𝑡^{(𝑛)}) + \frac{ℎ^{2}}{2} 𝑢^{″} (𝜉^{(𝑛)}), 𝜉^{(𝑛)} \in [𝑡^{(𝑛 - 1)}, 𝑡^{(𝑛)}] .

于是有

𝑅^{(𝑛)} = - \frac{ℎ}{2} 𝑢^{″} (𝜉^{(𝑛)}) = 𝑂 (ℎ), 𝜉^{(𝑛)} \in [𝑡^{(𝑛 - 1)}, 𝑡^{(𝑛)}] .

误差 $𝐸^{(𝑛)}$ 可表示为如下形式：

\frac{𝐸^{(𝑛)} - 𝐸^{(𝑛 - 1)}}{ℎ} - 𝑅^{(𝑛)} = 𝑓 (𝑡^{(𝑛)}, 𝑢 (𝑡^{(𝑛)})) - 𝑓 (𝑡^{(𝑛)}, 𝑢^{(𝑛)}) .

令 $𝑅 = {max}_{𝑛} | 𝑅^{(𝑛)} |$ ，假定 $ℎ 𝐿 < 1$ ，对误差 $𝐸^{(𝑛)}$ 有如下估计：

\begin{aligned} | 𝐸^{(𝑛)} | & \leq {(1 - ℎ 𝐿)}^{- 1} | 𝐸^{(𝑛 - 1)} | + ℎ 𝑅 \leq {(1 - ℎ 𝐿)}^{- 𝑛} | 𝐸^{(0)} | + \frac{𝑅}{𝐿} (1 - ℎ 𝐿) ({(1 - ℎ 𝐿)}^{- 𝑛} - 1) \\ < {(1 + \frac{ℎ 𝐿}{1 - ℎ 𝐿})}^{𝑛} | 𝐸^{(0)} | + \frac{𝑅}{𝐿} ({(1 + \frac{ℎ 𝐿}{1 - ℎ 𝐿})}^{𝑛} - 1) \\ \leq exp (\frac{𝐿 𝑇}{1 - ℎ 𝐿}) | 𝐸^{(0)} | + \frac{𝑅}{𝐿} (exp (\frac{𝐿 𝑇}{1 - ℎ 𝐿}) - 1) \leq 𝑂 (ℎ), 𝑛 = 1, 2, 3, \dots, 𝑀 . \end{aligned}

跃点格式

跃点格式的局部截断误差 $𝑅^{(𝑛)}$ 可表示为如下等式

𝑅^{(𝑛)} = \frac{𝑢 (𝑡^{(𝑛 + 1)}) - 𝑢 (𝑡^{(𝑛 - 1)})}{2 ℎ} - 𝑢^{'} (𝑡^{(𝑛)}) .

对 $𝑢 (𝑡^{(𝑛 - 1)})$ 和 $𝑢 (𝑡^{(𝑛 + 1)})$ 进行 Taylor 展开有

𝑢 (𝑡^{(𝑛 - 1)}) = 𝑢 (𝑡^{(𝑛)}) - ℎ 𝑢^{'} (𝑡^{(𝑛)}) + \frac{ℎ^{2}}{2} 𝑢^{″} (𝜉_{1}^{(𝑛)}), 𝜉_{1}^{(𝑛)} \in [𝑡^{(𝑛 - 1)}, 𝑡^{(𝑛)}] .

𝑢 (𝑡^{(𝑛 + 1)}) = 𝑢 (𝑡^{(𝑛)}) + ℎ 𝑢^{'} (𝑡^{(𝑛)}) + \frac{ℎ^{2}}{2} 𝑢^{″} (𝜉_{2}^{(𝑛)}), 𝜉_{2}^{(𝑛)} \in [𝑡^{(𝑛)}, 𝑡^{(𝑛 + 1)}] .

于是有

𝑅^{(𝑛)} = \frac{ℎ}{4} (𝑢^{″} (𝜉_{1}^{(𝑛)}) + 𝑢^{″} (𝜉_{2}^{(𝑛)})) = 𝑂 (ℎ), 𝜉_{1}^{(𝑛)} \in [𝑡^{(𝑛 - 1)}, 𝑡^{(𝑛)}], 𝜉_{2}^{(𝑛)} \in [𝑡^{(𝑛)}, 𝑡^{(𝑛 + 1)}] .

误差 $𝐸^{(𝑛)}$ 可表示为如下形式：

\frac{𝐸^{(𝑛 + 1)} - 𝐸^{(𝑛 - 1)}}{2 ℎ} - 𝑅^{(𝑛)} = 𝑓 (𝑡^{(𝑛)}, 𝑢 (𝑡^{(𝑛)})) - 𝑓 (𝑡^{(𝑛)}, 𝑢^{(𝑛)}) .

令 $𝑅 = {max}_{𝑛} | 𝑅^{(𝑛)} |$ ，对误差 $𝐸^{(𝑛)}$ 有如下估计：

| 𝐸^{(𝑛 + 1)} | \leq | 𝐸^{(𝑛 - 1)} | + 2 ℎ 𝐿 | 𝐸^{(𝑛)} | + 2 ℎ 𝑅 .

令 $𝜇 = \sqrt{1 + ℎ^{2} 𝐿^{2}} - ℎ 𝐿$ ， $𝜆 = \sqrt{1 + ℎ^{2} 𝐿^{2}} + ℎ 𝐿$ ，则 $𝜇 𝜆 = 1$ ，且 $𝜆 - 𝜇 = 2 ℎ 𝐿$ 。于是上式可改写为如下形式：

\begin{aligned} | 𝐸^{(𝑛 + 1)} | + 𝜇 | 𝐸^{(𝑛)} | & \leq 𝜆 (| 𝐸^{(𝑛)} | + 𝜇 | 𝐸^{(𝑛 - 1)} |) + 2 ℎ 𝑅 \leq \dots \leq 𝜆^{𝑛} (| 𝐸^{(1)} | + 𝜇 | 𝐸^{(0)} |) + \frac{𝜆^{𝑛} - 1}{𝜆 - 1} 2 ℎ 𝑅 \\ \leq {(1 + 2 ℎ 𝐿)}^{𝑛} (| 𝐸^{(1)} | + 𝜇 | 𝐸^{(0)} |) + \frac{{(1 + 2 ℎ 𝐿)}^{𝑛} - 1}{𝜆 - 1} 2 ℎ 𝑅 \\ \leq 𝑒^{2 𝐿 𝑇} (| 𝐸^{(1)} | + 𝜇 | 𝐸^{(0)} |) + \frac{𝑒^{2 𝐿 𝑇} - 1}{𝜆 - 1} 2 ℎ 𝑅 . \end{aligned}

不妨假定上式中的 $𝐸^{(1)} = 𝑢 (𝑡^{(1)}) - 𝑢^{(1)}$ 是由向前欧拉格式得到的，注意到 $ℎ \to 0$ 时 $\frac{ℎ}{𝜆 - 1} \to \frac{1}{𝐿}$ ，于是有

| 𝐸^{(𝑛 + 1)} | \leq | 𝐸^{(𝑛 + 1)} | + 𝜇 | 𝐸^{(𝑛)} | \leq 𝑂 (ℎ), 𝑛 = 1, 2, 3, \dots, 𝑀 - 1 .

梯形格式

梯形格式的局部截断误差 $𝑅^{(𝑛)}$ 可表示为如下等式

𝑅^{(𝑛)} = \frac{𝑢 (𝑡^{(𝑛 + 1)}) - 𝑢 (𝑡^{(𝑛)})}{ℎ} - \frac{𝑢^{'} (𝑡^{(𝑛 + 1)}) + 𝑢^{'} (𝑡^{(𝑛)})}{2} .

对 $𝑢 (𝑡^{(𝑛 + 1)})$ 和 $𝑢^{'} (𝑡^{(𝑛 + 1)})$ 进行 Taylor 展开有

𝑢 (𝑡^{(𝑛 + 1)}) = 𝑢 (𝑡^{(𝑛)}) + ℎ 𝑢^{'} (𝑡^{(𝑛)}) + \frac{ℎ^{2}}{2} 𝑢^{″} (𝑡^{(𝑛)}) + \frac{ℎ^{3}}{6} 𝑢^{‴} (𝜉_{1}^{(𝑛)}), 𝜉_{1}^{(𝑛)} \in [𝑡^{(𝑛)}, 𝑡^{(𝑛 + 1)}] .

𝑢^{'} (𝑡^{(𝑛 + 1)}) = 𝑢^{'} (𝑡^{(𝑛)}) + ℎ 𝑢^{″} (𝑡^{(𝑛)}) + \frac{ℎ^{2}}{2} 𝑢^{‴} (𝜉_{2}^{(𝑛)}), 𝜉_{2}^{(𝑛)} \in [𝑡^{(𝑛)}, 𝑡^{(𝑛 + 1)}] .

于是有

𝑅^{(𝑛)} = \frac{ℎ^{2}}{12} (2 𝑢^{‴} (𝜉_{1}^{(𝑛)}) - 3 𝑢^{‴} (𝜉_{2}^{(𝑛)})) = 𝑂 (ℎ^{2}), 𝜉_{1}^{(𝑛)} \in [𝑡^{(𝑛)}, 𝑡^{(𝑛 + 1)}], 𝜉_{2}^{(𝑛)} \in [𝑡^{(𝑛)}, 𝑡^{(𝑛 + 1)}] .

误差 $𝐸^{(𝑛)}$ 可表示为如下形式：

\frac{𝐸^{(𝑛 + 1)} - 𝐸^{(𝑛)}}{ℎ} - 𝑅^{(𝑛)} = \frac{𝑓 (𝑡^{(𝑛 + 1)}, 𝑢 (𝑡^{(𝑛 + 1)})) + 𝑓 (𝑡^{(𝑛)}, 𝑢 (𝑡^{(𝑛)}))}{2} - \frac{𝑓 (𝑡^{(𝑛 + 1)}, 𝑢^{(𝑛 + 1)}) + 𝑓 (𝑡^{(𝑛)}, 𝑢^{(𝑛)})}{2} .

令 $𝑅 = {max}_{𝑛} | 𝑅^{(𝑛)} |$ ，假定 $ℎ 𝐿 < 2$ ，对误差 $𝐸^{(𝑛)}$ 有如下估计：

| 𝐸^{(𝑛 + 1)} | \leq | 𝐸^{(𝑛)} | + \frac{ℎ 𝐿}{2} (| 𝐸^{(𝑛 + 1)} | + | 𝐸^{(𝑛)} |) + ℎ 𝑅 .

\begin{aligned} | 𝐸^{(𝑛)} | & \leq \frac{2 + ℎ 𝐿}{2 - ℎ 𝐿} | 𝐸^{(𝑛 - 1)} | + \frac{2 ℎ 𝑅}{2 - ℎ 𝐿} \leq \frac{{(2 + ℎ 𝐿)}^{𝑛}}{{(2 - ℎ 𝐿)}^{𝑛}} | 𝐸^{(0)} | + \frac{𝑅}{𝐿} (\frac{{(2 + ℎ 𝐿)}^{𝑛}}{{(2 - ℎ 𝐿)}^{𝑛}} - 1) \\ \leq exp (\frac{2 𝐿 𝑇}{2 - ℎ 𝐿}) | 𝐸^{(0)} | + \frac{𝑅}{𝐿} (exp (\frac{2 𝐿 𝑇}{2 - ℎ 𝐿}) - 1) \leq 𝑂 (ℎ^{2}), 𝑛 = 1, 2, 3, \dots, 𝑀 . \end{aligned}

稳定性推导

考虑模型问题

\frac{d 𝑢}{d 𝑡} = 𝜆 𝑢,

其中 $𝜆$ 为常数。定义 $𝜀^{(𝑛)}$ 为舍入误差， $𝑧 = ℎ 𝜆$ 。我们称差分格式的绝对稳定区域（Absolute Stability Region）为复平面上所有使差分格式的解稳定的 $𝑧$ 的集合。

向前欧拉格式

将格式代入模型问题可得

𝑢^{(𝑛 + 1)} = (1 + ℎ 𝜆) 𝑢^{(𝑛)} .

考虑到舍入误差，我们有 ${\bar{𝑢}}^{(𝑛)}$ 而非 $𝑢^{(𝑛)}$ ：

{\bar{𝑢}}^{(𝑛 + 1)} = {\bar{𝑢}}^{(𝑛)} + ℎ 𝜆 {\bar{𝑢}}^{(𝑛)} + 𝜀^{(𝑛)} .

误差 $𝐸^{(𝑛)} ≔ {\bar{𝑢}}^{(𝑛)} - 𝑢^{(𝑛)}$ 满足

𝐸^{(𝑛 + 1)} = 𝐸^{(𝑛)} + ℎ 𝜆 𝐸^{(𝑛)} + 𝜀^{(𝑛)} = \dots = {(1 + ℎ 𝜆)}^{𝑛 + 1} 𝐸^{(0)} + {(1 + ℎ 𝜆)}^{𝑛} 𝜀^{(0)} + \dots + 𝜀^{(𝑛)} .

假设 $𝜀^{(𝑛)} < 𝜀$ ，若 $| 1 + ℎ 𝜆 | < 1$ ，则有

| 𝐸^{(𝑛 + 1)} | \leq {| 1 + ℎ 𝜆 |}^{𝑛 + 1} | 𝐸^{(0)} | + \frac{{| 1 + ℎ 𝜆 |}^{𝑛 + 1} + 1}{| ℎ 𝜆 |} 𝜀 \to 𝑐 𝜀,

其中 $𝑐$ 为常数，所以本格式的绝对稳定区域为 $| 1 + 𝑧 | < 1$ 。

向后欧拉格式

将格式代入模型问题可得

𝑢^{(𝑛)} = (1 - ℎ 𝜆) 𝑢^{(𝑛 + 1)} .

考虑到舍入误差，我们有 ${\bar{𝑢}}^{(𝑛)}$ 而非 $𝑢^{(𝑛)}$ ：

{\bar{𝑢}}^{(𝑛 + 1)} = {\bar{𝑢}}^{(𝑛)} + ℎ 𝜆 {\bar{𝑢}}^{(𝑛 + 1)} + 𝜀^{(𝑛)} .

误差 $𝐸^{(𝑛)} ≔ {\bar{𝑢}}^{(𝑛)} - 𝑢^{(𝑛)}$ 满足

\begin{aligned} 𝐸^{(𝑛 + 1)} & = 𝐸^{(𝑛)} + ℎ 𝜆 𝐸^{(𝑛 + 1)} + 𝜀^{(𝑛)} = {(1 - ℎ 𝜆)}^{- 1} 𝐸^{(𝑛)} + {(1 - ℎ 𝜆)}^{- 1} 𝜀^{(𝑛)} = \dots \\ = {(1 - ℎ 𝜆)}^{- (𝑛 + 1)} 𝐸^{(0)} + {(1 - ℎ 𝜆)}^{- (𝑛 + 1)} 𝜀^{(0)} + \dots + {(1 - ℎ 𝜆)}^{- 1} 𝜀^{(𝑛)} . \end{aligned}

假设 $𝜀^{(𝑛)} < 𝜀$ ，若 $| 1 - ℎ 𝜆 | > 1$ ，则有

| 𝐸^{(𝑛 + 1)} | \leq {| 1 - ℎ 𝜆 |}^{- (𝑛 + 1)} | 𝐸^{(0)} | + \frac{{| 1 - ℎ 𝜆 |}^{- 𝑛} + 1}{| ℎ 𝜆 |} 𝜀 \to 𝑐 𝜀,

其中 $𝑐$ 为常数，所以本格式的绝对稳定区域为 $| 1 - 𝑧 | > 1$ 。

跃点格式

将格式代入模型问题可得

𝑢^{(𝑛 + 1)} = 𝑢^{(𝑛 - 1)} + 2 ℎ 𝜆 𝑢^{(𝑛)} .

考虑到舍入误差，我们有 ${\bar{𝑢}}^{(𝑛)}$ 而非 $𝑢^{(𝑛)}$ ：

{\bar{𝑢}}^{(𝑛 + 1)} = {\bar{𝑢}}^{(𝑛 - 1)} + 2 ℎ 𝜆 {\bar{𝑢}}^{(𝑛)} + 𝜀^{(𝑛 - 1)} .

误差 $𝐸^{(𝑛)} ≔ {\bar{𝑢}}^{(𝑛)} - 𝑢^{(𝑛)}$ 满足

\begin{aligned} 𝐸^{(𝑛 + 1)} + 𝛼 𝐸^{(𝑛)} & = \frac{1}{𝛼} (𝐸^{(𝑛)} + 𝛼 𝐸^{(𝑛 - 1)}) + 𝜀^{(𝑛 - 1)} = \dots \\ = \frac{1}{𝛼^{𝑛}} (𝐸^{(1)} + 𝛼 𝐸^{(0)}) + \frac{1}{𝛼^{𝑛 - 1}} 𝜀^{(0)} + \dots + 𝜀^{(𝑛 - 1)}, \end{aligned}

其中 $𝛼 = {(ℎ^{2} 𝜆^{2} + 1)}^{\frac{1}{2}} - ℎ 𝜆$ ，多值函数 $𝑧^{\frac{1}{2}}$ 取将正数映射到正数的解析分支。随后我们有

\begin{aligned} 𝐸^{(𝑛 + 1)} & = 𝛼^{- 𝑛} (𝐸^{(1)} + 𝛼 𝐸^{(0)}) + \sum_{𝑘 = 0}^{𝑛 - 1} 𝛼^{𝑘 - (𝑛 - 1)} 𝜀^{(𝑘)} - 𝛼 𝐸^{(𝑛)} = \dots \\ = \frac{{(- 1)}^{𝑛 - 1} 𝑎^{𝑛} + 𝑎^{- 𝑛}}{𝑎^{2} + 1} (𝐸^{(1)} + 𝛼 𝐸^{(0)}) + \sum_{𝑘 = 0}^{𝑛 - 1} \frac{{(- 1)}^{𝑛 - 1 - 𝑘} 𝑎^{𝑛 + 1 - 𝑘} + 𝛼^{𝑘 - (𝑛 - 1)}}{𝑎^{2} + 1} 𝜀^{(𝑘)} + {(- 1)}^{𝑛} 𝛼 𝐸^{(1)} . \end{aligned}

假设 $𝜀^{(𝑛)} < 𝜀$ ，若 $| 1 + ℎ 𝜆 | > 1$ ，则有

| 𝐸^{(𝑛 + 1)} | \leq \frac{{| 𝛼 |}^{𝑛} + {| 𝛼 |}^{- 𝑛}}{| 𝛼^{2} + 1 |} | 𝐸^{(1)} + 𝛼 𝐸^{(0)} | + | 𝛼 | | 𝐸^{(1)} | + \frac{| 𝑎 |}{| 𝛼^{2} + 1 |} (\frac{| 𝛼^{𝑛} + 1 |}{| 𝛼^{- 1} + 1 |} + \frac{| 𝛼^{- 𝑛} - 1 |}{| 𝛼 - 1 |}) 𝜀 .

上式同时出现 ${| 𝛼 |}^{𝑛}$ 与 ${| 𝛼 |}^{- 𝑛}$ ，所以仅 $| 𝛼 | = 1$ 时跃点格式稳定。令 $𝑧 = 𝑎 + 𝑏 𝑖$ 有

{| 𝛼 |}^{2} = {(\sqrt{𝑟} cos 𝜃 - 𝑎)}^{2} + {(\sqrt{𝑟} sin 𝜃 - 𝑏)}^{2} = 𝑟 + 𝑎^{2} + 𝑏^{2} - 2 \sqrt{𝑟} (𝑎 cos 𝜃 + 𝑏 sin 𝜃) = 1,

其中

𝑟 = \sqrt{{(𝑎^{2} - 𝑏^{2} + 1)}^{2} + 4 𝑎^{2} 𝑏^{2}}, 𝜃 = \frac{arctan (\frac{2 𝑎 𝑏}{𝑎^{2} - 𝑏^{2} + 1})}{2} .

注意到仅 $𝑎 = 0$ 且 $| 𝑏 | \leq 1$ 时上式成立，所以跃点格式的绝对稳定区域为 $𝑧 = 𝑏 𝑖, 𝑏 \in [- 1, 1]$ 。

梯形格式

将格式代入模型问题可得

𝑢^{(𝑛 + 1)} = 𝑢^{(𝑛)} + \frac{ℎ 𝜆}{2} (𝑢^{(𝑛)} + 𝑢^{(𝑛 + 1)}) .

考虑到舍入误差，我们有 ${\bar{𝑢}}^{(𝑛)}$ 而非 $𝑢^{(𝑛)}$ ：

{\bar{𝑢}}^{(𝑛 + 1)} = {\bar{𝑢}}^{(𝑛)} + \frac{ℎ 𝜆}{2} {\bar{𝑢}}^{(𝑛)} + \frac{ℎ 𝜆}{2} {\bar{𝑢}}^{(𝑛 + 1)} + 𝜀^{(𝑛)} .

误差 $𝐸^{(𝑛)} ≔ {\bar{𝑢}}^{(𝑛)} - 𝑢^{(𝑛)}$ 满足

\begin{aligned} 𝐸^{(𝑛 + 1)} & = 𝐸^{(𝑛)} + \frac{ℎ 𝜆}{2} 𝐸^{(𝑛)} + \frac{ℎ 𝜆}{2} 𝐸^{(𝑛 + 1)} + 𝜀^{(𝑛)} = \frac{2 + ℎ 𝜆}{2 - ℎ 𝜆} 𝐸^{(𝑛)} + \frac{2}{2 - ℎ 𝜆} 𝜀^{(𝑛)} = \dots \\ = \frac{{(2 + ℎ 𝜆)}^{𝑛 + 1}}{{(2 - ℎ 𝜆)}^{𝑛 + 1}} 𝐸^{(0)} + \frac{{(2 + ℎ 𝜆)}^{𝑛}}{{(2 - ℎ 𝜆)}^{𝑛 + 1}} 2 𝜀^{(0)} + \dots + \frac{2}{2 - ℎ 𝜆} 𝜀^{(𝑛)} . \end{aligned}

假设 $𝜀^{(𝑛)} < 𝜀$ ，令 $𝑧 = 𝑎 + 𝑏 𝑖$ ，若 $𝑎 < 0$ ，则有

| 𝐸^{(𝑛 + 1)} | \leq \frac{| {(2 + ℎ 𝜆)}^{𝑛 + 1} |}{| {(2 - ℎ 𝜆)}^{𝑛 + 1} |} | 𝐸^{(0)} | + (\frac{| {(2 + ℎ 𝜆)}^{𝑛 + 1} |}{| {(2 - ℎ 𝜆)}^{𝑛 + 1} |} - 1) \to 𝑐 𝜀,

其中 $𝑐$ 为常数，所以本格式的绝对稳定区域为 $Re (𝑧) < 0$ 。

计算复杂性

记 $𝑛 = \frac{𝑇}{ℎ}$ 。

对于递推式中不含 $𝑓 (𝑡^{(𝑛 + 1)}, 𝑢^{(𝑛 + 1)})$ 的差分格式，其求解全部 $𝑢^{(𝑛)}$ 只需进行 $𝑛$ 次计算，即计算复杂度为 $𝑂 (𝑛)$ 。

反之对于递推式中含有 $𝑓 (𝑡^{(𝑛 + 1)}, 𝑢^{(𝑛 + 1)})$ 的差分格式，其计算每一个 $𝑢^{(𝑛)}$ 时都需要求解隐函数零点，这一过程是 $𝑂 (log 𝑛)$ 的，因此整个求解过程的计算复杂度为 $𝑂 (𝑛 log 𝑛)$ 。

数值实验

以下面的例子对四种差分格式进行数值实验。

{\begin{cases} 𝑢^{'} (𝑡) = 𝑓 (𝑡, 𝑢 (𝑡)) = - cos (5 𝜋 𝑡) tan (5 𝜋 𝑢) \\ 𝑢 (0) = 1 / 60 . \end{cases}

其中 $𝑡 \in [0, 1]$ 。

上述初值问题的解析解为

𝑢 (𝑡) = \frac{1}{5 𝜋} arcsin (\frac{\sqrt{6} - \sqrt{2}}{4 exp (sin (5 𝜋 𝑡))}) .

网格宽度 $ℎ = 0.01$ 时，使用差分格式求解上述初值问题，可得如下图表：

对比五种数值方法（u(t)、FE、BE、LF、Tz）在 t=0 到 1 间计算结果的图像。BE 曲线升至 0.1 后保持水平；u(t)、LF、Tz 曲线基本重合，有两个波峰，LF 曲线呈锯齿状；FE 曲线走势相似但峰值较低。 — 图 1 数值实验结果

结论

由上述数值实验的实验结果，可以看出虽然向前 Euler 格式、向后 Euler 格式和跃点格式的准确性都是 $𝑂 (ℎ)$ 的，但是跃点格式的实际准确性明显优于向前 Euler 格式和向后 Euler 格式。

而准确性为 $𝑂 (ℎ^{2})$ 的梯形格式则明显优于所有准确性为 $𝑂 (ℎ)$ 的差分格式。

附录

向前欧拉格式：

1
function [x, result] = forward_euler(T, step, u0, f)
2
    x = linspace(0, T, step + 1);
3
    h = T / step;
4
    result = zeros(1, step + 1);
5
    result(1) = u0;
6
    for k = 2:step + 1
7
        result(k) = result(k - 1) + h * f(x(k - 1), result(k - 1));
8
    end
9
end

向后欧拉格式：

1
function [x, result] = backward_euler(T, step, u0, f)
2
    x = linspace(0, T, step + 1);
3
    h = T / step;
4
    result = zeros(1, step + 1);
5
    result(1) = u0;
6
    for k = 2:step + 1
7
        temp_func = @(uk) result(k - 1) + h * f(x(k), uk) - uk;
8
        result(k) = fzero(temp_func, result(k - 1));
9
    end
10
end

跃点格式：

1
function [x, result] = leap_frog(T, step, u0, u1, f)
2
    x = linspace(0, T, step + 1);
3
    h = T / step;
4
    result = zeros(1, step + 1);
5
    result(1) = u0;
6
    result(2) = u1;
7
    for k = 3:step + 1
8
        result(k) = result(k - 2) + 2 * h * f(x(k - 1), result(k - 1));
9
    end
10
end

梯形格式：

1
function [x, result] = trapezoidal(T, step, u0, f)
2
    x = linspace(0, T, step + 1);
3
    h = T / step;
4
    result = zeros(1, step + 1);
5
    result(1) = u0;
6
    for k = 2:step + 1
7
        temp_func = @(uk) result(k - 1) - uk + ...
8
            h * (f(x(k), uk) + f(x(k - 1), result(k - 1))) / 2;
9
        result(k) = fzero(temp_func, result(k - 1));
10
    end
11
end