Simulating orbit of planet around the sun with RK4

Question

I am trying to simulate a planet going around the sun with the RK4 algorithm.

This is my code that i tried:

import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt

def calcvec(x1,y1,x2,y2):
    r = np.array([0,0,0])
    
    r[0]=x2-x1
    r[1]=y2-y1
    r[2]= (r[0]**2 + r[1]**2)**(3/2)
    
    return r

def orbit():
    dt = 0.001
    sx = 0.0
    sy = 0.0
    
    t = np.arange(0,100,dt)
    rx = np.zeros(len(t))
    ry = np.zeros(len(t))

    vx = np.zeros(len(t))
    vy = np.zeros(len(t))
    
    rx[0]=15.0
    ry[0]=0.0
    
    vx[0]=1.0
    vy[0]=1.0
    
    ms = 1
    
    for i in range(0,len(t)-1):
        k1x = vx[i]
        r = calcvec(rx[i],ry[i],sx,sy)
        k1vx = - (ms*r[0]/r[2])
        
        k2x = vx[i] + (dt/2)*k1vx
        r = calcvec((rx[i]+(dt/2)*k1x),ry[i],sx,sy)
        k2vx = -(ms*r[0]/r[2])
        
        k3x = vx[i] + (dt/2)*k2vx
        r = calcvec((rx[i]+(dt/2)*k2x),ry[i],sx,sy)
        k3vx = -(ms*r[0]/r[2])
        
        k4x = vx[i] + dt*k3vx
        r = calcvec((rx[i]+(dt)*k3x),ry[i],sx,sy)
        k4vx = -(ms*r[0]/r[2])
        
        rx[i+1] = rx[i] + (dt/6)*(k1x + 2*k2x + 2*k3x + k4x)
        vx[i+1] = vx[i] + (dt/6)*(k1vx + 2*k2vx + 2*k3vx + k4vx)
    
        print(str(k1vx) + ", " +str(k2vx) + ", " +str(k3vx) + ", " +str(k4vx))
    
        k1y = vy[i]
        r = calcvec(rx[i],ry[i],sx,sy)
        k1vy = - (ms*r[1]/r[2])
        
        k2y = vy[i] + (dt/2)*k1vy
        r = calcvec(rx[i],(ry[i]+(dt/2)*k1y),sx,sy)
        k2vy = -(ms*r[1]/r[2])
        
        k3y = vy[i] + (dt/2)*k2vy
        r = calcvec(rx[i],(ry[i]+(dt/2)*k2y),sx,sy)
        k3vy = -(ms*r[1]/r[2])
        
        k4y = vy[i] + dt*k3vy
        r = calcvec(rx[i],(ry[i]+(dt)*k3y),sx,sy)
        k4vy = -(ms*r[1]/r[2])
        
        ry[i+1] = ry[i] + (dt/6)*(k1y + 2*k2y + 2*k3y + k4y)
        vy[i+1] = vy[i] + (dt/6)*(k1vy + 2*k2vy + 2*k3vy + k4vy)
    
    fig, ax = plt.subplots()
    ax.plot(rx,ry, label='x(t)')
    ax.scatter(sx,sy)
    plt.title("orbit")
    plt.xticks(fontsize=10)
    plt.grid(color='black', linestyle='-', linewidth=0.5)
    plt.xlabel(r'x', fontsize=15)
    plt.ylabel(r'y', fontsize=15)
    plt.savefig("testtwobody.pdf")
    plt.show()


if __name__=="__main__":
    orbit()

When running this code i receive an "orbit" like this

which is obviously wrong, because I would expect an elliptical orbit around the sun. Therefore, there must be a grave error or some sort of misunderstanding on my part.

Thanks for your help in advance! Yours sincerly, chwu

Answer 1

First good night. OK! first the star is fixed at the origin of the Cartesian coordinate system and the planet describes a flat orbit around the star due to the mutual iteration of the two. The equations of motion are obtained by applying Newton's laws of dynamics in conjunction with the Newtonian theory of gravitation. After the physical analysis of the problem and calculations, we have an initial value problem. Observation! The dynamic equations of motion used in the code were taken from computational Physics Nicholas J. Giordano and Hisao Nakanishi Chapter 4 page 94 second edition. As we are using the python language, we can use methods from the scipy package to integrate the system of ordinary differential equations. The initial conditions are the same as what you provided in your code.

# Author      : Carlos Eduardo da Silva Lima
# Theme       : Movement of a Plant around a fixed star
# Language    : Python
# date        : 11/19/2022
# Environment : Google Colab
# Bibliography: computational Physics Nicholas J. Giordano and Hisao Nakanishi Chapter 4 page 94 second edition.

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint, solve_ivp

# Initial conditions of the planet
x0 = 15.0 # componente da posição inicial - x
y0 = 0.0  # componente da posição inicial - y
vx0 = 1.0  # componente da velôcidade inicial - vx
vy0 = 1.0 # componente da velôcidade inicial - vy
t_initial = 0.0
t_final   = 55.0
N = 1000

# Dynamic equations of motion of the planet around the Sun
# for odeint
def edo(r,t):
  x,y,vx,vy = r
  r3 = np.power((x**2+y**2),3/2)
  return np.array([vx,vy,-4*np.power(np.pi,2)*(x/r3),-4*np.power(np.pi,2)*(y/r3)])

# for solve_ivp
def edo_(t,r):
  x,y,vx,vy = r
  r3 = np.power((x**2+y**2),3/2)
  return np.array([vx,vy,-4*np.power(np.pi,2)*(x/r3),-4*np.power(np.pi,2)*(y/r3)])

# Integration of dynamic equations using odeint
t = np.linspace(t_initial,t_final,N)
r0 = np.array([x0,y0,vx0,vy0])
sol_0 = odeint(edo,r0,t)
#sol_1 = solve_ivp(edo_, t_span = [t_initial,t_final], y0 = r0, method='DOP853')

# New variables
x  = sol_0[:,0]
y  = sol_0[:,1]
vx = sol_0[:,2]
vy = sol_0[:,3]

# Plot
plt.style.use('dark_background')
ax = plt.figure(figsize = (10,10)).add_subplot(projection='3d')
ax.plot(x,y,0,'bo',0,0,0,'yo',lw=0.5)
ax.set_xlabel("X", color = 'white')
ax.set_ylabel("Y", color = 'white')
ax.set_zlabel("Z", color = 'white')
ax.set_title("star and planet")
plt.show()

Plot: star and planet

Hope I helped, see you :).