# A Simple Method for Calibrating a Numerical Integrator

One of the most common tasks in astrodynamics is the numerical solution of the orbital equations of motion. For motion which is not Keplerian or two-body, this is usually accomplished with some type of numerical method which propagates a spacecraft’s motion from one time epoch to another (earlier or later) epoch.

For elliptical orbits, it is often necessary to use a variable step size integrator. These techniques change the integration step size as a function of a spacecraft’s location in its orbit. Classical methods include Runge-Kutta algorithms that solve the first-order form of the differential equations of orbital motion. My favorite algorithm is a Runge-Kutta-Fehlberg (RKF) method because the integration coefficients have been “tailored” to solve the types of differential equations found in astrodynamics and celestial mechanics.

Many of these numerical methods allow the user to control how well the differential equations are solved via a truncation error tolerance or perhaps absolute and relative error controls. Setting these algorithm controls to correct values determines the performance of the numerical method.

If we formulate the differential equations of un-perturbed orbital motion in first-order form, it is possible to study the effect of algorithm controls on the orbit propagation using the following simple method.

The initial conditions for the x, y and z components of the position and velocity vectors are

r_x(0) = 1 – ecc; r_y(0) = 0; r_z(0) = 0;

v_x(0) = 0; v_y(0) = sqrt[ (1 + ecc) / (1 – ecc) ]; v_z(0) = 0

where ecc is the orbital eccentricity.

The system of six first-order differential equations (y_dot) is given by

y_dot(1) = v_x; y_dot(2) = v_y; y_dot(3) = v_z

y_dot(4) = –r_x / r^3; y_dot(5) = –r_y / r^3; y_dot(6) = –r_z / r^3

For this formulation, the Keplerian or un-perturbed orbital period is exactly equal to 2 pi. In the last three equations r is the scalar magnitude of the position vector; r = sqrt(r_x^2 + r_y^2 + r_z^2).

To test an integration method, select an eccentricity value, say 0.75, set typical values for the algorithm controls, and run the propagator for an integer number of 2 pi “orbits”, printing the position and velocity vector components at “time” intervals of 2 pi. An examination of the state vector (position and velocity vector components) at the end of each complete orbit will indicate how well the integrator matches the initial conditions defined above as the propagation time increases. The goal is to select the “best” algorithm control for acceptable state vector differences for a user-defined propagation duration. A highly elliptical orbit exercises the step size features of a numerical method.

This method can also be used with dimensional values (perhaps kilometers for position and kilometers per second for velocity) for the initial state vector and Keplerian orbital period (seconds). This form also requires the value of the Earth’s gravitational constant (mu) in the three differential equations expressed as

y_dot(4) = –mu * r_x/r^3; y_dot(5) = –mu * r_y/r^3;

y_dot(6) = –mu * r_z /r^3

We have used this technique to calibrate an RKF method for (1) propagating orbital debris in the geosynchronous orbit region for one hundred years, (2) studying the lifetime of Earth satellites, (3) propagating lunar and interplanetary trajectories, (4) predicting impact locations for re-entering orbital debris, etc.