利用龙格库塔法求解质点运动方程

一、待解问题

讨论常微分方程的初值问题，边值问题的数值解法，最常用的基本方法就是龙格-库塔法使一些物理方程的计算简便，所得结果的精确提高。

下面是利用龙格-库塔法求解一个质点运动的常微分方程 以初速v08km/s自地球表面（半径R6000km）竖直向上发射一质量为 m 的火箭，如图所示。若不计空气阻力，火箭所受引力 F 大小与它到地心距离的平方成反比，求火箭所能到达的最大高度H。

火箭为我们分析的对象，对其做受力分析，火箭在任意位置 x 处仅受地球引力F 作用。由题意知，F 的大小与 x2 成反比，设为比例系数，则有：

F2 （1）

x当火箭处于地面时，即xR时F = m g ，由式(1)可得mgR2，于

mgR2是火箭在任意位置 x 处所受地球引力 F 的大小为 F2

x由于火箭作直线运动，火箭的直线运动微分方程式为：

d2xmgR2m22dtx分离变量积分

d2xdvdvdxdvv2dtdxdxdt dt2dvmgR

mv

dxx2dvgR22dxvx

二、物理机理

1.龙格-库塔法的基本思想及一般形式

设初值问题yy(x)[a,b]，由微分中值定理可知，必存在[xn,xn1]，使

y(xn1)y(xn)hy()y(xn)hf(,y())设yny(xn)，并记K*f(,y())，则

y(xn1)ynhK*

其中K*称为y(x)在[xn,xn1]上的平均斜率，只要对平均斜率K*提供一种算法，上式就给出了一种数值解公式，例如，用K1f(xn,yn)代替K*，就得到欧拉公式，用K2f(xn1,yn1)代替K*，就得到向后欧拉公式，如果用K1,K2的平均值来代替K*，则可得到二阶精度的梯形公式。可以设想，如果在[xn,xn1]上能多预测几个点的斜率值，用它们的加权平均值代替K*，就有望的到具有较高精度的数值解公式，这就是龙格-库塔（Runge-Kutta）法的基本思想。 龙格-库塔公式的一般形式：

yn1ynhcikii1r

K1f(xn,yn)Kif(xnih,ynhijKj)j1i1 （1）

其中Ki是yy(x)在xnih(0i1)点的斜率预测值；ci,i,ij均为常数，选取这些常数的原则是使（1）式有尽可能高的精度。 2. 二阶龙格-库塔法

r2的龙格-库塔公式为:

yn1ynh(c1K1c2K2)K1f(xn,yn)K2f(xn2h,ynh21K1) （2）

式中c1,c2,2,21均为待定常数，我们希望适当的选取这些系数使公式的精度尽可能高，先将K2按照二元函数泰勒级数：

1f(xn2h,ynh21K1)(2hh21)kf(xn,yn)...

xyk0k!n展开得到：

K2f(xn,yn)2hfx(xn,yn)21hfy(xn,yn)f(xn,yn)122[2hfxx(xn,yn)2221h2fxy(xn,yn)f(xn,yn)2!2221hfyy(xn,yn)f2(xn,yn)]...

为了叙述方便，把f(xn,yn)和偏导数里面的xn,yn省略不写，将K2代入（2）式的第一式，得：

yn1ynh(c1c2)fh2c2(2fx21fyf)c222h(2fxx2221fxyf21fyyf2)23再用泰勒公式将y(xn1)展开，先计算得到以下各式：

yf(x,y)yfxfyfyfxx2fxyffyyf2fy(fxfyf)再代入泰勒公式中：

h2h3y(xn1)y(xn)hy(xn)y(xn)y(xn)...2!3!h2h3ynhf(fxfyf)(fxx2fxyffyyf2

26fxfyfyf)...得到局部阶段误差为：

Rn1y(xn1)yn111h(1c1c2)fh2[(c22)fx(c221)fyf]2211231 h[(c22)fxx(c2221)fxyf623112(c221)fyyf2fy(fxfyf)]...62要使局部截断误差为O(h3)，f,fx,fyf的系数需为零，于是

c1c211 21c2212c22该方程组有4个未知数3个方程，所以有无穷多个解，它的每组解代入（2）式得到的近似公式都为二阶龙格-库塔方法。例如，取

1c1c2,2211，得到近似公式为：

2

hyn1yn(K1K2)2K1f(xn,yn)K2f(xnh,ynhK1)此处参数还有许多种不同的值，这里不再举例说明。

3、四阶龙格-库塔法

四阶龙格-库塔法的表达形式如下：

yn1ynh(c1K1c2K2c3K3c4K4)K1f(xn,yn)K2f(xna1h,ynb1K1)K3f(xna2h,ynb2K1b3K2)K4f(xna3h,ynb4K1b5K2b6K3)仿照二阶龙格-库塔法的解法，通过N=4阶的泰勒级数的系数匹配，使局部截断误差为O(h5)，

b1a1b2b3a2b4b5b6a3w1w2w3w4112122w2a12w3a2w4a33133w2a13w3a2w4a34w2a1w3a2w4a3w3a1b3w4(a1b5a2b6)1618

w3a1a2b3w4a3(a1b5a2b6)1wabw4(abab)121w4a1b3b6242313215226得到11个方程和13个未知量，因此需要补充两个条件

1a1,b20，求解得到。

2

111a2,a31,b1,b3,b40,b50,2221111b61,w1,w2,w3,w46336最后得出经典的四阶龙格-库塔公式：

hyn1yn(K12K22K3K4)6K1f(xn,yn)hhK2f(xn,ynK1)22hhK3f(xn,ynK2)22K4f(xnh,ynhK3)

三、数值计算方案

微分方程用积分可以解得火箭达到的高度，但是如果用数值方法解微分方程，所得的结果要精确得多。所以对于模型中的微分方程，用四阶龙格-库塔法对其求解所得结果精度更高，该方法还可以应用更为复杂的微分方程，如布朗动力学运动方程，但是对于高阶的微分方程，需对其降阶再。

上述模型中v08km/s,地球半径取x0R6000km千米,火箭到达最大高度处的速度为0，设当vi0时所对应的火箭距离地球球心的高度为xi，那么HxiR即为火箭能达到的最大高度。

对于以下微分方程及其初值

dvgR236107 2dxvxvx2,RxRHv8,xR6000  0 0

73610可知 f ( x , v )  ，利用上述所求的四阶龙格-库塔公式得到vx2该模型的龙格-库塔式子：

hvn1vn(K12K22K3K4)6K1f(xn,vn)hhK2f(xn,vnK1)22hhK3f(xn,vnK2)22K4f(xnh,vnhK3)

根据题目，有初值v06，x06000, 即可求解. 取h=0.4

36107K1f(x0,v0)1.25 28(6000)K2f(6000.2,7.75)1.2902

K3f(6000.2,7.74196)1.2905

K4f(6000.4,7.4838)1.3360

0.4v18(1.2521.290221.29051.3360)7.48351

6......

vn0

代入四阶龙格-库塔公式计算，直至vn0。步长选取为h0.4,

相应的xix0h,vi对应的值，直到vi0时，取到相应的xi，则最大高度HxiR.

四、流程图

初始值v08,x06000,h0.4求出vi0将vi,xi再次代入 K1,K2,vi,xi vi0得出HxiR