向量形式四阶龙格库塔法的仿真细节

向量形式四阶龙格库塔法的仿真细节"/>

向量形式四阶龙格库塔法的仿真细节

摘要 给出了四阶龙格库塔法（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)+2h​K1​)f(tn​+2h​,x(n)+2h​K2​)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} K 1​=K 2​=K 3​=K 4​=x n+1​=​f(tn​,x (n))f(tn​+2h​,x (n)+2h​K 1​)f(tn​+2h​,x (n)+2h​K 2​)f(tn​+h,x (n)+hK 3​)x (n)+6h​(K 1​+2K 2​+2K 3​+K 4​)​(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​)=xb​x˙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)+2h​K1b​xb​(n)+2h​K2b​xb​(n)+hK3b​​K1b​=K2b​=K3b​=K4b​=​f(tn​,xa​(n),xb​(n))f(tn​+2h​,xa​(n)+2h​K1a​,xb​(n)+2h​K1b​)f(tn​+2h​,xa​(n)+2h​K2a​,xb​(n)+2h​K2b​)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)+2h​xb​(n),xb​(n)+2h​K1b​)f(tn​+2h​,xa​(n)+2h​xb​(n)+4h2​,xb​(n)+2h​K2b​)f(tn​+h,xa​(n)+hxb​(n)+2h2​K2b​,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(λ1​K1​+λ2​K2​)​
只需要满足 λ 1 + λ 2 = 1 \lambda_1+\lambda_2=1 λ1​+λ2​=1 和 λ 2 p = 1 \lambda_2p=1 λ2​p=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)+2h​K1​)
而 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 + 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) x n+1​=x (n)+6h​(K 1​+2K 2​+2K 3​+K 4​)