数值分析复习（六）——常微分方程数值解法

六、常微分方程数值解法

本章讨论一阶常微分方程初值问题的数值解
{y′=f(x,y),a≤x≤by(a)=η\left\{ \begin{aligned} &y^{'} = f(x, y),\ a \le x \le b \\ &y(a) = \eta \end{aligned} \right. {y′=f(x,y), a≤x≤by(a)=η
假设

f(x,y)f(x,y)f(x,y)，∂f(x,y)∂y\frac{\partial f(x, y)}{\partial y}∂y∂f(x,y) 连续
上式存在唯一解 y(x)y(x)y(x) 且在 [a,b][a,b][a,b] 上充分光滑

离散化：将 [a,b][a,b][a,b] 作 nnn 等分，记 h=b−anh=\frac{b-a}{n}h=nb−a，xi=a+ih,(i=0,1,...,n)x_i=a+ih,\ (i=0,1,...,n)xi=a+ih, (i=0,1,...,n)。称 hhh 为步长，数值解为求初值问题的解 y(x)y(x)y(x) 在离散点 xix_ixi 处的近似值 yiy_iyi

计算 yi+1y_{i+1}yi+1 时，如果只用到前一步的值 yiy_iyi，称这类方法为单步法
计算 yi+1y_{i+1}yi+1 时，如果用到前 r 步的值 yi,yi−1,...,yi−r+1y_i,y_{i-1},...,y_{i-r+1}yi,yi−1,...,yi−r+1，这类方法称为 r 步方法

Euler 公式

Euler 公式：

yi+1=yi+hf(xi,yi),i=0,1,...,n−1y_{i+1}=y_i+hf(x_i,y_i),\ i=0,1,...,n-1yi+1=yi+hf(xi,yi), i=0,1,...,n−1（单步显式公式）
局部截断误差：Ri+1=y(xi+1)−[y(xi)+hf(xi,y(xi))]=12h2y′′(ξi),ξi∈(xi,xi+1)R_{i+1}=y(x_{i+1})-[y(x_i)+hf(x_i,y(x_i))]=\frac{1}{2}h^2y^{''}(\xi_i),\xi_i \in (x_i, x_{i+1})Ri+1=y(xi+1)−[y(xi)+hf(xi,y(xi))]=21h2y′′(ξi),ξi∈(xi,xi+1)
若 Ri+1=O(hp+1)R_{i+1}=O(h^{p+1})Ri+1=O(hp+1)，则称该公式是 p 阶的，或具有 p 阶精度。Euler 公式是 1 阶的

后退 Euler 公式：

yi+1=yi+hf(xi+1,yi+1),i=0,1,..,n−1y_{i+1}=y_i+hf(x_{i+1}, y_{i+1}),\ i=0,1,..,n-1yi+1=yi+hf(xi+1,yi+1), i=0,1,..,n−1（单步隐式公式）
局部截断误差：Ri+1=y(xi+1)−y(xi)−hf(xi+1,y(xi+1))=−h22y′′(ξi),ξi∈(xi,xi+1)R_{i+1}=y(x_{i+1})-y(x_i)-hf(x_{i+1},y(x_{i+1}))=-\frac{h^2}{2}y^{''}(\xi_i),\xi_i \in (x_i, x_{i+1})Ri+1=y(xi+1)−y(xi)−hf(xi+1,y(xi+1))=−2h2y′′(ξi),ξi∈(xi,xi+1)
后退 Euler 公式是 1 阶的

利用梯形公式的 Euler 公式：

yi+1=yi+h2[f(xi,yi)+f(xi+1,yi+1)],i=0,1,...,n−1y_{i+1}=y_i+\frac{h}{2}[f(x_i,y_i)+f(x_{i+1}, y_{i+1})],\ i=0,1,...,n-1yi+1=yi+2h[f(xi,yi)+f(xi+1,yi+1)], i=0,1,...,n−1（单步隐式公式）
局部截断误差：Ri+1=y(xi+1)−y(xi)−h2[f(xi,y(xi))+f(xi+1,y(xi+1))]=−112y′′′(ξi)h3,ξi∈(xi,xi+1)R_{i+1}=y(x_{i+1})-y(x_i)-\frac{h}{2}[f(x_i,y(x_i))+f(x_{i+1},y(x_{i+1}))]=-\frac{1}{12}y^{'''}(\xi_i)h^3,\xi_i \in (x_i, x_{i+1})Ri+1=y(xi+1)−y(xi)−2h[f(xi,y(xi))+f(xi+1,y(xi+1))]=−121y′′′(ξi)h3,ξi∈(xi,xi+1)
该 Euler 公式是 2 阶的

改进 Euler 公式：
{yi+1(p)=yi+hf(xi,yi)预测公式yi+1=yi+h2[f(xi,yi)+f(xi+1,yi+1(p))]校正公式\left\{ \begin{aligned} y_{i+1}^{(p)} &= 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}^{(p)})]\ 校正公式 \end{aligned} \right. ⎩⎪⎨⎪⎧yi+1(p)yi+1=yi+hf(xi,yi) 预测公式=yi+2h[f(xi,yi)+f(xi+1,yi+1(p))] 校正公式

该公式是单步显式公式，其阶为 2 阶
局部误差：Ri+1=y(xi+1)−y(xi)−{h2[f(xi,y(xi))+f(xi+1,y(xi)+hf(xi,y(xi)))]}R_{i+1}=y(x_{i+1})-y(x_i)-\{\frac{h}{2}[f(x_i,y(x_i))+f(x_{i+1},y(x_i)+hf(x_i, y(x_i)))]\}Ri+1=y(xi+1)−y(xi)−{2h[f(xi,y(xi))+f(xi+1,y(xi)+hf(xi,y(xi)))]}，计算可通过泰勒展开来处理

整体截断误差：

设当步长为 hhh 时某种数值方法求得的数值解为 y1[h],y2[h],...,yn[h]y_1^{[h]},y_2^{[h]},...,y_n^{[h]}y1[h],y2[h],...,yn[h]，其中 y(xi),yi[h],i=1,2,...,ny(x_i),y_i^{[h]},i=1,2,...,ny(xi),yi[h],i=1,2,...,n，分别为精确解和数值解，则称
E(h)=max⁡i≤i≤n∣y(xi)−yi[h]∣E(h) = \max_{i \le i \le n}|y(x_i)-y_i^{[h]}| E(h)=i≤i≤nmax∣y(xi)−yi[h]∣
为该数值方法的整体截断误差。若 lim⁡h→0E(h)=0\lim_{h\rightarrow 0}E(h)=0limh→0E(h)=0 则称该方法收敛

Runge-Kutta

构造思想：

由 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))dxy(xi+1)=y(xi)+∫xixi+1f(x,y(x))dx 可得
y(xi+1)=y(xi)+hf(xi+θh,y(xi+θh)),θ∈(0,1)y(x_{i+1})=y(x_i)+hf(x_i+\theta h,y(x_i+\theta h)),\ \theta \in (0, 1) y(xi+1)=y(xi)+hf(xi+θh,y(xi+θh)), θ∈(0,1)
称 f(xi+θh,y(xi+θh))f(x_i+\theta h,y(x_i+\theta h))f(xi+θh,y(xi+θh)) 为 y(x)y(x)y(x) 在 [xi,xi+1][x_i, x_{i+1}][xi,xi+1] 上的平均斜率，记为 k∗k^*k∗

记 k1=f(xi,yi)k_1=f(x_i,y_i)k1=f(xi,yi)，k2=f(xi+1,yi+hk1)k_2=f(x_{i+1},y_i+hk_1)k2=f(xi+1,yi+hk1)，若用 k1k_1k1 近似 k∗k^*k∗，则得一阶 Euler 公式；若用 k1+k22\frac{k_1+k_2}{2}2k1+k2 近似 k∗k^*k∗，则得二阶改进的 Euler 公式

一般 r 级 Runge-Kutta 公式：
{yi+1=yi+h∑j=1rαjkjk1=f(xi,yi)kj=f(xi+λjh,yi+h∑l=1j−1μjlkl),j=2,3,...,r\left\{ \begin{aligned} y_{i+1} &= y_i+h\sum_{j=1}^r\alpha_jk_j \\ k_1 &= f(x_i, y_i) \\ k_j &= f(x_i+\lambda_jh,y_i+h\sum_{l=1}^{j-1}\mu_{jl}k_l),\ j=2,3,...,r \end{aligned} \right. ⎩⎪⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎪⎧yi+1k1kj=yi+hj=1∑rαjkj=f(xi,yi)=f(xi+λjh,yi+hl=1∑j−1μjlkl), j=2,3,...,r
选择参数 αj\alpha_jαj，λj\lambda_jλj，μjl\mu_{jl}μjl 使其具有一定的阶数

2 阶 Runge-Kutta 公式：
(式0){yi+1=yi+h(α1k1+α2k2)k1=f(xi,yi)k2=f(xi+λ2h,yi+hμ21k1)(式0)\left\{ \begin{aligned} y_{i+1} &= y_i+h(\alpha_1k_1+\alpha_2k_2) \\ k_1 &= f(x_i, y_i) \\ k_2 &= f(x_i+\lambda_2h,y_i+h\mu_{21}k_1) \end{aligned} \right. (式0)⎩⎪⎨⎪⎧yi+1k1k2=yi+h(α1k1+α2k2)=f(xi,yi)=f(xi+λ2h,yi+hμ21k1)
其截断误差为：
{Ri+1=y(xi+1)−y(xi)−h(α1K1+α2K2)K1=f(xi,y(xi))K2=f(xi+λ2h,y(xi)+hμ21K1)\left\{ \begin{aligned} R_{i+1} &= y(x_{i+1})-y(x_i)-h(\alpha_1K_1+\alpha_2K_2) \\ K_1 &= f(x_i, y(x_i)) \\ K_2 &= f(x_i+\lambda_2h, y(x_i)+h\mu_{21}K_1) \end{aligned} \right. ⎩⎪⎨⎪⎧Ri+1K1K2=y(xi+1)−y(xi)−h(α1K1+α2K2)=f(xi,y(xi))=f(xi+λ2h,y(xi)+hμ21K1)
注：可通过泰勒展开来处理上述截断误差

要使式 0 具有 2 阶精度，则
{1−α1−α2=012−α2λ2=012−α2μ21=0\left\{ \begin{aligned} &1- \alpha_1 - \alpha_2 = 0 \\ &\frac{1}{2} - \alpha_2 \lambda_2 = 0 \\ &\frac{1}{2} - \alpha_2 \mu_{21} = 0 \end{aligned} \right. ⎩⎪⎪⎪⎪⎨⎪⎪⎪⎪⎧1−α1−α2=021−α2λ2=021−α2μ21=0
显然 α2\alpha_2α2 不能为 0，则
{α1=1−α2λ2=12α2μ21=12α2\left\{ \begin{aligned} &\alpha_1 = 1 - \alpha_2 \\ &\lambda_2 = \frac{1}{2\alpha_2} \\ &\mu_{21} = \frac{1}{2\alpha_2} \end{aligned} \right. ⎩⎪⎪⎪⎪⎨⎪⎪⎪⎪⎧α1=1−α2λ2=2α21μ21=2α21
从而可得一类 2 阶 Runge-Kutta 公式
{yi+1=yi+h[(1−α2)k1+α2k2]k1=f(xi,yi)k2=f(xi+12α2h,yi+12α2hk1)\left\{ \begin{aligned} &y_{i+1} = y_i + h[(1-\alpha_2)k_1+\alpha_2k_2] \\ &k_1 = f(x_i, y_i) \\ &k_2 = f(x_i+\frac{1}{2\alpha_2}h, y_i+\frac{1}{2\alpha_2}hk_1) \end{aligned} \right. ⎩⎪⎪⎪⎨⎪⎪⎪⎧yi+1=yi+h[(1−α2)k1+α2k2]k1=f(xi,yi)k2=f(xi+2α21h,yi+2α21hk1)
当 α2=12\alpha_2=\frac{1}{2}α2=21，得改进的 Euler 公式。当 α2=1\alpha_2=1α2=1，得变形的 Euler 公式

单步法的收敛性和稳定性

考虑单步显式公式
(式1){yi+1=yi+hφ(xi,yi,h),i=0,1,...,n−1y0=η(式1)\ \left\{ \begin{aligned} &y_{i+1} = y_i + h\varphi(x_i,y_i,h),\ i=0,1,...,n-1 \\ &y_0 = \eta \end{aligned} \right. (式1) {yi+1=yi+hφ(xi,yi,h), i=0,1,...,n−1y0=η

本章讨论一阶常微分方程初值问题的数值解
(式2){y′=f(x,y),a≤x≤by(a)=η(式2)\ \left\{ \begin{aligned} &y^{'} = f(x,y),\ a \le x \le b \\ &y(a) = \eta \end{aligned} \right. (式2) {y′=f(x,y), a≤x≤by(a)=η
收敛性：

设 y(x)y(x)y(x) 是上述微分方程的解（式2），{yi}i=0n\{y_i\}_{i=0}^n{yi}i=0n 为单步显式公式（式1）的解。如果

存在常数 c0>0c_0>0c0>0，使得 ∣Ri+1∣≤c0hp+1,i=0,1,...,n−1|R_{i+1}|\le c_0h^{p+1},\ i=0,1,...,n-1∣Ri+1∣≤c0hp+1, i=0,1,...,n−1 => 至少 ppp 阶
存在 h0>0,L>0h_0 \gt 0,L \gt 0h0>0,L>0，使得 max⁡(x,y)∈Dδ,0<h≤h0∣∂φ(x,y,h)∂y∣≤L\max_{(x,y)\in D_\delta,0 \lt h \le h_0}|\frac{\partial\varphi(x,y,h)}{\partial y}| \le Lmax(x,y)∈Dδ,0<h≤h0∣∂y∂φ(x,y,h)∣≤L => 有界

则当 h≤min⁡{h0,δcp}h \le \min\{h_0, \sqrt[p]{\frac{\delta}{c}}\}h≤min{h0,pcδ} 时，有 E(h)≤chpE(h) \le ch^pE(h)≤chp

其中 Dδ={(x,y)∣a≤x≤b,y(x)−δ≤y≤y(x)+δ}D_\delta = \{(x, y)|a \le x \le b,y(x)-\delta \le y \le y(x)+\delta \}Dδ={(x,y)∣a≤x≤b,y(x)−δ≤y≤y(x)+δ}，c=c0L[eL(b−a)−1]c=\frac{c_0}{L}[e^{L(b-a)}-1]c=Lc0[eL(b−a)−1]

稳定性：

对于上述初值问题（式2），设 {yi}i=0n\{y_i\}_{i=0}^n{yi}i=0n 是由上述单步法（式1）得到得近似解，{zi}i=0n\{z_i\}_{i=0}^n{zi}i=0n 是单步法（式1）扰动后的解，即满足
{zi+1=zi+h[φ(xi,yi,h)+δi+1],i=0,1,...,n−1z0=η+δ0\left\{ \begin{aligned} &z_{i+1} = z_i + h[\varphi(x_i,y_i,h)+\delta_{i+1}],\ i=0,1,...,n-1 \\ &z_0 = \eta + \delta_0 \end{aligned} \right. {zi+1=zi+h[φ(xi,yi,h)+δi+1], i=0,1,...,n−1z0=η+δ0
如果存在正常数 CCC，ϵ0\epsilon_0ϵ0，h0h_0h0，使得对所有 ϵ∈(0,ϵ0]\epsilon \in (0,\epsilon_0]ϵ∈(0,ϵ0]，h∈(0,h0]h \in (0, h_0]h∈(0,h0]，当 max⁡0≤i≤n∣δi∣≤ϵ\max_{0 \le i \le n}|\delta_i| \le \epsilonmax0≤i≤n∣δi∣≤ϵ 时，有
max⁡0≤i≤n∣yi−zi∣≤Cϵ\max_{0 \le i \le n}|y_i-z_i| \le C \epsilon 0≤i≤nmax∣yi−zi∣≤Cϵ
则称单步法（式1）稳定。在满足收敛性的条件下，单步公式（式1）是稳定的

线性多步法

一般的线性 k 步方法为
yi+1=∑j=0k−1ajyi−j+h∑j=−1k−1bjf(xi−j,yi−j)y_{i+1}=\sum_{j=0}^{k-1}a_jy_{i-j}+h\sum_{j=-1}^{k-1}b_jf(x_{i-j},y_{i-j}) yi+1=j=0∑k−1ajyi−j+hj=−1∑k−1bjf(xi−j,yi−j)
其中 ak−1a_{k-1}ak−1，bk−1b_{k-1}bk−1 不同时为零。当 b−1=0b_{-1}=0b−1=0 时为显式公式；当 b−1≠0b_{-1} \ne 0b−1=0 时为隐式公式

Adams 公式

Adams 显式公式
- yi+1=yi+h∑j=0rβrjf(xi−j,yi−j)y_{i+1} = y_i + h\sum_{j=0}^r \beta_{rj}f(x_{i-j}, y_{i-j})yi+1=yi+h∑j=0rβrjf(xi−j,yi−j)
- 截断误差：Ri+1=αr+1hr+2y(r+2)(ξi)R_{i+1}=\alpha_{r+1}h^{r+2}y^{(r+2)}(\xi_i)Ri+1=αr+1hr+2y(r+2)(ξi)
- (r+1) 步、(r+1) 阶显式的 Adams 公式
- 当 r=0r=0r=0 时得 Euler 公式 yi+1=yi+hf(xi,yi)y_{i+1} = y_i+hf(x_i,y_i)yi+1=yi+hf(xi,yi)

其中 βrj\beta_{rj}βrj 和 αr+1\alpha_{r+1}αr+1 定义如下
βrj=∫01∏l=0,l≠jrl+tl−jdt,j=0,1,...,rαr+1=1(r+1)!∫01∏j=0r(j+t)dt\begin{aligned} &\beta_{rj} = \int_0^1 \prod_{l=0,l \ne j}^r \frac{l+t}{l-j}dt,\ j=0,1,...,r \\ &\alpha_{r+1} = \frac{1}{(r+1)!}\int_0^1 \prod_{j=0}^r (j+t) dt \end{aligned} βrj=∫01l=0,l=j∏rl−jl+tdt, j=0,1,...,rαr+1=(r+1)!1∫01j=0∏r(j+t)dt

Adams 隐式公式
- yi+1=yi+h∑j=−1r−1β‾rjf(xi−j,yi−j)y_{i+1} = y_i + h\sum_{j=-1}^{r-1} \overline \beta_{rj}f(x_{i-j}, y_{i-j})yi+1=yi+h∑j=−1r−1βrjf(xi−j,yi−j)
- 截断误差：Ri+1=α‾r+1hr+2y(r+2)(ξ‾i)R_{i+1}=\overline \alpha_{r+1}h^{r+2}y^{(r+2)}(\overline \xi_i)Ri+1=αr+1hr+2y(r+2)(ξi)
- r 步、(r+1) 阶隐式的 Adams 公式
- 当 r=1r=1r=1 时得梯形公式 yi+1=yi+h2[f(xi+1,yi+1)+f(xi,yi)]y_{i+1} = y_i+\frac{h}{2}[f(x_{i+1}, y_{i+1})+f(x_i, y_i)]yi+1=yi+2h[f(xi+1,yi+1)+f(xi,yi)]

其中 β‾rj\overline{\beta}_{rj}βrj 和 α‾r+1\overline{\alpha}_{r+1}αr+1 定义如下
β‾rj=∫01∏l=−1,l≠jr−1l+tl−jdt,j=0,1,...,rα‾r+1=1(r+1)!∫01∏j=−1r−1(j+t)dt\begin{aligned} &\overline{\beta}_{rj} = \int_0^1 \prod_{l=-1,l \ne j}^{r-1} \frac{l+t}{l-j}dt,\ j=0,1,...,r \\ &\overline{\alpha}_{r+1} = \frac{1}{(r+1)!} \int_0^1 \prod_{j=-1}^{r-1}(j+t)dt \end{aligned} βrj=∫01l=−1,l=j∏r−1l−jl+tdt, j=0,1,...,rαr+1=(r+1)!1∫01j=−1∏r−1(j+t)dt

Adams 预测校正公式：将同阶（以下是 2 阶公式结合）的显式 Adams 公式和隐式 Adams 公式结合起来

{yi+1(p)=yi+h2[3f(xi,yi)−f(xi−1,yi−1)]yi+1=yi+h2[f(xi+1,yi+1(p))+f(xi,yi)]\left\{ \begin{aligned} y_{i+1}^{(p)} &= y_i + \frac{h}{2}[3f(x_i, y_i)-f(x_{i-1}, y_{i-1})] \\ y_{i+1} &= y_i + \frac{h}{2}[f(x_{i+1},y_{i+1}^{(p)})+f(x_i, y_i)] \end{aligned} \right. ⎩⎪⎪⎨⎪⎪⎧yi+1(p)yi+1=yi+2h[3f(xi,yi)−f(xi−1,yi−1)]=yi+2h[f(xi+1,yi+1(p))+f(xi,yi)]

基于 Taylor 展开得待定系数法

若有 f(x,y(x))f(x,y(x))f(x,y(x))，可考虑将其转化为 y′(x)y^{'}(x)y′(x)
利用泰勒展开，将式子展开，如下

y(xi−j)=y(xi−jh)=y(xi)+(−jh)y′(xi)+(−jh)22!y′′(xi)+...\begin{aligned} &y(x_{i-j}) = y(x_i - jh) = y(x_i) + (-jh)y^{'}(x_i) + \frac{(-jh)^2}{2!}y^{''}(x_i)+... \\ \end{aligned} y(xi−j)=y(xi−jh)=y(xi)+(−jh)y′(xi)+2!(−jh)2y′′(xi)+...

合并阶数相同的项
若要求得的阶为 ppp，则令阶小于 ppp 的项的系数为 0，然后求解相应的系数即可