RK4 Gravity Simulator

Question

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]]

enter image description here

Any help would be appreciated, thanks very much.

score 18 · Answer 1 · answered Jun 15 '14 at 22:57

18

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.

answered Jun 15 '14 at 22:57

auxsvr

517
2
11

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. – InquisitiveInquirer 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 – InquisitiveInquirer Jun 16 '14 at 19:15

score 8 · Accepted Answer · answered Jun 18 '14 at 10:25

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:

enter image description here

Nasser · Answer 3 · 2014-06-18T10:50:19.270

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)"]}
  }]

Mathematica graphics

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)"]}
  }]

Mathematica graphics

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.

RK4 Gravity Simulator

3 Answers3

NDSolve result

RK4 results

Linked

Related