Group Abstract

Message Boards

WOLFRAM COMMUNITY

13.7K Views

13 Replies

5 Total Likes

View groups...

Follow this post

Share this post:

GROUPS:

Astronomy Mathematics Physics Equation Solving Wolfram Language Modeling

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

Senlau Minto

Posted 7 years ago

I am using Mathematica 10.3. I want fo perform a computational analysis of the motion of the Earth around the sun based on Keplers 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. Please help me to fix my code.

I am using Mathematica 10.3. I want fo perform a computational analysis of the motion of the Earth around the sun based on Keplers 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.

POSTED BY: Senlau Minto

13 Replies

Sort By:

Senlau Minto

Posted 7 years ago

I tried it but the results of running the program are like this. 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.

POSTED BY: Senlau Minto

leonia dasilva

Posted 7 years ago

Please help...

POSTED BY: leonia dasilva

Neil Singer

Neil Singer, AC Kinetics, Inc.

Posted 7 years ago

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 BY: Neil Singer

Senlau Minto

Posted 7 years ago

Thank you very much..

POSTED BY: Senlau Minto

Senlau Minto

Posted 7 years ago

The right code to display 22 graphs. Like this: 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. Is my code correct?

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 BY: Senlau Minto

leonia dasilva

Posted 7 years ago

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

POSTED BY: leonia dasilva

Senlau Minto

Posted 7 years ago

thanks Oliver Seipel. I want to look like this: 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..

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 BY: Senlau Minto

Oliver Seipel

Posted 7 years ago

MultipleListPlot is replaced by ListPlot and ListLinePlot (since Version 6). Try this: ListLinePlot[xypts, PlotStyle -> {Hue[0], Hue[0.66]}, AspectRatio -> Automatic]

POSTED BY: Oliver Seipel

Senlau Minto

Posted 7 years 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 BY: Senlau Minto

Oliver Seipel

Posted 7 years 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] Regards OS

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 BY: Oliver Seipel

Senlau Minto

Posted 7 years ago

thanks Oliver Seipel.

POSTED BY: Senlau Minto

Senlau Minto

Posted 7 years ago

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

POSTED BY: Senlau Minto

Ox Clouding

Posted 7 years 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}] It's another story if you want to write your own solver though.

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 BY: Ox Clouding

Reply to this discussion

Reply Preview

Attachments

Remove Add a file to this post

Follow this discussion

or Discard

Feedback