方程组及高阶方程的数值方法

常微分方程的数值方法

张瑞
中国科学技术大学数学科学学院

一阶常微分方程组的初值问题

一阶常微分方程组可以写成如下形式

\[\begin{cases} \frac{dy_1}{dx}=f_1(x,y_1,\cdots,y_m) \\ \frac{dy_2}{dx}=f_2(x,y_1,y_2,\cdots,y_m) \\ \cdots \\ \frac{dy_m}{dx}=f_m(x,y_1,y_2,\cdots,y_m) \\ y_1(a)=\eta_1 \\ y_2(a)=\eta_2 \\ \cdots \\ y_m(a)=\eta_m \end{cases} , a\leq x\leq b \]

写成向量的形式：

\[\begin{cases} \dfrac{dY}{dx}=F(x,Y) \\ Y(a)=\eta \end{cases} , x\in[a,b] \]

其中$Y=(y_1,y_2,\cdots, y_m)^T$, $F=(f_1,f_2, \cdots, f_m)^T$。

把各种方法直接应用过来

例 1. 构造如下方程组的数值格式

\[\begin{cases} \frac{dy}{dx}=f(x,y,z) \\ \frac{dz}{dx}=g(x,y,z) \\ y(a)=y_0 \\ z(a)=z_0 \end{cases} \]

解. 记$Y=(y,z)^T$, $F(x,Y)=(f(x,Y),g(x,Y))^T$，则

向前Euler格式

\[Y_{n+1}=Y_n+hF(x_n,Y_n) \]

即

\[\begin{pmatrix} y_{n+1} \\z_{n+1} \end{pmatrix} =\begin{pmatrix}y_{n} \\z_{n} \end{pmatrix} +h\begin{pmatrix}f(x_n,y_{n},z_n) \\g(x_n,y_n,z_{n}) \end{pmatrix} \]

可以得到格式为

\[\begin{cases} y_{n+1}=y_n+hf(x_n,y_n,z_n) \\ z_{n+1}=z_n+hg(x_n,y_n,z_n) \end{cases} \]

Runge-Kutta格式

\[\begin{pmatrix} y_{n+1} \\z_{n+1} \end{pmatrix} =\begin{pmatrix}y_{n} \\z_{n} \end{pmatrix} +\frac{h}6(K_1+2K_2+2K_3+K_4) \]

其中

\[K_1=F(x_n+Y_n)=\begin{pmatrix} f(x_n,y_n,z_n) \\ g(x_n,y_n,z_n)\end{pmatrix} =\begin{pmatrix} K_1^{(1)} \\ K_1^{(2)} \end{pmatrix} \]
\[\begin{aligned} K_2=&F(x_n+\frac{h}2,Y_n+\frac{h}2K_1) \\ =&\begin{pmatrix} f(x_n+\frac{h}2,y_n+\frac{h}2K_1^{(1)},z_n+\frac{h}2K_1^{(2)} \\ g(x_n+\frac{h}2,y_n+\frac{h}2K_1^{(1)},z_n+\frac{h}2K_1^{(2)} \end{pmatrix} \end{aligned} \]
\[\begin{aligned} K_3=&F(x_n+\frac{h}2,Y_n+\frac{h}2K_2) \\ =&\begin{pmatrix} f(x_n+\frac{h}2,y_n+\frac{h}2K_2^{(1)},z_n+\frac{h}2K_2^{(2)} \\ g(x_n+\frac{h}2,y_n+\frac{h}2K_2^{(1)},z_n+\frac{h}2K_2^{(2)} \end{pmatrix} \end{aligned} \]
\[\begin{aligned} K_4=&F(x_n+{h},Y_n+{h}K_3) \\ =&\begin{pmatrix} f(x_n+{h},y_n+{h}K_3^{(1)},z_n+{h}K_3^{(2)} \\ g(x_n+{h},y_n+{h}K_3^{(1)},z_n+{h}K_3^{(2)} \end{pmatrix} \end{aligned} \]

例 2. 解方程组

\[\begin{cases} \frac{du}{dt}=0.05u(1-\frac{u}{20})-0.002uv \\ \frac{dv}{dt}=0.09v(1-\frac{v}{15})-0.15uv \\ u(0)=0.193 \\ v(0)=0.083 \end{cases} \]

解. 记

\[F(t,u,v)=\begin{pmatrix} 0.05u(1-\frac{u}{20})-0.002uv \\0.09v(1-\frac{v}{15})-0.15uv \end{pmatrix} \]

得到结果 

=========== 4th order Runge-Kutta 
value: [[0.193      0.083     ]
 [0.2027603  0.08811575]
 [0.21300665 0.09340366]
 [0.22376252 0.09884991]
 [0.23505244 0.10443747]
 [0.24690206 0.11014585]
 [0.2593382  0.11595088]
 [0.27238886 0.12182453]
 [0.28608334 0.12773475]
 [0.30045223 0.13364542]
 [0.31552748 0.13951627]]
================= Euler Forward
value: [[0.193      0.083     ]
 [0.20252484 0.08802582]
 [0.21251289 0.09322754]
 [0.222986   0.09859406]
 [0.23396702 0.10411143]
 [0.24547981 0.10976263]
 [0.25754926 0.11552731]
 [0.27020138 0.12138159]
 [0.28346334 0.12729791]
 [0.29736346 0.13324485]
 [0.31193132 0.13918704]]

高阶微分方程的初值问题

对于高阶的初值问题

\[\begin{cases} y^{(m)}(x)=f(x,y,y'(x),\cdots,y^{(m-1)}(x)) \\ y(a)=s_1,y'(a)=s_2,\cdots,y^{(m-1)}(a)=s_m \end{cases} \]

引入变量$y_1=y$, $y_2=y'$, $\cdots$, $y_m=y^{(m-1)}$后，

得到

\[\begin{cases} y'_1=y_2 \\ y'_2=y_3 \\ \cdots \\ y'_{m-1}=y_m \\ y'_m(x)=f(x,y_1,y_2,\cdots,y_m) \\ y_1(a)=s_1,y_2(a)=s_2,\cdots,y_m(a)=s_m \end{cases} \]

这是一个关于$(y_1,y_2, \cdots, y_m)$的微分方程组。

例 3. 求解如下的二阶微分方程的初值问题

\[\begin{cases} y''(x)-y'(x)=x, x\in[0,1] \\ y(0)=0 \\ y'(0)=1 \end{cases} \]

解. 作变换$z=y'(x)$，则有

\[\begin{cases} y'(x)=z \\ z'(x)=z+x \\ y(0)=0 \\ z(0)=1 \end{cases} \]

用经典的4阶Runge-Kutta方法求解。

取$h=0.1$, $x_0=0$, $z_0=1$, $y_0=0$。

Runge-Kutta公式中，记$K=(l,m)^T$，则有

\[\begin{cases} l_1=z_0=1 \\ m_1=z_0+x_0=1+0=1 \end{cases} \]
\[\begin{cases} l_2=z_0+\frac{h}2m_1=1.05 \\ m_2=(z_0+\frac{h}2m_1)+(x_0+\frac{h}2)=1.1 \end{cases} \]
\[\begin{cases} l_3=z_0+\frac{h}2m_2=1.055 \\ m_3=(z_0+\frac{h}2m_2)+(x_0+\frac{h}2)=1.105 \end{cases} \]
\[\begin{cases} l_4=z_0+{h}m_3=1.105 \\ m_4=(z_0+{h}m_3)+(x_0+{h})=1.2105 \end{cases} \]

则

\[\begin{cases} y_1=y_0+\frac{h}6(l_1+2l_2+2l_3+l_4)=0.1053 \\ z_2=z_0+\frac{h}6(m_1+2m_2+2m_3+m_4)=1.1104 \end{cases} \]

注. 方程的真解 $y(x)=x-\frac{x^2}2$，则$y(0.1)=0.1050$

目录

本节读完

例 4.

[#ex9-1-0].