Orbit of the Earth - numerical methods leapfrog

Question

edit: I get these plots now from earths orbit. However, total energy looks not right. Should it not be oscillating between 0?

Thanks again!

I am attempting this numerical methods homework exercise part b). Orbit of the Earth

My goal is to modify my code from part a) so it returns kinetic and potential energy of earths orbit but I get this error:

(8.12).py:31: RuntimeWarning: overflow encountered in multiply
  z_vec[i+1,:]=z_vec[i,:] + h*f(z_half_step, t_vec[i])

My code for finding kinetic and potential energy of earths orbit is obviously wrong. My problem is that I don't know how to modify my code in a) to find potential and kinetic energy.

This is my code for a:

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


"""a"""


G = 6.6738*10**-11
M = 1.9891*10**30
h = 3600
y = 1.4710*10**11
vx = 3.0287*10**4

def LeapFrog(f, t_start, t_stop, z0, h):

    t_vec = np.arange(t_start, t_stop, h) 
    n = len(t_vec)
    d = len(z0)  
    z_vec = np.zeros((n, d))
    z_vec[0,:] = z0
    z_half_step=z_vec[0 , :] + (1/2)*h*f(z0,t_vec[0]) 
    
    
    for i in range(0, n - 1):
        z_vec[i+1,:]=z_vec[i,:] + h*f(z_half_step, t_vec[i])
        z_half_step += h*f(z_vec[i+1,:], t_vec[i])   


    return t_vec, z_vec, 


def f(z,t):   #dz/dt
    
    x=z[0]
    y=z[1] 
    vx=z[2] 
    vy=z[3] 
    r=np.sqrt(x**2+y**2)

    dz=np.zeros(4)
    
    dz[0]=vx 
    dz[1]=vy
    dz[2]=-G*M*x/r**3
    dz[3]=-G*M*y/r**3

    return dz

t_start = 0
t_stop = 24*365*h*5 #5 years
z0 = np.array([0,y,vx,0])
t_vec, z_vec = LeapFrog(f, t_start, t_stop, z0, h)



figure, ax = plt.subplots()

plt.title("Orbit of the Earth")
plt.plot(0,0,'yo', label = 'Sun positon')  
plt.plot(z_vec[:,0],z_vec[:,1], 'g', markersize=1, label='Earth trajectory')
ax.set_aspect('equal')

ax.set(xlim=(-1.55*10**11, 1.55*10**11), ylim = (-1.55*10**11, 1.55*10**11))
a_circle = plt.Circle((0, 0), 1.4710*10**11, fill=False,color='r')
ax.add_artist(a_circle)
ax.set_aspect('equal')

legend = plt.legend(['Sun position','Earths 5-year trajectory','Perfect circle'],loc='center left', bbox_to_anchor=(1, 0.5))

plt.show()



"""

We can see that the trajectory of the earth for 5 years is very slightly non-circular.

"""

For part b) I modified "dz2" and "dz[3]". Obviously it didn't work:

"""b"""

    m = 5.9722*10**24
    def f(z,t):   #dz/dt
    
        x=z[0]
        y=z[1] 
        vx=z[2] 
        vy=z[3] 
        r=np.sqrt(x**2+y**2)
    
        dz=np.zeros(4)
    
        dz[0]=vx 
        dz[1]=vy
        dz[2]=-G*M*m/r
        dz[3]=0.5*m*y**2
    
        return dz
    
    t_start = 0
    t_stop = 24*365*h*5 #5 years
    z0 = np.array([0,y,vx,0])
    LeapFrog(f, t_start, t_stop, z0, h)

Thanks

Answer 1

There should be no integration in part b), you use the result from part a) to compute the energies at each point.

    x,y,vx,vy = z_vec.T
    r = np.hypot(x,y)
    E_kin = 0.5*m*(vx**2+vy**2)
    E_pot = -G*M*m/r
    E = E_kin+E_pot
    plt.plot(t,E)