Orbit of the Earth - numerical methods leapfrog

[Graham87]

Homework Statement

Plot gravitational potential energy and kinetic energy

Relevant Equations

-GMm/r, 0.5mv^2

I am attempting this homework exercise part b).

I've modified my code but I get error overflow message. My goal is to modify my code so it returns kinetic and potential energy of Earth's orbit.
I made a new f(z,t) and modified the rows 99 and 100 with dz[2]=-G*M*m/r, and dz[3]=0.5*m*y**2 but obviously it didn't work.

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

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.

"""

"""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)

 

[pbuk]

The equation of motion is exactly the same as in part (a) so you don't need to change the integration. What information do you think you need to calculate kinetic and potential energy?

 

[Graham87]

The equation of motion is exactly the same as in part (a) so you don't need to change the integration. What information do you think you need to calculate kinetic and potential energy?

Aha! So I use the result from part a) to compute the energies:

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

plt.title('Energies function of time')
plt.plot(t_vec,E_pot)
plt.plot(t_vec,E_kin)
plt.plot(t_vec,E)
plt.legend(['Potential energy','Kinetic energy', 'Total energy'])
plt.show()

plt.title('Total E(t)')
plt.plot(t_vec,E)
plt.show()

I get these plots now. Did I plot it correctly? Shouldn't total energy oscillate between 0?

 

[pbuk]

I get these plots now. Did I plot it correctly?

Looks about right: the key point is that total energy is more or less constant, (with a relatively small oscillation which is an artefact of the leapfrog method: you'll probably cover this in the next lecture).