写在前面
之所以写这个程序, 是因为某天晚上无聊, 室友正在学习 MATLAB, 于是提议写一个三体运动的物理模拟程序来练练手. 就此, 我也写一份该程序来为室友做一个参考标准, 希望可以帮助室友进步的更快.
做出来的效果图大概这样子
效果图
本系列所有代码均在我的 GitHub 中存有备份, 可下载后直接运行, 点击 GitHub: HanpuLiang/Three-Body-by-MATLAB https://github.com/HanpuLiang/Three-Body-by-MATLAB 即可进入.
三体简介
三体一般指的就是三个物体受到相互之间的引力作用的影响而运动. 一般来说, 因为其运动方程太过于复杂, 所以并没有解析解, 并且因为对初值的敏感性, 略微变化一点初始条件就会对未来长远的结果产生巨大的影响.
在没有解析解的情况下, 只能通过数值解的方法对微分方程组求解. 所以数值解的误差也受计算步长的影响, 计算步长越小越精确, 但是因为数据一定会有精度, 并不能真正的无穷小, 所以实际上在时间足够长以后依旧会产生很大的误差.
综合很多原因, 才会有了大刘《三体》的剧情, 不然凭借三体人那么厉害的科技水平还怎么还是选择来搞地球.
不过说到底, 解不开这样的问题还是目前人类的数学水平不行, 或许以后就有办法了呢?
但是我们这里并不用分析力学的方法求解, 因为手头没有演草纸, 推方程有点麻烦, 所以直接用经典力学的方法去模拟整个运动, 这样子相信有点物理基础的大家也是可以看懂的.
运动过程分析
我们首先需要思考:
三个小球到底是怎么运动的? 引力作用.
小球运动, 哪些量在变化? 位置改变导致引力大小改变, 引力导致加速度改变, 加速度导致速度改变, 速度导致位置改变.
也就是说, 我们只需要集中在三个物理量上面就好: 坐标, 速度(大小与方向), 加速度(大小与方向). 这就是我们所需要, 随着时间变化的, 计算的所有数据.
接下来就要开始引进物理公式了.
两个物体之间的加速度
首先, 两个物体之间的万有引力可以通过公式
来计算, 其中 是物体 1 的质量,是物体 2 的质量,是引力系数 (模拟中为方便可以设为 1), 为物体 1 与物体 2 之间的直线距离.
根据力与加速度的公式 就可以得到, 在 时刻, 物体之间相对距离为 时(用上一次时间的距离算), 加速度为
两个物体之间的速度
在单位时间 内, 两个小球之间的速度变化量 为, 所以可以得到 时刻的速度为
两个小球的位置
令两个小球的坐标分别为, 则相对距离为
同理, 在 时刻, 受到速度由 变化到 的影响, 可以简单的得到此时刻的距离为
其中
矢量的正交分解
众所周知, 位置, 速度, 加速度均为矢量, 即存在方向和大小这两种属性, 也就是说我们需要考虑矢量方向不一致时的情况.
将矢量正交分解为 轴与 轴是最简单方便的做法.
正交分解示意图
首先是坐标, 这个已经是分解到了 轴与 轴这个坐标系上了, 毕竟我们写出来的就是两个点的坐标, 如果谁还不会用坐标点绘图就可以点右上角退出界面了.
其次是速度与加速度. 我们以物体自身为原点建立坐标系, 速度大小为 , 方向相对 轴正方向为 度, 可以得到一个矢量如图所示. 根据高中知识, 就可以得到其在 轴与 轴上的分解为
同样的, 加速度也可以这样子分解, 得到
而且,轴上的加速度只会影响 轴上的速度, 所以我们分解后, 在计算时, 只需要分别计算 轴的坐标变化即可, 不需要再考虑方向, 即
这样子, 我们就将方向成功分解为 轴分解进行计算, 大大化简了繁琐的方向变化问题.
矢量的叠加
但是这只是两个物体之间的相互作用, 如果是三个物体的话, 其中一个物体就要受到两个力的作用.
实际上两个力是没有受到干扰的, 所以当其分解到 轴后, 直接将其对应轴上的加速度直接相加即可得到总的加速度, 也就是
其中 就是物体 1 在 轴上的总的加速度, 它由两个分加速度组成: 来自物体 2 对物体 1 的力的, 在 轴的加速度 和来自物体 3 对物体 1 的.
其他同理, 这样子就可以完美解决所有问题了.
代码思路
根据上面的公式分析, 加速度, 速度, 距离之间如何变化已经很清楚了, 三个物体之间的各个物理量的正交分解也很明确了, 已经可以转化为了代码可以实现的情况, 下面我们就需要将公式化成代码.
不过因为这一篇博客已经比较长了, 所以将本篇作为理论分析篇, 下一篇博客中我们再进行详细解释代码.
来源: http://www.jianshu.com/p/d1a56bf54a6c