文章介绍了常微分方程在工程中的应用以及解析解的局限性，重点讲述了数值求解方法，包括Euler方法和Runge-Kutta方法。Euler方法分为向前和向后两种，Runge-Kutta方法则是一种更精确的单步多级方法，通过配置不同的系数可以实现不同精度的求解。
常微分方程数值方法
常微分方程
微分方程是由自变量、未知函数和未知函数的导数一起建立的方程。常微分方程是微分方程中未知函数是自变量一元函数的类型。常微分方程在电路设计，化学反应，自动控制等工程中有着广泛应用。只有很少一部分常微分方程能够找到解析解，在绝大数实际应用中，数值求解占据着主导地位。
本章将主要介绍如何离散求解常微分方程的初值问题，
{
d
y
d
x
=
f
(
x
,
y
)
,
a
<
x
<
b
y
(
a
)
=
y
0
\begin{cases} \frac{d y }{dx}=f(x,y),a<x<b \\ y(a)=y_0 \end{cases}
{dxdy=f(x,y),a<x<by(a)=y0
初值问题解的存在和唯一性，请参考文献1中的证明。其中，
y
y
y可以是向量也可以是标量。
符号说明
常微分方程的解析解
y
(
x
)
y(x)
y(x)是区间
[
a
,
b
]
[a,b]
[a,b]上关于变量
x
x
x的连续函数，数值解时
y
y
y在区间
[
a
,
b
]
[a,b]
[a,b]的有限个离散点。通常来说对区间
[
a
,
b
]
[a,b]
[a,b]上等距剖分成
N
N
N个区间，即有
x
n
=
a
+
n
h
,
n
=
0
,
1
,
…
,
N
x_n=a+nh,n=0,1,\ldots,N
xn=a+nh,n=0,1,…,N
其中，
h
=
(
b
−
a
)
/
N
h=(b-a)/N
h=(b−a)/N.而数值解
y
n
≈
y
(
x
n
)
y_n\approx y(x_n)
yn≈y(xn).
常微分方程的数值格式可以通过数值微商对导数的近似得到，或者通过数值积分公式得到。
单步方法
Euler 方法
假设
y
(
x
)
y(x)
y(x)是常微分方程的唯一解，在
x
n
x_n
xn处泰勒展开可以得到，
y
(
x
n
+
1
)
=
y
n
+
h
y
′
(
x
n
)
+
1
2
y
′
′
(
ξ
n
)
(1)
y(x_{n+1})=y_n+hy'(x_n)+\frac{1}{2}y''(\xi_n) \tag{1}
y(xn+1)=yn+hy′(xn)+21y′′(ξn)(1)
其中，
ξ
n
∈
[
x
n
,
x
n
+
1
]
\xi_n \in [x_n,x_{n+1}]
ξn∈[xn,xn+1].忽略二阶项,代入
y
′
(
x
n
)
=
f
(
x
n
,
y
n
)
y'(x_n)=f(x_n,y_n)
y′(xn)=f(xn,yn)得到向前Euler格式，
y
n
+
1
=
y
n
+
h
f
(
x
n
,
y
n
)
y_{n+1}=y_n+hf(x_n,y_n)
yn+1=yn+hf(xn,yn)
该求解格式计算简单，已知
x
n
,
y
n
,
x
n
+
1
x_n,y_n,x_{n+1}
xn,yn,xn+1，可以递推求解
y
n
+
1
y_{n+1}
yn+1。
类似地，在
x
n
+
1
x_{n+1}
xn+1处泰勒展开得到，
y
(
x
n
)
=
y
n
+
1
−
h
y
′
(
x
n
+
1
)
+
1
2
y
′
′
(
ξ
n
)
y(x_{n})=y_{n+1}-hy'(x_{n+1})+\frac{1}{2}y''(\xi_n)
y(xn)=yn+1−hy′(xn+1)+21y′′(ξn)
其中，
ξ
n
∈
[
x
n
,
x
n
+
1
]
\xi_n \in [x_n,x_{n+1}]
ξn∈[xn,xn+1]，但是与公式
(
1
)
(1)
(1)中
ξ
\xi
ξ不是一个变量。忽略二阶项,代入
y
′
(
x
n
+
1
)
=
f
(
x
n
+
1
,
y
n
+
1
)
y'(x_{n+1})=f(x_{n+1},y_{n+1})
y′(xn+1)=f(xn+1,yn+1)得到后Euler格式，
y
n
+
1
=
y
n
+
h
f
(
x
n
+
1
,
y
n
+
1
)
(2)
y_{n+1}=y_n+hf(x_{n+1},y_{n+1}) \tag{2}
yn+1=yn+hf(xn+1,yn+1)(2)
值得注意的是，当
(
2
)
(2)
(2)中，
f
(
x
n
+
1
,
y
n
+
1
)
f(x_{n+1},y_{n+1})
f(xn+1,yn+1)比较复杂时,每求解一个时间步都需要求解线性代数方程组或非线性方程。
显式Runge-Kutta龙格-库塔方法
龙格-库塔方法(下面简称RK)是单步多级的数值方法。可以通过配置系数，达到单步高精度的结果。
设m是一个正整数，
a
i
,
b
j
a_i,b_j
ai,bj，
c
i
c_i
ci是一些权重因子，则下述方法是m级Runge-Kutta方法，
y
n
+
1
=
y
n
+
h
(
c
1
K
1
+
…
+
c
m
K
m
)
K
1
=
f
(
x
n
,
y
n
)
K
i
=
f
(
x
n
+
a
i
h
,
y
n
+
h
∑
j
b
i
j
K
j
)
,
i
=
2
,
…
,
m
.
(3)
y_{n+1}=y_{n}+h(c_1K_1+\ldots+c_mK_m) \\ K_1 =f(x_n,y_n) \\ K_i=f(x_n+a_ih,y_n+h\sum_j b_{ij}K_j), i=2,\ldots,m. \tag{3}
yn+1=yn+h(c1K1+…+cmKm)K1=f(xn,yn)Ki=f(xn+aih,yn+hj∑bijKj),i=2,…,m.(3)
- 如果 j j j的取值范围是 0 ≤ j ≤ i 0 \leq j\leq i 0≤j≤i,则该格式是显格式。
由显格式的相容性，系数满足，
∑
j
=
1
m
c
j
=
1
,
a
i
=
∑
j
=
1
i
−
1
b
i
j
,
i
=
2
,
3
,
…
,
m
\sum_{j=1}^m c_j=1, a_i =\sum_{j=1}^{i-1} b_{ij},i=2,3,\ldots,m
j=1∑mcj=1,ai=j=1∑i−1bij,i=2,3,…,m
参数 a , b , c a,b,c a,b,c的获取方法
从Taylor展开的角度
f
f
f是关于x,y的二元函数，以两步显RK格式为例. 对于精确解
y
y
y,在
x
x
x处的Taylor展开，得到
y
(
x
+
h
)
=
y
(
x
)
+
h
y
′
(
x
)
+
h
2
2
y
′
′
(
x
)
+
o
(
h
3
)
=
y
(
x
)
+
h
(
f
(
x
,
y
)
+
h
2
(
f
x
(
x
,
y
)
+
f
y
(
x
,
y
)
f
(
x
,
y
)
+
o
(
h
2
)
)
(4)
y(x+h)=y(x)+hy'(x)+\frac{h^2}{2}y''(x)+o(h^3)\\ =y(x)+h(f(x,y)+\frac{h}{2}(f_x(x,y)+f_y(x,y)f(x,y)+o(h^2)) \tag{4}
y(x+h)=y(x)+hy′(x)+2h2y′′(x)+o(h3)=y(x)+h(f(x,y)+2h(fx(x,y)+fy(x,y)f(x,y)+o(h2))(4)
其中，使用了二元函数的求导技巧，
d
f
(
x
,
y
(
x
)
)
d
x
=
f
x
(
x
,
y
)
+
f
y
(
x
,
y
)
f
(
x
,
y
)
\frac{df(x,y(x))}{dx}=f_x(x,y)+f_y(x,y)f(x,y)
dxdf(x,y(x))=fx(x,y)+fy(x,y)f(x,y)
对
K
2
=
f
(
x
+
a
2
h
,
y
+
b
21
h
f
(
x
,
y
)
)
K_2=f(x+a_2h,y+b_{21}hf(x,y))
K2=f(x+a2h,y+b21hf(x,y))进行Taylor展开，得到
K
2
=
f
(
x
,
y
)
+
a
2
h
f
x
(
x
,
y
)
+
b
21
h
f
(
x
,
y
)
f
y
(
x
,
y
)
+
o
(
h
2
)
(5)
K_2 = f(x,y)+a_2hf_x(x,y)+b_{21}hf(x,y)f_y(x,y)+o(h^2) \tag{5}
K2=f(x,y)+a2hfx(x,y)+b21hf(x,y)fy(x,y)+o(h2)(5)
将
(
5
)
(5)
(5)带回RK格式(3)中，得到，
y
n
+
1
=
y
n
+
h
(
c
1
f
(
x
n
,
y
n
)
+
c
2
(
f
(
x
,
y
)
+
a
2
h
f
x
(
x
,
y
)
+
b
21
h
f
(
x
,
y
)
f
y
(
x
,
y
)
)
+
o
(
h
2
)
)
(6)
y_{n+1}=y_{n}+h(c_1f(x_n,y_n)+c_2(f(x,y)+a_2hf_x(x,y)+b_{21}hf(x,y)f_y(x,y))+o(h^2)) \tag{6}
yn+1=yn+h(c1f(xn,yn)+c2(f(x,y)+a2hfx(x,y)+b21hf(x,y)fy(x,y))+o(h2))(6)
比较格式
(
4
)
(
6
)
(4)(6)
(4)(6)可以得到系数
a
,
b
,
c
a,b,c
a,b,c的方程，
{
c
1
+
c
2
=
1
,
c
2
a
2
=
1
/
2
,
c
2
b
21
=
1
/
2
,
\begin{cases} c_1+c_2=1, \\ c_2a_2=1/2, \\ c_2b_{21}=1/2, \end{cases}
⎩
⎨
⎧c1+c2=1,c2a2=1/2,c2b21=1/2,
该方程有无穷多解。
- 取 c 1 = c 2 = 1 / 2 c_1=c_2=1/2 c1=c2=1/2, a 2 = b 21 = 1 a_2=b_{21}=1 a2=b21=1可以得到改进的Euler公式，
{ y n + 1 = y n + h 2 ( K 1 + K 2 ) K 1 = f ( x n , y n ) , K 2 = f ( x n + 1 , y n + h K 1 ) \begin{cases} y_{n+1}=y_n+\frac{h}{2}(K_1+K_2) \\ K_1=f(x_n,y_n), \\ K_2=f(x_{n+1},y_{n}+hK_1) \end{cases} ⎩ ⎨ ⎧yn+1=yn+2h(K1+K2)K1=f(xn,yn),K2=f(xn+1,yn+hK1)
- 取 c 1 = 1 / 4 , c 2 = 3 / 4 , a 2 = b 21 = 2 / 3 c_1=1/4,c_2=3/4,a_2=b_{21}=2/3 c1=1/4,c2=3/4,a2=b21=2/3,
{ y n + 1 = y n + h ( 1 4 K 1 + 3 4 K 2 ) K 1 = f ( x n , y n ) , K 2 = f ( x n + 2 h 3 , y n + 2 3 h f ( x n , y n ) ) \begin{cases} y_{n+1}=y_n+h(\frac{1}{4}K_1+\frac{3}{4}K_2) \\ K_1=f(x_n,y_n), \\ K_2=f(x_{n}+\frac{2h}{3},y_{n}+\frac{2}{3}hf(x_n,y_n)) \end{cases} ⎩ ⎨ ⎧yn+1=yn+h(41K1+43K2)K1=f(xn,yn),K2=f(xn+32h,yn+32hf(xn,yn))
从数值积分的角度
首先通过数值积分公式确定 c c c的系数，然后通过f(x,y)的展开确定 a , b a,b a,b的系数。
还是以两步显式RK格式为例，对方程
(
1
)
(1)
(1)在区间
[
x
n
,
x
n
+
1
]
[x_n,x_{n+1}]
[xn,xn+1]上做积分得到，
y
n
+
1
=
y
n
+
∫
x
n
x
n
+
1
f
(
x
,
y
)
d
x
y_{n+1}=y_n+\int_{x_n}^{x_{n+1}} f(x,y)dx
yn+1=yn+∫xnxn+1f(x,y)dx
使用梯形积分公式对上式中积分项做数值积分得到，此时已经确定
c
1
=
c
2
=
1
/
2
c_1=c_2=1/2
c1=c2=1/2,
y
n
+
1
=
y
n
+
h
2
(
f
(
x
n
,
y
n
)
+
f
(
x
n
+
1
,
y
n
+
1
)
）
+
o
(
h
2
)
y_{n+1}=y_n+\frac{h}{2}(f(x_n,y_n)+f(x_{n+1},y_{n+1})）+o(h^2)
yn+1=yn+2h(f(xn,yn)+f(xn+1,yn+1)）+o(h2)
上式中的主项
f
(
x
n
+
1
,
y
n
+
1
)
f(x_{n+1},y_{n+1})
f(xn+1,yn+1)对
y
y
y做展开得到，
K
2
=
f
(
x
,
y
)
+
a
2
h
f
x
(
x
,
y
)
+
b
21
h
f
(
x
,
y
)
f
y
(
x
,
y
)
+
o
(
h
2
)
,
f
(
x
n
+
1
,
y
n
+
1
)
=
f
(
x
,
y
)
+
h
f
x
(
x
,
y
)
+
h
f
(
x
,
y
)
f
y
(
x
,
y
)
+
o
(
h
2
)
K_2=f(x,y)+a_2hf_x(x,y)+b_{21}hf(x,y)f_y(x,y)+o(h^2), \\ f(x_{n+1},y_{n+1})=f(x,y)+hf_x(x,y)+hf(x,y)f_y(x,y)+o(h^2)
K2=f(x,y)+a2hfx(x,y)+b21hf(x,y)fy(x,y)+o(h2),f(xn+1,yn+1)=f(x,y)+hfx(x,y)+hf(x,y)fy(x,y)+o(h2)
比较上面两个式子得到
a
1
=
b
21
=
1
a_1=b_{21}=1
a1=b21=1.
扫一扫
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
