Astronomy & Celestial Mechanics
the sky, computed
Orbits

The singularity that is a change of variable

The Kepler problem blows up at zero separation, and a fixed-step integrator falls apart long before it gets there. Divide time by the radius and the equations become a harmonic oscillator — exactly, for every conic at once.

Assumes Universal variables and The ellipse.

18 min read 8 figures Everything is a conicTwo bodies, and no more

An inverse-square force has a hole in it. At zero separation the acceleration is infinite, the speed is infinite, and every numerical method ever written produces nonsense there. That much is expected of a formula with r2r^2 in a denominator. What is not expected is that the trouble starts a very long way from the hole: an orbit at eccentricity 0.9 never comes closer to its primary than a tenth of its semi-major axis, and a well-behaved integrator taking equal steps of time will still lose the plot on it.

The reason is not the force law. It is the choice of independent variable. Time is a poor coordinate for a Kepler orbit because the body spends most of it doing almost nothing and then does everything in a few per cent of it — and the fix, which is exact and not an approximation, is to stop using time.

One squaring, and the singularity is gone. Left, a harmonic oscillator: an ellipse centred on the origin, marked at 24 equal steps of its own phase. Right, the same points after squaring as complex numbers, u ↦ u². The image is an ellipse with the origin at a focus — checked here by the focal property, r₁ + r₂ = 2a to nine figures at four points — with semi-axes (A²+B²)/2 = 2.5000 and AB = 1.5000 and the focus at (A²−B²)/2 = 2.0000, which is ae exactly. Three things follow at once. The angle doubles, so one turn of the oscillator is one whole orbit traversed twice as fast in phase; the equal phase steps on the left arrive as the eccentric anomaly on the right, which is why that anomaly and not the true one is what the equations want; and the collision at r = 0, where the inverse square is infinite, is the point u = 0, where the oscillator has a perfectly ordinary velocity. The singularity was a property of the coordinates.
Fig. 1 The whole of the trick, in one map. On the left is a harmonic oscillator: an ellipse centred on the origin, with a linear restoring force and no singularity anywhere, marked at twenty-four equal steps of its own phase. On the right is the same set of points after squaring as complex numbers, z=u2z = u^2. The image is an ellipse with the origin at a focus — the focal property r1+r2=2ar_1 + r_2 = 2a holds to nine figures at every point checked — and it is a Kepler orbit. The collision at r=0r = 0, where the force is infinite, is the point u=0u = 0, where the oscillator is moving at an entirely ordinary speed. The singularity is a property of the coordinates.

What the map actually says

Take the plane as the complex numbers and write the physical position as zz. Set z=u2z = u^2 and, at the same time, replace the time tt by a new independent variable ss through

dtds=∣z∣=r.\frac{dt}{ds} = |z| = r.

Under those two substitutions together the equation of motion for the inverse-square force becomes

d2uds2−E2 u=0,\frac{d^2u}{ds^2} - \frac{E}{2}\,u = 0,

with EE the total energy. That is Hooke’s law. For a bound orbit E<0E < 0 and it is a harmonic oscillator of angular frequency −E/2\sqrt{-E/2}; for a hyperbolic orbit it is the same equation with the sign flipped and the solutions are hyperbolic sines; for the parabola E=0E = 0 and the motion in uu is uniform in a straight line. One linear equation covers every conic the force allows, and the branch that the universal formulation removes from the time equation has been removed from the equation of motion as well.

The geometry in the opening figure is why. Squaring doubles angles, so one circuit of the oscillator is one circuit of the orbit; and for uu tracing a centred ellipse with semi-axes AA and BB, the image has semi-axes (A2+B2)/2(A^2+B^2)/2 and ABAB about a centre displaced by (A2−B2)/2(A^2-B^2)/2 — and the focal distance a2−b2\sqrt{a^2-b^2} works out to that same displacement exactly. The origin is a focus by construction rather than by arrangement.

The step that chooses itself

The second substitution is the one that pays for itself in practice, and it can be understood without any complex arithmetic at all. Since dt/ds=rdt/ds = r, a fixed step in ss is a time step proportional to the current radius: short near periapsis, long near apoapsis, and set by nothing but where the body is.

The effect on a drawing is immediate.

36 steps round one orbit, spaced two ways. The same orbit at e = 0.9, marked at 36 equal steps of time on the left and 36 equal steps of the fictitious time χ on the right. Equal time steps are what a fixed-step integrator takes, and they put the marks where the body is slow: the longest gap is 11.8 times the shortest, and the pericentre passage — the only part of the orbit where anything interesting happens to the trajectory — is crossed in 1 of them. Stepping in χ instead makes the time step proportional to r, and the marks come out 2.25-to-one even, with 3 inside the same pericentre region. Nobody chose that. dt/dχ = r/√μ is the substitution itself, and since vr equals the angular momentum at both apses, a fixed step in χ is very nearly a fixed step in arc length — an adaptive integrator's behaviour obtained from a change of variable rather than from a controller.
Fig. 2 Thirty-six steps round one orbit at e=0.9e = 0.9, spaced two ways. Equal steps of time are what a fixed-step integrator takes, and they cluster where the body is slow: the longest gap between marks is 11.8 times the shortest, and the pericentre passage — the only stretch where the trajectory is doing anything — gets a handful of them. Equal steps of the fictitious time ss come out 2.25-to-one even, and nobody chose that. It follows from dt/ds=rdt/ds = r, and it goes further than a proportional time step: since the product vrvr equals the angular momentum at both apses and exceeds it only modestly in between, a fixed step in ss is very nearly a fixed step in arc length along the path.

That is an adaptive integrator’s behaviour obtained without an adaptive integrator. A step-size controller measures a local error estimate, compares it against a tolerance, accepts or rejects, and adjusts — machinery with parameters in it, and parameters are where a method’s behaviour stops being a theorem. The regularised step needs none of it. The variable itself carries the adaptation, and the same fixed step is right at every point of the orbit because the equation being stepped is a harmonic oscillator, whose solutions are smooth on a scale that does not vary at all.

Energy error over 200 revolutions, at one step size. The relative error in total energy against revolution number, for three integrators run on the same Kepler orbit at e = 0.7 with the same step of 400 per revolution. The exact energy is a constant, so every curve here is the method rather than the problem. Euler climbs steadily: its energy at the end is 127.8% wrong, and the orbit it draws has spiralled outwards. Runge–Kutta 4 begins 1.1e+0 times more accurate than leapfrog and ends at 5.18e-4, having grown by a factor of 10 across the run: the error is SECULAR. Leapfrog oscillates inside a band and stays there — worst error 6.07e-3, and the second half of the run is no worse than the first, which is measured here rather than claimed. That is the property that decides whether a five-billion-year integration means anything, and it is not accuracy: a symplectic method is the exact solution of a Hamiltonian a step-size away from the intended one, so its energy cannot wander, while a more accurate non-symplectic method has no such constraint and eventually wanders further.
Fig. 3 What is at stake. Three integrators run on the same eccentric orbit for two hundred revolutions, with the total energy — a constant of the motion — plotted against time. Euler’s method gains energy steadily, Runge–Kutta loses it steadily, and the symplectic leapfrog does neither: it is wrong about where the body is and right about how much energy it has. All three are being asked to resolve the pericentre passage with a step chosen for the whole orbit, and all three do their damage there — the run is drawn at four hundred steps per orbit because at a hundred and sixty the leapfrog’s energy band stops being bounded at all, which is the resolution failure this essay is about arriving in the one method that is supposed to be immune to it. The regularised form does not appear on this plot because it does not belong on it: with ss as the independent variable the orbit is a sinusoid, and every one of these methods integrates a sinusoid to machine precision.

Where the eccentric anomaly comes from

There is a reward for the change of variable that has nothing to do with numerics, and it explains something that otherwise looks like an accident of history.

For a bound orbit the fictitious time ss is proportional to the eccentric anomaly EE — specifically s=E/μ/a3⋅μ/as = E/\sqrt{\mu/a^3}\cdot\sqrt{\mu}/a, or more usefully, equal steps of ss are equal steps of EE. So the angle that Kepler introduced as a construction on the auxiliary circle, and that everyone since has treated as an intermediate quantity to be got rid of, is the natural time of the regularised problem. That reframing is worth stating plainly because it inverts the usual order of explanation. The eccentric anomaly is not a trick for making Kepler’s equation writable. Kepler’s equation is the price of insisting on the wrong independent variable, and the eccentric anomaly is the right one showing through.

One Newton loop, four conics, 6–11 steps. The residual |√μΔt − f(χ)| after each Newton step of the universal Kepler equation, at e = 0.6, e = 1, e = 1.0001, e = 1.4, all through the same code, the same starting rule and the same convergence test. Each orbit has periapsis at 0.5 AU and the solve is for 315.0 days past it. The step counts are 6 at e = 0.6, 10 at e = 1, 11 at e = 1.0001, 9 at e = 1.4, and the pair that matters is the crossing: 10 steps at e = 1 against 11 at e = 1.0001, a difference of 1. The spread across the whole set is 1.83, and it is a property of the STARTING GUESS rather than of the conic class — the run that begins nearest its own answer finishes first, on every branch. What the figure is about is the absence of a step at e = 1: the derivative of the equation being solved is r itself, which is positive on every conic, so the iteration has nothing to notice there. The classical forms have the opposite property — Kepler's equation degenerates as e → 1 from below, Barker's cubic is exact only AT 1, and e sinh H − H loses digits to cancellation just above it.
Fig. 4 The same idea one rung down the ladder, where it is about arithmetic rather than about dynamics. Newton’s method solving the universal time equation at four eccentricities — an ellipse, a parabola, a curve a ten-thousandth away from parabolic, and a hyperbola — with the residual after each iteration. All four converge at the same rate on the same code path, and the near-parabolic case is not special. A method that needed a different number of steps either side of e=1e = 1 would not have removed the branch, only hidden it.

Three dimensions, and why it took until 1965

The plane is easy because the plane is the complex numbers, and squaring is a perfectly good map on them. Space is not so obliging. Squaring in three dimensions has no meaning, and there is no three-dimensional algebra in which it would.

The resolution, found by Kustaanheimo and Stiefel in 1965, is to go up rather than sideways: the physical position in R3\mathbb{R}^3 is written as a quadratic function of a point in R4\mathbb{R}^4, using what amounts to quaternion multiplication. The extra dimension is not physical and is fixed by a constraint, and with it the same result holds — the regularised equations of motion in four dimensions are four uncoupled harmonic oscillators, all of the same frequency, with the energy as that frequency.

This is a genuinely surprising place for the Kepler problem to end up. The oldest solved problem in dynamics turns out to be a linear oscillator in disguise, and the disguise is a dimension count. It is the same fact that shows up elsewhere as the hidden symmetry of the inverse square — the conserved vector that keeps a bound orbit from precessing is a consequence of a four-dimensional rotation group acting on exactly this space, and Bertrand’s theorem about which force laws close their orbits is the statement that only two force laws admit a transformation of this kind.

What was actually measured, and where it is used

Nothing in this essay is an observation of the sky, and it would be dishonest to dress it as one. What can be pointed at is where the method earns its keep, and the numbers there are real.

Every modern high-precision solar-system ephemeris integrates regularised or partially regularised equations, and the reason is close encounters. A comet passing Jupiter at a hundredth of an astronomical unit, an asteroid grazing the Earth, a satellite of Saturn passing another — each is a pericentre passage about a body it is not orbiting, and each would eat a fixed-step integrator alive. The same is true of the star-cluster codes that follow a few thousand stars for a few billion years: binaries form, harden and pass through pericentre at separations of a few stellar radii, and it is entirely routine for such a code to spend most of its time on the tightest pair in the system. Regularisation of those pairs — chain regularisation, which applies the same substitution to a whole hierarchy of close bodies at once — is what makes the calculation finish at all.

Sundman, and a series that converges to nothing useful

The transformation has a curious place in the history of the subject, because it was invented for a problem it did not solve.

In 1912 Karl Sundman published a series solution of the three-body problem. That sentence is usually met with disbelief, since the three-body problem is famously the thing that has no solution — and both statements are true, because they mean different things. Sundman produced a power series in t1/3t^{1/3}, convergent for all real time, giving the positions of three gravitating bodies with the sole exception of the case where the total angular momentum vanishes and a triple collision is possible. The machinery that made it work was regularisation: binary collisions had to be shown to be removable before any series could be continued through them, and Levi-Civita’s transformation, published in its modern form eight years later, is the two-body core of that argument.

The catch is the one every account gives and few make concrete. The series converges, and it converges so slowly that estimates of the number of terms needed for astronomical accuracy run to 108imes10610^{8 imes10^{6}}. It is a solution in the sense a mathematician means and in no sense anybody else means, and its real legacy is the technique rather than the result.

What a cluster code does with a binary

The transformation earns its keep in one particular place, and it is worth describing because it explains why an obscure change of variables is standard equipment.

Simulating a star cluster means integrating thousands of bodies for billions of years. The difficulty is not the number of bodies; it is that a few of them form very tight binaries, whose orbital periods are shorter than the cluster’s dynamical time by many orders of magnitude. An integrator that must resolve a one-day orbit cannot advance a ten-million-year simulation.

Worse, those binaries are where the interesting physics happens. A hard binary is an energy reservoir: it hardens further in three-body encounters and gives up energy to the cluster, which is what halts the core collapse a purely collisional system would otherwise undergo. Removing them is not an option.

Regularisation is how they are kept without paying for them. A tight pair is taken out of the direct summation and integrated in regularised variables, where the motion becomes a harmonic oscillator and the step size no longer has to shrink at pericentre. The pair’s centre of mass remains in the main integration as a single particle, and the internal motion is advanced separately and accurately.

The technique extends to more than two bodies. A hierarchical triple, or a chain of several bodies undergoing a close encounter, is handled by chain regularisation: the bodies are ordered along a chain of relative coordinates, each link regularised, so that the singularities of every close pair are removed at once. That is what allows a cluster simulation to follow a resonant three-body encounter — the ones that eject stars and harden binaries — without either failing or spending the entire computation on it.

The change of variables is therefore not a curiosity about the two-body problem. It is the reason a class of astrophysical simulation is possible at all, and the reason those codes carry a separate module for anything that comes too close.

The same transformation, in a different subject

There is one further destination for the four-dimensional version, and it belongs to a subject that has nothing to do with orbits.

The Kepler problem and the hydrogen atom have the same potential. So the transformation that turns a three-dimensional inverse-square problem into a four-dimensional harmonic oscillator does so in quantum mechanics as well, and the Coulomb problem in three dimensions becomes an oscillator problem in four.

That matters because the oscillator is the one problem everybody can solve. Its spectrum is a ladder of equally spaced levels; mapping that back through the transformation produces the hydrogen atom’s 1/n21/n^2 energies, and the mapping explains their degeneracy — the accidental degeneracy of hydrogen’s levels, which looks unexplained in three dimensions, is the ordinary degeneracy of a four-dimensional oscillator.

The classical statement of the same fact is the conserved vector along the major axis, which is why the orbit closes and does not precess. The extra symmetry that keeps a Kepler ellipse from turning and the extra degeneracy in hydrogen’s energy levels are one thing, seen twice, and the change of variables in this essay is the map between the two views.

Where the trick stops working

Regularising a two-body collision succeeds completely, and it is worth asking why the same treatment does not dispose of every singularity in the problem — because the answer is a theorem rather than a shortage of cleverness.

A binary collision is, in the regularised variables, a perfectly ordinary point of the motion: the solution passes through it analytically and continues on the other side, describing a body that falls in, bounces and comes back out along the same line. That is not physics — real bodies do not bounce — but it is a legitimate analytic continuation, and it means the singularity was never in the dynamics, only in the coordinates.

A triple collision is different. Sundman proved that a triple collision can occur only if the total angular momentum of the system is exactly zero, which is a set of measure zero among initial conditions — so in that sense it never happens. But when it does, the singularity is genuine: the solution cannot be continued through it, because the way the three bodies approach one another does not settle down to a single limiting shape. Different sequences of times give different limiting configurations, and there is nothing for the continuation to continue to.

There is a stranger result still. A singularity in the NN-body problem need not be a collision at all. For five or more bodies there exist solutions in which no two bodies ever touch and yet the motion ceases to be defined after a finite time, because some of the bodies are thrown to infinity in that finite time — driven by an oscillating tight binary that feeds them energy. Such non-collision singularities were conjectured in the nineteenth century and constructed only in 1992.

So the change of variable in this essay is exactly as general as the two-body problem is, and the singularities it cannot remove are the ones that are not artefacts of a coordinate choice. The distinction is the point: a singularity that can be transformed away was never there, and one that cannot be is saying something about the system.

It also explains why the cluster codes of the previous section regularise pairs and chains rather than the whole system at once. There is no transformation that removes every singularity of an NN-body problem, because some of them are not removable; what the codes do is remove the ones that are, which happen to be the ones that occur.

Which is a fair statement of the whole method’s standing: it is a complete answer to a problem that occurs constantly and a partial answer to the general one, and the boundary between the two is a theorem rather than a matter of effort.

One more orbit shows the map applied where there is no near-singularity to remove.

One squaring, and the singularity is gone. Left, a harmonic oscillator: an ellipse centred on the origin, marked at 24 equal steps of its own phase. Right, the same points after squaring as complex numbers, u ↦ u². The image is an ellipse with the origin at a focus — checked here by the focal property, r₁ + r₂ = 2a to nine figures at four points — with semi-axes (A²+B²)/2 = 1.7500 and AB = 1.4000 and the focus at (A²−B²)/2 = 1.0500, which is ae exactly. Three things follow at once. The angle doubles, so one turn of the oscillator is one whole orbit traversed twice as fast in phase; the equal phase steps on the left arrive as the eccentric anomaly on the right, which is why that anomaly and not the true one is what the equations want; and the collision at r = 0, where the inverse square is infinite, is the point u = 0, where the oscillator has a perfectly ordinary velocity. The singularity was a property of the coordinates.
Fig. 5 The squaring map on a moderately eccentric orbit with a distant periapsis. The transformed curve is smooth and so was the original, so the map has cost a change of variable and bought nothing — which is the honest statement of when regularisation is worth doing.

Where the picture stops

Two honest limitations, and neither is small.

The first is that regularisation removes the singularity of the two-body problem and no other. A body falling into a point mass is fine; three bodies arriving at one place at one instant is a triple collision, and no change of variable removes it, because the problem genuinely is not analytic there. That result is Siegel’s, and it is one of the reasons the three-body problem is unsolvable in the specific sense that word has rather than merely difficult.

The second is that the transformed problem is only linear when the force is exactly inverse-square. Add an oblateness term, a third body, atmospheric drag or a relativistic correction and the equation in uu acquires a forcing term. It is still enormously better conditioned than the original — the singularity is still gone, and the step is still adaptive for free — but it is no longer a harmonic oscillator, and the exactness that makes the picture beautiful is the first thing a real application loses.

One analytic function each, and a case split that is only in the writing. Left: the Stumpff functions c₂ and c₃ across z = αχ², from z = −24 to 24. Positive z is an ellipse, negative z a hyperbola, and z = 0 exactly a parabola — and there is nothing at the origin: c₂ passes through 0.5 and c₃ through 0.166667 with every derivative continuous, because each is the single power series Σ(−z)^j/(k+2j)!. Written as (1 − cos√z)/z and (√z − sin√z)/z^3/2 they need three cases; written as their series they need none, and the case split in the code above is about floating point rather than about the mathematics. Right: how much floating point. The closed form subtracts two numbers that agree ever more closely as z → 0, so it loses digits — 11.9 of the sixteen it started with at z = 10⁻¹², measured against the series rather than estimated. That is why the drawn curve switches to the series inside |z| < 0.001, and why an implementation that uses the closed form everywhere is worst exactly at the parabola it was written to include.
Fig. 6 The functions the universal formulation runs on, and the numerical trap in them. c2(z)c_2(z) and c3(z)c_3(z) are entire functions of the energy parameter, perfectly smooth through z=0z = 0, which is what makes one solve cover ellipse, parabola and hyperbola. Their closed forms are not: written with the trigonometric or hyperbolic expressions, both lose every significant digit near the parabolic case, because they subtract two nearly equal numbers. The right-hand panel counts the digits lost. The remedy is to evaluate the series instead when ∣z∣|z| is small — a switch made on a threshold, and the one place in the whole method where a number has to be chosen by hand.

The map and the functions it produces are the two objects the trick consists of, and each is worth reading at a setting more extreme than the ones the essay uses to introduce it.

One squaring, and the singularity is gone. Left, a harmonic oscillator: an ellipse centred on the origin, marked at 24 equal steps of its own phase. Right, the same points after squaring as complex numbers, u ↦ u². The image is an ellipse with the origin at a focus — checked here by the focal property, r₁ + r₂ = 2a to nine figures at four points — with semi-axes (A²+B²)/2 = 6.0000 and AB = 1.8735 and the focus at (A²−B²)/2 = 5.7000, which is ae exactly. Three things follow at once. The angle doubles, so one turn of the oscillator is one whole orbit traversed twice as fast in phase; the equal phase steps on the left arrive as the eccentric anomaly on the right, which is why that anomaly and not the true one is what the equations want; and the collision at r = 0, where the inverse square is infinite, is the point u = 0, where the oscillator has a perfectly ordinary velocity. The singularity was a property of the coordinates.
Fig. 7 The squaring map applied to an orbit of eccentricity 0.95 with a periapsis at three tenths of the unit distance. The collision-bound trajectory becomes a smooth arc and the near-singular part of the motion is stretched out along it — the closer the approach, the more the map has to do and the more it does.
One analytic function each, and a case split that is only in the writing. Left: the Stumpff functions c₂ and c₃ across z = αχ², from z = −60 to 60. Positive z is an ellipse, negative z a hyperbola, and z = 0 exactly a parabola — and there is nothing at the origin: c₂ passes through 0.5 and c₃ through 0.166667 with every derivative continuous, because each is the single power series Σ(−z)^j/(k+2j)!. Written as (1 − cos√z)/z and (√z − sin√z)/z^3/2 they need three cases; written as their series they need none, and the case split in the code above is about floating point rather than about the mathematics. Right: how much floating point. The closed form subtracts two numbers that agree ever more closely as z → 0, so it loses digits — 11.9 of the sixteen it started with at z = 10⁻¹², measured against the series rather than estimated. That is why the drawn curve switches to the series inside |z| < 0.001, and why an implementation that uses the closed form everywhere is worst exactly at the parabola it was written to include.
Fig. 8 The Stumpff functions over sixty units of their argument. Each is a single analytic function through the elliptic, parabolic and hyperbolic cases, which is what removes the case split that every classical treatment carries — the branch is in the sign of one variable and in nothing else.

Where this ladder goes next

The regularised form is a statement about conditioning: the same motion, described in variables that make its difficulty disappear. That is a theme with more rungs on it than this one. The next asks a related question about the elements rather than about the trajectory — the six numbers that fix an orbit include three angles that are undefined for perfectly ordinary orbits, and the cure is the same in kind: not a better solver, but a better set of coordinates.

Further along the same line sit the questions of what a set of variables costs rather than what it hides. An orbit determined from angles alone is badly conditioned for a reason no substitution removes, because the information is genuinely absent rather than merely awkwardly parameterised. Telling those two cases apart — a difficulty that is in the coordinates and a difficulty that is in the data — is most of what numerical celestial mechanics consists of, and the Levi-Civita map is the cleanest example of the first that the subject has.

Where this goes next

Four conics through one periapsis, drawn by one solve. Distance from the Sun against time for four orbits sharing a periapsis of 0.5 AU, at e = 0.6, e = 1, e = 1.0001, e = 1.4, over 900 days. Every point on every curve came from the same universal Kepler solve — no branch on the conic class anywhere in it — and each curve was then checked against the classical solution of its own kind at the midpoint: Kepler's equation at e = 0.6 agrees to machine precision; Barker's cubic at e = 1 agrees to machine precision; e sinh H − H at e = 1.0001 agrees to machine precision; e sinh H − H at e = 1.4 agrees to machine precision. The curve to read twice is e = 1.0001: over this arc it lies within 0.02% of the parabola and is indistinguishable from it, and it is the only one of the four whose fate the drawing cannot show. What the classical parameterisation costs there is arithmetic rather than impossibility: at the midpoint of this arc, e sinh H − H throws away 3.3 of its sixteen digits to cancellation, against 0.2 at e = 1.4 — enough to matter to an ephemeris and not enough to stop a plot. The one solve that does not ask which conic it is Kepler's equation is for ellipses, Barker's cubic for parabolas, and a hyperbolic sine for the rest — three parameterisations of one motion, each worst exactly where its neighbour takes over. The universal variable removes the question, and the removal is not a convenience. Energy error over 240 revolutions, at one step size. The relative error in total energy against revolution number, for three integrators run on the same Kepler orbit at e = 0.5 with the same step of 200 per revolution. The exact energy is a constant, so every curve here is the method rather than the problem. Euler climbs steadily: its energy at the end is 106.5% wrong, and the orbit it draws has spiralled outwards. Runge–Kutta 4 begins 4.7e+2 times more accurate than leapfrog and ends at 2.22e-4, having grown by a factor of 10 across the run: the error is SECULAR. Leapfrog oscillates inside a band and stays there — worst error 2.62e-3, and the second half of the run is no worse than the first, which is measured here rather than claimed. That is the property that decides whether a five-billion-year integration means anything, and it is not accuracy: a symplectic method is the exact solution of a Hamiltonian a step-size away from the intended one, so its energy cannot wander, while a more accurate non-symplectic method has no such constraint and eventually wanders further. Wrong about where, and right about how much Runge–Kutta is the more accurate method and loses energy steadily; leapfrog is cruder and its energy error never leaves a band. Over five billion years only one of those properties survives — and neither method knows where the planet is. The three anomalies at E = 1.15 rad. The auxiliary circle construction. The eccentric anomaly E is measured at the centre, the true anomaly ν at the focus, and the mean anomaly M is time expressed as an angle. Here E = 1.150, ν = 1.827 and M = 0.602 radians, related by M = E − e sin E. The position that has no formula, and is computed anyway Kepler's equation relates where a body is to when it is there. It cannot be solved in elementary functions, Kepler said so, and nobody has managed it since — which has stopped nothing.

About the same objects

Not linked from either essay — found by the objects both name.

What links here

Essays that link to this one from their own argument.

The objects this essay names

Each one links to every other essay that touches it.

Collision singularityConditioningEccentric anomalyEnergy driftFictitious timeHarmonic oscillatorLevi civita transformationRegularisationStep size controlStumpff functions