I've looked at this 1-D SE question. That question seems to have difficulties around diverging potential functions, and I don't have that. I've looked at this 2-D 3-body demonstration, but that starts with equations of motion and I want to start with Lagrangians and symbolically derive the equations of motion.

The big picture is that I am developing a generalized Lagrangian solver and sanity checking it with small, standard scenarios. I got the original idea for the solver from the Wolfram blog, in this article by Moylan. My version of Moylan's solver passes his scenario and passes a damped simple harmonic oscillator, but it fails on the 2-D "reduced" Kepler problem. I don't know whether I have bugs or whether I am just using Mathematica clumsily, specifically with NDSolve.

I have checked this as much as I can visually and with cross references. For instance, the equations of motion I symbolically derive match those in this MIT courseware document.

I show all details below and will be very grateful for advice.

I only need 2-D because all 2-body Kepler problems are planar.

The "reduced" Kepler problem is a problem in two degrees of freedom, $r$ and $\theta$, concerning motion of a particle of "reduced mass" $m_{Earth}m_{Sun}/(m_{Earth}+m_{Sun})$ (again, see MIT).

Following Moylan, we write:

enter image description here

I hope you can see the reason for my interest in Moylan's method: the conciseness is valuable and worth generalizing to other problems. The Kepler problem doesn't have any non-conservative forces, but I have preserved that term in the code just for the purposes of that generalization. It works for a damped harmonic oscillator, for example.

coordinates = {r[τ], θ[τ]};
velocities = D[coordinates, τ];
L = 1/2 μ r'[τ]^2 + 1/2 μ r[τ]^2 θ'[τ]^2 - G M μ/r[τ];
nonConservativeForces = -{0, 0}*velocities;
equations = 
  MapThread[{q, v, h} \[Function] D[D[L, v], τ] - D[L, q] == h,
   {coordinates, velocities, nonConservativeForces}];

For the rest of this post, I will paste images of a Notebook. Self-contained code is in this gist that anyone can paste into a Notebook and run.

Mathematica finds the following equations of motion, which match the MIT reference up to units and constants:

enter image description here

We need numerical values, which I supply with a list of substitution rules:

enter image description here

And sanity-check by comparing the gravitational force to the centripetal force with these constants:

enter image description here enter image description here

Should be close enough.

I then solve the numerical equations with the following code:

enter image description here

The results plot as follows (please see the gist for the lengthy but uninformative plotting code):

enter image description here

OK, pretty much nothing is right. I expected $r(\tau)$ to be more-or-less constant, but it blows up. I expected $\theta(\tau)$ to be more-or-less linearly increasing, but it's not linear. Angular momentum and energy are not even slightly conserved.

I hope this is just something dumb on my part, but I don't see it. Anyone spot a problem that's easy to fix? Or have I waded into deeper waters?

EDIT: just realized that my numerical value for day is off by a factor of 365.25 (oops) changing it does NOT solve the bigger problem, though. The angular momentum and energy are still not conserved. I updated the public gist.

  • $\begingroup$ It would help, if instead of screenshots of code, there would be code... $\endgroup$
    – gwr
    Sep 5 '17 at 17:52
  • $\begingroup$ @gwr code is in this gist gist.github.com/rebcabin/180b54ca8d7dafc25806c6d80baada39, as stated in the body of the question :) $\endgroup$
    – Reb.Cabin
    Sep 5 '17 at 17:54
  • $\begingroup$ @J. M. I'm curious about the source of the Greek symbols in your (fine) edit. Did you just paste in Unicode symbols manually, or is there an automated method of copy-paste from Mathematica? $\endgroup$
    – Reb.Cabin
    Sep 6 '17 at 1:27
  • 1
    $\begingroup$ Personally, I know the codepoints and enter them directly. ;) Other people either use this or type in the $\LaTeX$ version (e.g. $\tau\mu$), and then directly copy the output of MathJax. $\endgroup$
    – J. M.'s torpor
    Sep 6 '17 at 1:31
  • 2
    $\begingroup$ The M.SE editor toolbar buttons that J.M. linked to in his comment are very handy. $\endgroup$
    – Chris K
    Sep 6 '17 at 1:39

Remember in your action to subtract your potential energy from the kinetic energy ${\cal L}= T- V$.

Forget the angular stuff for a second. You want $\partial_t (M r') = - G M/r^2$. This isn't the equation of motion you were getting -- see the glitch? You had set up a repulsive system instead of an attractive one.

The potential (you want) is attractive, so your ${\cal L}$ should be:

 L = 1/2 μ (r'[τ]^2 + r[τ]^2 θ'[τ]^2) + G M μ/r[τ];

The sign in front of $G M /r$ is switched relative to yours.

This allows your integration to proceed merrily, and in the period of one year, you can plot your orbit:

ParametricPlot[{r[τ] Cos[θ[τ]], r[τ] Sin[θ[τ]]} /. soln[[1]] /.
               τ -> t // Evaluate, {t, 0, tlim}]



Here's the entire output of your approach applied to a year with the attractive action (and notice the correct equation of motion in r):


  • $\begingroup$ Excellent catch! Energy conservation is still fluctuating in the fifth decimal place, but $r$ and $theta$ are behaving as expected. Much obliged. $\endgroup$
    – Reb.Cabin
    Sep 5 '17 at 20:10
  • $\begingroup$ From looking at L, it may be that there are numerical problems. $\endgroup$
    – evanb
    Sep 5 '17 at 20:55
  • $\begingroup$ @John Joseph M. Carrasco --- energy conservation demonstrated by removing erronious division by degree symbol and changing the sign in the potential-energy term. Your eyes on this much appreciated: another instances of Linus's Law at work en.wikipedia.org/wiki/Linus%27s_Law $\endgroup$
    – Reb.Cabin
    Sep 6 '17 at 0:03
  • 2
    $\begingroup$ For physical systems like this, one would usually want to switch to a Hamiltonian formulation, which would then allow you to use the "SymplecticPartitionedRungeKutta" method of NDSolve[]. $\endgroup$
    – J. M.'s torpor
    Sep 6 '17 at 0:58
  • 2
    $\begingroup$ Few quick comments: a) I always find it easier to work and debug in units close to one [ in this case: {AU, years, SolarMass, radians}]. This way you don't have to worry so much about noise from large cancellations. Not a real prob here, but good practice. (Can always convert to PlanckUnits if want large exponents later) b) might as well throw in a check for conserved quantities and analytic partial solns -- easy to see here EoM indep of 1 of your vars -- only it's derivative (conservation of ang. momenta). c) define Kin and Pot separately and combine them as necessary for L, H etc. $\endgroup$ Sep 6 '17 at 7:44

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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