On An Ethereal Orrery

My fellow students and I have lately been wondering whether we might be able to employ Professor B------'s Experimental Clockwork Mathematical Apparatus[1] to fashion an ethereal orrery, making a model of the heavenly bodies with equations rather than brass.
In particular we have been curious as to whether we might construct such a model using nought but Sir N-----'s law of universal gravitation, which posits that those bodies are attracted to one another with a force that is proportional to the product of their masses divided by the square of the distance between them, and laws of motion, which posit that a body will remain at rest or move with constant velocity if no force acts upon it, that if a force acts upon it then it will be accelerated at a rate proportional to that force divided by its mass in the direction of that force and that it in return exerts a force of equal strength in the opposite direction.

The Two Body Problem

When there are but two bodies under consideration, Sir N----- proved that they will move along conic sections upon a plane, being the curves that are made by the intersections of planes with the surfaces of cones and which may be defined with the formulae
\begin{align*} x(\theta) &= \frac{p e \cos \theta}{1 + e \cos \theta}\\ y(\theta) &= \frac{p e \sin \theta}{1 + e \cos \theta} \end{align*}
where $$p$$ is known as the focal parameter and $$e$$ is known as the eccentricity. Note that these only hold true for positive eccentricities but we may also define such curves for an eccentricity of zero with the formulae
\begin{align*} x(\theta) &= \cos \theta\\ y(\theta) &= \sin \theta \end{align*}
which describe a unit circle centred at the origin and that, in both cases, the coordinates must be translated, rotated and scaled by appropriate constant factors.
We have put together deck 1 to plot points lying upon such conic sections for a focal parameter of one and eccentricities of up to two to illustrate the shape of the paths that Sir N----- deduced that one of a pair of bodies must traverse under the influence of the gravitational force of the other.

Deck 1: Conic Sections

Alas, the general case has proven utterly resistant to mathematical reckoning and we shall consequently be compelled to numerically approximate its behaviour.
Now we shall have no chance of doing so if we cannot well approximate the motion of a body that is subject to an inconstant acceleration, such as will arise from the gravitational forces of several others in motion about it.

Motion Under Acceleration

For a particle moving upon a straight line under an acceleration $$a(t)$$, its velocity $$v(t)$$ and position $$x(t)$$ are related to that acceleration by the derivatives
\begin{align*} v(t) &= \frac{\mathrm{d}x(t)}{\mathrm{d}t}\\ a(t) &= \frac{\mathrm{d}v(t)}{\mathrm{d}t} \end{align*}
or, equivalently, by the integrals
\begin{align*} v(t) &= v(0) + \int_0^t a(u) \, \mathrm{d}u\\ x(t) &= x(0) + \int_0^t v(u) \, \mathrm{d}u \end{align*}
Note that if the acceleration is constant then these yield
\begin{align*} v(t) &= v(0) + \int_0^t a \, \mathrm{d}u = v(0) + \big[a u\big]_0^t = v(0) + at\\ x(t) &= x(0) + \int_0^t v(0) + a u \, \mathrm{d}u = x(0) + \big[v(0) u + \tfrac12 a u^2\big]_0^t = x(0) + v(0) t + \tfrac12 a t^2 \end{align*}
which, assuming that
\begin{align*} x(0) &= 0\\ v(0) &= u \end{align*}
result in the familiar schoolroom formulae
\begin{align*} v(t) &= u + at\\ x(t) &= u t + \tfrac12 a t^2 \end{align*}

A First Approximation

Far and away the simplest manner by which we might approximate such motion is to replace the derivatives in the relationships between the acceleration, velocity and position with finite differences, being their differences at particular times divided by the differences between those times
\begin{align*} v_i &= \frac{\Delta x_i}{\Delta t_i} = \frac{x_{i+1}-x_i}{t_{i+1}-t_i}\\ a_i &= \frac{\Delta v_i}{\Delta t_i} = \frac{v_{i+1}-v_i}{t_{i+1}-t_i} \end{align*}
Rearranging these yields the approximate equations of motion
\begin{align*} x_{i+1} &= x_i + \left(t_{i+1}-t_i\right) v_i\\ v_{i+1} &= v_i + \left(t_{i+1}-t_i\right) a_i \end{align*}
To test the accuracy of these equations my fellow students and I supposed that a particle was moving exponentially, so that
\begin{align*} x(t) &= e^t\\ v(t) &= \frac{\mathrm{d}x(t)}{\mathrm{d}t} = e^t\\ a(t) &= \frac{\mathrm{d}v(t)}{\mathrm{d}t} = e^t \end{align*}
and we may therefore estimate $$a_i$$ with $$x_i$$, echoing the fact that we shall derive the gravitational forces within our ethereal orrery from the realised positions of its bodies.
We have assembled deck 2 to compare the approximate position that those equations yield to the exact position that we had supposed, plotting the former in red and the latter in black for steps in time of one tenth.

Deck 2: The First Approximation

These positions most evidently diverge which is made all the more clear by deck 3 which plots the difference between them.

Deck 3: The First Approximation's Error

In retrospect this inaccuracy shouldn't have come as a particular surprise since the equations are approximating the derivative with the forward finite difference, which is a notoriously poor substitute[2].

A Second Approximation

A markedly better approximation of the derivative is the symmetric finite difference, which is given by
$\frac{\mathrm{d}f(x)}{\mathrm{d}x} \approx \frac{f(x+\delta)-f(x-\delta)}{2 \delta}$
for some appropriate value of $$\delta$$. Using this to approximate the equations of motion yields
\begin{align*} v_i &= \frac{\Delta x_i}{\Delta t_i} = \frac{x_{i+1}-x_{i-1}}{t_{i+1}-t_{i-1}}\\ a_i &= \frac{\Delta v_i}{\Delta t_i} = \frac{v_{i+1}-v_{i-1}}{t_{i+1}-t_{i-1}} \end{align*}
and consequently
\begin{align*} x_{i+1} &= x_{i-1} + \left(t_{i+1}-t_{i-1}\right) v_i\\ v_{i+1} &= v_{i-1} + \left(t_{i+1}-t_{i-1}\right) a_i\\ a_{i+1} &= x_{i+1} \end{align*}
which we have put to test against exponential motion in deck 4.

Deck 4: The Second Approximation

This is a significant improvement upon our first attempt that deck 5 makes plain by plotting its deviation from the true motion at one tenth the scale of deck 3.

Deck 5: The Second Approximation's Error

A Third Approximation

If we assume that the acceleration of the particle is constant between times $$t_i$$ and $$t_{i+1}$$ then we might approximate its motion by employing our schoolroom formulae
\begin{align*} x_{i+1} &= x_i + v_i \left(t_{i+1}-t_i\right) + \tfrac12 a_i \left(t_{i+1}-t_i\right)^2\\ v_{i+1} &= v_i + a_i \left(t_{i+1}-t_i\right)\\ a_{i+1} &= x_{x+1} \end{align*}
which my fellow students and I put to the test with deck 6.

Deck 6: The Third Approximation

Whilst this isn't anything like as accurate as the symmetric finite difference, it's certainly an improvement upon the forward finite difference as evidenced by deck 7 which plots its error at the same scale.

Deck 7: The Third Approximation's Error

A Fourth Approximation

We may improve upon our calculation of the velocity by employing the trapezium rule, by which we might approximate an integral with
$\int_a^b f(x) \, \mathrm{d}x \approx \tfrac12 \left(f(a)+f(b)\right) \times (b-a)$
so that
\begin{align*} x_{i+1} &= x_i + v_i \left(t_{i+1}-t_i\right) + \tfrac12 a_i \left(t_{i+1}-t_i\right)^2\\ v_{i+1} &= v_i + \tfrac12\left(a_i+x_{i+1}\right) \times \left(t_{i+1}-t_i\right)\\ a_{i+1} &= x_{i+1} \end{align*}
Deck 8 compares the results of this approximation for an exponential motion with its correct position

Deck 8: The Fourth Approximation

and deck 9 plots the difference between them at the one tenth scale

Deck 9: The Fourth Approximation's Error

A Fifth Approximation

Now it transpires[3][4] that if we take several such steps with carefully chosen lengths, both forwards and backwards in time, then we can dramatically improve the accuracy of the approximation.
For example, given the coefficients
\begin{align*} c_0 &= \phantom{-}0.784513610477560\\ c_1 &= \phantom{-}0.235573213359357\\ c_2 &= -1.17767998417887\\ c_3 &= \phantom{-}1.31518632068391\\ c_4 &= \phantom{-}c_2\\ c_5 &= \phantom{-}c_1\\ c_6 &= \phantom{-}c_0 \end{align*}
defining the state of the motion at the commencement of those several steps as
\begin{align*} x_{i,0} &= x_i\\ v_{i,0} &= v_i\\ a_{i,0} &= a_i \end{align*}
and its evolution with
\begin{align*} \Delta t_{i,j} &= c_j \Delta t_i\\ x_{i,j+1} &= x_{i,j} + v_{i,j} \Delta t_{i,j} + \tfrac12 a_{i,j} \Delta t_{i,j}^2\\ v_{i,j+1} &= v_{i,j} + \tfrac12 \left(a_{i,j}+x_{i,j+1}\right) \Delta t_{i,j}\\ a_{i,j+1} &= x_{i,j+1} \end{align*}
then the state at the end will be
\begin{align*} x_{i+1} &= x_{i,7}\\ v_{i+1} &= v_{i,7}\\ a_{i+1} &= a_{i,7} \end{align*}
Now I must confess that my fellow students and I could not fathom why those particular coefficients should prove so effective, but deck 10 certainly lends support to claims that they do with its having no discernible deviation of the approximate from the true position.

Deck 10: The Fifth Approximation

To show the extent of the deviation it was necessary to plot it at a scale of one millionth that of deck 3.

Deck 11: The Fifth Approximation's Error

An important feature of this scheme is that its errors are proportional to the sixth power of $$\Delta t_i$$ so that if we halve it then the error should be reduced by a factor of roughly sixty four. Deck 12 demonstrates that this holds for our exponential motion, plotting the results at one hundredth the scale of deck 11.

Deck 12: Halving The Step Lengths

Whilst there is no doubt in our minds that greater effort should yield yet more accurate results, we believe that this shall be sufficient for our needs and plan to use it as the foundation of our orrery.
$$\Box$$

References

[1] On An Age Of Wonders, www.thusspakeak.com, 2014.

[2] You're Going To Have To Think! Why Finite Differences Won't Cure Your Calculus Blues, www.thusspakeak.com, 2013.

[3] Yoshida, H., Construction of higher order symplectic integrators, Physics Letters A, Vol. 150, No. 5, 6, 7, Elsevier, 1990.

[4] Hut, P. & Makino, J., The 2-Body Problem: Higher-Order Integrators, www.artcompsci.org, 2007.