Abstract:
Oscillatory systems are ubiquitous in physics: they arise in celestial and quantum mechanics, electrical circuits, molecular dynamics, and beyond. Yet even in the simplest case, when the frequency of oscillations changes slowly but is large, the vast majority of numerical methods struggle to solve such equations. Methods based on approximating the solution with polynomials are forced to take O(ω) timesteps, where ω is the characteristic frequency of oscillations. This scaling can generate unacceptable computational costs when the ODE in question needs to be solved billions of times, e.g. as the forward modeling step of Bayesian parameter estimation.
In this talk I will introduce an efficient method for solving 2nd order, linear, homogeneous ODEs whose solutions may vary between highly oscillatory and slowly changing over the solution interval.
The solver employs two methods: in regions where the solution varies slowly, it uses Chebyshev-grid based collocation with an adaptive stepsize, but in the highly oscillatory phase it automatically switches to constructing a local phase function. I propose a defect-correction iteration that gives an asymptotic series for such a phase function; this is numerically approximated on a Chebyshev grid. In the talk I will present how the method fits in the landscape of oscillatory solvers, details of the algorithm, results from numerical experiments, and a brief error analysis.