利用龙格-库塔法求解质点运动常微分方程
一、待解问题
讨论常微分方程的初值问题，边值问题的数值解法，最常用的基本方法就是龙格-库塔法使一些物理方程的计算简便，所得结果的精确提高。
下面是利用龙格-库塔法求解一个质点运动的常微分方程 以初速v 0=8km /s 自地球表面（半径R =6000km ）竖直向上发射一质量为 m 的火箭，如图所示。若不计空气阻力，火箭所受引力 F 大小与它到地心距离的平方成反比，求火箭所能到达的最大高度H 。
火箭为我们分析的对象，对其做受力分析，火箭在任意位置 x 处仅受地球引力F 作用。由题意知，F 的大小与 x2 成反比，设μ为比例系数，则有：
（1）
F =
μ
x 当火箭处于地面时，即x =R 时F = m g ，由式(1)可得μ=mgR 2，于
m gR 2
是火箭在任意位置 x 处所受地球引力 F 的大小为 F =2
x
由于火箭作直线运动，火箭的直线运动微分方程式为：
d 2x mgR 2m 2=-2dt x
分离变量积分
d 2x dv dv dx dv
==⋅=v 2
dt dx dx dt dt
dv mgR 2mv ⋅=-2
dx x
dv gR 2
=-2dx vx
二、物理机理
1. 龙格-库塔法的基本思想及一般形式
设初值问题y =y (x ) ∈[a , b ]，由微分中值定理可知，必存在ξ∈[x n , x n +1]，使
y (x n +1) =y (x n ) +h y '(ξ) =y (x n ) +hf (ξ, y (ξ))
设y n =y (x n ) ，并记K *=f (ξ, y (ξ)) ，则
y (x n +1) =y n +hK *
其中K *称为y (x ) 在[x n , x n +1]上的平均斜率，只要对平均斜率K *提供一种算法，上式就给出了一种数值解公式，例如，用K 1=f (x n , y n ) 代替K *，就得到欧拉公式，用K 2=f (x n +1, y n +1) 代替K *，就得到向后欧拉公式，如果用K 1, K 2的平均值来代替K *，则可得到二阶精度的梯形公式。可以设想，如果在[x n , x n +1]*库塔（Runge-Kutta ）法的基本思想。 龙格-库塔公式的一般形式：
y n +1=y n +h ∑c i k i
i =1
r
K 1=f (x n , y n )
K i =f (x n +λi h , y n +h ∑μij K j )
j =1i -1
（1）
其中K i 是y =y (x ) 在x n +λi h (0≤λi ≤1) 点的斜率预测值；c i , λi , μij 均为常数，选取这些常数的原则是使（1）式有尽可能高的精度。 2. 二阶龙格-库塔法
r =2的龙格-库塔公式为:
y n +1=y n +h (c 1K 1+c 2K 2) K 1=f (x n , y n )
K 2=f (x n +λ2h , y n +h μ21K 1)
（2）
式中c 1, c 2, λ2, μ21均为待定常数，我们希望适当的选取这些系数使公式的精度尽可能高，先将K 2按照二元函数泰勒级数：
1∂∂
f (x n +λ2h , y n +h μ21K 1) =∑(λ2h +h μ21) k f (x n , y n ) +...
∂x ∂y k =0k !
n
展开得到：
K 2=f (x n , y n ) +λ2hf x (x n , y n ) +μ21hf y (x n , y n ) f (x n , y n ) 122
[λ2h f xx (x n , y n ) +2λ2μ21h 2f xy (x n , y n ) f (x n , y n ) 2! 22+μ21h f yy (x n , y n ) f 2(x n , y n )]+... +
为了叙述方便，把f (x n , y n ) 和偏导数里面的x n , y n 省略不写，将K 2代入（2）式的第一式，得：
y n +1=y n +h (c 1+c 2) f +h 2c 2(λ2f x +μ21f y f ) c 222
+h (λ2f xx +2λ2μ21f xy f +μ21f yy f 2)
2
3
再用泰勒公式将y (x n +1) 展开，先计算得到以下各式：
y '=f (x , y ) y ''=f x +f y f
y
'''=f xx +2f xy f +f yy f 2+f y (f x +f y f )
再代入泰勒公式中：
h 2h 3
y (x n +1) =y (x n ) +h y '(x n ) +y ''(x n ) +y '''(x n ) +...
2! 3!
h 2h 3
=y n +hf +(f x +f y f ) +(f xx +2f xy f +f yy f 2
26
+f x f y +f y f ) +...
得到局部阶段误差为：
R n +1=y (x n +1) -y n +1
11
=h (1-c 1-c 2) f +h 2-c 2λ2) f x +(-c 2μ21) f y f ]
22
11231 +h -c 2λ2) f xx +(-c 2λ2μ21) f xy f 623112+(-c 2μ21) f yy f 2+f y (f x +f y f )]+... 62
要使局部截断误差为O (h 3) ，f , f x , f y f 的系数需为零，于是
c 1+c 2=11 21
c 2μ21=
2c 2λ2=
该方程组有4个未知数3个方程，所以有无穷多个解，它的每组解代入（2）式得到的近似公式都为二阶龙格-库塔方法。例如，取
1
c 1=c 2=, λ2=μ21=1，得到近似公式为：
2
h
y n +1=y n +(K 1+K 2)
2
K 1=f (x n , y n ) K 2=f (x n +h , y n +hK 1)
此处参数还有许多种不同的值，这里不再举例说明。
3、四阶龙格-库塔法
四阶龙格-库塔法的表达形式如下：
y n +1=y n +h (c 1K 1+c 2K 2+c 3K 3+c 4K 4) K 1=f (x n , y n )
K 2=f (x n +a 1h , y n +b 1K 1) K 3=f (x n +a 2h , y n +b 2K 1+b 3K 2) K 4=f (x n +a 3h , y n +b 4K 1+b 5K 2+b 6K 3)
仿照二阶龙格-库塔法的解法，通过N=4阶的泰勒级数的系数匹配，使局部截断误差为O (h 5) ，
b 1=a 1b 2+b 3=a 2b 4+b 5+b 6=a 3w 1+w 2+w 3+w 4=112122
w 2a 12+w 3a 2+w 4a 3=
3133
w 2a 13+w 3a 2+w 4a 3=
4w 2a 1+w 3a 2+w 4a 3=w 3a 1b 3+w 4(a 1b 5+a 2b 6) =
16
1
w 3a 1a 2b 3+w 4a 3(a 1b 5+a 2b 6) =
812
w 3a 12b 3+w 4(a 12b 5+a 2b 6) =
12
1
w 4a 1b 3b 6=
24
得到11个方程和13个未知量，因此需要补充两个条件
a 1=
1
, b 2=0，求解得到。 2
111
, a 3=1, b 1=, b 3=, b 4=0, b 5=0, 222
1111
b 6=1, w 1=, w 2=, w 3=, w 4=
6336a 2=
最后得出经典的四阶龙格-库塔公式：
h
y n +1=y n +(K 1+2K 2+2K 3+K 4)
6
K 1=f (x n , y n )
h h
K 2=f (x n +, y n +K 1)
22h h
K 3=f (x n +, y n +K 2)
22
K 4=f (x n +h , y n +hK 3)
三、数值计算方案
微分方程用积分可以解得火箭达到的高度，但是如果用数值方法解微分方程，所得的结果要精确得多。所以对于模型中的微分方程，用四阶龙格-库塔法对其求解所得结果精度更高，该方法还可以应用更为复杂的微分方程，如布朗动力学运动方程，但是对于高阶的微分方程，需对其降阶再。
km 千米, 火箭到上述模型中v 0=8km /s , 地球半径取x 0=R =6000
达最大高度处的速度为0，设当v i =0时所对应的火箭距离地球球心的高度为x i ，那么H =x i -R 即为火箭能达到的最大高度。 对于以下微分方程及其初值
⎧dv -gR 2-36⨯107⎪==2⎨dx vx vx 2, R ≤x ≤R +H ⎪⎩v 0=8, x 0=R =6000
7-36⨯10可知 f ( x , v ) = ，利用上述所求的四阶龙格-库塔公式得到vx 2
该模型的龙格-库塔式子：
h
v n +1=v n +(K 1+2K 2+2K 3+K 4)
6
K 1=f (x n , v n )
h h
K 2=f (x n +, v n +K 1)
22h h
K 3=f (x n +, v n +K 2)
22
K 4=f (x n +h , v n +hK 3)
根据题目，有初值v 0=6，x 0=6000, 即可求解. 取h=0.4
-36⨯107
K 1=f (x 0, v 0) ==-1. 25 2
8⨯(6000)
K 2==f (6000. 2, 7. 75) =-1. 2902 K 3=f (6000. 2, 7. 74196) =-1. 2905
K 4=f (6000. 4, 7. 4838) =-1. 3360
v 1=8-
......
0. 4
(1. 25+2⨯1. 2902+2⨯1. 2905+1. 3360) =7. 48351 6
v n =0
代入四阶龙格-库塔公式计算，直至v n
=0。步长选取为h =0. 4,
相应的x i =x 0+h , v i =对应的值，直到v i =0时，取到相应的x i ，则最大高度H =x i -R .
四、流程图
利用龙格-库塔法求解质点运动常微分方程
一、待解问题
讨论常微分方程的初值问题，边值问题的数值解法，最常用的基本方法就是龙格-库塔法使一些物理方程的计算简便，所得结果的精确提高。
下面是利用龙格-库塔法求解一个质点运动的常微分方程 以初速v 0=8km /s 自地球表面（半径R =6000km ）竖直向上发射一质量为 m 的火箭，如图所示。若不计空气阻力，火箭所受引力 F 大小与它到地心距离的平方成反比，求火箭所能到达的最大高度H 。
火箭为我们分析的对象，对其做受力分析，火箭在任意位置 x 处仅受地球引力F 作用。由题意知，F 的大小与 x2 成反比，设μ为比例系数，则有：
（1）
F =
μ
x 当火箭处于地面时，即x =R 时F = m g ，由式(1)可得μ=mgR 2，于
m gR 2
是火箭在任意位置 x 处所受地球引力 F 的大小为 F =2
x
由于火箭作直线运动，火箭的直线运动微分方程式为：
d 2x mgR 2m 2=-2dt x
分离变量积分
d 2x dv dv dx dv
==⋅=v 2
dt dx dx dt dt
dv mgR 2mv ⋅=-2
dx x
dv gR 2
=-2dx vx
二、物理机理
1. 龙格-库塔法的基本思想及一般形式
设初值问题y =y (x ) ∈[a , b ]，由微分中值定理可知，必存在ξ∈[x n , x n +1]，使
y (x n +1) =y (x n ) +h y '(ξ) =y (x n ) +hf (ξ, y (ξ))
设y n =y (x n ) ，并记K *=f (ξ, y (ξ)) ，则
y (x n +1) =y n +hK *
其中K *称为y (x ) 在[x n , x n +1]上的平均斜率，只要对平均斜率K *提供一种算法，上式就给出了一种数值解公式，例如，用K 1=f (x n , y n ) 代替K *，就得到欧拉公式，用K 2=f (x n +1, y n +1) 代替K *，就得到向后欧拉公式，如果用K 1, K 2的平均值来代替K *，则可得到二阶精度的梯形公式。可以设想，如果在[x n , x n +1]*库塔（Runge-Kutta ）法的基本思想。 龙格-库塔公式的一般形式：
y n +1=y n +h ∑c i k i
i =1
r
K 1=f (x n , y n )
K i =f (x n +λi h , y n +h ∑μij K j )
j =1i -1
（1）
其中K i 是y =y (x ) 在x n +λi h (0≤λi ≤1) 点的斜率预测值；c i , λi , μij 均为常数，选取这些常数的原则是使（1）式有尽可能高的精度。 2. 二阶龙格-库塔法
r =2的龙格-库塔公式为:
y n +1=y n +h (c 1K 1+c 2K 2) K 1=f (x n , y n )
K 2=f (x n +λ2h , y n +h μ21K 1)
（2）
式中c 1, c 2, λ2, μ21均为待定常数，我们希望适当的选取这些系数使公式的精度尽可能高，先将K 2按照二元函数泰勒级数：
1∂∂
f (x n +λ2h , y n +h μ21K 1) =∑(λ2h +h μ21) k f (x n , y n ) +...
∂x ∂y k =0k !
n
展开得到：
K 2=f (x n , y n ) +λ2hf x (x n , y n ) +μ21hf y (x n , y n ) f (x n , y n ) 122
[λ2h f xx (x n , y n ) +2λ2μ21h 2f xy (x n , y n ) f (x n , y n ) 2! 22+μ21h f yy (x n , y n ) f 2(x n , y n )]+... +
为了叙述方便，把f (x n , y n ) 和偏导数里面的x n , y n 省略不写，将K 2代入（2）式的第一式，得：
y n +1=y n +h (c 1+c 2) f +h 2c 2(λ2f x +μ21f y f ) c 222
+h (λ2f xx +2λ2μ21f xy f +μ21f yy f 2)
2
3
再用泰勒公式将y (x n +1) 展开，先计算得到以下各式：
y '=f (x , y ) y ''=f x +f y f
y
'''=f xx +2f xy f +f yy f 2+f y (f x +f y f )
再代入泰勒公式中：
h 2h 3
y (x n +1) =y (x n ) +h y '(x n ) +y ''(x n ) +y '''(x n ) +...
2! 3!
h 2h 3
=y n +hf +(f x +f y f ) +(f xx +2f xy f +f yy f 2
26
+f x f y +f y f ) +...
得到局部阶段误差为：
R n +1=y (x n +1) -y n +1
11
=h (1-c 1-c 2) f +h 2-c 2λ2) f x +(-c 2μ21) f y f ]
22
11231 +h -c 2λ2) f xx +(-c 2λ2μ21) f xy f 623112+(-c 2μ21) f yy f 2+f y (f x +f y f )]+... 62
要使局部截断误差为O (h 3) ，f , f x , f y f 的系数需为零，于是
c 1+c 2=11 21
c 2μ21=
2c 2λ2=
该方程组有4个未知数3个方程，所以有无穷多个解，它的每组解代入（2）式得到的近似公式都为二阶龙格-库塔方法。例如，取
1
c 1=c 2=, λ2=μ21=1，得到近似公式为：
2
h
y n +1=y n +(K 1+K 2)
2
K 1=f (x n , y n ) K 2=f (x n +h , y n +hK 1)
此处参数还有许多种不同的值，这里不再举例说明。
3、四阶龙格-库塔法
四阶龙格-库塔法的表达形式如下：
y n +1=y n +h (c 1K 1+c 2K 2+c 3K 3+c 4K 4) K 1=f (x n , y n )
K 2=f (x n +a 1h , y n +b 1K 1) K 3=f (x n +a 2h , y n +b 2K 1+b 3K 2) K 4=f (x n +a 3h , y n +b 4K 1+b 5K 2+b 6K 3)
仿照二阶龙格-库塔法的解法，通过N=4阶的泰勒级数的系数匹配，使局部截断误差为O (h 5) ，
b 1=a 1b 2+b 3=a 2b 4+b 5+b 6=a 3w 1+w 2+w 3+w 4=112122
w 2a 12+w 3a 2+w 4a 3=
3133
w 2a 13+w 3a 2+w 4a 3=
4w 2a 1+w 3a 2+w 4a 3=w 3a 1b 3+w 4(a 1b 5+a 2b 6) =
16
1
w 3a 1a 2b 3+w 4a 3(a 1b 5+a 2b 6) =
812
w 3a 12b 3+w 4(a 12b 5+a 2b 6) =
12
1
w 4a 1b 3b 6=
24
得到11个方程和13个未知量，因此需要补充两个条件
a 1=
1
, b 2=0，求解得到。 2
111
, a 3=1, b 1=, b 3=, b 4=0, b 5=0, 222
1111
b 6=1, w 1=, w 2=, w 3=, w 4=
6336a 2=
最后得出经典的四阶龙格-库塔公式：
h
y n +1=y n +(K 1+2K 2+2K 3+K 4)
6
K 1=f (x n , y n )
h h
K 2=f (x n +, y n +K 1)
22h h
K 3=f (x n +, y n +K 2)
22
K 4=f (x n +h , y n +hK 3)
三、数值计算方案
微分方程用积分可以解得火箭达到的高度，但是如果用数值方法解微分方程，所得的结果要精确得多。所以对于模型中的微分方程，用四阶龙格-库塔法对其求解所得结果精度更高，该方法还可以应用更为复杂的微分方程，如布朗动力学运动方程，但是对于高阶的微分方程，需对其降阶再。
km 千米, 火箭到上述模型中v 0=8km /s , 地球半径取x 0=R =6000
达最大高度处的速度为0，设当v i =0时所对应的火箭距离地球球心的高度为x i ，那么H =x i -R 即为火箭能达到的最大高度。 对于以下微分方程及其初值
⎧dv -gR 2-36⨯107⎪==2⎨dx vx vx 2, R ≤x ≤R +H ⎪⎩v 0=8, x 0=R =6000
7-36⨯10可知 f ( x , v ) = ，利用上述所求的四阶龙格-库塔公式得到vx 2
该模型的龙格-库塔式子：
h
v n +1=v n +(K 1+2K 2+2K 3+K 4)
6
K 1=f (x n , v n )
h h
K 2=f (x n +, v n +K 1)
22h h
K 3=f (x n +, v n +K 2)
22
K 4=f (x n +h , v n +hK 3)
根据题目，有初值v 0=6，x 0=6000, 即可求解. 取h=0.4
-36⨯107
K 1=f (x 0, v 0) ==-1. 25 2
8⨯(6000)
K 2==f (6000. 2, 7. 75) =-1. 2902 K 3=f (6000. 2, 7. 74196) =-1. 2905
K 4=f (6000. 4, 7. 4838) =-1. 3360
v 1=8-
......
0. 4
(1. 25+2⨯1. 2902+2⨯1. 2905+1. 3360) =7. 48351 6
v n =0
代入四阶龙格-库塔公式计算，直至v n
=0。步长选取为h =0. 4,
相应的x i =x 0+h , v i =对应的值，直到v i =0时，取到相应的x i ，则最大高度H =x i -R .
四、流程图