Symplectic Integrators

As result of issues with the Three Body Problem mathematicians have opted to create special algorithms for use in multi-body problems as a method of keeping the system stable. Usually associated with the Hamiltonian and KAM theories, Symplectic Integrators are special geometric integrators which preserve the geometry of an orbiting system.

Wikipedia defines Symplectic Integrators as:

“ In mathematics, a symplectic integrator (SI) is a numerical integration scheme for Hamiltonian systems. Symplectic integrators form the subclass of geometric integrators which, by definition, are canonical transformations. ”

Integrator Comparison

Sun-Earth-Moon System

Computing the long term evolution of the solar system with geometric numerical integrators (Archive)
By mathematicians Shaula Fiorelli Vilmart (bio) and Gilles Vilmart (bio)

Abstract

“ Simulating the dynamics of the Sun–Earth–Moon system with a standard algorithm yields a dramatically wrong solution, predicting that the Moon is ejected from its orbit. In contrast, a well chosen algorithm with the same initial data yields the correct behavior. We explain the main ideas of how the evolution of the solar system can be computed over long times by taking advantage of so-called geometric numerical methods [Ed: See Numerical Solutions]. Short sample codes are provided for the Sun–Earth–Moon system. ”

The standard algorithms produce 'wrong solution' because the Moon is ejected from its orbit. A different algorithm is necessary to keep it together, and produces the 'correct' behavior.

Figure 7 from the paper shows a comparison between a non-symplectic and symplectic integrator:

The Sun-Earth-Moon System

“ Figure 7: Comparison of the explicit Euler method (left) and the symplectic Euler method (right) for the Sun-Earth-Moon system simulated over one year. The distance between the Moon (blue trajectory) and the Earth (black trajectory) is scaled by a factor of 100 in the plots, to better distinguish the Earth and the Moon. ”

The paper describes that the algorithm which keeps it together is the symplectic integrator.

“ The symplectic method, an integrator that preserves the energy well over long times. ”

Initial Data

“ In Table 1 we provide masses, positions, and initial velocities of the Sun, the Earth, and the Moon at a given date (here 1st of January 2016), and a value for the gravitational constant G

Table 1: Initial data from ^[8] for Sun, Earth, Moon on 01/01/2016 at 0h00. Masses are given relative to the Sun mass M_☉.

The Sun is chosen as reference and located at the origin (0, 0, 0). Distances are expressed in astronomical units ua, a quantity based on the Earth–Sun distance (1 ua is about 150 million kilometers), and the time is in Earth days ”

References

^[8]Observatoire virtuel de l’IMCCE, Portail système solaire, Miriade ephemeris generator, Observatoire de Paris & CNRS (2016), http://vo.imcce.fr/webservices/miriade.

A history of the n-body problem is given:

History - Is the Solar System Stable?

“ A question closely related to the topic of this snapshot is the issue of the stability of the solar system. Soon after Newton proposed his universal law of gravitation, many researchers, including, amongst others, Pierre-Simon Laplace (1749–1827), Joseph-Louis Lagrange (1736–1813), and Siméon Denis Poisson (1781–1840), were studying the question whether the regular trajectories of the planets will continue nicely until the end of times, or if collisions or ejections will occur.

In 1885, King Oscar II of Sweden sponsored a competition about this question. The prize was awarded to Henri Poincaré (1854–1912), although he did not really solve the problem. His contribution, however, is at the origin of the theory of dynamical systems. It also led to important developments in “Hamiltonian perturbation theory” and gave rise to the so-called Kolmogorov–Arnold–Moser (KAM) theory, which deals with the persistence of quasi-periodic motions under small perturbations, see the survey [7]. Unfortunately, this beautiful theory does not apply to realistic solar system models. ”

Conclusion

“ We have seen that the energy, a key invariant of all mechanical systems, is well preserved by the symplectic Euler method. In contrast, the explicit Euler method, and more generally any standard explicit Runge–Kutta methods, do not preserve it and are thus not suitable for integration over long time intervals. A mathematical theory called “backward error analysis” permits to demonstrate that symplectic integrators have a good energy conservation for such mechanical systems. ”

A special integrator, called symplectic integrators, which preserves the energy and prevents the body from escaping is used, and is therefore 'more suitable'. Does this sound like a full clean simulation of gravity?

Source Code

The above paper includes the source code for the above Sun-Earth-Moon System, for use with the free open source software Scilab, where the integrators can be tested.

Sun-Jupiter-Saturn System

Another paper states the same for the Sun-Jupiter-Saturn system:

“ Following [8, 7], let us consider the Sun-Jupiter-Saturn system, where for simplicity we neglect the other bodies and influences in the solar system. Surprisingly, applying a standard numerical method yields a dramatically wrong solution, where one of the planets is ejected from its orbit. In contrast, a well chosen symplectic integrator with the same initial data yields the correct behavior. ”

Initial Data

“ We provide in Table 1 the positions and initial velocities for the Sun, Jupiter and Saturn at a given date (here September 5th 1994), expressed in astronomical units, based on the Earth-Sun distance (1 A.U. is about 150 million kilometers), and the time is in earth days. ”

Always Returns Stable Conditions

The following is a paper which looks at symplectic and non-symplectic simulations and says that symplectic integrators always give stable conditions regardless of the perturbations which affects the body:

Numerical Integration Techniques in Orbital Mechanics Applications (Archive)

“ Revisiting Table I, the non-symplectic integrators do not give a stable solution even at time step of 1e-4. The solutions from these integrators are chaotic and shoot off to infinity.

~

The symplectic integrators are very good at keeping the orbit stable for even long period integrations. For all the integrators that provide stable solution at time step of 1e-3, increasing the integration time to six orbital period yields the same result. Symplectic integrators are particularly good at this since it can keep the errors bounded.

~

From the two problems analyzed in this paper, one can clearly see the advantages of the symplectic integrator over the non-symplectic integrators. The symplectic integrators produce consistently better results with higher accuracy and slightly less run time. The fourth order symplectic integrators in particular are extremely successful in propagating both the restricted three body problem and the simple two body problem. As discussed earlier, the symplectic integrator is able to achieve stable solutions at lower integration tolerance and run time. For long duration integration, which is often the case for most orbital mechanics applications, the error on the non-symplectic integrators are unbounded and can grow. On the other hand, the symplectic integrator, by keeping the Hamiltonian constant, is able to bound the error and prevent it from growing substantially, always able to return to stable condition after perturbations. ”

We read that the non-sympletic integrators do not give stable solutions, while symplectic integrators are "always able to return to stable conditions after perturbations".

Does an algorithm which always returns a system to stable conditions after perturbations sound like a legitimate reflection of bodies operating under Newton's laws?

Description

The purpose of the symplectic integrator is to preserve the area or geometry of the phase space.

University of Rochester
PHY411 Lecture notes Part 7 – Integrators (Archive)

Figure 4 shows different trajectories:

Further down in the document:

...

We read above that symplectic integrators are designed to preserve the geometry of phase space.

Solving the Hamiltonian

A paper Symplectic integration of hierarchical stellar systems by H. Beust describes that the symplectic integrator does not actually solve the real Hamiltonian of the problem:

https://www.aanda.org/articles/aa/full/2003/12/aa3133/aa3133.right.html (Archive)

“ An integration of any Hamiltonian system is symplectic when it exactly preserves the generalized areas in phase-space. This is in particular the case if the integration exactly solves (within computer round-off errors) a Hamiltonian. Symplectic integrators do not actually strictly solve the real Hamiltonian of the problem, but do it for a surrogate Hamiltonian assumed to be close to the real one. Once this is achieved it can be shown that there is no long-term drift of the energy. This is not the case for more classical methods (Runge-Kutta, Burlish & Stoer...) used in celestial mechanics, that are far from being symplectic. ”

Phase Space vs. Position Space

It has been argued that these definitions apply only to "Phase Space," and that Phase Space is entirely different than the normal Position Space that we know which uses the x, y, and z coordinates. It has been argued that Phase Space does not contain the Position Space coordinates.

In truth, Phase Space is merely Position Space with additional dimensions. We read a definition from Physics for Degree Students B.Sc Second Year (Archive):

10.4 Phase Space

“ A combination of the position space and momentum space is known as phase space. Thus phase space has six dimensions. A point in phase space is, therefore, completely specified by six coordinates x, y, z, p_x, p_y, p_z. Complete information about any particle in a dynamic system can be obtained from a knowledge of these six co-ordinates which completely determine its position as well as momentum. As there are n particles a knowledge of 6n co-ordinates gives complete information regarding position and momentum of all then particles in the phase space for a dynamic system. The concept of phase space is very useful while dealing with dynamical systems actually existing in nature. ”

Scrolling up to section 10.2 we read the definitions of Position Space, verifying that it is the normal x, y, z position coordinates:

10.2 Position Space

“ The three dimensional space in which the location of a particle is completely given by the three position coordinates, is known as position space. ”

Thus, if the geometry of Phase Space is preserved with a Symplectic Integrator, the geometry of Position Space is also preserved.

Figure Eight Example

Compare this figure eight three body problem in Phase Space and Position Space:

https://demonstrations.wolfram.com/PlanarThreeBodyProblemInPhaseSpace/ (Archive)

Compare the above to the Figure Eight Three Body Problem solution in Position Space, plotted on the x and y.

https://www.researchgate.net/publication/253681041_Dynamics_and_stability_of_multiple_stars (Archive)

Plotted on the x and y, and looks surprisingly like the Phase Space version.

It looks the same because it is the same. Phase Space is merely Position Space with detail in additional dimensions which represents momentum. The shape of an orbit represents geometry, as well as energy loss/gain.

Usage

Is the solar system stable? (Archive)
Scott Tremaine
University of Toronto and Institute of Astronomy, Cambridge

In the above paper Scott Tremaine performs an analysis of the Solar System, concluding that the Solar System will be stable for millions of years. It is seen that this analysis uses Symplectic Integration techniques:

Related Examples

Symplectic Integrators For Solar System Dynamics by Scott Tremaine
High precision symplectic integrators for the Solar System by Jacques Laskar, et all.
Symplectic Maps for the N-body Problem with Applications to Solar System Dynamics by Matthew Jon Holman

New Scientist - The World of Symplectic Space

A New Scientist article titled The World of Symplectic Space (Archive) by Robert McLachlan, Ph.D. (bio), provides a background and history of this subject, admitting that Newton's Laws of Motion cannot simulate the Solar System or the n-Body problems. It is necessary to use a workaround using Symplectic integration, which preserves area and geometry.

“ Conventional calculations of the future of the Solar System quickly degenerate into disarray as computer errors build up. Symplectic integration could save the day

WHAT will the Solar System be like in the distant future? Will Pluto and Neptune collide? Will the Earth be thrown into a different orbit by the combined gravitational pull from all the other planets? You might think that the answers are easily calculated. Just program a computer with Newton’s laws of motion, tell it the positions of the planets now, and wait while it grinds out the future of the Solar System for the next billion years. Right?

Wrong. With a calculation as complicated as this, the computer is almost certain to come up with the wrong answer. It is not that the computer, or even that the person programming it, makes mistakes. The problem arises because computer replaces real time with a series of snapshots. Consequently, the calculations a computer makes are not absolutely precise, so it can only provide us with an approximate picture of what will happen in the real world. Normally the errors are so small that they go unnoticed, but when computers are set to work on the enormously long string of calculations needed to simulate the movement of the planets round the Sun, tiny errors in each step can build up to make the final result wildly inaccurate.

Confusion reigns

The errors are inevitable because the equations describing the Solar System are so complicated that precise solutions cannot even be attempted. Problems arise when errors build up systematically or, worse, when the errors become chaotic. If this happens, the error in the calculated result will not only be large, but unpredictable too. They can lead to results that defy the laws of physics—in an extreme case the planets could spiral into the Sun, for example, or gain energy from nowhere and spin off into space.

These errors affect every computer model although their effects go unnoticed because calculations usually run for a short time preventing wild variations from building up. Even a simple pendulum could start swinging like a propeller if the simulation were left to run for long enough.

Mathematicians have discovered they can get around these problems if the computer model is built not on the laws of motion that apply in our familiar three-dimensional space, but on the geometrical laws of a much larger mathematical world called symplectic space. What we understand as movement in our space can be represented as pure geometry in the very different world of symplectic space. This geometry provides a much more efficient way of representing movement mathematically: while it cannot prevent the computer from introducing errors, it ensures that, whatever the errors, the outcome is a physically reasonable one.

The laws of geometry in symplectic space are applied through mathematical tools known as symplectic integrators: simple formulae that a computer can use to create reliable simulations of chaotic and complex aspects of real systems. Symplectic integration is already helping scientists to model the forces between tens of thousands of atoms in a crystal lattice—a system far too complex for conventional methods to handle reliably—and successfully predict properties of the material such as its strength or the way it vibrates.

Not every system in our Universe is symplectic, however. Dissipative forces such as friction or viscosity do not obey geometric laws and cannot be simulated using symplectic integration. Applying the method to weather forecasting is complicated, for example, because air resistance is a dissipative force and has to be ignored.

One example where symplectic simulations can help is in designing circular particle accelerators. The Main Ring accelerator at Fermilab in Illinois, which was built in 1972 before symplectic motion was understood, cannot store particles for long because small variations in their paths rapidly build up into uncontrollably large deviations. Physicists now realise that with the new symplectic methods they would have been able to simulate these effects, and design the machine to cope.

Conventional computer models use position and time to deduce velocity, but position and velocity are treated on an equal footing in symplectic space. For instance, in a model of the Solar System each planet is defined by six dimensions—three for its velocity in each direction and three for its position. The dimensions are coupled by special kinds of angles known as symplectic angles. These angles cannot be measured with a protractor: two lines superimposed on each other, for example, have a normal angle between them of 0° but a symplectic angle of 90°.

The special features of this weird space are its laws of geometry. That space and geometry are closely linked can be seen by looking at the properties of a triangle in two different types of space. We are taught at school that the angles of a triangle add up to 180°, but this is not always true if the triangle is not on a flat plane. For example, imagine a triangle on the surface of the globe, which has one corner at the North Pole and the other two on the equator. No matter What the angle at the pole, both the angles at the equator will be exactly 90°, so the three angles are bound to add up to more than 180°. By choosing the space carefully, mathematicians can arrange for certain geometric laws to hold true. In the example of the triangle, the angles add up to 180° only if the space is flat.

If the space is complex, then the geometric laws can be complex too. As the planets move in symplectic space according to the symplectic laws of geometry, so they move in three-dimensional space according to Newton’s laws of motion. Symplectic geometry ensures that symplectic angles remain the same as the planets move. This geometry can describe all the complicated motion of the Solar System.

Symplectic space has been hard to explore because its geometry is so unlike that of three-dimensional space (the name comes from the Greek word symplegma, meaning tangled or plaited). After many decades of study, the breakthrough came during the 1950s when the Russian mathematician Vladimir Arnol’d at the Moscow State University, along with Andrei Kolmogorov and Jiirgen Moser, proved a theorem that explains some of the implications that these hidden geometrical laws hold for real motion.

In the case of a single planet in a circular orbit around the Sun, the motion is nonchaotic and well understood. But one problem that could not be by conventional means is what happens when there is a second planet exerting a small gravitational pull on the first. Does the circular orbit merely become slightly elliptical? Does it develop a chaotic wobble? Or will the first planet wander off entirely? The Kolmogorov-Arnol’d-Moser (KAM) team proved that all three are possible. Given certain initial conditions, the orbit can still be regular. But if the starting conditions are slightly different, it could be chaotic. Once in a chaotic state, the orbit might even “leak out” in a process known as Arnol’d diffusion, which causes the planet to wander away from its circular orbit.

Alternative orbits

The problem lies in determining which type of orbit a system will adopt. The KAM theorem shows that chaos and order are infinitely mixed. Between any two regular orbits lie chaotic ones, and the planet could adopt any one of an infinite number of each type of orbit. But a planet that behaves nonchaotically can never become chaotic.

Traditional computer models of planetary orbits produce outrageous results because the build-up of errors in the calculation leads to results that run counter to the laws of motion on which the model is based. Symplectic integration avoids these pitfalls by modelling not just the forces and accelerations as happens in conventional computer simulations, but by also keeping symplectic angles fixed in symplectic space. While the computation will still, inevitably, accumulate errors, the KAM theorem guarantees that they will not nudge a planet into an impossible orbit. The result does not predict the exact motion of our Solar System, but it does provide useful information about it. For example, astronomers can work out Pluto’s distance from the Sun in a billion years’ time, but not which side of the Sun it will then be on.

Although scientists have started to use symplectic integration only recently, it is not a new idea. In the 1950s René de Vogalaére, a mathematician then at the University of Notre Dame in Indiana, suggested rewriting formulae to preserve symplectic angles at each step. But his paper was rejected by a mathematical journal and the idea was forgotten. In 1983, it was rediscovered independently by Ronald Ruth at Lawrence Livermore National laboratory in California and Feng Kang at the Chinese Academy of Sciences in Peking. Since then, the study of symplectic integration has gone from strength to strength and the technique is giving scientists new insights into the workings of the Solar System. That chaotic motion exists in the Solar System was first suggested in 1988 by conventional computer calculations. But these were painfully and unnecessarily slow, and symplectic integrations can be carried out in a fraction of the time to simulate this chaos far into the future. Scientists now model the entire lifespan of the Solar System, thought to be about 10 billion years, and know that the model is reliable because it obeys the laws of symplectic geometry. ”

Solar System Integrator Comparison

Below we find a simulation of the outer planets of the Solar System, showcasing the Forward Euler method versus the Symplectic Euler method.

(Archive 1 2)

Description: “ On the left the solar system is evolved forward in time using the Forward Euler method while on the right the Symplectic Euler method is used. Both schemes may be evaluated explicitly; however, it should be noted that the Symplectic Euler method is defined implicitly and is only made explicit due to the form of the Hamiltonian being separable into functions of purely the positions and momenta. The Forward Euler scheme does not preserve any properties of the system and is only 1st order accurate. The Symplectic Euler scheme is also only 1st order accurate, but it preserves the structure of the elliptical orbits and Hamiltonian providing the time step is reasonably small. Both methods used the same time step of 200 days with the rest of the parameters being drawn from Hairer, Lubich, and Wanner's text on geometric integration as in my previous video ”

We see in the above video that with the Euler Forward method that Jupiter is ejected from the Solar System after a single orbit around the Sun. Saturn is thrown beyond Neptune, &c., as the Solar System quickly degenerates.

Forward Euler Descriptions

Wikipedia states that the 'Forward Euler method' is also known as the 'Euler method':

“ In mathematics and computational science, the Euler method (also called forward Euler method) is a first-order numerical procedure for solving ordinary differential equations (ODEs) with a given initial value. ”

Descriptions from Science Direct sources state:

“ Forward Euler is the simplest numerical integrator. ”

Source Files

The author provides steps for reproducing his plots with MATLAB in the description of his linked associated video, which shows a successful Störmer–Verlet simulation of the Solar System:

This is a demonstration of the Stormer-Verlet method applied to the outer solar system consisting of Jupiter through Pluto. The computation was carried out in Matlab with a time step of 200 days, a final time of 200,000 days, distances in astronomical units (AU), and masses normalized by the Sun's mass. The masses of Mercury through Mars were lumped together with the Sun. The initial conditions for position and velocity along with the scaled masses were obtained from pg. 13-14 of the 2nd Ed. of "Geometric Numerical Integration: Structure-Preserving Algorithms for Ordinary Differential Equations" by Hairer, Lubich, and Wanner.

To reproduce this plot, download the main script and functions linked below and run the main script with the functions either in the same folder as the main script or somewhere else in Matlab's search path (for Windows users C:\Users\YourUserName\Documents\MATLAB). The only variable that needs editing in the script is "method = methods{index}", where "index" should be the number corresponding to the numerical method you want to use from the list given below.

The numerical methods provided include:

1) Forward Euler

2) Runge-Kutta 23

3) Runge-Kutta 45

4) Symplectic Forward Euler

5) Stormer-Verlet

Main script:

Solar_System_Example.m

Functions:

solarSysDiffEQ.m

solarSysHam.m

Störmer–Verlet is Symplectic

It is seen that the Störmer–Verlet is a symplectic method.

Paper: Qualitative study of the symplectic Störmer–Verlet integrator

Abstract: “ Symplectic numerical integrators, such as the Störmer–Verlet method, are useful in preserving properties that are not preserved by conventional numerical integrators. ”

Quotes

Symplectic integrators are used in other areas in science. From the abstract of a particle physics paper (Archive) from the 6th International Particle Accelerator Conference we read that they are required to keep simulations stable:

“ It has been long understood that long time single particle tracking requires symplectic integrators to keep the simulations stable ”

Symplectic integratos provide approximations of mechanical systems. From the abstract of a 2016 paper from Physical Review:

“ Explicit symplectic integrators have been important tools for accurate and efficient approximations of mechanical systems with separable Hamiltonians. ”

From the abstract of Hamiltonian Approximants for Symplectic Integrators we see:

“ Symplectic integrators do not, in general, reproduce all the features of the dynamics of the Hamiltonian systems which they approximate. ”

Numerical Solutions

From the abstract of Symplectic Runge-Kutta and related methods: recent results (Archive):

“ Symplectic algorithms are numerical integrators for Hamiltonian systems that preserve the symplectic structure in phase space. In long time integrations these algorithms tend to perform better than their nonsymplectic counterparts. ”

See Numerical Solutions - The available solutions for celestial systems do not fully simulate gravity. For instance, the following is a symplectic integrator which breaks gravity down into multiple two body problems, rather than using a single multiple body problem:

https://academic.oup.com/mnras/article/452/2/1934/1069988 (Archive)

“ In this paper, we present a new symplectic integrator for collisional gravitational N-body dynamics. The integrator is inspired by the non-symplectic and non-reversible integrator in Gonçalves Ferrari et al. (2014), SAKURA, and makes use of Kepler solvers. Like SAKURA we decompose the N-body problem into two-body problems. In contrast to SAKURA, our two-body problems are not independent. The integrator is reversible and symplectic and conserves nine integrals of motion of the N-body problem to machine precision. ”

Addendum

As admitted, Newton's Laws cannot simulate the Solar System. Exotic methods are required to preserve the area and geometry of an orbiting system.

On KAM theory New Scientist states "While the computation will still, inevitably, accumulate errors, the KAM theorem guarantees that they will not nudge a planet into an impossible orbit.", while another researcher remarks regarding KAM theory - "Unfortunately, this beautiful theory does not apply to realistic solar system models" ^[1]. One paper concludes that "the symplectic integrator, by keeping the Hamiltonian constant, is able to bound the error and prevent it from growing substantially, always able to return to stable conditions after perturbations." ^[2]

New Scientist tells us that that the simulations are not built on the laws of motion that apply in our familiar three-dimensional space ^[3], that the conventional models "quickly degenerate into disarray"^[4], and that these alternative geometric methods discussed cannot be used for prediction of the positions of bodies: "The result does not predict the exact motion of our Solar System, but it does provide useful information about it. For example, astronomers can work out Pluto’s distance from the Sun in a billion years’ time, but not which side of the Sun it will then be on." ^[5] Which is interesting since any direct and real simulation of motion would inherently know that Pluto makes an orbit every 248 years around the Sun, regardless of the author's downplay of 'a billion years'. This is also a contradiction of the popular claims that the Solar System can be truly simulated with Newton's laws out to millions or billions of years, or even for a shorter period of time barring the special geometry preserving methods which keeps it together.

All of this tells us that these methods are not realistic models of orbiting systems which can be used for prediction, and are really niche academic frivolities — the end result of a submission to the Three Body Problem and science's desperate attempt to get invalid laws and an unworkable system to work.

Instead of a working model based on gravity and the laws of motion we are presented with a system based on "geometric laws". The truth is that if Newton's laws could keep the Solar System together he would have demonstrated that himself in the 1700's, rather than concluding that divine intervention keeps the Solar System stable^[6], and which this geometric approach is the ultimate manifestation of. Newton's own conclusion that the Solar System is unworkable under his laws was true then and remains true today. We are once again faced with the reality that the accepted laws of gravity and motion are not, and cannot, be used to simulate the Sun-Earth-Moon system or the Solar System.