Logo
ScienceStack
Table of Contents
Fetching ...
Sign In

Fast and stable computation of highly oscillatory and/or exponentially decaying integrals using a Clenshaw-Curtis product-integration rule

Victor Dominguez

TL;DR

The paper tackles the challenge of numerically evaluating $I_z(f)=\int_0^2 f(s) e^{zs}\,{\rm d}s$ for complex $z$ with potentially large negative real parts, where the integrand is oscillatory and/or exponentially decaying. It introduces a Clenshaw-Curtis product-integration rule that interpolates $f$ at Chebyshev nodes and integrates the interpolant exactly, with weights $\omega_\ell(z)=\int_0^2 T_\ell(s-1) e^{zs}\,{\rm d}s$ computed stably via a three-phase, near-linear-cost algorithm that remains robust across all $z$ and avoids special functions. The authors provide convergence and stability analyses, showing error bounds that decay with the node count $L$ and favorable behavior as $|z|$ grows, and they develop a practical implementation strategy that scales efficiently with $L$ while maintaining numerical stability. The approach is applicable to Laplace-transform-based solutions of fractional-order evolution equations, among other problems, enabling fast, accurate quadrature across a wide range of exponential parameters and offering a robust, function-free alternative to traditional oscillatory-quadrature techniques.

Abstract

We propose, analyze, and implement a quadrature method for evaluating integrals of the form $\int_0^2 f(s)\exp(zs)\, {\rm d}s$, where $z$ is a complex number with a possibly large negative real part. The integrand may exhibit exponential decay, highly oscillatory behavior, or both simultaneously, making standard quadrature rules computationally expensive. Our approach is based on a Clenshaw-Curtis product-integration rule: the smooth part of the integrand is interpolated using a polynomial at Chebyshev nodes, and the resulting integral is computed exactly. We analyze the convergence of the method with respect to both the number of nodes and the parameter $z$. Additionally, we provide a stable and efficient implementation whose computational cost is essentially independent of $z$ and scales linearly with $L$. Notably, our approach avoids the use of special functions, enhancing its numerical robustness.

Fast and stable computation of highly oscillatory and/or exponentially decaying integrals using a Clenshaw-Curtis product-integration rule

TL;DR

The paper tackles the challenge of numerically evaluating for complex with potentially large negative real parts, where the integrand is oscillatory and/or exponentially decaying. It introduces a Clenshaw-Curtis product-integration rule that interpolates at Chebyshev nodes and integrates the interpolant exactly, with weights computed stably via a three-phase, near-linear-cost algorithm that remains robust across all and avoids special functions. The authors provide convergence and stability analyses, showing error bounds that decay with the node count and favorable behavior as grows, and they develop a practical implementation strategy that scales efficiently with while maintaining numerical stability. The approach is applicable to Laplace-transform-based solutions of fractional-order evolution equations, among other problems, enabling fast, accurate quadrature across a wide range of exponential parameters and offering a robust, function-free alternative to traditional oscillatory-quadrature techniques.

Abstract

We propose, analyze, and implement a quadrature method for evaluating integrals of the form , where is a complex number with a possibly large negative real part. The integrand may exhibit exponential decay, highly oscillatory behavior, or both simultaneously, making standard quadrature rules computationally expensive. Our approach is based on a Clenshaw-Curtis product-integration rule: the smooth part of the integrand is interpolated using a polynomial at Chebyshev nodes, and the resulting integral is computed exactly. We analyze the convergence of the method with respect to both the number of nodes and the parameter . Additionally, we provide a stable and efficient implementation whose computational cost is essentially independent of and scales linearly with . Notably, our approach avoids the use of special functions, enhancing its numerical robustness.

Paper Structure

This paper contains 17 sections, 9 theorems, 162 equations, 5 figures, 4 tables, 4 algorithms.

Table of Contents

  1. Keywords:
  2. AMS subject classifications:
  3. Introduction
  4. Oscillatory and/or exponentially decaying integrals and product Clenshaw-Curtis rules
  5. Error for the product Clenshaw-Curtis rule
  6. Implementation
  7. Phase 1: Computation of $\omega_n(z)$ and $\rho_n(z)$ for $n=0,\ldots,L_0$, with $L_0 \le n_0$
  8. Phase 2: Computation of $\omega_n(z)$ and $\rho_n(z)$ for $n=n_0+1,\ldots..., L$
  9. Computation of $\rho_{m }(z)$ for some $m\ge L$
  10. Summary: Fast and stable computation of the weights for all $L$ and $z$
  11. Numerical experiments
  12. Implementation of the rule
  13. Stability in the computation of the weights
  14. Error of the quadrature rule
  15. Applications to the solutions of fractional order evolution equations
  16. ...and 2 more sections

Key Result

Lemma 3.1

For $n \ge m+1 \ge 1$, and for $f$ sufficiently smooth, the following identity holds:

Figures (5)

  • Figure 1: Modulus of the weights $\bm{\rho}^{\rm unstab}_n(z),$ (blue) and $\bm{\rho}^{\rm unstab,qp}_n(z),$ (orange) for $n=0,\ldots,250$ and different values of $z$. We observe that instabilities affect the computation, the worst scenario occurs when $z =-40\pi\approx 126$. Quadruple precision arithmetic only mitigates the problem
  • Figure 3: Error between $|\bm{\rho}_n(z)- \bm{\rho}^{\rm qp}_n(z)|$ for several values of $z$. The calculations remains stable in all the range
  • Figure 4: Absolute error in the computation of $I(n,z)$ cf. \ref{['eq:Inz']} for $|z|=250$. Notice that the error in exact arithmetic should be zero
  • Figure 5: Relative error in the computation of $I(n,z)$ cf. \ref{['eq:Inz']} for $|z|=250$. In this case the exact integral decay exponentially to zero which affects the computations due to cancellation error in the computation of the interpolating polynomial in the Chebyshev basis.
  • Figure 6: Domain of the problem

Theorems & Definitions (13)

  • Remark 2.1
  • Lemma 3.1
  • Proposition 3.2
  • Proposition 3.3
  • Lemma 3.4
  • Theorem 3.5
  • Proposition 4.1
  • Remark 4.2
  • Theorem 4.3
  • Remark 4.4
  • ...and 3 more