Lagrangian form of Earth-Sun Kepler blows up NDSolve?

Question

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.

Answer 1

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

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):