In an earlier post, I showed how to build a 2nd order symplectic integrator. To recap, we can think of the symplectic integrator as a small 'machine' that takes 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]
For details on this, and of the force function [MATH]\left\{F_{x},F_{y},F_{z}\right\}[/MATH], see the earlier post (http://www.orbiter-forum.com/showthread.php?t=35441).
Although 'only' second-order accurate, this is still quite a serviceable integrator - and very general (so long as one is using Cartesian coordinates). And with this second-order integrator, it is also possible to quickly build a fourth-order. Here's how:
Step 1. Let's introduce a short-hand notation such that [MATH]\left\{ \mathbf{Q}_{i},\mathbf{P_{i}}\right\} =\left\{ Q_{x,i},Q_{y,i},Q_{z,i},P_{x,i},P_{y,i},P_{z,i}\right\} [/MATH]. Then symbolically, we can represent the second-order symplectic integrator as a function which we'll call [MATH]G_2[/MATH] such that:
[MATH]\left\{ \mathbf{Q}_{2},\mathbf{P_{2}}\right\} \leftarrow G_{2}\left(\left\{ \mathbf{Q}_{0},\mathbf{P_{0}}\right\} ;\delta t\right)[/MATH]
[MATH]G_2[/MATH], then, is just the nine procedural steps given above which we have grouped together and called [MATH]G[/MATH]. The subscript '2' is there to remind us that collectively these steps represent a second order integrator. On the right-hand side, there is reference to [MATH]\delta t[/MATH]. This is to remind us that the new position [MATH]\left\{ \mathbf{Q}_{2},\mathbf{P_{2}}\right\}[/MATH] generated by applying [MATH]G_2[/MATH] is a function of the step-size [MATH]\delta t[/MATH]. Change the step-size, and we will get a different result.
Step 2. Now, with this preamble, we can define the fouth-order integrator which, for convenience, we will call [MATH]G_4[/MATH]. [MATH]G_4[/MATH] consists of just three steps:
[MATH]\left\{ \mathbf{Q}_{2},\mathbf{P_{2}}\right\} \leftarrow G_{2}\left(\left\{ \mathbf{Q}_{0},\mathbf{P_{0}}\right\} ;\omega_0\,\delta t\right)[/MATH]
[MATH]\left\{ \mathbf{Q}_{4},\mathbf{P_{4}}\right\} \leftarrow G_{2}\left(\left\{ \mathbf{Q}_{2},\mathbf{P_{2}}\right\} ;\omega_1\,\delta t\right)[/MATH]
[MATH]\left\{ \mathbf{Q}_{6},\mathbf{P_{6}}\right\} \leftarrow G_{2}\left(\left\{ \mathbf{Q}_{4},\mathbf{P_{4}}\right\} ;\omega_0\,\delta t\right)[/MATH]
where
[MATH]\omega_0 = \frac{1}{2-2^{1/3}}[/MATH]
[MATH]\omega_1 = -\frac{2^{1/3}}{2-2^{1/3}}[/MATH]
In other words, the fourth-order integrator, is just the second-order integrator applied three times, but where we have used the same time-step for the first and last step; and a different time-step for the middle step. Varying time-steps in this way leads to a magical cancellation of terms so that the end result is fourth order accurate. If you add up the step-size of the three steps, we still have a total overall step-size of just [MATH]\delta t[/MATH].
More succinctly, the fourth-order integrator can be written:
[MATH]\left\{ \mathbf{Q}_{6},\mathbf{P_{6}}\right\} \leftarrow G_{4}\left(\left\{ \mathbf{Q}_{0},\mathbf{P_{0}}\right\} ;\delta t\right)[/MATH]
Simple!
Once you have built the second-order integrator, it becomes extremely easy to build this fourth order-integrator. It is also possible to combined three of these fourth-order accurate steps to produce a 6th order accurate result. And similarly for an 8th order integrator. However, above fourth-order, this becomes an inefficient way of building higher-order integrators, so for the time being we'll stop at the fourth order. In another post, I'll show how to build a better a sixth-order symplectic integrator.
In any event, the fourth-order symplectic integrator described above is an effective tool. With a few tweaks to better control numerical rounding error, it was the main integrator I used to do the test of Orbiter's integration engine. And that seemed to work pretty well. (See http://www.orbiter-forum.com/showthread.php?t=35243)
Explicit versus Implicit integrators
Both the second-order and fourth-order integrators described above are called 'explicit' integrators. If you look closely at the algorithm you will find that aside from a need to recalculate the force function at a new location from time to time, the previous steps of the algorithm give you just enough information to allow you to calculate the next line. In other words, you can start at the first line of the algorithm and calculate the intermediate result. Then you can go on to the second line and calculate a new intermediate result. And then you move onto the third line and so on until you reach the end of the algorithm.
However, there is another class of integrators, known as 'implicit' integrators, where this step-by-step evaluation is not possible. Why? Because the algorithm is usually written in such a way that enumeration of the first step in the algorithm requires information from what will be the last step in the algorithm - i.e, the problem description is circular. (To know A, I need to know B. To know B, I need to know C. And to know C, I need to know A). Many symplectic Runge Kutta integrators are of this kind. The solution to this circularity is to use an algorithm such as Newton-Raphson to find the solution. Alternatively, one can try a simple fixed-point integrator. Neither method is guaranteed to always find a solution. The circularity means that finding a solution is slower and less reliable. But if/when a solution is found, it is typically more accurate than using an explicit integrator of the same order. For situations where speed and reliability is paramount, explicit integrators have the upper-hand.
[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]
For details on this, and of the force function [MATH]\left\{F_{x},F_{y},F_{z}\right\}[/MATH], see the earlier post (http://www.orbiter-forum.com/showthread.php?t=35441).
Although 'only' second-order accurate, this is still quite a serviceable integrator - and very general (so long as one is using Cartesian coordinates). And with this second-order integrator, it is also possible to quickly build a fourth-order. Here's how:
Step 1. Let's introduce a short-hand notation such that [MATH]\left\{ \mathbf{Q}_{i},\mathbf{P_{i}}\right\} =\left\{ Q_{x,i},Q_{y,i},Q_{z,i},P_{x,i},P_{y,i},P_{z,i}\right\} [/MATH]. Then symbolically, we can represent the second-order symplectic integrator as a function which we'll call [MATH]G_2[/MATH] such that:
[MATH]\left\{ \mathbf{Q}_{2},\mathbf{P_{2}}\right\} \leftarrow G_{2}\left(\left\{ \mathbf{Q}_{0},\mathbf{P_{0}}\right\} ;\delta t\right)[/MATH]
[MATH]G_2[/MATH], then, is just the nine procedural steps given above which we have grouped together and called [MATH]G[/MATH]. The subscript '2' is there to remind us that collectively these steps represent a second order integrator. On the right-hand side, there is reference to [MATH]\delta t[/MATH]. This is to remind us that the new position [MATH]\left\{ \mathbf{Q}_{2},\mathbf{P_{2}}\right\}[/MATH] generated by applying [MATH]G_2[/MATH] is a function of the step-size [MATH]\delta t[/MATH]. Change the step-size, and we will get a different result.
Step 2. Now, with this preamble, we can define the fouth-order integrator which, for convenience, we will call [MATH]G_4[/MATH]. [MATH]G_4[/MATH] consists of just three steps:
[MATH]\left\{ \mathbf{Q}_{2},\mathbf{P_{2}}\right\} \leftarrow G_{2}\left(\left\{ \mathbf{Q}_{0},\mathbf{P_{0}}\right\} ;\omega_0\,\delta t\right)[/MATH]
[MATH]\left\{ \mathbf{Q}_{4},\mathbf{P_{4}}\right\} \leftarrow G_{2}\left(\left\{ \mathbf{Q}_{2},\mathbf{P_{2}}\right\} ;\omega_1\,\delta t\right)[/MATH]
[MATH]\left\{ \mathbf{Q}_{6},\mathbf{P_{6}}\right\} \leftarrow G_{2}\left(\left\{ \mathbf{Q}_{4},\mathbf{P_{4}}\right\} ;\omega_0\,\delta t\right)[/MATH]
where
[MATH]\omega_0 = \frac{1}{2-2^{1/3}}[/MATH]
[MATH]\omega_1 = -\frac{2^{1/3}}{2-2^{1/3}}[/MATH]
In other words, the fourth-order integrator, is just the second-order integrator applied three times, but where we have used the same time-step for the first and last step; and a different time-step for the middle step. Varying time-steps in this way leads to a magical cancellation of terms so that the end result is fourth order accurate. If you add up the step-size of the three steps, we still have a total overall step-size of just [MATH]\delta t[/MATH].
More succinctly, the fourth-order integrator can be written:
[MATH]\left\{ \mathbf{Q}_{6},\mathbf{P_{6}}\right\} \leftarrow G_{4}\left(\left\{ \mathbf{Q}_{0},\mathbf{P_{0}}\right\} ;\delta t\right)[/MATH]
Simple!
Once you have built the second-order integrator, it becomes extremely easy to build this fourth order-integrator. It is also possible to combined three of these fourth-order accurate steps to produce a 6th order accurate result. And similarly for an 8th order integrator. However, above fourth-order, this becomes an inefficient way of building higher-order integrators, so for the time being we'll stop at the fourth order. In another post, I'll show how to build a better a sixth-order symplectic integrator.
In any event, the fourth-order symplectic integrator described above is an effective tool. With a few tweaks to better control numerical rounding error, it was the main integrator I used to do the test of Orbiter's integration engine. And that seemed to work pretty well. (See http://www.orbiter-forum.com/showthread.php?t=35243)
Explicit versus Implicit integrators
Both the second-order and fourth-order integrators described above are called 'explicit' integrators. If you look closely at the algorithm you will find that aside from a need to recalculate the force function at a new location from time to time, the previous steps of the algorithm give you just enough information to allow you to calculate the next line. In other words, you can start at the first line of the algorithm and calculate the intermediate result. Then you can go on to the second line and calculate a new intermediate result. And then you move onto the third line and so on until you reach the end of the algorithm.
However, there is another class of integrators, known as 'implicit' integrators, where this step-by-step evaluation is not possible. Why? Because the algorithm is usually written in such a way that enumeration of the first step in the algorithm requires information from what will be the last step in the algorithm - i.e, the problem description is circular. (To know A, I need to know B. To know B, I need to know C. And to know C, I need to know A). Many symplectic Runge Kutta integrators are of this kind. The solution to this circularity is to use an algorithm such as Newton-Raphson to find the solution. Alternatively, one can try a simple fixed-point integrator. Neither method is guaranteed to always find a solution. The circularity means that finding a solution is slower and less reliable. But if/when a solution is found, it is typically more accurate than using an explicit integrator of the same order. For situations where speed and reliability is paramount, explicit integrators have the upper-hand.
Last edited:
