nd是什么意思翻译成中文,牛顿是什么人格

2022-07-23 03:11:54 经济时评

　　全文是在7726,写的，预计学习时间是15分钟或更长。　　

　　图片来自flickr，凯文吉尔　　

　　中国作家刘的科幻小说《三体》描绘了一个虚构的外星文明，它存在于一个被三颗恒星包围的“三体”星球上。你能想象这个文明的存在因为三颗星而和我们大相径庭吗？耀眼的阳光？一个持续的夏天？事实证明，情况要糟糕得多。　　

　　生活在只有一颗主要恒星的太阳系中是幸运的，因为这使得这颗恒星(太阳)的轨道是可预测的。即使增加一颗星，系统也会保持稳定。系统有一个解叫解析解，就是描述解方程，得到一个能准确提供系统从一秒到一百万年演化的函数。　　

　　但是，当第三体加入后，就会出现一种特殊的情况。系统将变得混乱和不可预测。没有解析解(少数特定情况除外)，方程只能在计算机上数值求解。它们会突然由稳定变为不稳定，或者由不稳定变为稳定。在这样一个混乱的世界里，三人在混乱的时代里进化出了“脱水”和沉睡的能力，而在稳定的时代里却恢复并温和地生活着。　　

　　书中对恒星系统有趣的可视化描述，启发我们研究万有引力中的N体问题以及求解这类问题的数值算法。本文论述了理解问题所必需的引力的一些核心概念，以及描述系统方程所必需的数值算法的核心概念。　　

　　本文将描述以下工具和概念的实现方法：　　

　　利用Scipy模型中的odeint函数在Python中求解微分方程　　

　　无量纲方程　　

　　Matploylib中的三维渲染　　

　　万有引力概要 　　

　　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 scipyimport 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() #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 包含有关二体系统所有的信息，包括位置向量和速率向量。为了创建图和动画，只需要有能延伸到两个不同变量的位置向量就可以了。

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 formulav_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() #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>

可以把上述章节给出的代码稍作修改后运用到最后的绘图上。正如下图呈现的混乱，这个轨道没有预期的图案。

显示三颗星球轨道时间演变的Matplotlib图像