摘要 给出了四阶龙格库塔法(ODE4)的向量形式,推导了二阶积分器串联型系统的ODE4更新公式,解释了在使用ODE4仿真高阶系统和带外部输入系统时的各种注意事项,最后给出四阶龙格库塔法只能使用一次的重要结论。
标量和向量形式
对微分方程初值问题
x
˙
=
f
(
t
,
x
)
\dot x=f(t,x)
x˙=f(t,x)
常用的四阶龙格库塔法(下文简称ODE4)更新公式为
K
1
=
f
(
t
n
,
x
(
n
)
)
K
2
=
f
(
t
n
+
h
2
,
x
(
n
)
+
h
2
K
1
)
K
3
=
f
(
t
n
+
h
2
,
x
(
n
)
+
h
2
K
2
)
K
4
=
f
(
t
n
+
h
,
x
(
n
)
+
h
K
3
)
x
(
n
+
1
)
=
x
(
n
)
+
h
6
(
K
1
+
2
K
2
+
2
K
3
+
K
4
)
\begin{aligned} K_1 =& f(t_n,x(n)) \\ K_2 =& f(t_n+\frac h2,x(n)+\frac h2K_1) \\ K_3 =& f(t_n+\frac h2,x(n)+\frac h2K_2) \\ K_4 =& f(t_n+h,x(n)+hK_3) \\ x(n+1) =& x(n)+\frac h6(K_1+2K_2+2K_3+K_4) \end{aligned}
K1=K2=K3=K4=x(n+1)=f(tn,x(n))f(tn+2h,x(n)+2hK1)f(tn+2h,x(n)+2hK2)f(tn+h,x(n)+hK3)x(n)+6h(K1+2K2+2K3+K4)
向量形式为
x
⃗
˙
=
f
(
t
,
x
⃗
)
\dot{\vec x}=f(t,\vec x)
x˙=f(t,x)
K
⃗
1
=
f
(
t
n
,
x
⃗
(
n
)
)
K
⃗
2
=
f
(
t
n
+
h
2
,
x
⃗
(
n
)
+
h
2
K
⃗
1
)
K
⃗
3
=
f
(
t
n
+
h
2
,
x
⃗
(
n
)
+
h
2
K
⃗
2
)
K
⃗
4
=
f
(
t
n
+
h
,
x
⃗
(
n
)
+
h
K
⃗
3
)
x
⃗
n
+
1
=
x
⃗
(
n
)
+
h
6
(
K
⃗
1
+
2
K
⃗
2
+
2
K
⃗
3
+
K
⃗
4
)
(1)
\begin{aligned} \vec K_1 =& f(t_n,\vec x(n)) \\ \vec K_2 =& f(t_n+\frac h2,\vec x(n)+\frac h2\vec K_1) \\ \vec K_3 =& f(t_n+\frac h2,\vec x(n)+\frac h2\vec K_2) \\ \vec K_4 =& f(t_n+h,\vec x(n)+h\vec K_3) \\ \vec x_{n+1} =& \vec x(n)+\frac h6(\vec K_1+2\vec K_2+2\vec K_3+\vec K_4) \end{aligned} \tag{1}
K1=K2=K3=K4=xn+1=f(tn,x(n))f(tn+2h,x(n)+2hK1)f(tn+2h,x(n)+2hK2)f(tn+h,x(n)+hK3)x(n)+6h(K1+2K2+2K3+K4)(1)
例如对二阶积分器串联型系统
x
¨
=
f
(
t
,
x
,
x
˙
)
\ddot x=f(t,x,\dot x)
x¨=f(t,x,x˙)
可以写成下面的形式。为了便于区分,把
x
1
,
x
2
x_1,x_2
x1,x2 对应改成
x
a
,
x
b
x_a,x_b
xa,xb。
x
˙
a
=
f
a
(
t
,
x
a
,
x
b
)
=
x
b
x
˙
b
=
f
b
(
t
,
x
a
,
x
b
)
=
f
(
t
,
x
a
,
x
b
)
\begin{aligned} & \dot x_a=f_a(t,x_a,x_b)=x_b \\ & \dot x_b=f_b(t,x_a,x_b)=f(t,x_a,x_b) \\ \end{aligned}
x˙a=fa(t,xa,xb)=xbx˙b=fb(t,xa,xb)=f(t,xa,xb)
此时更新公式为
K
1
a
=
x
b
(
n
)
K
1
b
=
f
(
t
n
,
x
a
(
n
)
,
x
b
(
n
)
)
K
2
a
=
x
b
(
n
)
+
h
2
K
1
b
K
2
b
=
f
(
t
n
+
h
2
,
x
a
(
n
)
+
h
2
K
1
a
,
x
b
(
n
)
+
h
2
K
1
b
)
K
3
a
=
x
b
(
n
)
+
h
2
K
2
b
K
3
b
=
f
(
t
n
+
h
2
,
x
a
(
n
)
+
h
2
K
2
a
,
x
b
(
n
)
+
h
2
K
2
b
)
K
4
a
=
x
b
(
n
)
+
h
K
3
b
K
4
b
=
f
(
t
n
+
h
,
x
a
(
n
)
+
h
K
3
a
,
x
b
(
n
)
+
h
K
3
b
)
(2)
\begin{aligned} K_{1a} =& x_b(n)& \quad K_{1b} =& f(t_n,x_a(n),x_b(n)) \\ K_{2a} =& x_b(n)+\frac h2 K_{1b}& K_{2b} =& f(t_n+\frac h2,x_a(n)+\frac h2 K_{1a},x_b(n)+\frac h2 K_{1b}) \\ K_{3a} =& x_b(n)+\frac h2 K_{2b}& K_{3b} =& f(t_n+\frac h2,x_a(n)+\frac h2 K_{2a},x_b(n)+\frac h2 K_{2b}) \\ K_{4a} =& x_b(n)+hK_{3b}& K_{4b} =& f(t_n+h,x_a(n)+hK_{3a},x_b(n)+hK_{3b}) \\ \end{aligned} \tag{2}
K1a=K2a=K3a=K4a=xb(n)xb(n)+2hK1bxb(n)+2hK2bxb(n)+hK3bK1b=K2b=K3b=K4b=f(tn,xa(n),xb(n))f(tn+2h,xa(n)+2hK1a,xb(n)+2hK1b)f(tn+2h,xa(n)+2hK2a,xb(n)+2hK2b)f(tn+h,xa(n)+hK3a,xb(n)+hK3b)(2)
把
K
1
a
K_{1a}
K1a 代入
K
1
b
K_{1b}
K1b 可得
K
1
b
=
f
(
t
n
,
x
a
(
n
)
,
x
b
(
n
)
)
K
2
b
=
f
(
t
n
+
h
2
,
x
a
(
n
)
+
h
2
x
b
(
n
)
,
x
b
(
n
)
+
h
2
K
1
b
)
K
3
b
=
f
(
t
n
+
h
2
,
x
a
(
n
)
+
h
2
x
b
(
n
)
+
h
2
4
,
x
b
(
n
)
+
h
2
K
2
b
)
K
4
b
=
f
(
t
n
+
h
,
x
a
(
n
)
+
h
x
b
(
n
)
+
h
2
2
K
2
b
,
x
b
(
n
)
+
h
K
3
b
)
x
b
(
n
+
1
)
=
x
b
(
n
)
+
h
6
(
K
1
b
+
2
K
2
b
+
2
K
3
b
+
K
4
b
)
x
a
(
n
+
1
)
=
x
a
(
n
)
+
h
6
(
K
1
a
+
2
K
2
a
+
2
K
3
a
+
K
4
a
)
=
x
a
(
n
)
+
h
x
b
(
n
)
+
h
2
6
(
K
1
b
+
K
2
b
+
K
3
b
)
\begin{aligned} K_{1b} =& f(t_n,x_a(n),x_b(n)) \\ K_{2b} =& f(t_n+\frac h2,x_a(n)+\frac h2 x_b(n),x_b(n)+\frac h2 K_{1b}) \\ K_{3b} =& f(t_n+\frac h2,x_a(n)+\frac h2 x_b(n)+\frac{h^2}4, x_b(n)+\frac h2 K_{2b}) \\ K_{4b} =& f(t_n+h,x_a(n)+hx_b(n)+\frac{h^2}2K_{2b},x_b(n)+hK_{3b}) \\ x_b(n+1) =& x_b(n) + \frac h6(K_{1b}+2K_{2b}+2K_{3b}+K_{4b}) \\ x_a(n+1) =& x_a(n) + \frac h6(K_{1a}+2K_{2a}+2K_{3a}+K_{4a}) \\ =& x_a(n) + hx_b(n) + \frac {h^2}6(K_{1b}+K_{2b}+K_{3b}) \end{aligned}
K1b=K2b=K3b=K4b=xb(n+1)=xa(n+1)==f(tn,xa(n),xb(n))f(tn+2h,xa(n)+2hxb(n),xb(n)+2hK1b)f(tn+2h,xa(n)+2hxb(n)+4h2,xb(n)+2hK2b)f(tn+h,xa(n)+hxb(n)+2h2K2b,xb(n)+hK3b)xb(n)+6h(K1b+2K2b+2K3b+K4b)xa(n)+6h(K1a+2K2a+2K3a+K4a)xa(n)+hxb(n)+6h2(K1b+K2b+K3b)
需要注意的是,四阶龙格库塔法不止一种形式,
h
h
h 前的系数不固定。例如二阶龙格库塔法的通用格式为
K
1
=
f
(
t
n
,
x
(
n
)
)
K
2
=
f
(
t
n
+
p
,
x
(
n
)
+
p
h
K
1
)
x
(
n
+
1
)
=
x
(
n
)
+
h
(
λ
1
K
1
+
λ
2
K
2
)
\begin{aligned} K_1 =& f(t_n,x(n)) \\ K_2 =& f(t_{n+p},x(n)+phK_1) \\ x(n+1) =& x(n)+h(\lambda_1K_1+\lambda_2K_2) \end{aligned}
K1=K2=x(n+1)=f(tn,x(n))f(tn+p,x(n)+phK1)x(n)+h(λ1K1+λ2K2)
只需要满足
λ
1
+
λ
2
=
1
\lambda_1+\lambda_2=1
λ1+λ2=1 和
λ
2
p
=
1
\lambda_2p=1
λ2p=1 的约束即可,类似的四阶龙格库塔法中为了计算方便而让
h
h
h 前的系数取了
1
2
\frac 12
21。
含有外部输入时的细节扩展
如果一个用微分方程描述的系统含有外部输入信号,即
x
˙
=
f
(
t
,
x
)
+
u
(
t
)
\dot x=f(t,x)+u(t)
x˙=f(t,x)+u(t)
或向量形式
x
⃗
˙
=
f
(
t
,
x
⃗
)
+
u
⃗
(
t
)
\dot{\vec x}=f(t,\vec x)+\vec u(t)
x˙=f(t,x)+u(t)
可以看作是
x
⃗
˙
=
g
(
t
,
x
⃗
)
=
f
(
t
,
x
⃗
)
+
u
⃗
(
t
)
\dot{\vec x}=g(t,\vec x)=f(t,\vec x)+\vec u(t)
x˙=g(t,x)=f(t,x)+u(t)
这样就可以使用同样的方法。但另一方面,如果外部输入是另一个积分器的输出,例如
x
˙
1
=
f
1
(
t
,
x
1
,
x
2
)
x
˙
2
=
u
(
t
)
\begin{aligned} \dot x_1 =& f_1(t,x_1,x_2) \\ \dot x_2 =& u(t) \\ \end{aligned}
x˙1=x˙2=f1(t,x1,x2)u(t)
在计算
x
1
x_1
x1 时,可以使用式(1)给出的ODE4的向量形式,也可以像式(2)一样分开计算
K
i
a
K_{ia}
Kia 和
K
i
b
K_{ib}
Kib,但就是不能把
x
2
x_2
x2 看成外部输入,不能把
x
˙
1
=
f
1
(
t
,
x
1
,
x
2
)
\dot x_1 = f_1(t,x_1,x_2)
x˙1=f1(t,x1,x2) 看作
x
˙
1
=
f
1
(
t
,
x
1
,
u
2
(
t
)
)
=
g
(
t
,
x
1
)
\dot x_1 = f_1(t,x_1,u_2(t))=g(t,x_1)
x˙1=f1(t,x1,u2(t))=g(t,x1),因为此时计算
K
2
,
K
3
K_2,K_3
K2,K3 时就变成了
K
2
=
f
(
t
n
+
h
2
,
x
2
(
t
n
+
h
2
)
,
x
1
(
n
)
+
h
2
K
1
)
K_2=f\left(t_n+\frac h2,x_2(t_n+\frac h2),x_1(n)+\frac h2K_1\right)
K2=f(tn+2h,x2(tn+2h),x1(n)+2hK1)
而
x
2
(
t
n
+
h
2
)
x_2(t_n+\frac h2)
x2(tn+2h) 的值无法计算。
另外,像式(2)一样分开计算时需要注意,计算顺序必须是
K
1
a
/
K
1
b
→
K
2
a
/
K
2
b
→
K
3
a
/
K
3
b
→
K
4
a
/
K
4
b
K_{1a}/K_{1b}\rightarrow K_{2a}/K_{2b} \rightarrow K_{3a}/K_{3b}\rightarrow K_{4a}/K_{4b}
K1a/K1b→K2a/K2b→K3a/K3b→K4a/K4b
而不能是先算
K
1
a
K_{1a}
K1a 到
K
4
a
K_{4a}
K4a 再算
K
1
b
K_{1b}
K1b 到
K
4
b
K_{4b}
K4b。换句话说,在使用欧拉法离散化时可以使用下面的两个公式
x
1
(
n
+
1
)
=
x
1
(
n
)
+
h
x
2
(
n
)
x
2
(
n
+
1
)
=
x
2
(
n
)
+
h
u
(
n
)
\begin{aligned} & x_1(n+1)=x_1(n)+hx_2(n) \\ & x_2(n+1)=x_2(n)+hu(n) \\ \end{aligned}
x1(n+1)=x1(n)+hx2(n)x2(n+1)=x2(n)+hu(n)
而且两个公式可以不分先后顺序。但是在使用ODE4时不能这样用,也就是下面的两个式子
x
1
(
n
+
1
)
=
x
1
(
n
)
+
h
6
(
K
1
(
x
2
)
+
⋯
+
K
4
(
x
2
)
)
x
2
(
n
+
1
)
=
x
2
(
n
)
+
h
6
(
K
1
(
u
)
+
⋯
+
K
4
(
u
)
)
\begin{aligned} & x_1(n+1)=x_1(n)+\frac h6(K_1(x_2)+\cdots+K_4(x_2)) \\ & x_2(n+1)=x_2(n)+\frac h6(K_1(u)+\cdots+K_4(u)) \\ \end{aligned}
x1(n+1)=x1(n)+6h(K1(x2)+⋯+K4(x2))x2(n+1)=x2(n)+6h(K1(u)+⋯+K4(u))
不能先计算一个再计算另一个。
总结一下就是,在编程实现ODE4时需要注意,如果第一阶段使用ODE4编写并仿真了一个包含两个积分器的二阶系统,后来需要再串联两个积分器变成4阶系统,必须把每个二维向量
K
i
K_i
Ki 改成4维,而不是重新写另外4个二维向量
K
i
K_i
Ki。也就是说,不管系统是多少阶的,下面的ODE4更新公式只能有一个,其中的
x
⃗
(
n
)
\vec x(n)
x(n) 包含了系统中所有的需要积分的状态。
x
⃗
n
+
1
=
x
⃗
(
n
)
+
h
6
(
K
⃗
1
+
2
K
⃗
2
+
2
K
⃗
3
+
K
⃗
4
)
\vec x_{n+1} = \vec x(n)+\frac h6(\vec K_1+2\vec K_2+2\vec K_3+\vec K_4)
xn+1=x(n)+6h(K1+2K2+2K3+K4)