2.2. Integrators

Warning

The pages in this guide describe the software implementation of numerical integrators within Tudat, thus the reader is expected to be already familiar with the mathematics behind such tools.

Many applications require solving sets of first-order differential equations, where in the field of astrodynamics this is regularly done when using propagators, as described in Simulator Set-Up. The purpose of a numerical integrator is to find a solution to a system of differential equations that satisfies a set of initial conditions, a problem commonly referred to in literature as an Initial-Value Problem (IVP). Tudat includes a framework for using numerical integrators, where several types are offered depending on how the step-size is discretized:

• Single-step methods: These integrators use derivative information from a single step. Among these methods, Tudat includes the Euler method, the Runge-Kutta 4 method and several Runge-Kutta variable step-size methods.
• Multi-step methods: These integrators use derivative information from multiple steps. Note that at the moment multi-step methods have not been implemented in Tudat.

2.2.1. Setting up a numerical integrator

2.2.1.1. Creating an integrator object

All numerical integrators presented in this page follow the same framework, where the first step is to create the integrator object, which has the NumericalIntegratorXdPointer type:

NumericalIntegratorXdPointer integrator
                = boost::make_shared< EulerIntegratorXd >(
                    stateDerivativeFunction,
                    intervalStart,
                    initialState );

Although in this case NumericalIntegratorXdPointer is a typedef for a boost shared_ptr, the constructor can be created explicitly. To initialize an Euler integrator, the following arguments need to be specified:

• EulerIntegratorXd is the input argument for boost’s make_shared, which defines the derived integrator class that we want to use.
• stateDerivativeFunction passes the state derivative function that defines the differential equation to be integrated. See below for more details.
• intervalStart provides the start of the integration interval and must have the type of the independent variable.
• initialState provides the initial state and must have the type of the dependent variable.

It is important to emphasize that the provided state derivative function must have a particular structure for the numerical integration to work. A sample state derivative function is defined here:

static inline Eigen::VectorXd stateDerivativeFunction( const double time,
                                                       const Eigen::VectorXd& state )
{
    ...
}

Any user-defined state derivative function must satisfy the following requirements:

• The first input argument provides the independent variable, which must be a const double and in this case is named time.
• The second input argument provides the dependent variable, which type is free to be defined. In this case, the type to be integrated is a Eigen::VectorXd which is passed by reference.
• The output argument type must coincide with the dependent variable type.

The contents of the state derivative function can be freely defined according to the user’s needs, but it is important to respect the function’s input-output structure.

2.2.1.2. Using the integrator object

Once the integrator object has been created, the multiple functions available within the NumericalIntegrator class can be accessed. The following functions are available and they are exemplified using the integrator and the stateDerivativeFunction declared above:

• getNextStepSize( )

  Returns the step-size at the next integration step. The return type is defined by the type of the independent variable:

  double nextStepSize = integrator->getNextStepSize( );
  

• getCurrentState( )

  Returns the state at the current integration step. The return type is defined by the type of the dependent variable.

  Eigen::VectorXd currentState = integrator->getCurrentState( );
  

• getCurrentIndependentVariable( )

  Returns the value of the independent variable at the current integration step. The return type is defined by the type of the dependent variable.

  double currentIndependentVariable = integrator->getCurrentIndependentVariable( );
  

• performIntegrationStep( stepSize )

  Perform an integration step with the step-size fed to the first argument and return the integration result.

  Eigen::VectorXd stateAfterIntegrationStep = integrator->performIntegrationStep( nextStepSize );
  

• integrateTo( intervalEnd , initialTimeStep )

  Performs an integration until the specified intervalEnd, given the provided initialTimeStep. This function returns the state at the end of the interval.

  Eigen::VectorXd stateAtIntervalEnd = integrator->integrateTo( intervalEnd , initialTimeStep );
  

Note

The functions described above are virtual functions and thus redefined for each integrator method described in this page. Selection of the integrator method is made at the stage of creating the integrator object, where selection of the functions is taken care of by the implementation framework.

2.2.2. Selecting a numerical integrator

2.2.2.1. Euler integrator

The Euler integrator is the simplest integrator available but is also a first-order method, meaning that the global error is proportional to the step-size. Thus, its use for high-accuracy applications is not encouraged. Creating the Euler integrator object is done as follows:

NumericalIntegratorXdPointer integrator
                = boost::make_shared< EulerIntegratorXd >(
                    stateDerivativeFunction,
                    intervalStart,
                    initialState );

2.2.2.2. Runge-Kutta 4 integrator

The Runge-Kutta 4 (RK4) integrator is a fourth-order fixed step-size integrator, thus performing better than the Euler integrator. The RK4 integrator is created as follows:

NumericalIntegratorXdPointer integrator
                = boost::make_shared< RungeKutta4IntegratorXd >(
                    stateDerivativeFunction,
                    intervalStart,
                    initialState );

2.2.2.3. Runge-Kutta variable step-size integrator

The Runge-Kutta variable step-size integrator involves a number of methods where the step-size is adjusted throughout the integration interval to bound the numerical error. Creating such integrators differs from the Euler integrator and the Runge-Kutta 4 fixed step-size methods:

RungeKuttaVariableStepSizeIntegratorXd integrator(
            rungeKuttaCoefficients,
            stateDerivativeFunction,
            initialTime,
            initialState,
            minimumStepSize,
            maximumStepSize,
            relativeErrorTolerance,
            absoluteErrorTolerance )

where the following arguments are necessary:

• RungeKuttaVariableStepSizeIntegratorXd is the input argument for boost’s make_shared, which defines the derived integrator class that we want to use.

• rungeKuttaCoefficients provides the set of coefficients that define the particular variable step-size method being used. A number of Runge-Kutta coefficient sets are available in Tudat:

  // Runge-Kutta-Fehlberg 4(5)
  RungeKuttaCoefficients rungeKuttaCoefficients =
          RungeKuttaCoefficients::get( RungeKuttaCoefficients::rungeKuttaFehlberg45 );
  
  // Runge-Kutta-Fehlberg 5(6)
  RungeKuttaCoefficients rungeKuttaCoefficients =
          RungeKuttaCoefficients::get( RungeKuttaCoefficients::rungeKuttaFehlberg56 );
  
  // Runge-Kutta-Fehlberg 7(8)
  RungeKuttaCoefficients rungeKuttaCoefficients =
          RungeKuttaCoefficients::get( RungeKuttaCoefficients::rungeKuttaFehlberg78 );
  
  // Runge-Kutta-Fehlberg 8(7) Dormand-Prince
  RungeKuttaCoefficients rungeKuttaCoefficients =
          RungeKuttaCoefficients::get( RungeKuttaCoefficients::rungeKuttaFehlberg87 );
  

• stateDerivativeFunction passes the state derivative function that defines the differential equation to be integrated.

• initialTime provides the initial value of the independent variable.

• initialState provides the initial state and must have the type of the dependent variable.

• minimumStepSize defines the minimum step-size that the variable step-size integrator can take.

• maximumStepSize defines the maximum step-size that the variable step-size integrator can take.

• relativeErrorTolerance defines the relative error tolerance.

• absoluteErrorTolerance defines the absolute error tolerance.

2.2.3. Using a numerical integrator to propagate an orbit

The numerical integrators described in this page are commonly used to propagate the orbit of spacecraft and celestial bodies. The reader is referred to Integrator Settings, which discusses the how the numerical integrator fit in the simulator framework of Tudat.