An Adaptive Spectral Method for Oscillatory Second-Order Linear ODEs with Frequency-Independent Cost

SIAM Journal on Numerical Analysis, Volume 62, Issue 1, Page 295-321, February 2024.
Abstract. We introduce an efficient numerical method for second-order linear ODEs whose solution may vary between highly oscillatory and slowly changing over the solution interval. In oscillatory regions the solution is generated via a nonoscillatory phase function that obeys the nonlinear Riccati equation. We propose a defect correction iteration that gives an asymptotic series for such a phase function; this is numerically approximated on a Chebyshev grid with a small number of nodes. For analytic coefficients we prove that each iteration, up to a certain maximum number, reduces the residual by a factor of order of the local frequency. The algorithm adapts both the stepsize and the choice of method, switching to a conventional spectral collocation method away from oscillatory regions. In numerical experiments we find that our proposal outperforms other state-of-the-art oscillatory solvers, most significantly at low to intermediate frequencies and at low tolerances, where it may use up to [math] times fewer function evaluations. Even in high-frequency regimes, our implementation is on average 10 times faster than other specialized solvers.