计算方法（六）：常微分方程初值问题的数值解法

常微分方程初值问题的数值解法

本文着重讨论一阶常微分方程初值问题：
$$\{\begin{array}{l}{y}^{\prime }=f(x,y)\\ y(x_0)=y_0\end{array}$$
的数值解法。

初值问题简介：一阶常微方程初值问题：
$$\{\begin{array}{l}\frac{dy}{dx}=f(x,y),x\in [a,b]\\ y(a)=y_0\end{array}$$
解的存在唯一性问题：设 $x_0\in [a,b],f(x,y)$ 对 $x$ 连续且关于 $y$ 满足利普希茨条件：存在常数 $L$，使 $\forall x\in [a,b]$ 及任何实数 $y_1,y_2$，有
$$|f(x,y_1)-f(x,y_2)|\le L|y_1-y_2|$$
则上述初值问题在 $[a,b]$ 上存在唯一解。

解的适定性问题：

定义：称上述初值问题对初始值 $y_0$ 和函数 $f(x,y)$ 是适定的，如果存在常数 $K>0,\eta >0$，对任意 $0<\epsilon <\eta$，当
$$\begin{gathered} |y_0-{\tilde{y}}_0|<\epsilon \\ |f(x,y)-\tilde{f}(x,y)|<\epsilon \\ \forall (x,y)\in [x_0,b]\times (-\infty ,\infty ) \end{gathered}$$
时，初值问题
$$\{\begin{array}{l}{z}^{\prime }=\tilde{f}(x,z),x_0\le x\le b\\ z(x_0)={\tilde{y}}_0\end{array}$$
的解存在，且
$$|y(x)-z(x)|\le K\epsilon$$

定理：假设 $f(x,y)$ 在 $[x_0,b]\times (-\infty ,\infty )$ 上对 $y$ 满足利普希茨条件，则初值问题是适定的。

欧拉（Euler）方法与改进欧拉方法

欧拉方法

初值问题：
$$\{\begin{array}{l}{y}^{\prime }=f(x,y)\\ y(x_0)=y_0\end{array}$$
的解 $y=y(x)$ 代表通过点 $(x_0,y_0)$ 的一条曲线，称为微分方程的积分曲线。积分曲线上每一点 $(x,y)$ 的切线的斜率 $y^{\prime }(x)$ 等于函数 $f(x,y)$ 在这点的值。

几何意义：欧拉法是过点 $(x_0,y_0)$ 作曲线的切线与 $x_1$ 交于点 $(x_1,y_1)$，用 $y_1$ 作为曲线 $y(x)$ 上的点 $(x_1,y(x_1))$ 的纵坐标 $y(x_1)$ 的近似值。

过 $(x_0,y_0)$ 以 $f(x_0,y_0)$ 为斜率的切线方程：
$$y-y_0=f(x_0,y_0)(x-x_0)$$
当 $x=x_1$ 时，得
$$y_1=y_0+f(x_0,y_0)(x_1-x_0)$$
取 $y_1$ 作为解 $y(x_1)$ 的近似值，之后再过点 $(x_1,y_1)$ 以 $f(x_1,y_1)$ 为斜率作直线
$$y-y_1=f(x_1,y_1)(x-x_1)$$
以此类推，一般地，已求得点 $(x_i,y_i)$，以 $f(x_i,y_i)$ 为斜率作直线
$$y-y_i=f(x_i,y_i)(x-x_i)$$
当 $x=x_{i+1}$ 时，得
$$y_{i+1}=y_i+f(x_i,y_i)(x_{i+1}-x_i)$$
取 $y(x_{i+1})\approx y_{i+1}$，这样就可算出一系列数值解 $y_1,y_2,\cdots ,y_n$。

通常取 $x_{i+1}-x_i=h_i=h$，则计算格式为：
$$\{\begin{array}{l}{y}_{i+1}=y_i+hf(x_i,y_i)\\ x_i=x_0+ih,i=0,1,2,\cdots \end{array}$$
数值微分和数值积分等的角度：在 $x_i$ 点微分方程：
$$y^{\prime }(x_i)=f(x_i,y(x_i))$$
用差商 $\frac{y(x_{i+1})-y(x_i)}{h}$ 代替其中的导数项 $y^{\prime }(x_i)$，即
$$\begin{gathered} y^{\prime }(x_i)\approx \frac{y(x_{i+1})-y(x_i)}{h} \\ y(x_{i+1})\approx y(x_i)+h{y}^{\prime }(x_i)=y(x_i)+hf(x_i,y(x_i)) \end{gathered}$$
用 $y=y(x_i)$ 的近似值 $y_i$ 代入上式右端，记所得结果为 $y_{i+1}$，就有：
$$y_{i+1}=y_i+hf(x_i,y_i),i=0,1,\cdots ,n-1$$
上述也称显式欧拉公式或向前欧拉公式。

若将方程 $y^{\prime }=f(x,y)$ 的两端从 $x_n$ 到 $x_{n+1}$ 求积分
$$\begin{gathered} {\int }_{x_n}^{x_{n+1}}{y}^{\prime }dx={\int }_{x_n}^{x_{n+1}}f(x,y(x))dx \\ y(x_{n+1})=y(x_n)+{\int }_{x_n}^{x_{n+1}}f(x,y(x))dx \end{gathered}$$
用左矩形计算积分项：${\int }_{x_n}^{x_{n+1}}f(x,y(x))dx\approx hf(x_n,y(x_n))$，代入上式得
$$y(x_{n+1})\approx y(x_n)+hf(x_n,y(x_n))$$
泰勒展开：对 $y(x_{n+1})$ 在 $x_n$ 处按二阶泰勒展开有
$$\begin{gathered} y(x_{n+1})=y(x_n+h)=y(x_n)+h{y}^{\prime }(x_n)+\frac{1}{2!}{h}^2{y}^{\prime \prime }({\epsilon }_n),x_n\le {\epsilon }_n\le x_{n+1} \\ y(x_{n+1})\approx y(x_n)+h{y}^{\prime }(x_n)=y(x_n)+hf(x_n,y(x_n)) \end{gathered}$$
用近似值 $y_n$ 代替 $y(x_n)$，有
$$y_{n+1}=y_n+hf(x_n,y_n)$$

欧拉公式的局部截断误差与精度分析

定义1：在 $y_i$ 准确的前提下，即 $y_j=y(x_j)(j\le i)$ 时，用数值方法计算 $y_{i+1}$ 误差
$$R_{i+1}=y(x_{i+1})-y_{i+1}$$
称为该数值方法计算 $y_{i+1}$ 时的局部截断误差。

定义2：若一个数值方法的局部截断误差为 $O(h^{p+1})$，则称这种数值方法的阶数是 $p$。显然，步长 $h(<1)$ 越小，$p$ 值越高，则局部截断误差越小，计算精度越高。

欧拉法的局部截断误差：假定 $y_i=y(x_i)$，则有
$$y_{i+1}=y(x_i)+hf(x_i,y(x_i))=y(x_i)+h{y}^{\prime }(x_i)$$
二阶泰勒公式：
$$\begin{gathered} y(x_{i+1})=y(x_i+h)=y(x_i)+h{y}^{\prime }(x_i)+\frac{h^2}{2}{y}^{\prime \prime }({\zeta }_i),x_i<{\zeta }_i<x_{i+1} \\ {R}_{i+1}=y(x_{i+1})-y_{i+1}=\frac{h^2}{2}{y}^{\prime \prime }({\zeta }_i) \end{gathered}$$
如果 $y(x)$ 有三阶导数，则：
$$R_{i+1}=y(x_{i+1})-y_{i+1}=\frac{h^2}{2}{y}^{\prime \prime }(x_i)+\frac{h^3}{6}{y}^{\prime \prime \prime }({\zeta }_i)$$
欧拉法的局部截断误差为 $O(h^2)$，欧拉方法为一阶方法。

改进欧拉方法

设改用向后差商 $\frac{1}{h}(y(x_{i+1})-y(x_i))$ 替代方程 $y^{\prime }(x_{i+1})=f(x_{i+1},y(x_{i+1}))$ 中的导数项 $y^{\prime }(x_{i+1})$，再离散化，可得：
$$y_{i+1}=y_i+hf(x_{i+1},y_{i+1})$$
称为向后欧拉公式（又称隐式欧拉公式）。

向后欧拉公式的显式化（迭代）：
$$y_{i+1}^{(k+1)}=y_i+hf(x_{i+1},y_{i+1}^{(k)}),k=0,1,2,\cdots$$
隐式欧拉公式的局部截断误差：
$$R_{i+1}=-\frac{h^2}{2}{y}^{\prime \prime }(x_i)+O(h^3)$$
两步欧拉格式：为了改善精度，改用中心差商 $\frac{1}{2h}(y(x_{i+1})-y(x_{i-1}))$ 替代方程 $y^{\prime }(x_i)=f(x_i,y(x_i))$ 中的导数项，并取离散化得：
$$y_{i+1}=y_{i-1}+2hf(x_i,y_i)$$
设 $y_i=y(x_i),y_{i-1}=y(x_{i-1})$，前两步准确，则对两步欧拉格式，有
$$y_{i+1}=y(x_{i-1})+2hf(x_i,y(x_i))$$
将 $y(x_{i+1})$ 泰勒展开：
$$\begin{gathered} y(x_{i+1})=y(x_i)+h{y}^{\prime }(x_i)+\frac{h^2}{2}{y}^{\prime \prime }(x_i)+\frac{h^3}{6}{y}^{\prime \prime \prime }({\zeta }_i) \\ y(x_{i-1})=y(x_i)-h{y}^{\prime }(x_i)+\frac{h^2}{2}{y}^{\prime \prime }(x_i)-\frac{h^3}{6}{y}^{\prime \prime \prime }({\zeta }_i) \\ y(x_{i+1})-y(x_{i-1})=2h{y}^{\prime }(x_i)+\frac{h^3}{3}{y}^{\prime \prime \prime }({\zeta }_i),x_{i-1}<{\zeta }_i<x_{i+1} \end{gathered}$$
所以有
$$R_{i+1}=y(x_{i+1})-y_{i+1}=\frac{h^3}{3}{y}^{\prime \prime \prime }({\zeta }_i)=O(h^3)$$
所以这是一种二阶方法。

梯形公式和改进的欧拉方法：将方程 $y^{\prime }=f(x,y)$ 的两端从 $x_i,x_{i+1}$，求积分得
$$y(x_{i+1})=y(x_i)+{\int }_{x_i}^{x_{i+1}}f(x,y(x))dx$$
用梯形方法代替矩形方法：
$${\int }_{x_i}^{x_{i+1}}f(x,y(x))dx\approx \frac{h}{2}[f(x_i,y(x_i))+f(x_{i+1},y(x_{i+1}))]$$
离散化之后的结果为：
$$y_{i+1}=y_i+\frac{h}{2}[f(x_i,y_i)+f(x_{i+1},y_{i+1})]$$
梯形公式显式化：
$$\{\begin{array}{l}{y}_{i+1}^{(0)}=y_i+hf(x_i,y_i)\\ y_{i+1}^{(k+1)}=y_i+\frac{h}{2}[f(x_i,y_i)+f(x_{i+1},y_{i+1}^{(k)})],k=0,1,2,\cdots \end{array}$$
用 $|y_{i+1}^{(k+1)}-y_{i+1}^{(k)}|\le \epsilon$ 控制迭代次数，其中 $\epsilon$ 为允许误差，把满足误差要求的 $y_{i+1}^{(k+1)}$ 作为 $y(x_{i+1})$ 的近似值 $y_{i+1}$。

在实用上，当 $h$ 取值较小，先用欧拉格式得一个初步近似值 $y_{i+1}^{(0)}$，称之为预估值，预估值的精度不高，用它替代梯形法右端的 $y_{i+1}$，再直接计算出 $y_{i+1}$，并称之为校正值，得到改进的欧拉方法（预估—校正公式）：
$$\{\begin{array}{l}预估：y_{i+1}^{(0)}=y_i+hf(x_i,y_i)\\ 校正：y_{i+1}=y_i+\frac{h}{2}[f(x_i,y_i)+f(x_{i+1},y_{i+1}^{(0)})]\end{array}$$
上式可以表示为嵌套形式：
$$y_{i+1}=y_i+\frac{h}{2}[f(x_i,y_i)+f(x_{i+1},y_i+hf(x_i,y_i))]$$
或者表示成下列平均化形式：
$$\{\begin{array}{l}{y}_p=y_i+hf(x_i,y_i)\\ y_c=y_i+hf(x_{i+1},y_p)\\ y_{i+1}=\frac{1}{2}(y_p+y_c)\end{array}$$
改进的欧拉方法是二阶方法。

龙格-库塔(Runge-Kutta)法

构造原理

设 $y_i=y(x_i)$，将 $y(x_{i+1})$ 在 $x_i$ 处展开：
$$y(x_{i+1})=y(x_i+h)=y(x_i)+h{y}^{\prime }(x_i)+\frac{h^2}{2}{y}^{\prime \prime }(x_i)+\frac{h^3}{3!}{y}^{\prime \prime \prime }(x_i)+...$$
取上式前两项就是局部截断误差为 $O(h^2)$ 的欧拉格式。

取前三项，可得局部截断误差为 $O(h^3)$ 的公式：
$$\begin{gathered} y(x_{i+1})\approx y(x_i)+h{y}^{\prime }(x_i)+\frac{h^2}{2}{y}^{\prime \prime }(x_i) \\ =y(x_i)+hf(x_i,y(x_i))+\frac{h^2}{2}[f_x(x_i,y(x_i))+f(x_i,y(x_i))\cdot f_y(x_i,y(x_i))] \end{gathered}$$
可以构造：
$$y_{i+1}=y_i+hf(x_i,y_i)+\frac{h^2}{2}[f_x(x_i,y_i)+f(x_i,y_i)f_y(x_i,y_i)]$$
类似的，若取前 $p+1$ 项，可得到局部截断误差为 $O(h^{p+1})$ 的数值计算公式。

基本思路：用复合函数的计算来代替各阶偏导数，通过不同点的函数值组合间接使用泰勒展开来达到高阶局部截断误差的目的。

设 $m$ 为一个正整数，称方程：
$$\{\begin{array}{l}{y}_{i+1}=y_i+h(a_1{K}_1+a_2{K}_2+\cdots +a_m{K}_m)\\ K_1=f(x_i,y_i)\\ K_2=f(x_i+{\lambda }_{2}h,y_i+{\mu }_{2}h{K}_1)\\ \cdots \cdots \\ K_m=f(x_i+{\lambda }_{m}h,y_i+{\mu }_{m}h{K}_{m-1})\end{array}$$
为 $m$ 级龙格-库塔法，其中 $0\le {\lambda }_j\le 1,a_j,{\mu }_j$ 是待定常数，由局部截断误差阶等来确定。

基本选取原则：$y_{n+1}$ 的展开表达式：
$$y_{i+1}=y_i+a_{1}h+\frac{1}{2!}{a}_2{h}^2+\frac{1}{3!}{a}_3{h}^3+\cdots$$
$y(x_{i+1})$ 在 $(x_i,y_i)$ 处的泰勒展开式：
$$y(x_{i+1})=y(x_i)+h{y}^{\prime }(x_i)+\frac{h^2}{2!}{y}^{\prime \prime }(x_i)+\frac{h^3}{3!}{y}^{\prime \prime \prime }(x_i)+\cdots$$
上述两式要有尽可能多的项重合，以减小局部截断误差。

以 $m=2$ 为例：
$$\{\begin{array}{l}{y}_{i+1}=y_i+h(a_1{K}_1+a_2{K}_2)\\ K_1=f(x_i,y_i)\\ K_2=f(x_i+{\lambda }_{2}h,y_i+{\mu }_{2}h{K}_1)\end{array}$$
将 $K_2$ 在 $(x_i,y_i)$ 处泰勒展开，且 $y_i=y(x_i)$，有：
$$\begin{gathered} y_{i+1}=y_i+h(a_1{K}_1+a_2{K}_2) \\ =y_i+a_{1}hf(x_i,y_i)+a_{2}hf(x_i+{\lambda }_{2}h,y_i+{\mu }_{2}h{K}_1) \\ =y_i+a_{1}hf(x_i,y_i)+a_{2}h[f(x_i,y_i)+{\lambda }_{2}h{f}_x(x_i,y_i)+{\mu }_{2}h{K}_1{f}_y(x_i,y_i)+O(h^2)] \\ =y_i+(a_1+a_2)f(x_i,y_i)h+a_2[{\lambda }_2{f}_x(x_i,y_i)+{\mu }_{2}f(x_i,y_i)f_y(x_i,y_i)]h^2+O(h^3) \end{gathered}$$
又：
$$\begin{gathered} y(x_{i+1})=y(x_i+h)=y(x_i)+h{y}^{\prime }(x_i)+\frac{h^2}{2}{y}^{\prime \prime }(x_i)+O(h^3) \\ =y_i+f(x_i,y_i)h+\frac{1}{2}[f_x(x_i,y_i)+f(x_i,y_i)f_y(x_i,y_i)]h^2+O(h^3) \end{gathered}$$
比较 $y_{i+1}、y(x_{i+1})$ 的表达式，选取：
$$\{\begin{array}{l}{a}_1+a_2=1\\ a_2{\lambda }_2=\frac{1}{2}\\ a_2{\mu }_2=\frac{1}{2}\end{array}$$
得到一族二阶龙格-库塔法。

取 $a_1=0,a_2=1,{\lambda }_2={\mu }_2=\frac{1}{2}$，得到中点公式：
$$y_{i+1}=y_i+hf(x_i+\frac{1}{2}h,y_i+\frac{h}{2}f(x_i,y_i))$$
取 $a_1=\frac{1}{4},a_2=\frac{3}{4},{\lambda }_2={\mu }_2=\frac{2}{3}$，得到休恩公式：
$$y_{i+1}=y_i+\frac{1}{4}h[f(x_i,y_i)+3f(x_i+\frac{2}{3}h,y_i+\frac{2}{3}hf(x_i,y_i))]$$
取 $a_1=\frac{1}{2},a_2=\frac{1}{2},{\lambda }_2={\mu }_2=1$，得到改进欧拉公式：
$$y_{i+1}=y_i+\frac{h}{2}[f(x_i,y_i)+f(x_{i+1},y+hf(x_i,y_i))]$$

经典龙格-库塔法

当 $m=4$ 时，可得到四阶经典龙格-库塔法：
$$\{\begin{array}{l}{y}_{i+1}=y_i+\frac{h}{6}(K_1+2K_2+2K_3+K_4)\\ K_1=f(x_i,y_i)\\ K_2=f(x_i+\frac{h}{2},y_i+\frac{h}{2}{K}_1)\\ K_3=f(x_i+\frac{h}{2},y_i+\frac{h}{2}{K}_2)\\ K_4=f(x_i+h,y_i+h{K}_3)\end{array}$$
具有四阶精度，即局部截断误差为 $O(h^5)$。

一般地，$m$ 级龙格-库塔法所能达到的最大阶 $p$ 的关系如下：

步长的自动选择

以四阶经典龙格-库塔法为例。从结点 $x_i$ 出发，以 $h$ 为步长求一个近似值记为 $y_{i+1}^{(h)}$，由于局部截断误差为 $O(h^5)$ ，所以有
$$y(x_{i+1})-y_{i+1}^{(h)}\approx c{h}^5$$
假定系数 $c$ 变化很慢，近似常数，并且在 $h$ 很小时，$c$ 与 $h$ 无关。然后将步长折半，即取 $\frac{h}{2}$ 为步长，从 $x_i$ 跨两步到 $x_{i+1}$，求得一个近似值 $y_{i+1}^{(\frac{h}{2})}$，每跨一步的截断误差是 $c(\frac{h}{2}{)}^5$，因此有
$$y(x_{i+1})-y_{i+1}^{(\frac{h}{2})}\approx 2c(\frac{h}{2}{)}^5$$
步长折半后，误差大约减少了 $\frac{1}{16}$，即
$$\frac{y(x_{i+1})-y_{i+1}^{(\frac{h}{2})}}{y(x_{i+1})-y_{i+1}^{(h)}}\approx \frac{1}{16}$$
由此易得出下列误差估计式：
$$y(x_{i+1})\approx \frac{2^4{y}_{i+1}^{(\frac{h}{2})}-y_{i+1}^{(h)}}{2^4-1}$$
取：
$$y_{i+1}=\frac{2^4{y}_{i+1}^{(\frac{h}{2})}-y_{i+1}^{(h)}}{2^4-1}$$
作为 $y(x_{i+1})$ 的近似值，精度比 $y_{i+1}^{(h)}、y_{i+1}^{(\frac{h}{2})}$ 都高。（理查森外推方法）

也可以写成：
$$y(x_{i+1})-y_{i+1}^{(\frac{h}{2})}\approx \frac{y_{i+1}^{(\frac{h}{2})}-y_{i+1}^{(h)}}{2^4-1}$$
可以通过检查步长折半前后两次计算结果的偏差
$$\Delta =|y_{i+1}^{(\frac{h}{2})}-y_{i+1}^{(h)}|$$
来判断所选的步长是否合适：

1. 对于给定的精度 $\epsilon$，如果 $\Delta >\epsilon$，反复将步长折半进行计算，直至 $\Delta <\epsilon$，这时取步长折半后的“新值”作为结果；
2. 如果 $\Delta <\epsilon$，反复将步长加倍，直至 $\Delta >\epsilon$，这时取步长加倍后的“老值”作为结果。

上述为变步长方法。

收敛性与稳定性

收敛性

定义3：若一个数值方法对任意固定的点 $x_i=x_0+ih$，当 $h=\frac{x_i-x_0}{i}\to 0$（即 $i\to \infty$），都有 $y_i\to y(x_i)$ ，则该方法是收敛的。

收敛性与方法的整体截断误差有关，记整体截断误差 ${\epsilon }_i=y(x_i)-y_i$ 为在 $x_i$ 处的准确值 $y(x_i)$ 与数值方法得到的近似值 $y_i$ 之间的误差，则有：

定理1：设 $f(x,y)$ 关于 $y$ 满足利普希茨条件，即存在常数 $L$，使得：
$$|f(x,y_1)-f(x,y_2)|\le L|y_1-y_2|(\forall x\in [a,b])$$
且 $y^{\prime \prime }(x)$ 有界，记 $M={\max }_{x\in [a,b]}|y^{\prime \prime }(x)|$，则欧拉方法的整体截断误差有估计式：
$$|{\epsilon }_i|\le e^{L(b-a)}|{\epsilon }_0|+\frac{Mh}{2L}(e^{L(b-a)}-1)$$
其中 ${\epsilon }_0=y(x_0)-y_0$。

当 ${\epsilon }_0=0$ 时，有
$$|{\epsilon }_i|\le \frac{Mh}{2L}(e^{L(b-a)}-1)=O(h)$$
即具有一阶收敛速度。

整体截断误差与局部截断误差之间关系： 整体截断误差=$O(h^{-1}\ast$局部截断误差）。

稳定性

如果存在正常数 $C、h_0$，使得对任意初始值 $y_0、z_0$ 的欧拉方法的解 $y_i、z_i$ 满足估计式：
$$|y_i-z_i|\le C|y_0-z_0|,x_0\le x_0+ih\le b,h\le h_0$$
则称欧拉方法是稳定的。

定理^(*)：设 $f(x,y)$ 关于 $y$ 满足利普希茨条件，则欧拉方法是稳定的。

定义4：设用某一数值方法计算 $y_i$ 时，所得到的实际计算结果为 ${\tilde{y}}_i$，且由误差 ${\delta }_i=y_i-{\tilde{y}}_i$ 引起以后各结点处 $y_j(j>i)$ 的误差为 ${\delta }_j$，如果总有 $|{\delta }_j|\le |{\delta }_i|$，则称该方法是绝对稳定的。

一个数值方法的绝对稳定性与方法本身有关，也与 $f(x,y)$ 和步长 $h$ 有关。

稳定性问题比较复杂，为简化讨论，通常考虑如下试验方程，取 $f(x,y)=\lambda y$：
$$y^{\prime }=\lambda y$$
其中 $\lambda$ 是一个复常数，记 $\tilde{h}=\lambda h$。能使某一数值方法绝对稳定的 $\tilde{h}$ 的允许取值范围称为该方法的绝对稳定域。

欧拉方法的绝对稳定性：
$$y_{i+1}=y_i+hf(x_i,y_i)=y_i+h\lambda y_i=(1+\tilde{h})y_i$$
当 $y_i$ 有误差而变为 ${\tilde{y}}_i$ 时，有
$${\tilde{y}}_{i+1}=(1+\tilde{h}){\tilde{y}}_i$$
记 ${\delta }_i=y_i-{\tilde{y}}_i$，两式相减有
$${\delta }_{i+1}=(1+\tilde{h}){\delta }_i$$
要使误差不增加，需满足：
$$|1+\tilde{h}|\le 1$$
即欧拉方法的绝对稳定区域是以(-1,0)为中心、半径为1的圆形区域。欧拉方法是条件稳定的。

向后欧拉方法绝对稳域：$|1-\tilde{h}|>1$

改进欧拉方法的绝对稳域:
$$\begin{gathered} y_{i+1}=y_i+\frac{h}{2}[\lambda y_i+\lambda (y_i+h\lambda y_i)]=(1+\tilde{h}+\frac{1}{2}{\tilde{h}}^2)y_i \\ {\delta }_{i+1}=(1+\tilde{h}+\frac{1}{2}{\tilde{h}}^2){\delta }_i \\ |1+\tilde{h}+\frac{1}{2}{\tilde{h}}^2|\le 1 \end{gathered}$$
经典龙格-库塔法的绝对稳域：
$$|1+\tilde{h}+\frac{1}{2}{\tilde{h}}^2+\frac{1}{6}{\tilde{h}}^3+\frac{1}{24}{\tilde{h}}^4|\le 1$$
对于一般方程 $\frac{dy}{dx}=f(x,y)$，可近似地取：
$$\lambda \approx -|\frac{\partial f}{\partial y}{|}_{(x_i,y_i)}$$
以便判断绝对稳定性，并用以确定求 $y_{i+1}$ 时的步长 $h_{i+1}$。

一阶方程组与高阶方程的数值解法

一阶方程组初值问题的数值解法

$$\{\begin{array}{l}\frac{d{y}_i}{dx}=f_i(x,y_1,y_2,\cdots ,y_n)\\ y_i(x_0)=y_{i0}\end{array},i=1,2,\cdots ,n$$

若把其中的未知函数、方程右端都表成向量形式：
$$Y=(y_1,y_2,\cdots ,y_n{)}^T,F=(f_1,f_2,\cdots ,f_n{)}^T$$
初值条件表示为：
$$Y(x_0)=Y_0=(y_{10},y_{20},\cdots ,y_{n0}{)}^T$$
方程可表示为：
$$\{\begin{array}{l}\frac{dY}{dx}=F(x,Y)\\ Y(x_0)=Y_0\end{array}$$

例如：
$$\{\begin{array}{l}{y}^{\prime }=f(x,y,z),y(x_0)=y_0\\ z^{\prime }=g(x,y,z),z(x_0)=z_0\end{array}$$
令 $x_i=x_0+ih$，
$$Y(x)=\left[\begin{array}{c}y(x)\\ z(x)\end{array}\right],F(x,Y)=\left[\begin{array}{c}f(x,y,z)\\ g(x,y,z)\end{array}\right],Y(x_0)=\left[\begin{array}{c}y(x_0)\\ z(x_0)\end{array}\right]=\left[\begin{array}{c}{y}_0\\ z_0\end{array}\right]$$
改进的欧拉格式的预报公式：
$${\tilde{Y}}_{i+1}=Y_i+hF(x_i,Y_i),\left[\begin{array}{c}{\tilde{y}}_{i+1}\\ {\tilde{z}}_{i+1}\end{array}\right]=\left[\begin{array}{c}{y}_i\\ z_i\end{array}\right]+h\left[\begin{array}{c}f(x_i,y_i,z_i)\\ g(x_i,y_i,z_i)\end{array}\right]$$
校正公式：
$$\{\begin{array}{l}{y}_{i+1}=y_i+\frac{h}{2}[f(x_i,y_i,z_i)+f(x_{i+1},{\tilde{y}}_{i+1},{\tilde{z}}_{i+1})]\\ z_{i+1}=z_i+\frac{h}{2}[g(x_i,y_i,z_i)+g(x_{i+1},{\tilde{y}}_{i+1},{\tilde{z}}_{i+1})]\end{array}$$
四阶经典龙格—库塔公式：
$$\begin{gathered} \{\begin{array}{l}{y}_{i+1}=y_i+\frac{h}{6}(K_1+2K_2+2K_3+K_4)\\ z_{i+1}=z_i+\frac{h}{6}(L_1+2L_2+2L_3+L_4)\end{array} \\ {K}_1=f(x_i,y_i,z_i),L_1=g(x_i,y_i,z_i) \\ {K}_2=f(x_i+\frac{h}{2},y_i+\frac{h}{2}{K}_1,z_i+\frac{h}{2}{L}_1),L_2=g(x_i+\frac{h}{2},y_i+\frac{h}{2}{K}_1,z_i+\frac{h}{2}{L}_1) \\ {K}_3=f(x_i+\frac{h}{2},y_i+\frac{h}{2}{K}_2,z_i+\frac{h}{2}{L}_2),L_3=g(x_i+\frac{h}{2},y_i+\frac{h}{2}{K}_2,z_i+\frac{h}{2}{L}_2) \\ {K}_4=f(x_i+h,y_i+h{K}_3,z_i+h{L}_3),L_4=g(x_i+h,y_i+h{K}_3,z_i+h{L}_3) \end{gathered}$$
若引入向量记号：
$$\mathbf{y}=\left[\begin{array}{c}y\\ z\end{array}\right],\mathbf{f}=\left[\begin{array}{c}f\\ g\end{array}\right],{\mathbf{y}}_i=\left[\begin{array}{c}{y}_i\\ z_i\end{array}\right],{\mathbf{K}}_i=\left[\begin{array}{c}{K}_i\\ L_i\end{array}\right](i=1,2,3,4)$$
初值问题即为：
$$\begin{gathered} \{\begin{array}{l}{\mathbf{y}}^{\prime }=\mathbf{f}(x,\mathbf{y})\\ \mathbf{y}(x_0)={\mathbf{y}}_0\end{array} \\ f(x,\mathbf{y})=f(x,y,z) \end{gathered}$$
相应的经典R-K公式为：
$$\begin{gathered} {\mathbf{y}}_{i+1}={\mathbf{y}}_i+\frac{h}{6}({\mathbf{K}}_1+2{\mathbf{K}}_2+2{\mathbf{K}}_3+{\mathbf{K}}_4) \\ {\mathbf{K}}_1=\mathbf{f}(x_i,{\mathbf{y}}_i) \\ {\mathbf{K}}_2=\mathbf{f}(x_i+\frac{h}{2},{\mathbf{y}}_i+\frac{h}{2}{\mathbf{K}}_1) \\ {\mathbf{K}}_3=\mathbf{f}(x_i+\frac{h}{2},{\mathbf{y}}_i+\frac{h}{2}{\mathbf{K}}_2) \\ {\mathbf{K}}_4=\mathbf{f}(x_i+h,{\mathbf{y}}_i+h{\mathbf{K}}_3) \end{gathered}$$

高阶方程初值问题的数值解法

对于如下形式的高阶常微分方程：
$$y^{(n)}=f(x,y,y^{\prime },\cdots ,y^{(n-1)})$$
其初值问题应在初值点 $x=x_0$ 处给出 $n$ 个条件：
$$y(x_0)=y_0,y^{\prime }(x_0)=y_0^{\prime },\cdots ,y^{(n-1)}(x_0)=y_0^{(n-1)}$$
可引进新的未知函数：$y_1=y,y_2=y^{\prime },\cdots ,y_n=y^{(n-1)}$，把上述初值问题变成一个常微分方程组：
$$\{\begin{array}{l}{y}_1^{\prime }=y_2\\ y_2^{\prime }=y_3\\ ⋮\\ y_{n-1}^{\prime }=y_n\\ y_n^{\prime }=f(x,y_1,y_2,\cdots ,y_{n-1})\end{array}$$
初始条件为：
$$y_1(x_0)=y_0,y_2(x_0)=y_0^{\prime },\cdots ,y_n(x_0)=y_0^{(n-1)}$$

例如：
$$\{\begin{array}{l}{y}^{\prime \prime }=f(x,y,y^{\prime })\\ y(x_0)=y_0\\ y^{\prime }(x_0)=y_0^{\prime }\end{array}$$
引入新的变量 $z=y^{\prime }$，可化为一解方程组的初值问题：
$$\{\begin{array}{l}{y}^{\prime }=z,y(x_0)=y_0\\ z^{\prime }=f(x,y,z),z(x_0)=y_0^{\prime }\end{array}$$
四阶经典龙格—库塔公式：
$$\begin{gathered} \{\begin{array}{l}{y}_{i+1}=y_i+\frac{h}{6}(K_1+2K_2+2K_3+K_4)\\ z_{i+1}=z_i+\frac{h}{6}(L_1+2L_2+2L_3+L_4)\end{array} \\ {K}_1=z_i,L_1=f(x_i,y_i,z_i) \\ {K}_2=z_i+\frac{h}{2}{L}_1,L_2=f(x_i+\frac{h}{2},y_i+\frac{h}{2}{K}_1,z_i+\frac{h}{2}{L}_1) \\ {K}_3=z_i+\frac{h}{2}{L}_2,L_3=f(x_i+\frac{h}{2},y_i+\frac{h}{2}{K}_2,z_i+\frac{h}{2}{L}_2) \\ {K}_4=z_i+h{L}_3,L_4=f(x_i+h,y_i+h{K}_3,z_i+h{L}_3) \end{gathered}$$

边值问题的数值解法

求解二阶常微分方程：$y^{\prime \prime }=f(x,y,y^{\prime }),x\in [a,b]$，需两个特定条件，除初值条件外。也可以通过给定边界条件来确定。

边界条件的三种给法：

1. 第一边值条件：$y(a)=\alpha ,y(b)=\beta$；
2. 第二边值条件：$y^{\prime }(a)=\alpha ,y^{\prime }(b)=\beta$；
3. 第三边值条件：$\{\begin{array}{l}{y}^{\prime }(a)-{\omega }_{1}y(a)=\alpha \\ y^{\prime }(b)+{\omega }_{2}y(b)=\beta \end{array}$，其中 ${\omega }_1\ge 0,{\omega }_2\ge 0,{\omega }_1+{\omega }_2>0$。

打靶法

考虑二阶常微分方程第一边值问题：
$$\{\begin{array}{l}{y}^{\prime \prime }=f(x,y,y^{\prime }),a<x<b\\ y(a)=\alpha ,y(b)=\beta \end{array}$$
解记为 $y(x)$。

假设存在唯一解，选取适当的 $t$，构造初值问题：
$$\{\begin{array}{l}{y}^{\prime \prime }=f(x,y,y^{\prime }),a<x<b\\ y(a)=\alpha ,y^{\prime }(a)=t\end{array}$$
解记为 $y(x,t)$，$\epsilon$ 为给定精度要求，若能使：
$$|y(b,t)-\beta |<\epsilon$$
则可把 $y(x,t)$ 作为 $y(x)$ 的近似解，即表示近似“命中”精确解。

要通过不断试算和修正的方法来获得 $t$ 的值。先凭经验确定 $y^{\prime }(a)$ 的两个预测值 $t_1,t_2$ 进行试算，分别得到两个解 $y(x,t_1),y(x,t_2)$。

记 ${\beta }_1=y(b,t_1),{\beta }_2=y(b,t_2)$，若 $|{\beta }_1-\beta |<\epsilon$，则取 $y(x,t_1)$ 作为近似解，若 $|{\beta }_2-\beta |<\epsilon$，则取 $y(x,t_2)$。否则，通过线性插值或他法确定 $t$ 的新的预测值：
$$t_3=t_1+\frac{t_2-t_1}{{\beta }_2-{\beta }_1}(\beta -{\beta }_1)$$
重新试算可得 ${\beta }_3,\cdots$，直到满足 $|{\beta }_i-\beta |<\epsilon$ 为止，其中 ${\beta }_i=y(b,t_i)$：
$$t_i=t_{i-2}+\frac{t_{i-1}-t_{i-2}}{{\beta }_{i-1}-{\beta }_{i-2}}(\beta -{\beta }_{i-2}),i=3,4,\cdots$$

1. 也可以通过牛顿迭代法等方法来得到 $t$ 的值。
2. 对于第二、第三边值问题，类似地可以得到相应的解。

有限差分法

考虑两点边值问题：
$$\{\begin{array}{l}Ly=-\frac{d}{dx}(p(x)\frac{dy}{dx})+r(x)\frac{dy}{dx}+q(x)y=f(x),a<x<b\\ y(a)=\alpha ,y(b)=\beta \end{array}$$
$p(x)\ge p_0>0,q(x)\ge 0$。在 $x=a,x=b$ 两个端点给定两个边值条件，也可给定其他带一阶导数值的边值条件。

将区间 $[a,b]$ 均匀分成 $N$ 等份网格，即：
$$x_n=a+nh,h=\frac{b-a}{N},n=0,1,\cdots ,N$$
在每个内网格结点 $x_n$，令 $Ly{|}_{x_n}=f{|}_{x_n},n=1,2,\cdots ,N-1$

假定 $y(x)$ 充分光滑，由泰勒展开式，有：
$$\{\begin{array}{l}\frac{dy}{dx}{|}_{x_n}=\frac{y(x_{n+1})-y(x_{n-1})}{2h}+O(h^2)\\ \frac{d}{dx}(p(x)\frac{dy}{dx}){|}_{x_n}=\frac{1}{h}((p\frac{dy}{dx}{)}_{x_{n+\frac{1}{2}}}-(p\frac{dy}{dx}{)}_{x_{n-\frac{1}{2}}})+O(h^2)\\ =\frac{1}{h}[p_{n+\frac{1}{2}}\frac{y(x_{n+1})-y(x_n)}{h}-p_{n-\frac{1}{2}}\frac{y(x_n)-y(x_{n-1})}{h}]+O(h^2)\end{array}$$
这里 $x_{n+\frac{1}{2}}=\frac{1}{2}(x_{n+1}+x_n),p_{n\pm \frac{1}{2}}=p(x_{n\pm \frac{1}{2}})$。

于是有：
$$\begin{gathered} -\frac{1}{h}[p_{n+\frac{1}{2}}\frac{y(x_{n+1})-y(x_n)}{h}-p_{n-\frac{1}{2}}\frac{y(x_n)-y(x_{n-1})}{h}]+r_n\frac{y(x_{n+1})-y(x_{n-1})}{2h}+q_{n}y(x_n)+R_n=f_n \\ {R}_n=O(h^2) \end{gathered}$$
忽略高阶项 $R_n$，并用 $y_n$ 表示 $y(x_n)$ 的近似值，得：
$$\begin{gathered} -\frac{1}{h}[p_{n+\frac{1}{2}}\frac{y_{n+1}-y_n}{h}-p_{n-\frac{1}{2}}\frac{y_n-y_{n-1}}{h}]+r_n\frac{y_{n+1}-y_{n-1}}{2h}+q_n{y}_n=f_n,n=1,2,\cdots ,N-1 \\ 边值条件：y_0=\alpha ,y_N=\beta \end{gathered}$$
上述方程是关于未知量 $y_n(n=1,2,\cdots ,N-1)$ 的一个线性方程组，其系数矩阵为：
$$\begin{gathered} A={\left[\begin{array}{ccccc}{b}_1& -c_1& & & \\ -a_2& b_2& -c_2& & \\ & \ddots & \ddots & \ddots & \\ & & -a_{N-2}& b_{N-2}& -c_{N-2}\\ & & & -a_{N-1}& b_{N-1}\end{array}\right]}_{(N-1)\times (N-1)} \\ {b}_n=\frac{1}{h^2}(p_{n+\frac{1}{2}}+p_{n-\frac{1}{2}})+q_n \\ {c}_n=\frac{1}{h^2}{p}_{n+\frac{1}{2}}+\frac{1}{2h}{r}_n \\ {a}_n=\frac{1}{h^2}{p}_{n-\frac{1}{2}}-\frac{1}{2h}{r}_n \end{gathered}$$
当网格步长 $h$ 充分小时，$A$ 是一个对角占优的三对角矩阵，因此 $A$ 是非奇异矩阵，可用追赶法直接求解。

若 $p(x)\equiv 1,r(x)\equiv 0,f(x)=-\gamma (x)$，则两点边值问题为：
$$\begin{gathered} \{\begin{array}{l}Ly=-y^{\prime \prime }+q(x)y=-\gamma (x),q(x)\ge 0,a<x<b\\ y(a)=\alpha ,y(b)=\beta \end{array} \\ {y}^{\prime \prime }(x_i)-q(x_i)y(x_i)=\gamma (x_i) \\ \{\begin{array}{l}\frac{y_{i+1}-2y_i+y_{i-1}}{h^2}-q_i{y}_i={\gamma }_i\\ y_0=\alpha ,y_N=\beta \end{array},i=1,2,\cdots ,N-1 \end{gathered}$$
矩阵形式：
$$\begin{gathered} Ay=b \\ A={\left[\begin{array}{ccccc}-2-q_1{h}^2& 1& & & \\ 1& -2-q_2{h}^2& 1& & \\ & \ddots & \ddots & \ddots & \\ & & 1& -2-q_{N-2}{h}^2& 1\\ & & & 1& -2-q_{N-1}{h}^2\end{array}\right]}_{(N-1)\times (N-1)} \\ y=\left[\begin{array}{c}{y}_1\\ y_2\\ ⋮\\ y_{N-2}\\ y_{N-1}\end{array}\right],b=\left[\begin{array}{c}{\gamma }_1{h}^2-\alpha \\ {\gamma }_2{h}^2\\ ⋮\\ {\gamma }_{N-2}{h}^2\\ {\gamma }_{N-1}{h}^2-\beta \end{array}\right] \end{gathered}$$