On the adaptive Levin method

TL;DR

The paper addresses the numerical evaluation of oscillatory integrals $I=\int_a^b e^{i g(x)} f(x)\,dx$ and the issue of low-frequency breakdown in Levin-type methods. It proves that a Chebyshev spectral discretization of the Levin equation, solved with a truncated SVD, remains accurate regardless of the magnitude of $g'$ and even in the presence of stationary points, and it extends this framework to adaptive schemes. The authors show that the cost of representing the Levin solution via a polynomial expansion can be made independent of $|g'|$ (and decreases when $|g'|$ is small), enabling robust adaptive subdivision and integration. They further demonstrate that coupling the adaptive Levin method with phase-function representations for ODEs enables efficient evaluation of a wide class of oscillatory integrals involving special functions. Comprehensive numerical experiments across elementary functions, Bessel, Legendre, Hermite functions, and Green’s functions illustrate the method’s accuracy, efficiency, and versatility, highlighting its practical impact as a general-purpose tool for oscillatory integrals with stationary points or many oscillations.

Abstract

The Levin method is a well-known technique for evaluating oscillatory integrals, which operates by solving a certain ordinary differential equation in order to construct an antiderivative of the integrand. It was long believed that this approach suffers from "low-frequency breakdown," meaning that the accuracy of the calculated value of the integral deteriorates when the integrand is only slowly oscillating. Recently presented experimental evidence, however, suggests that if a Chebyshev spectral method is used to discretize the differential equation and the resulting linear system is solved via a truncated singular value decomposition, then no low-frequency breakdown occurs. Here, we provide a proof that this is the case, and our proof applies not only when the integrand is slowly oscillating, but even in the case of stationary points. Our result puts adaptive schemes based on the Levin method on a firm theoretical foundation and accounts for their behavior in the presence of stationary points. We go on to point out that by combining an adaptive Levin scheme with phase function methods for ordinary differential equations, a large class of oscillatory integrals involving special functions, including products of such functions and the compositions of such functions with slowly-varying functions, can be easily evaluated without the need for symbolic computations. Finally, we present the results of numerical experiments which illustrate the consequences of our analysis and demonstrate the properties of the adaptive Levin method.

Figures (11)

• Figure 1: The results of the first two experiments of Section \ref{['section:experiments:elem']}. The first row of plots pertains to the integral $I_1$ and the second to $I_2$. Each plot on the left gives the error in the value of the integral computed via the adaptive Levin method as a function of $\lambda$, while each plot on the right gives the running time in milliseconds as a function of $\lambda$.
• Figure 2: The results of the last two experiments of Section \ref{['section:experiments:elem']}. The first row of plots pertains to the integral $I_3$ and the second to $I_4$. Each plot on the left gives the error in the value of the integral computed via the adaptive Levin method as a function of $\lambda$, while each plot on the right gives the running time in milliseconds as a function of $\lambda$.
• Figure 3: The results of the experiments of Section \ref{['section:experiments:saddle']}. The plots in the first row give the number of subintervals in the adaptive discretization of $[-1,1]$ formed in the course of evaluating $I_9(\lambda,m)$ as a function of $\lambda$ for $m=2,3,4,5$; those in the second give the absolute error in the calculated error in the value of $I_9(\lambda,m)$ as a function of $\lambda$ for $m=2,3,4,5$; and the plots in the third row show the number of subintervals in the adaptive discretization of $[-1,1]$ as a function of $m$ for $\lambda=10^2, 10^3, 10^4, 10^5$. The plots in the column on the left concern experiments executed with the precision parameter $\epsilon$ taken to be $10^{-7}$, while those on the right correspond to $\epsilon = 10^{-12}$.
• Figure 4: The results of the first two experiments of Section \ref{['section:experiments:bes']}. In the first row, the plot on the left gives the absolute error in the calculation of the integral $I_{10}(\nu)$ as a function $\nu$ and the plot on the right shows the time take by the adaptive Levin method and the time taken to constuct the phase function as functions of $\nu$. The plot at bottom left gives the absolute error in the obtained value of $I_{11}(\lambda)$ as a function of $\lambda$, and the bottom-right plot gives the time required to compute $I_{11}(\lambda)$, including the time required to construct any necessary phase functions, as a function of $\lambda$.
• Figure 5: The results of the last two experiments of Section \ref{['section:experiments:bes']}. The first row pertains to $I_{12}(\lambda)$ while the second concerns $I_{13}(\lambda)$. In each row, the plot on the left gives the absolute error in the calculation of the integral as a function of $\lambda$ and the plot on the right shows the total time required to compute the integral via the adaptive Levin method and to construct any necessary phase functions as a function of $\lambda$.

Theorems & Definitions (22)

• Theorem 1
• proof
• Definition 1
• Corollary 1
• Lemma 1
• proof
• Lemma 2
• proof
• Lemma 3
• proof