全文共 7726 字，預計學習時長 15 分鐘或更長
圖片來自flickr, 凱文·吉爾
中國作家劉慈欣的科幻小說《三體》中描繪了存在于被三顆恒星環繞的“三體”星球上的一種虛構外星文明。能想象這種文明的存在因三顆恒星而和我們的文明大不相同嗎？炫目的陽光？持續的夏日？事實證明，情況要糟糕很多。
生活在僅有一顆主要恒星的太陽系是值得慶幸的，因為這使得這顆恒星（太陽）的軌道有可預測性。即使增加一顆恒星，這個系統仍能保持穩定。該系統有個被稱為分析解的解法，即描繪解方程式，并得到可以精確提供該系統從一秒到百萬年時間演變的函數。
然而，當加入第三體時，就會發生特殊情況。系統會變得混亂而不可預測。沒有分析解法（除了少數特定情況），并且只能在計算機上用數字解方程式。它們會突然從穩定變到不穩定，或者從不穩定變到穩定。在如此混亂的世界中，三體人進化出了在“混亂時代”自我“脫水”并休眠，而在“穩定時代”溫和地復蘇并生活的能力。
書中對恒星系統有趣的形象化描述激發了我們對萬有引力中n體類問題以及解決這類問題的數值算法的研究。此文涉及一些理解問題必要的萬有引力核心概念，以及描繪系統解方程式必要的數值算法的核心概念。
此文將講述以下工具和概念的實施方法：
· 使用Scipy模型中的odeint函數解Python中的微分方程
·?無量綱化方程式
· 在Matploylib中進行3D繪制
萬有引力概要
1. 牛頓的萬有引力定律
牛頓的萬有引力定律提出，任何兩個質點之間都存在相互吸引力（稱之為萬有引力），其大小與它們質量的乘積直接成正比，與它們之間距離的平方成反比。 以下方程式以向量的形式表示了這條定律：
在這里，G為萬有引力常數，m?和 m?為兩物體的質量，r為物體間的距離。單位向量從m?指向 m? 并且力的作用方向也是相同的。
2. 運動方程
?
根據牛頓第二運動定律，物體上的合力造成了物體動量的凈變化——簡而言之，力就是質量乘以加速度。因此，將上述方程式應用到質量m? 的物體中，就可以得到以下運動微分方程。
這里要注意的是，單位向量r被分解成向量r除以其大小|r|，因此要增加分母中r項的冪到3.
?
現在得到一個二階微分方程，它描繪了兩個物體間因重力而存在的相互作用。為簡化解法，可以將其分解成兩個一階微分方程。
?
物體的加速度是物體速率隨時間的變化，因此速率的一階微分可以替代位置的二階微分。類似地，速率可表示為位置的一階微分。
指數i 表示要計算位置和速率的物體，而指數j表示和物體i相互作用的其他物體。因此，對一個二體系統來說，就要解兩組方程式。
3. 質心
?
另一個值得記下來的實用概念是系統的質心。質心是一個點，在這個點上系統的所有質量矩總和為零——簡而言之，可以將其想象成整個系統質量平衡的點。
?
有個簡單的公式可以找到系統的質心和速率，該公式包括位置和速率向量的質量加權平均值。
在三體系統建模之前，先來為一個二體系統建模，觀測其行為然后將代碼延伸應用到三體系統中。
二體模型
1. 半人馬座α星恒星系統
一個著名的二體系統實例或許就是半人馬座α星恒星系統。它包括三顆恒星——半人馬座α星A，半人馬座α星B，半人馬座α星C（通常稱之為比鄰星）。然而，由于和其他兩顆恒星相比，比鄰星的質量小到可以忽視，半人馬座α星也被視作雙星系統。這兒有個需注意的要的點是，n體系統中考慮到的物體都有相似的質量。因此，日-地-月不是一個三體系統，因為它們沒有同等的質量，并且地球和月球對太陽的軌跡影響不大。
約翰·科洛西莫在智利的帕洛馬天文臺捕捉到的半人馬座α星雙星系統
2. 無量綱化
在解方程式之前，需要先對方程式進行無量綱化。這意味著什么？把方程式中（如位置、速率、質量等）所有具有維數（如分別為m, m/s, kg）的量轉換成大小接近單位的無因次量。這樣做的原因是：
· 在微分方程中，不同項有不同的數量級（從0.1到103?）。如此巨大的差異可能導致數值算法收斂變得緩慢。
· 如果所有項的大小都十分接近單位，從計算方面上講所有運算價格都將比數值大小不平衡時低廉。
· 將會得到一個有關大小的參考點。例如，如果給一個4×103? kg的量，可能無法衡量在宇宙尺度上它是大是小。然而，如果說是太陽質量的2倍，就很容易把握這個量的意義。
將每個量除以一個固定的參考量以對方程式進行無量綱化。例如，用質量項除以太陽的質量，半人馬座α星系統中兩星間距離的位置（或距離）項，半人馬座α星軌道周期的時間項，以及地球繞日公轉的相對速率。
當把每一項除以參考量的時候，也需要將其相乘以免改變方程式。所有這些量和G在一起就能變成一個常數，比如方程式1中的K?和方程式2中的 K?。因此，無量綱化的方程式即如下所示：
項上的橫杠表示該項是無因次的。因此這就是要用在模擬中的最終方程式。
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
接下來，定義無量綱化方程式中用到的常數和參考量，以及凈常數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
現在要定義一些參數，需要模擬兩顆恒星的質量、初始位置和初速度。需要注意的是這些參數是無因次的，所以定義半人馬座α星A的質量為1.1（意味著其質量為參考量太陽質量的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)
現在已經定義了大多數主要的模擬要求的量。可以為解方程組在scipy中準備odeint 求解器了。
為了解常微分方程，需要有方程式（肯定的！），一組初始條件以及解方程所需的時間跨度。Odeint求解器也要求具備這三樣基本要素。方程式通過函數的形式被定義。函數以其順序接受包含所有因變量（此處為位置和速率）的數組和包含所有因變量（此處為時間）的數組，函數返回數組中所有微分的值。
#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
從這段代碼中也許很容易區分出微分方程。那其他零碎的東西是什么？記得嗎？現在在解三維方程，所以每一位置和速率向量都會有3個分量。考慮一下之前章節給出的雙矢量微分方程，它們需要解向量的所有3個分量。因此，需要為一個物體解6標量的微分方程。同理，兩個物體就要解12標量的微分方程。所以要制作一個大小為12，含兩個物體的質量和速率的數組w。
在函數的最后，連結或加入所有不同的導數并返回大小為12的derivs。
最難的部分已經做完了！剩下的就是把函數，初始條件和和時間跨度輸入到odeint函數中去。
#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 包含有關二體系統所有的信息，包括位置向量和速率向量。為了創建圖和動畫，只需要有能延伸到兩個不同變量的位置向量就可以了。
r1_sol=two_body_sol[:,:3]
r2_sol=two_body_sol[:,3:6]
現在是繪圖！在這要用到Matplotlib的3D標圖功能。
#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)
最終的圖像清晰表明，軌道遵循可預測的模式，正如二體問題算法所期望的那樣。
顯示兩顆星球軌道時間演變的Matplotlib圖像
這里有一個展示軌道一步步演變的動畫。
Matplotlib中展示軌道一步步演變的動畫
還可以再制作一個可視化圖像，它來自質心的參考框架。上面的圖像來自空間中任一平穩點，但如果觀察來自質心系統的兩個物體的移動，就可以看到一個更加形象的圖案。
因此，先找到每一時間步上質心的位置，然后從兩個物體的位置向量上減去其向量以找到它們相對應質心的位置。
#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的顯示兩顆星球軌道時間演變的Matplotlib圖像
如果坐在COM觀察兩個物體，就可以看到以上軌道。從這個模擬來看它并不清晰，因為時間尺度非常小，然而即使是這些軌道也一直在輕微地旋轉。
現在就很清楚了，兩個物體嚴格遵循預測的路線，而且可以用函數——也許是個橢圓體方程——來描繪它們在空間中如二體系統所期望的移動。
三體模型
1.?代碼
現在為了把之前的代碼延伸到三體系統，需要給常數增加一些東西——增加第三體的質量、位置和速率向量。把第三恒星的質量視作和太陽的質量等同。
#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")
?
需要更新代碼中質心和質心速率的公式。
#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)
對一個三體系統來說，需要修改運動方程使之包括另一物體施加的額外引力。因此，需要在RHS上，對問題中每一對物體施加力的其他物體增加一個力項。在三體系統的情況下，一個物體會受到其余兩個物體施加的力的影響并因此在RHS上出現兩個力項。數學上可表示為：
為在代碼中反映這些變化，需要為odeint求解器創建一個新函數。
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
最后，調用odeint函數并向其提供上述函數連同初始條件。
#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))
?
正如二體模擬一樣，需要提取所有三體的位置進行繪圖。
r1_sol=three_body_sol[:,:3]
r2_sol=three_body_sol[:,3:6]
r3_sol=three_body_sol[:,6:9]
可以把上述章節給出的代碼稍作修改后運用到最后的繪圖上。 正如下圖呈現的混亂，這個軌道沒有預期的圖案。
顯示三顆星球軌道時間演變的Matplotlib圖像
?
動畫會使這張混亂的圖變得更容易理解。
Matplotlib中展示軌道一步步演變的動畫
這里有另一個初始配置的解法，其中可以觀察到，解法似乎在一開始是穩定的，但接著就突然變得不穩定了。
Matplotlib中展示軌道一步步演變的動畫
可以試著玩玩這些初始條件，看看不同種類的解法。近年來，由于計算機能力的增強，人們已經發現了很多關于三體問題的有趣的解法，其中一些是周期性的——如figure-8解法，在這里三個物體全都在平面的figure-8路線上運動。
留言 點贊 發個朋友圈
我們一起分享AI學習與發展的干貨
編譯組：鮑歆祎、李林虹
相關鏈接：
https://towardsdatascience.com/modelling-the-three-body-problem-in-classical-mechanics-using-python-9dc270ad7767
如需轉載，請后臺留言，遵守轉載規范
推薦文章閱讀
ACL2018論文集50篇解讀
EMNLP2017論文集28篇論文解讀
2018年AI三大頂會中國學術成果全鏈接
ACL2017 論文集：34篇解讀干貨全在這里
10篇AAAI2017經典論文回顧
長按識別二維碼可添加關注
讀芯君愛你
更多文章、技術交流、商務合作、聯系博主
微信掃碼或搜索：z360901061
微信掃一掃加我為好友
QQ號聯系： 360901061
您的支持是博主寫作最大的動力，如果您喜歡我的文章，感覺我的文章對您有幫助，請用微信掃描下面二維碼支持博主2元、5元、10元、20元等您想捐的金額吧，狠狠點擊下面給點支持吧，站長非常感激您！手機微信長按不能支付解決辦法：請將微信支付二維碼保存到相冊，切換到微信，然后點擊微信右上角掃一掃功能，選擇支付二維碼完成支付。
【本文對您有幫助就好】元
