N体问题的Python实现 [英] Python implementation of n-body problem issue
问题描述
我目前正在尝试使用欧拉方法求解微分方程来实现N体问题.但是，图形输出似乎不正确，并且我不确定经过一段时间的测试后代码中的问题在哪里.我目前正在使用半人马座A和B的近似值进行测试.这是我的代码:
将numpy导入为np导入matplotlib.pyplot作为plt从数学导入层#引力常数G = 6.67430e-11#天文单位au = 1.496e11sec_in_day = 60 * 60 * 24dt = 1 * sec_in_day类Body(object):def __init __(自身，名称，init_pos，init_vel，质量):self.name =名称self.p = init_posself.v = init_velself.m =质量def run_sim(body，t):mass = np.array([[b.m] for body in b]，dtype = float)#(n，1，1)vel = np.array([b.v for body in b]，dtype = float)#(n，1，3)pos = np.array([b.p for body in b]，dtype = float)#(n，1，3)名称= np.array([b中的b.b.name]，dtype = str)#保存绘图的位置和速度plt_pos = np.empty((floor(t/dt)，pos.shape [0]，pos.shape [1]))plt_vel = np.empty((floor(t/dt)，pos.shape [0]，pos.shape [1]))#重心com_p = np.sum(np.multiply(质量，pos)，轴= 0)/np.sum(质量，轴= 0)curr = 0我= 0而curr<t:dr = np.nan_to_num(pos [None ,:]-pos [:，None])r3 = np.nan_to_num(np.sum(np.abs(dr)** 2，axis = -1)**(0.5)).reshape((pos.shape [0]，pos.shape [0]，1))a = G * np.sum((np.nan_to_num(np.divide(dr，r3))* np.tile(mass，(pos.shape [0]，1)).reshape(pos.shape [0]，pos.shape [0]，1))，轴= 1)pos + = vel * dtplt_pos [i] = posvel + = a * dtplt_vel [i] = velcurr + = dt我+ = 1无花果= plt.figure(figsize =(15,15))斧= fig.add_subplot()对于列表中的我(范围(plt_pos.shape [1])):ax.plot(plt_pos [:，i，0]，plt_pos [:，i，1]，alpha = 0.5，label = names [i])ax.scatter(plt_pos [-1，i，0]，plt_pos [-1，i，1]，标记="o"，标签= f'{i}')plt.legend()plt.show()run_sim(身体= [身体('Alpha Centauri A'，[0，0，0]，[0,22345,0]，1.989e30 * 1.1)，正文('Alpha Centauri B'，[23 * au，0，0]，[0，-18100,0]，1.989e30 * 0.907)，]，t = 100 * 365 * sec_in_day) 这是结果图.我希望它们的轨道变通性较小，而更呈圆形，类似于维恩图式的形式.
正确绘制图有3个步骤.
首先也是最重要的是，正确实现物理模型的实现. r3 应该包含距离的三次方，因此平方根的三次方具有指数 1.5
r3 = np.nan_to_num(np.sum(np.abs(dr)** 2，axis = -1)**(1.5)).reshape((pos.shape [0]，pos.形状[0]，1)) 然后给出清理后的情节
请注意比例尺的差异，必须将图像水平压缩以在两个方向上获得相同的比例尺.
第二，这意味着初始速度过大，恒星彼此逃跑.这些速度可能恰好在恒星最靠近的位置.快速解决方法是将速度除以10.这样就可以得出图
可以通过评估和转换
对于欧拉方法来说，轨道向外螺旋形是典型的，因为各个步骤沿着切线移动到精确解的凸椭圆.
只有在RK4上使用5天的时间步长才能得到完美的椭圆形
对于RK4实施，最重要的步骤是将非平凡的导数计算打包到一个单独的子过程中
def run_sim(body，t，dt，method ="RK4"):...def acc(pos):dr = np.nan_to_num(pos [None ,:]-pos [:，None])r3 = np.nan_to_num(np.sum(np.abs(dr)** 2，axis = -1)**(1.5)).reshape((pos.shape [0]，pos.shape [0]，1))返回G * np.sum((np.nan_to_num(np.divide(dr，r3))* np.tile(mass，(pos.shape [0]，1)).reshape(pos.shape [0]，pos.shape [0]，1))，轴= 1) 然后可以使Euler跳出时间循环
def Euler_step(pos，vel，dt):a = acc(pos);返回pos + vel * dt，vel + a * dt 并类似地执行RK4步骤
def RK4_step(pos，vel，dt):v1 = vela1 = acc(pos)v2 = vel + a1 * 0.5 * dta2 = acc(pos + v1 * 0.5 * dt)v3 = vel + a2 * 0.5 * dta3 = acc(pos + v2 * 0.5 * dt)v4 = vel + a3 * dta4 = acc(pos + v3 * dt)返回pos +(v1 + 2 * v2 + 2 * v3 + v4)/6 * dt，vel +(a1 + 2 * a2 + 2 * a3 + a4)/6 * dt 选择类似方法
如果方法=="RK4"，则 stepper = RK4_step.其他Euler_step 然后时间循环采用通用形式
N =下限(t/dt)...对于范围(1，N + 1)中的i:pos，vel = stepper(pos，vel，dt)plt_pos [i] = posplt_vel [i] = vel I am currently trying to implement the N-body problem using Euler's method for solving differential equations. However, the graphical outputs do not seem correct, and I'm not sure where the issue in my code is after a while of testing. I'm currently using approximate values for Alpha Centauri A and B to test. This is my code:
import numpy as np
import matplotlib.pyplot as plt
from math import floor
# gravitation constant
G = 6.67430e-11
# astronomical units
au = 1.496e11
sec_in_day = 60 * 60 * 24
dt = 1 * sec_in_day
class Body(object):
def __init__(self, name, init_pos, init_vel, mass):
self.name = name
self.p = init_pos
self.v = init_vel
self.m = mass
def run_sim(bodies, t):
mass = np.array([[b.m] for b in bodies], dtype=float) # (n, 1, 1)
vel = np.array([b.v for b in bodies], dtype=float) # (n, 1, 3)
pos = np.array([b.p for b in bodies], dtype=float) # (n, 1, 3)
names = np.array([b.name for b in bodies], dtype=str)
# save positions and velocities for plotting
plt_pos = np.empty((floor(t/dt), pos.shape[0], pos.shape[1]))
plt_vel = np.empty((floor(t/dt), pos.shape[0], pos.shape[1]))
# center of mass
com_p = np.sum(np.multiply(mass, pos),axis=0) / np.sum(mass,axis=0)
curr = 0
i = 0
while curr < t:
dr = np.nan_to_num(pos[None,:] - pos[:,None])
r3 = np.nan_to_num(np.sum(np.abs(dr)**2, axis=-1)**(0.5)).reshape((pos.shape[0],pos.shape[0],1))
a = G * np.sum((np.nan_to_num(np.divide(dr, r3)) * np.tile(mass,(pos.shape[0],1)).reshape(pos.shape[0],pos.shape[0],1)), axis=1)
pos += vel * dt
plt_pos[i] = pos
vel += a * dt
plt_vel[i] = vel
curr += dt
i += 1
fig = plt.figure(figsize=(15,15))
ax = fig.add_subplot()
for i in list(range(plt_pos.shape[1])):
ax.plot(plt_pos[:,i,0], plt_pos[:,i,1], alpha=0.5, label=names[i])
ax.scatter(plt_pos[-1,i,0], plt_pos[-1,i,1], marker="o", label=f'{i}')
plt.legend()
plt.show()
run_sim(bodies = [ Body('Alpha Centauri A', [0, 0, 0], [0,22345,0], 1.989e30*1.1),
Body('Alpha Centauri B', [23 * au, 0, 0], [0,-18100,0], 1.989e30*0.907),
],
t = 100 * 365 * sec_in_day
)
And this is the resulting plot. I would expect their orbits to be less variant and more circular, sort of in a Venn diagram-esque form.
There are 3 steps to a correctly looking plot.
First and most important, get the implementation of the physical model right. r3 is supposed to contain the third powers of the distances, thus the third power of the square root has exponent 1.5
r3 = np.nan_to_num(np.sum(np.abs(dr)**2, axis=-1)**(1.5)).reshape((pos.shape[0],pos.shape[0],1))
This then give the cleaned up plot
Note the differences in the scales, one would have to horizontally compress the image to get the same scale in both directions.
Second, this means that the initial velocity is too large, the stars flee from each other. These velocities might be right in the position where the stars are closest together. As a quick fix, divide the velocities by 10. This gives the plot
Better initial values could be obtained by evaluating and transforming the supposedly more realistic data from https://towardsdatascience.com/modelling-the-three-body-problem-in-classical-mechanics-using-python-9dc270ad7767 or use the Kepler laws with the more general data from http://www.solstation.com/orbits/ac-absys.htm
Third, the mass center is not at the origin and has a non-zero velocity. Normalize the initial values for that
# center of mass
com_p = np.sum(np.multiply(mass, pos),axis=0) / np.sum(mass,axis=0)
com_v = np.sum(np.multiply(mass, vel),axis=0) / np.sum(mass,axis=0)
for p in pos: p -= com_p
for v in vel: v -= com_v
(or apply suitable broadcasting magic instead of the last two lines) to get the plot that you were probably expecting.
That the orbits spiral outwards is typical for the Euler method, as the individual steps move along the tangents to the convex ellipses of the exact solution.
The same only using RK4 with 5-day time steps gives prefect looking ellipses
For the RK4 implementation the most important step is to package the non-trivial derivatives computation into a separate sub-procedure
def run_sim(bodies, t, dt, method = "RK4"):
...
def acc(pos):
dr = np.nan_to_num(pos[None,:] - pos[:,None])
r3 = np.nan_to_num(np.sum(np.abs(dr)**2, axis=-1)**(1.5)).reshape((pos.shape[0],pos.shape[0],1))
return G * np.sum((np.nan_to_num(np.divide(dr, r3)) * np.tile(mass,(pos.shape[0],1)).reshape(pos.shape[0],pos.shape[0],1)), axis=1)
Then one can take the Euler step out of the time loop
def Euler_step(pos, vel, dt):
a = acc(pos);
return pos+vel*dt, vel+a*dt
and similarly implement the RK4 step
def RK4_step(pos, vel, dt):
v1 = vel
a1 = acc(pos)
v2 = vel + a1*0.5*dt
a2 = acc(pos+v1*0.5*dt)
v3 = vel + a2*0.5*dt
a3 = acc(pos+v2*0.5*dt)
v4 = vel + a3*dt
a4 = acc(pos+v3*dt)
return pos+(v1+2*v2+2*v3+v4)/6*dt, vel+(a1+2*a2+2*a3+a4)/6*dt
Select the method like
stepper = RK4_step if method == "RK4" else Euler_step
and then the time loop takes the generic form
N = floor(t/dt)
...
for i in range(1,N+1):
pos, vel = stepper(pos, vel, dt)
plt_pos[i] = pos
plt_vel[i] = vel
这篇关于N体问题的Python实现的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
