Editing User:MikeB/Symplectic Numeric Integration



== 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 second-order symplectic integrator ==
First, let's focus on a basic second-order symplectic integrator.  It is worthwhile spending a little time on this because, so long as the time-step of the integration is small enough, it actually works pretty well.  It is also an integrator from which it is possible to build higher-order 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 long-run 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 second-order because it has been constructed so as to 'kill' error terms up to the second-order only.

So, what is this integrator?  If we just focus on a one-dimensional 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 one-dimensional example, you can think of 'Q' as representing the x-coordinate 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 x-coordinate of the momentum of a particle.  But since momentum is just mass * velocity (in this coordinate system), we can think of 'P' as the x-coordinate 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 state-vector 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 one-dimension.)

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 one-dimensional 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}{(Q-Q^*)^{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 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 one-dimensional 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 hay-wire, 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 not 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 x-direction 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 y-direction 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 x-direction 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 second-order symplectic integrator.  To integrate, one chooses a suitable small step-size (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 time-steps. 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 second-order integrator.

And as a sequel to that I will talk about the distinction between 'implicit' and 'explicit' integrators.