Lagrangian form of Earth-Sun Kepler blows up NDSolve?

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: 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: We need numerical values, which I supply with a list of substitution rules: And sanity-check by comparing the gravitational force to the centripetal force with these constants:

Should be close enough.

I then solve the numerical equations with the following code: The results plot as follows (please see the gist for the lengthy but uninformative plotting code): 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.

• It would help, if instead of screenshots of code, there would be code...
– gwr
Sep 5 '17 at 17:52
• @gwr code is in this gist gist.github.com/rebcabin/180b54ca8d7dafc25806c6d80baada39, as stated in the body of the question :) Sep 5 '17 at 17:54
• @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? Sep 6 '17 at 1:27
• 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. Sep 6 '17 at 1:31
• The M.SE editor toolbar buttons that J.M. linked to in his comment are very handy. 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[] /.
τ -> t // Evaluate, {t, 0, tlim}]

yielding: 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): • Excellent catch! Energy conservation is still fluctuating in the fifth decimal place, but $r$ and $theta$ are behaving as expected. Much obliged. Sep 5 '17 at 20:10
• From looking at L, it may be that there are numerical problems. Sep 5 '17 at 20:55
• @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 Sep 6 '17 at 0:03
• 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[]. Sep 6 '17 at 0:58
• 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. Sep 6 '17 at 7:44