问题描述
正如我在标题中所说，即使代码正在运行且没有任何错误消息，代码也无法正常工作。我正在将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,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],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,m_mercury,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,m_earth,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',d_mars,m_mars,[v_mars,0])
m_jupiter = 318* EARTH_MASS
d_jupister= 5 * ASTRO_UNIT
v_jupiter = (28_000 * MILE)/ (1* HOUR)
jupiter = Body('Jupiter','red',[d_jupister,m_jupiter,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',d_saturn,m_saturn,[v_saturn,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,m_uranus,v_uranus,0])
m_neptune = 17 * EARTH_MASS
d_neptune = 30 * ASTRO_UNIT
v_neptune = (12_146* MILE)/(1*HOUR)
neptune = Body('Neptune',d_neptune,m_neptune,[v_neptune,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)
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,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()
解决方法
您的问题是典型的问题，当将以粒子为中心的哲学来实施欧拉/渐近欧拉/维莱特的结构化且适当的方法扩展到高阶方法时，就会发生您的问题。简而言之，它不起作用。
最大的区别在于，对于高阶方法，人们使用的临时状态（从时间步开始就保持原始状态）不属于数值结果，甚至略微偏离轨迹，请参见{{ 3}}。换句话说，您需要为所有粒子完成Runge-Kutta方法的阶段1，然后为阶段2的计算计算新的（虚拟）位置，然后为所有粒子完成阶段2，然后再开始阶段3的任何操作，等等。
总体目标应该是使模型的物理性质和数值积分方法分开。该模型仅应实现通用支持功能，例如，将自身克隆到因另一状态的缩放导数或一组导数的线性组合而移动的状态下。这样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,temp1)
temp3.compute_derivatives()
newstate = state.clone_shift_multiple([dt/6,dt/3,dt/6],[state,temp1,temp2,temp3])
该模型的另一个变体是实现以下功能的功能：以先前固定的顺序在平面向量之间来回复制粒子的状态和导数分量，从而可以通过{{1 }}数组之类。这意味着传递给求解器方法的派生函数（无论是Visualization of Third Order Runge-Kutta还是其他numpy求解器）都是模型函数的包装，
scipy.integrate.solve_ivp