# Simulating the motion of the Earth around the sun based on Kepler’s laws

I am using Mathematica 10.3. I want fo perform a computational analysis of the motion of the Earth around the sun based on Kepler’s laws.

Here is my code so far.

eulerStep[{t_, state_List}, h_, f_List] :=
{t + h, state + h Through[f[{t, state}]]}
solveSystemEuler [{t0_state0 _}, h_, n_Integer, f_List] :=
NestList[eulerStep[#, h, f] &, {t0, state0}, n]
midptStep[{t_, state_List}, h_, f_List] :=
{t + h, state + h Through[f[{t + 1/2 h, state + 1/2 h Through[f[{t, state}]]}]]}
solveSytemMidPt[{t0_, state0_}, h_, n_Integer, f_List] :=
NestList[midptStep[#, h, f] &, {t0, state0}, n]

L = 1/2 m (x'[t]^2 + y'[t]^2) + GMm/Sqrt[x[t]^2 + y[t]^2];

D[D[L, x'[t]], t] - D[L, x[t]] == 0
D[D[L, y'[t]], t] - D[L, y[t]] == 0

xdot[{t_, {x_, vx_, y_, vy_}}] := vx
vxdot[{t_, {x_, vx_, y_, vy_}}] := -x/(x^2 + y^2)^(3/2)
ydot[{t_, {x_, vx_, y_, vy_}}] := vy
vydot[{t_, {x_, vx_, y_, vy_}}] := -y/(x^2 + y^2)^(3/2)

start = {1, 0, 0, 1};
fcns = {xdot, vxdot, ydot, vydot};

orbit = solveSystemEuler[{0, start}, 0.01, 800, fcns];

<< StatisticsDataManipulation
xypts = Column[Column[orbit, 2], {1, 3}];
ListPlot[xypts, PlotJoined -> True];


Running the program gave the following error messages.

• The package you were trying to load is really old, and incompatible with the current system. Where did you get your code from? – J. M. will be back soon Oct 8 '18 at 4:23

Since today I corrected errors in this code in another question, I think that this question should also be answered

eulerStep[{t_, state_List}, h_, f_List] := {t + h,
state + h Through[f[{t, state}]]}
solveSystemEuler[{t0_, state0_}, h_, n_Integer, f_List] :=
NestList[eulerStep[#, h, f] &, {t0, state0}, n]
midptStep[{t_, state_List}, h_, f_List] := {t + h,
state + h Through[
f[{t + 1/2 h, state + 1/2 h Through[f[{t, state}]]}]]}
solveSytemMidPt[{t0_, state0_}, h_, n_Integer, f_List] :=
NestList[midptStep[#, h, f] &, {t0, state0}, n]

L = 1/2 m (x'[t]^2 + y'[t]^2) + GMm/Sqrt[x[t]^2 + y[t]^2];

D[D[L, x'[t]], t] - D[L, x[t]] == 0;
D[D[L, y'[t]], t] - D[L, y[t]] == 0;

xdot[{t_, {x_, vx_, y_, vy_}}] := vx
vxdot[{t_, {x_, vx_, y_, vy_}}] := -x/(x^2 + y^2)^(3/2)
ydot[{t_, {x_, vx_, y_, vy_}}] := vy
vydot[{t_, {x_, vx_, y_, vy_}}] := -y/(x^2 + y^2)^(3/2)

start = {1, 0, 0, 1};
fcns = {xdot, vxdot, ydot, vydot};

orbit = solveSystemEuler[{0, start}, 0.01, 800, fcns];

xypts = Table[{orbit[[i, 2]][[1]], orbit[[i, 2]][[3]]}, {i, 1,
Length[orbit]}];
ListPlot[xypts]
`

• Your answer is very good, very amazing. Thank you mr. alex Trounev for really helping me. :) – Teodora da Silva Nov 6 '18 at 5:29
• You're welcome! – Alex Trounev Nov 6 '18 at 10:53