ODE数值解

转载自 李森科在zhihu 链接

第一部分常微分方程的数值解法

常微分方程数值解法主要分为两大部分，初值问题与边值问题的数值方法。本部分只讨论初值问题。

第一章常微分方程初值问题

本章讨论一阶常微分方程初值问题 $\begin{cases} \frac{du(t)}{dt}=f(t,u(t))\\ u(t_0)=u_0\\ \end{cases}\ (1.0.1)$ 的数值求解，其中 $f$ 为 $t$ 和 $u$ 的已知函数， $u_0$ 为给定的初值。在基本假设下，初值问题 $(1.0.1)$ 在区间 $[t_0,T]$ 上有唯一解 $u(t)$ ，并且 $u(t)$ 是连续可微的。

数值解法是一种离散化方法，利用这种方法，可以在一系列离散点 $t_1,...,t_N$ 上求出未知函数 $u(t)$ 的近似值 $u_1,...,u_N$ ，这个近似值称为初值问题的一个数值解。

自变量 $t$ 的离散值 $t_1,...,t_N$ 是事先取定的，称之为节点，通常取成等距的，即 $t_1=t_0+h,...,t_N=t_0+Nh$ ，其中 $h>0$ 称为步长。

1.1 基本概念 Euler法与梯形法

1.1.1 Euler法介绍

本节通过Euler法及梯形法的讨论，说明常微分方程数值解法中的一些基本概念与主要研究的问题。

首先得到初值问题 $(1.0.1)$ 等价的积分形式 $u(t+h)=u(t)+\int^{t+h}_{t}f(\tau,u(\tau))d\tau \ (1.1.2)$

设 $t_1=t_0+h$ ，使用左矩形公式计算一步 $u_1=u_0+hf(t_0,u_0)$ ，类似地利用 $u_1$ 和 $f(t_1,u_1)$ 计算出 $u(t_2)$ 的近似值，那么可以得到Euler法的计算公式 $\begin{cases}u_0=u_0，\\ u_{n+1}=u_n+hf(t_n,u_n)\ (1.1.1),n\geq0\\ \end{cases}\$

Euler法的几何意义是十分清楚的：在这种方法中，实际上是用一条过 $(t_0,u_0)$ 的折现来近似替代过 $(t_0,u_0)$ 的积分曲线，因此这种方法又称为折线法。

1.1.2 收敛性，截断误差，稳定性问题的简单介绍

为考察由Euler法提供的数值解是否具有实用价值，首先应该知道，当步长 $h$ 取得充分小时，所得数值解 $u_n$ 能否足够精确地逼近初值问题 $(1.0.1)$ 的真解 $u(t_n)$ ，这就是收敛性问题。

其次，还必须估计数值解与真解之间的误差，以便于在实际计算中根据精度要求确定计算方案。在Euler法中产生的误差称为截断误差。

其次，计算过程中还会由于数值误差的舍入产生舍入误差。由于Euler法是一种步进方法，只有当最初产生的误差在以后各步的计算中不会无限制扩大时，即只有当 $(1.1.1)$ 的解对初值具有某种连续相依性质时，方法才具有实用价值。这种性质称为稳定性问题。

1.1.3 Euler法的截断误差估计

首先讨论Euler法的截断误差估计及收敛性问题。初值问题 $(1.0.1)$ 可以写成等价的积分形式 $u(t+h)=u(t)+\int^{t+h}_{t}f(\tau,u(\tau))d\tau\ (1.1.2)$

在上式中，令 $t=t_n$ ，并用左矩形公式计算右端积分，有

$u(t_{n+1})=u(t_n)+\int^{t_{n+1}}_{t_n}f(\tau,u(\tau))d\tau$ $=u(t_n)+\int^{t_{n+1}}_{t_n}[f(\tau,u(\tau)-f(t_n,u(t_n))]d\tau+\int^{t_{n+1}}_{t_n}f(t_n,u(t_n))d\tau$ $=u(t_n)+hf(t_n,u(t_n))+\int^{t_{n+1}}_{t_n}[f(\tau,u(\tau)-f(t_n,u(t_n))]d\tau$ $=u(t_n)+hf(t_n,u(t_n))+R_n\ (1.1.3)$

$R_n$ 称为Euler法的局部截断误差，它表示当 $u_n=u(t_n)$ 为精确值时，利用公式 $(1.1.1)$ 计算 $u(t_n+h)$ 的误差。在逐步计算的过程中，局部截断误差会传播和积累，因此还必须对这种误差的积累与传播作出估计。

设 $u_n$ 是在无舍入误差情况下用 $(1.1.1)$ 计算出的数值解，而 $u(t_n)$ 是真解， $\epsilon_n=u(t_n)-u_n$ 为Euler法 $(1.1.1)$ 的整体截断误差。

为估计 $\epsilon_n$ ，先给出 $R_n$ 的上界。从式 $(1.1.3)$ 知，

$R_n=\int^{t_{n+1}}_{t_n}[f(\tau,u(\tau)-f(t_n,u(t_n))]d\tau$ $=\int^{t_{n+1}}_{t_n}[f(\tau,u(\tau)-f(t_n,u(\tau))]d\tau+\int^{t_{n+1}}_{t_n}[f(t_n,u(\tau)-f(t_n,u(t_n))]d\tau$

若函数 $f(t,u)$ 关于 $t$ 也满足Lipschitz条件，以 $K$ 记相应的Lipschitz常数，则由上式推出 $|R_n|\leq\int^{t_{n+1}}_{t_n}\Big|f(\tau,u(\tau)-f(t_n,u(\tau))\Big|d\tau+\int^{t_{n+1}}_{t_n}\Big|f(t_n,u(\tau)-f(t_n,u(t_n))\Big|d\tau$ $\leq K\int^{t_{n+1}}_{t_n}|\tau-t_n|d\tau+L\int^{t_{n+1}}_{t_n}|u(\tau)-u(t_n)|d\tau$ $\leq K\int^{h}_{0}|u|du+hL\int^{t_{n+1}}_{t_n}|\frac{u(\tau)-u(t_n)}{h}|d\tau$ $\leq K\int^{h}_{0}udu+hL\int^{t_{n+1}}_{t_n}Md\tau=\frac{1}{2}h^2K+h^2LM\leq h^2(K+LM)$

其中 $M=\max_{t_0\leq t\leq T}|u'(t)|=\max_{t_0\leq t\leq T}|f(t,u(t))|$ ，于是 $|R_n|\leq R=h^2(K+LM)$

用式 $(1.1.3)$ 减去式 $(1.1.1)$ ，得到整体截断误差所满足的方程

$\epsilon_{n+1}=\epsilon_{n}+h[f(t_n,u(t_n))-f(t_n,u_n)]+R_n$ ，进一步可以得到

$|\epsilon_{n+1}|\leq|\epsilon_{n}|+h|f(t_n,u(t_n))-f(t_n,u_n)|+R$ $\leq|\epsilon_n|+hL|\epsilon_n|+R$ $=(1+hL)|\epsilon_n|+R\leq(1+(T-t_0)L)|\xi_n|+R$

这是一个递推公式，我们需要得到对任意 $n$ 均成立的上界。下面的引理能够做到这一点。

引理1.1.1 设数列 $\xi_i$ 满足不等式 $|\xi_{i+1}|\leq(1+\delta)|\xi_i|+B,$ ，其中 $\ \delta>0,\ B\geq0,\ i=0,1,2,...$ ，则有 $|\xi_i|\leq e^{i\delta}|\xi_0|+\frac{e^{i\delta}-1}{\delta}B,\ i>1$ 。

由此我们可以得到整体截断误差的估计 $|\epsilon_n|\leq e^{L(T-t_0)}|\epsilon_0|+\frac{h}{2}(e^{L(T-t_0)}-1)(\frac{LM+K}{L})\ (1.1.5)$

从 $u_0=u(t_0)$ 及式 $(1.1.5)$ ，又知 $\epsilon_n=O(h)，h\rightarrow0$ ，即Euler法的整体阶段误差与 $h$ 同阶，此时称Euler法为一阶方法。它的直观意义是，当 $u(t)$ 是一次多项式时，Euler法是精确的。

1.1.4 Euler法的稳定性

这里考虑仅初值 $u_0$ 有误差，而其他计算步骤无误差的简单情况。

…

上述结论所表述的Euler法的性质，及所谓关于初值的稳定性：当初始误差充分小时，以后各步的误差也可充分小。

1.1.5 梯形法

由式 $(1.1.2)$ ，取 $t=t_n$ ，并用梯形求积公式计算右端积分，舍去每一步梯形求积公式的余项，可以得到梯形法的计算公式 $\begin{cases} u_0=u_0\\ u_{n+1}=u_n+\frac{h}{2}[f_n+f_{n+1}],n\geq0\\ \end{cases}$

对于梯形法，可类似地建立其截断误差估计式及收敛性，且 $|u(t_n)-u_n|=O(h^2),h\rightarrow0$ ，此时称梯形法为二阶方法。

这里要注意梯形法与Euler法的一个重要区别，即式 $(1.1.8)$ 只给出一个关于 $u_{n+1}$ 的（非线性）方程，而Euler法给出的是 $u_{n+1}$ 的显式表达式，因此Euler法是一个显式方法，而梯形法是隐式方法。

当梯形法求 $u_{n+1}$ 时，可利用求解非线性方程的各种方法，最简单的是下述迭代方法： $u_{n+1}^{(m+1)}=u_n+\frac{h}{2}[f_n+f(t_{n+1},u_{n+1}^{(m)})],m\geq0$ 。

实际计算时，初始近似 $u_{n+1}^{(0)}$ 可用Euler法求出，然后再使用迭代法，这便是预测校正格式（又称为改进Euler法）：

$\begin{cases} u_{n+1}^{(0)}=u_n+hf(t_{n},u_n)\\ u_{n+1}^{(m+1)}=u_n+\frac{h}{2}[f_n+f(t_{n+1},u_{n+1}^{(m)})],m\geq0\\ \end{cases}$

1.1.6 总结

通过上述两个方法的讨论可知，数值方法的研究有如下基本问题：计算公式的构造，方法的收敛性，截断误差估计，稳定性及方法的实际应用等。

1.2 Runge-Kutta方法及一般单步方法

前节介绍的两个方法精度是很低的。本节从提高局部截断误差阶入手，来构造具有较高精度的方法。

1.2.1 Taylor级数法

设初值问题 $(1.0.1)$ 的解 $u(t)$ 具有 $q+1$ 次连续导数。考虑 $u(t)$ 在 $t=t_0$ 处的Taylor展开式，可得局部截断误差 $u(t_0+h)-u_1=O(h^{q+1})$ 。

…

这种方法称为Taylor级数法。虽然它的局部截断误差的阶可任意高，但各阶导数 $u^{(i)}(t)$ 的计算量非常大，不实用。现在介绍一类既可以达到较高精度，又可避免高阶导数计算的方法——Runge-Kutta方法。

1.2.2 RK方法

根据公式 $(1.1.2)$ 及积分中值定理 $\int^{t+h}_{t}f(\tau,u(\tau))d\tau =hf(t_n+\theta h,u(t_n+\theta h))$ ，可得 $u(t_n+h)=u(t_n)+hf(t_n+\theta h,u(t_n+\theta h)) \ (1.2.3)$ 。

但 $f(t_n+\theta h,u(t_n+\theta h))$ 是无法计算的。我们用 $f$ 在位于 $[t_n,t_n+h]$ 上的若干个点（例如 $s$ 个点）值的线性组合来近似它，并使之有尽可能高的精度。具体来讲，就是用下列公式代替式 $(1.2.3)$ $\begin{equation} u_{n+1}=u(t_n)+h\sum^{s}_{i=1}W_iK_i\ (1.2.4)\\ \begin{cases} K_1=f(t_n,u(t_n))，\\ K_i=f(t_n+h\alpha_i,u(t_n)+h\sum^{i-1}_{j=1}\beta_{ij}K_i)\ \\ \end{cases}(1.2.5)\\ \end{equation}\$

现在去选定系数 $\alpha_i,\beta_{ij}(i\geq2)$ ，以及 $W_i(i\geq1)$ ，使RK方法的整体截断误差具有尽可能高的阶 $q$ ，这个阶数 $q$ 当然与级数 $s$ 有关，可记为 $q(s)$ 。

构造RK方法的基本步骤是，选定上述系数，使 $u(t_n+h)-u(t_n)$ 与 $h\sum^{s}_{i=1}W_iK_i$ 两者的Taylor展开式中含 $h^r(r=1,...,q)$ 的项的系数对尽可能大的 $q$ 相等。

书上考虑的是 $s=4$ 的情况，在这里我们只考虑 $s=2$ 的情况。

-————————————————————————————————-

由一元函数Taylor展开可得

$u(t_n+h)-u(t_n)=u'(t_n)h+u''(t_n)\frac{h^2}{2}+O(h^3)$ $=f(t_n,u(t_n))h+\Big[f_t\big(t_n,u(t_n)\big)+f_u\big(t_n,u(t_n)\big)f\big(t_n,u(t_n)\big)\Big]\frac{h^2}{2}+O(h^3)$ $=hf(t_n,u(t_n))+\frac{h^2}{2}f_t\big(t_n,u(t_n)\big)+\frac{h^2}{2}f_u\big(t_n,u(t_n)\big)f\big(t_n,u(t_n)\big)+O(h^3)$

由二元函数的Taylor展开可得

$K_2=f(t_n+a_2h,u(t_n)+\beta_{21}K_1h)$

$=f(t_n,u(t_n))+f_t\big(t_n,u(t_n)\big)\alpha_2h+f_u(t_n,u(t_n))\beta_{21}K_1h$ $+O(h^2)$

$=f\big(t_n,u(t_n)\big)+f_t(t_n,u(t_n))\alpha_2h+f\big(t_n,u(t_n)\big)f_u(t_n,u(t_n))\beta_{21}h$ $+O(h^2)\ (RK.2)$

将 $(RK.2)$ 代入 $(1.2.4)$ 可得

$u_{n+1}-u(t_n)=h(W_1K_1+W_2K_2)=hW_1K_1+hW_2K_2$ $=hW_1f\big(t_n,u(t_n)\big)+hW_2\Big(f\big(t_n,u(t_n)\big)+f_t(t_n,u(t_n))\alpha_2h$ $+f\big(t_n,u(t_n)\big)f_u(t_n,u(t_n))\beta_{21}h+O(h^2)\Big)$ $=hW_1f\big(t_n,u(t_n)\big)+hW_2f\big(t_n,u(t_n)\big)+h^2W_2\alpha_2f_t(t_n,u(t_n))$ $+h^2W_2\beta_{12}f\big(t_n,u(t_n)\big)f_u(t_n,u(t_n))+O(h^3)\$

$=(W_1+W_2)hf\big(t_n,u(t_n)\big)+W_2\alpha_2h^2f_t(t_n,u(t_n))$ $+W_2\beta_{12}h^2f\big(t_n,u(t_n)\big)f_u(t_n,u(t_n))+O(h^3)$

比较 $(RK.1)$ 与 $(RK.3)$ ，使 $h^r(r=1,2)$ 的系数相等，得

$\begin{cases} W_1+W_2=1\\ W_2\alpha_2=\frac{1}{2}\\ W_2\beta_{21}=\frac{1}{2}\\ \end{cases}$ ，即 $\begin{cases} W_1+W_2=1\\ W_2\alpha_2=\frac{1}{2}\\ \beta_{21}=\alpha_2\\ \end{cases}$

取不同的 $\alpha_2$ ，可以得到不同的2级RK公式。

（以上 $s=2$ 情况的推导与书本无关）

-————————————————————————————————-

…

含 $h^r(r=1,2,3,4)$ 项的系数对应相等，且一般不能再使含 $h^5$ 项的系数相等，故所得方法的局部截断误差为 $O(h^5)$ ，即方法是四阶的。

…

从这些公式的讨论得知，当 $s\leq4$ 时， $s$ 级RK方法的阶 $q(s)=s$ ；当 $s=5,6,7$ ， $q(s)=s-1$ ；当 $s=8,9$ ， $q(s)=s-2$ ； $q(10)\leq8$ 。可见，对于四级以上的公式，计算函数值的计算量增加较快，而精度即阶数提高较慢，因此在实际问题中比较常用的是四级RK公式。

1.2.3 单步法的相容性与收敛性

无论是Euler法还是各级RK法，他们都是在已知 $u_n$ 的条件下算出 $u_{n+1}$ ，因此称为显式单步方法，简称为单步方法。这类方法的一般形式为 $u_{n+1}=u_n+h\varphi(t_n,u_n,h)\ (1.2.28)$ ，其中 $\varphi(t,u,h)$ 称为方法 $(1.2.28)$ 的增量函数，它是依赖于 $u$ 及 $f$ 的非线性函数（此处为简单计，略写了 $\varphi$ 对 $f$ 的依赖关系）。

在前面的讨论中，我们一再提到了每个具体方法的阶，现对单步方法 $(1.2.8)$ 给出阶的一般定义。

定义1.2.1 单步方法 $(1.2.28)$ 称为是 $q$ 阶的，如果对于真解 $u(t)$ ， $q$ 是使关系式 $u(t+h)-u(t)=h\varphi(t,u(t),h)+O(h^{q+1})$ 成立的最大整数。

…

1.3 线性多步方法

在前面所讨论的单步方法中，为提高截断误差的阶，每个时间布必须增加计算右端 $f(t,u)$ 的次数。当 $f$ 的结构比较复杂时，计算量会比较大。现在指出另一个提高截断误差阶的办法：构造这样的方法，在计算公式中 $u_{n+1}-u_{n}$ 不仅依赖于 $u_{n}$ ，且也直接依赖于 $u_{n-1},n_{n-2}$ 等已算出的值，但每进一步，只计算一次 $f(t,u)$ 的值。这样的方法称为多步方法。记 $f_n=f(t_n,u_n)$ ，若函数值 $f_n,f_{n-1}$ 等以线性组合的形式出现于公式中，则称方法为线性多步方法。

1.3.1 数值积分法

仍然考虑初值问题 $(1.0.1)$ 。易知 $u(t_{n+k})-u(t_{n-j})=\int^{t_{n+k}}_{t_{n-j}}f(t,u(t))dt\ (1.3.1)$

我们用被积函数 $f(t,u(t))$ 的 $q$ 次Lagrange插值多项式 $P_q(t)=\sum^{q}_{i=0}f(t_{n-i},u(t_{n-i}))L_i(t)$ 来近似替代 $(1.3.1)$ 中的被积函数，其中 $L_i(t)=\prod_{l=0,l\neq i}^{q}\frac{t-t_{n-l}}{t_{n-i}-t_{n-l}}=\prod_{l=0,l\neq i}^{q}\frac{t-t_n+t_n-t_{n-l}}{t_{n-i}-t_{n-l}}$ $=\prod_{l=0,l\neq i}^{q}\frac{t-t_n+lh}{(-i+l)h}$ $=\prod_{l=0,l\neq i}^{q}\frac{\frac{t-t_n}{h}+l}{(-i+l)}$ ，而有下面的恒等式成立 $\int^{t_{n+k}}_{t_{n-j}}L_i(t)dt=\prod_{l=0,l\neq i}^{q}\int^{t_{n}+kh}_{t_{n}-jh}\frac{\frac{t-t_n}{h}+l}{-i+l}dt =h\prod_{l=0,l\neq i}^{q}\int^{k}_{-j}\frac{s+l}{-i+l}ds,\$ $i=0,1,...,q$ ，插值节点是 $t_{n-q},...,t_n$ 。

于是得到近似公式

$u(t_{n+k})-u(t_{n-j})\approx\sum^{q}_{i=0}f\big(t_{n-i},u(t_{n-i})\big)\int^{t_{n+k}}_{t_{n-j}}L_i(t)dt$ $=h\sum^{q}_{i=0}\Bigg(\prod_{l=0,l\neq i}^{q}\int^{k}_{-j}\frac{s+l}{(-i+l)}ds\Bigg)f\big(t_{n-i},u(t_{n-i})\big)$ $=h\sum^{q}_{i=0}\beta_{qi}f\big(t_{n-i},u(t_{n-i})\big)\ (1.3.2)$

在式 $(1.3.2)$ 中，用 $u_n$ 代替 $u(t_n)$ ，仍用 $f_n$ 表示 $f(t_n,u_n)$ ，用等号代替 $\approx$ ，则得到线性多步法公式 $u_{n+k}=u_{n-j}+h\sum^{q}_{i=0}\beta_{qi}f_{n-i}\ (1.3.4)$

对 $k,j$ 和 $q$ 的不同选择，得到不同类型的具体公式。对 $k=1,j=0$ 和 $q=0,1,...$ ，我们得到Adams显式方法

$\begin{cases} u_{n+1}=u_n+h(\beta_{q0}f_n+\beta_{q1}f_{n-1}+...+\beta_{qq}f_{n-q})\\ \beta_{qi}=\int^{1}_{0}\prod^{q}_{l=0,l\neq 1}\frac{s+l}{-i+l}ds,\ i=0,1,...,q\\ \end{cases}\ (1.3.5)$

$q+1$ 步Adams显式方法和 $q$ 步Admas隐式方法的局部截断误差都与 $h^{q+2}$ 同阶，即他们是 $q+1$ 阶方法。

1.3.2 待定系数法

线性 $k$ 步方法的一般形式为 $\sum^{k}_{j=0}\alpha_ju_{n+j}=h\sum^{k}_{j=0}\beta_jf_{n+j}$ ，这里 $\alpha_j$ 和 $\beta_j$ 为实常数（注意这个 $h$ 先被提了出来），且最大步长系数 $\alpha_k\neq0$ ， $|\alpha_0|+|\beta_0|>0$ （其他系数都有可能为0）。为讨论方便，改为以 $u_{n+k}$ 为未知量， $u_{n+k-1},...,u_n$ 为已知量。我们用待定系数法来构造这种公式。

…

1.6 数值稳定性

前面各节，对一个数值方法进行理论分析时，引进了相容性，收敛性等有关概念。此时，我们作出过共同的假设，即 $h\rightarrow0$ 。这些概念在数值分析中虽然很重要，但还没有充分考虑到实际计算过程的主要特征：步长 $h$ 由于计算工具等条件的限制，不可能任意取小，更不要说趋于零。另一方面，由于每计算一步，总会产生舍入误差，而且这些舍入误差还要步步传播下去，步数增多就有可能使误差累积到影响数值解的精度，甚至于淹没了所求的近似解。因此，从实际计算的角度来看，反而要求步长 $h$ 尽可能大些。对有限的步长 $h$ ，考察一个方法对舍入误差的敏感性，这就是数值稳定性问题。这不但与计算公式有关，而且与微分方程本身的性质有关。

1.6.1 线性多步方法的绝对稳定性

首先考察一个数值例子。

…

从这个数值表看到，当用有限步长 $h=0.1$ 进行 $40$ 余步的计算之后，所得数值解与真解之值偏离相当大，已无法作为近似解提供使用。但我们知道，Milne方法 $(1.6.1)$ 是收敛的二步四阶方法，在所有的二步方法中精度阶最高，应该给出较好的数值计算结果。为什么会出现这种反常的不稳定现象？一个方法具有什么性质才能避免上述不稳定现象？这正是本节要讨论的问题。

首先假设所考虑的线性多步方法是收敛的。对于一般的线性 $k$ 步方法 $\sum^{k}_{j=0}\alpha_ju_{n+j}=h\sum^{k}_{j=0}\beta_jf_{n+j}\ (1.6.2)$

…

为了对 $e_n$ 的性质进行分析，我们作两个简化问题的假设：…；二是 $\frac{\partial f}{\partial u}=\mu=\rm const.$ ，即此时所考虑的是模型方程 $u'=\mu u$ 的求解问题。

…

记 $\overline{h}=\mu h$ ，用 $\xi_j=\xi_j(\overline{h})$ 表示稳定多项式 $\pi(\lambda;\overline{h})\equiv\rho(\lambda)-\overline{h}\sigma(\lambda)\ (1.6.7)$ 的零点，其中 $\rho(\lambda)=\lambda^k+\sum^{k-1}_{j=0}\alpha_j\lambda^j,\ \sigma(\lambda)=\sum^{k}_{j=0}\beta_j\lambda^j$ 。

…

我们感兴趣的不是对 $e_n$ 大小的估计，而是判定当 $n$ 增大时， $e_n$ 是随着增大还是减小或是振荡的问题。若 $e_n$ 是在减小，那就是说每步计算所产生的舍入误差对以后计算结果的影响在减弱，即可受到控制。据此，我们引进绝对稳定性的概念。

定义1.6.1 若对于给定的 $\overline{h}=\mu h$ ，稳定多项式 $\pi(\lambda;\overline{h})$ 的所有根 $\xi_j$ 都满足 $|\xi_j|<1,j=1,...,k$ ，则称方法 $(1.6.2)$ 关于 $\overline{h}$ 是绝对稳定的；若对实轴上的某区间（或复平面上某区域）内的任意 $\overline{h}$ ，方法都是绝对稳定的，则称此区间（或区域）为绝对稳定区间（或区域）。

Euler法

计算公式为 $u_{n+1}=u_n+hf_n$ 。对模型方程 $\frac{du}{dt}=f=\mu u$ ，可化为 $u_{n+1}-u_n=h\mu u_n$

这是线性 $1$ 步方法， $k=1$ ， $\alpha_0=-1,\alpha_1=1,$ $\beta_0=\frac{\mu}{h},\beta_1=0$ ， $\rho(\lambda)=\lambda+\alpha_0=\lambda-1,\sigma(\lambda)=\frac{\mu}{h}$

稳定多项式为 $\pi(\lambda;\overline{h})=\rho(\lambda)-\overline{h}\sigma(\lambda)=\lambda-1-\overline{h}\frac{\mu}{h}=\lambda-(1+\overline{h})$

绝对稳定区域满足 $|1+\overline{h}|<1$ 。当 $\mu$ 为复常数时，上式表示复平面上以 $-1$ 为中心的单位圆域内部。当 $\mu$ 为实常数时，绝对稳定区间为实轴上的线段 $(-2,0)$ 。

向后Euler法

计算公式为 $u_{n+1}=u_n+hf_{n+1}$ 。稳定多项式 $\pi(\lambda;\overline{h})=(1-\overline{h})\lambda-1$ 。绝对稳定区间满足 $\frac{1}{|1-\overline{h}|}<1$ ， $|1-\overline{h}|>1$ ，即方法的绝对稳定区域为复平面上以 $+1$ 为中心的单位圆域外部。特别地，当 $Re\ \mu<0$ 时，方法绝对稳定。若 $\mu$ 为实常数，绝对稳定区间为 $(-\infty,0)\cup(2,\infty)$ 。

梯形法

计算公式为 $u_{n+1}=u_n+\frac{1}{2}(f_n+f_{n+1})$

稳定多项式为 $\pi(\lambda;\overline{h})=(1-\frac{1}{2}\overline{h})\lambda-(1+\frac{1}{2}\overline{h})$

绝对稳定区间满足 $|\xi_1|=\Bigg|\frac{(1+\frac{1}{2}\overline{h})}{(1-\frac{1}{2}\overline{h})}\Bigg|<1$ 。当 $Re\ \mu<0$ 时，总有 $|\xi_1|<1$ ，故方法的绝对稳定区间为左半复平面 $Re\ \mu<0$ ，而绝对稳定区间为 $(-\infty,0)$ 。

改进Euler法

计算公式为 $\begin{cases} u_{n+1}^{(0)}=u_n+hf(t_{n},u_n)\\ u_{n+1}^{(m+1)}=u_n+\frac{h}{2}[f_n+f(t_{n+1},u_{n+1}^{(m)})],m\geq0\\ \end{cases}$

稳定多项式为 $\pi(\lambda;\overline{h})=\lambda-(1+\overline{h}+\overline{h}^2)$

绝对稳定区间满足 $|\xi_1|=|1+\overline{h}+\frac{1}{2}\overline{h}^2|=|(\frac{1}{\sqrt{2}}\overline{h}+\sqrt{2})^2-1|<1$ ，即 $(\frac{1}{\sqrt{2}}\overline{h}+\sqrt{2})^2\leq2, \overline{h}\in C$ ，这是一个复椭圆。绝对稳定区间为 $(-2,0)$ 。

二级二阶RK方法

绝对稳定区间为 $(-2,0)$ 。

四级四阶RK方法

绝对稳定区间为 $(-2.78,0)$ 。

通过本节的讨论可看到，从绝对稳定性的要求来看，对于 $Re\ \mu>0$ 的方程，没有什么可靠的算法，因为误差函数随着计算步数 $n$ 在增长。若从另一个角度来讨论问题，由于方法的收敛性，此时数值解必然随着真解的增长而增长，但只要误差函数的增长速度低于数值解的增长速度，而不影响数值解的精度，所得的计算结果也具有使用价值，这就是所谓的相对稳定性问题。

第一部分 常微分方程的数值解法

第一章 常微分方程初值问题