?[1]??采用牛頓力學(xué)法，在直角坐標(biāo)系下，建立空間三體混沌系統(tǒng)的運(yùn)動(dòng)微分方程；[2] 將其轉(zhuǎn)換為一階運(yùn)動(dòng)微分方程組；[3] 采用MATLAB?ode45函數(shù)（一種基于顯式Runge-Kutta的單步法求解器）求解。

1 關(guān)于三 體

1.1 什么是三體系統(tǒng) ？

三體是天體力學(xué)中的基本力學(xué)模型。研究三個(gè)可視為質(zhì)點(diǎn)的天體在相互之間萬(wàn)有引力作用下的運(yùn)動(dòng)規(guī)律。這三個(gè)天體的質(zhì)量、初位置和初速度都是任意的。此時(shí)不考慮外力，整個(gè)系統(tǒng)質(zhì)心保持慣性運(yùn)動(dòng)，系統(tǒng)質(zhì)心動(dòng)量守恒。

1.2 三體運(yùn)動(dòng)方程的“階” ？

在一般三體問(wèn)題中，每個(gè)質(zhì)點(diǎn)具有3個(gè)方向的位置分量，3個(gè)方向的速度分量，因此系統(tǒng)共計(jì)(3+3）X 3 = 18個(gè)分量，三體是1個(gè)十八階的常微分方程（實(shí)際上不大可能顯示地寫(xiě)出這個(gè)方程），或者理解成18個(gè)一階常微分方程組成的微分方程組，它們是等效的，并且后者對(duì)求解才有意義。

1.3 如何理解“混沌“ ？

首先“混沌“不等同于“混亂“，混沌是有序的，只不過(guò)受制于現(xiàn)有技術(shù)無(wú)法捕捉到一段時(shí)間后的真實(shí)狀態(tài)。“混沌“的亂和我們小學(xué)三年級(jí)就知道的那個(gè)沒(méi)有規(guī)律的無(wú)理數(shù) \pi 完全不一樣!

由于三體問(wèn)題沒(méi)有顯示解，即不能寫(xiě)成u= f(t)這樣的形式，所以無(wú)法給定一個(gè)時(shí)間，就計(jì)算出準(zhǔn)確的狀態(tài)。所以需要采用數(shù)值算法，建立遞推公式，根據(jù)某個(gè)時(shí)刻的狀態(tài)推測(cè)出下一個(gè)時(shí)刻的狀態(tài)（該過(guò)程從零時(shí)刻開(kāi)始），在這個(gè)遞推過(guò)程中誤差會(huì)被指數(shù)級(jí)放大，因此步長(zhǎng)設(shè)置地再小，計(jì)算出的結(jié)果對(duì)未來(lái)很遠(yuǎn)的某個(gè)時(shí)刻都是沒(méi)有意義的。

總結(jié)一下，就是混沌系統(tǒng)對(duì)初值極為敏感，只要有任何一點(diǎn)點(diǎn)誤差，指數(shù)倍放大后都是一個(gè)天文數(shù)字，因此時(shí)間長(zhǎng)了后結(jié)果不可信。比如天氣系統(tǒng)就是一個(gè)混沌系統(tǒng)，所以長(zhǎng)期的天氣具有不可預(yù)測(cè)性，比如給你一個(gè)30天的天氣預(yù)報(bào)，那基本就是鬧著玩的。

2 建立運(yùn)動(dòng)微分方程

由于不受外力，質(zhì)心可看作慣性系統(tǒng)，以其為參考點(diǎn)建立方程，這里為了簡(jiǎn)化求解，將參考點(diǎn)設(shè)置在原點(diǎn)。

在下面的推導(dǎo)中，粗體表示矢量，變量上帶‘點(diǎn)’表示對(duì)時(shí)間求導(dǎo)，有幾個(gè)‘點(diǎn)’就是幾階導(dǎo)數(shù)。

X1表示質(zhì)點(diǎn)1的位置，F(xiàn)21表示質(zhì)點(diǎn)2對(duì)1的萬(wàn)有引力，以此類(lèi)推。

結(jié)合萬(wàn)有引力公式，得到質(zhì)點(diǎn)1的運(yùn)動(dòng)微分方程：

同理可以得到其它質(zhì)點(diǎn)的運(yùn)動(dòng)微分方程：

由于系統(tǒng)動(dòng)量守恒（這里直接初始設(shè)置為0矢量），需要對(duì)初始條件加一個(gè)約束：

上文提到將質(zhì)心和原點(diǎn)重合，便于計(jì)算，需要再加一個(gè)約束：

最終，得到建立好的運(yùn)動(dòng)方程：

為了采用數(shù)值解法，我們需要一個(gè)小技巧，將方程組轉(zhuǎn)化為1階的：

至此，我們建立了可以用于數(shù)值求解的運(yùn)動(dòng)微分方程組，是不是很簡(jiǎn)單？對(duì)于更多的星體，格式完全一樣！

3 數(shù)值求解運(yùn)動(dòng)微分方程組

定義好系統(tǒng)參數(shù)和初值條件：

定義好迭代格式:

選擇求解器：

這里采用匿名函數(shù)定義可以將更多的參數(shù)（本例中為: m1,m2,m3）傳入ode45求解器。

至此，我們建立并且數(shù)值求解了三體運(yùn)動(dòng)的微分方程，如果你還不是太明白，可以點(diǎn)擊閱讀原文，查看視頻教程。

接下來(lái)，圖圖將在極坐標(biāo)系下推導(dǎo)單星與雙星的運(yùn)動(dòng)微分方程 | 介紹龍格庫(kù)塔迭代格式...

?完整代碼：關(guān)注公眾號(hào)“圖通道”回復(fù)“三體”下載；

?MATLAB交流群：1129425848；