|
圖片來自flickr, 凱文·吉爾 中國作家劉慈欣的科幻小說《三體》中描繪了存在于被三顆恒星環(huán)繞的“三體”星球上的一種虛構外星文明。能想象這種文明的存在因三顆恒星而和我們的文明大不相同嗎,？炫目的陽光？持續(xù)的夏日？事實證明,，情況要糟糕很多,。 生活在僅有一顆主要恒星的太陽系是值得慶幸的，因為這使得這顆恒星（太陽）的軌道有可預測性,。即使增加一顆恒星,，這個系統(tǒng)仍能保持穩(wěn)定。該系統(tǒng)有個被稱為分析解的解法,，即描繪解方程式,，并得到可以精確提供該系統(tǒng)從一秒到百萬年時間演變的函數(shù)。 然而,，當加入第三體時,，就會發(fā)生特殊情況。系統(tǒng)會變得混亂而不可預測,。沒有分析解法（除了少數(shù)特定情況）,，并且只能在計算機上用數(shù)字解方程式。它們會突然從穩(wěn)定變到不穩(wěn)定,，或者從不穩(wěn)定變到穩(wěn)定,。在如此混亂的世界中，三體人進化出了在“混亂時代”自我“脫水”并休眠,，而在“穩(wěn)定時代”溫和地復蘇并生活的能力,。 書中對恒星系統(tǒng)有趣的形象化描述激發(fā)了我們對萬有引力中n體類問題以及解決這類問題的數(shù)值算法的研究。此文涉及一些理解問題必要的萬有引力核心概念,，以及描繪系統(tǒng)解方程式必要的數(shù)值算法的核心概念,。 此文將講述以下工具和概念的實施方法： · 使用Scipy模型中的odeint函數(shù)解Python中的微分方程 · 無量綱化方程式 · 在Matploylib中進行3D繪制 萬有引力概要 1. 牛頓的萬有引力定律 牛頓的萬有引力定律提出,，任何兩個質(zhì)點之間都存在相互吸引力（稱之為萬有引力）,，其大小與它們質(zhì)量的乘積直接成正比，與它們之間距離的平方成反比,。以下方程式以向量的形式表示了這條定律： 在這里，G為萬有引力常數(shù),，m?和 m?為兩物體的質(zhì)量,，r為物體間的距離。單位向量從m?指向 m? 并且力的作用方向也是相同的,。 2. 運動方程 根據(jù)牛頓第二運動定律,，物體上的合力造成了物體動量的凈變化——簡而言之，力就是質(zhì)量乘以加速度,。因此,，將上述方程式應用到質(zhì)量m? 的物體中，就可以得到以下運動微分方程。 這里要注意的是,，單位向量r被分解成向量r除以其大小|r|，因此要增加分母中r項的冪到3. 現(xiàn)在得到一個二階微分方程,，它描繪了兩個物體間因重力而存在的相互作用,。為簡化解法，可以將其分解成兩個一階微分方程,。 物體的加速度是物體速率隨時間的變化,，因此速率的一階微分可以替代位置的二階微分。類似地,，速率可表示為位置的一階微分,。 指數(shù)i 表示要計算位置和速率的物體,，而指數(shù)j表示和物體i相互作用的其他物體,。因此，對一個二體系統(tǒng)來說,，就要解兩組方程式,。 3. 質(zhì)心 另一個值得記下來的實用概念是系統(tǒng)的質(zhì)心。質(zhì)心是一個點,，在這個點上系統(tǒng)的所有質(zhì)量矩總和為零——簡而言之,，可以將其想象成整個系統(tǒng)質(zhì)量平衡的點。 有個簡單的公式可以找到系統(tǒng)的質(zhì)心和速率,，該公式包括位置和速率向量的質(zhì)量加權平均值,。 在三體系統(tǒng)建模之前,，先來為一個二體系統(tǒng)建模,，觀測其行為然后將代碼延伸應用到三體系統(tǒng)中。 二體模型 1. 半人馬座α星恒星系統(tǒng) 一個著名的二體系統(tǒng)實例或許就是半人馬座α星恒星系統(tǒng),。它包括三顆恒星——半人馬座α星A，半人馬座α星B,，半人馬座α星C（通常稱之為比鄰星）,。然而，由于和其他兩顆恒星相比,，比鄰星的質(zhì)量小到可以忽視,，半人馬座α星也被視作雙星系統(tǒng)。這兒有個需注意的要的點是,，n體系統(tǒng)中考慮到的物體都有相似的質(zhì)量,。因此,，日-地-月不是一個三體系統(tǒng),，因為它們沒有同等的質(zhì)量,，并且地球和月球對太陽的軌跡影響不大。 約翰·科洛西莫在智利的帕洛馬天文臺捕捉到的半人馬座α星雙星系統(tǒng) 2. 無量綱化 在解方程式之前,，需要先對方程式進行無量綱化。這意味著什么,？把方程式中（如位置,、速率、質(zhì)量等）所有具有維數(shù)（如分別為m, m/s, kg）的量轉換成大小接近單位的無因次量,。這樣做的原因是： · 在微分方程中,，不同項有不同的數(shù)量級（從0.1到103?）。如此巨大的差異可能導致數(shù)值算法收斂變得緩慢,。 · 如果所有項的大小都十分接近單位,，從計算方面上講所有運算價格都將比數(shù)值大小不平衡時低廉。 · 將會得到一個有關大小的參考點,。例如,，如果給一個4×103? kg的量，可能無法衡量在宇宙尺度上它是大是小,。然而,，如果說是太陽質(zhì)量的2倍，就很容易把握這個量的意義,。 將每個量除以一個固定的參考量以對方程式進行無量綱化,。例如，用質(zhì)量項除以太陽的質(zhì)量,，半人馬座α星系統(tǒng)中兩星間距離的位置（或距離）項,，半人馬座α星軌道周期的時間項，以及地球繞日公轉的相對速率,。 當把每一項除以參考量的時候,，也需要將其相乘以免改變方程式。所有這些量和G在一起就能變成一個常數(shù),，比如方程式1中的K?和方程式2中的 K?,。因此，無量綱化的方程式即如下所示： 項上的橫杠表示該項是無因次的,。因此這就是要用在模擬中的最終方程式。 3. 代碼 從為模擬輸入所要求的模塊開始,。 #Import scipyimport scipy as sci#Import matplotlib and associated modules for 3D andanimationsimport matplotlib.pyplot as pltfrom mpl_toolkits.mplot3d import Axes3Dfrom matplotlib import animation 接下來,，定義無量綱化方程式中用到的常數(shù)和參考量，以及凈常數(shù)K? 和 K?。 #Define universal gravitation constantG=6.67408e-11 #N-m2/kg2#Reference quantitiesm_nd=1.989e+30 #kg #mass of the sunr_nd=5.326e+12 #m #distance between stars in Alpha Centauriv_nd=30000 #m/s #relative velocity of earth around the sunt_nd=79.91*365*24*3600*0.51 #s #orbital period of Alpha Centauri#Net constantsK1=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 massesm1=1.1 #Alpha Centauri Am2=0.907 #Alpha Centauri B#Define initial position vectorsr1=[-0.5,0,0] #mr2=[0.5,0,0] #m#Convert pos vectors to arraysr1=sci.array(r1,dtype='float64')r2=sci.array(r2,dtype='float64')#Find Centre of Massr_com=(m1*r1+m2*r2)/(m1+m2)#Define initial velocitiesv1=[0.01,0.01,0] #m/sv2=[-0.05,0,-0.1] #m/s#Convert velocity vectors to arraysv1=sci.array(v1,dtype='float64')v2=sci.array(v2,dtype='float64')#Find velocity of COMv_com=(m1*v1+m2*v2)/(m1+m2) 現(xiàn)在已經(jīng)定義了大多數(shù)主要的模擬要求的量,。可以為解方程組在scipy中準備odeint 求解器了,。 為了解常微分方程,，需要有方程式（肯定的！）,，一組初始條件以及解方程所需的時間跨度,。Odeint求解器也要求具備這三樣基本要素。方程式通過函數(shù)的形式被定義,。函數(shù)以其順序接受包含所有因變量（此處為位置和速率）的數(shù)組和包含所有因變量（此處為時間）的數(shù)組,，函數(shù)返回數(shù)組中所有微分的值。 #A function defining the equations of motiondef 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)在在解三維方程,，所以每一位置和速率向量都會有3個分量,。考慮一下之前章節(jié)給出的雙矢量微分方程,，它們需要解向量的所有3個分量,。因此，需要為一個物體解6標量的微分方程,。同理,，兩個物體就要解12標量的微分方程。所以要制作一個大小為12,，含兩個物體的質(zhì)量和速率的數(shù)組w,。 在函數(shù)的最后，連結或加入所有不同的導數(shù)并返回大小為12的derivs,。 最難的部分已經(jīng)做完了,！剩下的就是把函數(shù)，初始條件和和時間跨度輸入到odeint函數(shù)中去,。 #Package initial parametersinit_params=sci.array([r1,r2,v1,v2]) #create array of initial paramsinit_params=init_params.flatten() #flatten array to make it 1Dtime_span=sci.linspace(0,8,500) #8 orbital periods and 500 points#Run the ODE solverimport scipy.integratetwo_body_sol=sci.integrate.odeint(TwoBodyEquations,init_params,time_span,args=(G,m1,m2)) 變量two_body_sol 包含有關二體系統(tǒng)所有的信息,，包括位置向量和速率向量,。為了創(chuàng)建圖和動畫，只需要有能延伸到兩個不同變量的位置向量就可以了,。 r1_sol=two_body_sol[:,:3]r2_sol=two_body_sol[:,3:6] 現(xiàn)在是繪圖,！在這要用到Matplotlib的3D標圖功能。 #Create figurefig=plt.figure(figsize=(15,15))#Create 3D axesax=fig.add_subplot(111,projection='3d')#Plot the orbitsax.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 starsax.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 whistlesax.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中展示軌道一步步演變的動畫 還可以再制作一個可視化圖像,，它來自質(zhì)心的參考框架,。上面的圖像來自空間中任一平穩(wěn)點，但如果觀察來自質(zhì)心系統(tǒng)的兩個物體的移動,，就可以看到一個更加形象的圖案,。 因此，先找到每一時間步上質(zhì)心的位置,，然后從兩個物體的位置向量上減去其向量以找到它們相對應質(zhì)心的位置,。 #Find location of COMrcom_sol=(m1*r1_sol+m2*r2_sol)/(m1+m2)#Find location of Alpha Centauri A w.r.t COMr1com_sol=r1_sol-rcom_sol#Find location of Alpha Centauri B w.r.t COMr2com_sol=r2_sol-rcom_sol 最后，可以使用更改變量后的繪制之前圖像使用的代碼,，繪制以下圖像,。
來自COM的顯示兩顆星球軌道時間演變的Matplotlib圖像 如果坐在COM觀察兩個物體,，就可以看到以上軌道,。從這個模擬來看它并不清晰，因為時間尺度非常小,，然而即使是這些軌道也一直在輕微地旋轉,。 現(xiàn)在就很清楚了，兩個物體嚴格遵循預測的路線,，而且可以用函數(shù)——也許是個橢圓體方程——來描繪它們在空間中如二體系統(tǒng)所期望的移動,。
三體模型 1. 代碼 現(xiàn)在為了把之前的代碼延伸到三體系統(tǒng),，需要給常數(shù)增加一些東西——增加第三體的質(zhì)量,、位置和速率向量。把第三恒星的質(zhì)量視作和太陽的質(zhì)量等同,。 #Mass of the Third Starm3=1.0 #Third Star#Position of the Third Starr3=[0,1,0] #mr3=sci.array(r3,dtype='float64')#Velocity of the Third Starv3=[0,-0.01,0]v3=sci.array(v3,dtype='float64') 需要更新代碼中質(zhì)心和質(zhì)心速率的公式,。 #Update COM formular_com=(m1*r1+m2*r2+m3*r3)/(m1+m2+m3)#Update velocity of COM formulav_com=(m1*v1+m2*v2+m3*v3)/(m1+m2+m3) 對一個三體系統(tǒng)來說，需要修改運動方程使之包括另一物體施加的額外引力,。因此,，需要在RHS上,，對問題中每一對物體施加力的其他物體增加一個力項。在三體系統(tǒng)的情況下,，一個物體會受到其余兩個物體施加的力的影響并因此在RHS上出現(xiàn)兩個力項,。數(shù)學上可表示為：
為在代碼中反映這些變化,，需要為odeint求解器創(chuàng)建一個新函數(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 parametersinit_params=sci.array([r1,r2,r3,v1,v2,v3]) #Initial parametersinit_params=init_params.flatten() #Flatten to make 1D arraytime_span=sci.linspace(0,20,500) #20 orbital periods and 500 points#Run the ODE solverimport scipy.integratethree_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] 可以把上述章節(jié)給出的代碼稍作修改后運用到最后的繪圖上,。正如下圖呈現(xiàn)的混亂,，這個軌道沒有預期的圖案。
顯示三顆星球軌道時間演變的Matplotlib圖像 動畫會使這張混亂的圖變得更容易理解,。 Matplotlib中展示軌道一步步演變的動畫 這里有另一個初始配置的解法，其中可以觀察到,，解法似乎在一開始是穩(wěn)定的,，但接著就突然變得不穩(wěn)定了。 Matplotlib中展示軌道一步步演變的動畫 可以試著玩玩這些初始條件,，看看不同種類的解法,。近年來，由于計算機能力的增強,，人們已經(jīng)發(fā)現(xiàn)了很多關于三體問題的有趣的解法,，其中一些是周期性的——如figure-8解法，在這里三個物體全都在平面的figure-8路線上運動,。
|
|
|
