Computation of highly oscillatory integrals
using a Fourier extension approximation

Akash Anand Akash Anand, Department of Mathematics and Statistics, Indian Institute of Technology, Kanpur, UP 208016 [email protected] and Damini Dhiman Damini Dhiman, Department of Mathematics and Statistics, Indian Institute of Technology, Kanpur, UP 208016 [email protected]

Abstract.

The numerical evaluation of integrals of the form

\displaystyle\int_{a}^{b}f(x)e^{ikg(x)}\,dx

is an important problem in scientific computing with significant applications in many branches of applied mathematics, science and engineering. The numerical approximation of such integrals using classical quadratures can be prohibitively expensive at high oscillation frequency ( $k\gg 1$ ) as the number of quadrature points needed for achieving a reasonable accuracy must grow proportionally to $k$ . To address this significant computational challenge, starting with Filon in 1930, several specialized quadratures have been developed to compute such oscillatory integrals efficiently. A crucial element in such Filon-type quadrature is the accurate evaluation of certain moments which poses a significant challenge when non-linear phase functions $g$ are involved.

In this paper, we propose an equispaced-grid Filon-type quadrature for computing such highly oscillatory integrals that utilizes a Fourier extension of the slowly varying envelope $f$ . This strategy is primarily aimed at significantly simplifying the moment calculations, even when the phase function has stationary points. Moreover, the proposed approach can also handle certain integrable singularities in the integrand. We analyze the scheme to theoretically establish high-order convergence rates. We also include a wide variety of numerical experiments, including oscillatory integrals with algebraic and logarithmic singularities, to demonstrate the performance of the quadrature.

Key words and phrases:

oscillatory integrals, integrable singularities, Fourier extension, Filon quadrature

1. Introduction

The evaluation of highly oscillatory integrals plays a significant role in many branches of applied mathematics. The accuracy of their numerical approximation using a classical quadrature suffers significantly when the oscillation frequency is much more than the number of quadrature points. Since the initial attempt by Filon [1] to solve this difficulty, numerous approaches have been developed [2, 3, 4, 5, 6]. The Filon quadrature for approximating the integral of the form

(1)

\displaystyle\int_{a}^{b}f(x)e^{ikg(x)}\,dx,

relies on the interpolating polynomial $p$ satisfying $p(x_{j})=f(x_{j})$ for interpolation nodes $x_{0},x_{1},\ldots,x_{n},$ to obtain the approximation

(2)

\displaystyle\int_{a}^{b}f(x)e^{ikg(x)}\,dx\approx\int_{a}^{b}p(x)e^{ikg(x)}\,dx.

Both the Filon method and its generalisation [2, 3, 4] where $f$ is approximated by splines, are shown to be efficient for suitably smooth functions provided the moments $\int_{a}^{b}x^{j}e^{ikg(x)}\,dx$ can be explicitly calculated. While obtaining these moments is easy for linear oscillators, that is, $g(x)=\alpha x+\beta,\ \alpha,\beta\in\mathbb{R}$ , in general, the accurate moment calculations pose a significant challenge. For example, The Filon-Clenshaw-Curtis quadrature [7], that uses $\cos(j\pi/N),j=0,1,\dots,N$ as nodes and Chebyshev approximations for $f$ , obtain the underlying moments recursively using multiple recurrence relations. The Filon-Clenshaw-Curtis and its variations [8, 9, 10, 11] have been developed to efficiently handle integrable algebraic and logarithmic singularities in the integrand.

Several Filon-like approaches use the change of variable $y=g(x)$ to transform the integral to the one with linear oscillations

(3)

\displaystyle\int_{a}^{b}f(x)e^{ikg(x)}\,dx=\int_{g(a)}^{g(b)}\frac{1}{g^{\prime}(g^{-1}(y))}f(g^{-1}(y))e^{iky}\,dy

under the assumption that $g^{\prime}(x)>0$ for all $x\in(a,b)$ or $g^{\prime}(x)<0$ for all $x\in(a,b)$ . Indeed, if $g^{\prime}(x)$ has finitely many zeros in $(a,b)$ , say $x_{j},j=1,2,\ldots,J$ , then, by setting $x_{0}=a$ and $x_{J+1}=b$ , we have

\displaystyle\int_{a}^{b}f(x)e^{ikg(x)}\,dx=\sum_{j=0}^{J}\int_{x_{j}}^{x_{j+1}}f(x)e^{ikg(x)}\,dx

where, for each $j$ , the integral over $[x_{j},x_{j+1}]$ can be transformed into the one with linear oscillations.

In this paper, we propose a Filon-type quadrature for oscillatory integrals with equispaced quadrature nodes. We show that this approach not only handles algebraic and logarithmic singularities in the integrand but also significantly simplifies the solution of the underlying moment problem. In fact, for several special singular weights, the proposed approach yields moment problems that can be solved analytically. Moreover, our approach requires only discrete functional data on a uniform grid. The fact that approximating the integral to high order does not require explicit knowledge of the derivatives is a significant strength of this approach.

The method to yield a more tractable moment problem relies on the interpolating trigonometric polynomial approximation of the smooth non-oscillatory part of the integrand. However, in general, such an approximation is inflicted with Gibb’s oscillations that manifest due to the discontinuous periodic extension of the function being approximated. To overcome this difficulty, we first construct a periodic extension of the function, which is then approximated by the interpolating trigonometric polynomial using FFT. We discuss the periodic extension employed in the paper in section 2.2. In section 2.4, we delve into the moment problem and discuss the solution of several important instances. In section 3, we present a convergence analysis. A variety of numerical experiments are presented in section 4 to demonstrate the performance of the proposed quadrature.

2. Methodology

Motivated by the need to compute the integrals such as those that appear on the right-hand side of the eq. 3, we consider the problem of approximation of the integrals

(4)

\displaystyle I^{[a,b]}_{k}(w,f)=\int_{a}^{b}w(x)f(x)e^{ikx}\,dx,

where $w\in L^{1}(a,b)$ is the weight function and $f$ is a smooth function on $[a,b]$ .

2.1. The basic idea

For clarity, we introduce the main ideas assuming that the $f$ is known everywhere in the interval $[a,b]$ and its derivatives are available at $x=a$ and $x=b$ . At the heart of this procedure is the step where we approximate $f$ by a trigonometric polynomial

\displaystyle\sum_{\ell=-n}^{n-1}c_{\ell}(f)e^{2\pi i\ell(x-a)/\xi}

of period $\xi\geq b-a$ . The primary motivation for doing so is to simplify the underlying moment calculations. Indeed, we see that

(5)

\displaystyle I^{[a,b]}_{k}(w,f)\approx\int_{a}^{b}w(x)\left(\sum_{\ell=-n}^{n-1}c_{\ell}(f)e^{2\pi i\ell(x-a)/\xi}\right)e^{ikx}\,dx=e^{ika}\sum_{\ell=-n}^{n-1}c_{\ell}(f)W^{[a,b]}_{k,\ell}(w)

where

(6)

\displaystyle c_{\ell}(f)=\frac{1}{\xi}\int_{a}^{a+\xi}f(x)e^{-2\pi i\ell(x-a)/\xi}\,dx\quad\text{and}\quad W^{[a,b]}_{k,\ell}(w)=\int_{a}^{b}w(x)e^{i(k\xi+2\pi\ell)(x-a)/\xi}\,dx.

For many important weights $w$ , the moments $W^{[a,b]}_{k,\ell}(w)$ can be obtained analytically as we see later in section 2.4.

The natural choice for the approximating trigonometric polynomial of $f$ is to use its truncated Fourier series. However, if we use the period $\xi=b-a$ , it fails to converge uniformly if the $(b-a)$ -periodic extension of $f$ is discontinuous, that is, if $f(a)\neq f(b)$ , due to Gibb’s phenomenon. Instead of directly approximating $f$ , we utilize the Fourier series of a periodic function with a larger period that coincides with $f$ on the interval $[a,b]$ to eliminate Gibb’s oscillations. In the next sub-section, we discuss one such construction that we use for numerical experiments in this paper.

2.2. A periodic extension:

Among several possible periodic extensions for a given function, for simplicity of implementation and analysis, we utilize the construction based on the two-point Hermite polynomial [12, 13]. Given $f$ on $[a,b]$ , we introduce a periodic function $f^{e}$ of period $\xi=2(b-a)$ given by

(7)

\displaystyle f^{e}(x)=\begin{dcases}f(x),&x\in[a,b],\\ p(x),&x\in[b,2b-a),\end{dcases}

where the polynomial $p$ is chosen as

(8)

\displaystyle p(x)=\sum_{m=0}^{r}f^{(m)}(b)p_{m,r}^{[b,2b-a]}(x)+\sum_{m=0}^{r}f^{(m)}(a)p_{m,r}^{[2b-a,b]}(x)

with

\displaystyle p_{m,r}^{[t_{1},t_{2}]}(t)

\displaystyle=\frac{1}{m!}(t-t_{1})^{m}\left(\frac{t-t_{2}}{t_{1}-t_{2}}\right)^{r+1}\sum_{n=0}^{r-m}\binom{r+n}{r}\left(\frac{t-t_{1}}{t_{2}-t_{1}}\right)^{n}.

2.3. The Discrete Problem:

In many practical applications, only discrete functional data is available. In this regard, we consider the equispaced grid of size $n+1$ on the interval $[a,b]$ with grid points $x_{j}=a+(b-a)j/n$ where the functional data $f_{j},j=0,1,\dots,n$ is assumed to be available. We introduce the grid function $f_{r,q}^{e}$ given as

(9)

\displaystyle(f_{r,q}^{e})_{j}=\begin{dcases}f_{j},&j=0,\dots,n,\\ p_{r,q}(a+(b-a)j/n),&j=n+1,\dots,2n,\\ \end{dcases}

where

(10)

\displaystyle p_{r,q}(x)=\sum_{m=0}^{r}\mathcal{D}^{m,-}_{n,q}(f)(x_{n})p_{m,r}^{[b,2b-a]}(x)+\sum_{m=0}^{r}\mathcal{D}^{m,+}_{n,q}(f)(x_{0})p_{m,r}^{[2b-a,b]}(x),

for $1\leq m\leq r$ , the operators $\mathcal{D}^{m,+}_{n,q}$ and $\mathcal{D}^{m,-}_{n,q}$ respectively denote the forward and backward finite difference derivative operators for approximating $f^{(m)}$ of the form

\displaystyle\mathcal{D}^{m,\pm}_{n,q}(f)(x_{j})=(\pm n)^{m}\left(\sum_{\ell=0}^{m+q-1}(a_{q}^{m})_{\ell}f_{j\pm\ell}\right)

with appropriately chosen constants $(a_{q}^{m})_{\ell}$ so that the order of accuracy is $q$ . The discrete Fourier coefficients are obtained as

(11)

\displaystyle d_{\ell,n}(f_{r,q}^{e})=\frac{1}{2n}\sum_{j=0}^{2n-1}(f_{r,q}^{e})_{j}e^{-\pi i\ell j/n}

for $\ell=-n,\dots,n-1$ . Finally, the integration scheme for approximation of $I_{k}^{[a,b]}(w,f)$ is given by

(12)

\displaystyle I_{k}^{[a,b]}(w,f)\approx I_{k,n}^{[a,b]}(w,f):=e^{ika}\sum_{\ell=-n}^{n-1}d_{\ell,n}(f_{r,q}^{e})W_{k,\ell}^{[a,b]}(w).

2.4. The moment problem:

The numerical integration scheme requires obtaining the moments given by

	$\displaystyle W^{[a,b]}_{k,\ell}(w)$	$\displaystyle=\int_{a}^{b}w(x)e^{\pi i\ell(x-a)/(b-a)}e^{ik(x-a)}\,dx=(b-a)\int_{0}^{1}w(a+(b-a)u)e^{i((b-a)k+\pi\ell)u}\,du$
		$\displaystyle=\frac{b-a}{2}e^{i((b-a)k+\pi\ell)/2}\int_{-1}^{1}w\left(\frac{a+b}{2}+\frac{b-a}{2}v\right)e^{i((b-a)k+\pi\ell)v/2}\,dv$

The following are solutions to some instances of the moment problem corresponding to important classes of weights of interest.

Example 2.1.

For the simplest case where $w(x)=1$ , we have

\displaystyle W^{[a,b]}_{k,\ell}

\displaystyle=\begin{cases}b-a,&\ell=-(b-a)k/\pi,\\ \dfrac{2(b-a)}{(b-a)k+\pi\ell}e^{i((b-a)k+\pi\ell)/2}\sin(((b-a)k+\pi\ell)/2),&\ell\neq-(b-a)k/\pi.\end{cases}

Example 2.2.

For one sided algebraic singularity where $w(x)=(x-a)^{\beta},$ $\beta\in(-1,0)$ , the moments are given by

\displaystyle W^{[a,b]}_{k,\ell}

\displaystyle=\begin{cases}\dfrac{(b-a)^{1+\beta}}{1+\beta},&\ell=-(b-a)k/\pi,\\ \left(\dfrac{i(b-a)}{(b-a)k+\pi\ell}\right)^{1+\beta}(\Gamma(1+\beta)-\Gamma(1+\beta,-i((b-a)k+\pi\ell))),&\ell\neq-(b-a)k/\pi.\end{cases}

Example 2.3.

For $w(x)=(x-a)^{\beta}(b-x)^{\beta}$ with $\beta\in(-1,0)$ , the weak algebraic singularity at both end-points of the same order $\beta$ , we have

\displaystyle W^{[a,b]}_{k,\ell}=\sqrt{\pi}\ \Gamma(1+\beta)(b-a)^{1+2\beta}e^{i((b-a)k+\pi\ell)/2}((b-a)k+\pi\ell)^{-1/2-\beta}J_{1/2+\beta}(((b-a)k+\pi\ell)/2)

provided $\ell\neq-(b-a)k/\pi$ , whereas

\displaystyle W^{[a,b]}_{k,\ell}=\left(\frac{b-a}{2}\right)^{1+2\beta}\frac{\sqrt{\pi}\hskip 2.84544pt\Gamma(1+\beta)}{\Gamma\left(\frac{3}{2}+\beta\right)},

if $\ell=-(b-a)k/\pi$ . In particular, for $\beta=-1/2$ , we get

\displaystyle W^{[a,b]}_{k,\ell}=\pi e^{i((b-a)k+\pi\ell)/2}J_{0}(((b-a)k+\pi\ell)/2).

Example 2.4.

For $w(x)=(x-a)^{\alpha}(b-x)^{\beta}$ , for $\alpha,\beta\in(-1,0)$ , when the singularity is present at both end-points with different exponents, then

	$\displaystyle W^{[a,b]}_{k,\ell}$	$\displaystyle=\frac{b-a}{2}e^{i((b-a)k+\pi\ell)/2}\int_{-1}^{1}\left(\frac{b-a}{2}\right)^{\alpha+\beta}(1+v)^{\alpha}(1-v)^{\beta}e^{i((b-a)k+\pi\ell)v/2}\,dv$
		$\displaystyle=\frac{(b-a)^{\alpha+\beta+1}}{\Gamma(2+\alpha+\beta)}\Gamma(1+\alpha)\Gamma(1+\beta)_{1}F_{1}(1+\alpha,2+\alpha+\beta,i((b-a)k+\pi\ell)).$

Example 2.5.

For $w(x)=\log(x-a)$ , when the logarithmic singularity is present at the left end-point, then

	$\displaystyle W^{[a,b]}_{k,\ell}=\frac{i(b-a)}{(a-b)k-\pi\ell}($	$\displaystyle\gamma+\Gamma(0,-i((b-a)k+\pi\ell))+(-1+e^{i((b-a)k+\pi\ell)})\log(b-a)$
		$\displaystyle+\log(-i((b-a)k+\pi\ell)),$

if $(b-a)k+\pi\ell\neq 0$ and $W^{[a,b]}_{k,\ell}=(b-a)(\log(b-a)-1)$ otherwise.

3. Error Analysis

We begin by noting that, using the Fourier series of the periodic extension of $f$ , that is,

\displaystyle f(x)=\sum_{\ell=-\infty}^{\infty}c_{\ell}(f^{e})e^{\pi i\ell(x-a)/(b-a)},

where

\displaystyle c_{\ell}(f^{e})=\frac{1}{\xi}\int_{a}^{a+\xi}f^{e}(x)e^{-2\pi i\ell(x-a)/\xi}\,dx=\frac{1}{2(b-a)}\int_{a}^{2b-a}f^{e}(x)e^{-\pi i\ell(x-a)/(b-a)}\,dx,

the exact integral can be expressed as

\displaystyle I^{[a,b]}_{k}(w,f)=e^{ika}\sum_{\ell=-\infty}^{\infty}c_{\ell}(f^{e})W^{[a,b]}_{k,\ell}(w).

The error in the numerical integration, therefore, is given by

	$\displaystyle I^{[a,b]}_{k}(w,f)-I^{[a,b]}_{k,n}(w,f)=e^{ika}\left(\sum_{\ell=-\infty}^{\infty}c_{\ell}(f^{e})W_{k,\ell}^{[a,b]}(w)-\sum_{\ell=-n}^{n-1}d_{\ell,n}(f_{r,q}^{e})W_{k,\ell}^{[a,b]}(w)\right)=$
(13)		$\displaystyle e^{ika}\left(\sum_{\ell=-n}^{n-1}\sum_{\begin{subarray}{c}m=-\infty\\ m\neq 0\end{subarray}}^{\infty}c_{\ell+2mn}(f^{e})W_{k,\ell+2mn}^{[a,b]}(w)+\sum_{\ell=-n}^{n-1}\left(c_{\ell}(f^{e})-d_{\ell,n}(f_{r,q}^{e})\right)W_{k,\ell}^{[a,b]}(w)\right)$

Using

\displaystyle c_{\ell}(f^{e})-d_{\ell,n}(f_{r,q}^{e})=c_{\ell}(f^{e}-f_{r,q}^{e})-\sum_{\begin{subarray}{c}m=-\infty\\ m\neq 0\end{subarray}}^{\infty}c_{\ell+2mn}(f_{r,q}^{e}),

we rewrite the second summation in eq. 13 as

(14)

\displaystyle\sum_{\ell=-n}^{n-1}c_{\ell}(f^{e}-f_{r,q}^{e})W_{k,\ell}^{[a,b]}(w)-\sum_{\ell=-n}^{n-1}\sum_{\begin{subarray}{c}m=-\infty\\ m\neq 0\end{subarray}}^{\infty}c_{\ell+2mn}(f_{r,q}^{e})W_{k,\ell}^{[a,b]}(w),

and, therefore, combining the expressions in eq. 13 and eq. 14, the error reads

	$\displaystyle I^{[a,b]}_{k}(w,f)-I^{[a,b]}_{k,n}(w,f)=e^{ika}\left(\sum_{\ell=-n}^{n-1}\sum_{\begin{subarray}{c}m=-\infty\\ m\neq 0\end{subarray}}^{\infty}c_{\ell+2mn}(f^{e})\left(W_{k,\ell+2mn}^{[a,b]}(w)-W_{k,\ell}^{[a,b]}(w)\right)+\right.$
(15)		$\displaystyle\left.\sum_{\ell=-n}^{n-1}\sum_{\begin{subarray}{c}m=-\infty\\ m\neq 0\end{subarray}}^{\infty}c_{\ell+2mn}(f^{e}-f_{r,q}^{e})W_{k,\ell}^{[a,b]}(w)+\sum_{\ell=-n}^{n-1}c_{\ell}(f^{e}-f_{r,q}^{e})W_{k,\ell}^{[a,b]}(w)\right)$

To see that

(16)

\displaystyle\sum_{\ell=-n}^{n-1}c_{\ell}(f^{e}-f_{r,q}^{e})W_{k,\ell}^{[a,b]}(w)=-\sum_{\ell=-n}^{n-1}\sum_{\begin{subarray}{c}m=-\infty\\ m\neq 0\end{subarray}}^{\infty}c_{\ell+2mn}(f^{e}-f_{r,q}^{e})W_{k,\ell+2mn}^{[a,b]}(w)

we introduce

\displaystyle A(g)(z):=\int_{a}^{2b-a}w_{0}(y)g(z+y)e^{ik(y-a)}\,dy

where $w_{0}$ is the $2(b-a)$ -periodic function with $w_{0}(x)=w(x)$ for $x\in[a,b]$ and $w_{0}(x)=0$ for $x\in(b,2b-a)$ . As

\displaystyle A(g)(z)=\sum_{\ell=-\infty}^{\infty}c_{\ell}(g)e^{\pi i\ell z/(b-a)}W^{[a,b]}_{k,\ell}(w),

we have

\displaystyle\sum_{\ell=-\infty}^{\infty}c_{\ell}(f^{e}-f_{r,q}^{e})W^{[a,b]}_{k,\ell}(w)=A(f^{e}-f_{r,q}^{e})(0)=0.

Therefore, using eq. 16, and after simplification, the error in section 3 is expressed as

(17)

\displaystyle I^{[a,b]}_{k}(w,f)-I^{[a,b]}_{k,n}(w,f)=e^{ika}\sum_{\ell=-n}^{n-1}\sum_{\begin{subarray}{c}m=-\infty\\ m\neq 0\end{subarray}}^{\infty}c_{\ell+2mn}(f_{r,q}^{e})\left(W_{k,\ell+2mn}^{[a,b]}(w)-W_{k,\ell}^{[a,b]}(w)\right).

Thus, to obtain an error bound, we need an estimate on $\left|c_{\ell}(f_{r,q}^{e})\right|$ . In this direction, it is important to note that $f_{r,q}^{e}(x)=f^{e}(x)+n^{-q}g(x)$ where $g(x)=0$ for $x\in[a,b]$ , $g\in C^{2}(b,2b-a)$ and $g(b)=g(2b-a)=0$ . Thus, for $\ell\neq 0$ , we have

(18)

\displaystyle c_{\ell}(f_{r,q}^{e})=c_{\ell}(f^{e})+n^{-q}\frac{b-a}{2(\pi i\ell)^{2}}\left((-1)^{\ell}g^{\prime}(b)-g^{\prime}(2b-a)+\int_{b}^{2b-a}g^{\prime\prime}(x)e^{-\pi i\ell(x-a)/(b-a)}\,dx\right)

and

(19)			$\displaystyle c_{\ell}(f^{e})=\frac{(b-a)^{r+1}}{2(\pi i\ell)^{r+2}}\Bigg{(}(-1)^{\ell}\left(p^{(r+1)}(b)-f^{(r+1)}(b)\right)+f^{(r+1)}(a)-p^{(r+1)}(2b-a)$
		$\displaystyle\quad\quad\quad\quad\quad+\int_{a}^{b}f^{(r+2)}(x)e^{-\pi i\ell(x-a)/(b-a)}\,dx+\int_{b}^{2b-a}p^{(r+2)}(x)e^{-\pi i\ell(x-a)/(b-a)}\,dx\Bigg{)}.$

Therefore, using eq. 18 and eq. 19, we conclude that there exits a positive constant $C$ such that

(20)

\displaystyle\left|c_{\ell}(f_{r,q}^{e})\right|\leq C\left(|\ell|^{-(r+2)}+n^{-q}|\ell|^{-2}\right),\quad\ell\in\mathbb{Z}\setminus\{0\}.

Finally, combining eq. 17 and eq. 20, we have

(21)

\displaystyle\left|I^{[a,b]}_{k}(w,f)-I^{[a,b]}_{k,n}(w,f)\right|\leq Cn^{-(\min\{r,q\}+2)}\sum_{\ell=-n}^{n-1}\sum_{m=-\infty}^{\infty}\frac{\left|W_{k,\ell+2mn}^{[a,b]}(w)-W_{k,\ell}^{[a,b]}(w)\right|}{\left(2|m|-1\right)^{2}}.

We, therefore, see that the convergence rate for the approximations does depend on the weight function $w$ . In this context, using the notation

\displaystyle\mathcal{W}^{[a,b]}_{k,n,M}(w):=\sum_{\ell=-n}^{n-1}\sum_{m=-M}^{M}\frac{\left|W_{k,\ell+2mn}^{[a,b]}(w)-W_{k,\ell}^{[a,b]}(w)\right|}{\left(2|m|-1\right)^{2}}

we introduce a class of weight functions defined as

\displaystyle\mathbb{W}_{\gamma}(a,b):=\left\{w\in L^{1}(a,b):\exists\ n_{0}\in\mathbb{N},\ C>0,\mathcal{W}^{[a,b]}_{k,n,\infty}(w)\leq Cn^{\gamma},\text{ for all }n\geq n_{0}\right\}

for which the approximations converge at the same rate for a given $f$ . Obviously, $\mathbb{W}_{\gamma_{1}}(a,b)\subseteq\mathbb{W}_{\gamma_{2}}(a,b)$ for $\gamma_{1}\leq\gamma_{2}$ . The numerical calculations plotted in fig. 1 indicate that the the weight $w=(x-a)^{\beta},\beta\in(-1,0)$ belongs to $\mathbb{W}_{-\beta}(a,b)$ whereas $w=(x-a)^{\alpha}(b-x)^{\beta},\alpha,\beta\in(-1,0)$ belongs to $\mathbb{W}_{\gamma}(a,b)$ with $\gamma=\max\{-\alpha,-\beta\}$ .

The following theorem captures the essence of the error analysis of the proposed integration scheme.

Refer to caption — (a) $w(x)=x^{\beta},\ k=10$

Theorem 3.1.

Let $w\in\mathbb{W}_{\gamma}(a,b)$ and $f\in C^{\infty}([a,b])$ . Let

\displaystyle I_{k,n}^{[a,b]}(w,f):=e^{ika}\sum_{\ell=-n}^{n-1}d_{\ell,n}(f_{r,q}^{e})W_{k,\ell}^{[0,1]}(w),\quad q,n\in\mathbb{N},\ \ r\in\mathbb{N}\cup\{0\},

where

	$\displaystyle f_{r,q}^{e}(x)=\begin{dcases}f(x),&x\in[a,b],\\ \sum_{m=0}^{r}\mathcal{D}^{m,-}_{n,q}(f)(b)p_{m,r}^{[b,2b-a]}(x)+\sum_{m=0}^{r}\mathcal{D}^{m,+}_{n,q}(f)(a)p_{m,r}^{[2b-a,b]}(x),&x\in(b,2b-a),\end{dcases}$
	$\displaystyle d_{\ell,n}(f_{r,q}^{e})=\frac{1}{2n}\sum_{j=0}^{2n-1}f_{r,q}^{e}(a+(b-a)j/n)e^{-\pi i\ell j/n},\ W^{[a,b]}_{k,\ell}(w)=\int_{a}^{b}w(x)e^{i(k(b-a)+\pi\ell)(x-a)/(b-a)}\,dx,$

be the approximation to the oscillatory integral $I_{k}^{[a,b]}(w,f)$ . Then, there exist $n_{0}\in\mathbb{N}$ and a constant $C>0$ such that

\displaystyle\left|I_{k}^{[a,b]}(w,f)-I_{k,n}^{[a,b]}(w,f)\right|\leq Cn^{-(\min\{r,q\}+2-\gamma)},

for all $n\geq n_{0}$ .

The optimal choice $q=r$ is, therefore, recommended in practice for approximating the derivatives that results in order $r+2-\gamma$ rate of convergence for integral approximations.

4. Numerical Results

In this section, we present a variety of numerical experiments to validate the performance of the proposed method. Unless otherwise stated, we use $q=r$ in quadrature approximations.

Example 4.1.

To begin with, we test our method on an example used by David Levin in [14] and consider the integral

\displaystyle I_{k}=\int_{0}^{1}\sin(t)e^{ik(t+t^{2})}\,dt.

The change of variable $x=t+t^{2}$ yields $I_{k}=\int_{0}^{2}f(x)e^{ikx}\,dx$ where

\displaystyle f(x)=\sin\left(\frac{-1+\sqrt{4x+1}}{2}\right)\frac{1}{\sqrt{4x+1}}

is smooth on the interval $[0,2]$ . Thus, the quadrature with weight function $w(x)=1$ produces rapidly convergent approximations $I_{k,n}$ of $I_{k}$ where the convergence rate can be raised arbitrarily by increasing the parameter $r$ . We studied the convergence at $k=100$ , $k=500$ and $k=1000$ for $r=1,2,3,4,$ and presented the results in fig. 2 where we plot $\log_{10}|I_{k,n}-I_{k}|$ against $\log_{10}(1/n)$ to gauge the speed of convergence through the slope of respective lines. The plots clearly show that approximations converge with the rate $r+2$ .

Example 4.2.

For studying the convergence behaviour of the numerical integration scheme on a second example without any singularity, that is, $w(x)=1$ , we consider the integral

\displaystyle I_{k}=\int_{0}^{1}\sin(\cos t)\sin t\,e^{ik(\cos t)}\,dt

The results in fig. 3 show that the order of convergence $r+2$ does not degrade with increasing oscillation frequency.

Example 4.3.

For the next example, we consider the integral

\displaystyle I_{k}^{p}=\int_{0}^{1}e^{ikt^{p}}\,dt\quad=\frac{1}{p}\int_{0}^{1}x^{-(p-1)/p}e^{ikx}\,dx

that has been studied in literature for measuring the performance of quadratures in dealing with end-point algebraic singularities. For instance, the results for this example is discussed in [15] where, for $k=10000$ and $p=2$ , a relative error of the order $10^{-7}$ is achieved with $n=100$ . For $k=3000$ and $p=10$ , a relative error of $10^{-12}$ is obtained with $n=100$ . Approximations to these integrals using Filon-Clenshaw-Curtis have also been reported in [8] for $p=2/3$ and $p=4/3$ under experiment 1. Our approach, by construction, captures these integrals exactly upto the machine precision through analytic moment computations as seen in table 1 where we report the errors $|I^{p}_{k,2}-I^{p}_{k}|$ for $k=10^{m},\ m=3,4,5,6,7$ and $p=2/3,4/3,2,10$ . The results shown are for $n=2$ and $r=0$ .

$k$	$p=2/3$	$p=4/3$	$p=2$	$p=10$
$10^{3}$	$7.4325\times 10^{-17}$	$1.7110\times 10^{-16}$	$1.2337\times 10^{-16}$	$4.6653\times 10^{-16}$
$10^{4}$	$0.0000\times 10^{+00}$	$2.7730\times 10^{-16}$	$9.7618\times 10^{-17}$	$5.8885\times 10^{-16}$
$10^{5}$	$2.2818\times 10^{-17}$	$2.2485\times 10^{-16}$	$1.5455\times 10^{-16}$	$5.5786\times 10^{-16}$
$10^{6}$	$0.0000\times 10^{+00}$	$2.9916\times 10^{-16}$	$0.0000\times 10^{+00}$	$8.1510\times 10^{-16}$
$10^{7}$	$0.0000\times 10^{+00}$	$6.9014\times 10^{-16}$	$1.3676\times 10^{-16}$	$4.4208\times 10^{-16}$

Table 1. The approximation errors

|I_{k,n}^{p}-I_{k}^{p}|

for the integral in 4.3 with

n=2

and

r=0

Example 4.4.

Next, we study the performance of our method in approximating the integral

\displaystyle I_{k}=\int_{0}^{\pi}e^{ik\sqrt{(\cos t-1)^{2}+\sin^{2}t}}\,dt,

Using the change of variable $x=\sqrt{(\cos t-1)^{2}+\sin^{2}t}$ , we have

\displaystyle I_{k}=\int_{0}^{2}\frac{2}{\sqrt{4-x^{2}}}e^{ikx}\,dx=\int_{0}^{2}\frac{1}{\sqrt{2-x}}\frac{2}{\sqrt{2+x}}e^{ikx}\,dx,

where we have an algebraic singularity of order $\beta=-1/2$ at the upper end point of the integration interval. Moreover, the function $f(x)=1/\sqrt{2+x}$ is smooth on $[0,2]$ . The results in fig. 4 show that the approximations convergence at the expected rate $r+3/2$ .

Example 4.5.

For the following integral

\displaystyle I_{k}=\int_{0}^{1}x^{\beta}\left|x-\frac{1}{2}\right|e^{ikx}\,dx

where $f(x)=|x-1/2|$ is only piecewise smooth, high-order approximation can still be obtained by breaking the integral appropriately as seen in the results reported in fig. 8 where numerical integrals exhibit the convergence rate $r+2+\beta$ .

Example 4.6.

We now consider an integral with a weight of the form $w(x)=(x-a)^{\alpha}(b-x)^{\beta}$ . In particular, we approximate the following integral

\displaystyle I_{k}^{1}=\int_{0}^{1}x^{-1/2}(1-x)^{-1/3}e^{x}e^{ikx}\,dx

using the proposed scheme for $k=10$ and $k=100$ . The convergence rates for $r=1,2,3,4$ , shown in fig. 6, are consistent with the theoretical rate $r+2-\max\{-\alpha,-\beta\}$ .

The the results for a second integral of this type

\displaystyle I_{k}^{2}=\int_{2}^{3}(x-2)^{-1/4}(3-x)^{-2/3}\sin(x)e^{ikx}\,dx,

presented in fig. 7, confirm the convergence rates.

Example 4.7.

Finally, we study the performance of the proposed method in approximating the following integral

\displaystyle I_{k}=\int_{0}^{\pi/2}\log\left(\sqrt{(\cos t-1)^{2}+\sin^{2}t}\right)e^{ik\sqrt{(\cos t-1)^{2}+\sin^{2}t}}\,dt,

that has a logarithmic singularity at $t=0$ . As in 4.4, using the change of variable $x=\sqrt{(\cos t-1)^{2}+\sin^{2}t}$ , yields a linear oscillator

\displaystyle I_{k}=\int_{0}^{\sqrt{2}}\log(x)\frac{2}{\sqrt{4-x^{2}}}e^{ikx}\,dx,

where $w(x)=\log(x)$ and $f(x)=2/\sqrt{4-x^{2}}$ is smooth on $[0,\sqrt{2}]$ .

5. Conclusion

In this paper, we introduce a new Filon-type quadrature for numerical of highly oscillatory integrals of the form eq. 4 that can be employed to approximate more general oscillatory integrals of the form eq. 1 with algebraic singularities and stationary points. The proposed method utilises an interpolating trigonometric polynomial to approximate the smooth and non-oscillatory part of the integrand. Normally, when a function is not periodic, such trigonometric polynomial approximations suffer from the Gibb’s phenomenon. In order that the unwanted Gibb’s oscillations do not adversely impact the accuracy of numerical integration, the method constructs a periodic extension of the underlying function before its approximation. The main advantage of using trigonometric interpolation lies in obtaining a relatively simpler moment problem that can be solved analytically for many important classes of integrands. This approach can easily handle one and two sided integrable algebraic singularities as well as logarithmic singularities, The accompanying numerical analysis of the method reveals that the approximations converge rapidly to the exact integral. The scheme is tested on a variety of highly oscillatory integrals and these numerical experiments confirm that the approximations converge at the rates predicted by the theory.

References

[1] L. N. G. Filon, Iii.—on a quadrature formula for trigonometric integrals, Proceedings of the Royal Society of Edinburgh 49 (1930) 38–47.
[2] Y. L. Luke, On the computation of oscillatory integrals, in: Mathematical Proceedings of the Cambridge Philosophical Society, Vol. 50, Cambridge University Press, 1954, pp. 269–277.
[3] E. Flinn, A modification of filon’s method of numerical integration, Journal of the ACM (JACM) 7 (2) (1960) 181–184.
[4] E. Tuck, A simple “filon-trapezoidal” rule, Mathematics of Computation 21 (98) (1967) 239–241.
[5] A. Iserles, On the numerical quadrature of highly-oscillating integrals i: Fourier transforms, IMA Journal of Numerical Analysis 24 (3) (2004) 365–391.
[6] A. Iserles, S. Nørsett, On quadrature methods for highly oscillatory integrals and their implementations, Tech. rep. (2004).
[7] V. Dominguez, I. Graham, V. Smyshlyaev, Stability and error estimates for filon–clenshaw–curtis rules for highly oscillatory integrals, IMA journal of numerical analysis 31 (4) (2011) 1253–1280.
[8] V. Dominguez, I. G. Graham, T. Kim, Filon–clenshaw–curtis rules for highly oscillatory integrals with algebraic singularities and stationary points, SIAM Journal on Numerical Analysis 51 (3) (2013) 1542–1566.
[9] V. Domínguez, Filon–clenshaw–curtis rules for a class of highly-oscillatory integrals with logarithmic singularities, Journal of computational and applied mathematics 261 (2014) 299–319.
[10] B. Li, S. Xiang, Efficient methods for highly oscillatory integrals with weakly singular and hypersingular kernels, Applied Mathematics and Computation 362 (2019) 124499.
[11] H. Kang, S. Xiang, Efficient integration for a class of highly oscillatory integrals, Applied Mathematics and Computation 218 (7) (2011) 3553–3564.
[12] A. Anand, A. K. Tiwari, A fourier extension based numerical integration scheme for fast and high-order approximation of convolutions with weakly singular kernels, SIAM Journal on Scientific Computing 41 (5) (2019) A2772–A2794.
[13] D. Huybrechs, S. Vandewalle, On the evaluation of highly oscillatory integrals by analytic continuation, SIAM Journal on Numerical Analysis 44 (3) (2006) 1026–1048.
[14] D. Levin, Procedures for computing one-and two-dimensional integrals of functions with rapid irregular oscillations, Mathematics of Computation 38 (158) (1982) 531–538.
[15] S. Zaman, et al., New quadrature rules for highly oscillatory integrals with stationary points, Journal of Computational and Applied Mathematics 278 (2015) 75–89.

Computation of highly oscillatory integrals using a Fourier extension approximation

Abstract.

Key words and phrases:

1. Introduction

2. Methodology

2.1. The basic idea

2.2. A periodic extension:

2.3. The Discrete Problem:

2.4. The moment problem:

Example 2.1.

Example 2.2.

Example 2.3.

Example 2.4.

Example 2.5.

3. Error Analysis

Theorem 3.1.

4. Numerical Results

Example 4.1.

Example 4.2.

Example 4.3.

Example 4.4.

Example 4.5.

Example 4.6.

Example 4.7.

5. Conclusion

References

Computation of highly oscillatory integrals
using a Fourier extension approximation