（建议阅读原文）预备知识 中点法解常微分方程（组）龙格库塔法是一类数值解微分方程的算法， 其中较常见的是四阶龙格库塔法． 这里不进行推导， 仅仅给出公式如下（
![]() 的定义类比式 5 ）
其中
由以上两式， 不难把该算法拓展到方程组的情况． 对于
![]() 元微分方程组
我们可以把该式记为矢量函数的形式
现在我们仅需要把式 1 和 式 2 中的所有
![]() 和
![]() 都变为
![]() 维列矢量
![]() 和
![]() 即可将微分方程拓展为微分方程组．
例程
到目前为止， 我们每求一个微分方程的数值解都要重新写一次程序， 对于一些较为复杂的算法这样做效率较低． 我们这里不妨把四阶龙格库塔法写到一个单独的函数文件 odeRK4.m 中， 当我们要解某个特定的方程时， 只需把式 2 中的
![]() 作为自变量输入即可解出
![]() ．
代码 1：odeRK4.m
我们先来看第 1 行的函数声明， 输入变量中，f 是式 4 中
![]() 的函数句柄，
tspan 是一个
2×1 的列矢量，
tspan(1) 是初始时间，
tspan(2) 是终止时间，
Y0 是一个列矢量，
Y0(ii) 是第
ii 个因变量的初始值，
Nt 是
![]() 的个数，
tspan 定义的时间区间被等分为
Nt - 1 个小区间． 因变量中，
Y 的行数是因变量的个数， 列数是
Nt，
t 是一个行矢量， 由第 6 行定义，
Y(ii, jj) 就是第
ii 个变量在
t(jj) 时刻的值． 第 5 行把初值
Y0 赋给
Y 的第 1 列， 第 8-14 行的循环根据式 1 和式 2 的矢量形式由
Y 的第
ii 列（
![]() ）求第
ii+1 列（
![]() ）．
我们先来用这个函数来计算 “天体运动的简单数值计算” 中的问题． 我们令因变量
![]() 的四个分量依次为一阶方程组（式 4 ）
中的
![]() ． 程序代码如下
代码 2：keplerRK4.m
运行结果如图 1 所示．
图 1：运行结果
我们先来看函数 fun （20 行）， 这个函数就相当于式 5 ． 第一个输入变量 Y 是一个列矢量， 是
![]() 的值， 第二个输入变量是
![]() ， 但由于式 5 中没有出现
![]() ， 我们用波浪线代替． 第三个输入变量是参数
![]() ， 即万有引力常数和中心天体质量之积． 输出变量
Y1 是一个列矢量， 是
![]() 的值．
再来看主函数
KeplerRK4 （第 1 行）， 参数设定中除了步数从 4000 变为了 100， 其他都和 “天体运动的简单数值计算” 中的程序一样， 然而这里运行结果却精确得多（曲线几乎闭合）， 可见这种算法的优越性．
主函数第 10 行中将
fun(Y, t, GM) 变为函数句柄
f(Y, t)， 这样
GM 就可以在 “参数设定” 中设置， 而不用在
fun 函数内部设置． 第 11 行调用了上文中的
odeRK4 函数解方程组， 由于我们在画图时不需要用
t， 所以把第二个输出变量改为波浪线．
| 
转播
