公式(5)、(6)所表示的解有很大的实用价值,在后面的例题和练习中希望大家不断的学习和体会,直到能够熟练运用这一方法。我们可以通过MATLAB编写一个自定义函数rk4( )来实现以上内容,并保存为rk4.m(注 ),具体脚本语句如下:
2.2 建立计算模型
前面提到我们将用四阶龙格-库塔方法来求解行星的运动微分方程,在MATLAB程序中我们首先要给出初始位置和初始速度以及在程序中要用到的一些常数:
从以上程序中可以看到,我们首先设定了初始位置和速度,为了简单方便,我们将初始位置定在x轴上,由于速度方向总是沿轨道的切线方向,所以相应的初始速度只有y分量。同时我们还把初始速度的大小设定为2的h倍,这样我们就能通过控制h来得到不同的轨道形状。例如,若设定r0=1AU,则当h取1时,行星轨道为圆形(参见3.2.4节内容),如果h取其他值,我们就能对应得到椭圆、抛物线或者双曲线轨道。