(py)oscode: fast solutions of oscillatory ODEs

The frequency ω(x) and damping γ(x) terms do not need to be explicit functions of x (they can instead be e.g. the result of another numerical solution of an ODE), as they are supplied as sequences ωj , γj evaluated at xi ≤ xj ≤ xf , where (xi, xf ) is the integration range. (py)oscode is written in C++, but comes with a Python wrapper. Its Python interface was designed to be similar to those included in SciPy’s (Virtanen et al., 2020) numerical ODE solution modules. This is demonstrated in the example below whose output is shown in Figure 1.

We present (py)oscode, a general-purpose numerical routine for solving a class of highly oscillatory ordinary differential equations (ODEs) efficiently. The package has been designed to solve equations which describe a single harmonic oscillator with a time-dependent frequency and damping term, i.e. are of the form y ′′ + 2γ(x)y ′ + ω 2 (x)y = 0. (1) The frequency ω(x) and damping γ(x) terms do not need to be explicit functions of x (they can instead be e.g. the result of another numerical solution of an ODE), as they are supplied as sequences ω j , γ j evaluated at (py)oscode is written in C++, but comes with a Python wrapper. Its Python interface was designed to be similar to those included in SciPy's (Virtanen et al., 2020) numerical ODE solution modules. This is demonstrated in the example below whose output is shown in Figure 1.
import numpy as np import scipy.special as sp import pyoscode The increase in step-size of pyoscode's internal steps (orange dots) is due to the algorithm switching from using the RK method to the WKB approximation in the presence of high-frequency oscillations. The orange segment shows dense output, the solution at these points was computed at no additional evaluations of terms in the differential equation.

Statement of need
Even if the terms in Equation 1 change slowly, if the frequency of oscillations in the solution is high enough, standard numerical methods struggle to solve such equations quickly. Traditional methods have to trace every oscillation in the solution, taking many steps in x at an enormous computational cost. The algorithm underlying (py)oscode, published in  and based on W. J. Handley et al. (2016), can detect when the solution is oscillatory and switch to a method based on an analytic approximation (Wentzel-Kramers-Brillouin, WKB) suited for oscillatory functions, otherwise using a Runge-Kutta (RK) method. Using the WKB approximation allows the algorithm to skip over several wavelengths of oscillation in a single step, reducing the number of steps taken drastically. It adaptively updates its step-size to keep the local numerical error within a user-specified tolerance. (py)oscode is capable of producing a solution estimate at an arbitrary value of x, not just at its internal steps, therefore it can be used to generate a "continuous" solution, or dense output .

Related research and software
(py)oscode's development was motivated by the need for a significantly more efficient solver for the evolution of early-universe quantum fluctuations. These perturbations are thought to have been stretched to macroscopic scales by a phase of accelerated expansion of the universe (cosmic inflation), to later become the large-scale structure we see today. To understand the origins of structure it is therefore essential to model the perturbations and understand the physics involved in inflation. (py)oscode has been used to speed up the numerical evolution of quantum fluctuations in the early universe, enabling the exploration of models beyond the standard model of cosmology (W. Handley, 2019). It served as inspiration for other numerical methods aiming to extend the range of oscillatory ODEs to solve (Bamber & Handley, 2020).
The efficient solution of oscillatory ODEs is a long-standing numerical analysis problem with many existing methods to handle certain sub-classes of equations. Examples include successful methods by Petzold (Petzold, 1981), reviewed in Petzold et al. (1997) with many references therein, Iserles et al. (Condon et al., 2009(Condon et al., , 2011Deaño et al., 2017), and Bremer (Bremer, 2018a), with code available from Bremer (2018b).