Modeling the trajectory of the planet around the sun. Ellipse is not closing

Question

I was trying to create the projection of planets around the sun. Using the RungeKutta I'm trying to project and create the matplotlib graph. However, the out ellipse is not closing. Could you please help me, where exactly I'm doing the mistake?

Used Runge Kutta to find a vector R having the position as a function of time. And when I use the plot to draw the trajectory it doesn't make an ellipse which is what I am supposed to find.

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

#unités de normalisation 
UA=149.59787e6                       #distance moyenne Terre-Soleil
MASSE=6.0e24                         #masse Terre
JOUR=24*3600                      

#données : 
k=0.01720209895  
G=(k**2)                             # constante de gravitation en km^3/kg/s²
r0= 1                                # distance initiale terre soleil en m    
Ms = 2.0e30/MASSE                    #masse Soleil
Mt = 1                               #masse Terre

w0 = 30/(UA/JOUR)                    #vitesse angulaire en Km/s
C = r0*(w0**2)
m = (Ms*Mt/Ms+Mt)                    #masse réduite

K = G*m                              #on choisit K > 0 

b = 2 #beta


def RK4(F, h, r, theta, t, *args):
    k1=F(t,r,theta,)[0]
    l1=F(t,r,theta,)[1]
    k2=F(t+h/2,r+h*k1/2,theta+h*l1/2)[0]
    l2=(t+h/2,r+h*k1/2,theta+h*l1/2)[1]
    k3=F(t+h/2,r+h*k2/2,theta+h*l2/2)[0]
    l3=F(t+h/2,r+h*k2/2,theta+h*l2/2)[1]
    k4=F(t+h,r+h*k3,theta+h*l3/2)[0]
    l4=F(t+h,r+h*k3,theta+h*l3/2)[1]

    return [r +(h/6)*(k1+2*k2+2*k3+k4),theta +(h/6)*(l1+2*l2+2*l3+l4)]

def F1(t,r,theta):
    return np.array([r[1],r[0]*theta[1]**2-K/r[0]**b]),np.array([theta[1],-2*r[1]*theta[1]/r[0]])

def RK4N(F1,N,r0,r1,theta0,theta1,ta,tb):
    h=0.05
    ri=np.array([r0,r1])
    thetai=np.array([theta0,theta1])
    ti=ta
    R=np.zeros((N,2))
    THETA=np.zeros((N,2))
    lt=np.zeros(N)
    lt[0], R[0],THETA[0] = ta , ri ,thetai
    for i in range (1, N):
        a = ri
        ri = RK4(F1,h,ri,thetai,ti)[0]
        thetai=RK4(F1,h,a,thetai,ti)[1]
        ti=ti+h
        R[i]=ri
        THETA[i]=thetai
        lt[i]=ti
    return R,THETA,lt
def trace_position(F1,N,r0,r1,theta0,theta1,ta,tb):

    R,THETA,lt=RK4N(F1,N,r0,r1,theta0,theta1,ta,tb)
    lx,ly=[],[]
    for i in range(N):
        lx.append(R[i][0]*np.cos(THETA[i][0]))
        ly.append(R[i][0]*np.sin(THETA[i][0]))
    plt.plot(lx,ly)
    plt.plot(0,0)
    plt.show()

    # rapport a/b    
    max_x,min_x,max_y,min_y= max(lx),min(lx),max(ly),min(ly)
    rapport = (max_x-min_x)/(max_y-min_y)
    print("rapport a/b = ",rapport)
    if abs(rapport-1)< 10e-2:
        print("le mouvement est presque circulaire")
    else:
        print("le mouvement est elliptique")

def trace_Ep(F1,N,r0,r1,theta0,theta1,ta,tb):
    R,THETA,lt=RK4N(F1,N,r0,r1,theta0,theta1,ta,tb)
    lEp = []
    for i in range(N):
        lEp.append(-K/R[i][0]**(b-1))

    #fig, (ax1, ax2) = plt.subplots(1, 2)
    #ax1.plot(lt,lEp)
    #ax2.plot(R[:,0],lEp)
    plt.plot(lt,lEp)
    plt.show()

trace_position(F1,380,r0,0,0,w0,0,1)

Output:

Lutz Lehmann · Accepted Answer · 2020-01-27T15:16:22.913

You made a not so uncommon copy-paste error. There is an erroneous division by two left over in

k4=F(t+h,r+h*k3,theta+h*l3/2)[0]
l4=F(t+h,r+h*k3,theta+h*l3/2)[1]

Note also the missing F in the l2 line.

You can shorten these computation and avoid redundant or duplicated computations, and reduce places for errors by one half, by using tuple assignments like in

k4,l4 = F(t+h,r+h*k3,theta+h*l3)

and later

ri, thetai = RK4(F1,h,ri,thetai,ti)

Changing the step size computation to h=(tb-ta)/N as probably initially intended, and using tb=150 to get a closed loop, for a selection of subdivisions one gets the increasingly closed orbits via

for k in range(4):
    N = [16,19,25,120][k]
    plt.subplot(2,2,k+1,aspect='equal')
    trace_position(F1,N,r0,0,0,w0,0,150)

Stratubas · Answer 2 · 2020-01-27T14:22:47.963

3

A Runge-Kutta method can't give you a perfectly closed trajectory, because there is energy drift (energy is gradually decreasing).

If you need energy to remain close to the initial value, you can use a symplectic integrator. You can read more about this in Chris Rackauckas's answer on "What does 'symplectic' mean[...]".

However, as Lutz Lehmann pointed out, your time step is small enough, so that's not why your trajectory isn't visibly closed. This info might be interesting for others to know, so I'm leaving this answer here.

edited Jan 27 '20 at 14:22

answered Jan 27 '20 at 13:08

Stratubas

2,939
1
13
18

1

Your claim is quite wrong, to get that kind of gap in the elliptical orbit with a correctly implemented RK4 you need to reduce the number of points in one cycle dramatically to `N=19` (using ta,tb=0,150 and h=(tb-ta)/N). In general, the error also in the energy is O(h^4*t), slowly degrading in time but rather small even for modestly small step sizes. – Lutz Lehmann Jan 27 '20 at 14:13

Modeling the trajectory of the planet around the sun. Ellipse is not closing

2 Answers2