Planetary orbit shown as linear graph using rk4

Question

I am trying to simulate the orbit of a planet around a star using the Runge-Kutta 4 method. After speaking to tutors my code should be correct. However, I am not generating my expected 2D orbital plot but instead a linear plot. This is my first time using solve_ivp to solve a second order differential. Can anyone explain why my plots are wrong?

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
# %% Define derivative function
def f(t, z):
    
    x = z[0] # Position x
    y = z[1] # Position y
    dx = z[2] # Velocity x
    dy = z[3] # Velocity y
    
    G = 6.674 * 10**-11 # Gravitational constant
    M = 2 # Mass of binary stars in solar masses
    c = 2*G*M 
    r = np.sqrt(y**2 + x**2) # Distance of planet from stars
    
    zdot = np.empty(6) # Array for integration solutions
    zdot[0] = x
    zdot[1] = y
    zdot[2] = dx # Velocity x
    zdot[3] = dy #Velocity y
    zdot[4] = (-c/(r**3))*(x) # Acceleration x
    zdot[5] = (-c/(r**3))*(y) # Acceleration y
    
    return zdot
# %% Define time spans, initial values, and constants
tspan = np.linspace(0., 10000., 100000000)
xy0 = [0.03, -0.2, 0.008, 0.046, 0.2, 0.3] # Initial positions x,y in R and velocities
# %% Solve differential equation
sol = solve_ivp(lambda t, z: f(t, z), [tspan[0], tspan[-1]], xy0, t_eval=tspan)
# %% Plot 
#plot
plt.grid()
plt.subplot(2, 2, 1)
plt.plot(sol.y[0],sol.y[1], color='b')
plt.subplot(2, 2, 2)
plt.plot(sol.t,sol.y[2], color='g')
plt.subplot(2, 2, 3)
plt.plot(sol.t,sol.y[4], color='r')
plt.show()

Answer 1

With the ODE function as given, you are solving in the first components the system

xdot = x
ydot = y

which has well-known exponential solutions. As the exponential factor is the same long both solutions, the xy-plot will move along a line through the origin.

The solution is of course to fill zdot[0:2] with dx,dy, and zdot[2:4] with ax,ay or ddx,ddy or however you want to name the components of the acceleration. Then the initial state also has only 4 components. Or you need to make and treat position and velocity as 3-dimensional.

---

You need to put units to your constants and care that all use the same units. G as cited is in m^3/kg/s^2, so that any M you define will be in kg, any length is in m and any velocity in m/s. Your constants might appear ridiculously small in that context.

It does not matter what the comment in the code says, there will be no magical conversion. You need to use actual conversion computations to get realistic numbers. For instance using the numbers

G       = 6.67408e-11 # m^3 s^-2 kg^-1
AU      = 149.597e9 # m
Msun    = 1.988435e30 # kg
hour = 60*60 # seconds in an hour
day = hour * 24 # seconds in one day
year = 365.25*day # seconds in a year (not very astronomical)

one could guess that for a sensible binary system of two stars of equal mass one has

M = 2*Msun # now actually 2 sun masses
x0 = 0.03*AU
y0 = -0.2*AU
vx0 = 0.008*AU/day
vy0 = 0.046*AU/day

For the position only AU makes sense as unit, the speed could also be in AU/hour. By https://math.stackexchange.com/questions/4033996/developing-keplers-first-law and Cannot get RK4 to solve for position of orbiting body in Python the speed for a circular orbit of radius R=0.2AU around a combined mass of 2*M is

sqrt(2*M*G/R)=sqrt(4*Msun*G/(0.2*AU)) = 0.00320 * AU/hour = 0.07693 AU/day

which is ... not too unreasonable if the given speeds are actually in AU/day. Invoke the computations from https://math.stackexchange.com/questions/4050575/application-of-the-principle-of-conservation to compute if the Kepler ellipse would look sensible

r0 = (x0**2+y0**2)**0.5
dotr0 = (x0*vx0+y0*vy0)/r0
L = x0*vy0-y0*vx0           # r^2*dotphi = L constant, L^2 = G*M_center*R
dotphi0 = L/r0**2

R = L**2/(G*2*M)
wx = R/r0-1; wy = -dotr0*(R/(G*2*M))**0.5
E = (wx*wx+wy*wy)**0.5; psi = m.atan2(wy,wx)
print(f"half-axis R={R/AU} AU, eccentr. E={E}, init. angle psi={psi}")
print(f"min. rad. = {R/(1+E)/AU} AU, max. rad. = {R/(1-E)/AU} AU")

which returns

half-axis R=0.00750258 AU, eccentr. E=0.96934113, init. angle psi=3.02626659
min. rad. = 0.00380969 AU, max. rad. = 0.24471174 AU

This gives an extremely thin ellipse, which is not that astonishing as the initial velocity points almost directly to the gravity center.

orbit variants with half-day steps marked, lengths in AU

If the velocity components were swapped one would get

half-axis R=0.07528741 AU, eccentr. E=0.62778767, init. angle psi=3.12777251
min. rad. = 0.04625137 AU, max. rad. = 0.20227006 AU

This is a little more balanced.