我尝试实现RK4到N体的模拟.但从图表上看，代码工作不正常

Question

正如我在标题中所说，即使代码运行时没有任何错误消息，代码也不能正常工作。我正在将RK4方法应用于N体模拟代码，特别是太阳系。这张图看起来很奇怪，好像行星之间没有引力。我的RK4代码运行正常吗？我认为错误是存在的，因为EULER方法运行良好。

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
class Body():
"""
This class contains adjustable parameters as attributes. These
parameters include current and previous positions, velocities, and
accelerations.
"""
def __init__(self,
id, facecolor, pos,
mass=1, vel=None, acc=None):
self.id = id
self.facecolor = facecolor
self.pos = np.array(pos, dtype=float)
self.mass = mass
self.vel = self.autocorrect_parameter(vel)
self.acc = self.autocorrect_parameter(acc)
self.vector_pos = [self.pos]
self.vector_vel = [self.vel]
self.vector_acc = [self.acc]

def autocorrect_parameter(self, args):
if args is None:
return np.zeros(self.pos.shape, dtype=float)
return np.array(args, dtype=float)
def scalar_pos(self):
return np.sqrt(np.sum(np.square(self.vector_pos), axis=1))
def partial_step(point1, point2, dt):
ret=0
ret=point1+point2*dt
return ret   
class computinggravity():
"""
This class contains methods to run a simulation of N bodies that interact
gravitationally over some time. Each body in the simulation keeps a record
of parameters (pos, vel, and acc) as time progresses.
"""
def __init__(self, bodies, t=0, gravitational_constant=6.67e-11):
self.bodies = bodies
self.nbodies = len(bodies)
self.ndim = len(bodies[0].pos)
self.t = t
self.moments = [t]
self.gravitational_constant = gravitational_constant
def get_acc(self, ibody, jbody):
dpos = ibody.pos - jbody.pos
r = np.sum(dpos**2)
acc = -1 * self.gravitational_constant * jbody.mass*dpos / np.sqrt(r**3)
return acc
def update_accelerations(self,dt):
for ith_body, ibody in enumerate(self.bodies):
ibody.acc *= 0
k1 =0
k2 =0
k3 =0
k4 =0
acc_pos =0
acc_vel =0
for jth_body, jbody in enumerate(self.bodies):
if ith_body != jth_body:
acc = self.get_acc(ibody, jbody)
k1=acc*(jbody.pos-ibody.pos)
#acc_vel=partial_step(ibody.vel,k1,0.5)
acc_vel=ibody.vel+k1*0.5
#acc_pos=partial_step(ibody.pos,acc_vel,0.5)
acc_pos=ibody.pos+acc_vel*0.5

k2=(acc_pos-(acc_pos+acc_vel*0.5*dt))*acc
#acc_vel=partial_step(ibody.vel, k2, 0.5)
acc_vel=ibody.vel+k2*0.5
k3=(acc_pos-(acc_pos+acc_vel*0.5*dt))*acc
#acc_vel=partial_step(ibody.vel, k3, 1)
acc_vel=ibody.vel+k3*1
#acc_pos=partial_step(ibody.pos, acc_vel, 0.5)
acc_pos=ibody.pos+acc_vel*0.5

k4=(jbody.pos-(acc_pos+acc_vel*dt))*acc
#ibody.acc=ibody.acc+acc
ibody.acc=ibody.acc+(k1+k2*2+k3*2+k4)/6
ibody.vector_acc.append(np.copy(ibody.acc))

def updatingvelocity(self,dt):
for ibody in self.bodies:
#ibody.acc=self.update_accelerations(self)
ibody.vel=ibody.vel+ibody.acc*dt
ibody.vector_vel.append(np.copy(ibody.vel))
def updatingposition(self,dt):
for ibody in self.bodies:
ibody.pos=ibody.pos+ibody.vel * dt
ibody.vector_pos.append(np.copy(ibody.pos))
def simulation(self, dt, duration):
nsteps = int(duration / dt)
for ith_step in range(nsteps):
self.update_accelerations(dt)
self.updatingvelocity(dt)
self.updatingposition(dt)
self.t= self.t+ dt
self.moments.append(self.t)


# Masses
SOLAR_MASS = 1.988e30
EARTH_MASS = 5.9722e24
# Distances
ASTRO_UNIT = 1.496e11
MILE = 1609
# Durations
HOUR = 3600
DAY = 24 * HOUR
YEAR = 365 * DAY
def main():
m_sun = 1 * SOLAR_MASS
sun = Body('Sun', 'yellow', [0, 0, 0], m_sun)
m_mercury = 0.05227 * EARTH_MASS
d_mercury = 0.4392 * ASTRO_UNIT
v_mercury = (106_000 * MILE) / (1 * HOUR)
mercury = Body('Mercury', 'gray',
[d_mercury, 0, 0], m_mercury,
[0, v_mercury, 0])
m_earth = 1 * EARTH_MASS
d_earth = 1 * ASTRO_UNIT
v_earth = (66_600 * MILE) / (1 * HOUR)
earth = Body('Earth', 'blue', [d_earth, 0, 0], m_earth, [0, v_earth, 0])
m_mars = 0.1704 * EARTH_MASS
d_mars = 1.653 * ASTRO_UNIT
v_mars = (53_900 * MILE) / (1 * HOUR)
mars = Body('Mars', 'darkred', [0, d_mars, 0], m_mars, [v_mars, 0, 0])

m_jupiter = 318* EARTH_MASS
d_jupister= 5 * ASTRO_UNIT
v_jupiter = (28_000 * MILE)/ (1* HOUR)
jupiter = Body('Jupiter', 'red',[d_jupister,0,0], m_jupiter, [0,v_jupiter,0])

m_saturn = 95* EARTH_MASS
d_saturn= 9.5 * ASTRO_UNIT
v_saturn = (21_675 * MILE)/ (1* HOUR)
saturn = Body('Saturn', 'brown',[0,d_saturn,0], m_saturn, [v_saturn,0,0])

m_uranus = 14.5 * EARTH_MASS
d_uranus = 19.8 * ASTRO_UNIT
v_uranus = (15_233* MILE)/(1 * HOUR)
uranus = Body('Uranus', 'cyan', [d_uranus,0,0], m_uranus, [0,v_uranus,0])

m_neptune = 17 * EARTH_MASS 
d_neptune = 30 * ASTRO_UNIT
v_neptune = (12_146* MILE)/(1*HOUR)
neptune = Body('Neptune', 'blue', [0, d_neptune,0], m_neptune, [v_neptune,0,0])
bodies = [sun, mercury, earth, mars, jupiter,saturn, uranus, neptune]
dt = 1 * DAY
duration = 40 * YEAR
gd = computinggravity(bodies)
gd.simulation(dt, duration)
fig = plt.figure(figsize=(11, 7))
ax_left = fig.add_subplot(1, 2, 1, projection='3d')
ax_left.set_title('3-D Position')
ax_right = fig.add_subplot(1, 2, 2)
ax_right.set_title('Displacement vs Time')
ts = np.array(gd.moments) / YEAR
xticks = np.arange(max(ts)+1).astype(int)
for body in gd.bodies:
vector = np.array(body.vector_pos)
vector_coordinates = vector / ASTRO_UNIT
scalar = body.scalar_pos()
scalar_coordinates = scalar / ASTRO_UNIT
ax_left.scatter(*vector_coordinates.T, marker='.',
c=body.facecolor, label=body.id)
ax_right.scatter(ts, scalar_coordinates, marker='.',
c=body.facecolor, label=body.id)
ax_right.set_xticks(xticks)
ax_right.grid(color='k', linestyle=':', alpha=0.3)
fig.subplots_adjust(bottom=0.3)
fig.legend(loc='lower center', mode='expand', ncol=len(gd.bodies))
plt.show()
# plt.close(fig)
if __name__ == '__main__':
main()

Answer 1

您的问题是一个典型的问题，当在以粒子为中心的哲学上实现Euler/辛Euler/Verlet的结构化和适当的方法扩展到更高阶方法时，就会发生该问题。简而言之，它不起作用。

最大的区别是，对于高阶方法，使用不属于数值结果的临时状态(同时从时间步长开始保持原始状态(，甚至稍微偏离轨迹，请参阅三阶龙格-库塔可视化，了解如何想象这一点。或者换句话说，你需要为所有粒子完成龙格-库塔方法的第一阶段，然后为第二阶段的计算计算新的(虚拟(位置，然后在开始第三阶段等之前为所有粒子结束第二阶段。

总体目标应该是将模型的物理性质和数值积分方法分开。模型应该只实现一般的支持功能，例如，在一个状态下克隆自己，该状态被另一个状态的缩放导数或一组导数的线性组合所移动。因此RK4被实现为类似的东西

state.compute_derivatives()
temp1 = state.clone_shift_single(0.5*dt,state)
temp1.compute_derivatives()
temp2 = state.clone_shift_single(0.5*dt,temp1)
temp2.compute_derivatives()
temp3 = state.clone_shift_single(1.0*dt,temp2)
temp3.compute_derivatives()
newstate = state.clone_shift_multiple([dt/6,dt/3,dt/3,dt/6],[state,temp1,temp2,temp3])

另一种变体是，该模型实现以先前固定的顺序将粒子的状态和导数分量复制到平面向量和从平面向量复制粒子的状态分量和导数分量的函数，使得RK方法内部的向量运算可以通过numpy阵列等来实现。这意味着传递给求解器方法的导数函数；烹饪书；实现或scipy.integrate.solve_ivp之类的其他求解器，是看起来像的模型函数的包装器

def derivs(t,u):
model.copy_state_from(u)
model.compute_derivatives(t) # if there is some time dependent component
du_dt = np.zeros_like(u)
model.copy_derivatives_to(du_dt)
return du_dt