# Plotting elliptical orbits using Verlet method

I'm looking to solve a 2 body system consisting of the Earth and a body of mass in orbit.

I'm using the Verlet method to numerically integrate Newton's equation of motion in order to plot an elliptical orbit.

m = 10(*body of mass*)
e = 0.705(*eccentricity*)
M = 354.4223(*Mean Anomoly*)
Mo = 12.30336608(*Mean Motion*)
n = Mo*(2*Pi/86400)
μ = 3.986004418*10^14(*Grav parameter of Earth*)
ψ = M + 2*e*Sin[M] + 1.25*e^2*Sin[2*M](*True Anomoly*)
a = ((μ^(1/3))/(n^(2/3)))(*semi major axis*)
(*r = a (1 - e^2)/(1 + e*Cos[ψ])(*equation of ellitpical motion*)*)

Fg = (μ*m)/r^2(*gravitational force*)
acc = Fg/m(*initial aceleration*)

(*Implement Verlet Algorithm:*)

orb = {};
For[t = 0; rini = 1.15792*10^7(*initial position*);
vini = 4307.79(*initialvelocity*);
accini = 6.3168×10^6 (*initial acceleration*), t < 86400,
t += dt; j = accini;
d = 0.01;
r = rini + vini*dt + (accini*dt^2)/2;
AppendTo[orb, {vini, accini}];
v = vini + (j + accini)*dt/2];


This is what I have so far, but I'm stuck what do next to actually plot an elliptical orbit. Any guidance on if I'm implementing the Verlet method correctly would help.

• Welcome Luke. Right now, you just have two states, r and r’, making this a 1D problem. Since it is an orbit, you are going to need expressions the theta too, or do it in Cartesian coordinates with X and Y. Also, your acceleration should tied to the force, but I’m not seeing it. – MikeY Mar 25 at 22:28
• Also, consider changing your units (e.g. AU instead of kilometers) so you're not bandying about large magnitudes. – J. M. is away Mar 26 at 1:22

There are several problems with the OP's code. It's not clear what differential equation is to be solved, in what coordinate system, and so forth. The OP hasn't seen fit to clarify, so here's some code cobbed from an exercise I did to write an NDSolve method. It seems a good idea to separated the numerical method from the particular ode.

(* Stormer-Verlet method to solve q''[x] == f[q] with step size h *)
ClearAll[stormerVerlet];
stormerVerlet[f_, h_][ {q0_, p0_}] := stormerVerlet[f, h][N@ {q0, p0, f[q0]}];
stormerVerlet[f_, h_][ {q0_, p0_, f0_}] := Module[{pp, p1, q1, f1},
pp = p0 + h/2. f0;
q1 = q0 + h pp;
f1 = f[q1];
p1 = pp + h/2. f1;
{q1, p1, f1}
];


Here's a vector form of the gravitational field ODE $$q'' = - q\,/\, \| q \|^3$$ in an arbitrary dimensional space.. The function rhs[q] represents the right-hand side of the ODE.

ClearAll[rhs];
rhs[q_] := -q/(Norm@q)^3;


The first example shows a solution in the plane:

stepsize = 0.01;
pos0 = {1., 0.};
vel0 = {0., 0.5};
orb = NestList[stormerVerlet[rhs, stepsize], {pos0, vel0}, 600];
Graphics[{Line[orb[[All, 1]],
VertexColors -> ColorData["Rainbow"] /@ Array[# &, Length@orb, {0., 1.}]]},
Frame -> True]


Here are three orbits in space:

stepsize = 0.01;
pos1 = {1., 0., 0};
vel1 = {0., 0.5, 0.1};
orb1 = NestList[stormerVerlet[rhs, stepsize], {pos1, vel1}, 300];
pos2 = {1.5, 0., 0.2};
vel2 = {0., 0.8, -0.1};
orb2 = NestList[stormerVerlet[rhs, stepsize], {pos2, vel2}, 1200];
pos3 = {1.25, 0., -0.1};
vel3 = {0., 0.8, 0.2};
orb3 = NestList[stormerVerlet[rhs, stepsize], {pos3, vel3}, 1200];
Graphics3D[{Thick,
Line[orb1[[All, 1]],
VertexColors -> ColorData["Rainbow"] /@ Array[# &, Length@orb1, {0., 1.}]],
Line[orb2[[All, 1]],
VertexColors -> ColorData["Rainbow"] /@ Array[# &, Length@orb2, {0., 1.}]],
Line[orb3[[All, 1]],
VertexColors -> ColorData["Rainbow"] /@ Array[# &, Length@orb3, {0., 1.}]]}]


You can construct an InterpolatingFunction solution by putting the solution data in the form $$\{\{\{t_0\}, q_0, q'_0, q''_0\}, \{\{t_1\}, q_1, q'_1, q''_1\},\ldots\}$$:

orbFN = Interpolation[
Join[
ArrayReshape[
Subdivide[0., (Length[orb] - 1) stepsize, Length[orb]],
{Length[orb], 1, 1}],
orb,
2],
InterpolationOrder -> 5];

ParametricPlot[orbFN[t],
{t, 0, (Length[orb] - 1) stepsize},
ColorFunction -> (ColorData["Rainbow"][#3] &)]