Cowell's formulation

Cowell's method. Forces from all perturbing bodies (black and gray) are summed to form the total force on body ${\displaystyle \ i\ }$ (red), and this is numerically integrated starting from the initial position (the epoch of osculation).

Cowell's formulation (so named for Philip H. Cowell, who, with A.C.D. Cromellin, used a similar method to predict the return of Halley's comet) is perhaps the simplest of the special perturbation methods.^[8] In a system of ${\displaystyle \ n\ }$ mutually interacting bodies, this method mathematically solves for the Newtonian forces on body ${\displaystyle \ i\ }$ by summing the individual interactions from the other ${\displaystyle j}$ bodies:

${\displaystyle \mathbf {\ddot {r}} _{i}=\sum _{\underset {j\neq i}{j=1}}^{n}\ G\ m_{j}{\frac {\ (\mathbf {r} _{j}-\mathbf {r} _{i})\ }{\ \|\mathbf {r} _{j}-\mathbf {r} _{i}\|^{3}}}}$

where ${\displaystyle \ \mathbf {\ddot {r}} _{i}\ }$ is the acceleration vector of body ${\displaystyle i}$, ${\displaystyle G}$ is the gravitational constant, ${\displaystyle \ m_{j}\ }$ is the mass of body ${\displaystyle j}$, ${\displaystyle \ \mathbf {r} _{i}\ }$ and ${\displaystyle \ \mathbf {r} _{j}\ }$ are the position vectors of objects ${\displaystyle \ i\ }$ and ${\displaystyle \ j\ }$ respectively, and ${\displaystyle \ r_{ij}\equiv \|\mathbf {r} _{j}-\mathbf {r} _{i}\|\ }$ is the distance from object ${\displaystyle i}$ to object ${\displaystyle \ j\ }$, all vectors being referred to the barycenter of the system. This equation is resolved into components in ${\displaystyle \ x\ ,}$ ${\displaystyle \ y\ ,}$ and ${\displaystyle \ z\ ,}$ and these are integrated numerically to form the new velocity and position vectors. This process is repeated as many times as necessary. The advantage of Cowell's method is ease of application and programming. A disadvantage is that when perturbations become large in magnitude (as when an object makes a close approach to another) the errors of the method also become large.^[9] However, for many problems in celestial mechanics, this is never the case. Another disadvantage is that in systems with a dominant central body, such as the Sun, it is necessary to carry many significant digits in the arithmetic because of the large difference in the forces of the central body and the perturbing bodies, although with high precision numbers built into modern computers this is not as much of a limitation as it once was.^[10]

Encke's method

Encke's method. Greatly exaggerated here, the small difference δr (blue) between the osculating, unperturbed orbit (black) and the perturbed orbit (red), is numerically integrated starting from the initial position (the epoch of osculation).

Encke's method begins with the osculating orbit as a reference and integrates numerically to solve for the variation from the reference as a function of time.^[11] Its advantages are that perturbations are generally small in magnitude, so the integration can proceed in larger steps (with resulting lesser errors), and the method is much less affected by extreme perturbations. Its disadvantage is complexity; it cannot be used indefinitely without occasionally updating the osculating orbit and continuing from there, a process known as rectification.^[9] Encke's method is similar to the general perturbation method of variation of the elements, except the rectification is performed at discrete intervals rather than continuously.^[12]

Letting ${\displaystyle {\boldsymbol {\rho }}}$ be the radius vector of the osculating orbit, ${\displaystyle \mathbf {r} }$ the radius vector of the perturbed orbit, and ${\displaystyle \delta \mathbf {r} }$ the variation from the osculating orbit,

${\displaystyle \delta \mathbf {r} =\mathbf {r} -{\boldsymbol {\rho }}}$, and the equation of motion of ${\displaystyle \delta \mathbf {r} }$ is simply

(1)

${\displaystyle {\ddot {\delta \mathbf {r} }}=\mathbf {\ddot {r}} -{\boldsymbol {\ddot {\rho }}}}$.

(2)

${\displaystyle \mathbf {\ddot {r}} }$ and ${\displaystyle {\boldsymbol {\ddot {\rho }}}}$ are just the equations of motion of ${\displaystyle \mathbf {r} }$ and ${\displaystyle {\boldsymbol {\rho }},}$

${\displaystyle \mathbf {\ddot {r}} =\mathbf {a} _{\text{per}}-{\mu \over r^{3}}\mathbf {r} }$ for the perturbed orbit and

(3)

${\displaystyle {\boldsymbol {\ddot {\rho }}}=-{\mu \over \rho ^{3}}{\boldsymbol {\rho }}}$ for the unperturbed orbit,

(4)

where ${\displaystyle \mu =G(M+m)}$ is the gravitational parameter with ${\displaystyle M}$ and ${\displaystyle m}$ the masses of the central body and the perturbed body, ${\displaystyle \mathbf {a} _{\text{per}}}$ is the perturbing acceleration, and ${\displaystyle r}$ and ${\displaystyle \rho }$ are the magnitudes of ${\displaystyle \mathbf {r} }$ and ${\displaystyle {\boldsymbol {\rho }}}$.

Substituting from equations (3) and (4) into equation (2),

${\displaystyle {\ddot {\delta \mathbf {r} }}=\mathbf {a} _{\text{per}}+\mu \left({{\boldsymbol {\rho }} \over \rho ^{3}}-{\mathbf {r} \over r^{3}}\right),}$

(5)

which, in theory, could be integrated twice to find ${\displaystyle \delta \mathbf {r} }$. Since the osculating orbit is easily calculated by two-body methods, ${\displaystyle {\boldsymbol {\rho }}}$ and ${\displaystyle \delta \mathbf {r} }$ are accounted for and ${\displaystyle \mathbf {r} }$ can be solved. In practice, the quantity in the brackets, ${\displaystyle {{\boldsymbol {\rho }} \over \rho ^{3}}-{\mathbf {r} \over r^{3}}}$, is the difference of two nearly equal vectors, and further manipulation is necessary to avoid the need for extra significant digits.^[13][14] Encke's method was more widely used before the advent of modern computers, when much orbit computation was performed on mechanical calculating machines.