I tried implementing RK4 to N-body simulation. But the look from the graph the code is not working properly

Question

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()

Answer 1

Your problem is a typical one and occurs when the structured and appropriate approach to implement Euler/symplectic Euler/Verlet on a particle-centered philosophy is extended to higher order methods. In short, it does not work.

The big difference is that for higher-order methods one uses temporary states (while keeping the original state from the start of the time step) that are not part of the numerical result and even further slightly aside the trajectory, see Visualization of Third Order Runge-Kutta for how one could imagine that. Or in other words, you need to complete stage 1 of the Runge-Kutta method for all particles, then compute the new (virtual) positions for the computations of stage 2, then complete stage 2 for all particles before starting anything of stage 3 etc.

The overall aim should be to keep the physics of the model and the numerical integration method apart. The model should only implement the general support functionality, for instance cloning itself at a state that is shifted by the scaled derivatives of another state, or by the linear combinations of a set of derivatives. So that RK4 is implemented as something like

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])

Another variant is for the model to implement functions that copy the state and derivative components of the particles to and from a flat vector in a previously fixed order, so that the vector operations inside the RK method can be implemented by numpy arrays or the likes. This means that the derivatives function that is passed to the solver method, be it RK4 in a "cook book" implementation or some other solver as scipy.integrate.solve_ivp, is a wrapper for the model functions that can look like

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