文章目录

常微分方程初值问题的数值解法
- 欧拉（Euler）方法与改进欧拉方法
- - 欧拉方法
  - 欧拉公式的局部截断误差与精度分析
  - 改进欧拉方法
- 龙格-库塔(Runge-Kutta)法
- - 构造原理
  - 经典龙格-库塔法
  - 步长的自动选择
- 收敛性与稳定性
- - 收敛性
  - 稳定性
- 一阶方程组与高阶方程的数值解法
- - 一阶方程组初值问题的数值解法
  - 高阶方程初值问题的数值解法
- 边值问题的数值解法
- - 打靶法
  - 有限差分法

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

本文着重讨论一阶常微分方程初值问题：
{y′=f(x,y)y(x0)=y0\begin{cases}y'=f(x,y)\\y(x_0)=y_0\end{cases} {y′=f(x,y)y(x0)=y0
的数值解法。

初值问题简介：一阶常微方程初值问题：
{dydx=f(x,y),x∈[a,b]y(a)=y0\begin{cases}\frac{dy}{dx}=f(x,y),x\in[a,b]\\y(a)=y_0\end{cases} {dxdy=f(x,y),x∈[a,b]y(a)=y0
解的存在唯一性问题：设 x0∈[a,b],f(x,y)x_0\in[a,b],f(x,y)x0∈[a,b],f(x,y) 对 xxx 连续且关于 yyy 满足利普希茨条件：存在常数 LLL，使 ∀x∈[a,b]\forall x\in[a,b]∀x∈[a,b] 及任何实数 y1,y2y_1,y_2y1,y2，有
∣f(x,y1)−f(x,y2)∣≤L∣y1−y2∣|f(x,y_1)-f(x,y_2)|\leq L|y_1-y_2| ∣f(x,y1)−f(x,y2)∣≤L∣y1−y2∣
则上述初值问题在 [a,b][a,b][a,b] 上存在唯一解。

解的适定性问题：

定义：称上述初值问题对初始值 y0y_0y0 和函数 f(x,y)f(x,y)f(x,y) 是适定的，如果存在常数 K>0,η>0K>0,\eta>0K>0,η>0，对任意 0<ε<η0<\varepsilon<\eta0<ε<η，当
∣y0−y~0∣<ε∣f(x,y)−f~(x,y)∣<ε∀(x,y)∈[x0,b]×(−∞,∞)|y_0-\tilde y_0|<\varepsilon\\|f(x,y)-\tilde f(x,y)|<\varepsilon\\\forall(x,y)\in[x_0,b]\times(-\infin,\infin) ∣y0−y~0∣<ε∣f(x,y)−f~(x,y)∣<ε∀(x,y)∈[x0,b]×(−∞,∞)
时，初值问题
{z′=f~(x,z),x0≤x≤bz(x0)=y~0\begin{cases}z'=\tilde f(x,z),x_0\leq x\leq b\\z(x_0)=\tilde y_0\end{cases} {z′=f~(x,z),x0≤x≤bz(x0)=y~0
的解存在，且
∣y(x)−z(x)∣≤Kε|y(x)-z(x)|\leq K\varepsilon ∣y(x)−z(x)∣≤Kε

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

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

欧拉方法

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

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

过 (x0,y0)(x_0,y_0)(x0,y0) 以 f(x0,y0)f(x_0,y_0)f(x0,y0) 为斜率的切线方程：
y−y0=f(x0,y0)(x−x0)y-y_0=f(x_0,y_0)(x-x_0) y−y0=f(x0,y0)(x−x0)
当 x=x1x=x_1x=x1 时，得
y1=y0+f(x0,y0)(x1−x0)y_1=y_0+f(x_0,y_0)(x_1-x_0) y1=y0+f(x0,y0)(x1−x0)
取 y1y_1y1 作为解 y(x1)y(x_1)y(x1) 的近似值，之后再过点 (x1,y1)(x_1,y_1)(x1,y1) 以 f(x1,y1)f(x_1,y_1)f(x1,y1) 为斜率作直线
y−y1=f(x1,y1)(x−x1)y-y_1=f(x_1,y_1)(x-x_1) y−y1=f(x1,y1)(x−x1)
以此类推，一般地，已求得点 (xi,yi)(x_i,y_i)(xi,yi)，以 f(xi,yi)f(x_i,y_i)f(xi,yi) 为斜率作直线
y−yi=f(xi,yi)(x−xi)y-y_i=f(x_i,y_i)(x-x_i) y−yi=f(xi,yi)(x−xi)
当 x=xi+1x=x_{i+1}x=xi+1 时，得
yi+1=yi+f(xi,yi)(xi+1−xi)y_{i+1}=y_i+f(x_i,y_i)(x_{i+1}-x_i) yi+1=yi+f(xi,yi)(xi+1−xi)
取 y(xi+1)≈yi+1y(x_{i+1})\approx y_{i+1}y(xi+1)≈yi+1，这样就可算出一系列数值解 y1,y2,⋯,yny_1,y_2,\cdots,y_ny1,y2,⋯,yn。

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

若将方程 y′=f(x,y)y'=f(x,y)y′=f(x,y) 的两端从 xnx_nxn 到 xn+1x_{n+1}xn+1 求积分
∫xnxn+1y′dx=∫xnxn+1f(x,y(x))dxy(xn+1)=y(xn)+∫xnxn+1f(x,y(x))dx\int_{x_n}^{x_{n+1}}y'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 ∫xnxn+1y′dx=∫xnxn+1f(x,y(x))dxy(xn+1)=y(xn)+∫xnxn+1f(x,y(x))dx
用左矩形计算积分项：∫xnxn+1f(x,y(x))dx≈hf(xn,y(xn))\int_{x_n}^{x_{n+1}}f(x,y(x))dx\approx hf(x_n,y(x_n))∫xnxn+1f(x,y(x))dx≈hf(xn,y(xn))，代入上式得
y(xn+1)≈y(xn)+hf(xn,y(xn))y(x_{n+1})\approx y(x_n)+hf(x_n,y(x_n)) y(xn+1)≈y(xn)+hf(xn,y(xn))
泰勒展开：对 y(xn+1)y(x_{n+1})y(xn+1) 在 xnx_nxn 处按二阶泰勒展开有
y(xn+1)=y(xn+h)=y(xn)+hy′(xn)+12!h2y′′(εn),xn≤εn≤xn+1y(xn+1)≈y(xn)+hy′(xn)=y(xn)+hf(xn,y(xn))y(x_{n+1})=y(x_n+h)=y(x_n)+hy'(x_n)+\frac{1}{2!}h^2y''(\varepsilon_n),x_n\leq\varepsilon_n\leq x_{n+1}\\y(x_{n+1})\approx y(x_n)+hy'(x_n)=y(x_n)+hf(x_n,y(x_n)) y(xn+1)=y(xn+h)=y(xn)+hy′(xn)+2!1h2y′′(εn),xn≤εn≤xn+1y(xn+1)≈y(xn)+hy′(xn)=y(xn)+hf(xn,y(xn))
用近似值 yny_nyn 代替 y(xn)y(x_n)y(xn)，有
yn+1=yn+hf(xn,yn)y_{n+1}=y_n+hf(x_n,y_n) yn+1=yn+hf(xn,yn)

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

定义1：在 yiy_iyi 准确的前提下，即 yj=y(xj)(j≤i)y_j=y(x_j)(j\leq i)yj=y(xj)(j≤i) 时，用数值方法计算 yi+1y_{i+1}yi+1 误差
Ri+1=y(xi+1)−yi+1R_{i+1}=y(x_{i+1})-y_{i+1} Ri+1=y(xi+1)−yi+1
称为该数值方法计算 yi+1y_{i+1}yi+1 时的局部截断误差。

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

欧拉法的局部截断误差：假定 yi=y(xi)y_i=y(x_i)yi=y(xi)，则有
yi+1=y(xi)+hf(xi,y(xi))=y(xi)+hy′(xi)y_{i+1}=y(x_i)+hf(x_i,y(x_i))=y(x_i)+hy'(x_i) yi+1=y(xi)+hf(xi,y(xi))=y(xi)+hy′(xi)
二阶泰勒公式：
y(xi+1)=y(xi+h)=y(xi)+hy′(xi)+h22y′′(ζi),xi<ζi<xi+1Ri+1=y(xi+1)−yi+1=h22y′′(ζi)y(x_{i+1})=y(x_i+h)=y(x_i)+hy'(x_i)+\frac{h^2}{2}y''(\zeta_i),x_i<\zeta_i<x_{i+1}\\R_{i+1}=y(x_{i+1})-y_{i+1}=\frac{h^2}{2}y''(\zeta_i) y(xi+1)=y(xi+h)=y(xi)+hy′(xi)+2h2y′′(ζi),xi<ζi<xi+1Ri+1=y(xi+1)−yi+1=2h2y′′(ζi)
如果 y(x)y(x)y(x) 有三阶导数，则：
Ri+1=y(xi+1)−yi+1=h22y′′(xi)+h36y′′′(ζi)R_{i+1}=y(x_{i+1})-y_{i+1}=\frac{h^2}{2}y''(x_i)+\frac{h^3}{6}y'''(\zeta_i) Ri+1=y(xi+1)−yi+1=2h2y′′(xi)+6h3y′′′(ζi)
欧拉法的局部截断误差为 O(h2)O(h^2)O(h2)，欧拉方法为一阶方法。

改进欧拉方法

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

向后欧拉公式的显式化（迭代）：
yi+1(k+1)=yi+hf(xi+1,yi+1(k)),k=0,1,2,⋯y_{i+1}^{(k+1)}=y_i+hf(x_{i+1},y_{i+1}^{(k)}),k=0,1,2,\cdots yi+1(k+1)=yi+hf(xi+1,yi+1(k)),k=0,1,2,⋯
隐式欧拉公式的局部截断误差：
Ri+1=−h22y′′(xi)+O(h3)R_{i+1}=-\frac{h^2}{2}y''(x_i)+O(h^3) Ri+1=−2h2y′′(xi)+O(h3)
两步欧拉格式：为了改善精度，改用中心差商 12h(y(xi+1)−y(xi−1))\frac{1}{2h}(y(x_{i+1})-y(x_{i-1}))2h1(y(xi+1)−y(xi−1)) 替代方程 y′(xi)=f(xi,y(xi))y'(x_i)=f(x_i,y(x_i))y′(xi)=f(xi,y(xi)) 中的导数项，并取离散化得：
yi+1=yi−1+2hf(xi,yi)y_{i+1}=y_{i-1}+2hf(x_i,y_i) yi+1=yi−1+2hf(xi,yi)
设 yi=y(xi),yi−1=y(xi−1)y_i=y(x_i),y_{i-1}=y(x_{i-1})yi=y(xi),yi−1=y(xi−1)，前两步准确，则对两步欧拉格式，有
yi+1=y(xi−1)+2hf(xi,y(xi))y_{i+1}=y(x_{i-1})+2hf(x_i,y(x_i)) yi+1=y(xi−1)+2hf(xi,y(xi))
将 y(xi+1)y(x_{i+1})y(xi+1) 泰勒展开：
y(xi+1)=y(xi)+hy′(xi)+h22y′′(xi)+h36y′′′(ζi)y(xi−1)=y(xi)−hy′(xi)+h22y′′(xi)−h36y′′′(ζi)y(xi+1)−y(xi−1)=2hy′(xi)+h33y′′′(ζi),xi−1<ζi<xi+1y(x_{i+1})=y(x_i)+hy'(x_i)+\frac{h^2}{2}y''(x_i)+\frac{h^3}{6}y'''(\zeta_i)\\y(x_{i-1})=y(x_i)-hy'(x_i)+\frac{h^2}{2}y''(x_i)-\frac{h^3}{6}y'''(\zeta_i)\\y(x_{i+1})-y(x_{i-1})=2hy'(x_i)+\frac{h^3}{3}y'''(\zeta_i),x_{i-1}<\zeta_i<x_{i+1} y(xi+1)=y(xi)+hy′(xi)+2h2y′′(xi)+6h3y′′′(ζi)y(xi−1)=y(xi)−hy′(xi)+2h2y′′(xi)−6h3y′′′(ζi)y(xi+1)−y(xi−1)=2hy′(xi)+3h3y′′′(ζi),xi−1<ζi<xi+1
所以有
Ri+1=y(xi+1)−yi+1=h33y′′′(ζi)=O(h3)R_{i+1}=y(x_{i+1})-y_{i+1}=\frac{h^3}{3}y'''(\zeta_i)=O(h^3) Ri+1=y(xi+1)−yi+1=3h3y′′′(ζi)=O(h3)
所以这是一种二阶方法。

梯形公式和改进的欧拉方法：将方程 y′=f(x,y)y'=f(x,y)y′=f(x,y) 的两端从 xi,xi+1x_i,x_{i+1}xi,xi+1，求积分得
y(xi+1)=y(xi)+∫xixi+1f(x,y(x))dxy(x_{i+1})=y(x_i)+\int_{x_i}^{x_{i+1}}f(x,y(x))dx y(xi+1)=y(xi)+∫xixi+1f(x,y(x))dx
用梯形方法代替矩形方法：
∫xixi+1f(x,y(x))dx≈h2[f(xi,y(xi))+f(xi+1,y(xi+1))]\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}))] ∫xixi+1f(x,y(x))dx≈2h[f(xi,y(xi))+f(xi+1,y(xi+1))]
离散化之后的结果为：
yi+1=yi+h2[f(xi,yi)+f(xi+1,yi+1)]y_{i+1}=y_i+\frac{h}{2}[f(x_i,y_i)+f(x_{i+1},y_{i+1})] yi+1=yi+2h[f(xi,yi)+f(xi+1,yi+1)]
梯形公式显式化：
{yi+1(0)=yi+hf(xi,yi)yi+1(k+1)=yi+h2[f(xi,yi)+f(xi+1,yi+1(k))],k=0,1,2,⋯\begin{cases}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{cases} {yi+1(0)=yi+hf(xi,yi)yi+1(k+1)=yi+2h[f(xi,yi)+f(xi+1,yi+1(k))],k=0,1,2,⋯
用 ∣yi+1(k+1)−yi+1(k)∣≤ε|y_{i+1}^{(k+1)}-y_{i+1}^{(k)}|\leq\varepsilon∣yi+1(k+1)−yi+1(k)∣≤ε 控制迭代次数，其中 ε\varepsilonε 为允许误差，把满足误差要求的 yi+1(k+1)y_{i+1}^{(k+1)}yi+1(k+1) 作为 y(xi+1)y(x_{i+1})y(xi+1) 的近似值 yi+1y_{i+1}yi+1。

在实用上，当 hhh 取值较小，先用欧拉格式得一个初步近似值 yi+1(0)y_{i+1}^{(0)}yi+1(0)，称之为预估值，预估值的精度不高，用它替代梯形法右端的 yi+1y_{i+1}yi+1，再直接计算出 yi+1y_{i+1}yi+1，并称之为校正值，得到改进的欧拉方法（预估—校正公式）：
{预估：yi+1(0)=yi+hf(xi,yi)校正：yi+1=yi+h2[f(xi,yi)+f(xi+1,yi+1(0))]\begin{cases}预估：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{cases} {预估：yi+1(0)=yi+hf(xi,yi)校正：yi+1=yi+2h[f(xi,yi)+f(xi+1,yi+1(0))]
上式可以表示为嵌套形式：
yi+1=yi+h2[f(xi,yi)+f(xi+1,yi+hf(xi,yi))]y_{i+1}=y_i+\frac{h}{2}[f(x_i,y_i)+f(x_{i+1},y_i+hf(x_i,y_i))] yi+1=yi+2h[f(xi,yi)+f(xi+1,yi+hf(xi,yi))]
或者表示成下列平均化形式：
{yp=yi+hf(xi,yi)yc=yi+hf(xi+1,yp)yi+1=12(yp+yc)\begin{cases}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{cases} ⎩⎪⎨⎪⎧yp=yi+hf(xi,yi)yc=yi+hf(xi+1,yp)yi+1=21(yp+yc)
改进的欧拉方法是二阶方法。

龙格-库塔(Runge-Kutta)法

构造原理

设 yi=y(xi)y_i=y(x_i)yi=y(xi)，将 y(xi+1)y(x_{i+1})y(xi+1) 在 xix_ixi 处展开：
y(xi+1)=y(xi+h)=y(xi)+hy′(xi)+h22y′′(xi)+h33!y′′′(xi)+...y(x_{i+1})=y(x_i+h)=y(x_i)+hy'(x_i)+\frac{h^2}{2}y''(x_i)+\frac{h^3}{3!}y'''(x_i)+... y(xi+1)=y(xi+h)=y(xi)+hy′(xi)+2h2y′′(xi)+3!h3y′′′(xi)+...
取上式前两项就是局部截断误差为 O(h2)O(h^2)O(h2) 的欧拉格式。

取前三项，可得局部截断误差为 O(h3)O(h^3)O(h3) 的公式：
y(xi+1)≈y(xi)+hy′(xi)+h22y′′(xi)=y(xi)+hf(xi,y(xi))+h22[fx(xi,y(xi))+f(xi,y(xi))⋅fy(xi,y(xi))]y(x_{i+1})\approx y(x_i)+hy'(x_i)+\frac{h^2}{2}y''(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))·f_y(x_i,y(x_i))] y(xi+1)≈y(xi)+hy′(xi)+2h2y′′(xi)=y(xi)+hf(xi,y(xi))+2h2[fx(xi,y(xi))+f(xi,y(xi))⋅fy(xi,y(xi))]
可以构造：
yi+1=yi+hf(xi,yi)+h22[fx(xi,yi)+f(xi,yi)fy(xi,yi)]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)] yi+1=yi+hf(xi,yi)+2h2[fx(xi,yi)+f(xi,yi)fy(xi,yi)]
类似的，若取前 p+1p+1p+1 项，可得到局部截断误差为 O(hp+1)O(h^{p+1})O(hp+1) 的数值计算公式。

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

设 mmm 为一个正整数，称方程：
{yi+1=yi+h(a1K1+a2K2+⋯+amKm)K1=f(xi,yi)K2=f(xi+λ2h,yi+μ2hK1)⋯⋯Km=f(xi+λmh,yi+μmhKm−1)\begin{cases}y_{i+1}=y_i+h(a_1K_1+a_2K_2+\cdots+a_mK_m)\\K_1=f(x_i,y_i)\\K_2=f(x_i+\lambda_2h,y_i+\mu_2hK_1)\\\cdots\cdots\\K_m=f(x_i+\lambda_mh,y_i+\mu_mhK_{m-1})\end{cases} ⎩⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎧yi+1=yi+h(a1K1+a2K2+⋯+amKm)K1=f(xi,yi)K2=f(xi+λ2h,yi+μ2hK1)⋯⋯Km=f(xi+λmh,yi+μmhKm−1)
为 mmm 级龙格-库塔法，其中 0≤λj≤1,aj,μj0\leq\lambda_j\leq1,a_j,\mu_j0≤λj≤1,aj,μj 是待定常数，由局部截断误差阶等来确定。

基本选取原则：yn+1y_{n+1}yn+1 的展开表达式：
yi+1=yi+a1h+12!a2h2+13!a3h3+⋯y_{i+1}=y_i+a_1h+\frac{1}{2!}a_2h^2+\frac{1}{3!}a_3h^3+\cdots yi+1=yi+a1h+2!1a2h2+3!1a3h3+⋯
y(xi+1)y(x_{i+1})y(xi+1) 在 (xi,yi)(x_i,y_i)(xi,yi) 处的泰勒展开式：
y(xi+1)=y(xi)+hy′(xi)+h22!y′′(xi)+h33!y′′′(xi)+⋯y(x_{i+1})=y(x_i)+hy'(x_i)+\frac{h^2}{2!}y''(x_i)+\frac{h^3}{3!}y'''(x_i)+\cdots y(xi+1)=y(xi)+hy′(xi)+2!h2y′′(xi)+3!h3y′′′(xi)+⋯
上述两式要有尽可能多的项重合，以减小局部截断误差。

以 m=2m=2m=2 为例：
{yi+1=yi+h(a1K1+a2K2)K1=f(xi,yi)K2=f(xi+λ2h,yi+μ2hK1)\begin{cases}y_{i+1}=y_i+h(a_1K_1+a_2K_2)\\K_1=f(x_i,y_i)\\K_2=f(x_i+\lambda_2h,y_i+\mu_2hK_1)\end{cases} ⎩⎪⎨⎪⎧yi+1=yi+h(a1K1+a2K2)K1=f(xi,yi)K2=f(xi+λ2h,yi+μ2hK1)
将 K2K_2K2 在 (xi,yi)(x_i,y_i)(xi,yi) 处泰勒展开，且 yi=y(xi)y_i=y(x_i)yi=y(xi)，有：
yi+1=yi+h(a1K1+a2K2)=yi+a1hf(xi,yi)+a2hf(xi+λ2h,yi+μ2hK1)=yi+a1hf(xi,yi)+a2h[f(xi,yi)+λ2hfx(xi,yi)+μ2hK1fy(xi,yi)+O(h2)]=yi+(a1+a2)f(xi,yi)h+a2[λ2fx(xi,yi)+μ2f(xi,yi)fy(xi,yi)]h2+O(h3)y_{i+1}=y_i+h(a_1K_1+a_2K_2)\\=y_i+a_1hf(x_i,y_i)+a_2hf(x_i+\lambda_2h,y_i+\mu_2hK_1)\\=y_i+a_1hf(x_i,y_i)+a_2h[f(x_i,y_i)+\lambda_2hf_x(x_i,y_i)+\mu_2hK_1f_y(x_i,y_i)+O(h^2)]\\=y_i+(a_1+a_2)f(x_i,y_i)h+a_2[\lambda_2f_x(x_i,y_i)+\mu_2f(x_i,y_i)f_y(x_i,y_i)]h^2+O(h^3) yi+1=yi+h(a1K1+a2K2)=yi+a1hf(xi,yi)+a2hf(xi+λ2h,yi+μ2hK1)=yi+a1hf(xi,yi)+a2h[f(xi,yi)+λ2hfx(xi,yi)+μ2hK1fy(xi,yi)+O(h2)]=yi+(a1+a2)f(xi,yi)h+a2[λ2fx(xi,yi)+μ2f(xi,yi)fy(xi,yi)]h2+O(h3)
又：
y(xi+1)=y(xi+h)=y(xi)+hy′(xi)+h22y′′(xi)+O(h3)=yi+f(xi,yi)h+12[fx(xi,yi)+f(xi,yi)fy(xi,yi)]h2+O(h3)y(x_{i+1})=y(x_i+h)=y(x_i)+hy'(x_i)+\frac{h^2}{2}y''(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) y(xi+1)=y(xi+h)=y(xi)+hy′(xi)+2h2y′′(xi)+O(h3)=yi+f(xi,yi)h+21[fx(xi,yi)+f(xi,yi)fy(xi,yi)]h2+O(h3)
比较 yi+1、y(xi+1)y_{i+1}、y(x_{i+1})yi+1、y(xi+1) 的表达式，选取：
{a1+a2=1a2λ2=12a2μ2=12\begin{cases}a_1+a_2=1\\a_2\lambda_2=\dfrac{1}{2}\\a_2\mu_2=\dfrac{1}{2}\end{cases} ⎩⎪⎪⎪⎨⎪⎪⎪⎧a1+a2=1a2λ2=21a2μ2=21
得到一族二阶龙格-库塔法。

取 a1=0,a2=1,λ2=μ2=12a_1=0,a_2=1,\lambda_2=\mu_2=\dfrac{1}{2}a1=0,a2=1,λ2=μ2=21，得到中点公式：
yi+1=yi+hf(xi+12h,yi+h2f(xi,yi))y_{i+1}=y_i+hf(x_i+\dfrac{1}{2}h,y_i+\dfrac{h}{2}f(x_i,y_i)) yi+1=yi+hf(xi+21h,yi+2hf(xi,yi))
取 a1=14,a2=34,λ2=μ2=23a_1=\dfrac{1}{4},a_2=\dfrac{3}{4},\lambda_2=\mu_2=\dfrac{2}{3}a1=41,a2=43,λ2=μ2=32，得到休恩公式：
yi+1=yi+14h[f(xi,yi)+3f(xi+23h,yi+23hf(xi,yi))]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))] yi+1=yi+41h[f(xi,yi)+3f(xi+32h,yi+32hf(xi,yi))]
取 a1=12,a2=12,λ2=μ2=1a_1=\dfrac{1}{2},a_2=\dfrac{1}{2},\lambda_2=\mu_2=1a1=21,a2=21,λ2=μ2=1，得到改进欧拉公式：
yi+1=yi+h2[f(xi,yi)+f(xi+1,y+hf(xi,yi))]y_{i+1}=y_i+\frac{h}{2}[f(x_i,y_i)+f(x_{i+1},y+hf(x_i,y_i))] yi+1=yi+2h[f(xi,yi)+f(xi+1,y+hf(xi,yi))]

经典龙格-库塔法

当 m=4m=4m=4 时，可得到四阶经典龙格-库塔法：
{yi+1=yi+h6(K1+2K2+2K3+K4)K1=f(xi,yi)K2=f(xi+h2,yi+h2K1)K3=f(xi+h2,yi+h2K2)K4=f(xi+h,yi+hK3)\begin{cases}y_{i+1}=y_i+\dfrac{h}{6}(K_1+2K_2+2K_3+K_4)\\K_1=f(x_i,y_i)\\K_2=f(x_i+\dfrac{h}{2},y_i+\dfrac{h}{2}K_1)\\K_3=f(x_i+\dfrac{h}{2},y_i+\dfrac{h}{2}K_2)\\K_4=f(x_i+h,y_i+hK_3)\end{cases} ⎩⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎧yi+1=yi+6h(K1+2K2+2K3+K4)K1=f(xi,yi)K2=f(xi+2h,yi+2hK1)K3=f(xi+2h,yi+2hK2)K4=f(xi+h,yi+hK3)
具有四阶精度，即局部截断误差为 O(h5)O(h^5)O(h5)。

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

步长的自动选择

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

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

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

上述为变步长方法。

收敛性与稳定性

收敛性

定义3：若一个数值方法对任意固定的点 xi=x0+ihx_i=x_0+ihxi=x0+ih，当 h=xi−x0i→0h=\dfrac{x_i-x_0}{i}\to0h=ixi−x0→0（即 i→∞i\to\infini→∞），都有 yi→y(xi)y_i\to y(x_i)yi→y(xi) ，则该方法是收敛的。

收敛性与方法的整体截断误差有关，记整体截断误差 εi=y(xi)−yi\varepsilon_i=y(x_i)-y_iεi=y(xi)−yi 为在 xix_ixi 处的准确值 y(xi)y(x_i)y(xi) 与数值方法得到的近似值 yiy_iyi 之间的误差，则有：

定理1：设 f(x,y)f(x,y)f(x,y) 关于 yyy 满足利普希茨条件，即存在常数 LLL，使得：
∣f(x,y1)−f(x,y2)∣≤L∣y1−y2∣(∀x∈[a,b])|f(x,y_1)-f(x,y_2)|\leq L|y_1-y_2|\ \ \ (\forall x\in[a,b]) ∣f(x,y1)−f(x,y2)∣≤L∣y1−y2∣ (∀x∈[a,b])
且 y′′(x)y''(x)y′′(x) 有界，记 M=max⁡x∈[a,b]∣y′′(x)∣M=\max_{x\in[a,b]}|y''(x)|M=maxx∈[a,b]∣y′′(x)∣，则欧拉方法的整体截断误差有估计式：
∣εi∣≤eL(b−a)∣ε0∣+Mh2L(eL(b−a)−1)|\varepsilon_i|\leq e^{L(b-a)}|\varepsilon_0|+\frac{Mh}{2L}(e^{L(b-a)}-1) ∣εi∣≤eL(b−a)∣ε0∣+2LMh(eL(b−a)−1)
其中 ε0=y(x0)−y0\varepsilon_0=y(x_0)-y_0ε0=y(x0)−y0。

当 ε0=0\varepsilon_0=0ε0=0 时，有
∣εi∣≤Mh2L(eL(b−a)−1)=O(h)|\varepsilon_i|\leq\frac{Mh}{2L}(e^{L(b-a)}-1)=O(h) ∣εi∣≤2LMh(eL(b−a)−1)=O(h)
即具有一阶收敛速度。

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

稳定性

如果存在正常数 C、h0C、h_0C、h0，使得对任意初始值 y0、z0y_0、z_0y0、z0 的欧拉方法的解 yi、ziy_i、z_iyi、zi 满足估计式：
∣yi−zi∣≤C∣y0−z0∣,x0≤x0+ih≤b,h≤h0|y_i-z_i|\leq C|y_0-z_0|, \ x_0\leq x_0+ih\leq b,\ \ h\leq h_0 ∣yi−zi∣≤C∣y0−z0∣, x0≤x0+ih≤b, h≤h0
则称欧拉方法是稳定的。

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

定义4：设用某一数值方法计算 yiy_iyi 时，所得到的实际计算结果为 y~i\tilde y_iy~i，且由误差 δi=yi−y~i\delta_i=y_i-\tilde y_iδi=yi−y~i 引起以后各结点处 yj(j>i)y_j(j>i)yj(j>i) 的误差为 δj\delta_jδj，如果总有 ∣δj∣≤∣δi∣|\delta_j|\leq|\delta_i|∣δj∣≤∣δi∣，则称该方法是绝对稳定的。

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

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

欧拉方法的绝对稳定性：
yi+1=yi+hf(xi,yi)=yi+hλyi=(1+h~)yiy_{i+1}=y_i+hf(x_i,y_i)=y_i+h\lambda y_i=(1+\tilde h)y_i yi+1=yi+hf(xi,yi)=yi+hλyi=(1+h~)yi
当 yiy_iyi 有误差而变为 y~i\tilde y_iy~i 时，有
y~i+1=(1+h~)y~i\tilde y_{i+1}=(1+\tilde h)\tilde y_i y~i+1=(1+h~)y~i
记 δi=yi−y~i\delta_i=y_i-\tilde y_iδi=yi−y~i，两式相减有
δi+1=(1+h~)δi\delta_{i+1}=(1+\tilde h)\delta_i δi+1=(1+h~)δi
要使误差不增加，需满足：
∣1+h~∣≤1|1+\tilde h|\leq1 ∣1+h~∣≤1
即欧拉方法的绝对稳定区域是以(-1,0)为中心、半径为1的圆形区域。欧拉方法是条件稳定的。

向后欧拉方法绝对稳域：∣1−h~∣>1|1-\tilde h|>1∣1−h~∣>1

改进欧拉方法的绝对稳域:
yi+1=yi+h2[λyi+λ(yi+hλyi)]=(1+h~+12h~2)yiδi+1=(1+h~+12h~2)δi∣1+h~+12h~2∣≤1y_{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|\leq1 yi+1=yi+2h[λyi+λ(yi+hλyi)]=(1+h~+21h~2)yiδi+1=(1+h~+21h~2)δi∣1+h~+21h~2∣≤1
经典龙格-库塔法的绝对稳域：
∣1+h~+12h~2+16h~3+124h~4∣≤1|1+\tilde h+\frac{1}{2}\tilde h^2+\frac{1}{6}\tilde h^3+\frac{1}{24}\tilde h^4|\leq1 ∣1+h~+21h~2+61h~3+241h~4∣≤1
对于一般方程 dydx=f(x,y)\dfrac{dy}{dx}=f(x,y)dxdy=f(x,y)，可近似地取：
λ≈−∣∂f∂y∣(xi,yi)\lambda\approx-|\frac{\partial f}{\partial y}|_{(x_i,y_i)} λ≈−∣∂y∂f∣(xi,yi)
以便判断绝对稳定性，并用以确定求 yi+1y_{i+1}yi+1 时的步长 hi+1h_{i+1}hi+1。

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

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

{dyidx=fi(x,y1,y2,⋯,yn)yi(x0)=yi0,i=1,2,⋯,n\begin{cases}\dfrac{dy_i}{dx}=f_i(x,y_1,y_2,\cdots,y_n)\\y_i(x_0)=y_{i0}\end{cases},i=1,2,\cdots,n ⎩⎨⎧dxdyi=fi(x,y1,y2,⋯,yn)yi(x0)=yi0,i=1,2,⋯,n

若把其中的未知函数、方程右端都表成向量形式：
Y=(y1,y2,⋯,yn)T,F=(f1,f2,⋯,fn)TY=(y_1,y_2,\cdots,y_n)^T,F=(f_1,f_2,\cdots,f_n)^T Y=(y1,y2,⋯,yn)T,F=(f1,f2,⋯,fn)T
初值条件表示为：
Y(x0)=Y0=(y10,y20,⋯,yn0)TY(x_0)=Y_0=(y_{10},y_{20},\cdots,y_{n0})^T Y(x0)=Y0=(y10,y20,⋯,yn0)T
方程可表示为：
{dYdx=F(x,Y)Y(x0)=Y0\begin{cases}\dfrac{dY}{dx}=F(x,Y)\\Y(x_0)=Y_0\end{cases} ⎩⎨⎧dxdY=F(x,Y)Y(x0)=Y0

例如：
{y′=f(x,y,z),y(x0)=y0z′=g(x,y,z),z(x0)=z0\begin{cases}y'=f(x,y,z),y(x_0)=y_0\\z'=g(x,y,z),z(x_0)=z_0\end{cases}\\ {y′=f(x,y,z),y(x0)=y0z′=g(x,y,z),z(x0)=z0
令 xi=x0+ihx_i=x_0+ihxi=x0+ih，
Y(x)=[y(x)z(x)],F(x,Y)=[f(x,y,z)g(x,y,z)],Y(x0)=[y(x0)z(x0)]=[y0z0]Y(x)=\left[\begin{matrix}y(x)\\z(x)\end{matrix}\right],F(x,Y)=\left[\begin{matrix}f(x,y,z)\\g(x,y,z)\end{matrix}\right],Y(x_0)=\left[\begin{matrix}y(x_0)\\z(x_0)\end{matrix}\right]=\left[\begin{matrix}y_0\\z_0\end{matrix}\right] Y(x)=[y(x)z(x)],F(x,Y)=[f(x,y,z)g(x,y,z)],Y(x0)=[y(x0)z(x0)]=[y0z0]
改进的欧拉格式的预报公式：
Y~i+1=Yi+hF(xi,Yi),[y~i+1z~i+1]=[yizi]+h[f(xi,yi,zi)g(xi,yi,zi)]\tilde Y_{i+1}=Y_i+hF(x_i,Y_i),\left[\begin{matrix}\tilde y_{i+1}\\\tilde z_{i+1}\end{matrix}\right]=\left[\begin{matrix}y_i\\z_i\end{matrix}\right]+h\left[\begin{matrix}f(x_i,y_i,z_i)\\g(x_i,y_i,z_i)\end{matrix}\right] Y~i+1=Yi+hF(xi,Yi),[y~i+1z~i+1]=[yizi]+h[f(xi,yi,zi)g(xi,yi,zi)]
校正公式：
{yi+1=yi+h2[f(xi,yi,zi)+f(xi+1,y~i+1,z~i+1)]zi+1=zi+h2[g(xi,yi,zi)+g(xi+1,y~i+1,z~i+1)]\begin{cases}y_{i+1}=y_i+\dfrac{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+\dfrac{h}{2}[g(x_i,y_i,z_i)+g(x_{i+1},\tilde y_{i+1},\tilde z_{i+1})]\end{cases} ⎩⎪⎨⎪⎧yi+1=yi+2h[f(xi,yi,zi)+f(xi+1,y~i+1,z~i+1)]zi+1=zi+2h[g(xi,yi,zi)+g(xi+1,y~i+1,z~i+1)]
四阶经典龙格—库塔公式：
{yi+1=yi+h6(K1+2K2+2K3+K4)zi+1=zi+h6(L1+2L2+2L3+L4)K1=f(xi,yi,zi),L1=g(xi,yi,zi)K2=f(xi+h2,yi+h2K1,zi+h2L1),L2=g(xi+h2,yi+h2K1,zi+h2L1)K3=f(xi+h2,yi+h2K2,zi+h2L2),L3=g(xi+h2,yi+h2K2,zi+h2L2)K4=f(xi+h,yi+hK3,zi+hL3),L4=g(xi+h,yi+hK3,zi+hL3)\begin{cases}y_{i+1}=y_i+\dfrac{h}{6}(K_1+2K_2+2K_3+K_4)\\z_{i+1}=z_i+\dfrac{h}{6}(L_1+2L_2+2L_3+L_4)\end{cases}\\K_1=f(x_i,y_i,z_i),L_1=g(x_i,y_i,z_i)\\K_2=f(x_i+\dfrac{h}{2},y_i+\dfrac{h}{2}K_1,z_i+\dfrac{h}{2}L_1),L_2=g(x_i+\dfrac{h}{2},y_i+\dfrac{h}{2}K_1,z_i+\dfrac{h}{2}L_1)\\K_3=f(x_i+\dfrac{h}{2},y_i+\dfrac{h}{2}K_2,z_i+\dfrac{h}{2}L_2),L_3=g(x_i+\dfrac{h}{2},y_i+\dfrac{h}{2}K_2,z_i+\dfrac{h}{2}L_2)\\K_4=f(x_i+h,y_i+hK_3,z_i+hL_3),L_4=g(x_i+h,y_i+hK_3,z_i+hL_3) ⎩⎪⎨⎪⎧yi+1=yi+6h(K1+2K2+2K3+K4)zi+1=zi+6h(L1+2L2+2L3+L4)K1=f(xi,yi,zi),L1=g(xi,yi,zi)K2=f(xi+2h,yi+2hK1,zi+2hL1),L2=g(xi+2h,yi+2hK1,zi+2hL1)K3=f(xi+2h,yi+2hK2,zi+2hL2),L3=g(xi+2h,yi+2hK2,zi+2hL2)K4=f(xi+h,yi+hK3,zi+hL3),L4=g(xi+h,yi+hK3,zi+hL3)
若引入向量记号：
y=[yz],f=[fg],yi=[yizi],Ki=[KiLi](i=1,2,3,4)\boldsymbol y=\left[\begin{matrix}y\\z\end{matrix}\right],\boldsymbol f=\left[\begin{matrix}f\\g\end{matrix}\right],\boldsymbol y_i=\left[\begin{matrix}y_i\\z_i\end{matrix}\right],\boldsymbol K_i=\left[\begin{matrix}K_i\\L_i\end{matrix}\right]\ \ (i=1,2,3,4) y=[yz],f=[fg],yi=[yizi],Ki=[KiLi] (i=1,2,3,4)
初值问题即为：
{y′=f(x,y)y(x0)=y0f(x,y)=f(x,y,z)\begin{cases}\boldsymbol y'=\boldsymbol f(x,\boldsymbol y)\\\boldsymbol y(x_0)=\boldsymbol y_0\end{cases}\\f(x,\boldsymbol y)=f(x,y,z) {y′=f(x,y)y(x0)=y0f(x,y)=f(x,y,z)
相应的经典R-K公式为：
yi+1=yi+h6(K1+2K2+2K3+K4)K1=f(xi,yi)K2=f(xi+h2,yi+h2K1)K3=f(xi+h2,yi+h2K2)K4=f(xi+h,yi+hK3)\boldsymbol y_{i+1}=\boldsymbol y_i+\frac{h}{6}(\boldsymbol K_1+2\boldsymbol K_2+2\boldsymbol K_3+\boldsymbol K_4)\\\boldsymbol K_1=\boldsymbol f(x_i,\boldsymbol y_i)\\\boldsymbol K_2=\boldsymbol f(x_i+\frac{h}{2},\boldsymbol y_i+\frac{h}{2}\boldsymbol K_1)\\\boldsymbol K_3=\boldsymbol f(x_i+\frac{h}{2},\boldsymbol y_i+\frac{h}{2}\boldsymbol K_2)\\\boldsymbol K_4=\boldsymbol f(x_i+h,\boldsymbol y_i+h\boldsymbol K_3) yi+1=yi+6h(K1+2K2+2K3+K4)K1=f(xi,yi)K2=f(xi+2h,yi+2hK1)K3=f(xi+2h,yi+2hK2)K4=f(xi+h,yi+hK3)

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

对于如下形式的高阶常微分方程：
y(n)=f(x,y,y′,⋯,y(n−1))y^{(n)}=f(x,y,y',\cdots,y^{(n-1)}) y(n)=f(x,y,y′,⋯,y(n−1))
其初值问题应在初值点 x=x0x=x_0x=x0 处给出 nnn 个条件：
y(x0)=y0,y′(x0)=y0′,⋯,y(n−1)(x0)=y0(n−1)y(x_0)=y_0,y'(x_0)=y_0',\cdots,y^{(n-1)}(x_0)=y^{(n-1)}_0 y(x0)=y0,y′(x0)=y0′,⋯,y(n−1)(x0)=y0(n−1)
可引进新的未知函数：y1=y,y2=y′,⋯,yn=y(n−1)y_1=y,y_2=y',\cdots,y_n=y^{(n-1)}y1=y,y2=y′,⋯,yn=y(n−1)，把上述初值问题变成一个常微分方程组：
{y1′=y2y2′=y3⋮yn−1′=ynyn′=f(x,y1,y2,⋯,yn−1)\begin{cases}y_1'=y_2\\y_2'=y_3\\\vdots\\y'_{n-1}=y_n\\y_n'=f(x,y_1,y_2,\cdots,y_{n-1})\end{cases} ⎩⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎧y1′=y2y2′=y3⋮yn−1′=ynyn′=f(x,y1,y2,⋯,yn−1)
初始条件为：
y1(x0)=y0,y2(x0)=y0′,⋯,yn(x0)=y0(n−1)y_1(x_0)=y_0,y_2(x_0)=y_0',\cdots,y_n(x_0)=y_0^{(n-1)} y1(x0)=y0,y2(x0)=y0′,⋯,yn(x0)=y0(n−1)

例如：
{y′′=f(x,y,y′)y(x0)=y0y′(x0)=y0′\begin{cases}y''=f(x,y,y')\\y(x_0)=y_0\\y'(x_0)=y_0'\end{cases} ⎩⎪⎨⎪⎧y′′=f(x,y,y′)y(x0)=y0y′(x0)=y0′
引入新的变量 z=y′z=y'z=y′，可化为一解方程组的初值问题：
{y′=z,y(x0)=y0z′=f(x,y,z),z(x0)=y0′\begin{cases}y'=z,y(x_0)=y_0\\z'=f(x,y,z),z(x_0)=y_0'\end{cases} {y′=z,y(x0)=y0z′=f(x,y,z),z(x0)=y0′
四阶经典龙格—库塔公式：
{yi+1=yi+h6(K1+2K2+2K3+K4)zi+1=zi+h6(L1+2L2+2L3+L4)K1=zi,L1=f(xi,yi,zi)K2=zi+h2L1,L2=f(xi+h2,yi+h2K1,zi+h2L1)K3=zi+h2L2,L3=f(xi+h2,yi+h2K2,zi+h2L2)K4=zi+hL3,L4=f(xi+h,yi+hK3,zi+hL3)\begin{cases}y_{i+1}=y_i+\dfrac{h}{6}(K_1+2K_2+2K_3+K_4)\\z_{i+1}=z_i+\dfrac{h}{6}(L_1+2L_2+2L_3+L_4)\end{cases}\\K_1=z_i,L_1=f(x_i,y_i,z_i)\\K_2=z_i+\dfrac{h}{2}L_1,L_2=f(x_i+\dfrac{h}{2},y_i+\dfrac{h}{2}K_1,z_i+\dfrac{h}{2}L_1)\\K_3=z_i+\dfrac{h}{2}L_2,L_3=f(x_i+\dfrac{h}{2},y_i+\dfrac{h}{2}K_2,z_i+\dfrac{h}{2}L_2)\\K_4=z_i+hL_3,L_4=f(x_i+h,y_i+hK_3,z_i+hL_3) ⎩⎪⎨⎪⎧yi+1=yi+6h(K1+2K2+2K3+K4)zi+1=zi+6h(L1+2L2+2L3+L4)K1=zi,L1=f(xi,yi,zi)K2=zi+2hL1,L2=f(xi+2h,yi+2hK1,zi+2hL1)K3=zi+2hL2,L3=f(xi+2h,yi+2hK2,zi+2hL2)K4=zi+hL3,L4=f(xi+h,yi+hK3,zi+hL3)

边值问题的数值解法

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

边界条件的三种给法：

第一边值条件：y(a)=α,y(b)=βy(a)=\alpha,y(b)=\betay(a)=α,y(b)=β；
第二边值条件：y′(a)=α,y′(b)=βy'(a)=\alpha,y'(b)=\betay′(a)=α,y′(b)=β；
第三边值条件：{y′(a)−ω1y(a)=αy′(b)+ω2y(b)=β\begin{cases}y'(a)-\omega_1y(a)=\alpha\\y'(b)+\omega_2y(b)=\beta\end{cases}{y′(a)−ω1y(a)=αy′(b)+ω2y(b)=β，其中 ω1≥0,ω2≥0,ω1+ω2>0\omega_1\geq0,\omega_2\geq0,\omega_1+\omega_2>0ω1≥0,ω2≥0,ω1+ω2>0。

打靶法

考虑二阶常微分方程第一边值问题：
{y′′=f(x,y,y′),a<x<by(a)=α,y(b)=β\begin{cases}y''=f(x,y,y'),a<x<b\\y(a)=\alpha,y(b)=\beta\end{cases} {y′′=f(x,y,y′),a<x<by(a)=α,y(b)=β
解记为 y(x)y(x)y(x)。

假设存在唯一解，选取适当的 ttt，构造初值问题：
{y′′=f(x,y,y′),a<x<by(a)=α,y′(a)=t\begin{cases}y''=f(x,y,y'),a<x<b\\y(a)=\alpha,y'(a)=t\end{cases} {y′′=f(x,y,y′),a<x<by(a)=α,y′(a)=t
解记为 y(x,t)y(x,t)y(x,t)，ε\varepsilonε 为给定精度要求，若能使：
∣y(b,t)−β∣<ε|y(b,t)-\beta|<\varepsilon ∣y(b,t)−β∣<ε
则可把 y(x,t)y(x,t)y(x,t) 作为 y(x)y(x)y(x) 的近似解，即表示近似“命中”精确解。

要通过不断试算和修正的方法来获得 ttt 的值。先凭经验确定 y′(a)y'(a)y′(a) 的两个预测值 t1,t2t_1,t_2t1,t2 进行试算，分别得到两个解 y(x,t1),y(x,t2)y(x,t_1),y(x,t_2)y(x,t1),y(x,t2)。

记 β1=y(b,t1),β2=y(b,t2)\beta_1=y(b,t_1),\beta_2=y(b,t_2)β1=y(b,t1),β2=y(b,t2)，若 ∣β1−β∣<ε|\beta_1-\beta|<\varepsilon∣β1−β∣<ε，则取 y(x,t1)y(x,t_1)y(x,t1) 作为近似解，若 ∣β2−β∣<ε|\beta_2-\beta|<\varepsilon∣β2−β∣<ε，则取 y(x,t2)y(x,t_2)y(x,t2)。否则，通过线性插值或他法确定 ttt 的新的预测值：
t3=t1+t2−t1β2−β1(β−β1)t_3=t_1+\frac{t_2-t_1}{\beta_2-\beta_1}(\beta-\beta_1) t3=t1+β2−β1t2−t1(β−β1)
重新试算可得 β3,⋯\beta_3,\cdotsβ3,⋯，直到满足 ∣βi−β∣<ε|\beta_i-\beta|<\varepsilon∣βi−β∣<ε 为止，其中 βi=y(b,ti)\beta_i=y(b,t_i)βi=y(b,ti)：
ti=ti−2+ti−1−ti−2βi−1−βi−2(β−βi−2),i=3,4,⋯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 ti=ti−2+βi−1−βi−2ti−1−ti−2(β−βi−2),i=3,4,⋯

也可以通过牛顿迭代法等方法来得到 ttt 的值。

对于第二、第三边值问题，类似地可以得到相应的解。

有限差分法

考虑两点边值问题：
{Ly=−ddx(p(x)dydx)+r(x)dydx+q(x)y=f(x),a<x<by(a)=α,y(b)=β\begin{cases}Ly=-\dfrac{d}{dx}(p(x)\dfrac{dy}{dx})+r(x)\dfrac{dy}{dx}+q(x)y=f(x),a<x<b\\y(a)=\alpha,y(b)=\beta\end{cases} ⎩⎨⎧Ly=−dxd(p(x)dxdy)+r(x)dxdy+q(x)y=f(x),a<x<by(a)=α,y(b)=β
p(x)≥p0>0,q(x)≥0p(x)\geq p_0>0,q(x)\geq0p(x)≥p0>0,q(x)≥0。在 x=a,x=bx=a,x=bx=a,x=b 两个端点给定两个边值条件，也可给定其他带一阶导数值的边值条件。

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

假定 y(x)y(x)y(x) 充分光滑，由泰勒展开式，有：
{dydx∣xn=y(xn+1)−y(xn−1)2h+O(h2)ddx(p(x)dydx)∣xn=1h((pdydx)xn+12−(pdydx)xn−12)+O(h2)=1h[pn+12y(xn+1)−y(xn)h−pn−12y(xn)−y(xn−1)h]+O(h2)\begin{cases}\dfrac{dy}{dx}|_{x_n}=\dfrac{y(x_{n+1})-y(x_{n-1})}{2h}+O(h^2)\\\dfrac{d}{dx}(p(x)\dfrac{dy}{dx})|_{x_n}=\dfrac{1}{h}((p\dfrac{dy}{dx})_{x_{n+\frac{1}{2}}}-(p\dfrac{dy}{dx})_{x_{n-\frac{1}{2}}})+O(h^2)\\=\dfrac{1}{h}[p_{n+\frac{1}{2}}\dfrac{y(x_{n+1})-y(x_n)}{h}-p_{n-\frac{1}{2}}\dfrac{y(x_n)-y(x_{n-1})}{h}]+O(h^2)\end{cases} ⎩⎪⎪⎪⎪⎨⎪⎪⎪⎪⎧dxdy∣xn=2hy(xn+1)−y(xn−1)+O(h2)dxd(p(x)dxdy)∣xn=h1((pdxdy)xn+21−(pdxdy)xn−21)+O(h2)=h1[pn+21hy(xn+1)−y(xn)−pn−21hy(xn)−y(xn−1)]+O(h2)
这里 xn+12=12(xn+1+xn),pn±12=p(xn±12)x_{n+\frac{1}{2}}=\dfrac{1}{2}(x_{n+1}+x_n),p_{n\pm\frac{1}{2}}=p(x_{n\pm\frac{1}{2}})xn+21=21(xn+1+xn),pn±21=p(xn±21)。

于是有：
−1h[pn+12y(xn+1)−y(xn)h−pn−12y(xn)−y(xn−1)h]+rny(xn+1)−y(xn−1)2h+qny(xn)+Rn=fnRn=O(h2)-\dfrac{1}{h}[p_{n+\frac{1}{2}}\dfrac{y(x_{n+1})-y(x_n)}{h}-p_{n-\frac{1}{2}}\dfrac{y(x_n)-y(x_{n-1})}{h}]+r_n\dfrac{y(x_{n+1})-y(x_{n-1})}{2h}+q_ny(x_n)+R_n=f_n\\R_n=O(h^2) −h1[pn+21hy(xn+1)−y(xn)−pn−21hy(xn)−y(xn−1)]+rn2hy(xn+1)−y(xn−1)+qny(xn)+Rn=fnRn=O(h2)
忽略高阶项 RnR_nRn，并用 yny_nyn 表示 y(xn)y(x_n)y(xn) 的近似值，得：
−1h[pn+12yn+1−ynh−pn−12yn−yn−1h]+rnyn+1−yn−12h+qnyn=fn,n=1,2,⋯,N−1边值条件：y0=α,yN=β-\dfrac{1}{h}[p_{n+\frac{1}{2}}\dfrac{y_{n+1}-y_n}{h}-p_{n-\frac{1}{2}}\dfrac{y_n-y_{n-1}}{h}]+r_n\dfrac{y_{n+1}-y_{n-1}}{2h}+q_ny_n=f_n,n=1,2,\cdots,N-1\\边值条件：y_0=\alpha,y_N=\beta −h1[pn+21hyn+1−yn−pn−21hyn−yn−1]+rn2hyn+1−yn−1+qnyn=fn,n=1,2,⋯,N−1边值条件：y0=α,yN=β
上述方程是关于未知量 yn(n=1,2,⋯,N−1)y_n(n=1,2,\cdots,N-1)yn(n=1,2,⋯,N−1) 的一个线性方程组，其系数矩阵为：
A=[b1−c1−a2b2−c2⋱⋱⋱−aN−2bN−2−cN−2−aN−1bN−1](N−1)×(N−1)bn=1h2(pn+12+pn−12)+qncn=1h2pn+12+12hrnan=1h2pn−12−12hrnA=\left[\begin{matrix}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{matrix}\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 A=⎣⎢⎢⎢⎢⎡b1−a2−c1b2⋱−c2⋱−aN−2⋱bN−2−aN−1−cN−2bN−1⎦⎥⎥⎥⎥⎤(N−1)×(N−1)bn=h21(pn+21+pn−21)+qncn=h21pn+21+2h1rnan=h21pn−21−2h1rn
当网格步长 hhh 充分小时，AAA 是一个对角占优的三对角矩阵，因此 AAA 是非奇异矩阵，可用追赶法直接求解。

若 p(x)≡1,r(x)≡0,f(x)=−γ(x)p(x)\equiv1,r(x)\equiv0,f(x)=-\gamma(x)p(x)≡1,r(x)≡0,f(x)=−γ(x)，则两点边值问题为：
{Ly=−y′′+q(x)y=−γ(x),q(x)≥0,a<x<by(a)=α,y(b)=βy′′(xi)−q(xi)y(xi)=γ(xi){yi+1−2yi+yi−1h2−qiyi=γiy0=α,yN=β,i=1,2,⋯,N−1\begin{cases}Ly=-y''+q(x)y=-\gamma(x),q(x)\geq0,a<x<b\\y(a)=\alpha,y(b)=\beta\end{cases}\\y''(x_i)-q(x_i)y(x_i)=\gamma(x_i)\\\begin{cases}\dfrac{y_{i+1}-2y_i+y_{i-1}}{h^2}-q_iy_i=\gamma_i\\y_0=\alpha,y_N=\beta\end{cases},i=1,2,\cdots,N-1 {Ly=−y′′+q(x)y=−γ(x),q(x)≥0,a<x<by(a)=α,y(b)=βy′′(xi)−q(xi)y(xi)=γ(xi)⎩⎨⎧h2yi+1−2yi+yi−1−qiyi=γiy0=α,yN=β,i=1,2,⋯,N−1
矩阵形式：
Ay=bA=[−2−q1h211−2−q2h21⋱⋱⋱1−2−qN−2h211−2−qN−1h2](N−1)×(N−1)y=[y1y2⋮yN−2yN−1],b=[γ1h2−αγ2h2⋮γN−2h2γN−1h2−β]Ay=b\\A=\left[\begin{matrix}-2-q_1h^2&1\\1&-2-q_2h^2&1\\&\ddots&\ddots&\ddots\\&&1&-2-q_{N-2}h^2&1\\&&&1&-2-q_{N-1}h^2\end{matrix}\right]_{(N-1)\times(N-1)}\\y=\left[\begin{matrix}y_1\\y_2\\\vdots\\y_{N-2}\\y_{N-1}\end{matrix}\right],b=\left[\begin{matrix}\gamma_1h^2-\alpha\\\gamma_2h^2\\\vdots\\\gamma_{N-2}h^2\\\gamma_{N-1}h^2-\beta\end{matrix}\right] Ay=bA=⎣⎢⎢⎢⎢⎡−2−q1h211−2−q2h2⋱1⋱1⋱−2−qN−2h211−2−qN−1h2⎦⎥⎥⎥⎥⎤(N−1)×(N−1)y=⎣⎢⎢⎢⎢⎢⎡y1y2⋮yN−2yN−1⎦⎥⎥⎥⎥⎥⎤,b=⎣⎢⎢⎢⎢⎢⎡γ1h2−αγ2h2⋮γN−2h2γN−1h2−β⎦⎥⎥⎥⎥⎥⎤

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

计算物理学（数值分析）上机实验答案5、常微分方程初值问题的数值解法
实验五.常微分方程初值问题的数值解法常微分方程的求解问题在实践中经常遇到,因此研究常微分方程的数值解法就很有必要.欧拉方法是最简单.最基本的方法,利用差商代替微商,就可得到一系列欧拉公式.这些公 ...
二阶边值问题的数值解matlab,《二阶常微分方程边值问题的数值解法》-毕业论文.doc...
w 摘要本文主要研究二阶常微分方程边值问题的数值解法.对线性边值问题,我们总结了两类常用的数值方法,即打靶法和有限差分方法,对每种方法都列出了详细的计算步骤和Matlab程序代码,通过具体的算例对 ...
matlab初值微分方程,常微分方程初值问题的MATLAB解法
大连圣亚海洋世界官网-2021年2月7日发(作者:转身之后还是你) 用 Matlab 求常微分方程 < br>(ODE) 的初值问题 (IVP) 本节考虑一阶常微分方程  u   f ...
数值计算方法上机c语言编程,数值计算方法上机实验报告.doc-资源下载在线文库www.lddoc.cn...
<数值计算方法>上机实验报告.doc 华北电力大学实验报告实验名称数值计算方法上机实验课程名称数值计算方法专业班级电力实 08 学生姓名李超然学号20080100 ...
线性代数方程组数值解法
线性代数方程组的数值解法线性代数方程组数值解法一.向量范数与矩阵范数 1.1 向量范数 1.1.1 满足三个条件(向量范数公理) 1.1.2 常用的向量范数 1.2 矩阵范数 1.2.1 矩阵范数 ...
【数理知识】《数值分析》李庆扬老师-第9章-常微分方程初值问题数值解法
第8章回到目录无第9章-常微分方程初值问题数值解法 9.1 引言利普希茨 (Lipschitz) 条件 / 利普希茨常数定理1 解的存在唯一性定理定理2 解对初值依赖的敏感性 9.2 简单 ...
常微分方程初值问题数值解法[完整公式]（Python）
目录 1.概述 (1)常微分方程初值问题数值解法 (2)解题步骤 (3)数值微分解法 (4)数值积分解法 2.所有知识点代码 3.结果---以三阶Runge-Kutta公式为例(其他的类似) 1.概述 ...
常微分方程数值解的c语言程序,常微分方程的数值解法一阶常微分方程数值解的C语言编程实现.doc...
常微分方程的数值解法一阶常微分方程数值解的C语言编程实现导读:就爱阅读网友为您分享以下"一阶常微分方程数值解的C语言编程实现"资讯,希望对您有所帮助,感谢您对92的支持! 一阶 ...
【matlab】常微分方程的数值解法
实验任务 (一)常微分方程的符号计算和数值解法基本操作 1.课上例题:1.2.3.4 (二)专题实验(梯形格式) 编写梯形公式的程序,要求: 程序要有通用性,例如: functio ...

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