Logo
ScienceStack
Table of Contents
Fetching ...
Sign In

An accelerated frequency-independent solver for oscillatory differential equations

Tara Stojimirovic, James Bremer

TL;DR

This work introduces a frequency-independent solver for oscillatory second order linear ODEs of the form $y''(t) + \omega^2 q(t,\omega) y(t)=0$. It builds a slowly varying phase function via a Riccati-based representation, discrete it with Chebyshev spectral collocation, and solve with Newton-Kantorovich iterations, complemented by Appell's equation to propagate the phase across low-frequency regions. Theoretical results establish the existence of a nonoscillatory phase and quadratic convergence of the Newton-Kantorovich scheme, while the four-stage algorithm adaptively constructs and integrates the phase across the domain. Numerical experiments show orders-of-magnitude speedups over state-of-the-art methods (including Modified Magnus, ARDC, and BremerPhase) while maintaining high accuracy, demonstrating practical impact for evaluating special functions and solving Sturm–Liouville-type problems. Limitations in the low-frequency regime near turning points are discussed, with plans to address them via alternative phase representations and parallelization opportunities.

Abstract

Oscillatory second order linear ordinary differential equations arise in many scientific calculations. Because the running times of standard solvers increase linearly with frequency when they are applied to such problems, a variety of specialized methods, most of them quite complicated, have been proposed. Here, we point out that one of the simplest approaches not only works, but yields a scheme for solving oscillatory second order linear ordinary differential equations which is significantly faster than current state-of-the-art techniques. Our method, which operates by constructing a slowly varying phase function representing a basis of solutions of the differential equation, runs in time independent of the frequency and can be applied to second order equations whose solutions are oscillatory in some regions and slowly varying in others. In the high-frequency regime, our algorithm discretizes the nonlinear Riccati equation satisfied by the derivative of the phase function via a Chebyshev spectral collocation method and applies the Newton-Kantorovich method to the resulting system of nonlinear algebraic equations. We prove that the iterates converge quadratically to a nonoscillatory solution of the Riccati equation. The quadratic convergence of the Newton-Kantorovich method and the simple form of the linearized equations ensure that this procedure is extremely efficient. Our algorithm then extends the slowly varying phase function calculated in the high-frequency regime throughout the solution domain by solving a certain third order linear ordinary differential equation related to the Riccati equation. We describe the results of numerical experiments showing that our algorithm is orders of magnitude faster than existing schemes, including the modified Magnus method [18], the current state-of-the-art approach [7] and the recently introduced ARDC method [1].

An accelerated frequency-independent solver for oscillatory differential equations

TL;DR

This work introduces a frequency-independent solver for oscillatory second order linear ODEs of the form . It builds a slowly varying phase function via a Riccati-based representation, discrete it with Chebyshev spectral collocation, and solve with Newton-Kantorovich iterations, complemented by Appell's equation to propagate the phase across low-frequency regions. Theoretical results establish the existence of a nonoscillatory phase and quadratic convergence of the Newton-Kantorovich scheme, while the four-stage algorithm adaptively constructs and integrates the phase across the domain. Numerical experiments show orders-of-magnitude speedups over state-of-the-art methods (including Modified Magnus, ARDC, and BremerPhase) while maintaining high accuracy, demonstrating practical impact for evaluating special functions and solving Sturm–Liouville-type problems. Limitations in the low-frequency regime near turning points are discussed, with plans to address them via alternative phase representations and parallelization opportunities.

Abstract

Oscillatory second order linear ordinary differential equations arise in many scientific calculations. Because the running times of standard solvers increase linearly with frequency when they are applied to such problems, a variety of specialized methods, most of them quite complicated, have been proposed. Here, we point out that one of the simplest approaches not only works, but yields a scheme for solving oscillatory second order linear ordinary differential equations which is significantly faster than current state-of-the-art techniques. Our method, which operates by constructing a slowly varying phase function representing a basis of solutions of the differential equation, runs in time independent of the frequency and can be applied to second order equations whose solutions are oscillatory in some regions and slowly varying in others. In the high-frequency regime, our algorithm discretizes the nonlinear Riccati equation satisfied by the derivative of the phase function via a Chebyshev spectral collocation method and applies the Newton-Kantorovich method to the resulting system of nonlinear algebraic equations. We prove that the iterates converge quadratically to a nonoscillatory solution of the Riccati equation. The quadratic convergence of the Newton-Kantorovich method and the simple form of the linearized equations ensure that this procedure is extremely efficient. Our algorithm then extends the slowly varying phase function calculated in the high-frequency regime throughout the solution domain by solving a certain third order linear ordinary differential equation related to the Riccati equation. We describe the results of numerical experiments showing that our algorithm is orders of magnitude faster than existing schemes, including the modified Magnus method [18], the current state-of-the-art approach [7] and the recently introduced ARDC method [1].
Paper Structure (24 sections, 4 theorems, 111 equations, 5 figures)

This paper contains 24 sections, 4 theorems, 111 equations, 5 figures.

Table of Contents

  1. Introduction
  2. Mathematical and Numerical Preliminaries
  3. Notation and conventions
  4. The Newton-Kantorovich theorem
  5. The Banach fixed point theorem
  6. The Riccati equation and its variants
  7. A bound on the magnitude of Chebyshev coefficients
  8. Existence of Nonoscillatory Phase Functions
  9. Elementary proof of the existence of a nonoscillatory phase function
  10. A stronger results regarding nonoscillatory solutions of the Riccati equation
  11. Analysis of the Discretized Riccati Equation
  12. Numerical Algorithm
  13. Stage one: adaptive discretization of the coefficient $q$
  14. Stage two: solution of the Riccati equation and left-to-right sweep
  15. Stage three: right-to-left sweep
  16. ...and 9 more sections

Key Result

Theorem 1

Suppose that $\Omega$ is an open subset of the Banach space $X$, that $Y$ is a Banach space and that $F: \Omega \subset X \to Y$ is continuously differentiable. Suppose also that there exist a point $x_0 \in \Omega$ and constants $\lambda$ and $\eta$ such that Then, $F'(x)$ has a bounded inverse $F'(x)^{-1} \in L(Y,X)$ for each $x \in B_\eta(x_0)$, the sequence $\{x_n\}$ defined by is contained

Figures (5)

  • Figure 1: The results of the experiment of Section \ref{['section:experiments1']} in which the modified Magnus method of ModMagnus and the algorithm of this paper were used to evaluate the Legendre function $L_n(t)$ defined in (\ref{['experiments:lnu']}) on the interval $[0,0.9]$. The plot on the left gives the time in seconds required by each method as a function of the degree $n$ of the Legendre function, while the plot on the right gives the maximum relative error observed in the course of evaluating the Legendre function using each method. The plot on the right also shows the accuracy predicted by the condition number of evaluation of the Legendre function. We give the results for the modified Magnus method when the step size $h$ was taken to be a constant multiple of $n^{-1/2}$, and when it was taken to be a constant multiple of $n^{-3/4}$.
  • Figure 2: The results of the experiment of Subsection \ref{['section:experiments2']} in which the method of this paper, the smooth deformation scheme of BremerPhase and the ARDC method of ARDC were used to evaluate the Legendre function $L_n$ defined in (\ref{['experiments:lnu']}) on the interval $[0,0.999]$. The plot on the left gives the time required by each algorithm as a function of the degree $n$ of the Legendre function, while the plot on the right shows the maximum relative error observed when $L_n$ was evaluated at $100$ equispaced points on the interval $[0,.999]$ using each approach. The plot on the right also gives the maximum relative error predicted by the condition number of evaluation of this function.
  • Figure 3: The results of the experiment of Subsection \ref{['section:experiments3']} in which the accuracy of the phase functions computed by our algorithm was measured. The plot on the left gives the time taken by our algorithm as a function of the frequency parameter $n$. The plot on the right gives the maximum relative error observed while evaluating the solution of Kummer's equation calculated by our algorithm at 1,000 points in the interval $\left[0,1-10^{-7}\right]$, including the points $0$ and $1-10^{-7}$.
  • Figure 4: The results of the experiment of Subsection \ref{['section:experiments4']} in which our algorithm was used to evaluate Gegenbauer polynomials of orders $\alpha=-0.499,0.25,1.0$. The plot on the left gives the time required to construct each phase function. The plot on the right gives the maximum absolute accuracy observed while evaluating the Gegenbauer polynomials at 1,000 equispaced points on the interval $[0,0.999]$.
  • Figure 5: The results of the experiment of Subsection \ref{['section:experiments5']} in which our algorithm was applied to boundary value problem for the differential equation (\ref{['experiments:bvp']}) and the accuracy of the resulting solution was tested using a standard adaptive Chebyshev spectral solver. The plot on the left gives the time required to solve the problem using each method as a function of the frequency $\omega$, while the plot on the right gives the maximum absolute error observed while comparing the solution obtained using our method with the reference solution computed via the standard solver.

Theorems & Definitions (4)

  • Theorem 1: Newton-Kantorovich
  • Theorem 2: Banach fixed point theorem
  • Theorem 3
  • Theorem 4