A Symplectic Integrator
Contents
Introduction
On Orbiter Forum, orbinaut Keithth G has described some of his results comparing the Orbiter numerical integrator with one of his own. In response to a user question, he provided the following description of the principles underlying a symplectic numeric integrator.
A secondorder symplectic integrator
First, let's focus on a basic secondorder symplectic integrator. It is worthwhile spending a little time on this because, so long as the timestep of the integration is small enough, it actually works pretty well. It is also an integrator from which it is possible to build higherorder syplectic integrators.
A couple of introductory points:
 The integrator is called 'symplectic' because (up to second order at least) it preserves the symplectic character of the physical system. This means that the longrun properties of the integrator tend to be much more stable than other kinds of integrators and, largely because of this, they have become very popular amongst physicists in the last couple of decades.
 The integrator is secondorder because it has been constructed so as to 'kill' error terms up to the secondorder only.
So, what is this integrator? If we just focus on a onedimensional system for a moment, then the integrator maps a pair of points <MATH>\left\{Q_0,P_0\right\}</MATH> to a new pair of points <MATH>\left\{Q_2,P_2\right\}</MATH> at a time <MATH>\delta t </MATH> later. We can think of this as taking a pair of numbers, <MATH>\left\{Q_0,P_0\right\}</MATH>, that describe the state of some object at some start time <MATH>t = t_0</MATH> and updating this to a new pair of numbers, <MATH>\left\{Q_2,P_2\right\}</MATH>, that describes the new state of the same object at time <MATH>t = t_0 + \delta t</MATH>.
There are three steps to this updating rule:
<MATH>Q_{1} \leftarrow Q_{0}+\frac{1}{2}\,\delta t\, P_{0}</MATH>
<MATH>P_{2} \leftarrow P_{0}+\delta t\, F\left(Q_{1}\right)</MATH>
<MATH>Q_{2} \leftarrow Q_{1}+\frac{1}{2}\,\delta t\, P_{2}</MATH>
OK, so what do all of these symbols mean? Let's start with 'Q' and 'P'. In physics, largely because of a longstanding convention in nomenclature, 'Q' is often used to denote a spatial coordinate of something. Here, in our onedimensional example, you can think of 'Q' as representing the xcoordinate of some object moving in some gravitational field. In the same convention, 'P' is used to denote the momentum of the same object. (Formally, it is the generalised momentum conjugate to 'Q' but that's a nuance that we don't need to worry about here.). Here, we can think of 'P' as representing the xcoordinate of the momentum of a particle. But since momentum is just mass * velocity (in this coordinate system), we can think of 'P' as the xcoordinate of velocity (multiplied by the mass of the object). So, the pair of numbers <MATH>\left\{Q_0,P_0\right\}</MATH> just represents the position and velocity of the object at some initial time. And this is just the statevector of the object written in cartesian coordinates. In other words, the symplectic integrator is an updating rule that takes an object's state vector at some initial time, and returns a new state vector at some time <MATH>\delta t </MATH> later. (Now, you can extend this idea to three dimensions but for the time being we'll just stick with onedimension.)
Next, let's focus on the 'F' term. Again by convention, 'F' is used to denote the force acting on an object. Specifically, <MATH>F\left(Q_{1}\right)</MATH> represents the force on the object being modelled when it is at position <MATH>Q_1</MATH>. So, to carry out the three steps of the integration updating rule, we need to provide some force function, <MATH>F(Q)</MATH> which allows us to calculate the force on an object at a point in our onedimensional space (i.e., at any 'Q'). For object moving subject to a gravitational force from a single body, we know that:
<MATH> F(Q) = \frac{\mu\, m}{(QQ^*)^{2}}</MATH>
(provided that <MATH>Q > Q^*</MATH>)
where <MATH>\mu</MATH> is the gravitational constant for the body in question; 'm' is the mass of the object that we are modelling, and <MATH>Q^*</MATH> is the location of the source of the gravitational field (e.g., the centre of the Sun or the Earth). If we work in Gaussian units where we measure distances in AU (Astronomical Units); velocity in AU/day; and mass in units of the mass of the Sun, then for an object moving in the gravitational potential of the Sun, <MATH>\mu = 0.00029591220828559115</MATH> and <MATH>m = 1</MATH>. For an object moving in the gravitational field of the Earth, <MATH>m = 1/354710</MATH>. To convert from AU and days to metres and seconds, one needs to know that 1 AU = 149597870700 metres; and that 1 day = 86400 seconds. Unless the mass of the object that we are modelling is changing, it is convenient to set <MATH>m = 1</MATH>.
In this onedimensional system, the object moving in a gravitational potential is a bit limited in terms of directional options. It can either go up, or it can go down. Clearly, as the object moves towards the location of the gravitational source, the force acting on the object is going to increase without limit and the integration is going to going haywire, but so long as we are a reasonable distance away from <MATH>Q^*</MATH> then the integration scheme given above will provide a reasonable description of the object's motion.
In three dimensions
In three dimensions, the integration scheme looks much the same. But now we have to apply it to three spatial coordinates and three velocity components  i.e., the integration scheme becomes one that updates one set of 6 numbers <MATH>\left\{Q_{x,0},Q_{y,0},Q_{z,0},P_{x,0},P_{y,0},P_{z,0}\right\}</MATH> to a new set of six numbers <MATH>\left\{Q_{x,2},Q_{y,2},Q_{z,2},P_{x,2},P_{y,2},P_{z,2}\right\}</MATH>. And the second order integration scheme that does this is as follows:
<MATH>Q_{x,1} \leftarrow Q_{x,0}+\frac{1}{2}\,\delta t\, P_{x,0}</MATH>
<MATH>Q_{y,1} \leftarrow Q_{y,0}+\frac{1}{2}\,\delta t\, P_{y,0}</MATH>
<MATH>Q_{z,1} \leftarrow Q_{z,0}+\frac{1}{2}\,\delta t\, P_{z,0}</MATH>
<MATH>P_{x,2} \leftarrow P_{x,0}+\delta t\, F_x\left(Q_{x,1},Q_{y,1},Q_{z,1}\right)</MATH>
<MATH>P_{y,2} \leftarrow P_{y,0}+\delta t\, F_y\left(Q_{x,1},Q_{y,1},Q_{z,1}\right)</MATH>
<MATH>P_{z,2} \leftarrow P_{z,0}+\delta t\, F_z\left(Q_{x,1},Q_{y,1},Q_{z,1}\right)</MATH>
<MATH>Q_{x,2} \leftarrow Q_{x,1}+\frac{1}{2}\,\delta t\, P_{x,2}</MATH>
<MATH>Q_{y,2} \leftarrow Q_{y,1}+\frac{1}{2}\,\delta t\, P_{y,2}</MATH>
<MATH>Q_{z,2} \leftarrow Q_{z,1}+\frac{1}{2}\,\delta t\, P_{z,2}</MATH>
Each step of the integrator is the same  except now we apply it three times  once for each spatial coordinate.
You should note that we now have three force functions rather than just one. This is because force is really a 'vector' rather than 'scalar' quantity. <MATH>F_x\left(Q_{x},Q_{y},Q_{z}\right)</MATH> is the force acting in the xdirection on the object located at position <MATH>{Q_x, Q_y, Q_z}</MATH>; <MATH>F_y\left(Q_{x},Q_{y},Q_{z}\right)</MATH> is the force acting in the ydirection on the object located at position <MATH>{Q_x, Q_y, Q_z}</MATH>; and <MATH>F_z\left(Q_{x},Q_{y},Q_{z}\right)</MATH> is the force in the xdirection acting on the object located at position <MATH>{Q_x, Q_y, Q_z}</MATH>. And if we are in the gravitational field of a single body, then these functions become:
<MATH>F_x\left(Q_{x},Q_{y},Q_{z}\right) =  \frac{\mu\,m\,\left(Q_{x}Q_x^*\right)}{\left(\left(Q_{x}Q_x^*\right)^2+\left(Q_{y}Q_y^*\right)^2+\left(Q_{z}Q_z^*\right)^2\right)^{3/2}}</MATH>
<MATH>F_y\left(Q_{x},Q_{y},Q_{z}\right) =  \frac{\mu\,m\,\left(Q_{y}Q_y^*\right)}{\left(\left(Q_{x}Q_x^*\right)^2+\left(Q_{y}Q_y^*\right)^2+\left(Q_{z}Q_z^*\right)^2\right)^{3/2}}</MATH>
<MATH>F_z\left(Q_{x},Q_{y},Q_{z}\right) =  \frac{\mu\,m\,\left(Q_{z}Q_z^*\right)}{\left(\left(Q_{x}Q_x^*\right)^2+\left(Q_{y}Q_y^*\right)^2+\left(Q_{z}Q_z^*\right)^2\right)^{3/2}}</MATH>
where <MATH>\left(Q_x^*,Q_y^*,Q_z^*\right)</MATH> are the coordinates of the gravitating body (when the object whose motion you are modelling is at <MATH>\left(Q_x,Q_y,Q_z\right)</MATH>)
And that's about all there is to this this secondorder symplectic integrator. To integrate, one chooses a suitable small stepsize (so that the overall errors are as small as some tolerance that you require for your calculations); one specifies the initial conditions of the object  its position and velocity  and then one repeatedly applies the integration step for as long as you wish.
Now, for those interested, it is worthwhile setting up this integrator in, say, an spreadsheet and seeing how it performs under various sizes of timesteps. and initial conditions. Of course, if there is more than one gravitating body, then you have additional terms in the force functions, but the basic scheme of the updating rule remains the same.
As a sequel to this post, I will sketch how you can quickly build fourth and sixth order symplectic integrators from this simple secondorder integrator.
And as a sequel to that I will talk about the distinction between 'implicit' and 'explicit' integrators.
Precis
This article has a precis and appears in the Random addon or Random article section on the Main Page. The precis can be found at A Symplectic Integrator/precis and is displayed below.
A Symplectic Integrator. An explanation of how to set up a certain type of numeric integrator, a key aspect of celestial mechanics calculations, such as those performed by Orbiter. (More...)
