Equation of motion using Runge-kutta 4 and Verlet algorithm

[FahdEl]

Homework Statement

Hey, i used Runge-kutta 4 algorithm to solve the equation of motion (earth/sun) but the graph i get it s wrong,i thought my algorithm was wrong so i used Verlet algorithm but i had the same probleme

Relevant Equations

Equation of motion

import numpy as np
import matplotlib.pyplot as plt
G=6.67408e-11
M=1.989e30
m=5.972e24
X0=-147095000000
Y0=0
VX0=0
VY0=-30300
T=365*24*60
def rk(ax,ay,x,y,vx,vy,h):
    t=0
    n=int(T/h)
    A=[[t],[x],[y],[vx],[vy]]
    for i in range(1,n):
        k1x=vx
        k1y=vy
        q1x=ax(x,y)
        q1y=ay(x,y)
        k2x=vx+h*q1x/2
        k2y=vy+h*q1y/2
        q2x=ax(x+h*k1x/2,y+h*k1y/2)
        q2y=ay(x+h*k1x/2,y+h*k1y/2)
        k3x=vx+h*q2x/2
        k3y=vy+h*q2y/2       
        q3x=ax(x+h*k2x/2,y+h*k2y/2)
        q3y=ay(x+h*k2x/2,y+h*k2y/2)
        k4x=vx+h*q3x
        k4y=vy+h*q3y
        q4x=ax(x+h*k4x,y+h*k4y)
        q4y=ay(x+h*k4x,y+h*k4y)
        x+=h*(k1x+2*k2x+2*k3x+k4x)/6
        A[1].append(x)
        y+=h*(k1y+2*k2y+2*k3y+k4y)/6
        A[2].append(y)
        vx+=h*(q1x+2*q2x+2*q3x+q4x)/6
        A[3].append(vx)
        vy+=h*(q1y+2*q2y+2*q3y+q4y)/6
        A[4].append(vy)
        t+=h
        A[0].append(t)
    return A
#acceleration Newton
def ax(x,y):
    r=(x**2+y**2)**0.5
    return -G*M*x/(r**3)
def ay(x,y):
    r=(x**2+y**2)**0.5
    return -G*M*y/(r**3)

[CODE lang="python" title="Verlet algorithm"]import numpy as np
import matplotlib.pyplot as plt
G=6.67408e-11
M=1.989e30
m=5.972e24
X0=-147095000000
Y0=0
VX0=0
VY0=-30300
T=365
#acceleration Newton
def ax(x,y):
r=(x**2+y**2)**0.5
return -G*M*x/(r**3)
def ay(x,y):
r=(x**2+y**2)**0.5
return -G*M*y/(r**3)
#Verlet methode
def verlet(x,y,vx,vy,dt,Ax,Ay):
x_new = x+ vx*dt + Ax(x,y)*(dt**2)/2
y_new = y+ vy*dt + Ay(x,y)*(dt**2)/2
vx_new = vx + dt*(Ax(x,y) + Ax(x_new,y_new))/2
vy_new = vy + dt*(Ay(x,y) + Ay(x_new,y_new))/2
return (x_new,y_new,vx_new,vy_new)
#simulation
def coor(x,y,vx,vy,dt):
n=int(T/dt)
coord=[[0],[x],[y],[vx],[vy]]
t=0
for i in range(1,n):
x,y,vx,vy=verlet(x,y,vx,vy,dt,ax,ay)
t+=dt
coord[0].append(t)
coord[1].append(x)
coord[2].append(y)
coord[3].append(vx)
coord[4].append(vy)
return coord
[/CODE]

 

[Dr. Courtney]

Label your axes with quantity and units.

 

[FahdEl]

does it change the results ? i don't think so?what do u mean?
the probleme i must get an elliptique shape not this

 

[FahdEl]

this graph is for x axes and y axes