- #1
FahdEl
- 3
- 0
- 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
Python:
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]