Solving 2nd Order DEs with 4th Order RK Method

[Arman777]

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 .

[CODE lang="python" title="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()

[/CODE]

or

[CODE lang="python" title="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()[/CODE]

I am open to suggestions

 

[mfb]

I'm not sure what you mean by "couldnt uptade the distance".
One problem I see in the second method: You update rx and then calculate the change in y using the old ry but the new rx.

 

[Arman777]

I'm not sure what you mean by "couldnt uptade the distance".
One problem I see in the second method: You update rx and then calculate the change in y using the old ry but the new rx.

So I Should change the position of the rx. So you think that my code is correct ?

I changed my version 2 code.

 

[mfb]

Does it give the correct result for test cases?

 

[Arman777]

Does it give the correct result for test cases?

Not really, I guess it also does not matter at this point. It seems that the RK method is not the best choice to calculate a trajectory for an object.