On the Extension of Adams--Bashforth--Moulton Methods for Numerical Integration of Delay Differential Equations and Application to the Moon's Orbit

TL;DR

This paper extends Adams--Bashforth--Moulton methods to enable accurate forward and backward numerical integration of delay differential equations, with applications to modeling the Moon's orbit considering tidal delays.

Contribution

It introduces a modified Adams--Bashforth--Moulton method capable of integrating DDEs both forwards and backwards in time, addressing a key challenge in celestial mechanics.

Findings

Successful backward and forward integration of Moon's DDEs demonstrated

Enhanced numerical methods for delay differential equations developed

Application to lunar orbit modeling with improved accuracy

Abstract

One of the problems arising in modern celestial mechanics is the need of precise numerical integration of dynamical equations of motion of the Moon. The action of tidal forces is modeled with a time delay and the motion of the Moon is therefore described by a functional differential equation (FDE) called delay differential equation (DDE). Numerical integration of the orbit is normally being performed in both directions (forwards and backwards in time) starting from some epoch (moment in time). While the theory of normal forwards-in-time numerical integration of DDEs is developed and well-known, integrating a DDE backwards in time is equivalent to solving a different kind of FDE called advanced differential equation, where the derivative of the function depends on not yet known future states of the function. We examine a modification of Adams--Bashforth--Moulton method allowing to perform integration of the Moon's DDE forwards and backwards in time and the results of such integration.

Equations

$$\dot{\omega }={\left(\frac{I}{m}\right)}^{-1}\left[\frac{\mathbf{N}}{m}-\frac{\dot{I}}{m}\omega -\omega \times \left(\frac{I}{m}\omega \right)\right],$$

$$\frac{I}{m}$$

$$\displaystyle+k_{2}\frac{R_{\mathrm{M}}^{5}}{3\mu_{\mathrm{M}}}\left[\begin{array}[]{ccc}\omega_{x}^{2}-\frac{1}{3}(\omega^{2}-n^{2})&\omega_{x}\omega_{y}&\omega_{x}\omega_{z}\\ \omega_{x}\omega_{y}&\omega_{y}^{2}-\frac{1}{3}(\omega^{2}-n^{2})&\omega_{y}\omega_{z}\\ \omega_{x}\omega_{z}&\omega_{y}\omega_{z}&\omega_{z}^{2}-\frac{1}{3}(\omega^{2}+2n^{2})\\ \end{array}\right],$$

$$\dot{\mathbf{x}}(t)=\mathbf{f}(t,\mathbf{x}(t),\mathbf{x}(t-\tau ),\dot{\mathbf{x}}(t-\tau )),$$

$$\mathbf{x}(t_0)={\mathbf{x}}_0$$

$$\dot{\mathbf{x}}(t)=\mathbf{f}(t,\mathbf{x}(t),\mathbf{x}(t+\tau ),\dot{\mathbf{x}}(t+\tau )).$$

$$\{\begin{array}{l}\dot{\mathbf{x}}(t)=\mathbf{f}(t,\mathbf{x}(t),\mathbf{x}(t\pm \tau ),\dot{\mathbf{x}}(t\pm \tau )),\\ \mathbf{x}(t_0)={\mathbf{x}}_0.\end{array}$$

$$\dot{\mathbf{y}}(t)=\mathbf{g}(t,\mathbf{y}(t))=\mathbf{f}(t,\mathbf{y}(t),\mathbf{y}(t),0)$$

$$\{\begin{array}{l}\dot{\mathbf{y}}(t)=\mathbf{g}(t,\mathbf{y}(t)),\\ \mathbf{y}(t_n)={\mathbf{x}}_n.\end{array}$$

$${\mathbf{x}}^{(p-1)}(t)=\frac{\sum _{j=0}^p\frac{w_j}{t-t_j}{\mathbf{x}}_j}{\sum _{j=0}^p\frac{w_j}{t-t_j}},$$

$$w_j=\frac{1}{{\prod }_{k\ne j}(t_j-t_k)},j=0,...,p.$$

$${\mathbf{x}}_{n+1}^{(\mathrm{p})}={\mathbf{x}}_n+h\sum _{j=0}^{k-1}{\gamma }_j{\nabla }^j{\mathbf{f}}_n,$$

$${\mathbf{x}}_{n+1}={\mathbf{x}}_n+h\sum _{j=0}^k{c}_j{\nabla }^j{\mathbf{f}}_{n+1},$$

$$c_j=-\sum _{i=0}^{j-1}\frac{c_i}{j+1-i},$$

$${\gamma }_j=\sum _{k=0}^j{c}_k.$$

$${\mathbf{x}}_{n+1}={\mathbf{x}}_{n+1}^{(\mathrm{p})}+h\left(\sum _{j=0}^k{c}_j{\nabla }^j{\mathbf{f}}_{n+1}-\sum _{j=0}^{k-1}{\gamma }_j{\nabla }^j{\mathbf{f}}_n\right).$$

$$\mathcal{S}=\sum _{j=0}^k{c}_j{\nabla }^j{\mathbf{f}}_{n+1}-\sum _{j=0}^{k-1}{\gamma }_j{\nabla }^j{\mathbf{f}}_n.$$

$$\begin{array}{rl}\mathcal{S}& =\sum _{j=0}^k({\gamma }_j-{\gamma }_{j-1}){\nabla }^j{\mathbf{f}}_{n+1}-\sum _{j=0}^{k-1}{\gamma }_j{\nabla }^j{\mathbf{f}}_n=\\ & ={\gamma }_k{\nabla }^k{\mathbf{f}}_{n+1}+\sum _{j=0}^{k-1}{\gamma }_j\left({\nabla }^j{\mathbf{f}}_{n+1}-{\nabla }^j{\mathbf{f}}_n\right)-\sum _{j=0}^k{\gamma }_{j-1}{\nabla }^j{\mathbf{f}}_{n+1}.\end{array}$$

$$\begin{array}{rl}\mathcal{S}& ={\gamma }_k{\nabla }^k{\mathbf{f}}_{n+1}+\sum _{j=0}^{k-1}{\gamma }_j{\nabla }^{j+1}{\mathbf{f}}_{n+1}-\sum _{j=1}^k{\gamma }_{j-1}{\nabla }^j{\mathbf{f}}_{n+1}=\\ & ={\gamma }_k{\nabla }^k{\mathbf{f}}_{n+1}+\sum _{j=0}^{k-1}{\gamma }_j{\nabla }^{j+1}{\mathbf{f}}_{n+1}-\sum _{i=0}^{k-1}{\gamma }_i{\nabla }^{i+1}{\mathbf{f}}_{n+1}=\\ & ={\gamma }_k{\nabla }^k{\mathbf{f}}_{n+1}.\end{array}$$

$${\mathbf{x}}_{n+1}={\mathbf{x}}_n^{(\mathrm{p})}+h{\gamma }_k{\nabla }^k{\mathbf{f}}_{n+1},$$

$$\nabla {\mathbf{f}}_i$$

$${\nabla }^k{\mathbf{f}}_i$$

Full text

Numerical integration of the orbit is normally being performed in both directions (forwards and backwards in time) starting from some epoch (moment in time). While the theory of normal forwards-in-time numerical integration of DDEs is developed and well-known, integrating a DDE backwards in time is equivalent to solving a different kind of FDE called advanced differential equation, where the derivative of the function depends on not yet known future states of the function.

We examine a modification of Adams–Bashforth–Moulton method allowing to perform integration of the Moon’s DDE forwards and backwards in time and the results of such integration.

1 Introduction

In 1960’s, two things happened simultaneously: the space era began and computers made their way into research labs. In years that followed, requirements to accuracy of the coordinates of solar system bodies (ephemeris) for space missions grew. Eventually, the so-called “analytical” ephemeris, calculated via sophisticated series expansion, gave way to more precise numerical ephemeris, obtained by numerical integration of dynamical equations on a computer.

As the precision of astronomical observations improved, the dynamical models underlying the ephemeris became more complex over time. The laser ranging of the lunar retroreflectors, which began in the end of 1969, allowed to determine the presence of a liquid core inside the Moon three decades later [1]. Such findings require careful analysis of different causes of the observed motion; specifically in the case of the Moon, it is important to model the dissipation of energy that is caused by tides and the rotation of the elastic Moon.

Redistribution of mass in a large rocky object does not happen instantly, but rather takes hours to happen. Thus, the corresponding tidal terms in the dynamical equations are naturally written with a time delay.

Let us move on to the equation for which the method descibed in this article was developed. This is the equation for the rotational motion of the Moon [2, 3], written in a Moon-fixed frame:

[TABLE]

where $\bm{\mathbf{\omega}}$ is the angular velocity, $m$ is the mass of the Moon, $\bm{\mathbf{N}}$ is the external torque, and $I(t)$ is the inertia tensor of the lunar mantle. The inertia tensor is defined as

[TABLE]

where $I_{0}$ is the inertia tensor of the undistorted lunar mantle, $\mu_{\mathrm{E}}$ and $\mu_{\mathrm{M}}$ are the gravitational parameters of Earth and the Moon respectively, $R_{\mathrm{M}}$ is the equatorial radius of the Moon, $k_{2}$ is the degree-2 Love number of the Moon, $\bm{r}=(x,y,z)^{T}$ is the position of the Moon relative to Earth, $n$ is the lunar mean motion.

Most importantly, the inertia tensor of the Moon is subject to delayed tidal distortion from Earth and delayed spin distortion. Tidal and spin distortions are evaluated with a delayed argument, so calculation of $I(t)$ involves not $\bm{r}(t)$ and $\bm{\omega}(t)$, but $\bm{r}(t-\tau)$ and $\bm{\omega}(t-\tau)$, and $\tau$ is a constant time delay (in [3], $\tau=0.096\ \mathrm{d}$). Thus, the equation 1 is a delay differential equation of neutral type with constant delays [4], the general form for which would be

[TABLE]

where $\bm{\mathbf{x}}=(x_{1},x_{2},\dots,x_{n})$ and $\tau>0$.

The initial condition for numerical integration of the Moon’s orbit is of the form

[TABLE]

and is defined for $t_{0}=\mathrm{JD}\ 2446000.5$ (1984-10-27 midnight) [3]. To cover the whole interval of the lunar laser ranging observations, numerical integration has to be performed both forward and backward in time.

It is obvious that backward-in-time integration of DDE (3) is equivalent to forward-in-time integration of an advanced differential equation

[TABLE]

Therefore, what must be invented is a numerical method able to solve an initial value problem

[TABLE]

In other words, the method must be able to numerically integrate both delay and advanced differential equations given an initial condition (4).

2 Known approaches

Most of existing numerical integrators for DDEs, like [5, 6], work only in forward direction in time, and require a so-called “history function” that unambiguously defines the $\bm{\mathbf{x}}(t-\tau)$ for $t\in[0,\tau)$. That requirement conflicts with the nature of motion of celestial bodies, which is usually quasiperiodic and does not “take off” from some initially known trajectory.

Particularly for eq. (1), some dedicated methods are known that approximate the $\bm{\mathbf{x}}(t\pm\tau)$ and $\dot{\bm{\mathbf{x}}}(t\pm\tau)$. In [7], a quadratic approximation (based on known derivatives of the equations) is used to calculate the delayed terms. In [2], the approximation is done via high degree polynomials that are pre-fit before the main integration [8].

3 The “nested integration” method

Suppose that for an arbitrary point $t_{n}=t_{0}+nh$ the function value $\bm{\mathbf{x}}_{n}=\bm{\mathbf{x}}(t_{n})$ is known. Calculation of $\bm{\mathbf{f}}_{n}=\bm{\mathbf{f}}(t_{n},\,\bm{\mathbf{x}}_{n},\,\bm{\mathbf{x}}(t_{n}\pm\tau),\,\dot{\bm{\mathbf{x}}}(t_{n}\pm\tau))$ then requires knowledge of delayed function and derivative values $\bm{\mathbf{x}}(t_{n}\pm\tau)$ and $\dot{\bm{\mathbf{x}}}(t_{n}\pm\tau)$.

The “nested integration” method comes down to introducing a new equation

[TABLE]

and then finding an approximation to the delayed function value $\bm{\mathbf{x}}(t_{n}\pm\tau)$ as the numerical solution $\bm{\mathbf{y}}(t_{n}\pm\tau)$ of the initial value problem involving eq. (7)

[TABLE]

The delayed derivative $\dot{\bm{\mathbf{x}}}(t_{n}\pm\tau)$ in this case is simply calculated as $\bm{\mathbf{g}}(t_{n}\pm\tau,\,\bm{\mathbf{y}}(t_{n}\pm\tau))$. The value for $\dot{\bm{\mathbf{x}}}(t_{n}-\tau)$ in eq. 7 is chosen arbitrarily.

Both outer (of problem (6)) and inner (of problems (8) occurring at each outer integration step) integration can be performed by any numerical method, single- or multistep.

4 The interpolation method

Another way to calculate the delayed values of $\bm{\mathbf{x}}$ and $\dot{\bm{\mathbf{x}}}$ is to use polynomial interpolation of data computed on previous integration steps. If the function values at nodes $[t_{0},...,t_{p}]$ are $[\bm{\mathbf{x}}_{0},...,\bm{\mathbf{x}}_{p}]$, then for arbitrary point $t$ the value $\bm{\mathbf{x}}^{(p-1)}(t)$ of interpolating polynomial of degree $(p-1)$ is given by the barycentric Lagrange formula [9]

[TABLE]

where $w_{j}$ are the barycentric weights defined by

[TABLE]

The weights $\{w_{j}\}$ depend only on values of $\{t_{j}\}$, or, for an equidistant grid with $h=t_{n}-t_{n-1}$, on step size $h$. Calculating the weights with (10) requires $O(p^{2})$ operations and has to be done once for one step size $h$. Once the weights are known, evaluation of the Lagrange polynomial is performed with (9) in $O(p\cdot n)$ operations.

5 The ABMD Integrator

The ABMD integrator, which is currently used in ERA-8 software system [10], uses the Adams–Bashforth–Moulton predictor–corrector methods [11].

The Adams–Bashforth predictor and the Adams–Moulton corrector are given by the formulas:

[TABLE]

[TABLE]

where $k$ is the method order, $\bm{\mathbf{x}}^{\mathrm{(p)}}$ is the predicted solution value, and the coefficients $\gamma_{j}$ and $c_{j}$ are defined as

[TABLE]

[TABLE]

The classic Adams–Moulton corrector formula (12) can be rewritten in a more computationally efficient form. First, we will extract $\bm{\mathbf{x}}_{n}$ from (11) and insert it into 12:

[TABLE]

For convenience, let’s denote

[TABLE]

Now, substituting $c_{j}$ with $\gamma_{j}-\gamma_{j-1}$, taking out the last term of the first sum and regrouping the sums:

[TABLE]

As $\nabla^{j}\bm{\mathbf{f}}_{n+1}-\nabla^{j}\bm{\mathbf{f}}_{n}=\nabla^{j+1}\bm{\mathbf{f}}_{n+1}$ and $\gamma_{-1}=0$,

[TABLE]

Hence, eq. 12 simplifies to the modified Adams–Moulton corrector formula

[TABLE]

which, provided that the difference $\nabla^{k}\bm{\mathbf{f}}_{n+1}$ is known, requires $O(n)$ operations instead of $O(kn)$ required by the formula (12).

So, the three-stage predictor–evaluator–corrector (PEC) algorithm comes down to calculating a predicted solution value $\bm{\mathbf{x}}_{n+1}^{\mathrm{(p)}}$ with eq. (11), evaluating the right-hand side $\bm{\mathbf{f}}_{n+1}=\bm{\mathbf{f}}\left(t_{n+1},\ \bm{\mathbf{x}}_{n+1}^{\mathrm{(p)}}\right)$, and then correcting the predicted solution value with eq. (19). The evaluator and corrector steps may be repeated multiple times: such methods would be referred to as PECE, PECEC, PECECE etc. The ABMD integrator uses a PECEC scheme.

By definition of the backward difference operator $\nabla$,

[TABLE]

Therefore, calculating $k$ differences (from $\nabla\bm{\mathbf{f}}$ to $\nabla^{k}\bm{\mathbf{f}}$) requires $O(k^{2}n)$ operations. If, however, differences $\nabla\bm{\mathbf{f}}_{i}\dots\nabla^{k}\bm{\mathbf{f}}_{i}$ are known, then successively updating them to $\nabla\bm{\mathbf{f}}_{i+1}\dots\nabla^{k}\bm{\mathbf{f}}_{i+1}$ with the definition (20) requires only $O(kn)$ operations, which makes it the preferred method of calculating backward differences during integration, especially for problems with large values of $n$.

As follows from eqs. 11 and 19, to compute a new solution value, a method of order $k$ requires knowledge of $k$ previous RHS values. At the start of the integration only one solution value is known, thus, the other $k-1$ values have to be found with a startup procedure using a single-step method.

In the ABMD integrator, startup procedure uses the nested integration method. After the startup is finished, delayed solution values are retrieved with Lagrange interpolation (9), which turns to extrapolation when integration is performed backward in time.

6 Results

We did not find DDEs that are sensible and have known solutions or invariants to check against. So the primary criterion to validate the proposed method will be the numerical integration of the “real world” eqs. 1 and 2.

We used the ABMD integrator with Adams–Bashforth–Moulton PECEC scheme of order 13 with step size and with step size $h_{\mathrm{ABM}}=1/16\ \mathrm{d}$.

The startup uses the 8th order Dormand–Prince method [11] for the outer integration method and the 4th order Runge–Kutta for the inner integration. To match the precision of the 13th order ABM, step size for Dormand–Prince startup is reduced: $h_{\mathrm{DP}}=h_{\mathrm{ABM}}/8$.

Two tests were done. The first one is the usual forward–backward test. The full dynamical equations of the Moon from [3] were integrated from 1984 to 2014, and the solution at the last point was taken as the initial condition for integrating backwards, from 2014 to 1984. Differences between the two results are shown in fig. 1. It is evident that the inaccuracy inherent to the method is sufficient at present level of best lunar laser ranging observations (2 mm in Earth–Moon distance).

Another test was done to compare the delayed values $\bm{\mathbf{r}}(t-\tau)$ that are obtained during integration by interpolation with the same values of the final integrated orbit. Fig. 2 shows that values calculated with polynomial interpolation—forward from the epoch—are two orders of magnitude closer to the final orbit than those that are calculated with nested Runge–Kutta integration. For the values calculated with extrapolation—backward from the epoch—the improvement is one order.

7 Conclusion and further work

A modification of the Adams–Bashforth–Moulton multistep numerical integration method was proposed, and its implementation, ABMD, was developed that allows to integrate delay differential equations. Unlike most other methods for DDEs, the developed method does not require a “history function” and can integrate forwards and backwards in time.

There has been no theoretical analysis of the method’s correction or convergence. Most probably the applicability of ABMD is limited to relatively stable systems in which the delay term is not the main term of the equation. It has been shown that ABMD works well for one particular equation of the Earth–Moon system.

The proposed method is trivially generalised to DDEs with multiple constant delays with possibly different signs. The present model of the Earth–Moon system already includes delays in other equations besides the one that was described in this paper: these are equations for the Earth tides [12].

For spacecraft dynamics (as opposed to dynamics of celestial bodies), more precise and complicated tidal models are used, see e.g. [13] for lunar tides. Such models do not contain the delays directly in the equations, but instead use empirical Fourier series whose arguments are linear combinations of Delaunay variables linked to the state of the Earth–Moon–Sun system.

We believe that the new tool developed for integration of arbitrary DDEs can help to build more natural models for tides on Earth and the Moon. Such new models could be a step to a more adequate description of the mechanical processes that are going at the interior of the celestial bodies.

8 Acknowledgements

Authors are grateful to Sergey Kurdubov (IAA RAS) and John Chandler (Harvard CfA) for useful discussions related to this work. Steve Moshier’s implementation of Adams–Bashforth–Moulton method (http://moshier.net/ssystem.html) was a helpful reference.