超级会员免费看
本文介绍了龙格-库塔方法，包括显式和隐式形式，以及它们的Butcher阵列表示。显式方法计算简单但受CFL条件限制，隐式方法虽计算复杂但稳定性好，而半隐式方法则是两者的折衷。高阶精度的Runge-Kutta方法如四阶经典公式和六阶五阶精度的Nyström、Lawson公式被提及。对于不同的常微分方程，选择合适的方法至关重要。
1.1.4 龙格一库塔法
将向前欧拉法写成式 (1-37) 的形式, 可以看出它实际上利用了
f
(
x
,
u
)
f(x, u)
f(x,u) 在
x
n
x_n
xn 处的取值来计算
u
n
+
1
u_{n+1}
un+1 的值。类似地, 向后欧拉法使用了
f
(
x
,
u
)
f(x, u)
f(x,u) 在
x
n
+
1
x_{n+1}
xn+1 处的近似值来计算
u
n
+
1
u_{n+1}
un+1 的值, 而向前欧拉法和向后欧拉法的精度都是 1 阶。
{
u
n
+
1
=
u
n
+
k
1
k
1
=
h
f
(
x
n
,
u
n
)
(1-37)
\left\{\begin{array}{l} u_{n+1}=u_n+k_1 \\ k_1=h f\left(x_n, u_n\right) \end{array}\right.\tag{1-37}
{un+1=un+k1k1=hf(xn,un)(1-37)
同样, 将预测-校正法写成如下形式, 不难发现它其实是用
f
(
x
,
u
)
f(x, u)
f(x,u) 在
x
n
x_n
xn 和
x
n
+
1
x_{n+1}
xn+1 两处的取值来计算
u
n
+
1
u_{n+1}
un+1 的值, 并达到了 2 阶的精度。
{
u
n
+
1
=
u
n
+
(
k
1
+
k
2
)
/
2
k
1
=
h
f
(
x
n
,
u
n
)
k
2
=
h
f
(
x
n
+
h
,
u
n
+
k
1
)
(1-38)
\left\{\begin{array}{l} u_{n+1}=u_n+\left(k_1+k_2\right) / 2 \\ k_1=h f\left(x_n, u_n\right) \\ k_2=h f\left(x_n+h, u_n+k_1\right) \end{array}\right.\tag{1-38}
⎩
⎨
⎧un+1=un+(k1+k2)/2k1=hf(xn,un)k2=hf(xn+h,un+k1)(1-38)
上述方法都是通过
f
(
x
,
u
)
f(x, u)
f(x,u) 在不同位置的线性组合来计算
u
n
+
1
u_{n+1}
un+1 的值, 所考虑的位置越多,精度也就越高。顺着这样的思路想下去, 如果用
f
(
x
,
u
)
f(x, u)
f(x,u) 在更多位置的线性组合来构造递推公式, 将会得到更高的精度, 这就是龙格-库塔法（Runge-Kutta method）的思想。这样, 递推公式将有如下形式:
{
u
n
+
1
=
u
n
+
h
∑
i
=
1
r
R
i
k
i
k
i
=
f
(
x
n
+
a
i
h
,
u
n
+
∑
j
=
1
i
−
1
b
i
j
k
j
)
,
i
=
1
,
2
,
3
,
…
,
r
(1-39)
\left\{\begin{array}{l} u_{n+1}=u_n+h \sum\limits_{i=1}^r R_i k_i \\ k_i=f\left(x_n+a_i h, u_n+\sum\limits_{j=1}^{i-1} b_{i j} k_j\right), \quad i=1,2,3, \ldots, r \end{array}\right.\tag{1-39}
⎩
⎨
⎧un+1=un+hi=1∑rRikiki=f(xn+aih,un+j=1∑i−1bijkj),i=1,2,3,…,r(1-39)
其中,
R
i
、
a
i
、
b
i
j
R_i 、 a_i 、 b_{i j}
Ri、ai、bij 为待定常数。
以上 Runge-Kutta 公式又可以写成如下既直观又简洁的阵列形式 (Butcher表, Butcher, 1964):
a 1 b 1 , 1 ⋯ b 1 , r ⋮ ⋮ ⋱ ⋮ a r b r , 1 ⋯ b r , r R 1 ⋯ R r \begin{array}{c|ccc} a_1 & b_{1,1} & \cdots & b_{1, r} \\ \vdots & \vdots & \ddots & \vdots \\ a_r & b_{r, 1} & \cdots & b_{r, r} \\ \hline & R_1 & \cdots & R_r \end{array} a1⋮arb1,1⋮br,1R1⋯⋱⋯⋯b1,r⋮br,rRr
由于使用了
f
(
x
,
u
)
f(x, u)
f(x,u) 在
r
r
r 个位置的取值, 称递推公式 (1-39) 为
r
r
r 阶龙格-库塔法。为了确定其在某一精度下的待定系数, 需要将式中的
u
n
+
1
u_{n+1}
un+1 在
x
n
x_n
xn 处泰勒展开, 并与
u
(
x
n
+
1
)
u\left(x_{n+1}\right)
u(xn+1) 在
x
n
x_n
xn 处的泰勒展开式比较:
u
(
x
n
+
1
)
=
u
(
x
n
)
+
h
1
!
u
′
(
x
n
)
+
h
2
2
!
u
′
′
(
x
n
)
+
…
+
h
k
k
!
u
(
k
)
(
x
n
)
+
O
(
h
k
+
1
)
(1-40)
u\left(x_{n+1}\right)=u\left(x_n\right)+\frac{h}{1 !} u^{\prime}\left(x_n\right)+\frac{h^2}{2 !} u^{\prime \prime}\left(x_n\right)+\ldots+\frac{h^k}{k !} u^{(k)}\left(x_n\right)+O\left(h^{k+1}\right)\tag{1-40}
u(xn+1)=u(xn)+1!hu′(xn)+2!h2u′′(xn)+…+k!hku(k)(xn)+O(hk+1)(1-40)
若待定常数
R
i
、
a
i
、
b
i
j
R_i 、 a_i 、 b_{i j}
Ri、ai、bij 的取值使得两式最多有前
k
+
1
k+1
k+1 项相同，那么:
u
(
x
n
+
1
)
−
u
n
+
1
=
O
(
h
k
+
1
)
(1-41)
u\left(x_{n+1}\right)-u_{n+1}=O\left(h^{k+1}\right)\tag{1-41}
u(xn+1)−un+1=O(hk+1)(1-41)
这即是说, 此时
R
i
、
a
i
、
b
i
j
R_i 、 a_i 、 b_{i j}
Ri、ai、bij 的取值使得式 (1-39) 具有
k
k
k 阶的精度, 下面举例说明。 当
r
=
2
r=2
r=2 时, 龙格一库塔法递推公式为:
{
u
n
+
1
=
u
n
+
R
1
k
1
+
R
2
k
2
k
1
=
h
f
(
x
n
,
u
n
)
k
2
=
h
f
(
x
n
+
a
h
,
u
n
+
b
k
1
)
(1-42)
\left\{\begin{array}{l} u_{n+1}=u_n+R_1 k_1+R_2 k_2 \\ k_1=h f\left(x_n, u_n\right) \\ k_2=h f\left(x_n+a h, u_n+b k_1\right) \end{array}\right.\tag{1-42}
⎩
⎨
⎧un+1=un+R1k1+R2k2k1=hf(xn,un)k2=hf(xn+ah,un+bk1)(1-42)
将
k
2
k_2
k2 的表达式在
x
n
x_n
xn 处展开:
k
2
=
h
f
(
x
n
+
a
h
,
u
n
+
b
k
1
)
=
h
[
f
(
x
n
,
u
n
)
+
a
h
f
x
(
x
n
,
u
n
)
+
b
k
1
f
u
(
x
n
,
u
n
)
+
O
(
h
2
)
]
=
h
f
(
x
n
,
u
n
)
+
h
2
[
a
f
x
(
x
n
,
u
n
)
+
b
f
u
(
x
n
,
u
n
)
f
(
x
n
,
u
n
)
]
+
O
(
h
3
)
(1-43)
\begin{aligned} k_2 & =h f\left(x_n+a h, u_n+b k_1\right) \\ & =h\left[f\left(x_n, u_n\right)+a h f_x\left(x_n, u_n\right)+b k_1 f_u\left(x_n, u_n\right)+O\left(h^2\right)\right] \\ & =h f\left(x_n, u_n\right)+h^2\left[a f_x\left(x_n, u_n\right)+b f_u\left(x_n, u_n\right) f\left(x_n, u_n\right)\right]+O\left(h^3\right) \end{aligned}\tag{1-43}
k2=hf(xn+ah,un+bk1)=h[f(xn,un)+ahfx(xn,un)+bk1fu(xn,un)+O(h2)]=hf(xn,un)+h2[afx(xn,un)+bfu(xn,un)f(xn,un)]+O(h3)(1-43)
将
k
1
k_1
k1 和
k
2
k_2
k2 代入
u
n
+
1
u_{n+1}
un+1 的表达式:
u
n
+
1
=
u
n
+
h
(
R
1
+
R
2
)
f
(
x
n
,
u
n
)
+
h
2
[
a
R
2
f
x
(
x
n
,
u
n
)
+
b
R
2
f
u
(
x
n
,
u
n
)
f
(
x
n
,
u
n
)
]
+
O
(
h
3
)
(1-44)
u_{n+1}=u_n+h\left(R_1+R_2\right) f\left(x_n, u_n\right)+h^2\left[a R_2 f_x\left(x_n, u_n\right)+b R_2 f_u\left(x_n, u_n\right) f\left(x_n, u_n\right)\right]+O\left(h^3\right)\tag{1-44}
un+1=un+h(R1+R2)f(xn,un)+h2[aR2fx(xn,un)+bR2fu(xn,un)f(xn,un)]+O(h3)(1-44)
假设
u
n
=
u
(
x
n
)
u_n=u\left(x_n\right)
un=u(xn), 并根据式 (1-1) 中的第 1 个式子, 使式 (1-40) 和式（1-44）的前 3 项相同，得到:
{
R
1
+
R
2
=
1
a
R
2
=
1
/
2
b
R
2
=
1
/
2
(1-45)
\left\{\begin{array}{l} R_1+R_2=1 \\ a R_2=1 / 2 \\ b R_2=1 / 2 \end{array}\right.\tag{1-45}
⎩
⎨
⎧R1+R2=1aR2=1/2bR2=1/2(1-45)
满足上式的
R
1
、
R
2
、
a
、
b
R_1 、 R_2 、 a 、 b
R1、R2、a、b 可以有无穷多种情况, 每种情况的局部截断误差都是
O
(
h
3
)
O\left(h^3\right)
O(h3), 都被称为二阶龙格-库塔法。如果取
R
1
=
R
2
=
1
/
2
、
a
=
b
=
1
R_1=R_2=1 / 2 、 a=b=1
R1=R2=1/2、a=b=1, 就是前面提到的预测-校正法。 如果取
R
1
=
0
、
R
2
=
1
、
a
=
b
=
1
/
2
R_1=0 、 R_2=1 、 a=b=1 / 2
R1=0、R2=1、a=b=1/2, 就得到中点公式 (midpoint method):
{
u
n
+
1
=
u
n
+
k
2
k
1
=
h
f
(
x
n
,
u
n
)
k
2
=
h
f
(
x
n
+
1
2
h
,
u
n
+
1
2
k
1
)
(1-46)
\left\{\begin{array}{l} u_{n+1}=u_n+k_2 \\ k_1=h f\left(x_n, u_n\right) \\ k_2=h f\left(x_n+\frac{1}{2} h, u_n+\frac{1}{2} k_1\right) \end{array}\right.\tag{1-46}
⎩
⎨
⎧un+1=un+k2k1=hf(xn,un)k2=hf(xn+21h,un+21k1)(1-46)
二阶龙格一库塔法不能保证局部截断误差为
O
(
h
4
)
O\left(h^4\right)
O(h4), 所以接下来考虑
r
=
3
r=3
r=3 的情况, 即:
{
u
n
+
1
=
u
n
+
R
1
k
1
+
R
2
k
2
+
R
3
k
3
k
1
=
h
f
(
x
n
,
u
n
)
k
2
=
h
f
(
x
n
+
a
2
h
,
u
n
+
b
21
k
1
)
k
3
=
h
f
(
x
n
+
a
3
h
,
u
n
+
b
31
k
1
+
b
32
k
2
)
(1-47)
\left\{\begin{array}{l} u_{n+1}=u_n+R_1 k_1+R_2 k_2+R_3 k_3 \\ k_1=h f\left(x_n, u_n\right) \\ k_2=h f\left(x_n+a_2 h, u_n+b_{21} k_1\right) \\ k_3=h f\left(x_n+a_3 h, u_n+b_{31} k_1+b_{32} k_2\right) \end{array}\right.\tag{1-47}
⎩
⎨
⎧un+1=un+R1k1+R2k2+R3k3k1=hf(xn,un)k2=hf(xn+a2h,un+b21k1)k3=hf(xn+a3h,un+b31k1+b32k2)(1-47)
与
r
=
2
r=2
r=2 的情况类似, 先写出上式中的
u
n
+
1
u_{n+1}
un+1 在
x
n
x_n
xn 处的泰勒展开式, 并假设
u
n
=
u
(
x
n
)
u_n=u\left(x_n\right)
un=u(xn), 与式 (1-40) 比较, 要保证局部截断误差为
u
(
x
n
+
1
)
−
u
n
+
1
=
O
(
h
4
)
u\left(x_{n+1}\right)-u_{n+1}=O\left(h^4\right)
u(xn+1)−un+1=O(h4), 得:
{
R
1
+
R
2
+
R
3
=
1
a
2
=
b
21
a
3
=
b
31
+
b
32
a
2
R
2
+
a
3
R
3
=
1
/
2
a
2
2
R
2
+
a
3
2
R
3
=
1
/
3
a
2
b
32
R
3
=
1
/
6
(1-48)
\left\{\begin{array}{l} R_1+R_2+R_3=1 \\ a_2=b_{21} \\ a_3=b_{31}+b_{32} \\ a_2 R_2+a_3 R_3=1 / 2 \\ a_2^2 R_2+a_3^2 R_3=1 / 3 \\ a_2 b_{32} R_3=1 / 6 \end{array}\right.\tag{1-48}
⎩
⎨
⎧R1+R2+R3=1a2=b21a3=b31+b32a2R2+a3R3=1/2a22R2+a32R3=1/3a2b32R3=1/6(1-48)
方程组 (1-48) 的解有无穷多个, 每组解都确定一个三阶龙格-库塔法的公式。比较简单且重要的一组解是
R
1
=
1
/
6
、
R
2
=
4
/
6
、
R
3
=
1
/
6
、
a
2
=
1
/
2
、
a
3
=
1
、
b
21
=
1
/
2
、
b
31
=
−
1
、
b
32
=
2
R_1=1 / 6 、 R_2=4 / 6 、 R_3=1 / 6 、 a_2=1 / 2 、 a_3=1 、 b_{21}=1 / 2 、 b_{31}=-1 、 b_{32}=2
R1=1/6、R2=4/6、R3=1/6、a2=1/2、a3=1、b21=1/2、b31=−1、b32=2, 它对应的是三级三阶精度的显式库塔公式 (Kutta’s third-order method), 如下:
{
u
n
+
1
=
u
n
+
1
6
(
k
1
+
4
k
2
+
k
3
)
k
1
=
h
f
(
x
n
,
u
n
)
k
2
=
h
f
(
x
n
+
1
2
h
,
u
n
+
1
2
k
1
)
k
3
=
h
f
(
x
n
+
h
,
u
n
−
k
1
+
2
k
2
)
(1-49)
\left\{\begin{array}{l} u_{n+1}=u_n+\frac{1}{6}\left(k_1+4 k_2+k_3\right) \\ k_1=h f\left(x_n, u_n\right) \\ k_2=h f\left(x_n+\frac{1}{2} h, u_n+\frac{1}{2} k_1\right) \\ k_3=h f\left(x_n+h, u_n-k_1+2 k_2\right) \end{array}\right.\tag{1-49}
⎩
⎨
⎧un+1=un+61(k1+4k2+k3)k1=hf(xn,un)k2=hf(xn+21h,un+21k1)k3=hf(xn+h,un−k1+2k2)(1-49)
类似地,
r
=
4
r=4
r=4 时, 龙格-库塔法的递推公式的形式为:
{
u
n
+
1
=
u
n
+
R
1
k
1
+
R
2
k
2
+
R
3
k
3
+
R
4
k
4
k
1
=
h
f
(
x
n
,
u
n
)
k
2
=
h
f
(
x
n
+
a
2
h
,
u
n
+
b
21
k
1
)
k
3
=
h
f
(
x
n
+
a
3
h
,
u
n
+
b
31
k
1
+
b
32
k
2
)
k
4
=
h
f
(
x
n
+
a
4
h
,
u
n
+
b
41
k
1
+
b
42
k
2
+
b
43
k
3
)
(1-50)
\left\{\begin{array}{l} u_{n+1}=u_n+R_1 k_1+R_2 k_2+R_3 k_3+R_4 k_4 \\ k_1=h f\left(x_n, u_n\right) \\ k_2=h f\left(x_n+a_2 h, u_n+b_{21} k_1\right) \\ k_3=h f\left(x_n+a_3 h, u_n+b_{31} k_1+b_{32} k_2\right) \\ k_4=h f\left(x_n+a_4 h, u_n+b_{41} k_1+b_{42} k_2+b_{43} k_3\right) \end{array}\right.\tag{1-50}
⎩
⎨
⎧un+1=un+R1k1+R2k2+R3k3+R4k4k1=hf(xn,un)k2=hf(xn+a2h,un+b21k1)k3=hf(xn+a3h,un+b31k1+b32k2)k4=hf(xn+a4h,un+b41k1+b42k2+b43k3)(1-50)
在局部截断误差为
O
(
h
5
)
O\left(h^5\right)
O(h5) 的前提下, 得到关于上式中待定常数的方程组的过程比较复杂, 这里从略。下面给出两组常用的四阶龙格一库塔法的公式。
四级四阶精度的古典显式龙格-库塔公式 (classic fourth-order method):
{
u
n
+
1
=
u
n
+
1
6
(
k
1
+
2
k
2
+
2
k
3
+
k
4
)
k
1
=
h
f
(
x
n
,
u
n
)
k
2
=
h
f
(
x
n
+
1
2
h
,
u
n
+
1
2
k
1
)
k
3
=
h
f
(
x
n
+
1
2
h
,
u
n
+
1
2
k
2
)
k
4
=
h
f
(
x
n
+
h
,
u
n
+
k
3
)
(1-51)
\left\{\begin{array}{l} u_{n+1}=u_n+\frac{1}{6}\left(k_1+2 k_2+2 k_3+k_4\right) \\ k_1=h f\left(x_n, u_n\right) \\ k_2=h f\left(x_n+\frac{1}{2} h, u_n+\frac{1}{2} k_1\right) \\ k_3=h f\left(x_n+\frac{1}{2} h, u_n+\frac{1}{2} k_2\right) \\ k_4=h f\left(x_n+h, u_n+k_3\right) \end{array}\right.\tag{1-51}
⎩
⎨
⎧un+1=un+61(k1+2k2+2k3+k4)k1=hf(xn,un)k2=hf(xn+21h,un+21k1)k3=hf(xn+21h,un+21k2)k4=hf(xn+h,un+k3)(1-51)
四级四阶精度的显式吉尔 (Gill) 公式:
{
u
n
+
1
=
u
n
+
1
6
[
k
1
+
(
2
−
2
)
k
2
+
(
2
+
2
)
k
3
+
k
4
]
k
1
=
h
f
(
x
n
,
u
n
)
k
2
=
h
f
(
x
n
+
1
2
h
,
u
n
+
1
2
k
1
)
k
3
=
h
f
(
x
n
+
1
2
h
,
u
n
+
2
−
1
2
k
1
+
2
−
2
2
k
2
)
k
4
=
h
f
(
x
n
+
h
,
u
n
−
2
2
k
2
+
2
+
2
2
k
3
)
(1-52)
\left\{\begin{array}{l} u_{n+1}=u_n+\frac{1}{6}\left[k_1+(2-\sqrt{2}) k_2+(2+\sqrt{2}) k_3+k_4\right] \\ k_1=h f\left(x_n, u_n\right) \\ k_2=h f\left(x_n+\frac{1}{2} h, u_n+\frac{1}{2} k_1\right) \\ k_3=h f\left(x_n+\frac{1}{2} h, u_n+\frac{\sqrt{2}-1}{2} k_1+\frac{2-\sqrt{2}}{2} k_2\right) \\ k_4=h f\left(x_n+h, u_n-\frac{\sqrt{2}}{2} k_2+\frac{2+\sqrt{2}}{2} k_3\right) \end{array}\right.\tag{1-52}
⎩
⎨
⎧un+1=un+61[k1+(2−2)k2+(2+2)k3+k4]k1=hf(xn,un)k2=hf(xn+21h,un+21k1)k3=hf(xn+21h,un+22−1k1+22−2k2)k4=hf(xn+h,un−22k2+22+2k3)(1-52)
Gill公式有减少舍入误差的优点。
一般地, 方程 (1-39) 中的系数满足
∑
j
=
1
r
R
j
=
1
,
a
i
=
∑
j
=
1
i
−
1
b
i
,
j
,
i
=
2
,
3
,
⋯
r
.
\sum_{j=1}^r R_j=1, a_i=\sum_{j=1}^{i-1} b_{i, j}, i=2,3, \cdots r .
j=1∑rRj=1,ai=j=1∑i−1bi,j,i=2,3,⋯r.
待定参数
{
R
i
}
,
{
a
i
}
\left\{R_i\right\},\left\{a_i\right\}
{Ri},{ai}, 和
{
b
i
,
j
}
\left\{b_{i, j}\right\}
{bi,j} 的出现为构造高精度的数值方法创造了条件. 需要强调的是,
r
r
r 的取值并不总是和精度的阶数相同, 二者关系如下表所示。随着
r
r
r 的增加, 每递推一步的运算量也在增加, 对于多数实际问题, 四阶龙格-库塔法的精度就已经足够了, 因此它也是最常用的龙格-库塔法。当
r
=
r=
r=
2
,
3
,
4
2,3,4
2,3,4 时, 存在
r
r
r 阶精度的显式 Runge-Kutta 方法; 但当
r
=
5
r=5
r=5 时, 至多能获得四阶精度的显式 Runge-Kutta 方法; 为了构造五阶精度的显式 Runge-Kutta 方法,
r
r
r 至少应为 6; 当
r
=
7
r=7
r=7 时, 可得到六阶精度的显式 Runge-Kutta 方法; 而当
r
⩾
8
r \geqslant 8
r⩾8 时, 至多可构造出一个
(
r
−
2
)
(r-2)
(r−2) 阶精度的显式 Runge-Kutta 方法.
| r 取值 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | ⩾8 |
|---|---|---|---|---|---|---|---|---|
| 精度阶数 | 1 | 2 | 3 | 4 | 4 | 5 | 6 | r-2 |
下面是一些高阶精度的显式 Runge-Kutta 方法例子 (以 Butcher 阵列表示).
六级五阶精度的显式 Nyström 公式
0
1
3
1
3
2
5
4
25
6
25
1
1
4
−
3
15
4
2
3
6
81
90
81
−
50
81
8
81
4
5
6
75
36
75
10
75
8
75
0
23
192
0
125
192
0
−
81
192
125
192
\begin{aligned} &\begin{array}{c|cccccc} 0 & & & & & & \\ \frac{1}{3} & \frac{1}{3} & & & & & \\ \frac{2}{5} & \frac{4}{25} & \frac{6}{25} & & & & \\ 1 & \frac{1}{4} & -3 & \frac{15}{4} & & & \\ \frac{2}{3} & \frac{6}{81} & \frac{90}{81} & -\frac{50}{81} & \frac{8}{81} & & \\ \frac{4}{5} & \frac{6}{75} & \frac{36}{75} & \frac{10}{75} & \frac{8}{75} & 0 & \\ \hline & \frac{23}{192} & 0 & \frac{125}{192} & 0 & -\frac{81}{192} & \frac{125}{192} \end{array} \end{aligned}
0315213254312544181675619223256−3819075360415−8150751019212581875800−19281192125
六级五阶精度的显式 Lawson 公式
0
1
2
1
2
1
4
3
16
1
16
1
2
0
0
1
2
3
4
0
−
3
16
6
16
9
16
1
1
7
4
7
6
7
−
12
7
8
7
7
90
0
32
90
12
90
32
90
7
90
\begin{array}{c|cccccc} 0 & & & & & & \\ \frac{1}{2} & \frac{1}{2} & & & & & \\ \frac{1}{4} & \frac{3}{16} & \frac{1}{16} & & & & \\ \frac{1}{2} & 0 & 0 & \frac{1}{2} & & & \\ \frac{3}{4} & 0 & -\frac{3}{16} & \frac{6}{16} & \frac{9}{16} & & \\ 1 & \frac{1}{7} & \frac{4}{7} & \frac{6}{7} & -\frac{12}{7} & \frac{8}{7} & \\ \hline & \frac{7}{90} & 0 & \frac{32}{90} & \frac{12}{90} & \frac{32}{90} & \frac{7}{90} \end{array}
02141214312116300719071610−16374021166769032169−7129012789032907
七级六阶精度的显式 Butcher 公式
0
1
3
1
3
2
3
0
2
3
1
3
1
12
1
3
−
1
12
1
2
−
1
16
9
8
−
3
16
−
3
8
1
2
0
9
8
−
3
8
−
3
4
1
2
1
2
9
44
−
9
11
63
44
18
11
0
−
16
11
1
11
120
0
27
40
27
40
−
4
15
−
4
15
11
120
\begin{array}{c|ccccccc} 0 & & & & & & & \\ \frac{1}{3} & \frac{1}{3} & & & & & & \\ \frac{2}{3} & 0 & \frac{2}{3} & & & & \\ \frac{1}{3} & \frac{1}{12} & \frac{1}{3} & -\frac{1}{12} & & & \\ \frac{1}{2} & -\frac{1}{16} & \frac{9}{8} & -\frac{3}{16} & -\frac{3}{8} & & & \\ \frac{1}{2} & 0 & \frac{9}{8} & -\frac{3}{8} & -\frac{3}{4} & \frac{1}{2} & & \\ \frac{1}{2} & \frac{9}{44} & -\frac{9}{11} & \frac{63}{44} & \frac{18}{11} & 0 & -\frac{16}{11} & \\ \hline 1 & \frac{11}{120} & 0 & \frac{27}{40} & \frac{27}{40} & -\frac{4}{15} & -\frac{4}{15} & \frac{11}{120} \end{array}
03132312121211310121−16104491201132318989−1190−121−163−8344634027−83−4311184027210−154−1116−15412011
值得注意的是, 显式龙格库塔方法是一种数值求解常微分方程的方法，它使用固定的步长来计算下一个时间步的解。这种方法在计算上非常高效，因为它不需要解决非线性方程，但它的稳定性却受到严重的限制。
具体而言，显式龙格库塔方法的稳定性限制来自于一个称为"CFL条件"的限制，这个条件限制了时间步长和空间步长之间的关系。如果时间步长过长，那么计算结果就会出现震荡和不稳定的情况。因此，为了确保计算结果的稳定性，需要选择非常小的时间步长，这会导致计算的时间复杂度变高。
为了克服显式龙格库塔方法的这个不足，人们引入了一种称为"隐式龙格库塔方法"的方法。这种方法使用迭代的方式来解决非线性方程，因此可以使用更大的时间步长，从而提高计算效率。虽然隐式龙格库塔方法计算时的时间复杂度更高，但它却能够获得更好的稳定性和准确性，这使得它在实际应用中得到广泛的应用。
这里不再进一步介绍隐式 Runge-Kutta 的推导，而仅仅(以 Butcher 表列举三个典型的 Gauss-Legendre 方法（基于 Gauss 求积公式得到的 Runge-Kutta 方法称为 Gauss-Legendre 方法）
一级二阶精度的隐式中点公式
1
2
1
2
1
\begin{array}{l|l} \frac{1}{2} & \frac{1}{2} \\ \hline & 1 \end{array}
21211
二级四阶精度的隐式 Hammer 和 Hollingsworth 公式
3
−
3
6
1
4
3
−
2
3
12
3
+
3
6
3
+
2
3
12
1
4
1
2
1
2
\begin{array}{c|cc} \frac{3-\sqrt{3}}{6} & \frac{1}{4} & \frac{3-2 \sqrt{3}}{12} \\ \frac{3+\sqrt{3}}{6} & \frac{3+2 \sqrt{3}}{12} & \frac{1}{4} \\ \hline & \frac{1}{2} & \frac{1}{2} \end{array}
63−363+341123+2321123−234121
三级六阶精度的隐式 Kuntzmann 和 Butcher 公式
5
−
15
10
5
36
10
−
3
15
45
25
−
6
15
180
1
2
10
+
3
15
72
2
9
10
−
3
15
72
5
+
15
10
25
+
6
15
180
10
+
3
15
45
5
36
5
18
4
9
5
18
\begin{array}{c|ccc} \frac{5-\sqrt{15}}{10} & \frac{5}{36} & \frac{10-3 \sqrt{15}}{45} & \frac{25-6 \sqrt{15}}{180} \\ \frac{1}{2} & \frac{10+3 \sqrt{15}}{72} & \frac{2}{9} & \frac{10-3 \sqrt{15}}{72} \\ \frac{5+\sqrt{15}}{10} & \frac{25+6 \sqrt{15}}{180} & \frac{10+3 \sqrt{15}}{45} & \frac{5}{36} \\ \hline & \frac{5}{18} & \frac{4}{9} & \frac{5}{18} \end{array}
105−1521105+153657210+31518025+6151854510−315924510+3159418025−6157210−315365185
在计算相同个数的函数 f f f 的点值情况下, 隐式 Runge-Kutta 公式的精确度要 比显式 Runge-Kutta 公式的高, 且又有较好的数值稳定性. 这些优点对于解刚性方程组是非常有用的, 但是在隐式 Runge-Kutta 公式中, 由于每个 k i k_i ki 的 表达式中都含有 k 1 , k 2 , ⋯ , k s k_1, k_2, \cdots, k_s k1,k2,⋯,ks, 因而在使用隐式 Runge-Kutta 公式时, 往往需要解非线性方程 (组), 这会给计算带来不便和大的开销.为优化这个缺点，有学者提出半隐式的 Runge-Kutta 方法. 对单个方程, 它既保持隐式 Runge-Kutta 方法的好的稳定性, 又不需要解非线性方程 (组). 这种方法是由 Rosenbrock(1963) 通过将较高阶导数与 Runge-Kutta 方法相结合而得到的, 具体理论此处不再介绍。
显式龙格库塔方法、隐式龙格库塔方法和半隐式龙格库塔方法都是常微分方程数值解法中常用的方法，它们在不同的情况下有不同的适用性和优缺点。
- 显式龙格库塔方法
显式龙格库塔方法计算简单，计算速度快，是常用的数值方法之一。它的计算公式只需要考虑当前时间步的解，因此它的内存占用较小。但是，它的稳定性受到严格限制，步长必须小于一定的值，否则计算结果会出现震荡和不稳定的情况，这样就会导致时间复杂度增加。
- 隐式龙格库塔方法
隐式龙格库塔方法的计算公式需要同时考虑当前时间步和下一个时间步的解，需要使用迭代方法解决非线性方程，这导致计算量相对显式方法更大，时间复杂度更高。但是，由于隐式方法不受到稳定性限制，可以使用更大的步长来计算，从而减少计算时间和内存开销，并且能够处理更复杂的问题，如刚体动力学和流体动力学等。
- 半隐式龙格库塔方法
半隐式龙格库塔方法是显式和隐式方法的折中方案，它的计算公式需要同时考虑当前时间步和下一个时间步的解，其中一部分使用显式方法，另一部分使用隐式方法。半隐式方法的计算复杂度介于显式和隐式方法之间，同时具有较高的稳定性和更大的时间步长。这使得半隐式方法成为一些特殊问题的首选方法，例如多体动力学、弹性体力学和分子动力学等。
综上所述，不同的龙格库塔方法具有不同的适用范围和优缺点。在实际应用中，需要根据问题的特点和求解要求，选择适当的方法来求解常微分方程，以获得更准确和高效的结果。
扫一扫
586
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
