Fast Computation of Integrals with Fourier-Type Oscillator Involving Stationary Point

Abstract

An adaptive splitting algorithm was implemented for numerical evaluation of Fourier-type highly oscillatory integrals involving stationary point. Accordingly, a modified Levin collocation method was coupled with multi-resolution quadratures in order to tackle the stationary point and irregular oscillations of the integrand caused by $\omega$. Some test problems are included to verify the accuracy of the proposed methods.

1. Introduction

Accurate and efficient evaluation of highly oscillatory integrals is challenging as the analytical and classical computational methods fail to compute the integrals. The interest of computational scientists to evaluate these integrals quickly and accurately is developed due to the wide range of applications of these integrals in different fields of science and engineering such as optics, acoustics, quantum mechanics, seismology image processing, and electromagnetic [1,2,3,4]. Generally, these integrals can be written as

$$I[r,\omega ]={\int }_{-1}^{1}r\left(x\right)e^{i\omega \mathsf{\Theta }\left(x\right)}\mathrm{d}x,$$

(1)

where $r,\mathsf{\Theta }$ are both smooth and non-oscillatory functions over the interval $[-1,1]$, $\omega \gg 1$ and the function $\mathsf{\Theta }\left(x\right)$ has a stationary point of order k at $x=x_0$ i.e., ${\mathsf{\Theta }}^{\prime }\left(x_0\right)={\mathsf{\Theta }}^{″}\left(x_0\right)=\dots ={\mathsf{\Theta }}^{(k-1)}\left(x_0\right)=0$ and ${\mathsf{\Theta }}^{\left(k\right)}\left(x_0\right)\ne 0$. Large values of the frequency $\omega$ cause a highly oscillations of the integrand, which make the existing quadratures and softwares like MATHEMATICA in some cases, inaccurate. New algorithms need to be developed for accurate and efficient computation of these integrals.

In the last two decades, a number of accurate and efficient methods have been designed for numerical evaluation of one-dimensional highly oscillatory integrals, which include: the asymptotic method [5,6,7,8], the numerical steepest descent method [9], the Filon(-type) methods [10,11], and the Levin(-type) methods [12,13,14,15,16]. Among which, the Levin method has attracted much attention as it can handle highly oscillatory integrals with complicated phase functions. The asymptotic expansion theory which is considered in [10] is also one of the best methods for evaluation of highly oscillatory integrals. Since then, many authors [11,17,18] have implemented this method.

The Filon method [11] gives better approximation when the frequency $\omega ⟶\infty$, but the method has its limitations, that is, it is only applicable for linear phase functions. Compared to the Filon method, Levin’s method has the advantage of solving oscillatory integrals with complicated phase functions. According to the Levin method, a one-dimensional oscillatory integral is transformed into an ordinary differential equation (ODE) and then the ODE can be solved by the collocation method. Subsequently, a solution of such an integral is obtained.

In [13], the authors presented an improved Levin method with Chebyshev polynomials as a basis function and better accuracy was obtained. In [15], the authors evaluated non-oscillatory, mildly oscillatory, and highly oscillatory multi-dimensional integrals by using a meshless method based on a multi-quadric radial basis function and multi-resolution quadratures. In [14], the authors presented a meshless method based on a multi-quadric radial basis function (MQ-RBF) and multi-resolution quadrature to solve highly oscillatory integrals with and without critical points. In [19], the author split the integration domain into two or more sub-domains in order to separate the critical point and a two-point Gauss–Legendre quadrature was used to evaluate the integrals with critical points. The remaining integrals were computed by the modified Levin collocation method.

In [16], the authors developed an algorithm in which the meshless method with a multi-quadric radial basis function was coupled with the multi-resolution quadrature based on hybrid functions and the Haar wavelet to evaluate highly oscillatory integrals with critical point(s). The multi-resolution quadrature based on hybrid function of order 8 to evaluate integral of the form $I\left[f\right]={\int }_a^{b}f\left(x\right)\mathrm{d}x$, is given by

$$\begin{array}{c}{Q}_h^8\left[f\right]=\frac{h}{241920}\sum _{k=1}^n\left[295627f\left(a+\frac{h}{2}(16k-15)\right)+71329f\left(a+\frac{h}{2}(16k-13)\right)\right.\hfill \\ +471771f\left(a+\frac{h}{2}(16k-11)\right)+128953f\left(a+\frac{h}{2}(16k-9)\right)\\ +128953f\left(a+\frac{h}{2}(16k-7)\right)+471771f\left(a+\frac{h}{2}(16k-5)\right)\\ \hfill +71329f\left.\left(a+\frac{h}{2}(16k-3)\right)+295627f\left(a+\frac{h}{2}(16k-1)\right)\right],\end{array}$$

(2)

where $h=\frac{b-a}{8n}$.

For a = 0, b = 1, and n = 4, the truncation error bound of the hybrid function $Q_h^8\left[f\right]$ is given

$$E=\frac{-3194621\times h^9}{14515200}f{\left({\xi }_1\right)}^{\left(8\right)},$$

(3)

for some ${\xi }_1$ in (a, b) [20].

The formula of the Haar wavelet based quadrature for computing the same integral is given as

$$Q_H^w\left[f\right]=h\sum _{i=1}^{n}f\left(x_i\right)=h\sum _{i=1}^{n}f\left(a+h(i-0.5)\right),$$

(4)

where $h=\frac{b-a}{n}$ and $n=2M$. The truncation error of the quadrature rule $Q_H^w\left[f\right]$ is given by

$$\left|E\right|=\frac{h^3}{6}{f}^{″}\left({\varsigma }_1\right),$$

for some ${\varsigma }_1\in (a,b)$ [20].

The current work is the extension of the work reported in [16,19]. According to the proposed procedure, the interval of integration is divided into two or more subintervals in order to isolate the stationary point. Integrals defined over a stationary point oriented interval can be computed by the multi-resolution quadratures even for higher frequencies, while the remaining integrals with high oscillations can be approximated by the Levin collocation quadrature. Theoretical error analysis of the individual methods and the proposed algorithm was performed. Numerical results met the theoretical proofs of the proposed algorithm.

2. Evaluation Procedure

Levin type quadratures fail to compute Fourier-type oscillatory integrals (1) containing a stationary point. We propose a splitting algorithm, which couples two types of quadrature rules to tackle the stationary point. In this algorithm, we subdivide the interval $[-1,1]$ in such a manner that the stationary point is separated in a very small interval and can be handled by the multi-resolution quadratures. The evaluation procedure is discussed as follows.

2.1. Levin Quadrature

The Levin quadrature is applicable to evaluate Fourier-type oscillatory integrals of the form

$$I[r,\omega ]={\int }_a^{b}r\left(x\right)e^{i\omega \mathsf{\Theta }\left(x\right)}\mathrm{d}x,$$

(5)

where ${\mathsf{\Theta }}^{\prime }\left(x\right)\ne 0$. To find the approximate solution of the integral (5), an approximate function $\widehat{S}\left(x\right)$ is assumed to satisfy the following differential equation

$$\mathcal{L}\left(S\left(x\right)e^{i\omega \mathsf{\Theta }\left(x\right)}\right)=r\left(x\right)e^{i\omega \mathsf{\Theta }\left(x\right)},x\in {\mathbb{R}}^n,n=1,2,3,\dots ,$$

(6)

where $\mathcal{L}$ is defined as

$$\mathcal{L}=\left\{\begin{array}{c}\frac{d}{dx}\left[\widehat{S}\left(x\right)e^{i\omega \mathsf{\Theta }\left(x\right)}\right],\mathrm{for}n=1\hfill \\ \frac{{\partial }^n\left[\widehat{S}\left(x\right)e^{i\omega \mathsf{\Theta }\left(x\right)}\right]}{\underset{\mathrm{n}\text{-}\mathrm{variable}}{\underset{⏝}{\partial x_1\partial x_2\dots }}},\mathrm{for}n\ge 2.\hfill \end{array}\right.$$

On applying the interpolation condition to (6), one can find the approximate solution $\widehat{S}\left(x\right)$. Consequently, the desired solution of integral (5) can be obtained as

$$\begin{array}{cc}\hfill I(r,\omega )& ={\int }_a^b\frac{d}{dx}\left[\widehat{S}\left(x\right)e^{i\omega \mathsf{\Theta }\left(x\right)}\right]\hfill \\ & =\widehat{S}\left(b\right)e^{i\omega \mathsf{\Theta }\left(b\right)}-\widehat{S}\left(a\right)e^{i\omega \mathsf{\Theta }\left(a\right)}.\hfill \end{array}$$

2.2. Chebyshev Differentiation Matrix and Its Approximation

The Chebyshev differentiation matrix D [13] is used to approximate the derivative of a function and can be explicitly defined as

$${\mathbf{D}}_{kj}=\left\{\begin{array}{c}\frac{v_k}{v_j}.\frac{{(-1)}^{k+j}}{(x_k-x_j)},k,j=0,\dots ,N-1,k\ne j\hfill \\ -{\mathsf{\Sigma }}_{n=0,n\ne k}^{N-1},D_{kn},k=j,\hfill \end{array}\right.$$

(7)

where

$$v_j=\left\{\begin{array}{c}2,j=0,N-1\hfill \\ 1,j=1,\dots ,N-2,\hfill \end{array}\right.$$

(8)

and $x_j=cos\left(\frac{\pi j}{N-1}\right),j=0,1,2,\dots ,N-1$ are the N Chebyshev–Gauss–Lobatto nodes. If an integral is defined over any domain $[a,b]$, then, the Chebyshev–Gauss–Lobatto nodes $x_j$ can be obtained by the transformation equation $x_j=\frac{b-a}{2}cos\left(\frac{\pi j}{N-1}\right)+\frac{a+b}{2},j=0,1,\dots ,N-1$.

If $\mathbf{u}$ is a vector of function values of $f\left(x\right)$ at Chebyshev–Gauss–Lobatto nodes, then the approximate values of $f^{\prime }\left(x\right)$ at Chebyshev–Gauss–Lobatto nodes can be obtained by $\mathbf{Du}$. The procedure is applied to compute the first derivative of the following oscillatory type functions. The results in terms of absolute errors are shown in Figure 1.

$$f_1\left(x\right)=cos\left(\omega x^{10}\right),f_2\left(x\right)=sinh\left(5x\right).$$

(9)

From Figure 1, it is shown that the accuracy for the derivative of the functions approximated by the Chebyshev differentiation matrix $\mathbf{D}\mathbf{u}$ at Chebyshev–Gauss–Lobatto nodes improves on increasing N, and decaying for large frequency $\omega$. In the Levin procedure, we are concerned with the derivative of a smooth function $S\left(x\right)$, which is independent of $\omega$.

2.3. Chebyshev–Levin Quadrature

The Chebyshev differentiation matrix is used to find the approximate solution $\widehat{S}\left(x\right)$ of ODE (6). If $\mathsf{\Theta }\left(x\right)$ is any phase function of integral (1) such that ${\mathsf{\Theta }}^{\prime }\left(x\right)\ne 0$ for all $x\in [a,b]$, then the discretized form of the ODE (6) at Chebyshev–Gauss–Lobatto nodes is given by

$${\widehat{\mathbf{S}}}^{\prime }+i\omega {\mathsf{\Theta }}^{\prime }\odot \widehat{\mathbf{S}}=\mathbf{R},$$

(10)

where ⊙ represents Hadamard product of the real valued function ${\mathsf{\Theta }}^{\prime }$ and the vector $\widehat{\mathbf{S}}$. If ${\widehat{\mathbf{S}}}^{\prime }=\mathbf{D}\widehat{\mathbf{S}}$, then the descretized form ODE (10) is given as

$$(\mathbf{D}+i\omega \sum )\widehat{\mathbf{S}}=\mathbf{R},$$

(11)

where $\sum =diag\left({\mathsf{\Theta }}^{\prime }\right)$ and can be written in Matlab built-in form as $\mathsf{\Gamma }={\mathsf{\Theta }}^{\prime }.\ast \mathbf{E}$. Particularly if $\mathsf{\Theta }\left(x\right)=x$, Equation (11) can be written as

$$(\mathbf{D}+i\omega \mathbf{I})\widehat{\mathbf{S}}=\mathbf{R},$$

(12)

where $\mathbf{I}$ is the identity matrix and each of $\widehat{\mathbf{S}}$ and $\mathbf{R}$ are column matrices of order $N\times 1$. The vectors $\widehat{\mathbf{S}}$, $\mathsf{\Theta }$ and $\mathbf{R}$ are obtained from the real valued functions $\widehat{S}\left(x\right)$, $\mathsf{\Theta }\left(x\right)$ and $r\left(x\right)$, respectively, at Chebyshev–Gauss–Lobatto nodes as

$$\begin{array}{c}\widehat{\mathbf{S}}={[s\left(x_0\right)s\left(x_1\right)\dots s\left(x_{N-1}\right)]}^T,\mathsf{\Theta }={[\mathsf{\Theta }\left(x_0\right)\mathsf{\Theta }\left(x_1\right)\dots \mathsf{\Theta }\left(x_{N-1}\right)]}^T,\hfill \\ \mathrm{and}\mathbf{R}={[r\left(x_0\right)r\left(x_1\right)\dots r\left(x_{N-1}\right)]}^T.\hfill \end{array}$$

Equation (12) can be written in matrix form as

$$\mathbf{A}\widehat{\mathbf{S}}=\mathbf{R},$$

(13)

where $\mathbf{A}$ is a square matrix of order $N\times N$ and both column vectors $\widehat{\mathbf{S}}$ and $\mathbf{R}$ have length N. Then the system of Equation (13) for the unknown vector $\widehat{\mathbf{S}}\left(x\right)$ is solved. Subsequently, the desired numerical solution of the oscillatory integral (5) can be obtained by the Levin’s procedure as

$$Q_{C-L}\left[r\right]=\widehat{\mathbf{S}}\left(x_{N-1}\right)e^{i\omega \mathsf{\Theta }\left(x_{N-1}\right)}-\widehat{\mathbf{S}}\left(x_0\right)e^{i\omega \mathsf{\Theta }\left(x_0\right)}.$$

2.4. Adaptive Splitting

The Chebyshev–Levin quadrature $Q_{C-L}\left[r\right]$ fails to evaluate the Fourier-type oscillatory integral (1) with a stationary point at $x=x_0$. In this section, an adaptive splitting is used to evaluate the integral (1). According to this splitting, the domain interval is bifurcated into two or more sub-intervals in such manner that the stationary point is isolated in a very small interval. For this, a parameter $\zeta ={\left(\frac{n}{10\omega }\right)}^{\frac{1}{k}}$, a point in the vicinity of stationary point $x_0$, is defined such that for fixed n, $\zeta \to 0$ as $\omega \to \infty$.

Case-I When the stationary point $x_0=a$ lies at the left end point of the domain, then integral (1) can be split as

$$I[r,\omega ]={\int }_{x_0}^{\zeta }r\left(x\right)e^{i\omega \mathsf{\Theta }\left(x\right)}\mathrm{d}x+{\int }_{\zeta }^{b}r\left(x\right)e^{i\omega \mathsf{\Theta }\left(x\right)}\mathrm{d}x=I_1[r,\omega ]+I_2[r,\omega ]$$

(14)

The integral $I_1[r,\omega ]$ is computed by the quadrature rule $Q_h^8\left[f\right]$ or $Q_H^w\left[f\right]$ and the second integral $I_2[r,\omega ]$ is approximated by the proposed method $Q_{C-L}\left[r\right]$. The following transformation equation is used to discretize the interval $[\zeta ,b]$ into Chebyshev–Gauss–Lobatto nodes.

$$x_j=\frac{(b-\zeta )cos\left(\frac{j\pi }{N-1}\right)}{2}+\frac{\zeta +b}{2},j=0,1,\dots ,N-1.$$

(15)

As the stationary point occurs anywhere in the interval $[a,b]$ and hence for this, we have to discuss the following case as case-II.

Case-II When $a<x_0<b$, then the split form of the integral (1) is given as

$$\begin{array}{cc}\hfill I[r,\omega ]& ={\int }_a^{x_0-\zeta }r\left(x\right)e^{i\omega \mathsf{\Theta }\left(x\right)}\mathrm{d}x+{\int }_{x_0-\zeta }^{x_0+\zeta }r\left(x\right)e^{i\omega \mathsf{\Theta }\left(x\right)}\mathrm{d}x+{\int }_{x_0+\zeta }^{b}r\left(x\right)e^{i\omega \mathsf{\Theta }\left(x\right)}\mathrm{d}x\hfill \\ & =I_1[r,\omega ]+I_2[r,\omega ]+I_3[r,\omega ].\hfill \end{array}$$

(16)

The integral $I_2[r,\omega ]$ contains the stationary point and is computed by the multi-resolution quadratures $Q_h^8\left[f\right]$ or $Q_H^w\left[f\right]$ with n quadrature points. The integrals containing no stationary point are computed by $Q_{C-L}\left[r\right]$ with N collocation points. The following transformation Equation (15) for integral $I_3[r,\omega ]$ is used to descretize the interval in Chebyshev–Gauss–Lobatto nodes.

$$x_j=\frac{(b-x_0-\zeta )cos\left(\frac{j\pi }{N-1}\right)}{2}+\frac{x_0+\zeta +b}{2},j=0,1,\dots ,N-1,$$

and similarly we modify the transformation Equation (15) for $I_1[r,\omega ]$. The final values of the component integrals are combined as

• If $Q_h^8\left[f\right]$ is used to compute the integral having a stationary point, then Chebyshev–hybrid quadrature is given by

  $$ChQ=Q_h^8\left[f\right]+Q_{C-L}\left[r\right].$$
• If $Q_H^w\left[f\right]$ is used for computing the integral having stationary point, then Chebyshev–Haar quadrature can be written as

  $$CHQ=Q_H^w\left[f\right]+Q_{C-L}\left[r\right].$$

Algorithm-I:

• $x_j=cos\left(\frac{i\ast \pi }{N-1}\right);j=0,1,\dots ,N-1$;
• $\xi ={\left(\frac{N_0}{10\ast \omega }\right)}^{1/k}$; (Slitting parameter.)
• $\left[\mathbf{D}\right]=Cheb.diff.matrix(N,1,2)$;
• $\mathbf{R}=r\left(x_j\right)$, for $j=1,2,3,\dots ,N$; (where r(x) is the amplitude function of (1).)
• $\mathbf{B}=\mathbf{D}+i\ast \omega \ast \mathbf{I}(N,N)$; (The coefficient matrix of the ODE (8) and I is the identity matrix.)
• $\widehat{\mathbf{p}}=inv\left(\mathbf{B}\right)\ast {\mathbf{R}}^{\prime }$; (The approximate solution of the ODE (8).)
• $App1={\widehat{\mathbf{p}}}_{N-1}\ast exp(i\ast \omega \ast \mathsf{\Theta }\left(x_{N-1}\right))-{\widehat{\mathbf{p}}}_{x_0}\ast exp(i\ast \omega \ast \mathsf{\Theta }\left(x_0\right))$; (Approximate Chebshev solution of integral having no stationary point.)
• $App2=hybrid8(a,b,\omega ,n)$; (Approximate hybrid solution of the integral having stationary point.)
• ChQ = App1 + App2; (Solution of integral (1)).

3. Error Bounds

In this section, error bounds of the individual methods as well as the proposed Algorithm-I in the inverse power of $\omega$ are obtained to ensure the asymptotic convergence rate of the new algorithm.

Lemma 1

([19]). Suppose that $\mathsf{\Theta }\left(x\right)$ is a real valued and smooth function in $(a,b)$ and $|{\mathsf{\Theta }}^{\left(k\right)}\left(x\right)|\ge 1$ for all $x\in (a,b)$ and a fixed value of k. Then $\left|{\int }_a^b{e}^{i\omega \mathsf{\Theta }\left(x\right)}dx\right|\le c\left(k\right){\omega }^{-1/k}$ holds when:

i. 

$k\ge 2$ or

ii. 

$k=1$ and ${\mathsf{\Theta }}^{\prime }\left(x\right)$ is monotonic.

The bound $c\left(k\right)$ is independent of Θ and ω, and $c\left(k\right)={5.2}^{k-1}-2.$

Lemma 2

([3,19]). Under the assumption on $\mathsf{\Theta }\left(x\right)$ in Lemma 1, the following result is obtained

$$\left|{\int }_a^b{e}^{i\omega \mathsf{\Theta }\left(x\right)}F\left(x\right)dx\right|\le c\left(k\right){\omega }^{\frac{-1}{k}}\left[\mid F\left(b\right)\mid +{\int }_a^b\mid F^{\prime }\left(x\right)\mid dx\right].$$

Lemma 3.

Suppose that the integral $I_2(r,\omega )$ of (14) has a unique stationary point at $x=x_0$ which lies at the left end of the domain. Then the error bound of the Chebyshev–Levin quadrature $Q_{C-L}\left[r\right]$ with N collocation points for computing $I_2$ is given by

$$|{\mathcal{E}}_2|=|I_2(r,\omega )-Q_{C-L}\left[r\right]|=O\left(\frac{{(b-\zeta )}^N}{{\omega }^{1/k}}\right),$$

where k is the order of stationary point.

Proof. 

Let $\widehat{p}\left(x\right)={\sum }_{i=0}^{N-1}{\lambda }_i{T}_i\left(x\right)$ be the approximate Chebyshev interpolating polynomial of the first kind of order $N-1$ that satisfies the following ODE

$$p^{\prime }\left(x\right)+i\omega {\mathsf{\Theta }}^{\prime }\left(x\right)p\left(x\right)=r\left(x\right).$$

Here ${\lambda }_i,i=0,1,\dots N-1$ are the N unknown coefficients, which can be determined by the following interpolation condition

$$\psi \left(x_k\right)=r\left(x_k\right),k=0,1,\dots N-1,$$

where $\psi \left(x\right)={\widehat{p}}^{\prime }\left(x\right)+i\omega {\mathsf{\Theta }}^{\prime }\left(x\right)\widehat{p}\left(x\right)$. By applying the recurrence relation $T_i\left(x\right)=2x{T}_{i-1}\left(x\right)-T_{i-2}\left(x\right),i=2,3,\dots$ [20] with $T_0=1$ and $T_1=x$, it follows that the function $\psi \left(x\right)$ can be written as

$$\psi \left(x\right)=\sum _{i=2}^{N-1}{\lambda }_i\left[2x{T}_{i-1}^{\prime }\left(x\right)-T_{i-2}^{\prime }\left(x\right)+(2+i\omega x{\mathsf{\Theta }}^{\prime }\left(x\right))T_{i-1}\left(x\right)-i\omega {\mathsf{\Theta }}^{\prime }\left(x\right)T_{i-2}\left(x\right)\right].$$

Let $\mathsf{\Phi }\left(x\right)=Ma{x}_{x\in [a,b]}(r\left(x\right)-\psi \left(x\right))$. Using Lemma 2 and [16], the error bound of the integral $I_2(r,\omega )$ of (14) computed by the $Q_{C-L}\left[r\right]$ is given as

$$\begin{array}{cc}\hfill |{\mathcal{E}}_2|& =|I_2(r,\omega )-Q_{C-L}\left[r\right]|={\int }_{\zeta }^b(r\left(x\right)-\psi \left(x\right))e^{i\omega \mathsf{\Theta }\left(x\right)}\mathrm{d}x\le {\int }_{\zeta }^b\mathsf{\Phi }\left(x\right)e^{i\omega \mathsf{\Theta }\left(x\right)}\mathrm{d}x\hfill \\ & \le \frac{3(b-\zeta +N)\parallel {\mathsf{\Phi }}^{\left(m\right)}{\parallel }_{\infty }{(b-\zeta )}^N}{N!\omega \delta }.\hfill \end{array}$$

(17)

Following the procedure [16] and (17), we can find the error bound of the proposed method $Q_{C-L}\left[r\right]$ as

$$|{\mathcal{E}}_2|=|I_2(r,\omega )-Q_{C-L}\left[r\right]|=O\left(\frac{{(b-\zeta )}^N}{{\omega }^{1/k}}\right).$$

□

Lemma 4.

Suppose that the integrand of $I_1(r,\omega )$ of (14) has a unique stationary point at $x=x_0$ that lies at the left endpoint of the domain interval. Then the error bound of $Q_h^8\left[f\right]$ for computing $I_1(r,\omega )={\int }_{x_0}^{\zeta }r\left(x\right)e^{i\omega \mathsf{\Theta }\left(x\right)}\mathrm{d}x$ with n quadrature points is given by

$$|{\mathcal{E}}_1|=\left|I_1(r,\omega )-Q_h^8\left[f\right]\right|=O\left({\omega }^{-(1+1/k)}\right),$$

where $n=\frac{\xi -x_0}{8h}$.

Proof. 

Since the error bound of the hybrid function $Q_h^8\left[f\right]$ with n quadrature points is given by (See Theorem 1, [16])

$$|{\mathcal{E}}_1|\le \frac{C({\left(\frac{n}{10\omega }\right)}^{1/k}-a)}{4.54\times {10}^8},\omega \gg 1,$$

where a is the stationary point, C is a constant independent of $\omega$, h and k is the order of stationary point. Since $Q_h^8\left[f\right]$ can compute the integral $I_1(r,\omega )$ having stationary point at $x=x_0$. Therefore, the error bound of the hybrid function $Q_h^8\left[f\right]$ for computing $I_1(r,\omega )={\int }_{x_0}^{\zeta }r\left(x\right)e^{i\omega \mathsf{\Theta }\left(x\right)}\mathrm{d}x$, (where $\zeta$ is the nearest point of $x_0$) is given as

$$|{\mathcal{E}}_1|=\left|I_1(r,\omega )-Q_h^8\left[f\right]\right|\le \frac{C_1(\frac{n}{10\omega }-x_0)}{4.54\times {10}^8},\omega \gg 1.$$

(18)

If $x_0=0$, then (18) becomes

$$|{\mathcal{E}}_1|=\left|I_1(r,\omega )-Q_h^8\left[f\right]\right|\le C_1\left(\frac{n}{45.4\times {10}^8\omega }\right),$$

(19)

where $C_1$ is independent of h and $\omega$.

Since $n=\frac{\zeta }{8h}$, the inequality in (19) can be written as

$$|{\mathcal{E}}_1|\le \frac{C_1\zeta }{363.2\times {10}^8\omega h}.$$

(20)

Using $\zeta ={\left(\frac{n}{10\omega }\right)}^{\frac{1}{k}}$, Equation (20) can be written as

$$\begin{array}{cc}\hfill |{\mathcal{E}}_1|& =\frac{C_2{n}^{1/k}}{h{\omega }^{1+1/k}},\mathrm{where}{C}_2=\frac{C_1}{363.2\times {10}^{8+1/k}}\hfill \\ & =O\left({\omega }^{-(1+1/k)}\right).\hfill \end{array}$$

□

As the Chebyshev–hybrid quadrature splitting method (ChQ) is obtained as

$$ChQ=Q_h^8\left[f\right]+Q_{C-L}\left[r\right],$$

therefore, the error bound of the splitting method ChQ can be obtained as

$$|{\mathcal{E}}_1|=|I[r,\omega ]-ChQ|\le \min \left(\frac{C_2{n}^{1/k}}{h{\omega }^{1+1/k}},\frac{C{(b-\zeta )}^N}{{\omega }^{1/k}}\right),$$

(21)

where $C_2$ and C are independent of h and $\omega$. For each subinterval of size h, Equation (21) can be written as

$$|{\mathcal{E}}_1|=|I[r,\omega ]-ChQ|\le O\left(\min \left\{\frac{n^{1/k}}{h{\omega }^{1+1/k}},\frac{h^{N-1}}{{\omega }^{1/k}}\right\}\right).$$

4. Numerical Examples and Discussion

In this section, the proposed method was tested on solving some benchmark problems [11,19]. The reference solution was obtained by MAPLE 18. The absolute errors $L_{abs}$ and scaled absolute errors were computed in each test problem. Results of the new methods were compared with methods reported in [11,19]. MATLAB 2009a platform was used for computation of numerical results.

Example 1.

Consider the following integral [11]

$$I_1[r,\omega ]={\int }_0^{1}sinhx{e}^{i\omega (x^3+x^2+x)}\mathrm{d}x.$$

(22)

Integral (22) is highly oscillatory having no stationary point and can be computed by the Chebyshev–Levin method $Q_{C-L}\left[r\right]$. The absolute errors are shown in Figure 2 and Figure 3 and

Table 1. $L_{abs}$ produced by the $Q_{C-L}\left[r\right]$, $(n,N=10)$ for test Problem 1.

$\mathit{\omega }$	${\mathit{Q}}_{\mathit{C}-\mathit{L}}\left[\mathit{r}\right]$	${\mathit{Q}}_{\mathit{h}}^8\left[\mathit{f}\right]$	${\mathit{Q}}_{\mathit{H}}^{\mathit{w}}\left[\mathit{f}\right]$
${10}^5$	$4.37\times {10}^{-16}$	$7.14\times {10}^{-2}$	$7.12\times {10}^{-2}$
${10}^6$	$7.39\times {10}^{-18}$	$5.49\times {10}^{-2}$	$1.11\times {10}^{-1}$
${10}^7$	$5.86\times {10}^{-20}$	$5.48\times {10}^{-2}$	$2.02\times {10}^{-1}$
${10}^8$	$8.25\times {10}^{-22}$	$2.23\times {10}^{-2}$	$2.77\times {10}^{-2}$
${10}^9$	$6.10\times {10}^{-24}$	$2.80\times {10}^{-2}$	$3.61\times {10}^{-2}$

. Figure 2a shows that the asymptotic order of convergence of the proposed method reached to $O\left({\omega }^{-3.5}\right)$, while the asymptotic order of convergence of the methods reported in [11] is $O\left({\omega }^{-2}\right)$ as shown in Figure 2b. The proposed method was tested for higher frequencies. The absolute errors are shown in Table 1. It is shown in the table that the proposed method $Q_{C-L}\left[r\right]$ improved accuracy on increasing $\omega$.

Integral (22) was computed by the proposed methods $Q_{C-L}\left[r\right]$ for fixed frequency and varying nodal points and the results were compared with those obtained by $Q_h^8\left[f\right]$ and $Q_H^w\left[f\right]$. The results are shown in Figure 3a, while in Figure 3b, the nodal points are fixed and the frequencies vary. It is shown in the figure that the proposed method was better than all the other methods. The main advantage of the proposed method $Q_{C-L}\left[r\right]$ was that it improved accuracy on increasing $\omega$ as well as N.

Example 2.

Consider the following integral [19]

$$I_2[r,\omega ]={\int }_0^1{e}^{i\omega x^{10}}\mathrm{d}x.$$

(23)

Oscillatory behavior of the real part of the integrand (23) is shown in Figure 4a. The integral has a stationary point of order $k=10$ at $x=0$, which lies at the extreme left position of the domain. The Chebyshev–Levin quadrature $Q_{C-L}\left[r\right]$ failed to evaluate the integral due to the existence of stationary point in the domain interval. Multi-resolution quadrature $Q_h^8\left[f\right]$ or $Q_H^w\left[f\right]$ needed dense nodes to give the desired accuracy for high frequencies, which was impractical. The integral was computed by the new splitting procedure ChQ accurately. Results in the form of scaled absolute errors and relative errors were computed and are shown in Figure 5. In Figure 5a, it has been shown that the asymptotic order of convergence of the splitting method ChQ reached to $O\left({\omega }^{-3}\right)$. The new splitting method improved accuracy on increasing nodal points as shown in Figure 5b.

The two splitting methods ChQ and Chebyshev–Haar quadrature (CHQ) were implemented for fixed nodes and varying frequencies. The results are shown in Figure 6a. The new methods were compared with the result of the meshless splitting method [16] in Figure 6b. From the figures, it is shown that the new splitting methods were better than the methods reported in [16]. From all figures, it is evident that the new splitting methods were more accurate than the existing methods even for smaller nodes.

Example 3.

Consider the following integral [11]

$$I_3[r,\omega ]={\int }_{-1}^1\frac{1}{x+2}{e}^{i\omega (1-cosx-\frac{1}{2}{x}^2+x^3)}\mathrm{d}x.$$

(24)

The integral in (24) is highly oscillatory. Irregular high oscillations of real part of the integrand are shown in Figure 4b. The integral has a stationary point of order 2 at $x=0\in [-1,1]$. The integral was split according to the case-II of the splitting methods and was computed by the two new splitting methods ChQ and CHQ. Absolute errors and scaled absolute errors are analyzed in

Table 2. $L_{abs}$ produced by ChQ for test Problem 3.

y	$\mathit{m}=10$	$\mathit{m}=20$	$\mathit{m}=30$	$\mathit{m}=40$
${10}^5$	$6.9845\times {10}^{-8}$	$2.5559\times {10}^{-10}$	$1.8400\times {10}^{-13}$	$2.0019\times {10}^{-15}$
${10}^6$	$2.7704\times {10}^{-10}$	$1.3929\times {10}^{-12}$	$4.0003\times {10}^{-15}$	mp
${10}^7$	$6.1799\times {10}^{-12}$	$3.0001\times {10}^{-14}$	mp	mp

and Figure 7. Figure 7a indicates that the new splitting method ChQ retained the same asymptotic order of convergence $O\left({\omega }^{-3}\right)$ for this problem as well. In Figure 7b, results of the new splitting methods were compared with multi-resolution quadratures for fixed frequency and varying nodal points. The new splitting methods reduced the error up to $O\left({10}^{-15}\right)$, which was a higher rate of convergence than the method reported in [11].

Method ChQ was tested for higher frequencies. The results are shown in Table 2. It is shown in the table that the new methods improved the results on increasing $\omega$ as well as N.

Example 4.

Consider the following integral [11]

$$I_4[r,\omega ]={\int }_0^1{e}^x{e}^{i\omega (x^{3/2}+x^{5/2})}\mathrm{d}x.$$

(25)

The integral in (25) has a stationary point at $x=0$. The integral was computed by the new splitting algorithms ChQ and CHQ. Results are analyzed in Figure 8 and Figure 9. The result in terms of scaled absolute errors is shown in Figure 8a. In the same figure, results of the proposed methods are compared with the results of Filon-type method [11]. It has been shown that the asymptotic order of convergence of the new method ChQ reached to $O\left({\omega }^{-3}\right)$, while the method reported in [11] was of order $O\left({\omega }^{-2}\right)$.

Comparison of the proposed methods was performed for fixed $\omega$ and varying nodal points and vice versa as shown in Figure 9. It has been shown that the proposed splitting algorithms ChQ and CHQ improved the accuracy on increasing $\omega$ and N.

5. Conclusions

An adaptive splitting procedure which couples Chebyshev–Levin quadrature $Q_{C-L}\left[r\right]$ with multi-resolution quadratures was implemented to evaluate highly oscillatory integrals with a stationary point. Theoretical error analysis in terms of inverse powers of $\omega$ of the proposed procedure was performed and numerically verified by solving a few benchmark problems. The asymptotic order of convergence of both, the Chebyshev–Levin quadrature $Q_{C-L}\left[r\right]$ and the splitting algorithm reached to $O\left({\omega }^{-3}\right)$. Results of the proposed algorithm were compared with some existing methods reported in the literature. Moreover, the improvement in the accuracy could easily be witnessed from the results of new algorithm.