HOODESolver.jl: A Julia package for highly oscillatory problems

Highly oscillatory ordinary differential equations (ODEs) have a long history since they are ubiquitous to describe dynamical multi-scale physical phenomena in physics or chemistry. They can be obtained by an appropriate spatial discretization of partial differential equations or can directly describe the behavior of dynamical quantities. In addition to the standard difficulties relating to their numerical resolution, highly oscillatory ODEs involve a stiffness (characterized by a parameter ε ∈ ]0, 1]) which gives rise to high oscillations in the solution. Hence, to capture these small scales (or high oscillations), conventional methods have to consider a time step smaller than ε leading to unacceptable computational cost.


Summary
Highly oscillatory ordinary differential equations (ODEs) have a long history since they are ubiquitous to describe dynamical multi-scale physical phenomena in physics or chemistry. They can be obtained by an appropriate spatial discretization of partial differential equations or can directly describe the behavior of dynamical quantities. In addition to the standard difficulties relating to their numerical resolution, highly oscillatory ODEs involve a stiffness (characterized by a parameter ε ∈ ]0, 1]) which gives rise to high oscillations in the solution. Hence, to capture these small scales (or high oscillations), conventional methods have to consider a time step smaller than ε leading to unacceptable computational cost.
We present here HOODESolver.jl 1 , a general-purpose library written in Julia dedicated to the efficient resolution of highly oscillatory ODEs. In the documentation 2 details are given to explain how to simulate highly oscillatory ODEs using a Uniformly Accurate (UA) method ie a method which is able to capture the solution while keeping the time step (and then the computational cost) independent of the degree of stiffness ε.

Statement of need
The Julia package DifferentialEquations.jl (Rackauckas & Nie, 2017) efficiently solves many ODE problems using recent and advanced numerical techniques. However the available algorithms do not easily solve the stiff problems discussed above, because they do not take into account the highly oscillatory character of the solutions. Indeed, the solution presents oscillations whose period is proportional to ε. If ε is small, conventional methods struggle to solve such multi-scale phenomena since they require to use tiny time steps to capture high oscillations and become computationally very costly. On the one side, specific methods inspired by the averaging theory have been designed to deal with the regime ε ≪ 1. On the other side, when ε ∼ 1 the problem ceases to be stiff and a classical integrator gives accurate result in reasonable time. The true difficulty emerges for intermediate values of ε, for which averaging techniques are not accurate enough and, due to computational cost, standard methods are inefficient. Thus, a new paradigm has been recently introduced, the so-called uniformly accurate (UA) methods which are indeed able to solve the original highly oscillatory problem with a precision and a computational cost that are independent of the value ε. In particular, these methods allows skipping several oscillations per time step, reducing the number of iterations (and then the cost of the simulation) drastically.
HOODESolver.jl intends to gather and unify recent research around highly oscillatory problems (Bao & Zhao, 2019;P. Chartier et al., 2015P. Chartier et al., , 2020Crouseilles et al., 2017); its development has been motivated by these research needs and it has already been used in some papers (Philippe Chartier et al., 2021). HOODESolver.jl provides software implementations of several theoretical ideas contained in the recent literature around the so-called two-scale method. In particular, a very recent extension proposed in (Philippe Chartier et al., 2021) enables to reach high order accuracy. The implementation focuses on a multistep method (namely Adams-Bashforth method) coupled with a spectral method for the discretization of the additional variable representing the fast scale. Hence, HOODESolver.jl provides an efficient way for researchers to solve a highly oscillatory ODE system, and as such it can be used by the scientific community: • researchers interested in solving highly oscillatory problems arising in their research field (electromagnetic waves, quantum mechanics, plasma physics, molecular dynamics, . . .), • it can guide some future potential numerical or theoretical developments, • it will serve as a reference to benchmark a new method designed by researchers.

Features
HOODESolver.jl is designed to solve the following highly oscillatory ordinary differential systemu where The numerical solution of Equation 1 is computed by simply entering the different components of the equation (A, f , ε, t start , t end , u in ) following the required format. The user simply chooses an order order of the Adams-Bashforth time integrator and the time step h = (t start − t end )/nb_t. The result is given as a function object which can be evaluated in an arbitrary time t, not just at the discrete times. In addition to the methodology introduced in HOODESolver.jl, the package includes: 1. Arbitrary precision arithmetic via BigFloats, 2. New technique to compute the first iterations required for the initialization of the Adams-Bashforth method, this requires that f has to be order times differentiable on [t start − order h, t end ], 3. Extension of the two-scale method to non-homogeneous problems.
The package has been thought to be in close connection to DifferentialEquation.jl. We offer a common interface with it by extending the SplitODE problem type 3 . Users can use our package more easily and it facilitates the cross comparisons with other methods.
The function LinearHOODEOperator has been written following the interface AbstractDif fEqLinearOperator 4 in order to solve a SplitODEProblem using the HOODEAB algorithm. It defines the stiff operator 1 ε A with both ε and A from the studied Equation 1.

Example
The following is an example with the system of Hénon-Heiles 5 : using

Related research and software
The development of the HOODESolver.jl package was initially motivated by the need of efficient multiscale solvers for the charged particles dynamics in an external strong magnetic field. Indeed, due to the Lorentz force, charged particles undergo rapid circular motion around the magnetic field lines. This constitutes the basis of the magnetic confinement of a plasma in a chamber. Obviously, computing a highly oscillatory dynamics is a long-standing problem occurring in many relevant applications (Engquist et al., 2009;Hairer et al., 2006). However, we are not aware of other software packages with similar purpose, excepting the very recent (py)oscode package (Agocs, 2020) which combines WKB techniques and standard integration methods to ensure a user-specified tolerance.