[1]  采用牛顿力学法，在直角坐标系下，建立空间三体混沌系统的运动微分方程；[2] 将其转换为一阶运动微分方程组；[3] 采用MATLAB ode45函数（一种基于显式Runge-Kutta的单步法求解器）求解。

1 关于三 体

1.1 什么是三体系统 ？

三体是天体力学中的基本力学模型。研究三个可视为质点的天体在相互之间万有引力作用下的运动规律。这三个天体的质量、初位置和初速度都是任意的。此时不考虑外力，整个系统质心保持惯性运动，系统质心动量守恒。

1.2 三体运动方程的“阶” ？

在一般三体问题中，每个质点具有3个方向的位置分量，3个方向的速度分量，因此系统共计(3+3）X 3 = 18个分量，三体是1个十八阶的常微分方程（实际上不大可能显示地写出这个方程），或者理解成18个一阶常微分方程组成的微分方程组，它们是等效的，并且后者对求解才有意义。

1.3 如何理解“混沌“ ？

首先“混沌“不等同于“混乱“，混沌是有序的，只不过受制于现有技术无法捕捉到一段时间后的真实状态。“混沌“的乱和我们小学三年级就知道的那个没有规律的无理数 \pi 完全不一样!

由于三体问题没有显示解，即不能写成u= f(t)这样的形式，所以无法给定一个时间，就计算出准确的状态。所以需要采用数值算法，建立递推公式，根据某个时刻的状态推测出下一个时刻的状态（该过程从零时刻开始），在这个递推过程中误差会被指数级放大，因此步长设置地再小，计算出的结果对未来很远的某个时刻都是没有意义的。

总结一下，就是混沌系统对初值极为敏感，只要有任何一点点误差，指数倍放大后都是一个天文数字，因此时间长了后结果不可信。比如天气系统就是一个混沌系统，所以长期的天气具有不可预测性，比如给你一个30天的天气预报，那基本就是闹着玩的。

2 建立运动微分方程

由于不受外力，质心可看作惯性系统，以其为参考点建立方程，这里为了简化求解，将参考点设置在原点。

在下面的推导中，粗体表示矢量，变量上带‘点’表示对时间求导，有几个‘点’就是几阶导数。

X1表示质点1的位置，F21表示质点2对1的万有引力，以此类推。

结合万有引力公式，得到质点1的运动微分方程：

同理可以得到其它质点的运动微分方程：

由于系统动量守恒（这里直接初始设置为0矢量），需要对初始条件加一个约束：

上文提到将质心和原点重合，便于计算，需要再加一个约束：

最终，得到建立好的运动方程：

为了采用数值解法，我们需要一个小技巧，将方程组转化为1阶的：

至此，我们建立了可以用于数值求解的运动微分方程组，是不是很简单？对于更多的星体，格式完全一样！

3 数值求解运动微分方程组

定义好系统参数和初值条件：

定义好迭代格式:

选择求解器：

这里采用匿名函数定义可以将更多的参数（本例中为: m1,m2,m3）传入ode45求解器。

至此，我们建立并且数值求解了三体运动的微分方程，如果你还不是太明白，可以点击阅读原文，查看视频教程。

接下来，图图将在极坐标系下推导单星与双星的运动微分方程 | 介绍龙格库塔迭代格式...

•完整代码：关注公众号“图通道”回复“三体”下载；

•MATLAB交流群：1129425848；