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.
$\tau\mu$
), and then directly copy the output of MathJax. $\endgroup$ – J. M.'s ennui♦ Sep 6 '17 at 1:31