常微分方程的数值解法

基本原理

符号定义

符号 含义
y(xn)y(x_n) xnx_n处精确值，无误差
yny_n xnx_n处理论计算值，只有截断误差影响
yˉn\bar{y}_n xnx_n处实际计算值，还有舍入误差影响
En=y(xn)−ynE_n = y(x_n)-y_n 局部截断误差，En=O(hP+1)E_n = O(h^{P+1})有PP阶精度。

一阶常微分方程

{dydx=f(x,y)a≤x≤by(x0)=y0\begin{cases} \dfrac{dy}{dx}=f(x,y) & a≤x≤b \\ y(x_0)=y_0 \end{cases}

李普希茨（Lipschitz）条件

如果函数f(x,y)f(x,y)在区域{a≤x≤b,m≤y≤M}\{ a≤x≤b, m≤y≤M \}上连续，且关于yy满足Lipschitz条件：

∣f(x,y1)−f(x,y2)∣≤L∣y1−y2∣|f(x,y_1)-f(x,y_2)| ≤ L|y_1-y_2|

其中LL为一个正数称为Lipschitz常数。若满足Lipschitz条件，则一阶常微分方程的解存在且唯一。若不满足，则未必存在唯一解。

基本差分公式

y(xn+1)−y(xn)=∫xnxn+1f(x,y(x))dx   ①y(x_{n+1})-y(x_n) = \int_{x_n}^{x_{n+1}} f(x, y(x)) {\rm d}x { \ \ \ } ①

y(xn+1)−y(xn−1)=∫xn−1xn+1f(x,y(x))dx   ②y(x_{n+1})-y(x_{n-1}) = \int_{x_{n-1}}^{x_{n+1}} f(x, y(x)) {\rm d}x { \ \ \ } ②

单步法收敛性

对于一阶常微分方程，其单步方法可写成：yn+1=yn+hΦ(xn,yn,h)y_{n+1} = y_n + hΦ(x_n, y_n, h)

定义1

lim⁡h→0yn=y(xn)\lim_{h→0} y_n = y(x_n)

定义2

lim⁡k→∞yn+1[k]=yn+1\lim_{k→∞} y_{n+1}^{[k]} = y_{n+1}

单步法稳定性

定义1

对于给定的初值问题，取定步长hh，用差分方法进行计算时，假设只在一个节点值yny_n上产生计算误差σσ，如果这个误差引起的以后各节点计算值ym,(m>n)y_m, (m>n)的变化均不超过σσ，则称此差分方法是绝对稳定的。

通常仅限于实验方程

dxdy=λy\dfrac{ {\rm d}x }{ {\rm d}y } = λy

其中λλ是复数，且Re(λ)<0Re(λ)<0（即λλ实部<0<0）

定义2

将单步方法用于解试验方程，假设得到

yn+1yn=E(λh)\dfrac{y_{n+1}}{y_{n}} = E(λh)

若满足条件∣E(λh)∣<1|E(λh)|<1则称此单步方法是绝对稳定的。否则是不稳定的。

在复平面上，变量λhλh满足∣E(λh)∣<1|E(λh)|<1的区域称为此方法的绝对稳定域，它与实轴的交集称为绝对稳定区间。

计算稳定性

若满足

∣ρn+1∣∣ρn∣<1\dfrac{ |ρ_{n+1}| }{ |ρ_{n}| } < 1

则说明计算中摄入误差可以得到控制，数值方法是稳定的。

一般仅考虑模型方程y′=λyy' = λy，要求λ<0λ<0或Re(λ)<0Re(λ)<0。

例

y′=−10yy'=-10y，y(0)=1y(0)=1，当用Euler方法求解时，hh取值范围为多少才能使计算稳定。

解

f(x,y)=−10yf(x,y) = -10y，构造Euler公式得yn+1=yn+h(−10yn)=(1−10h)yny_{n+1} = y_n + h(-10y_n) = (1-10h)y_n

则

∣ρn+1ρn∣=∣yn+1∗−yn+1yn∗−yn∣=∣1−10h∣<1\left| \dfrac{ρ_{n+1}}{ρ_{n}} \right| = \left| \dfrac{ y_{n+1}^* - y_{n+1} }{ y_{n}^* - y_{n} } \right| = |1-10h| < 1

解得0<h<0.20<h<0.2

欧拉（Euler）方法

欧拉公式

由基本差分公式①和左矩形近似（∫abf(x)dx≈(b−a)f(a)\int_{a}^{b} f(x) {\rm d}x ≈ (b-a)f(a)）得

yn+1=yn+hf(xn,yn)   n=0,1,2,⋯y_{n+1} = y_n + hf(x_n, y_n) { \ \ \ } n=0,1,2,\cdots

1阶精度。

局部截断误差

y(xn+1)−yn+1=y′′(xn)2!h2+O(h3)=O(h2)y(x_{n+1}) - y_{n+1} = \dfrac{y''(x_n)}{2!}h^2 + O(h^3) = O(h^2)，故欧拉公式为一阶方法。

例

一阶常微分方程：{dydx=f(x,y)a≤x≤by(x0)=y0\begin{cases} \dfrac{dy}{dx}=f(x,y) & a≤x≤b \\ y(x_0)=y_0 \end{cases}

求{y′=−y+x+1x≥0y(0)=1\begin{cases} y'=-y+x+1 & x≥0 \\ y(0)=1 \end{cases}的初值，取h=0.1h=0.1，计算到x=0.5x=0.5

解

f(x)=−y+x+1f(x)=-y+x+1 =x+e−x=x+e^{-x}，由欧拉公式得yn+1=yn+h(−yn+xn+1)y_{n+1} = y_n + h(-y_n + x_n + 1)，即

{xn+1=xn+hyn+1=0.9yn+0.1(xn+1)\begin{cases} x_{n+1} = x_n + h \\ y_{n+1} = 0.9y_n + 0.1(x_n + 1) \end{cases}

nn xnx_n yny_n
0 0 1
1 0.1 1.0
2 0.2 1.01
3 0.3 1.029
4 0.4 1.0561
5 0.5 1.09049

梯形公式

由基本差分公式①和梯形近似（∫abf(x)dx≈b−a2[f(a)+f(b)]\int_{a}^{b} f(x) {\rm d}x ≈ \dfrac{b-a}{2}\left[ {f(a) + f(b)} \right]）得

yn+1=yn+h2[f(xn,yn)+f(xn+1,yn+1)]   n=0,1,2,⋯y_{n+1} = y_n + \dfrac{h}{2}[ { f(x_n, y_n) + f(x_{n+1}, y_{n+1}) }] { \ \ \ } n=0,1,2,\cdots

中点欧拉公式

由基本差分公式②和中矩形近似（∫abf(x)dx≈(b−a)f(a+b2)\int_{a}^{b} f(x) {\rm d}x ≈ (b-a)f(\dfrac{a+b}{2})）得

yn+1=yn−1+2hf(xn,yn)   n=1,2,⋯y_{n+1} = y_{n-1} +2hf(x_n, y_n) { \ \ \ } n=1,2,\cdots

改进的欧拉公式

yn+1=yn+h2[f(xn,yn)+f(xn+1,yn+hf(xn,yn))]   n=0,1,2,⋯y_{n+1} = y_n + \dfrac{h}{2}[ f(x_n, y_n) + f(x_{n+1}, y_n + hf(x_n,y_n)) ] { \ \ \ } n=0,1,2,\cdots

2阶精度。

局部截断误差

y(xn+1)−yn+1=O(h3)y(x_{n+1}) - y_{n+1} = O(h^3)，故改进的欧拉公式为二阶方法。

龙格库塔（Runge-Kutta）方法

pp级方法

{yn+1=yn+h(∑i=1pλiKi)K1=f(xn,yn)K2=f(xn+α2h,yn+hβ21K1)⋯ ,⋯Kp=f(xn+aph,yn+h(∑i=1p−1βpiKi))\begin{cases} y_{n+1} = y_{n} + h(\sum\limits_{i=1}^{p} λ_iK_i) \\ K_1 = f(x_n, y_n) \\ K_2 = f(x_n + α_2h, y_n + hβ_{21}K_1) \\ \cdots, \cdots \\ K_p = f(x_n + a_ph, y_n + h(\sum\limits_{i=1}^{p-1} β_{pi}K_i)) \end{cases}

其中{λi,αi,βij}\{λ_i, α_i, β_{ij}\}为待定参数，此公式称为：pp级Runge-Kutta方法。

R-K方法

若pp级Runge-Kutta方法的局部截断误差为O(hp+1)O(h^{p+1})，则称为pp级R-K方法。

当p=2p=2时

{yn+1=yn+h(λ1K1+λ2K2)K1=f(xn,yn)K2=f(xn+αh,yn+hβK1)\begin{cases} y_{n+1} = y_{n} + h(λ_1K_1 + λ_2K_2) \\ K_1 = f(x_n, y_n) \\ K_2 = f(x_n + αh, y_n + hβK_1) \end{cases}

若取α=1,λ1=λ2=12,β=1α=1, λ_1=λ_2=\frac{1}{2}, β=1，则此时该公式为：改进的Euler方法。

若取α=β=12,λ1=0，λ2=1α=β=\frac{1}{2}, λ_1=0， λ_2=1，则此时该公式为：中点公式。

线性多步方法

k+1k+1步线性多部方法的一般形式为

yn+1=∑i=0kαiyn−i+h∑i=−1kβifn−iy_{n+1} = \sum_{i=0}^{k} α_iy_{n-i} + h\sum_{i=-1}^{k}β_if_{n-i}

β−1=0β_{-1}=0 β−1≠0β_{-1}≠0
显式公式 隐式公式

局部截断误差

y(xn+1)−yn+1=O(hk+2)y(x_{n+1})-y_{n+1} = O(h^{k+2})

外插公式

阿达姆斯（Adams）显式公式、四阶方法

设已求得精确解y(x)y(x)在步长为hh的前k+1k+1个等距节点xn−k,⋯ ,xnx_{n-k},\cdots,x_n上的近似值yn−k,⋯ ,yny_{n-k},\cdots,y_n，构造kk次Lagrange插值多项式

Pk(x)=∑i=0kln−i(x)fn−iP_k(x) = \sum_{i=0}^{k} l_{n-i}(x)f_{n-i}

由基本差分公式①和f(x,y)=Pk(x)+Rk(x)f(x,y) = P_k(x) + R_k(x)得k+1k+1步Adams显式公式

yn+1=yn+h∑i=0kβkifn−iy_{n+1} = y_n + h\sum_{i=0}^{k} β_{ki}f_{n-i}

当k=3k=3时，即四阶Adams显式公式为

yn+1=yn+h24(55fn−59fn−1+37fn−2−9fn−3)y_{n+1} = y_n + \dfrac{h}{24}( 55f_{n} - 59f_{n-1} + 37f_{n-2} - 9f_{n-3} )

内插公式

阿达姆斯（Adams）隐式公式、四阶方法

如果利用k+1k+1个数据(xn−k+1,fn−k+1),⋯ ,(xn+1,fn+1)(x_{n-k+1},f_{n-k+1}),\cdots,(x_{n+1},f_{n+1})构造kk次Lagrange插值多项式Pk(x)P_k(x)，则可导出数值稳定性好的Adams隐式公式

yn+1=yn+h∑i=0kβki∗fn−i+1y_{n+1} = y_n + h\sum_{i=0}^{k} β_{ki}^*f_{n-i+1}

当k=3k=3时，即四阶Adams隐式公式

yn+1=yn+h24(9fn+1+19fn−5fn−1+fn−2)y_{n+1} = y_n + \dfrac{h}{24}( 9f_{n+1} + 19f_{n} - 5f_{n-1} + f_{n-2} )