9
$\begingroup$

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.

Improve this question
$\endgroup$
0

3 Answers 3

Reset to default
8
+50
$\begingroup$

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

Improve this answer
$\endgroup$
18
$\begingroup$

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.

Improve this answer
$\endgroup$
3
  • $\begingroup$ 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. $\endgroup$
    – InquisitiveInquirer
    Commented Jun 16, 2014 at 12:09
  • $\begingroup$ 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. $\endgroup$
    – auxsvr
    Commented Jun 16, 2014 at 13:06
  • $\begingroup$ 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 $\endgroup$
    – InquisitiveInquirer
    Commented Jun 16, 2014 at 19:15
6
$\begingroup$

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.

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.