Plotting orbits in python using integrate.solve_ivp

Question

Im trying to plot the orbit of jupiter around the sun (as well as 2 astroid clusters in the Lagrange points) using integrate.solve_ivp, but when i plot a graph of x position and y, I'm getting a spiral, rather than a stable orbit. Can anyone help?

import numpy as np
import scipy.integrate
import matplotlib.pyplot as plt

#takes the position of two masses and outputs the vector difference, and the modulus
def vectorise(x1,y1,z1,x2,y2,z2):
    v = np.array([x2-x1,y2-y1,z2-z1])
    return v, np.linalg.norm(v)

def derivatives(t, y):
    G =4*np.pi**2
    Mj = 0.001
    Ms = 1

    #Vij gives the vector pointing from i to j (leading to force on j from i)
    Vjs = vectorise(y[3],y[4],y[5],y[0],y[1],y[2])
    Vsg = vectorise(y[0],y[1],y[2],y[6],y[7],y[8])
    Vjg = vectorise(y[3],y[4],y[5],y[6],y[7],y[8])
    Vst = vectorise(y[0],y[1],y[2],y[9],y[10],y[11])
    Vjt = vectorise(y[3],y[4],y[5],y[9],y[10],y[11])

    return [y[12],y[13],y[14],#first differentials of sun position
    y[15],y[16],y[17],#first differentials of Jupiter position
    y[18],y[19],y[20],#first differentials of Greek position
    y[21],y[22],y[23], #first differentials of Trojan position
    -G*Mj*1/(Vjs[1]**3) *Vjs[0][0], #second differentail of y[12] (sun x)
    -G*Mj*1/(Vjs[1]**3) *Vjs[0][1], #second differentail of y[13] (sun y)
    -G*Mj*1/(Vjs[1]**3) *Vjs[0][2], #second differentail of y[14] (sun z)
    G*Ms*1/(Vjs[1]**3) *Vjs[0][0], #second differentail of y[15] (Jupiter x)
    G*Ms*1/(Vjs[1]**3) *Vjs[0][1], #second differentail of y[16] (Jupiter y)
    G*Ms*1/(Vjs[1]**3) *Vjs[0][2], #second differentail of y[17] (Jupiter z)
    -G*(Ms*1/(Vsg[1]**3) * Vsg[0][0] + Mj*1/(Vjg[1]**3) * Vjg[0][0]), #second differentail of y[18] (Greek x)
    -G*(Ms*1/(Vsg[1]**3) * Vsg[0][1] + Mj*1/(Vjg[1]**3) * Vjg[0][1]), #second differentail of y[19] (Greek y)
    -G*(Ms*1/(Vsg[1]**3) * Vsg[0][2] + Mj*1/(Vjg[1]**3) * Vjg[0][2]), #second differentail of y[20] (Greek z)
    -G*(Ms*1/(Vst[1]**3) * Vst[0][0] + Mj*1/(Vjt[1]**3) * Vjt[0][0]), #second differentail of y[21] (Trojan x)
    -G*(Ms*1/(Vst[1]**3) * Vst[0][1] + Mj*1/(Vjt[1]**3) * Vjt[0][1]), #second differentail of y[22] (Trojan y)
    -G*(Ms*1/(Vst[1]**3) * Vst[0][2] + Mj*1/(Vjt[1]**3) * Vjt[0][2])] #second differentail of y[23] (Trojan z)

def solver():
    solution = scipy.integrate.solve_ivp(
        fun = derivatives,
        t_span=(0,118.6),
        # initial velocities are given in sun frame, and are in AU/year, distances in AU (1AU = 1.496 *10^11 m)
        # jupiter has average orbital velocity of 13.07 km/s = 2.755 au/year
        y0 = (0.0, 0.0, 0.0, #initial values of y[0-2]; sun position
        0.0, 5.2, 0.0, #initial values of y[3-5] jupiter position
        -5.2*np.sqrt(3/4), 5.2*0.5, 0,#initial values of y[6-8]; greek position
        5.2*np.sqrt(3/4), 5.2*0.5 , 0,#initial values of y[9-11]; trojan position
        0,0,0,#initial values of y[12-14]; sun velocity
        2.755, 0, 0,#initial values of y[15-17]; jupiter velocity in AU per year
        2.755*0.5, 2.755*np.sqrt(3/4), 0,#initial values of y[18-20]; greek velocity in AU per year
        2.755*0.5, -2.755*np.sqrt(3/4), 0),#initial values of y[21-23]; trojan velocity in AU per year
        t_eval = np.linspace(0,118.6,  1000),
    )
    return solution

plt.plot(solver().y[3],solver().y[4])
plt.show()

It looks like the numerical method accuracy problem.. – Jan Stránský Aug 23 '20 at 14:27 — Jan Stránský, Aug 23 '20 at 14:27

score 1 · Accepted Answer · answered Aug 23 '20 at 14:34

1

These are typical symptoms of a wrong numerical method or wrong parameters to the method.

Reading the documentation, you can use several methods. For the default "RK45" I got what you described. However, using

scipy.integrate.solve_ivp(...,method="RK23",...)

I got nice ellipses.

answered Aug 23 '20 at 14:34

Jan Stránský

1,592
1
10
14

Thank you so much! Thats much much better. After a while, it still spirals, if I run it for, say 100000 years, is there a method that would make it stay stable for longer / indefinitely? – Jeremy Aug 24 '20 at 09:11
No. Theoretically you can make the numerical method result with infinitesimally small error (e.g. using infinitesimally small time step), but 1) it would be computationally too expensive and 2) there still would be computer round-off errors. If you have two bodies, you can use analytical solution. If you have more than two bodies, you are [doomed to not-exactly-accurate numerical solution](https://en.wikipedia.org/wiki/Three-body_problem).. – Jan Stránský Aug 24 '20 at 10:57
I've been reading about the RK methods, and it seems like RK45 is a more accurate method, do you know why in this case RK23 worked better? – Jeremy Sep 11 '20 at 17:18
1

I am not mathematician an know only roughly the background of the methods. As this is not really related to programming but is really a mathematical question, maybe you can try to ask on [math stack exchange](https://math.stackexchange.com/) or [math overflow](https://mathoverflow.net/) or similar? – Jan Stránský Sep 13 '20 at 19:01

Plotting orbits in python using integrate.solve_ivp

1 Answers1