odeint: Ordinary Differential Equations Integrator#

The integration of the equations of motions is taken care by the odeint module. It is based on boost::numeric::odeint but restricted to the needs of the project.

How does it work?#

The problems are described in the form

\[\begin{split}\begin{align} x'(t) &= f(x, t)\\ x(0) &= x_0. \end{align}\end{split}\]

The numerical integration methods calculate the solution to \(x(t)\) iteratively, in general there is a formula that obtains \(x(t + \Delta t)\) as a function of \(x(t)\), \(f\) and \(\Delta t\). For example, the Euler method approximates

\[x(t + \Delta t) = x(t) + x'(t) \Delta t.\]

In odeint the methods of integration are implemented in stepper classes that are only required to have a method do_step(system, state, current_time, delta_time) that takes the equation system \(f\), the current state \(x(t)\), current time \(t\) and the time delta \(\Delta t\), then applies the method in order to approximate the state \(x(t + \Delta t)\) and outputs to the out parameter state.

The do_step method is called iteratively in order to integrate the equation to the desired time using the following function:

template<typename stepper_type, typename system_type, typename state_type, typename scalar_type, typename observer_type>
size_t integrate(stepper_type &stepper, system_type &sys, state_type &x, scalar_type t0, scalar_type dt, size_t Nsteps, observer_type &obs, size_t obs_skip_steps = 0)#

Integrate a system with fixed size steps and method given by the stepper.

Parameters

  • stepper – Stepper that implements the method. Has a do_step method that takes the arguments do_step(system, state, current_time, delta_time); calculates the state at current_time + delta_time and sets state to it.

  • sys – Equation system. It’s called with arguments (x, dxdy, t), x is the current state and t the time; sets dxdy to the derivative of x with respect to t at current state.

  • x – State that will be updated during integration.

  • t0 – Initial time

  • dt – Time delta

  • Nsteps – Amount of steps of integration

  • obs – Observes the state during the integration

  • obs_skip_steps – Amount of steps to skip between observations (default = 0)

This function accepts templated arguments so that it can use different integration methods, equation systems, and observers. The function loops Nsteps times and evaluates the stepper.do_step() with sys, the current time and state and then updates the state and the time (+= dt). obs is called every obs_skip_steps + 1 to observe (output to file or save) the current state.

Note

In the future it might be useful or necessary to implement other integration function that uses adaptative time steps.

Implemented steppers#

The currently implemented steppers are

template<typename system_type, typename state_type, typename scalar_type>
class EulerStepper#

Euler integrator order 1. Just for testing purposes as it is quite unprecise.

template<typename system_type, typename state_type, typename scalar_type>
class RK46NL#

Low-dissipation and low-dispersion fourth-order Runge–Kutta

template<typename system_type, typename state_type = State, typename scalar_type = double>
class Boris#

Boris algorithm for integrating charged particles subjected to Lorentz force [1, 2]

[1] https://www.particleincell.com/2011/vxb-rotation/

[2] Qin, H., Zhang, S., Xiao, J., Liu, J., Sun, Y., & Tang, W. M. (2013). Why is Boris algorithm so good?. Physics of Plasmas, 20(8), 084503. https://doi.org/10.1063/1.4818428