如何寫出三體的MATLAB程序-理論分析篇

如何寫出三體的MATLAB程序-理論分析篇

寫在前面

之所以寫這個程序，是因為某天晚上無聊，室友正在學習MATLAB，于是提議寫一個三體運動的物理模擬程序來練練手。就此，我也寫一份該程序來為室友做一個參考標準，希望可以幫助室友進步的更快。

做出來的效果圖大概這樣子

效果圖

本系列所有代碼均在我的Github中存有備份，可下載后直接運行，點擊Github: HanpuLiang/Three-Body-by-MATLAB即可進入。

三體簡介

三體一般指的就是三個物體受到相互之間的引力作用的影響而運動。一般來說，因為其運動方程太過于復雜，所以并沒有解析解，并且因為對初值的敏感性，略微變化一點初始條件就會對未來長遠的結(jié)果產(chǎn)生巨大的影響。

在沒有解析解的情況下，只能通過數(shù)值解的方法對微分方程組求解。所以數(shù)值解的誤差也受計算步長的影響，計算步長越小越精確，但是因為數(shù)據(jù)一定會有精度，并不能真正的無窮小，所以實際上在時間足夠長以后依舊會產(chǎn)生很大的誤差。

綜合很多原因，才會有了大劉《三體》的劇情，不然憑借三體人那么厲害的科技水平還怎么還是選擇來搞地球。

不過說到底，解不開這樣的問題還是目前人類的數(shù)學水平不行，或許以后就有辦法了呢？

但是我們這里并不用分析力學的方法求解，因為手頭沒有演草紙，推方程有點麻煩，所以直接用經(jīng)典力學的方法去模擬整個運動，這樣子相信有點物理基礎的大家也是可以看懂的。

運動過程分析

我們首先需要思考：

  • 三個小球到底是怎么運動的？引力作用。
  • 小球運動，哪些量在變化？位置改變導致引力大小改變，引力導致加速度改變，加速度導致速度改變，速度導致位置改變。

也就是說，我們只需要集中在三個物理量上面就好：坐標，速度(大小與方向)，加速度(大小與方向)。這就是我們所需要，隨著時間變化的，計算的所有數(shù)據(jù)。

接下來就要開始引進物理公式了。

兩個物體之間的加速度

首先，兩個物體之間的萬有引力可以通過公式

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

來計算，其中m_1是物體1的質(zhì)量，m_2是物體2的質(zhì)量，G是引力系數(shù)(模擬中為方便可以設為1)，r為物體1與物體2之間的直線距離。

根據(jù)力與加速度的公式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內(nèi)，兩個小球之間的速度變化量\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軸是最簡單方便的做法。

正交分解示意圖

首先是坐標，這個已經(jīng)是分解到了x軸與y軸這個坐標系上了，畢竟我們寫出來的就是兩個點的坐標，如果誰還不會用坐標點繪圖就可以點右上角退出界面了。

其次是速度與加速度。我們以物體自身為原點建立坐標系，速度大小為v，方向相對x軸正方向為\theta度，可以得到一個矢量如圖所示。根據(jù)高中知識，就可以得到其在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}。

其他同理，這樣子就可以完美解決所有問題了。

代碼思路

根據(jù)上面的公式分析，加速度、速度、距離之間如何變化已經(jīng)很清楚了，三個物體之間的各個物理量的正交分解也很明確了，已經(jīng)可以轉(zhuǎn)化為了代碼可以實現(xiàn)的情況，下面我們就需要將公式化成代碼。

不過因為這一篇博客已經(jīng)比較長了，所以將本篇作為理論分析篇，下一篇博客中我們再進行詳細解釋代碼。

如果這一篇我講的比較不錯的話，還希望可以點個贊、加個收藏、來個關(guān)注噢。

最后編輯于 ：
?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請聯(lián)系作者
【社區(qū)內(nèi)容提示】社區(qū)部分內(nèi)容疑似由AI輔助生成，瀏覽時請結(jié)合常識與多方信息審慎甄別。
平臺聲明：文章內(nèi)容（如有圖片或視頻亦包括在內(nèi)）由作者上傳并發(fā)布，文章內(nèi)容僅代表作者本人觀點，簡書系信息發(fā)布平臺，僅提供信息存儲服務。

相關(guān)閱讀更多精彩內(nèi)容

友情鏈接更多精彩內(nèi)容