New Algorithms to Compute Hypersingular Integrals with Rapidly Oscillating Kernels

The paper is dedicated to the study of two new algorithms for approximating hypersingular integrals with rapidly oscillating kernels. The methods use an interpolatory procedure at zeros of orthogonal polynomials. Estimates of the error are given, as well as of the amplification factors. Some numerical examples show the coherence of the theoretical results with the numerical ones.

Introduction

In [1] the authors propose an approximate method for evaluating Cauchy singular integral with rapidly oscillating kernel

(1) $$\begin{array}{r}{I}^{\omega }(f;t):={\int }_{-1}^1\frac{f\left(x\right)}{x-t}{e}^{i\omega x}dx,\omega >0,-1<t<1,\end{array}$$

and compare this method with some other numerical approximations available in previous literature.

The present paper is devoted to the issue of construction of two new algorithms for evaluating the finite-part hypersingular integral

(2) $$\begin{array}{r}{J}^{\omega }(f;t):={\int }_{-1}^1\frac{f\left(x\right)}{{(x-t)}^2}{e}^{i\omega x}dx,\omega >0,-1<t<1,\end{array}$$

taking into account the results in [1].

The mathematical modeling of wave processes, electromagnetic scattering and fracture mechanics in many areas of physics and technology bring importance into the evaluation of singular and hypersingular integrals with rapidly oscillating kernel (cf. [2]–[5]), and in the last years many papers are devoted to numerical methods for approximating the integrals in (2) (see for example [6]–[8] and the references given therein). Here, we follow the idea presented first in [1] and in [9], where quadrature formulas are considered based on interpolation processes that are convergent, stable and can be implemented with small attempt.

Defined the integral (1) as Cauchy principal value

$$\begin{array}{r}{I}^{\mathrm{\omega }}(f;t)=\underset{\mathrm{\epsilon }\to 0^{+}}{lim}({\int }_{-1}^{t-\mathrm{\epsilon }}+{\int }_{t+\mathrm{\epsilon }}^1)\frac{f\left(x\right)}{x-t}{e}^{i\omega x}dx,\mathrm{\omega }>0,-1<t<1,\end{array}$$

the integral (2) can be written as the derivative of (1)

$$\begin{array}{r}{J}^{\mathrm{\omega }}(f;t)=\frac{d}{dt}{I}^{\mathrm{\omega }}(f;t),\mathrm{\omega }>0,-1<t<1,\end{array}$$

see ([10]). Moreover, from Lemma 6.1, Ch. II in [11], we can write

$$\begin{array}{r}{J}^{\mathrm{\omega }}(f;t)={\int }_{-1}^1\frac{\frac{d}{dx}\left(f\left(x\right)e^{i\mathrm{\omega }x}\right)}{x-t}dx-\frac{e^{i\mathrm{\omega }}f\left(1\right)}{1-t}-\frac{e^{-i\mathrm{\omega }}f(-1)}{1+t}=\end{array}$$

(3) $$\begin{array}{r}={\int }_{-1}^1\frac{f^{\prime }\left(x\right)e^{i\mathrm{\omega }x}}{x-t}dx+i\mathrm{\omega }{\int }_{-1}^1\frac{f\left(x\right)e^{i\mathrm{\omega }x}}{x-t}dx-\frac{e^{i\mathrm{\omega }}f\left(1\right)}{1-t}-\frac{e^{-i\mathrm{\omega }}f(-1)}{1+t},\mathrm{\omega }>0,-1<t<1.\end{array}$$

Recalling that $e^{i\mathrm{\omega }x}=\cos \mathrm{\omega }x+isin\mathrm{\omega }x$, finally we obtain

$$\begin{array}{r}{J}^{\mathrm{\omega }}(f;t)={\int }_{-1}^1\frac{\cos \mathrm{\omega }x(f^{\prime }\left(x\right)+i\mathrm{\omega }f\left(x\right))}{x-t}dx+\end{array}$$

(4) $$\begin{array}{r}{\int }_{-1}^1\frac{\sin \mathrm{\omega }x(f^{\prime }\left(x\right)+i\mathrm{\omega }f\left(x\right))}{x-t}dx-\frac{e^{i\mathrm{\omega }}f\left(1\right)}{1-t}-\frac{e^{-i\mathrm{\omega }}f(-1)}{1+t},\mathrm{\omega }>0,-1<t<1.\end{array}$$

At first, we remember that if the functions $f$ and $f^{\prime }$ are Hölder continuous then, we obtain the existence of the integrals (1)–(3) (see [10]).

Very recently, in [9] the authors have presented a method for evaluating integral in (3) making use of the values of f and f’ at the first kind Chebyshev zeros, proving bounds of the error and of the amplification factor. Although the deduced formula has the drawback of using twice as many functional evaluations, it has the advantage of having recourse only to the weights of the quadrature sum proposed in [1] to approximate (1).

In the present paper we give two new algorithms to compute (3) which although use new weights, have the advantage of making use only of the values of the function $f$ at the same points, and are convergent in a weaker condition on the smoothness of the density function $f$. Therefore, these latter formulas are useful, for instance, in quadrature method for solving hypersingular integral equations with rapidly oscillating kernel.

The paper is organized as follows: In Section 2 as well as presenting the algorithm, bounds of the error and of the amplification factor are provided; in Section 3 we present an alternative approximative method having the advantages of that ones in Section 2 about the computation and having a better convergence behavior. Finally, in Section 4 we present some numerical examples that show the coherence of the theoretical results with the numerical ones.

An Algorithm to Evaluate Integral (3)

From (4) we observe that the numerical approximation of (3) is strictly related to the quadrature of the following integrals:

(5) $$\begin{array}{r}{I}_S^{\mathrm{\omega }}(f^{\prime };t):={\int }_{-1}^1\frac{f^{\prime }\left(x\right)}{x-t}\sin \mathrm{\omega }xdx,I_C^{\mathrm{\omega }}(f^{\prime };t):={\int }_{-1}^1\frac{f^{\prime }\left(x\right)}{x-t}\cos \mathrm{\omega }xdx,\end{array}$$

and $I_S^{\mathrm{\omega }}(f;t),I_C^{\mathrm{\omega }}(f;t)$ that can be approximated by using the methods in [1]. In [9] the authors have suggested the use of

$$\begin{array}{r}{\hat{I}}_{S,m}^{\mathrm{\omega }}(v^{\mathrm{\alpha },\mathrm{\beta }};f^{\prime };t):={\int }_{-1}^1\frac{{ℒ}_m(v^{\mathrm{\alpha },\mathrm{\beta }};f^{\prime };x)}{x-t}\sin \mathrm{\omega }xdx\end{array}$$

to approximate $I_S^{\mathrm{\omega }}(f^{\prime };t)$.

Here, in order to establish a formula that does not depend on the samples of the derivative function $f^{\prime }$, we approximate $f^{\prime }$ in (5) with the Lagrange interpolant polynomial ${ℒ}_m^{\prime }$ of the function f instead of the Lagrange interpolant polynomial ${ℒ}_m$ of the function $f^{\prime }$. In what follow we will consider quadrature formulas of $I_S^{\mathrm{\omega }}(f^{\prime };t)$, since the integral $I_C^{\mathrm{\omega }}(f^{\prime };t)$ can be treated similarly. In particular, we consider the following approximation

(6) $$\begin{array}{r}{I}_{S,m}^{\mathrm{\omega }}(v^{\mathrm{\alpha },\mathrm{\beta }};f^{\prime };t):={\int }_{-1}^1\frac{{ℒ}_m^{\prime }(v^{\mathrm{\alpha },\mathrm{\beta }};f;x)}{x-t}\sin \mathrm{\omega }xdx,\end{array}$$

where the Lagrange polynomial ${ℒ}_m(v^{\mathrm{\alpha },\mathrm{\beta }};g)$ interpolates a given function $g$ at the points $x_{m,k}^{\mathrm{\alpha },\mathrm{\beta }},k=1,\dots ,m$ zeros of the $m$th Jacobi polynomial $p_m^{\mathrm{\alpha },\mathrm{\beta }},m\in N$ with respect to the exponent $\mathrm{\alpha },\mathrm{\beta }\succ -1$.

First, let us introduce some notations.

Let us denote by ${\mathrm{\omega }}_{\mathrm{\phi }}(f;\mathrm{\delta })$ the modulus of smoothness of a given function $g$, defined as

$$\begin{array}{r}{\mathrm{\omega }}_{\mathrm{\phi }}:={\text{Sup}}_{h\le \mathrm{\delta }}\underset{\left|x\right|\le 1}{max}\left|{△}_{h\mathrm{\phi }}g\left(x\right)\right|,\end{array}$$

where $\mathrm{\phi }\left(x\right)=\sqrt{1-x^2}$ and ${△}_{h\mathrm{\phi }}g\left(x\right)=g(x+h/2\mathrm{\phi }\left(x\right))-g(x-h/2\mathrm{\phi }\left(x\right))$, (cf. [12]). Further, we denote by ${\Vert g\Vert }_{\mathrm{\infty }}=\underset{\left|x\right|\le 1}{max}\left|g\left(x\right)\right|$ the usual uniform norm and by ${\mathrm{\Lambda }}_m^{\mathrm{\alpha },\mathrm{\beta }},m\in N$ the $m$th Lebesgue constant corresponding to the weight function $v^{\mathrm{\alpha },\mathrm{\beta }}\left(x\right)={(1-x)}^{\mathrm{\alpha }}{(1+x)}^{\mathrm{\beta }},\mathrm{\alpha },\mathrm{\beta }\succ 1,\left|x\right|\le 1$.

We shall study the convergence of the sequence $I_{S,m}^{\omega }(v^{-1/2,-1/2};f^{\prime };t)=I_{S,m}^{\omega }(f^{\prime };t)$ in (6) with $\mathrm{\alpha }=\mathrm{\beta }=-1/2$ to $I_S^{\mathrm{\omega }}(f^{\prime };t)$. For this purpose, we state a theorem showing the behavior of the function $I_S^{\mathrm{\omega }}(f^{\prime };t)$.

Theorem 1 Let $f\in C^1$ and $\omega \ge 0$. Then for $\left|t\right|<1$,

$$\begin{array}{r}\left|I_S^{\mathrm{\omega }}(f^{\prime };t)\right|\le c\log \frac{e}{1-t^2}\{(1+\mathrm{\omega })|f^{\prime }{|}_{\mathrm{\infty }}+{\int }_0^1\frac{{\mathrm{\omega }}_{\mathrm{\phi }}(f^{\prime };\mathrm{\delta })}{\mathrm{\delta }}d\mathrm{\delta }\},\end{array}$$

where $c$ denotes a positive constant independent of $t,f$ and $\omega$. □

Proof. See Theorem 3.1 in [1].

We recall the next result on the simultaneous polynomial approximation of a given function $g$ (see [13], Theorem p. 113).

Lemma 2 For every function $g\in C^k$, there exists an algebraic polynomial $q_m$ of degree $m\ge 4k+5$ such that

$$\begin{array}{r}|g^{\left(i\right)}\left(x\right)-q_m^{\left(i\right)}\left(x\right)|\le c{\left(\frac{\sqrt{1-x^2}}{m}\right)}^{k-i}\mathrm{\omega }(f^{\left(k\right)};\frac{\sqrt{1-x^2}}{m}),i=0,1,\dots ,k,\end{array}$$

where $\left|x\right|\le 1$ and $c$ is a positive constant independent of $m,g$ and $x$.

In the previous lemma and in what follows $\mathrm{\omega }(g;.)$ denotes the ordinary modulus of continuity of a given function $g$.

Then, for the quadrature rule (6) with $\mathrm{\alpha }=\mathrm{\beta }=-1/2$, the next result holds true.

Theorem 3 For every function $f\in C^k$, $k\ge 0$ and $\omega \ge 0$, we have

(7) $$\begin{array}{r}\left|I_{S,m}^{\mathrm{\omega }}(f^{\prime };t)\right|\le c\log \frac{e}{1-t^2}{m}^2\log m(2+\mathrm{\omega }+\log m){\Vert f\Vert }_{\mathrm{\infty }},\end{array}$$

and

(8) $$\begin{array}{r}|I_S^{\mathrm{\omega }}(f^{\prime };t)-I_{S,m}^{\mathrm{\omega }}(f^{\prime };t)|\le c\log \frac{e}{1-t^2}\{\frac{\mathrm{\omega }(f^{\left(k\right)};\frac{1}{m})}{m^{k-2}}\log m(2+\mathrm{\omega }+\log m)+{\int }_0^{\frac{1}{m}}\frac{{\omega }_{\phi }(f^{\prime };\delta )}{\delta }d\delta \},\end{array}$$

where $c$ denotes a positive constant independent of $m,f,\omega$ and $t\in (-1,1)$.

Proof. In view of Theorem 1 we can write

$$\begin{array}{rcl}& & \left|I_{S,m}^{\mathrm{\omega }}(f^{\prime };t)\right|\\ & & \le c\log \frac{e}{1-t^2}\{(1+\omega ){\Vert {ℒ}_m^{\prime }\left(f\right)\Vert }_{\mathrm{\infty }}+{\int }_0^{1/m}\frac{{\omega }_{\phi }({ℒ}_m^{\prime };\delta )}{\delta }d\delta +{\int }_{1/m}^1\frac{{\omega }_{\phi }({ℒ}_m^{\prime };\delta )}{\delta }d\delta \}\\ & & \le c\log \frac{e}{1-t^2}\{(1+\mathrm{\omega }+\log m){\Vert {ℒ}_m^{\prime }\left(f\right)\Vert }_{\mathrm{\infty }}+\frac{1}{m}{\Vert {ℒ}_m^{″}\left(f\right)\Vert }_{\mathrm{\infty }}\}\\ & & \le c\log \frac{e}{1-t^2}(2+\mathrm{\omega }+\log m){\Vert {ℒ}_m^{\prime }\left(f\right)\Vert }_{\mathrm{\infty }}\\ & & \le c\log \frac{e}{1-t^2}{m}^2(2+\mathrm{\omega }+\log m){\Vert {ℒ}_m\left(f\right)\Vert }_{\mathrm{\infty }},\end{array}$$

by using Bernstein inequality. Hence (7) follows from ${\Vert {ℒ}_m\left(f\right)\Vert }_{\mathrm{\infty }}=\mathrm{????}(\log m)$ when $\mathrm{\alpha }=\mathrm{\beta }=-1/2$.

In order to prove (8) we remark that in view of Theorem 1

$$\begin{array}{r}|I_S^{\mathrm{\omega }}(f^{\prime };t)-I_{S,m}^{\mathrm{\omega }}(f^{\prime };t)|=\left|I_S^{\mathrm{\omega }}(f^{\prime }-{ℒ}_m^{\prime }\left(f\right);t)\right|\end{array}$$

(9) $$\begin{array}{r}\le c\log \frac{e}{1-t^2}\{(1+\omega ){\Vert f^{\prime }-{ℒ}_m^{\prime }\left(f\right)\Vert }_{\mathrm{\infty }}+{\int }_0^1\frac{{\omega }_{\phi }(f^{\prime }-{ℒ}_m^{\prime };\delta )}{\delta }d\delta \}.\end{array}$$

Now let $q_{m-1}$ be the polynomial of Lemma 2. Thus

(10) $$\begin{array}{rcl}& & {\Vert f^{\prime }-{ℒ}_m^{\prime }\left(f\right)\Vert }_{\mathrm{\infty }}\le {\Vert f^{\prime }-q_{m-1}^{\prime }\Vert }_{\mathrm{\infty }}+{\Vert {ℒ}_m^{\prime }(q_{m-1}-f)\Vert }_{\mathrm{\infty }}\\ & & \le c\{\frac{1}{m^{k-1}}\mathrm{\omega }(f^{(k-1)};\frac{1}{m})+\frac{1}{m^k}\mathrm{\omega }(f^{\left(k\right)};\frac{1}{m}){\Vert {ℒ}_m^{\prime }\Vert }_{\mathrm{\infty }}\}\\ & & \le c\frac{1}{m^{k-2}}\mathrm{\omega }(f^{\left(k\right)};\frac{1}{m})logm.\end{array}$$

Further

$$\begin{array}{l}{\displaystyle {\int }_0^1\frac{{\omega }_{\phi }(f^{\prime }-{ℒ}_m^{\prime };\delta )}{\delta }d\delta }\\ \le {\int }_0^{\frac{1}{m}}\frac{{\omega }_{\phi }(f^{\prime }-q_{m-1}^{\prime };\delta )}{\delta }d\delta +{\int }_{\frac{1}{m}}^1\frac{{\omega }_{\phi }(f^{\prime }-q_{m-1}^{\prime };\delta )}{\delta }d\delta +{\int }_0^{\frac{1}{m}}\frac{{\omega }_{\phi }({ℒ}_m^{\prime }(q_{m-1}-f);\delta )}{\delta }d\delta +{\int }_{\frac{1}{m}}^1\frac{{\omega }_{\phi }({ℒ}_m^{\prime }(q_{m-1}-f);\delta )}{\delta }d\delta \end{array}$$

(11) $$\begin{array}{rcl}& & \le c\{{\int }_0^{\frac{1}{m}}\frac{{\omega }_{\phi }(f^{\prime };\delta )}{\delta }d\delta +logm{\Vert f^{\prime }-q_{m-1}^{\prime }\Vert }_{\mathrm{\infty }}+\frac{1}{m}{\Vert {ℒ}_m^{″}(q_{m-1}-f)\phi \Vert }_{\mathrm{\infty }}\\ & & +logm{\Vert {ℒ}_m^{\prime }(q_{m-1}-f)\Vert }_{\mathrm{\infty }}\}\\ & & \le c\{{\int }_0^{\frac{1}{m}}\frac{{\omega }_{\phi }(f^{\prime };\delta )}{\delta }d\delta +\frac{logm}{m^{k-1}}\mathrm{\omega }(f^{\left(k\right)};\frac{1}{m})\\ & & +m^2(1+\log m){\Vert {ℒ}_m\Vert }_{\mathrm{\infty }}{\Vert q_{m-1}-f\Vert }_{\mathrm{\infty }}\}\end{array}$$

$$\begin{array}{r}\le c\{{\int }_0^{\frac{1}{m}}\frac{{\omega }_{\phi }(f^{\prime };\delta )}{\delta }d\delta +\frac{logm}{m^{k-2}}(1+\log m)\mathrm{\omega }(f^{\left(k\right)};\frac{1}{m})\},\end{array}$$

having used again Bernstein inequality, Lemma 2 and ${\Vert {ℒ}_m\left(f\right)\Vert }_{\mathrm{\infty }}=\mathrm{????}(\log m)$ for $\mathrm{\alpha }=\mathrm{\beta }=-1/2$. Combining (10) and (11) with (9) we deduce (8). □

We want to emphasize that from (7) we can deduce the following bound

$$\begin{array}{r}\Vert I_{S,m}^{\omega }\left(f^{\prime }\right)lo{g}^{-1}\frac{e}{1-{(.)}^2}\Vert \le c{m}^{2}logm(2+\omega +logm),\end{array}$$

that provides the behavior of the weighted amplification factor.

We go on to see how to compute $I_{S,m}^{\omega }(f^{\prime };t)$ in (6) with $\mathrm{\alpha }=\mathrm{\beta }=-1/2$. So we denote by

$$\begin{array}{r}{p}_0^{-\frac{1}{2},-\frac{1}{2}}\left(x\right)=p_0\left(x\right)=\frac{1}{\sqrt{\pi }},p_m^{-\frac{1}{2},-\frac{1}{2}}\left(x\right)=p_m\left(x\right)=\sqrt{\frac{2}{\pi }}{T}_m\left(x\right),m\ge 1\end{array}$$

the m-th Chebyshev orthonormal polynomial of the first kind and let $x_{m,k}^{-1/2,-1/2}=x_{m,k}=cos\left((2k-1)\pi /2m\right)$ be the zeros of the orthogonal polynomial $T_m$. Since

$$\begin{array}{r}{ℒ}_m^{\prime }(v^{-\frac{1}{2},-\frac{1}{2}};f;x)={ℒ}_m^{\prime }(f;x)=\sum _{k=1}^m{l}_{m,k}^{\prime }\left(x\right)f\left(x_{m,k}\right),\end{array}$$

where $l_{m,k},k=1,\dots ,m$ are the fundamental Lagrange polynomials with respect to the points $x_{m,k},k=1,\dots ,m$, we have

$$\begin{array}{r}{I}_{S,m}^{\mathrm{\omega }}(f^{\prime };t)=\sum _{k=1}^m[{\int }_{-1}^1\frac{l_{m,k}^{\prime }\left(x\right)}{x-t}\sin \omega xdx]f\left(x_{m,k}\right),\end{array}$$

with

$$\begin{array}{r}{l}_{m,k}^{\prime }\left(x\right)=\sum _{i=1}^{m-1}{a}_i\frac{d}{dx}{p}_i\left(x\right),k=1,\dots ,m,\end{array}$$

and

$$\begin{array}{r}{a}_i={\int }_{-1}^1\frac{l_{m,k}\left(x\right)p_i\left(x\right)}{\sqrt{1-x^2}}dx=\frac{\pi }{m}\sum _{j=1}^m{l}_{m,k}\left(x_{m,j}\right)p_i\left(x_{m,j}\right)=\frac{\pi }{m}{p}_i\left(x_{m,k}\right),i=1,\dots ,m-1.\end{array}$$

Thus,

$$\begin{array}{r}{l}_{m,k}^{\prime }\left(x\right)=\frac{2}{m}\sum _{i=1}^{m-1}{T}_i\left(x_{m,k}\right)T_i^{\prime }\left(x\right),k=1,\dots ,m,\end{array}$$

and

(12) $$\begin{array}{r}{I}_{S,m}^{\mathrm{\omega }}(f^{\prime };t)=\frac{2}{m}\sum _{k=1}^m\left[\sum _{i=1}^{m-1}{T}_i\left(x_{m,k}\right)q_i^{\prime \omega }\left(t\right)\right]f(x_{m,k}),\end{array}$$

where

$$\begin{array}{r}{q}_i^{\prime \omega }\left(t\right)={\int }_{-1}^1\frac{T_i^{\prime }\left(x\right)}{x-t}\sin \omega xdx,i=1,2,\dots .\end{array}$$

Denoting by ${\left\{U_n\right\}}_{n\in \mathbb{N}}$ the sequence of the Chebyshev orthogonal polynomials of the second kind, we have

$$\begin{array}{r}{T}_n^{\prime }\left(x\right)=n{U}_{n-1}\left(x\right),n=1,2,\dots ,\end{array}$$

and

$$\begin{array}{r}{U}_0\left(x\right)\equiv 1,U_1\left(x\right)=2x,U_{n+1}\left(x\right)=2x{U}_n\left(x\right)-U_{n-1}\left(x\right),n=1,2,\dots .\end{array}$$

Therefore

$$\begin{array}{r}{q}_i^{\prime \omega }\left(t\right)=i{\overline{q}}_{i-1}^{\omega }\left(t\right),i=1,2,\dots ,\end{array}$$

where

$$\begin{array}{r}{\overline{q}}_{i-1}^{\omega }\left(t\right)={\int }_{-1}^1\frac{U_i\left(x\right)}{x-t}\sin \omega xdx,i=0,1,\dots .,\end{array}$$

(13) $$\begin{array}{r}{\overline{q}}_{n+1}^{\omega }\left(t\right)=2t{\overline{q}}_n^{\omega }\left(t\right)-{\overline{q}}_{n-1}^{\omega }\left(t\right)+2{\overline{M}}_n^{\omega },n=1,2,\dots ,\end{array}$$

and

$$\begin{array}{r}{\overline{M}}_n^{\omega }={\int }_{-1}^1{U}_n\left(x\right)\sin \omega xdx,n=0,1,\dots .\end{array}$$

The accurate evaluation of ${\overline{M}}_n^{\omega }$ in (13) allows us to compute ${\overline{q}}_n^{\omega }\left(t\right)$ for $n=0,1,\dots ,$ together with

$$\begin{array}{r}{\overline{q}}_0^{\omega }\left(t\right)={\int }_{-1}^1\frac{\sin \omega x}{x-t}dx=\sin \omega x[Ci\left({\tau }_1\right)-Ci\left(\left|{\tau }_2\right|\right)]+\cos \omega x[Si\left({\tau }_1\right)+Si\left(\left|{\tau }_2\right|\right)],\end{array}$$

and

$$\begin{array}{r}{\overline{q}}_1^{\omega }\left(t\right)=2{\int }_{-1}^1\frac{x\sin \omega x}{x-t}dx=2t{\overline{q}}_0^{\omega }\left(t\right),\end{array}$$

where

$$\begin{array}{r}Si\left(\tau \right)={\int }_0^{\tau }\frac{\sin x}{x}dx,Ci\left(\tau \right)={\int }_0^{\tau }\frac{\cos x-1}{x}dx+log\tau +C,\tau >0,\end{array}$$

are the sine and cosine integral, respectively; ${\tau }_1=\omega (1-t),{\tau }_2=-\omega (1+t)$ and $C$ is the Euler constant. The starting values of (13) require the evaluation of the sine and cosine integrals that can be computed by some mathematical software like Mathematica [14]. Finally, we remark that (12) can be rewritten

(14) $$\begin{array}{r}{I}_{S,m}^{\mathrm{\omega }}(f^{\prime };t)=\frac{2}{m}\sum _{i=0}^{m-2}{c}_{m,i}(f){\overline{q}}_i^{\omega }\left(t\right),\end{array}$$

with respect to the coefficients

$$\begin{array}{r}{c}_{m,i}\left(f\right)=\sum _{k=1}^m(i+1)T_{i+1}\left(x_{m,k}\right)f\left(x_{m,k}\right),i=0,2,\dots ,m-2,\end{array}$$

which are not influenced by the value t and the oscillatory factor ω. Thus the evaluation of $I_{S,m}^{\mathrm{\omega }}(f^{\prime };t)$ in (14) can be done following Clenshaw type algorithm of this kind:

$$\begin{array}{r}{z}_m=z_{m-1}=0,w_{m-1}=0,\end{array}$$

$$\begin{array}{r}{z}_k=2t{z}_{k+1}-z_{k+2}+c_{m,k}\left(f\right),k=m-2,m-3,\dots ,0\end{array}$$

$$\begin{array}{r}{w}_k=2z_{k+1}{\overline{M}}_k^{\omega }+w_{k+1},k=m-2,m-3,\dots ,0\end{array}$$

$$\begin{array}{r}{I}_{S,m}^{\mathrm{\omega }}(f^{\prime };t)=\frac{2}{m}({\overline{q}}_0^{\omega }\left(t\right)z_0+w_0).\end{array}$$

We want point out that even if the quadrature $I_{S,m}^{\mathrm{\omega }}(f^{\prime };t)$ is to preferred to use formula (6) in [9], because it does not use the values $f^{\prime }\left(x_{m,k}\right),k=1,\dots ,m$, from a convergence point of view it performs worse (cf. the previous theorem and Theorem 1 in [9]). This is due to the fact that $I_{S,m}^{\mathrm{\omega }}(f^{\prime };t)$ uses the operator ${ℒ}_m^{\prime }$ which performs worse that the operator ${ℒ}_m$ used in [9]. Indeed, the Lagrange operator is not good for the simultaneous approximation even if we start from an optimal choice of interpolation points as in the case $\mathrm{\alpha }=\mathrm{\beta }=-1/2$.

Another Algorithm to Evaluate Integral (3)

In the sequel we shall propose a new formula to compute $I_S^{\mathrm{\omega }}(f^{\prime };t)$ having the advantages of $I_{S,m}^{\mathrm{\omega }}(f^{\prime };t)$ about the computation and having the same convergence behavior of formula (2.6) in [1].

We introduce the polynomial ${ℒ}_{m,1,1}\left(g\right)$ interpolating a given function $g$ at the points $x_{m,k},k=1,\dots ,m$, zeros of $T_m$ and at $x_{m,0}=-1,x_{m,m+1}=1$

${ℒ}_{m,1,1}(g;x)=(1-x^2){ℒ}_m({\displaystyle \frac{g}{1-{(.)}^2}};x)+[{(-1)}^m(1-x)g(-1)+(1+x)g\left(1\right)]{\displaystyle \frac{T_m(x)}{2}}$ .

Then we consider the new quadrature rule

$$\begin{array}{r}{I}_{S,m,1,1}^{\omega }(f^{\prime };t)=I_S^{\omega }({ℒ}_{m,1,1}^{\prime }\left(f\right);t)\end{array}$$

to approximate $I_S^{\mathrm{\omega }}(f^{\prime };t)$. For it we can prove the following result.

Theorem 4. For every function $f\in C^k,k\ge 0$ and $\omega \ge 0$, we have

(15) $$\begin{array}{r}\left|I_{S,m,1,1}^{\mathrm{\omega }}(f^{\prime };t)\right|\le clog\frac{e}{1-t^2}logm(2+\omega +logm){\Vert f\Vert }_{\mathrm{\infty }},\end{array}$$

and

(16) $$\begin{array}{r}|I_S^{\mathrm{\omega }}(f^{\prime };t)-I_{S,m,1,1}^{\mathrm{\omega }}(f^{\prime };t)|\le clog\frac{e}{1-t^2}\{\frac{\omega (f^{\left(k\right)};\frac{1}{m})}{m^{k-1}}logm(2+\omega +logm)+{\int }_0^{\frac{1}{m}}\frac{{\omega }_{\phi }(f^{\prime };\delta )}{\delta }d\delta \},\end{array}$$

where $c$ denotes a positive constant independent of $m,f,\omega$ and $t\in (-1,1)$.

Proof. To prove (15) and (16) we can follow the same steps to prove (7) and (8), respectively. Thus, the proof follows recalling that $\Vert {ℒ}_{m,1,1}^{\prime }\Vert =\mathrm{????}\left(logm\right)$ in view of Corollary 3.2 in [15].

We remark that the assumption of the existence of the Hölder continuous derivative $f^{\prime }$ besides ensuring the existence of$I_S^{\mathrm{\omega }}(f^{\prime };t)$ it provides the convergence of $I_{S,m,1,1}^{\mathrm{\omega }}(f^{\prime };t)$ (cf. Theorem 4). Instead, the same hypothesis on the smoothness on $f$ does not ensure the convergence of $I_{S,m}^{\mathrm{\omega }}(f;t)$ (cf. Theorem 3).

Finally, by standard computations we try

$$\begin{array}{rcl}{I}_{S,m,1,1}^{\mathrm{\omega }}(f^{\prime };t)& =& \frac{2}{m}\sum _{k=1}^m\{\sum _{j=2}^{m-1}{T}_j\left(x_{m,k}\right)[-2(t{q}_j^{\omega }\left(t\right)+M_j^{\omega })+(1-t^2)j{\overline{q}}_{j-1}^{\omega }\left(t\right)-jt{\overline{M}}_{j-1}^{\omega }\\ & & -j({\overline{M}}_j^{\omega }-M_j^{\omega })]-t{q}_0^{\omega }\left(t\right)+x_{m,k}[\frac{6}{{\omega }^2}[\omega \cos \omega -\sin \omega ]-2t{q}_1^{\omega }\left(t\right)\\ & & +(1-t^2){\overline{q}}_0^{\omega }\left(t\right)]\}\frac{f\left(x_{m,k}\right)}{1-x_{m,k}^2}+\frac{{(-1)}^{m+1}}{2}\{q_m^{\omega }\left(t\right)+m[{\overline{M}}_{m-1}^{\omega }-(1-t){\overline{q}}_{m-1}^{\omega }\left(t\right)]\}\end{array}$$

$$\begin{array}{r}f(-1)+\frac{1}{2}\{q_m^{\omega }\left(t\right)+m[{\overline{M}}_{m-1}^{\omega }+(1+t){\overline{q}}_{m-1}^{\omega }\left(t\right)]\}f\left(1\right),\end{array}$$

where ${\overline{q}}_n^{\omega }\left(t\right)$ and ${\overline{M}}_n^{\omega },n=0,1,2\dots$ are the same defined before and where

$$\begin{array}{r}{q}_n^{\omega }\left(t\right)={\int }_{-1}^1\frac{T_n\left(x\right)}{x-t}\sin \omega xdx,n=0,1,\dots ,\end{array}$$

$$\begin{array}{r}{M}_n^{\omega }={\int }_{-1}^1{T}_n\left(x\right)\sin \omega xdx,n=0,1,\dots .\end{array}$$

Recalling that the polynomials $T_n,n\in \mathbb{N}$, satisfy

$$\begin{array}{r}{T}_0\left(x\right)\equiv 1,T_1\left(x\right)=x,T_{n+1}\left(x\right)=2x{T}_n\left(x\right)-T_{n-1}\left(x\right),n=1,2,\dots ,\end{array}$$

we try

(17) $$\begin{array}{r}{q}_{n+1}^{\omega }(t)=2t{q}_n^{\omega }(t)-q_{n-1}^{\omega }(t)+2M_n^{\omega },n=1,2,\dots .\end{array}$$

The evaluation of the integral $M_n^{\omega }$ in (17) allows us to compute $q_n^{\omega }(t)$ for $n=1,2,\dots$, together with

$q_0^{\omega }\left(t\right)={\overline{q}}_0^{\omega }\left(t\right)$ and $q_1^{\omega }\left(t\right)=1/2{\overline{q}}_1^{\omega }\left(t\right)$.

Numerical Examples

In this Section we consider two numerical examples with the aim of showing the correspondence between the numerical results and the theoretical ones by applying the algorithms presented in the previous Sections. All the computations have been performed in double precision arithmetic.

In the first test we choose the function $f\left(x\right)=\exp \left(x\right)$, so that $f^{\prime }(x)=exp(x)$. In this case we know the exact solution, therefore we compare this with the numerical solutions obtained using the two proposed methods. We denote by

$E_{S,m}^{\mathrm{\omega }}(f^{\prime };t)=I_S^{\mathrm{\omega }}(f^{\prime };t)-I_{S,m}^{\mathrm{\omega }}(f^{\prime };t)$ and $E_{S,m,1,1}^{\mathrm{\omega }}(f^{\prime };t)=I_S^{\mathrm{\omega }}(f^{\prime };t)-I_{S,m,1,1}^{\mathrm{\omega }}(f^{\prime };t),$

the errors obtained with two different methods. In Tables I–III we show the errors of the two methods, compared with the exact solution of the integral evaluated with Mathematica package, for three different values of $t\in (-1,1),$ $\mathrm{\omega }$ and for increasing values of $m\in \mathbb{N}$.

Table I. $f\left(x\right)=exp\left(x\right),\mathrm{\omega }=10$

m	${\text{E}}_{\text{S},\text{m}}^{\mathrm{\omega }}({\text{f}}^{\prime };0.1)$	${\text{E}}_{\text{S},\text{m},1,1}^{\mathrm{\omega }}({\text{f}}^{\prime };0.1)$	${\text{E}}_{\text{S},\text{m}}^{\mathrm{\omega }}({\text{f}}^{\prime };0.5)$	${\text{E}}_{\text{S},\text{m},1,1}^{\mathrm{\omega }}({\text{f}}^{\prime };0.5)$	${\text{E}}_{\text{S},\text{m}}^{\mathrm{\omega }}({\text{f}}^{\prime };0.9)$	${\text{E}}_{\text{S},\text{m},1,1}^{\mathrm{\omega }}({\text{f}}^{\prime };0.9)$
4	0.1D-01	0.4D-03	0.1D-01	0.3D-03	0.9D-01	0.1D-02
8	0.2D-05	0.2D-07	0.5D-06	0.1D-07	0.6D-05	0.5D-07
$\ge 16$	0.1D-14	0.5D-15	0.2D-14	0.3D-15	0.1D-10	0.9D-14

Table II. $f\left(x\right)=exp\left(x\right),\mathrm{\omega }=50$

m	${\text{E}}_{\text{S},\text{m}}^{\mathrm{\omega }}({\text{f}}^{\prime };0.1)$	${\text{E}}_{\text{S},\text{m},1,1}^{\mathrm{\omega }}({\text{f}}^{\prime };0.1)$	${\text{E}}_{\text{S},\text{m}}^{\mathrm{\omega }}({\text{f}}^{\prime };0.5)$	${\text{E}}_{\text{S},\text{m},1,1}^{\mathrm{\omega }}({\text{f}}^{\prime };0.5)$	${\text{E}}_{\text{S},\text{m}}^{\mathrm{\omega }}({\text{f}}^{\prime };0.9)$	${\text{E}}_{\text{S},\text{m},1,1}^{\mathrm{\omega }}({\text{f}}^{\prime };0.9)$
4	0.8D-02	0.2D-03	0.8D-01	0.2D-02	0.8D-01	0.8D-03
8	0.1D-05	0.1D-07	0.5D-05	0.5D-07	0.6D-05	0.9D-09
$\ge 16$	0.1D-14	0.3D-15	0.2D-14	0.2D-15	0.3D-14	0.8D-15

Table III. $f\left(x\right)=exp\left(x\right),\mathrm{\omega }=100$

m	${\text{E}}_{\text{S},\text{m}}^{\mathrm{\omega }}({\text{f}}^{\prime };0.1)$	${\text{E}}_{\text{S},\text{m},1,1}^{\mathrm{\omega }}({\text{f}}^{\prime };0.1)$	${\text{E}}_{\text{S},\text{m}}^{\mathrm{\omega }}({\text{f}}^{\prime };0.5)$	${\text{E}}_{\text{S},\text{m},1,1}^{\mathrm{\omega }}({\text{f}}^{\prime };0.5)$	${\text{E}}_{\text{S},\text{m}}^{\mathrm{\omega }}({\text{f}}^{\prime };0.9)$	${\text{E}}_{\text{S},\text{m},1,1}^{\mathrm{\omega }}({\text{f}}^{\prime };0.9)$
4	0.2D-01	0.7D-03	0.8D-01	0.2D-02	0.9D-01	0.6D-03
8	0.3D-05	0.3D-07	0.5D-05	0.5D-07	0.1D-05	0.2D-08
$\ge 16$	0.9D-15	0.5D-15	0.8D-14	0.4D-15	0.1D-13	0.4D-15

In the second example we choose the function $f\left(x\right)=1/2(x\sqrt{1-x^2}+\arcsin x)$, so that $f^{\prime }\left(x\right)=\sqrt{1-x^2}$. In such case we don’t know the exact solution, and in Tables IV–VI we consider only the correct digits obtained using the two proposed methods and, as in the first example, for three different values of $t\in (-1,1)$, $\mathrm{\omega }$ and for increasing values of $m\in \mathbb{N}$. Taking into account the regularity of the functions considered in the examples, Theorems (3) and (4) give the same order of convergence, as we can see looking at the results presented in Tables I–III. Better results are obtained with the second method when the function is less regular.

Table IV. $f\left(x\right)=1/2(x\sqrt{1-x^2}+\arcsin x),\mathrm{\omega }=10$

m	$I_{\text{S},\text{m}}^{\mathrm{\omega }}({\text{f}}^{\prime };0.1)$	$I_{\text{S},\text{m},1,1}^{\mathrm{\omega }}({\text{f}}^{\prime };0.1)$	$I_{\text{S},\text{m}}^{\mathrm{\omega }}({\text{f}}^{\prime };0.5)$	$I_{\text{S},\text{m},1,1}^{\mathrm{\omega }}({\text{f}}^{\prime };0.5)$	$I_{\text{S},\text{m}}^{\mathrm{\omega }}({\text{f}}^{\prime };0.9)$	$I_{\text{S},\text{m},1,1}^{\mathrm{\omega }}({\text{f}}^{\prime };0.9)$
4	1.7	1.7	0.8	0.8	-1.	−1.0
8	1.7	1.7	0.8	0.8	−1.07	−1.0
16	1.763	1.763	0.86	0.86	−1.07	−1.07
$\ge 32$	1.7634	1.7634	0.8641	0.8641	−1.073	−1.073

Table V. $f\left(x\right)=1/2(x\sqrt{1-x^2}+\arcsin x),\mathrm{\omega }=50$

m	${\text{I}}_{\text{S},\text{m}}^{\mathrm{\omega }}({\text{f}}^{\prime };0.1)$	$I_{\text{S},\text{m},1,1}^{\mathrm{\omega }}({\text{f}}^{\prime };0.1)$	$I_{\text{S},\text{m}}^{\mathrm{\omega }}({\text{f}}^{\prime };0.5)$	$I_{\text{S},\text{m},1,1}^{\mathrm{\omega }}({\text{f}}^{\prime };0.5)$	$I_{\text{S},\text{m}}^{\mathrm{\omega }}({\text{f}}^{\prime };0.9)$	$I_{\text{S},\text{m},1,1}^{\mathrm{\omega }}({\text{f}}^{\prime };0.9)$
4	0.8	0.8	2.6	2.6	0.7	0.6
8	0.8	0.883	2.69	2.6	0.7	0.7
16	0.8832	0.8832	2.692	2.692	0.7	0.70
$\ge 32$	0.88327	0.88327	2.6924	2.6924	0.7097	0.7097

Table VI. $f\left(x\right)=1/2(x\sqrt{1-x^2}+\arcsin x),\mathrm{\omega }=100$

m	${\text{I}}_{\text{S},\text{m}}^{\mathrm{\omega }}({\text{f}}^{\prime };0.1)$	$I_{\text{S},\text{m},1,1}^{\mathrm{\omega }}({\text{f}}^{\prime };0.1)$	$I_{\text{S},\text{m}}^{\mathrm{\omega }}({\text{f}}^{\prime };0.5)$	$I_{\text{S},\text{m},1,1}^{\mathrm{\omega }}({\text{f}}^{\prime };0.5)$	$I_{\text{S},\text{m}}^{\mathrm{\omega }}({\text{f}}^{\prime };0.9)$	$I_{\text{S},\text{m},1,1}^{\mathrm{\omega }}({\text{f}}^{\prime };0.9)$
4	−2.6	−2.6	2.	2.6	−0.7	−0.57
8	−2.62	−2.62	2.624	2.62	−0.61	−0.6
16	−2.623	−2.623	2.624	2.624	−0.61	−0.61
$\ge 32$	−2.6234	−2.6234	2.6246	2.6246	−0.615	−0.615