# RK4 Gravity Simulator

I have the following RK4 solver which splits the two 2nd order ODEs, used to calculate x and y positions under the influence of a gravitating body where $$x''(t)=\frac{G m x(t)}{(x(t)^2+y(t)^2)^{3/2}}$$ and $$y''(t)=\frac{G m y(t)}{(x(t)^2+y(t)^2)^{3/2}}$$into four 1st order ODEs using the following code, and have attempted to plot the results which I had hoped would be a "spacecraft" in a 300km circular orbit about Earth:

Remove["Global*"]
(*dx/dt=*)f[t_, x_, y_, dx_, dy_] := -((G m x)/(x^2 + y^2)^(3/2));
(*dy/dt=*)g[t_, x_, y_, dx_, dy_] := -((G m y)/(x^2 + y^2)^(3/2));
(*d^2x/dt^2=*)u[t_, x_, y_, dx_, dy_] := dx;
(*d^2y/dt^2=*)v[t_, x_, y_, dx_, dy_] := dy;

G = 6.672*10^-11; (*Gravitational constant*)
m = 5.97219*10^24; (*Mass of Earth*)
r = 6.37101 *10^6; (*Earth mean equatorial radius*)
t[0] = 0; (*Initial time*)
x[0] = 0; (*Initial x-position*)
y[0] = r + 300000; (*Initial y-position*)
dx[0] = Sqrt[(G m)/(r + 300000)]; (*Initial x-velocity*)
dy[0] = 0; (*Initial y-velocity)
h = 0.005; (*Step size*)
tmax = 10000; (*Maximum running time*)

(*RK4 solver*)
Do[
{t[n] = t[0] + h n,

a1 = h f[t[n], x[n], y[n], dx[n], dy[n]];
b1 = h g[t[n], x[n], y[n], dx[n], dy[n]];
c1 = h u[t[n], x[n], y[n], dx[n], dy[n]];
d1 = h v[t[n], x[n], y[n], dx[n], dy[n]];

a2 = h f[t[n] + h/2, x[n] + a1/2, y[n] + b1/2, dx[n] + c1/2,
dy[n] + d1/2];
b2 = h g[t[n] + h/2, x[n] + a1/2, y[n] + b1/2, dx[n] + c1/2,
dy[n] + d1/2];
c2 = h u[t[n] + h/2, x[n] + a1/2, y[n] + b1/2, dx[n] + c1/2,
dy[n] + d1/2];
d2 = h v[t[n] + h/2, x[n] + a1/2, y[n] + b1/2, dx[n] + c1/2,
dy[n] + d1/2];

a3 = h f[t[n] + h/2, x[n] + a2/2, y[n] + b2/2, dx[n] + c2/2,
dy[n] + d2/2];
b3 = h g[t[n] + h/2, x[n] + a2/2, y[n] + b2/2, dx[n] + c2/2,
dy[n] + d2/2];
c3 = h u[t[n] + h/2, x[n] + a2/2, y[n] + b2/2, dx[n] + c2/2,
dy[n] + d2/2];
d3 = h v[t[n] + h/2, x[n] + a2/2, y[n] + b2/2, dx[n] + c2/2,
dy[n] + d2/2];

a4 = h f[t[n] + h, x[n] + a3, y[n] + b3, dx[n] + c3, dy[n] + d3];
b4 = h g[t[n] + h, x[n] + a3, y[n] + b3, dx[n] + c3, dy[n] + d3];
c4 = h u[t[n] + h, x[n] + a3, y[n] + b3, dx[n] + c3, dy[n] + d3];
d4 = h v[t[n] + h, x[n] + a3, y[n] + b3, dx[n] + c3, dy[n] + d3];

x[n + 1] = x[n] + 1/6 (a1 + 2 a2 + 2 a3 + a4);
y[n + 1] = y[n] + 1/6 (b1 + 2 b2 + 2 b3 + b4);
dx[n + 1] = dx[n] + 1/6 (c1 + 2 c2 + 2 c3 + c4);
dy[n + 1] = dy[n] + 1/6 (d1 + 2 d2 + 2 d3 + d4);
}, {n, 0, tmax}]

T1 = Table[{t[i], x[i]}, {i, 0, tmax}];
T2 = Table[{t[i], y[i]}, {i, 0, tmax}];
T3 = Table[{t[i], dx[i]}, {i, 0, tmax}];
T4 = Table[{t[i], dy[i]}, {i, 0, tmax}];

(*Graphical output*)
ListLinePlot[T1]
ListLinePlot[T2]
ListLinePlot[T3]
ListLinePlot[T4]

ListLinePlot[Table[{x[i], y[i]}, {i, 0, tmax}]]


When I compare the above output results to the following NDSolve code (where here $e=x$ and $p=y$), I see that my above code is incorrect and I do not get a circular orbit as is given below. Does anyone know where I might have gone wrong?

soln = NDSolve[{
e''[t] == -((G m   e[t])/(e[t]^2 + p[t]^2)^(3/2)),
p''[t] == -((G m   p[t])/(e[t]^2 + p[t]^2)^(3/2)),
e[0] == 0, p[0] == r + 300000, e'[0] == 7700, p'[0] == 0}, {e[t],
p[t]}, {t, 0, tmax}, Method -> "StiffnessSwitching",
MaxSteps -> 10000000]

Show[ParametricPlot[Evaluate[{{e[t], p[t]}} /. soln], {t, 0, tmax},
AxesLabel -> {e, p}, PlotStyle -> Automatic, PlotRange -> Full,
ImageSize -> Large]]


Any help would be appreciated, thanks very much.

-

## 3 Answers

Your problem is caused by a simple mistake. To solve a set of 2nd order ODEs with RK4, we need to rewrite the ODEs into a set of 1st order ODEs, like:

$$u'(t)=\frac{G m x(t)}{(x(t)^2+y(t)^2)^{3/2}}$$ $$v'(t)=\frac{G m y(t)}{(x(t)^2+y(t)^2)^{3/2}}$$ $$x'(t)=u(t)$$ $$y'(t)=v(t)$$

Then compare the above equations with the corresponding part in your code, it's easy to see that you've mixed up f, g and u, v, so the easiest way to fix your code is to exchange the position of f, g and u, v in the definition of them i.e.

(*du/dt=d^2x/dt^2=*)u[t_, x_, y_, dx_, dy_] := -((G m x)/(x^2 + y^2)^(3/2));
(*dv/dt=d^2y/dt^2=*)v[t_, x_, y_, dx_, dy_] := -((G m y)/(x^2 + y^2)^(3/2));
(*dx/dt=*)f[t_, x_, y_, dx_, dy_] := dx;
(*dy/dt=*)g[t_, x_, y_, dx_, dy_] := dy;


BTW, you've chosen too small a step size, something like h = 100; tmax = 50; is quite enough, the following is the corresponding result:

-

You can write your own algorithm and use it from NDSolve. For example, for RK4:

CRK4[]["Step"[rhs_, t_, h_, y_, yp_]] := Module[{k0, k1, k2, k3 },
k0 = h yp;
k1 = h rhs[t + h/2, y + k0/2];
k2 = h rhs[t + h/2, y + k1/2];
k3 = h rhs[t + h/2, y + k2];
{h, (k0 + 2 k1 + 2 k2 + k3)/6}]

CRK4[___]["DifferenceOrder"] := 4

CRK4[___]["StepMode"] := Fixed


Then call NDSolve[..., Method-> CRK4] after setting the options as you wish.

-
Thanks very much auxsvr, but I was hoping to figure out how to get the ODE system solved without having to make use of NDSolve, and was hoping someone would be able to show me where I went wrong in my RK4 algorithm. –  user7388 Jun 16 '14 at 12:09
Since this method is much simpler, it will be easier to spot the mistake this way. I'm afraid I can't help you any more than this, because it's been a long time since I used something other than NDSolve, and I copied this code from someone else. I've verified that it works, though. –  auxsvr Jun 16 '14 at 13:06
Damn, thanks anyway auxsvr; I use NDSolve a lot, so I wanted to make my own ODE solver without having to go back to replying on NDSolve again :P –  user7388 Jun 16 '14 at 19:15

Data

Remove["Global*"]
G = 6.672*10^-11;(*Gravitational constant*)
m = 5.97219*10^24;(*Mass of Earth*)
r = 6.37101*10^6;(*Earth mean equatorial radius*)
delT = 0.005;(*Step size*)
nSteps = 10000;


## NDSolve result

Clear[x, y, t];
soln = First@NDSolve[{
x''[t] == -((G m x[t])/(x[t]^2 + y[t]^2)^(3/2)),
y''[t] == -((G m y[t])/(x[t]^2 + y[t]^2)^(3/2)),
x[0] == 0,
y[0] == r + 300000,
x'[0] == Sqrt[(G m)/(r + 300000)],
y'[0] == 0},
{x[t], y[t], x'[t], y'[t]}, {t, 0, nSteps*delT},
Method -> "StiffnessSwitching", MaxSteps -> 10000000];

Grid[{{
Plot[Evaluate[x[t] /. soln], {t, 0, nSteps*delT}, PlotLabel -> "x(t)"],
Plot[Evaluate[x'[t] /. soln], {t, 0, nSteps*delT}, PlotLabel -> "x'(t)"]},
{Plot[Evaluate[y[t] /. soln], {t, 0, nSteps*delT}, PlotLabel -> "y(t)"],
Plot[Evaluate[y'[t] /. soln], {t, 0, nSteps*delT}, PlotLabel -> "y'(t)"]}
}]


## RK4 results

x = Table[0, {i, nSteps}, {j, 4}];
x[[1, 1]] = 0;(*Initial x-position*)
x[[1, 2]] = Sqrt[(G m)/(r + 300000)];(*Initial x-velocity*)
x[[1, 3]] = r + 300000;(*Initial y-position*)
x[[1, 4]] = 0;(*Initial y-velocity*)

f[{x1_, x2_, x3_, x4_}] := {x2, -(G m x1)/(x1^2 + x3^2)^(3/2), x4,
-(G m x3)/(x1^2 + x3^2)^(3/2)};

Do[k1 = f[x[[n]]];
k2 = f[  x[[n]] + delT  k1/2 ];
k3 = f[ x[[n]] + delT  k2/2];
k4 = f[  x[[n]] + delT  k3 ];
x[[n + 1]] = x[[n]] + delT/6 (k1 + 2 k2 + 2 k3 + k4), {n, 1, nSteps - 1}] ;


Plot it

t = Table[i*delT, {i, 0, nSteps - 1}];
Grid[{{
ListPlot[Transpose[{t, x[[All, 1]]}], PlotLabel -> "x(t)"],
ListPlot[Transpose[{t, x[[All, 2]]}], PlotLabel -> "x'(t)"]},
{ListPlot[Transpose[{t, x[[All, 3]]}], PlotLabel -> "Y(t)"],
ListPlot[Transpose[{t, x[[All, 4]]}], PlotLabel -> "Y'(t)"]}
}]


Notes: There is no need to break the system as you did. You can treat it as state space formulation. $x'(t) = f(x,t)$ where $x'(t)$ is vector, and $f(x,t)$ is the RHS. Hence you solve it as if it was one linear equation or 20 linear equations. It is the same method. only need to use $k1,k2,k3,k4$. I used $x_1$ for $x(t)$, and $x_2$ for $x'(t)$, and $x_3$ for $y(t)$ and $x_4$ for $y'(t)$ since this is very common way to do it in control systems.

-