本博客来自北理工《计算机仿真与matlab》课程的大作业，有任何问题欢迎批评指正。 源代码请参见如下链接： 链接：https://pan.baidu.com/s/1hkAoi9c16_IoUYJSxX92hg 提取码：4gr6

我使用matlab和simulink建立天体系统的模型，并使用不同的数值积分方法仿真求解，模拟推演天体运动轨迹。

下文第3、4节将演示双星系统的仿真模拟，第5节将演示三体系统的仿真模拟。

行星所受的万有引力为：

F ⃗ = G M m ∣ R ⃗ − r ⃗ ∣ 2 ⋅ R ⃗ − r ⃗ ∣ R ⃗ − r ⃗ ∣ \vec{F}=\frac{GMm}{\left|\vec{R}-\vec{r}\right|^2}·\frac{\vec{R}-\vec{r}}{\left|\vec{R}-\vec{r}\right|} F =∣∣∣​R −r ∣∣∣​2GMm​⋅∣∣∣​R −r ∣∣∣​R −r ​

其中 G G G、 M M M、 m m m 分别为万有引力常量、恒星质量、行星质量， R ⃗ \vec{R} R 是恒星位置， r ⃗ \vec{r} r 是行星位置。

由牛二律建立行星运动微分方程：

d 2 r ⃗ d t 2 = G M ( R ⃗ − r ⃗ ) ∣ R ⃗ − r ⃗ ∣ 3 \frac{d^2\vec{r}}{dt^2}=\frac{GM(\vec{R}-\vec{r})}{\left|\vec{R}-\vec{r}\right|^3} dt2d2r ​=∣∣∣​R −r ∣∣∣​3GM(R −r )​

考虑到恒星的质量远大于行星的质量，假设恒星在系统中静止，即 R ⃗ \vec{R} R 为常矢量。因此上式为关于行星位置 r ⃗ \vec{r} r 的二阶微分方程，求解时将 r ⃗ \vec{r} r 分解到x、y、z轴，得到与上式等价的三个微分方程组。

函数pre_load_twobody.m，预定义仿真参数（仿真时间间隔h），天体参数（万有引力常量G和天体质量）和初始条件（天体的初始位置和初始速度）。

twobody.m，使用梯形法求解微分方程并绘制天体轨迹。（其中调用的count_a函数用于计算加速度，此处略，请参见源代码）