- #1

- 2,168

- 193

- Homework Statement
- Create a trajectory for the comets orbit

- Relevant Equations
- ##\ddot{x} = -GMx/r^3## and ##\ddot{y} = -GMy/r^3## for ##r = \sqrt{x^2 + y^2}##

In second order case we should rewrite the equation in terms of 2 first order DE's. So I wrote,

$$dx/dt = wx$$ $$dwx/dt = -GMx/r^3$$ and $$dy/dt = wy$$, $$dwy/dt = -GMy/r^3$$

Now I guess there's two ways to do it in 4th order RK method. I would either do it component by component or just in once. However when I do it component by component I couldn't uptade the distance .

or

I am open to suggestions

$$dx/dt = wx$$ $$dwx/dt = -GMx/r^3$$ and $$dy/dt = wy$$, $$dwy/dt = -GMy/r^3$$

Now I guess there's two ways to do it in 4th order RK method. I would either do it component by component or just in once. However when I do it component by component I couldn't uptade the distance .

Version1:

```
from numpy import arange, array
from math import sqrt
from pylab import plot, show, xlabel, ylabel
G = 6.67408 * 10 ** -11 #m^3kg^^-1s^-2
M = 1.989 * 10 ** 30 # kg
t_i = 0
t_f = 100
N = 10000
h = (t_f - t_i) / N
def f(r, t):
x = r[0] # dx/dt = wx
wx = r[1] #dwx/dt = -GMx/r^3
y = r[2] # dx/dt = wx
wy = r[3] #dwx/dt = -GMx/r^3
f_x = wx
f_wx = (- G * M * x) / sqrt(x**2 + y**2)**3
f_y = wy
f_wy = (- G * M * y) / sqrt(x**2 + y**2)**3
return array([f_x, f_wx, f_y, f_wy], float)
t_data = arange(t_i, t_f, h)
x_points = []
y_points = []
r = array([4*10**12, 0, 0, 500], float) # [x, vx, y, vy]for t in t_data:
x_points.append(r[0])
y_points.append(r[2])
k1 = h * f(r, t)
k2 = h * f(r + 0.5 * k1, t + 0.5 * h)
k3 = h * f(r + 0.5 * k2, t + 0.5 * h)
k4 = h * f(r + k3, t + h)
r += 1/6 * (k1 + 2 * k2 + 2 * k3 + k4)plot(x_points, y_points)
show()
```

or

Version2:

```
from numpy import arange, array
from pylab import plot, show, xlabel, ylabel
from math import sqrt
G = 6.67408 * 10 ** -11 #m^3kg^^-1s^-2
M = 1.989 * 10 ** 30 # kg
t_i = 0
t_f = 2000
N = 10000
h = (t_i - t_f) / N
def fx(rx, ry, t):
x = rx[0] # dx/dt = wx
wx = rx[1] #dwx/dt = -GMx/r^3
y = ry[0]
f_x = wx
f_wx = (- G * M * x) / sqrt(x**2 + y**2) **3
return array([f_x, f_wx], float)
def fy(rx, ry, t):
y = ry[0] # dx/dt = wx
wy = ry[1] #dwx/dt = -GMx/r^3
x = rx[0]
f_y = wy
f_wy = (- G * M * y) / sqrt(x**2 + y**2) **3
return array([f_y, f_wy], float)
t_data = arange(t_i, t_f, N)
x_points = []
y_points = []
rx = array([4*10**12, 0], float) #position and velocity respectively
ry = array([0, 500], float)
# for the x-component
for t in t_data:
x_points.append(rx[0])
y_points.append(ry[0])
k1 = h * fx(rx, ry, t)
k2 = h * fx(rx + 0.5 * k1, ry, t + 0.5 * h)
k3 = h * fx(rx + 0.5 * k2, ry, t + 0.5 * h)
k4 = h * fx(rx + k3, ry, t + h)
k1 = h * fy(rx, ry, t)
k2 = h * fy(rx, ry + 0.5 * k1, t + 0.5 * h)
k3 = h * fy(rx, ry + 0.5 * k2, t + 0.5 * h)
k4 = h * fy(rx, ry+ k3, t + h)
rx += 1/6 * (k1 + 2 * k2 + 2 * k3 + k4)
ry += 1/6 * (k1 + 2 * k2 + 2 * k3 + k4)
plot(x_points, y_points)
show()
```

I am open to suggestions

Last edited: