用 Python 模拟一个微型太阳系

使用 Matplotlib 中的真实质量、距离和速度数据模拟一个包含太阳、地球、火星和未知彗星的微型太阳系——我对 PCA(主成分分析)感到惊讶,并考虑使用动画来展示机器学习过程.当我设法为图表制作动画时,我很想为一些很酷的东西制作动画。我想到的第一个很酷的东西是太阳系。这篇文章也是……

用 Python 模拟一个微型太阳系

使用 Matplotlib 中的真实质量、距离和速度数据模拟具有太阳、地球、火星和未知彗星的微型太阳系

用 Python 模拟一个微型太阳系

我对 PCA(主成分分析)感到惊讶,并考虑使用动画来展示机器学习过程。当我设法为图表制作动画时,我很想为一些很酷的东西制作动画。我想到的第一个很酷的东西是太阳系。这篇文章也是我回忆牛顿物理、向量空间、matplotlib知识的一次旅程。

为了模拟行星的轨道并为其设置动画,我需要解决两个关键问题:

  1. 生成轨道数据。
  2. 使用 Matplotlib 对数据进行动画处理。

在本文中,我将讨论上述两个问题。如果你赶时间不想看我啰嗦的八卦,你可以在 Github 上查看代码宿主。[0]

作者的所有图像和GIF。

牛顿万有引力定律

根据牛顿万有引力定律,我们可以计算两个物体之间的力。只考虑地球和太阳。太阳对地球的重力影响可以用下式表示:

用 Python 模拟一个微型太阳系

M 表示太阳的质量,m 表示地球的质量,r 表示两个物体中心之间的距离。 G 是一个常数值。

根据牛顿第二定律:

用 Python 模拟一个微型太阳系

我们可以得到加速度a:

用 Python 模拟一个微型太阳系

嗯,加速度a只与太阳的质量和日地距离r有关。现在有了一个已知的。我们可以计算增量时间后地球的速度 – dt。

用 Python 模拟一个微型太阳系

最后,得到位移的变化——dd:

用 Python 模拟一个微型太阳系

(我知道在一篇文章中使用/阅读这么多方程很烦人,而这可能是解释引力轨道模拟基础的最简单、最简洁的方法)

向量空间中的万有引力定律

由于我要使用 matplotlib 绘制地球的轨道,我需要计算 x,y 的值。为简单起见,我将使用 2D 坐标轴而不是 3D。

用 Python 模拟一个微型太阳系

我在上图中绘制了方程的矢量版本。如果这是您第一次阅读方程式,您可能会想到几个问题。

  1. 为什么引力方程中有减号(-)?因为力的方向与距离矢量r相反。当您计算从地球施加在太阳上的力时,您需要删除减号 (-)。在矢量坐标中,我们需要非常小心方向。
  2. 离心力呢,为什么不处理离心力呢?因为离心力实际上是万有引力的结果。如果我们在地球上进行力分析。只有一个引力作用在它上面。当地球围绕太阳旋转时,地球正在向太阳下落。然而,地球的初速度垂直于 G 力将其拉离太阳,这就是地球可以绕太阳公转的原因。

在 Python 代码中,假设太阳在位置 (0,0),地球在位置 (6,6):

sx,sy = 0,0
ex,ey = 6,6

距离 r 将是:

rx = ex - sx
ry = ey - sy

定义 GMm,即日地重力常数为 gravconst_e:

Ms          = 2.0e30                    # mass of sun in kg unit
Me = 5.972e24 # mass of earth in kg unit
G = 6.67e-11
gravconst = G*Ms*Me

模拟器将每 dt 时间计算一个数据点。我们可以将 dt 设置为一天:

daysec = 24.0*60*60
dt = 1*daysec

r (|r|³) 的立方是:

modr3_e = (rx**2 + ry**2)**1.5

要计算施加在地球方向上的力:

# get g force in x and y direction
fx_e = gravconst_e*rx/modr3_e
fy_e = gravconst_e*ry/modr3_e

基于牛顿第二定律的向量形式:

用 Python 模拟一个微型太阳系

结合上述方程:

用 Python 模拟一个微型太阳系

由于我已经计算了力向量,我可以在 dt 时间后得到新的速度向量:

xve += fx_e*dt/Me
yve += fy_e*dt/Me

最后,得到地球的新位移(速度乘以时间就是距离:d=v*t):

xe += xve*dt
ye += yve*dt

这是模拟地球轨道一年的完整代码。

G           = 6.67e-11                  # constant G
Ms = 2.0e30 # sun
Me = 5.972e24 # earth
AU = 1.5e11 # earth sun distance
daysec = 24.0*60*60 # seconds of a day
e_ap_v = 29290 # earth velocity at aphelion
gravconst_e = G*Me*Ms
# setup the starting conditions
# earth
xe,ye,ze = 1.0167*AU,0,0
xve,yve,zve = 0,e_ap_v,0
# sun
xs,ys,zs = 0,0,0
xvs,yvs,zvs = 0,0,0
t = 0.0
dt = 1*daysec # every frame move this time
xelist,yelist,zelist = [],[],[]
xslist,yslist,zslist = [],[],[]
# start simulation
while t<1*365*daysec:
################ earth #############
# compute G force on earth
rx,ry,rz = xe - xs, ye - ys, ze - zs
modr3_e = (rx**2+ry**2+rz**2)**1.5
fx_e = -gravconst_e*rx/modr3_e
fy_e = -gravconst_e*ry/modr3_e
fz_e = -gravconst_e*rz/modr3_e

# update quantities how is this calculated? F = ma -> a = F/m
xve += fx_e*dt/Me
yve += fy_e*dt/Me
zve += fz_e*dt/Me

# update position
xe += xve*dt
ye += yve*dt
ze += zve*dt

# save the position in list
xelist.append(xe)
yelist.append(ye)
zelist.append(ze)

################ the sun ###########
# update quantities how is this calculated? F = ma -> a = F/m
xvs += -fx_e*dt/Ms
yvs += -fy_e*dt/Ms
zvs += -fz_e*dt/Ms

# update position
xs += xvs*dt
ys += yvs*dt
zs += zvs*dt
xslist.append(xs)
yslist.append(ys)
zslist.append(zs)

# update dt
t +=dt

要查看轨道的样子:

import matplotlib.pyplot as plt
plt.plot(xelist,yelist,'-g',lw=2)
plt.axis('equal')
plt.show()

使用 matplotlib 为轨道设置动画

现在,我将把静止的地球轨道变成一个移动的物体,并在平面上保留一条尾迹。以下是我在构建动画过程中解决的一些问题:

  1. 所有的魔法都发生在更新函数内部。基本上,该函数将由 FuncAnimation 每一帧调用。您可以调用 ax.plot(…) 在平面中绘制新像素(所有像素将保留在下一帧中)或调用 set_data 或 set_position 来更新可视对象的数据(matplotlib 将根据数据重绘所有内容对于当前帧)。
  2. 表达式“line_e, = ax.plot(…)”中的逗号一开始看起来很奇怪。但请注意,构成元组的是逗号。 “ax.plot(…)” 返回一个元组对象。这就是为什么您需要在变量名之后添加一个逗号。
  3. 在地球之后绘制追踪尾随。我初始化空列表(exdata,eydata = [],[])以保存所有先前的跟踪点。每次调用更新函数时,将跟踪点附加到列表中。

我不打算把这部分变成动画教程。如果您以前从未尝试过 matplotlib 动画。 Jake VanderPlas 的 Matplotlib 动画教程非常适合阅读。[0]

这是动画轨道的代码:

import matplotlib.pyplot as plt
from matplotlib import animation

fig, ax = plt.subplots(figsize=(10,10))
ax.set_aspect('equal')
ax.grid()

line_e, = ax.plot([],[],'-g',lw=1,c='blue')
point_e, = ax.plot([AU], [0], marker="o"
, markersize=4
, markeredgecolor="blue"
, markerfacecolor="blue")
text_e = ax.text(AU,0,'Earth')

point_s, = ax.plot([0], [0], marker="o"
, markersize=7
, markeredgecolor="yellow"
, markerfacecolor="yellow")
text_s = ax.text(0,0,'Sun')

exdata,eydata = [],[] # earth track
sxdata,sydata = [],[] # sun track

print(len(xelist))

def update(i):
exdata.append(xelist[i])
eydata.append(yelist[i])

line_e.set_data(exdata,eydata)
point_e.set_data(xelist[i],yelist[i])
text_e.set_position((xelist[i],yelist[i]))

point_s.set_data(xslist[i],yslist[i])
text_s.set_position((xslist[i],yslist[i]))
ax.axis('equal')
ax.set_xlim(-3*AU,3*AU)
ax.set_ylim(-3*AU,3*AU)

return line_e,point_s,point_e,text_e,text_s

anim = animation.FuncAnimation(fig
,func=update
,frames=len(xelist)
,interval=1
,blit=True)
plt.show()

Result:

用 Python 模拟一个微型太阳系

在此处查看完整代码。[0]

What if…

改变初始地球速度,比如将地球速度降低到 19,290 m/s(从 29,290 m/s):

用 Python 模拟一个微型太阳系

改变地球的质量,将其增加到 5.972e29 kg(从 5.972e24 kg)。太阳在动!

用 Python 模拟一个微型太阳系

将火星和一颗未知彗星加入轨道系统

自从我成功模拟了地球和太阳的轨道。将更多对象添加到系统中并不难。例如,火星和一颗未知的彗星。

用 Python 模拟一个微型太阳系

在这里查看它的代码。[0]

Matplotlib3D 中的动画

通过对绘图代码进行一些调整,我什至可以在 3D 中为轨道设置动画!在上面的模拟代码中,我已经计算好了z轴数据。

用 Python 模拟一个微型太阳系

Check out the code here.[0]

Summary

在这篇文章中,我介绍了牛顿万有引力定律的基础,如何将方程转换为向量空间,并逐步计算模拟数据。手头有仿真数据,我也讲解了动画制作路上的三个关键问题。

感谢您阅读这里。通过对模拟过程的透彻了解,您甚至可以构建著名的三体问题模拟器。 (是的,我也是同名小说的粉丝)[0]

用 Python 模拟一个微型太阳系

如果您有任何问题,请告诉我,我会尽力回答。并且毫不犹豫地指出我在文章中犯的任何错误。期待与您的进一步联系。

让我停在这里,祝你享受它并快乐编程。

Reference Links

  • 牛顿万有引力定律[0]
  • Earth’s orbit[0]
  • 使用 Tex 渲染数学方程[0]
  • MV3.03 两体问题[0]
  • MV3.04 两体:Python[0]

附录——代码生成矢量坐标图

在运行代码之前,你需要在你的机器上安装 LaTex。

import numpy as np
import matplotlib.pyplot as plt
plt.rcParams['text.usetex'] = True
font = {'family' : 'normal',
'weight' : 'bold',
'size' : 18}
plt.rc('font',**font)

# initialize the stage
fig,ax = plt.subplots(figsize=(8,8))

# set x, and y axis,and remove top and right
ax.spines[['top','right']].set_visible(False)
ax.spines[['bottom','left']].set_position(('data',0.0))
ax.axis('equal')
ax.axes.xaxis.set_visible(False)
ax.axes.yaxis.set_visible(False)
ax.set_xlim(-3,10)
ax.set_ylim(-3,10)

# draw arrows
ax.plot(1, 0, ">k", transform=ax.get_yaxis_transform(), clip_on=False)
ax.plot(0, 1, "^k", transform=ax.get_xaxis_transform(), clip_on=False)

# sun
a = np.linspace(0,2*np.pi,360)
xs = np.sin(a)
ys = np.cos(a)
ax.plot(xs,ys,c='r')
ax.text(-1,1,'sun')

# earth
xec,yec = 6,6
xe = xec+ 0.5*xs
ye = yec+ 0.5*ys
ax.plot(xe,ye,c='b')
ax.text(xec+0.5,yec+0.5,'earth')
ax.text(xec-1,yec+1.1,r"$\vec{v}$")

# r vector
ax.quiver([0,xec,xec-0.3],[0,yec,yec+0.1]
,[xec,-2,-3],[yec,2,-3]
,units='xy',scale=1,width=0.07)
ax.text(3,2,r"$\vec{r} = \vec{earth}-\vec{sun}$")
f_eq = (r"$\vec{F}=-G\frac{Mm}{|r|^2}\frac{\vec{r}}{|r|}\\$"
r"$=-G\frac{Mm}{|r|^3}\vec{r}$")
ax.text(0.5,5.5,f_eq)

# plot data
plt.show()

Result:

用 Python 模拟一个微型太阳系

文章出处登录后可见!

已经登录?立即刷新

共计人评分,平均

到目前为止还没有投票!成为第一位评论此文章。

(0)
上一篇 2022年6月14日 下午12:53
下一篇 2022年6月14日 下午2:07
本站注重文章个人版权,不会主动收集付费或者带有商业版权的文章,如果出现侵权情况只可能是作者后期更改了版权声明,如果出现这种情况请主动联系我们,我们看到会在第一时间删除!本站专注于人工智能高质量优质文章收集,方便各位学者快速找到学习资源,本站收集的文章都会附上文章出处,如果不愿意分享到本平台,我们会第一时间删除!