[Paper Review] On the numerical solution of second order differential equations in the high-frequency regime

James Bremer|arXiv (Cornell University)|Sep 21, 2014

Numerical methods for differential equations21 references4 citations

TL;DR

This paper presents a fast, frequency-independent algorithm for solving high-frequency second-order linear ordinary differential equations by leveraging nonoscillatory phase functions. It reformulates the problem using a logarithmic form of Kummer's equation to compute phase functions that avoid costly oscillation resolution, enabling solutions with runtime independent of the high frequency parameter λ, as demonstrated on Bessel, Legendre, and prolate spheroidal wave functions with high accuracy and efficiency.

ABSTRACT

We describe an algorithm for the numerical solution of second order linear differential equations in the highly-oscillatory regime. It is founded on the recent observation that the solutions of equations of this type can be accurately represented using nonoscillatory phase functions. Unlike standard solvers for ordinary differential equations, the running time of our algorithm is independent of the frequency of oscillation of the solutions. We illustrate the performance of the method with several numerical experiments.

Motivation & Objective

• To address the inefficiency of standard ODE solvers in the high-frequency regime, where solution oscillations require excessive computational cost proportional to λ.
• To develop a general-purpose algorithm applicable to a wide class of second-order ODEs, including those defining special functions.
• To achieve solution accuracy independent of λ by using nonoscillatory phase functions derived from a logarithmic form of Kummer’s equation.
• To enable fast evaluation of special functions like Bessel, Legendre, and prolate spheroidal wave functions with high precision and minimal computational overhead.

Proposed method

• The method constructs a nonoscillatory phase function α(t) such that solutions of the ODE y'' + λ²q(t)y = 0 are represented via u(t) = cos(α(t))/|α'(t)|¹ᐟ² and v(t) = sin(α(t))/|α'(t)|¹ᐟ².
• It transforms Kummer’s third-order nonlinear ODE into a more tractable logarithmic form via r'(t) = λ exp(r(t)/2), leading to the equation r'' - ¼(r')² + 4λ²(exp(r) - q) = 0.
• The algorithm solves this logarithmic form numerically using a graded mesh near singularities and high-order Chebyshev spectral collocation on subintervals.
• Solutions are represented using high-order Chebyshev expansions on each subinterval, enabling efficient evaluation at arbitrary points.
• The method avoids resolving rapid oscillations directly by working with the phase function, ensuring computational cost is independent of λ.
• It is validated on equations for Bessel functions, Legendre functions, and prolate spheroidal wave functions, with comparisons to existing solvers like Osipov-Rokhlin.

Experimental results

Research questions

• RQ1 Can a general-purpose numerical method be developed for high-frequency second-order ODEs that is independent of the oscillation frequency λ?
• RQ2 Do nonoscillatory phase functions exist for a broad class of ODEs, including those defining special functions, and can they be computed efficiently?
• RQ3 Can the logarithmic form of Kummer’s equation be solved numerically with high accuracy to enable fast evaluation of solutions?
• RQ4 How does the performance and accuracy of this method compare to existing solvers like Osipov-Rokhlin for prolate spheroidal wave functions?
• RQ5 Can the method achieve machine precision accuracy for special functions such as Legendre functions of high order?

Key findings

• The algorithm achieves solution accuracy on the order of (μλ)⁻¹exp(−μλ), where μ depends on q but not λ, ensuring high precision even at large λ.
• The running time of the algorithm is independent of λ, in contrast to standard solvers whose cost scales linearly with λ.
• For Legendre functions of order ν = π×10⁵, the method computed solutions with maximum absolute error below 10⁻¹⁴ using 16th-order Chebyshev expansions on a graded mesh.
• The method outperformed the Osipov-Rokhlin algorithm in evaluation speed for prolate spheroidal wave functions while maintaining comparable accuracy.
• The nonoscillatory phase function for Bessel’s equation was computed accurately without elementary closed-form expression, demonstrating the method’s generality.
• Numerical experiments confirmed that the method maintains high accuracy and efficiency across diverse ODEs, including those with singular coefficients at endpoints.