（建议阅读全文）
预备知识 万有引力， 弹簧振子受迫运动的简单数值计算 下面我们来用一种极其简单的算法对单个天体在中心天体的万有引力作用下的运动进行数值计算． 事实上该问题存在解析解（见开普勒三定律）， 所以以下的算法只是用于演示数值解常微分方程的大致原理． 这种方法可以轻易地拓展到多个天体的情况， 而多体情况没有一般的解析解． 一种效率更高的常见常微分方程数值算法可参考 “四阶龙格库塔法”． 直角坐标系中， 设中心天体质量为
其中
用
假设已知初值条件
2．设经过一段极微小的时间步长
3．把
% 参数设定
GM = 1; % 万有引力常数乘以中心天体质量
x0 = 1; y0 = 0; % 初始位置
vx0 = 0; vy0 = 0.7; % 初始速度
T = 4; Nstep = 4000; % 总时间和步数
dt = T/Nstep; % 步长% 矩阵预赋值
x = nan(Nstep,1); y = x;
x1 = x; y1 = x;
x2 = x; y2 = x;% 初始位置，速度，加速度
x(1) = x0; y(1) = y0; % 初位置
x1(1) = vx0; y1(1) = vy0; % 初速度
x2(1) = -GM*x(1)/(x(1)^2+y(1)^2)^(3/2); % 代入方程得到 x''(0)
y2(1) = -GM*y(1)/(x(1)^2+y(1)^2)^(3/2); % 代入方程得到 y''(0)% 迭代循环
for ii = 2:Nstepx(ii) = x(ii-1)+x1(ii-1)*dt; % x的微分y(ii) = y(ii-1)+y1(ii-1)*dt; % y的微分x1(ii) = x1(ii-1)+x2(ii-1)*dt; % x' 的微分y1(ii) = y1(ii-1)+y2(ii-1)*dt; % y' 的微分x2(ii) = -GM*x(ii)/(x(ii)^2+y(ii)^2)^(3/2); % 代入微分方程求出 x''y2(ii) = -GM*y(ii)/(x(ii)^2+y(ii)^2)^(3/2); % 代入微分方程求出 y''
end% 画图
plot(x,y); % 画行星轨道
axis equal; % xy坐标长度一致
hold on; % 继续画图
scatter(0,0); % 标出中心天体程序运行结果如图 1 所示． 注意行星轨道并不是一个闭合的椭圆， 这是由于这种算法误差较大， 为了减小误差， 可以增加程序中 Nstep 的值． 
