全文共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到10³⁰)。如此巨大的差异可能导致数值算法收敛变得缓慢。
· 如果所有项的大小都十分接近单位，从计算方面上讲所有运算价格都将比数值大小不平衡时低廉。
· 将会得到一个有关大小的参考点。例如，如果给一个4×10³⁰ kg的量，可能无法衡量在宇宙尺度上它是大是小。然而，如果说是太阳质量的2倍，就很容易把握这个量的意义。
将每个量除以一个固定的参考量以对方程式进行无量纲化。例如，用质量项除以太阳的质量，半人马座α星系统中两星间距离的位置(或距离)项，半人马座α星轨道周期的时间项，以及地球绕日公转的相对速率。
当把每一项除以参考量的时候，也需要将其相乘以免改变方程式。所有这些量和G在一起就能变成一个常数，比如方程式1中的K₁和方程式2中的 K₂。因此，无量纲化的方程式即如下所示：
项上的横杠表示该项是无因次的。因此这就是要用在模拟中的最终方程式。
3. 代码
从为模拟输入所要求的模块开始。
#Import scipy
import scipy as sci
#Import matplotlib and associated modules for 3D andanimationsimport matplotlib.pyplot as pltfrom mpl_toolkits.mplot3d import Axes3Dfrom matplotlib import animation
接下来，定义无量纲化方程式中用到的常数和参考量，以及净常数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
现在要定义一些参数，需要模拟两颗恒星的质量、初始位置和初速度。需要注意的是这些参数是无因次的，所以定义半人马座α星A的质量为1.1(意味着其质量为参考量太阳质量的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)
现在已经定义了大多数主要的模拟要求的量。可以为解方程组在scipy中准备odeint 求解器了。
为了解常微分方程，需要有方程式(肯定的！)，一组初始条件以及解方程所需的时间跨度。Odeint求解器也要求具备这三样基本要素。方程式通过函数的形式被定义。函数以其顺序接受包含所有因变量(此处为位置和速率)的数组和包含所有因变量(此处为时间)的数组，函数返回数组中所有微分的值。
#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
从这段代码中也许很容易区分出微分方程。那其他零碎的东西是什么？记得吗？现在在解三维方程，所以每一位置和速率向量都会有3个分量。考虑一下之前章节给出的双矢量微分方程，它们需要解向量的所有3个分量。因此，需要为一个物体解6标量的微分方程。同理，两个物体就要解12标量的微分方程。所以要制作一个大小为12，含两个物体的质量和速率的数组w。
在函数的最后，连结或加入所有不同的导数并返回大小为12的derivs。
最难的部分已经做完了！剩下的就是把函数，初始条件和和时间跨度输入到odeint函数中去。
#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.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 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中展示轨道一步步演变的动画
还可以再制作一个可视化图像，它来自质心的参考框架。上面的图像来自空间中任一平稳点，但如果观察来自质心系统的两个物体的移动，就可以看到一个更加形象的图案。
因此，先找到每一时间步上质心的位置，然后从两个物体的位置向量上减去其向量以找到它们相对应质心的位置。
#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观察两个物体，就可以看到以上轨道。从这个模拟来看它并不清晰，因为时间尺度非常小，然而即使是这些轨道也一直在轻微地旋转。
现在就很清楚了，两个物体严格遵循预测的路线，而且可以用函数——也许是个椭圆体方程——来描绘它们在空间中如二体系统所期望的移动。
三体模型
1. 代码
现在为了把之前的代码延伸到三体系统，需要给常数增加一些东西——增加第三体的质量、位置和速率向量。把第三恒星的质量视作和太阳的质量等同。
#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")
需要更新代码中质心和质心速率的公式。
#Update COM formular_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 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.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学习与发展的干货
欢迎关注全平台AI垂类自媒体 “读芯术”
771
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
