微分方程数值解 第一章（1）

Euler法

对于常微分方程（ODE）问题：

$$\begin{cases} u' = f(t, u), & 0 \leq t \leq T \\ u(0) = u_0 \end{cases}$$

其中 $f(t, u)$ 关于 $u$ 满足 Lipschitz 条件

欧拉法推导：

连续问题离散化

设真解为 $u(t)$。将区间 $[0, T]$ 离散化，令 $\displaystyle \frac{T}{N}=h$，其中 $h$ 为步长，且 $0 = t_0 < t_1 < \cdots < t_N = T$，$t_n=nh$

已知 $u(t_0)=u_0$，$u’(t_0)=f(t_0, u_0)$

利用 Taylor 展开推导近似值

对 $u(t_1)$ 进行 Taylor 展开：$\displaystyle u(t_1)=u(t_0 + h)=u(t_0)+u’(t_0)h+\frac{u’’(\xi)}{2!}h^2$ ，其中 $\displaystyle \xi\in(t_0,t_1)$ ，记为 $\displaystyle u(t_0)+u’(t_0)h + R_0$ ，$\displaystyle R_0 = O(h^2)$ 是二阶小量

略去二阶小量 $R_0$ ，得到 $u(t_1)$ 的近似值 $\displaystyle u_1 = u_0+hf(t_0, u_0)\approx u(t_1)$

类似地，可得到 $u(t_2)$ 的近似值 $u_2 = u_1+hf(t_1, u_1)$

以此类推，得到递推公式：$u_{n + 1}=u_n+hf(t_n, u_n)$ ，这就是 Euler 法

显格式与隐格式

显格式：

$u_{n+1} = u_n + h f(t_n, u_n)$，变形为 $\displaystyle \frac{u_{n + 1}-u_n}{h}=f(t_n, u_n)$ ，即用差商 $\displaystyle \frac{u_{n + 1}-u_n}{h}$（向前差分）代替导数 $u’(t_n)$

隐格式：

$\displaystyle \frac{u_n - u_{n - 1}}{h}=f(t_n, u_n)$ ，也可写成 $\displaystyle u_{n + 1}=u_n+hf(t_{n + 1}, u_{n + 1})$ ，这里用向后差分 $\displaystyle \frac{u_n - u_{n - 1}}{h}$ 近似 $u’(t_n)$ ，隐格式一般需要迭代求解

迭代公式为

$$\displaystyle u\_{n + 1}^{[k + 1]}=u\_n+hf(t\_{n + 1}, u\_{n + 1}^{[k]})$$

数值积分法：ODE 数值解的等价积分形式

常微分方程 $u’(t)=f(t, u(t))$ ，$u(t_0)$ 已知，其等价积分形式为：

$$u(t)=u(t_0)+\int_{t_0}^{t}f(\tau, u(\tau))d\tau$$

通过对积分项进行不同的数值近似，可得到不同的数值解法

矩形公式

左矩形公式（显式 Euler 法）

公式：$\displaystyle u_1 = u_0 + f(t_0, u_0)h$

原理：用左端点的函数值 $f(t_0, u_0)$ 近似积分区间 $[t_0, t_0 + h]$ 上 $\displaystyle f(\tau, u(\tau))$ 的平均值，来计算积分 $\displaystyle \int_{t_0}^{t_0 + h}f(\tau, u(\tau))d\tau\approx hf(t_0, u_0)$ ，进而得到 $u(t_1)$（$t_1=t_0 + h$ ）的近似值 $u_1$

右矩形公式（隐式 Euler 法）

公式：$\displaystyle u_1 = u_0 + f(t_1, u_1)h$

原理：用右端点的函数值 $f(t_1, u_1)$ 近似积分区间 $[t_0, t_0 + h]$ 上 $f(\tau, u(\tau))$ 的平均值 ，但由于 $u_1$ 未知，该方法是隐式的，一般需迭代求解

梯形公式（改进 Euler 法）

公式：$\displaystyle u_{n + 1}=u_n+\frac{h}{2}[f(t_n, u_n)+f(t_{n + 1}, u_{n + 1})]$

原理：用梯形面积公式近似积分，即把积分区间 $[t_n, t_{n + 1}]$ 上 $f(\tau, u(\tau))$ 的平均值近似为 $\displaystyle \frac{1}{2}[f(t_n, u_n)+f(t_{n + 1}, u_{n + 1})]$ ，从而计算积分

$$\displaystyle \int_{t\_n}^{t\_{n + 1}}f(\tau, u(\tau))d\tau\approx\frac{h}{2}[f(t\_n, u\_n)+f(t\_{n + 1}, u\_{n + 1})]$$

Euler 误差分析

数值解与精确解

数值解：显式 Euler 法的迭代公式为 $\displaystyle u_{n + 1}=u_n+hf(t_n, u_n)$

精确解：对 $u(t_{n + 1})$ 在 $t_n$ 处进行 Taylor 展开，$\displaystyle u(t_{n + 1})=u(t_n)+hf(t_n, u(t_n))+\frac{u’’(t_n)}{2!}h^2+\frac{u’’’(\xi)}{3!}h^3$ ，$\xi\in(t_n, t_{n + 1})$

局部截断误差

定义：把精确解 $u(t)$ 代入差分格式得到的余项。差分格式 $\displaystyle L[u_n, h]=u_{n + 1}-u_n - hf(t_n, u_n)$ ，将真解代入得 $\displaystyle L[u(t_n), h]=u(t_{n + 1})-u(t_n)-hf(t_n, u(t_n))\neq0$ ，余项 $\displaystyle R_n=\frac{u’’(t_n)}{2!}h^2+\frac{u’’’(\xi)}{3!}h^3$ 。当 $h\to0$ 时，Euler 法的局部截断误差为 $O(h^2)$ ，意味着随着步长 $h$ 趋于 0 ，误差主要由 $h^2$ 阶项决定

改进 Euler 法局部截断误差

利用积分中值定理相关结论

$$\displaystyle \int_{a}^{b}v(x)dx=\frac{1}{2}(b - a)(v(a)+v(b))-\frac{(b - a)^3}{12}v''(\xi)$$

将真解代入改进 Euler 法格式 $\displaystyle u_{n + 1}-u_n-\frac{h}{2}[f(t_n, u_n)+f(t_{n + 1}, u_{n + 1})]=0$ ，得到余项为 $\displaystyle -\frac{h^3}{12}u’’’(\xi)$ ，$\xi\in(t_n, t_{n + 1})$ ，所以改进 Euler 法的局部截断误差为 $O(h^3)$ ，比 Euler 法精度更高

整体误差

记 $\displaystyle e_n = u(t_n)-u_n$

Euler 法：

$$\begin{align*} u_{n + 1}&=u_n+hf(t_n,u_n)\\ u(t_{n+1})&=u(t_n)+hf(t_n,u(t_n)) + R_n \end{align*}$$

作差：

$$\displaystyle e_{n + 1}=e_n+h[f(t_n,u_n)-f(t_n,u(t_n))]+R_n$$

记 $R=\max|R_n|$，由于 $f$ 满足 Lipschitz 条件得：

$$\begin{align*} |e_{n+1}|&\leq|e_n|+L\cdot h\cdot|e_n| + R\\ |e_n|&\leq(1 + Lh)|e_{n - 1}|+R\\ &\leq(1 + Lh)^2|e_{n - 2}|+(1 + Lh)R+R\\ &\vdots\\ &\leq(1 + Lh)^n|e_0|+R\sum_{j = 0}^{n - 1}(1 + Lh)^j\\ &=(1 + Lh)^n|e_0|+R\cdot\frac{(1 + Lh)^n-1}{Lh} \end{align*}$$

由 $\ln(1 + Lh)\leq Lh$，$\Rightarrow n\ln(1 + Lh)\leq n\cdot Lh\leq L(T - t_0)$，$\Rightarrow(1 + Lh)^n\leq e^{L(T - t_0)}$

于是

$$|e_n|\leq e^{L(T - t_0)}|e_0|+\frac{R}{Lh}(e^{L(T - t_0)}-1)，n = 1,2,\cdots,N$$

对 Euler 法，可取 $e_0 = 0$，$R = ch^2$

则

$$|e_n|\leq\frac{ch}{L}(e^{L(T - t_0)})\leq ch$$

即 $e_n = o(h)$

注意比较局部截断误差与整体误差的关系

稳定性

定义：中途误差是否连续地依赖于初始值误差。初始误差包括测量误差、舍入误差

Euler 法稳定性分析：

从初值 $u_0$，$v_0$ 算出的节点值分别记为 ${u_n}$，${v_n}$。

$$\begin{align*} u_n&=u_{n - 1}+hf(t_{n - 1},u_{n - 1})\\ v_n&=v_{n - 1}+hf(t_{n - 1},v_{n - 1}) \end{align*}$$

两式相减，记 $e_n = u_n - v_n$，则

$$e_n=e_{n - 1}+h[f(t_{n - 1},u_{n - 1})-f(t_{n - 1},v_{n - 1})]$$

由 Lipschitz 条件，有

$$|e_n|\leq|e_{n - 1}|+hL|e_{n - 1}|=(1 + Lh)|e_{n - 1}|$$

递推可得：

$$\begin{align*} |e_n|&\leq(1 + Lh)^n|e_0|\\ &\leq e^{L(T - t_0)}|e_0| \end{align*}$$

这表明 $e_n$ 连续地依赖于初始误差 $e_0$，即 Euler 法稳定

线性多步法

数值积分法

数值积分公式：

$$u(t_{n + 1})=u(t_n)+\int_{t_n}^{t_{n + 1}}f(t,u(t))dt$$

思想：适当地选取节点，用 Lagrange 插值多项式代替 $f(t,u(t))$

Adams 外插法（显式）

节点：$t_{n - k},\cdots,t_n$（$k + 1$ 个点），对应 $k + 1$ 步法

使用 $k$ 次 Lagrange 插值多项式

Adams 内插法（隐式）

节点：$\displaystyle t_{n - k},\cdots,t_n,t_{n + 1}$（未知 $t_{n + 1}$ 故仍为 $k + 1$ 步法）

使用 $k + 1$ 次 Lagrange 插值多项式

Adams 外插法

插值分解

将 $f(t, u(t))$ 分解为 $\displaystyle f(t, u(t)) = L_{n,k}(t) + r_{n,k}(t)$ ，其中 $\displaystyle L_{n,k}(t)$ 是插值多项式，$\displaystyle r_{n,k}(t)$ 是插值余项

根据积分的线性性质，有

$$\displaystyle u(t_{n + 1}) = u(t_n)+\int_{t_n}^{t_{n + 1}}L_{n,k}(t)dt+\int_{t_n}^{t_{n + 1}}r_{n,k}(t)dt$$

记为式 ①

离散格式

$$\displaystyle u_{n + 1}=u_n+\int_{t_n}^{t_{n + 1}}L_{n,k}(t)dt$$

是通过用插值多项式 $L_{n,k}(t)$ 近似 $f(t, u(t))$ 并积分得到的数值求解公式

局部截断误差

局部截断误差 $\displaystyle R_{n,k}=\int_{t_n}^{t_{n + 1}}r_{n,k}(t)dt$ ，它衡量了用插值多项式近似 $f(t, u(t))$ 并积分时产生的误差

Lagrange 插值公式回顾

Lagrange 插值多项式

对于给定节点 $x_0, x_1,\cdots, x_k$ ，函数 $f(x)$ 的 $k$ 次 Lagrange 插值多项式为 $\displaystyle L_k(x)=\sum_{i = 0}^{k}f(x_i)l_i(x)$ ，其中 $l_i(x)$ 是 Lagrange 基函数

Newton 型插值多项式

向前形式：$\displaystyle f(x_0)+f[x_0,x_1](x - x_0)+\cdots$

向后形式：$\displaystyle f(x_k)+f[x_{k - 1},x_k](x - x_k)+\cdots + f[x_0,\cdots,x_k](x - x_k)\cdots(x - x_0)$ ，其中 $\displaystyle f[x_i,\cdots,x_j]$ 是 $i$ 阶差商

差分

向前差分：$\displaystyle \Delta_+ f_j = f_{j + 1}-f_j$

向后差分：$\displaystyle \Delta_- f_j = f_j - f_{j - 1}$ ，且 $\displaystyle \Delta_-^k f_j=\Delta_+^k f_{-k}$

采用 Newton 向后插值公式

Newton 向后插值公式 $L_{n,k}(t)$ 是基于节点处函数值 $g(t_i)$ 及其差商，通过差商的线性组合来近似 $g(t)$ ，进而用于数值积分求解 $u(t_{n + 1})$ ，展开形式是利用差商和节点信息逐步构建，以逼近 $g(t)$ 在区间 $[t_n, t_{n + 1}]$ 上的积分

设 $t_{n - k},\cdots,t_{n - 1},t_n$ 为节点，令 $f(t, u(t)) = g(t)$ ，$t = t_n+\tau h$，$\tau\in[0,1]$

$$\begin{align*} L_{n,k}(t)&=L_{n,k}(t)=g(t_n)+g(t_n,t_{n-1})(t-t_n)+\cdots+g(t_n,\cdots,t_{n-k})(t-t_n)\cdots(t-t_{n-k+1})\\ &=g(t_n)+\frac{g(t_n-g_{n-1})}{t_n-t_{n-1}}+\cdots\cdots \end{align*}$$

Newton 向后插值多项式展开

把 $L_{n,k}(t)$ 在 $t = t_n+\tau h$ 处展开，$\displaystyle g_{n - j}=f_{n - j}=f(t_{n - j},u_{n - j})$ 是节点处的函数值。利用向前差分 $\displaystyle \Delta_+^i g_{n - i}$ 来构建插值多项式，$\displaystyle \binom{s}{j}=\frac{s(s - 1)\cdots(s - j + 1)}{j!}$ ，$\displaystyle \binom{s}{0}=1$ ，得到 $L_{n,k}(t)$ 的具体形式

$$L_{n,k}(t)=g(t_n)+\frac{\Delta_+^i g_{n - 1}}{h}\tau h+\frac{\Delta_+^i g_{n - 2}}{2!}\tau(\tau+1)+ \cdots+\frac{\tau(\tau+1)\cdots(\tau+k-1)}{k!}\Delta_+^i g_{n - k}$$

则

$$\begin{align*} L_{n,k}(t)&=L_{n,k}(t_n+\tau h) \\ &=\sum_{j=0}^k (-1)^j\binom{-\tau}{j}\Delta_+^i f_{n - j}\\ \end{align*}$$

代入积分并化简

将 $L_{n,k}(t)$ 代入 $\displaystyle u_{n + 1}=u_n+\int_{t_n}^{t_{n + 1}}L_{n,k}(t)dt$ ，通过换元 $t = t_n+\tau h$ ，$dt = h d\tau$ ，把积分区间从 $[t_n,t_{n + 1}]$ 转换为 $[0,1]$ 。然后对积分进行计算，提出与 $j$ 无关的项 $h$ 和 $\Delta_+^j f_{n - j}$ ，得到 $\displaystyle u_{n + 1}=u_n + h\sum_{j = 0}^{k}a_j\Delta_+^j f_{n - j}$ ，其中记 $\displaystyle a_j = (-1)^j\int_{0}^{1}\binom{-\tau}{j}d\tau$

具体步骤如下：

$$\begin{align*} u_{n+1}&=u_n+\int_{t_n}^{t_{n+1}} L_{n,k}(t)dt\\ &=u_n+\int_{0}^{1} L_{n,k}(t_n+\tau h)hd\tau\\ &=\displaystyle u_n + h\sum_{j = 0}^{k}(-1)^j\int_{0}^{1}\binom{-\tau}{j}d\tau\Delta_+^j f_{n - j}\\ &=\displaystyle u_n + h\sum_{j = 0}^{k}a_j\Delta_+^j f_{n - j} \end{align*}$$

计算系数（生成函数法）

把 $(-1)^j(-\tau)(-\tau - 1)\cdots(-\tau - j + 1)$ 看成是 $(1 - t)^{-\tau}$ 关于 $t$ 在 $t = 0$ 处的导数

利用泰勒级数展开

$$(1 - t)^{-\tau}=\sum_{j = 0}^{\infty}\frac{(-1)^j(-\tau)(-\tau - 1)\cdots(-\tau - j + 1)}{j!}t^j$$

对上述展开式的左右两端在 $[0,1]$ 上积分并进行一系列变形

得

$$\displaystyle -\frac{t}{(1 - t)\ln(1 - t)}=\sum_{j = 0}^{\infty}a_j t^{j + 1}$$

变形得

$$-\frac{\ln(1 - t)}{t}\sum_{j = 0}^{\infty}a_j t^{j + 1}=\frac{1}{1-t}=1+t+t^2+\cdots+t^n$$

由

$$-\frac{\ln(1 - t)}{t}=1+\frac{t}{2}+\frac{t^2}{3}+\cdots+\frac{t^{n-1}}{n}$$

代入得

$$(1+\frac{t}{2}+\frac{t^2}{3}+\cdots+\frac{t^{n-1}}{n})(a_0+a_1+\cdots+a_nt^n)=(1+t+t^2+\cdots+t^n)$$

比较系数得到 $a_j$

$$\begin{equation} \begin{cases} a_0=1 \\ a_1+\frac{1}{2}a_0= 1 \\ \quad \quad \vdots \\ a_n + \frac{1}{2}a_{n-1}+\frac{1}{3}a_{n-2}+\cdots+\frac{1}{n+1}a_{0} = 1 \end{cases} \end{equation}$$

计算 $\Delta_+^j f_{n-j}$

原理：利用向前差分算子 $\Delta_+$ 的性质，$\Delta_+=(S - I)$（$S$ 为移位算子，$\displaystyle Sf_n = f_{n + 1}$ ，$I$ 为恒等算子 ），通过二项式定理 $\displaystyle (a + b)^j=\sum_{l = 0}^{j}C_{j}^{l}a^{j - l}b^{l}$ 展开 $(E - I)^j$ 来计算 $\Delta_+^j f_{n - j}$

具体实现步骤：

$$\Delta_+ f_n=f_{n+1}-f_{n}=(S-I)f_n \\$$ $$\begin{align*} \Delta_+ f_n &=(S-I)^jf_{n-j} \\ &=\sum_{l = 0}^{j}C_{j}^{l}S^{J - l}(-I)^{l}f_{n-j}\\ &=\sum_{l = 0}^{j}(-1)^l\binom{j}{l}f_{n-l} \end{align*}$$

因而

$$\begin{align*} u_{n + 1}&=u_n + h\sum_{j = 0}^{k}a_j\Delta_+^j f_{n - j} \\ &=u_n + h\sum_{j = 0}^{k}a_j\sum_{l=0}^{j}(-1)^l\binom{j}{l}f_{n - l}\\ &=u_n + h\sum_{l=0}^{j}\sum_{j = 0}^{k}(-1)^la_j\binom{j}{l}f_{n - l}\\ &=u_n+h\sum_{j = 0}^{k}b_{kl}f_{n - l} \end{align*}$$

其中记 $\displaystyle b_{kl}=\sum_{j = 0}^{k}(-1)^la_j\binom{j}{l}$

$k=0$：此时公式为 $\displaystyle u_{n + 1}=u_n+hf(t_n,u_n)$ ，这就是显式 Euler 法，它是 Adams 外插法的最简单形式，基于当前点的函数值来预估下一点的值

$k=1$：$\displaystyle u_{n + 1}=u_n+\frac{h}{2}(3f_n - f_{n - 1})$ ，利用前两个节点的函数值来计算下一个节点的数值解，相比 Euler 法考虑了更多历史信息

$k=2$：$\displaystyle u_{n + 1}=u_n+\frac{h}{12}(23f_n - 16f_{n - 1}+5f_{n - 2})$ ，使用了前三个节点的函数值，精度进一步提高

$k=3$：$\displaystyle u_{n + 1}=u_n+\frac{h}{24}(55f_n - 59f_{n - 1}+37f_{n - 2}-9f_{n - 3})$ ，借助更多历史节点信息来逼近真实解

Adams 外插法局部截断误差

插值余项公式：对于函数 $f(x)$ 的 $n$ 次插值多项式 $P_n(x)$

插值余项

$$\displaystyle E(f,x)=f(x)-P_n(x)=\frac{(x - x_0)(x - x_1)\cdots(x - x_n)}{(n + 1)!}f^{(n + 1)}(\xi)\quad(\xi 介于 x 与节点之间)$$

应用到 Adams 外插法：这里 $g(t)=f(t,u(t))$ ，$L_{n,k}(t)$ 是插值多项式，$\displaystyle r_{n,k}(t)=g(t)-L_{n,k}(t)$ 为插值余项

根据插值余项公式

$$\displaystyle r_{n,k}(t)=\frac{(t - t_{n - k})(t - t_{n - k + 1})\cdots(t - t_n)}{(k + 1)!}g^{(k + 1)}(\xi)$$

$\xi$ 与 $t$ 有关，局部截断误差就是对 $\displaystyle r_{n,k}(t)$ 在积分区间上积分产生的误差

将 $t = t_n+\tau h$ 代入插值余项表达式 $\displaystyle r_{n,k}(t)$ 。根据插值余项公式 $\displaystyle r_{n,k}(t)=\frac{(t - t_{n - k})(t - t_{n - k + 1})\cdots(t - t_n)}{(k + 1)!}u^{(k + 2)}(\xi)$ ，当 $t = t_n+\tau h$ 时，$\displaystyle t - t_{n - i}=\tau h + ih = (\tau + i)h$ （$i = 0,1,\cdots,k$ ），

代入后经过整理得到

$$\displaystyle r_{n,k}(t_n+\tau h)=\frac{(\tau + k)(\tau + k - 1)\cdots(\tau + 1)\tau h^{k + 1}}{(k + 1)!}u^{(k + 2)}(\xi)$$

进一步变形为

$$\displaystyle r_{n,k}(t_n+\tau h)=(-1)^{k + 1}\binom{-\tau}{k + 1}h^{k + 1}u^{(k + 2)}(\xi)$$

积分变换：已知 $\displaystyle R_{n,k}=\int_{t_n}^{t_{n + 1}}r_{n,k}(t)dt$ ，通过换元 $t = t_n+\tau h$ ，$dt = h d\tau$ ，积分区间从 $[t_n,t_{n + 1}]$ 变为 $[0,1]$

则

$$R_{n,k}=\int_{0}^{1}r_{n,k}(t_n+\tau h)\cdot h d\tau$$

代入化简：将 $\displaystyle r_{n,k}(t_n+\tau h)=(-1)^{k + 1}\binom{-\tau}{k + 1}h^{k + 1}u^{(k + 2)}(\xi)$ 代入上式

得到

$$R_{n,k}=\int_{0}^{1}(-1)^{k + 1}\binom{-\tau}{k + 1}u^{(k + 2)}(\xi)d\tau h^{k + 2}$$

结果分析：由于 $\displaystyle u^{(k + 2)}(\xi)$ 在积分区间内有界（在一定条件下 ），令 $\displaystyle \alpha_{k + 1}=\int_{0}^{1}(-1)^{k + 1}\binom{-\tau}{k + 1}d\tau$

则

$$\displaystyle R_{n,k}=\alpha_{k + 1}u^{(k + 2)}(\overline{\xi})h^{k + 2}$$

可知 $\displaystyle R_{n,k}=O(h^{k + 2})$ ，这表明 Adams 外插法的局部截断误差阶为 $h^{k + 2}$ ，即随着步长 $h$ 减小，误差以 $\displaystyle h^{k + 2}$ 的速度趋近于 0

Adams 内插法

Adams 内插法为隐格式

局部截断误差比外插法高一阶，为 $O(h^{k + 3})$

$k = 0, 1, 2, 3$的内插公式：

$k = 0$：$\displaystyle u_{n + 1} = u_n + hf_{n + 1}$ → 隐式 Euler

$k = 1$：$\displaystyle u_{n + 1} = u_n + \frac{h}{2}(f_{n + 1} + f_n)$ → 改进 Euler

$k = 2$：$\displaystyle u_{n + 1} = u_n + \frac{h}{12}(5f_{n + 1} + 8f_n - f_{n - 1})$

$k = 3$：$\displaystyle u_{n + 1} = u_n + \frac{h}{24}(9f_{n + 1} + 19f_n - 5f_{n - 1} + f_{n - 2})$

多步法一般形式：

$$\sum_{j = 0}^{k} \alpha_j u_{n + j} = h\sum_{j = 0}^{k} \beta_j f_{n + j}$$

由 $\displaystyle u_n, u_{n + 1}, \cdots, u_{n + k - 1}$ → $u_{n + k}$ （共 $k$ 个，即$k$步法）

Adams 方法：

$$\alpha_k = 1，\alpha_{k - 1} = -1，\displaystyle \alpha_{k - 2} = \cdots = \alpha_0 = 0$$

注：$\alpha_k$ 总可取成 1，显方法 $\beta_k = 0$，隐方法 $\beta_k \neq 0$

实际计算时，初始 $u_0, u_1, \cdots, u_{k - 1}$，可由单步法算

待定系数法

将真解 $u(t)$ 代入格式：$\displaystyle L[u(t_n); h] = \sum_{j = 0}^{k} [\alpha_j u(t_{n + j}) - h\beta_j f(t_{n + j}, u(t_{n + j}))]$

为简便且讨论一般，在 $t_n$ 处作 Taylor 展开：

$$\begin{align} &u(t_{n + j}) = u(t_n + jh) = u(t_n) + u'(t_n)jh + \frac{u''(t_n)}{2!}(jh)^2 + \cdots + \frac{u^{(q)}(t_n)}{q!}(jh)^q + \cdots \\ &f(t_{n + j}, u(t_{n + j})) = u'(t_{n + j}) = u'(t_n + jh) = u'(t_n) + u''(t_n)jh + \frac{u'''(t_n)}{2!}(jh)^2 + \cdots + \frac{u^{(q + 1)}(t_n)}{q!}(jh)^q + \cdots \end{align}$$

则：

$$\begin{align*} &\sum_{j = 0}^{k} [\alpha_j u(t_{n + j}) - h\beta_j f(t_{n + j}, u(t_{n + j}))] \\ &=\sum_{j = 0}^{k} [\alpha_j (u(t_n) + u'(t_n)jh + \frac{u''(t_n)}{2!}(jh)^2 + \cdots + \frac{u^{(q)}(t_n)}{q!}(jh)^q + \cdots) - h\beta_j (u'(t_n) + u''(t_n)jh + \frac{u'''(t_n)}{2!}(jh)^2 + \cdots + \frac{u^{(q + 1)}(t_n)}{q!}(jh)^q)] \\ &=\sum_{j = 0}^{k} \alpha_j u(t_n) + \sum_{j = 0}^{k} (j\alpha_j - \beta_j)u'(t_n) \cdot h + \cdots + \sum_{j = 0}^{k} (\alpha_j \frac{j^q}{q!} - \beta_j \frac{j^{q - 1}}{(q - 1)!}) \cdot h^q \\ &= C_0 u(t_n) + C_1 u'(t_n) + \cdots + C_q u^{(q)}(t_n) \cdot h^q + \cdots \end{align*}$$

其中：

$$\begin{equation} \begin{cases} C_0 = \alpha_0 + \alpha_1 + \cdots + \alpha_k \\ C_1 = \alpha_1 + 2\alpha_2 + \cdots + k\alpha_k - (\beta_0 + \beta_1 + \cdots + \beta_k) \\ \quad \quad \vdots \\ C_q = \frac{1}{q!}(\alpha_1 + 2^q\alpha_2 + \cdots + k^q\alpha_k) - \frac{1}{(q - 1)!}(\beta_1 + 2^{q - 1}\beta_2 + \cdots + k^{q - 1}\beta_k) \end{cases} \end{equation}$$

若 $u(t)$ 有 $p + 2$ 次连续导数，则可选 $\alpha_j, \beta_j$

使得：

$$C_0 = C_1 = \cdots = C_p = 0，C_{p + 1} \neq 0$$

此时：

$$L[u(t_n); h] = C_{p + 1} u^{(p + 1)}(t_n)h^{p + 1} + O(h^{p + 2})$$

得到了 $p$ 阶 $k$ 步法：

$$\sum_{j = 0}^{k} \alpha_j u_{n + j} = h\sum_{j = 0}^{k} \beta_j f_{n + j}$$

局部截断误差为 $O(h^{p + 1})$
整体误差为 $O(h^p)$

相容性、稳定性和误差估计

k 步法：

$$\sum_{j = 0}^{k} \alpha_j u_{n + j} = h\sum_{j = 0}^{k} \beta_j f_{n + j}$$

相容性

相容性：差分算子逼近微分算子

以 Euler 法为例：

$$u_{n + 1} = u_n + hf(t_n, u_n)，变形为 \frac{u_{n + 1} - u_n}{h} = f(t_n, u_n)$$

差分算子：

$$L[u(t_n); h] = u(t_{n + 1}) - u(t_n) - hf(t_n, u(t_n))$$

差商算子：

$$\frac{1}{h}L[u(t_n); h] = \frac{u(t_{n + 1}) - u(t_n)}{h} - f(t_n, u(t_n))$$

微分算子：

$$u'(t) - f(t, u(t))$$

Euler 法相容性：

$$\frac{1}{h}L[u(t_n); h] - [u'(t_n) - f(t_n, u(t_n))] \to 0，当 h \to 0$$

k 步法相容性：

$$\frac{1}{h} [\sum_{j = 0}^{k} \alpha_j u(t_{n + j}) - h\sum_{j = 0}^{k} \beta_j f(t_{n + j}, u(t_{n + j}))] - [u'(t_n) - f(t_n, u(t_n))] = o(1)，当 h \to 0$$

即

$$\sum_{j = 0}^{k} \alpha_j u(t_{n + j}) - h\sum_{j = 0}^{k} \beta_j f(t_{n + j}, u(t_{n + j})) = o(h) \ \Leftrightarrow (C_0 = 0, C_1 = 0)$$

其中

$$\sum_{j = 0}^{k} \alpha_j u(t_{n + j}) - h\sum_{j = 0}^{k} \beta_j f(t_{n + j}, u(t_{n + j})) = C_0u(t) + C_1u'(t)h + C_2u''(t)h^2 + \cdots$$

特征多项式

第一特征多项式：

$$\rho(\lambda) = \sum_{j = 0}^{k} \alpha_j \lambda^j$$

第二特征多项式：

$$\sigma(\lambda) = \sum_{j = 0}^{k} \beta_j \lambda^j$$

定理3.1

$$k步法相容 \Leftrightarrow C_0 = 0, C_1 = 0 \Leftrightarrow \rho(1) = 0, \rho'(1) = \sigma(1)$$

其中：

$$\begin{align*} &C_0 = \alpha_0 + \alpha_1 + \cdots + \alpha_k = \rho(1)\\ &C_1 = \alpha_1 + 2\alpha_2 + \cdots + k\alpha_k - (\beta_0 + \beta_1 + \cdots + \beta_k) = \rho'(1) - \sigma(1) \end{align*}$$

稳定性

稳定性：解对初值的连续依赖性

k 步法：由初值 $u_0, u_1, \cdots, u_{k - 1}$ 得到 $u_n$ 。若 $u_0, u_1, \cdots, u_{k - 1}$ 有误差，$u_n$ 的误差是否会无限增长？

定义 3.2：若 $\exists \tau > 0, h_0 > 0$ ，使得当网格步长 $h$ 满足 $0 < h < h_0$ 时，对 $k$ 步法的任意两个解 ${u_n}$， ${v_n}$

都有

$$\max_{n \leq \frac{T}{h}} |u_n - v_n| \leq \tau \max_{0 \leq j \leq k} |u_j - v_j|$$

则称 $k$ 步法稳定

具体方法的误差递推关系

显式 Euler：

$$|e_{n + 1}| \leq |e_n| + Lh|e_n| + R$$

此为误差递推关系

隐式 Euler：

$$\begin{align} &u\_{n + 1} = u\_n + hf(t\_{n + 1}, u\_{n + 1}) \\ &v\_{n + 1} = v\_n + hf(t\_{n + 1}, v\_{n + 1}) \\ &\rightarrow|e\_{n + 1}| \leq |e\_n| + Lh|e\_{n + 1}| \\ &\rightarrow|e\_{n + 1}| \leq \frac{1}{1 - Lh}|e\_n| \end{align}$$

k 步法：

已知：由 $u_0, u_1, \cdots, u_{k - 1}$ 得到 $u_n$ ，满足

$$\sum_{j = 0}^{k} \alpha_j u_{n + j} = h\sum_{j = 0}^{k} \beta_j f(t_{n + j}, u_{n + j})$$

由 $v_0, v_1, \cdots, v_{k - 1}$ 得到 $v_n$ ，满足

$$\sum_{j = 0}^{k} \alpha_j v_{n + j} = h\sum_{j = 0}^{k} \beta_j f(t_{n + j}, v_{n + j})$$

记 $e_n = u_n - v_n$ ，则

$$\sum_{j = 0}^{k} \alpha_j e_{n + j} = h\sum_{j = 0}^{k} \beta_j [f(t_{n + j}, u_{n + j}) - f(t_{n + j}, v_{n + j})]$$

记 $hb_n$ 为：

$$hb_n=h\sum_{j = 0}^{k} \beta_j [f(t_{n + j}, u_{n + j}) - f(t_{n + j}, v_{n + j})]$$

进一步得到

$$e_{n + k} = \frac{1}{\alpha_k}(-\alpha_{k - 1}e_{n + k - 1} - \cdots - \alpha_0 e_n + hb_n)$$

将 $k$ 个连续的 $e_i$ 构成向量形式：

$$\begin{bmatrix}e_{n + k}\\e_{n + k - 1}\\\vdots\\e_{n + 1}\end{bmatrix}= \begin{bmatrix} -\frac{\alpha_{k - 1}}{\alpha_k} &-\frac{\alpha_{k - 2}}{\alpha_k} &\cdots& -\frac{\alpha_0}{\alpha_k} \\ 1&0&\cdots&0\\ \vdots&\ddots&\ddots&\vdots\\0&\cdots&1&0\end{bmatrix} \begin{bmatrix}e_{n + k - 1}\\e_{n + k - 2}\\\vdots\\e_{n}\end{bmatrix}+\begin{bmatrix}h\frac{b_n}{\alpha_k}\\0\\\vdots\\0\end{bmatrix}$$

即

$$E_{n + 1} = CE_n + B_n$$

类似数值代数中迭代法

且

$$E_n = C^nE_0 + \sum_{l = 0}^{n - 1} C^l B_{n - 1 - l}$$

关于 $C^n$ 有界性分析

特征方程

计算 $\det(C - \lambda I_k)$ ，其结果为

$$(\alpha_k\lambda^k+\cdots+\alpha_1\lambda+\alpha_0)\cdot(-1)^k\alpha_k^{-1}$$

记

$$\rho(\lambda)=\alpha_k\lambda^k+\cdots+\alpha_1\lambda+\alpha_0$$

利用 Jordan 标准型计算

设 $\lambda_j$ 是 $C$ 的特征值，$\displaystyle x_j=(d_{k - 1},\cdots,d_0)^T$ 是相应特征向量 ，则有 $Cx_j=\lambda_jx_j$ ，由此推导出一系列等式：

$$\begin{align} -\alpha_k^{-1}&(\alpha_{k - 1}d_{k - 1}+\cdots+\alpha_1d_1+\alpha_0d_0)=\lambda_jd_{k - 1} \\ &d_{k - 2}=\lambda_jd_{k - 3} ，\cdots ，d_1=\lambda_jd_0 \end{align}$$

进而得到

$$d_{k - 1}=\lambda_j^{k - 1}d_0\\ x_j = d_0(\lambda_j^{k - 1},\cdots,\lambda_j,1)$$

因而 $C$ 的每个特征值的特征子空间维数都是 1，即每个特征值几何重数为 1，每个特征值的 Jordan 块只有 1 块。记 $C$ 的 Jordan 标准型为 $J$

则

$$C = PJP^{-1}$$

Jordan 块形式：$J$ 的 Jordan 块要么是 $\lambda$ ，要么是 $J_r(\lambda)$，$r$ 是 $\lambda$ 的重数

$$J_r(\lambda)=\begin{bmatrix} \lambda & 1 & \cdots & 0 \\ 0 & \lambda & \ddots & \vdots \\ \vdots & \vdots & \ddots & 1 \\ 0 & 0 & \cdots & \lambda \end{bmatrix}$$

$\displaystyle C^n = PJ^nP^{-1}$ ，且 $J^n$ 也是分块矩阵，每一块形式为 $\lambda^n$ 或 $J_r^n$

$$J_r^n(\lambda)=\begin{bmatrix} \lambda & c_n^1\lambda^{n-1}-\cdots c_n^{r-1}\lambda^{n-r+1} & \cdots & 0 \\ 0 & \lambda & \ddots & \vdots \\ \vdots & \vdots & \ddots & c_n^1\lambda^{n-1} \\ 0 & 0 & \cdots & \lambda \end{bmatrix}$$

对于 $J_r=\lambda I + S$ ，（其中 $S$ 是幂零矩阵）

$$S=\begin{bmatrix} 0& 1 & \cdots & 0 \\ 0 & 0 & \ddots & \vdots \\ \vdots & \vdots & \ddots & 1 \\ 0 & 0 & \cdots & 0 \end{bmatrix} \quad S^2=\begin{bmatrix} 0 & 0 &1 & \cdots & 0 \\ 0 & 0& \ddots & \ddots &\vdots\\ 0 & 0 &\ddots &\ddots &1\\ \vdots & \vdots & \ddots & 0 &0\\ 0 & 0 & \cdots & 0&0 \end{bmatrix} \quad S^r=0$$

利用二项式定理计算

$$\begin{align*} J_r^n &= (\lambda I+S)^n \\ &=\sum_{j = 0}^{n}C_{n}^{j}\lambda^{n - j}S^j \\ &=\sum_{j = 0}^{r - 1}C_{n}^{j}\lambda^{n - j}S^j \end{align*}$$

这里 $\displaystyle C_{n}^{j}=\frac{n!}{j!(n - j)!}$

引理 1.1

$C^n$有界 $\Leftrightarrow |\lambda_j| \leq 1$ ，且若 $|\lambda_j| = 1$ ，则 $\lambda_j$ 是单根

$C^n$有极限 $\Leftrightarrow |\lambda_j| \leq 1$ ，且若 $|\lambda_j| = 1$ ，则 $\lambda_j$ 不仅是单根且 $\lambda_j = 1$

$C^n \to 0 \Leftrightarrow |\lambda_j| < 1$

定理：线性多步法稳定的充要条件

线性 $k$ 步法稳定 $\Leftrightarrow$ 第一特征多项式 $\rho(\lambda)$ 的根满足 $|\lambda| \leq 1$ ，且位于单位圆周上的根都是单根（即 $\rho(\lambda)$ 满足根条件

证明

必要性证明

将 $k$ 步法用于方程 $u’ = 0$ ，此时 $f(t, u) \equiv 0$ ，则 $\displaystyle \sum_{j = 0}^{k} \alpha_j u_{n + j} = 0$ ，${u_n = 0}$ 是上述方程的平凡解。

因为 $k$ 步法稳定，所以 $\exists \tau, h_0$ ，使得当 $0 < h < h_0$ 时，$\displaystyle \max_{n \leq \frac{T}{h}} |u_n| \leq \tau \max_{0 \leq j \leq k} |u_j|$ ，即 ${u_n}$ 关于 $n, h$ 一致有界。

又因为 $u_n = C^n u_0$ （这里 $\displaystyle u_0 = [u_{k - 1}, \cdots, u_0]^T$ ，$\displaystyle u_n = [u_{n + k - 1}, \cdots, u_n]^T$ ），所以 ${C^n}$ 一致有界，从而 $\rho(\lambda)$ 满足根条件。

充分性证明

设 ${u_n}, {v_n}$ 是 $k$ 步法的两个解，令 $e_n = u_n - v_n$ ，则

$$\displaystyle \sum_{j = 0}^{k} \alpha_j e_{n + j} = h\sum_{j = 0}^{k} \beta_j [f(t_{n + j}, u_{n + j}) - f(t_{n + j}, v_{n + j})] = hb_n$$

设 $B = \max{|\beta_0|, \cdots, |\beta_k|}$ ，由 $f$ 关于 $u$ 满足 Lipschitz 条件，可得

$$|b_n| \leq B\cdot L \sum_{j = 0}^{k} |e_{n + j}|$$

记 $\displaystyle E_n = (e_{n + k - 1}, \cdots, e_n)^T$ ，$\displaystyle B_n = (h\alpha_k^{-1}b_n, 0, \cdots, 0)^T$ ，则 $E_{n + 1} = CE_n + B_n$ ，进而有

$$E_n = C^n E_0 + \sum_{l = 0}^{n - 1} C^l B_{n - 1 - l}$$

由于 $\rho(\lambda)$ 满足根条件，所以 $|C^n| \leq M, \forall n$ ，且

$$\|B_n| \leq h|\alpha_k^{-1}| |b_n| \leq h|\alpha_k^{-1}| \cdot B\cdot L \sum{j = 0}^{k} |e{n + j}|$$

经过推导可得

$$\begin{align} \|E_n\| &\leq M\|E_0\| + \sum_{l = 0}^{n - 1} M\cdot h|\alpha_k^{-1}| \cdot B\cdot L \sum_{j = 0}^{k} |e_{n - 1 - l + j}| \\&\leq M\|E_0\| + M_1 \sum_{l = 0}^{n - 1} \sum_{j = 0}^{k} |e_{n - 1 - l + j}| \end{align}$$

对 $\displaystyle \sum_{l = 0}^{n - 1} \sum_{j = 0}^{k} |e_{n - l - 1 + j}|$ 进行放缩：

$$\sum_{l = 0}^{n - 1} \sum_{j = 0}^{k} |e_{n - l - 1 + j}|=\sum_{l = 0}^{n - 1} |e_{n + k - 1 - l}|+\sum_{l = 0}^{n - 1} \sum_{j = 0}^{k} |e_{n - l - 1 + j}|$$

利用向量范数性质，有

$$\sum_{l = 0}^{n - 1} |e_{n + k - 1 - l}|+\sum_{l = 0}^{n - 1} \sum_{j = 0}^{k} |e_{n - l - 1 + j}| \leq \sum_{l = 0}^{n - 1} \|E_{n - l - 1}\|+\sum_{l = 0}^{n - 1} \sqrt{k} \|E_{n - l - 1}\|$$

进一步得到

$$\sum_{l = 0}^{n - 1} \|E_{n - l - 1}\|+\sum_{l = 0}^{n - 1} \sqrt{k} \|E_{n - l - 1}\|=\|E_n\|+(1 + \sqrt{k})\sum_{i = 0}^{n - 1} \|E_i\|$$

结合前面 $|E_n|$ 的表达式

$$\|E_n\| \leq M\|E_0\| + M B L|\alpha_k^{-1}| h \|E_n\|+M B L|\alpha_k^{-1}| h\cdot(\sqrt{k}+1)\sum_{l = 0}^{n - 1} \|E_l\|$$

在 $\displaystyle MBL|\alpha_k^{-1}| h_0 < 1$ 条件下，令 $\displaystyle K_1=\frac{M}{1 - M B L|\alpha_k^{-1}| h_0}$ ，$\displaystyle K_2 = K_1 M B L|\alpha_k^{-1}|(\sqrt{k}+1)$ ，可得

$$\|E_n\| \leq M\|E_0\| + M B L|\alpha_k^{-1}| h \|E_n\|+M B L|\alpha_k^{-1}| h\cdot(\sqrt{k}+1)\sum_{l = 0}^{n - 1} \|E_l\|$$

再由 Gronwall 不等式，得到

$$\|E_n\| \leq e^{K_2 T}(K_1 + K_2 h_0)\|E_0\|$$

进而有

$$\max_{n \leq \frac{T}{h}} |u_n - v_n| \leq \overline{c} \max_{0 \leq j \leq k} |u_j - v_j|$$

说明方法稳定

应用：Euler 法稳定性、Adams 外插和内插法稳定性

收敛性与误差估计

定义与设定：设真解 $u(t_n)$ 和数值解 $u_n$ 满足 ：

$$\begin{align} \sum_{j = 0}^{k} \alpha_j u(t_{n + j}) &= h\sum_{j = 0}^{k} \beta_j f(t_{n + j}, u(t_{n + j})) + R_n \\ \sum_{j = 0}^{k} \alpha_j u_{n + j} &= h\sum_{j = 0}^{k} \beta_j f(t_{n + j}, u_{n + j}) \end{align}$$

令 $\displaystyle e_n = u(t_n) - u_n$

则

$$\begin{align} \sum_{j = 0}^{k} \alpha_j e_{n + j} &= h\sum_{j = 0}^{k} \beta_j [f(t_{n + j}, u(t_{n + j})) - f(t_{n + j}, u_{n + j})]+R_n \\&= hb_n + R_n \end{align}$$

设 $\displaystyle R_n = C_{p + 1}h^{p + 1}$ ，记 $\displaystyle E_{n + 1} = CE_n + B^{\ast}_n$ ，其中 $\displaystyle B^{\ast}_n=(h\alpha_k^{-1}b^{\ast}_n, 0, \cdots, 0)^T$ ，且

$$\displaystyle \|B_n^*\| \leq h B L|\alpha_k^{-1}| \sum_{j = 0}^{k} |e_{n + j}|+M_{p + 1}h^{p + 1}$$

定理

相容性 + 稳定性 $\Rightarrow$ 收敛。且若 $E_0 = 0$ ，则 $|E_n| \leq C h^p$ 。即当数值方法满足相容性（差分算子能逼近微分算子）和稳定性（解对初值连续依赖）时，随着步长 $h$ 趋于 0 ，数值解会收敛到真解 ，且在初始误差为 0 时，误差阶为 $O(h^p)$