全文共 7726 字，預(yù)計(jì)學(xué)習(xí)時(shí)長 15 分鐘或更長
圖片來自flickr, 凱文·吉爾
中國作家劉慈欣的科幻小說《三體》中描繪了存在于被三顆恒星環(huán)繞的“三體”星球上的一種虛構(gòu)外星文明。能想象這種文明的存在因三顆恒星而和我們的文明大不相同嗎？炫目的陽光？持續(xù)的夏日？事實(shí)證明，情況要糟糕很多。
生活在僅有一顆主要恒星的太陽系是值得慶幸的，因?yàn)檫@使得這顆恒星（太陽）的軌道有可預(yù)測(cè)性。即使增加一顆恒星，這個(gè)系統(tǒng)仍能保持穩(wěn)定。該系統(tǒng)有個(gè)被稱為分析解的解法，即描繪解方程式，并得到可以精確提供該系統(tǒng)從一秒到百萬年時(shí)間演變的函數(shù)。
然而，當(dāng)加入第三體時(shí)，就會(huì)發(fā)生特殊情況。系統(tǒng)會(huì)變得混亂而不可預(yù)測(cè)。沒有分析解法（除了少數(shù)特定情況），并且只能在計(jì)算機(jī)上用數(shù)字解方程式。它們會(huì)突然從穩(wěn)定變到不穩(wěn)定，或者從不穩(wěn)定變到穩(wěn)定。在如此混亂的世界中，三體人進(jìn)化出了在“混亂時(shí)代”自我“脫水”并休眠，而在“穩(wěn)定時(shí)代”溫和地復(fù)蘇并生活的能力。
書中對(duì)恒星系統(tǒng)有趣的形象化描述激發(fā)了我們對(duì)萬有引力中n體類問題以及解決這類問題的數(shù)值算法的研究。此文涉及一些理解問題必要的萬有引力核心概念，以及描繪系統(tǒng)解方程式必要的數(shù)值算法的核心概念。
此文將講述以下工具和概念的實(shí)施方法：
· 使用Scipy模型中的odeint函數(shù)解Python中的微分方程
·?無量綱化方程式
· 在Matploylib中進(jìn)行3D繪制
萬有引力概要
1. 牛頓的萬有引力定律
牛頓的萬有引力定律提出，任何兩個(gè)質(zhì)點(diǎn)之間都存在相互吸引力（稱之為萬有引力），其大小與它們質(zhì)量的乘積直接成正比，與它們之間距離的平方成反比。 以下方程式以向量的形式表示了這條定律：
在這里，G為萬有引力常數(shù)，m?和 m?為兩物體的質(zhì)量，r為物體間的距離。單位向量從m?指向 m? 并且力的作用方向也是相同的。
2. 運(yùn)動(dòng)方程
?
根據(jù)牛頓第二運(yùn)動(dòng)定律，物體上的合力造成了物體動(dòng)量的凈變化——簡(jiǎn)而言之，力就是質(zhì)量乘以加速度。因此，將上述方程式應(yīng)用到質(zhì)量m? 的物體中，就可以得到以下運(yùn)動(dòng)微分方程。
這里要注意的是，單位向量r被分解成向量r除以其大小|r|，因此要增加分母中r項(xiàng)的冪到3.
?
現(xiàn)在得到一個(gè)二階微分方程，它描繪了兩個(gè)物體間因重力而存在的相互作用。為簡(jiǎn)化解法，可以將其分解成兩個(gè)一階微分方程。
?
物體的加速度是物體速率隨時(shí)間的變化，因此速率的一階微分可以替代位置的二階微分。類似地，速率可表示為位置的一階微分。
指數(shù)i 表示要計(jì)算位置和速率的物體，而指數(shù)j表示和物體i相互作用的其他物體。因此，對(duì)一個(gè)二體系統(tǒng)來說，就要解兩組方程式。
3. 質(zhì)心
?
另一個(gè)值得記下來的實(shí)用概念是系統(tǒng)的質(zhì)心。質(zhì)心是一個(gè)點(diǎn)，在這個(gè)點(diǎn)上系統(tǒng)的所有質(zhì)量矩總和為零——簡(jiǎn)而言之，可以將其想象成整個(gè)系統(tǒng)質(zhì)量平衡的點(diǎn)。
?
有個(gè)簡(jiǎn)單的公式可以找到系統(tǒng)的質(zhì)心和速率，該公式包括位置和速率向量的質(zhì)量加權(quán)平均值。
在三體系統(tǒng)建模之前，先來為一個(gè)二體系統(tǒng)建模，觀測(cè)其行為然后將代碼延伸應(yīng)用到三體系統(tǒng)中。
二體模型
1. 半人馬座α星恒星系統(tǒng)
一個(gè)著名的二體系統(tǒng)實(shí)例或許就是半人馬座α星恒星系統(tǒng)。它包括三顆恒星——半人馬座α星A，半人馬座α星B，半人馬座α星C（通常稱之為比鄰星）。然而，由于和其他兩顆恒星相比，比鄰星的質(zhì)量小到可以忽視，半人馬座α星也被視作雙星系統(tǒng)。這兒有個(gè)需注意的要的點(diǎn)是，n體系統(tǒng)中考慮到的物體都有相似的質(zhì)量。因此，日-地-月不是一個(gè)三體系統(tǒng)，因?yàn)樗鼈儧]有同等的質(zhì)量，并且地球和月球?qū)μ柕能壽E影響不大。
約翰·科洛西莫在智利的帕洛馬天文臺(tái)捕捉到的半人馬座α星雙星系統(tǒng)
2. 無量綱化
在解方程式之前，需要先對(duì)方程式進(jìn)行無量綱化。這意味著什么？把方程式中（如位置、速率、質(zhì)量等）所有具有維數(shù)（如分別為m, m/s, kg）的量轉(zhuǎn)換成大小接近單位的無因次量。這樣做的原因是：
· 在微分方程中，不同項(xiàng)有不同的數(shù)量級(jí)（從0.1到103?）。如此巨大的差異可能導(dǎo)致數(shù)值算法收斂變得緩慢。
· 如果所有項(xiàng)的大小都十分接近單位，從計(jì)算方面上講所有運(yùn)算價(jià)格都將比數(shù)值大小不平衡時(shí)低廉。
· 將會(huì)得到一個(gè)有關(guān)大小的參考點(diǎn)。例如，如果給一個(gè)4×103? kg的量，可能無法衡量在宇宙尺度上它是大是小。然而，如果說是太陽質(zhì)量的2倍，就很容易把握這個(gè)量的意義。
將每個(gè)量除以一個(gè)固定的參考量以對(duì)方程式進(jìn)行無量綱化。例如，用質(zhì)量項(xiàng)除以太陽的質(zhì)量，半人馬座α星系統(tǒng)中兩星間距離的位置（或距離）項(xiàng)，半人馬座α星軌道周期的時(shí)間項(xiàng)，以及地球繞日公轉(zhuǎn)的相對(duì)速率。
當(dāng)把每一項(xiàng)除以參考量的時(shí)候，也需要將其相乘以免改變方程式。所有這些量和G在一起就能變成一個(gè)常數(shù)，比如方程式1中的K?和方程式2中的 K?。因此，無量綱化的方程式即如下所示：
項(xiàng)上的橫杠表示該項(xiàng)是無因次的。因此這就是要用在模擬中的最終方程式。
3. 代碼
?
從為模擬輸入所要求的模塊開始。
#Import scipy
import scipy as sci
#Import matplotlib and associated modules for 3D andanimations
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import animation
接下來，定義無量綱化方程式中用到的常數(shù)和參考量，以及凈常數(shù)K??和 K?。
#Define universal gravitation constant
G=6.67408e-11 #N-m2/kg2
#Reference quantities
m_nd=1.989e+30 #kg #mass of the sun
r_nd=5.326e+12 #m #distance between stars in Alpha Centauri
v_nd=30000 #m/s #relative velocity of earth around the sun
t_nd=79.91*365*24*3600*0.51 #s #orbital period of Alpha Centauri
#Net constants
K1=G*t_nd*m_nd/(r_nd**2*v_nd)
K2=v_nd*t_nd/r_nd
現(xiàn)在要定義一些參數(shù)，需要模擬兩顆恒星的質(zhì)量、初始位置和初速度。需要注意的是這些參數(shù)是無因次的，所以定義半人馬座α星A的質(zhì)量為1.1（意味著其質(zhì)量為參考量太陽質(zhì)量的1.1倍）。任意定義速率，則沒有物體能夠逃脫彼此的引力。
#Define masses
m1=1.1 #Alpha Centauri A
m2=0.907 #Alpha Centauri B
#Define initial position vectors
r1=[-0.5,0,0] #m
r2=[0.5,0,0] #m
#Convert pos vectors to arrays
r1=sci.array(r1,dtype="float64")
r2=sci.array(r2,dtype="float64")
#Find Centre of Mass
r_com=(m1*r1+m2*r2)/(m1+m2)
#Define initial velocities
v1=[0.01,0.01,0] #m/s
v2=[-0.05,0,-0.1] #m/s
#Convert velocity vectors to arrays
v1=sci.array(v1,dtype="float64")
v2=sci.array(v2,dtype="float64")
#Find velocity of COM
v_com=(m1*v1+m2*v2)/(m1+m2)
現(xiàn)在已經(jīng)定義了大多數(shù)主要的模擬要求的量。可以為解方程組在scipy中準(zhǔn)備odeint 求解器了。
為了解常微分方程，需要有方程式（肯定的！），一組初始條件以及解方程所需的時(shí)間跨度。Odeint求解器也要求具備這三樣基本要素。方程式通過函數(shù)的形式被定義。函數(shù)以其順序接受包含所有因變量（此處為位置和速率）的數(shù)組和包含所有因變量（此處為時(shí)間）的數(shù)組，函數(shù)返回?cái)?shù)組中所有微分的值。
#A function defining the equations of motion
def TwoBodyEquations(w,t,G,m1,m2):
??? r1=w[:3]
??? r2=w[3:6]
??? v1=w[6:9]
??? v2=w[9:12]
? ? r=sci.linalg.norm(r2-r1) #Calculate magnitude or norm of vector
? ? dv1bydt=K1*m2*(r2-r1)/r**3
??? dv2bydt=K1*m1*(r1-r2)/r**3
??? dr1bydt=K2*v1
??? dr2bydt=K2*v2
? ? r_derivs=sci.concatenate((dr1bydt,dr2bydt))
? ? derivs=sci.concatenate((r_derivs,dv1bydt,dv2bydt))
??? return derivs
從這段代碼中也許很容易區(qū)分出微分方程。那其他零碎的東西是什么？記得嗎？現(xiàn)在在解三維方程，所以每一位置和速率向量都會(huì)有3個(gè)分量。考慮一下之前章節(jié)給出的雙矢量微分方程，它們需要解向量的所有3個(gè)分量。因此，需要為一個(gè)物體解6標(biāo)量的微分方程。同理，兩個(gè)物體就要解12標(biāo)量的微分方程。所以要制作一個(gè)大小為12，含兩個(gè)物體的質(zhì)量和速率的數(shù)組w。
在函數(shù)的最后，連結(jié)或加入所有不同的導(dǎo)數(shù)并返回大小為12的derivs。
最難的部分已經(jīng)做完了！剩下的就是把函數(shù)，初始條件和和時(shí)間跨度輸入到odeint函數(shù)中去。
#Package initial parameters
init_params=sci.array([r1,r2,v1,v2]) #create array of initial params
init_params=init_params.flatten() #flatten array to make it 1D
time_span=sci.linspace(0,8,500) #8 orbital periods and 500 points
#Run the ODE solver
import scipy.integrate
two_body_sol=sci.integrate.odeint(TwoBodyEquations,init_params,time_span,args=(G,m1,m2))
變量two_body_sol 包含有關(guān)二體系統(tǒng)所有的信息，包括位置向量和速率向量。為了創(chuàng)建圖和動(dòng)畫，只需要有能延伸到兩個(gè)不同變量的位置向量就可以了。
r1_sol=two_body_sol[:,:3]
r2_sol=two_body_sol[:,3:6]
現(xiàn)在是繪圖！在這要用到Matplotlib的3D標(biāo)圖功能。
#Create figure
fig=plt.figure(figsize=(15,15))
#Create 3D axes
ax=fig.add_subplot(111,projection="3d")
#Plot the orbits
ax.plot(r1_sol[:,0],r1_sol[:,1],r1_sol[:,2],color="darkblue")
ax.plot(r2_sol[:,0],r2_sol[:,1],r2_sol[:,2],color="tab:red")
#Plot the final positions of the stars
ax.scatter(r1_sol[-1,0],r1_sol[-1,1],r1_sol[-1,2],color="darkblue",marker="o",s=100,label="AlphaCentauri A")
ax.scatter(r2_sol[-1,0],r2_sol[-1,1],r2_sol[-1,2],color="tab:red",marker="o",s=100,label="AlphaCentauri B")
#Add a few more bells and whistles
ax.set_xlabel("x-coordinate",fontsize=14)
ax.set_ylabel("y-coordinate",fontsize=14)
ax.set_zlabel("z-coordinate",fontsize=14)
ax.set_title("Visualization of orbits of stars in a two-bodysystem\n",fontsize=14)
ax.legend(loc="upper left",fontsize=14)
最終的圖像清晰表明，軌道遵循可預(yù)測(cè)的模式，正如二體問題算法所期望的那樣。
顯示兩顆星球軌道時(shí)間演變的Matplotlib圖像
這里有一個(gè)展示軌道一步步演變的動(dòng)畫。
Matplotlib中展示軌道一步步演變的動(dòng)畫
還可以再制作一個(gè)可視化圖像，它來自質(zhì)心的參考框架。上面的圖像來自空間中任一平穩(wěn)點(diǎn)，但如果觀察來自質(zhì)心系統(tǒng)的兩個(gè)物體的移動(dòng)，就可以看到一個(gè)更加形象的圖案。
因此，先找到每一時(shí)間步上質(zhì)心的位置，然后從兩個(gè)物體的位置向量上減去其向量以找到它們相對(duì)應(yīng)質(zhì)心的位置。
#Find location of COM
rcom_sol=(m1*r1_sol+m2*r2_sol)/(m1+m2)
#Find location of Alpha Centauri A w.r.t COM
r1com_sol=r1_sol-rcom_sol
#Find location of Alpha Centauri B w.r.t COM
r2com_sol=r2_sol-rcom_sol
最后，可以使用更改變量后的繪制之前圖像使用的代碼，繪制以下圖像。
來自COM的顯示兩顆星球軌道時(shí)間演變的Matplotlib圖像
如果坐在COM觀察兩個(gè)物體，就可以看到以上軌道。從這個(gè)模擬來看它并不清晰，因?yàn)闀r(shí)間尺度非常小，然而即使是這些軌道也一直在輕微地旋轉(zhuǎn)。
現(xiàn)在就很清楚了，兩個(gè)物體嚴(yán)格遵循預(yù)測(cè)的路線，而且可以用函數(shù)——也許是個(gè)橢圓體方程——來描繪它們?cè)诳臻g中如二體系統(tǒng)所期望的移動(dòng)。
三體模型
1.?代碼
現(xiàn)在為了把之前的代碼延伸到三體系統(tǒng)，需要給常數(shù)增加一些東西——增加第三體的質(zhì)量、位置和速率向量。把第三恒星的質(zhì)量視作和太陽的質(zhì)量等同。
#Mass of the Third Star
m3=1.0 #Third Star
#Position of the Third Star
r3=[0,1,0] #m
r3=sci.array(r3,dtype="float64")
#Velocity of the Third Star
v3=[0,-0.01,0]
v3=sci.array(v3,dtype="float64")
?
需要更新代碼中質(zhì)心和質(zhì)心速率的公式。
#Update COM formula
r_com=(m1*r1+m2*r2+m3*r3)/(m1+m2+m3)
#Update velocity of COM formula
v_com=(m1*v1+m2*v2+m3*v3)/(m1+m2+m3)
對(duì)一個(gè)三體系統(tǒng)來說，需要修改運(yùn)動(dòng)方程使之包括另一物體施加的額外引力。因此，需要在RHS上，對(duì)問題中每一對(duì)物體施加力的其他物體增加一個(gè)力項(xiàng)。在三體系統(tǒng)的情況下，一個(gè)物體會(huì)受到其余兩個(gè)物體施加的力的影響并因此在RHS上出現(xiàn)兩個(gè)力項(xiàng)。數(shù)學(xué)上可表示為：
為在代碼中反映這些變化，需要為odeint求解器創(chuàng)建一個(gè)新函數(shù)。
def ThreeBodyEquations(w,t,G,m1,m2,m3):
??? r1=w[:3]
??? r2=w[3:6]
??? r3=w[6:9]
??? v1=w[9:12]
??? v2=w[12:15]
??? v3=w[15:18]
? ? r12=sci.linalg.norm(r2-r1)
??? r13=sci.linalg.norm(r3-r1)
??? r23=sci.linalg.norm(r3-r2)
???
? ? dv1bydt=K1*m2*(r2-r1)/r12**3+K1*m3*(r3-r1)/r13**3
? ? dv2bydt=K1*m1*(r1-r2)/r12**3+K1*m3*(r3-r2)/r23**3
??? dv3bydt=K1*m1*(r1-r3)/r13**3+K1*m2*(r2-r3)/r23**3
??? dr1bydt=K2*v1
??? dr2bydt=K2*v2
??? dr3bydt=K2*v3
? ??r12_derivs=sci.concatenate((dr1bydt,dr2bydt))
? ? r_derivs=sci.concatenate((r12_derivs,dr3bydt))
? ? v12_derivs=sci.concatenate((dv1bydt,dv2bydt))
??? v_derivs=sci.concatenate((v12_derivs,dv3bydt))
? ? derivs=sci.concatenate((r_derivs,v_derivs))
??? return derivs
最后，調(diào)用odeint函數(shù)并向其提供上述函數(shù)連同初始條件。
#Package initial parameters
init_params=sci.array([r1,r2,r3,v1,v2,v3]) #Initial parameters
init_params=init_params.flatten() #Flatten to make 1D array
time_span=sci.linspace(0,20,500) #20 orbital periods and 500 points
#Run the ODE solver
import scipy.integrate
three_body_sol=sci.integrate.odeint(ThreeBodyEquations,init_params,time_span,args=(G,m1,m2,m3))
?
正如二體模擬一樣，需要提取所有三體的位置進(jìn)行繪圖。
r1_sol=three_body_sol[:,:3]
r2_sol=three_body_sol[:,3:6]
r3_sol=three_body_sol[:,6:9]
可以把上述章節(jié)給出的代碼稍作修改后運(yùn)用到最后的繪圖上。 正如下圖呈現(xiàn)的混亂，這個(gè)軌道沒有預(yù)期的圖案。
顯示三顆星球軌道時(shí)間演變的Matplotlib圖像
?
動(dòng)畫會(huì)使這張混亂的圖變得更容易理解。
Matplotlib中展示軌道一步步演變的動(dòng)畫
這里有另一個(gè)初始配置的解法，其中可以觀察到，解法似乎在一開始是穩(wěn)定的，但接著就突然變得不穩(wěn)定了。
Matplotlib中展示軌道一步步演變的動(dòng)畫
可以試著玩玩這些初始條件，看看不同種類的解法。近年來，由于計(jì)算機(jī)能力的增強(qiáng)，人們已經(jīng)發(fā)現(xiàn)了很多關(guān)于三體問題的有趣的解法，其中一些是周期性的——如figure-8解法，在這里三個(gè)物體全都在平面的figure-8路線上運(yùn)動(dòng)。
留言 點(diǎn)贊 發(fā)個(gè)朋友圈
我們一起分享AI學(xué)習(xí)與發(fā)展的干貨
編譯組：鮑歆祎、李林虹
相關(guān)鏈接：
https://towardsdatascience.com/modelling-the-three-body-problem-in-classical-mechanics-using-python-9dc270ad7767
如需轉(zhuǎn)載，請(qǐng)后臺(tái)留言，遵守轉(zhuǎn)載規(guī)范
推薦文章閱讀
ACL2018論文集50篇解讀
EMNLP2017論文集28篇論文解讀
2018年AI三大頂會(huì)中國學(xué)術(shù)成果全鏈接
ACL2017 論文集：34篇解讀干貨全在這里
10篇AAAI2017經(jīng)典論文回顧
長按識(shí)別二維碼可添加關(guān)注
讀芯君愛你
更多文章、技術(shù)交流、商務(wù)合作、聯(lián)系博主
微信掃碼或搜索：z360901061
微信掃一掃加我為好友
QQ號(hào)聯(lián)系： 360901061
您的支持是博主寫作最大的動(dòng)力，如果您喜歡我的文章，感覺我的文章對(duì)您有幫助，請(qǐng)用微信掃描下面二維碼支持博主2元、5元、10元、20元等您想捐的金額吧，狠狠點(diǎn)擊下面給點(diǎn)支持吧，站長非常感激您！手機(jī)微信長按不能支付解決辦法：請(qǐng)將微信支付二維碼保存到相冊(cè)，切換到微信，然后點(diǎn)擊微信右上角掃一掃功能，選擇支付二維碼完成支付。
【本文對(duì)您有幫助就好】元
