量子物理·宇宙天文·电磁场·光·万有引力物理程序员

如何写出三体的MATLAB程序-理论分析篇

2019-03-07 本文已影响5人老梁家的风子

如何写出三体的MATLAB程序-理论分析篇

写在前面

之所以写这个程序，是因为某天晚上无聊，室友正在学习MATLAB，于是提议写一个三体运动的物理模拟程序来练练手。就此，我也写一份该程序来为室友做一个参考标准，希望可以帮助室友进步的更快。

做出来的效果图大概这样子

效果图

本系列所有代码均在我的Github中存有备份，可下载后直接运行，点击Github: HanpuLiang/Three-Body-by-MATLAB即可进入。

三体简介

三体一般指的就是三个物体受到相互之间的引力作用的影响而运动。一般来说，因为其运动方程太过于复杂，所以并没有解析解，并且因为对初值的敏感性，略微变化一点初始条件就会对未来长远的结果产生巨大的影响。

在没有解析解的情况下，只能通过数值解的方法对微分方程组求解。所以数值解的误差也受计算步长的影响，计算步长越小越精确，但是因为数据一定会有精度，并不能真正的无穷小，所以实际上在时间足够长以后依旧会产生很大的误差。

综合很多原因，才会有了大刘《三体》的剧情，不然凭借三体人那么厉害的科技水平还怎么还是选择来搞地球。

不过说到底，解不开这样的问题还是目前人类的数学水平不行，或许以后就有办法了呢？

但是我们这里并不用分析力学的方法求解，因为手头没有演草纸，推方程有点麻烦，所以直接用经典力学的方法去模拟整个运动，这样子相信有点物理基础的大家也是可以看懂的。

运动过程分析

我们首先需要思考：

三个小球到底是怎么运动的？引力作用。
小球运动，哪些量在变化？位置改变导致引力大小改变，引力导致加速度改变，加速度导致速度改变，速度导致位置改变。

也就是说，我们只需要集中在三个物理量上面就好：坐标，速度(大小与方向)，加速度(大小与方向)。这就是我们所需要，随着时间变化的，计算的所有数据。

接下来就要开始引进物理公式了。

两个物体之间的加速度

首先，两个物体之间的万有引力可以通过公式

$F=G\dfrac{m_1m_2}{r^2}$

来计算，其中 $m_1$ 是物体1的质量， $m_2$ 是物体2的质量， $G$ 是引力系数(模拟中为方便可以设为1)， $r$ 为物体1与物体2之间的直线距离。

根据力与加速度的公式 $F=ma$ 就可以得到，在 $t$ 时刻，物体之间相对距离为 $r(t-1)$ 时(用上一次时间的距离算)，加速度为

$\text{物体1的加速度: }\quad a_1(t)=G\dfrac{m_2}{r(t-1)^2}$

$\text{物体2的加速度: }\quad a_2(t)=G\dfrac{m_1}{r(t-1)^2}$

两个物体之间的速度

在单位时间 $\Delta t$ 内，两个小球之间的速度变化量 $\Delta v$ 为 $\Delta v=a\Delta t$ ，所以可以得到 $t$ 时刻的速度为

$\text{物体1的速度: }\quad v_1(t)=v_1(t-1) + a_1(t)\Delta t$

$\text{物体2的速度: }\quad v_2(t)=v_2(t-1) + a_2(t)\Delta t$

两个小球的位置

令两个小球的坐标分别为 $(x_1, y_1), (x_2, y_2)$ ，则相对距离为

$r = \sqrt{(x_1-x_2)^2+(y_1-y_2)^2}$

同理，在 $t$ 时刻，受到速度由 $v(t-1)$ 变化到 $v(t)$ 的影响，可以简单的得到此时刻的距离为

$r(t) = r(t-1) + \Delta r$

其中

$\Delta r = \dfrac12(v(t-1) + v(t))\Delta t.$

矢量的正交分解

众所周知，位置、速度、加速度均为矢量，即存在方向和大小这两种属性，也就是说我们需要考虑矢量方向不一致时的情况。

将矢量正交分解为 $x$ 轴与 $y$ 轴是最简单方便的做法。

正交分解示意图

首先是坐标，这个已经是分解到了 $x$ 轴与 $y$ 轴这个坐标系上了，毕竟我们写出来的就是两个点的坐标，如果谁还不会用坐标点绘图就可以点右上角退出界面了。

其次是速度与加速度。我们以物体自身为原点建立坐标系，速度大小为 $v$ ，方向相对 $x$ 轴正方向为 $\theta$ 度，可以得到一个矢量如图所示。根据高中知识，就可以得到其在 $x$ 轴与 $y$ 轴上的分解为

$v_x(t) = v(t)\cos\theta$

$v_y(t) = v(t)\sin\theta,$

同样的，加速度也可以这样子分解，得到

$a_x(t) = a(t)\cos\theta$

$a_y(t) = a(t)\sin\theta.$

而且， $x$ 轴上的加速度只会影响 $x$ 轴上的速度，所以我们分解后，在计算时，只需要分别计算 $xy$ 轴的坐标变化即可，不需要再考虑方向，即

$v_x(t) = v_x(t-1) + a_x(t)\Delta t.$

这样子，我们就将方向成功分解为 $xy$ 轴分解进行计算，大大化简了繁琐的方向变化问题。

矢量的叠加

但是这只是两个物体之间的相互作用，如果是三个物体的话，其中一个物体就要受到两个力的作用。

实际上两个力是没有受到干扰的，所以当其分解到 $xy$ 轴后，直接将其对应轴上的加速度直接相加即可得到总的加速度，也就是

$a_{1x}(t) = a_{12x}(t) + a_{13x}(t),$

其中 $a_{1x}$ 就是物体1在 $x$ 轴上的总的加速度，它由两个分加速度组成:来自物体2对物体1的力的、在 $x$ 轴的加速度 $a_{12x}$ 和来自物体3对物体1的 $x_{13x}$ 。

其他同理，这样子就可以完美解决所有问题了。

代码思路

根据上面的公式分析，加速度、速度、距离之间如何变化已经很清楚了，三个物体之间的各个物理量的正交分解也很明确了，已经可以转化为了代码可以实现的情况，下面我们就需要将公式化成代码。

不过因为这一篇博客已经比较长了，所以将本篇作为理论分析篇，下一篇博客中我们再进行详细解释代码。

如果这一篇我讲的比较不错的话，还希望可以点个赞、加个收藏、来个关注噢。

上一篇下一篇

猜你喜欢

热点阅读