as I said in the title, even though the code is running without any error messages, the code is not working properly. I'm applying RK4 method to the N-body simulation code, particularly solar system. The graph looks strange like there is no gravity between planets. Is my RK4 code running properly? I think the error is there because EULER method is just working fine.
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()
