Message Boards Message Boards

GROUPS:

Simulate the motion of the Earth around the Sun based on Kepler's Law?

Posted 2 months ago
553 Views
|
13 Replies
|
5 Total Likes
|

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

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

Running the program gave the following error messages. enter image description here

Please help me to fix my code.

13 Replies
Posted 2 months ago

Hi,

If you just want to solve the equations, why not use the built-in NDSolve function?

Here is how I would do this:

lagrangian = 1/2 (x'[t]^2 + y'[t]^2) + 1/Sqrt[x[t]^2 + y[t]^2];
eq = Table[
  D[lagrangian, x1] - D[D[lagrangian, D[x1, t]], t] == 
   0, {x1, {x[t], y[t]}}];
sol = NDSolve[
  Join[eq, {x[0] == 1, x'[0] == 0, y[0] == 0, y'[0] == 1}], {x[t], 
   y[t]}, {t, 0, 5}];
ParametricPlot[{x[t], y[t]} /. sol[[1]], {t, 0, 5}]

enter image description here

It's another story if you want to write your own solver though.

Posted 2 months ago

I want to analyze how completing the simulation of the earth's motion around the sun is solved by the euler method

Posted 2 months ago

With minor code fixes it works:

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 (x'[t]^2 + y'[t]^2) + 1/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 = orbit\[Transpose][[2]]\[Transpose][[{1, 3}]]\[Transpose];
ListPlot[xypts]

Graphical Output

Regards OS

Posted 1 month ago

thanks Oliver Seipel.

Posted 1 month ago

Oliver Seipel please help me:

MultipleListPlot[xypts, PlotJoined -> True, 
  SymbolShape -> None, PlotStyle -> {Hue[0], Hue[0.66]}, 
  AspectRatio -> Automatic];

Is my code correct? Please help.

Posted 1 month ago

MultipleListPlot is replaced by ListPlot and ListLinePlot (since Version 6).

Try this:

ListLinePlot[xypts, PlotStyle -> {Hue[0], Hue[0.66]}, 
 AspectRatio -> Automatic]

enter image description here

Posted 1 month ago

thanks Oliver Seipel.

I want to look like this: enter image description here

what is the right code?

earthorbit = xypts;

spacestart = {1, 0, 0, 1.25}; orbit = 
 solveSytemMidPt[{0, spacestart}, 0.01, 2200, fcns];
spacehiporbit = 
  orbit\[Transpose][[2]]\[Transpose][[{1, 3}]]\[Transpose];
ListPlot[spacehiporbit]

ListLinePlot[earthorbit, spacehiporbit, 
 PlotStyle -> {Hue[0], Hue[0.66]}, AspectRatio -> Automatic]

I tried it but did not get the graphics I wanted. Please help me..

Posted 1 month ago

Please help me.. i really need your help Oliver Siepel.. thank you so much..

Use {} around your data for ListLinePlot -- it expects a list of data (which are in the form of a list as well).

ListLinePlot[{earthorbit, spacehiporbit}, 
 PlotStyle -> {Hue[0], Hue[0.66]}, AspectRatio -> Automatic]
Posted 1 month ago

Thank you very much..

Posted 1 month ago

The right code to display 22 graphs. Like this:

enter image description here

The graph is obtained from the following old code:

Table[MultipleListPlot[Take[earthorbit, {n, n + 100}], 
Take[spaceshiporbit, {n, n + 100}], PlotJoined -> True, 
SymbolShape -> None, PlotStyle -> {Hue[0], Hue[0.66]}, 
PlotRange -> {{-3.6, 1.2}, {-2, 2}}, 
AspectRatio -> Automatic];, {n, 1, 2200, 100}];

then I change the code as follows:

Table[ListLinePlot[Take[earthorbit, {n, n + 100}], 
Take[spaceshiporbit, {n, n + 100}], 
PlotStyle -> {Hue[0], Hue[0.66]}, 
PlotRange -> {{-3.6, 1.2}, {-2, 2}}, 
AspectRatio -> Automatic];, {n, 1, 2200, 100}];

Running the program gave the following messages. enter image description here

Is my code correct?

Posted 1 month ago

Please help...

Posted 1 month ago

I tried it but the results of running the program are like this.

enter image description here

I want to show is a few points along a hundred-long space orbit. notice from space how to get a shorter trace when it's further: a description of Kepler's II law. Is my code correct? what's the problem?

I hope someone can help me. Thank You.

Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract