Computation of integrals with oscillatory and singular integrands using Chebyshev expansions

Journal of Computational and Applied Mathematics 242 (2013) 141–156

Contents lists available at SciVerse ScienceDirect

Journal of Computational and Applied Mathematics journal homepage: www.elsevier.com/locate/cam

Computation of integrals with oscillatory and singular integrands using Chebyshev expansions✩ Hongchao Kang, Shuhuang Xiang ∗ , Guo He School of Mathematics and Statistics, Central South University, Changsha, Hunan 410083, PR China

article

info

Article history: Received 20 March 2012 Received in revised form 19 October 2012

abstract We present a general method for computing oscillatory integrals of the form



1

f (x)G(x)eiωx dx,

−1

MSC: 65D32 65D30 Keywords: Oscillatory integrals Singularity Chebyshev expansion Modified moments Recurrence relations

where f is sufficiently smooth on [−1, 1], ω is a positive parameter and G is a product of singular factors of algebraic or logarithmic type. Based on a Chebyshev expansion of f and the properties of Chebyshev polynomials, the proposed method for such integrals is constructed with the help of the expansion of the oscillatory factor eiωx . Furthermore, due to numerically stable recurrence relations for the modified moments, the devised scheme can be employed to compute oscillatory integrals with algebraic or logarithmic singularities at the end or interior points of the interval of integration. Numerical examples are provided to confirm our analysis. Crown Copyright © 2012 Published by Elsevier B.V. All rights reserved.

1. Introduction In this paper we are concerned with the numerical evaluation of oscillatory integrals of the form



1

f (x)G(x)eiωx dx,

(1.1)

−1

where f is a sufficiently smooth function on [−1, 1], ω is a positive parameter and G is a product of singular factors of algebraic or logarithmic type. Here, weight functions G(x) and parameters are listed in Table 1. 1 The first case of Table 1 corresponds to −1 f (x)eiωx dx and the integrals of the forms 2–5, 6, 8, 10–17 belong to the special cases of the form 18. The integrals 2–18 with the weak singularities arise in the numerical approximations of solutions to Volterra integral equations of the first kind involving highly oscillatory kernels with weak singularities [1,2]. In addition, it is well-known that the Radon transform, which plays an important role in the CT, PET and SPECT technology of medical sciences and is widely applicable to tomography, the creation of an image from the scattering data associated to crosssectional scans of an object, is closely related to this form of oscillatory singular integrals [3,4]. The singularities in (1.1) are also called singularities of the Radon transform in medical tomography. Moreover, the integrals 2–18 are also used to the solution of the singular integral equation for classical crack problems in plane and antiplane elasticity [5]. Further, they can be taken as model integrals appearing in boundary integral equations for high-frequency acoustic scattering (e.g., highfrequency Helmholtz equation in two dimensions), where the kernels have algebraic or logarithmic singularities on the diagonal ([6,7] and references therein), which is also our main target application. ✩ This paper is supported by NSF of China (No. 11071260).

∗

Corresponding author. E-mail addresses: [email protected] (H. Kang), [email protected] (S. Xiang), [email protected] (G. He).

0377-0427/$ – see front matter Crown Copyright © 2012 Published by Elsevier B.V. All rights reserved. doi:10.1016/j.cam.2012.10.016

142

H. Kang et al. / Journal of Computational and Applied Mathematics 242 (2013) 141–156 Table 1 Weight functions G(x) and parameters. Number c

G(x)(G[c ] (x))

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18

1 (1 − x )α (1 + x )β ln(1 − x)(1 − x)α (1 + x)β (1 − x)α (1 + x)β ln(1 + x) ln(1 − x)(1 − x)α (1 + x)β ln(1 + x) |x − a|α |x − a|α sgn(x − a) |x − a|α ln |x − a| |x − a|α ln |x − a|sgn(x − a) (1 − x)α (1 + x)β |x − a|γ (1 − x)α (1 + x)β |x − a|γ ln |x − a| ln(1 − x)(1 − x)α (1 + x)β |x − a|γ (1 − x)α (1 + x)β ln(1 + x)|x − a|γ ln(1 − x)(1 − x)α (1 + x)β ln(1 + x)|x − a|γ (1 − x)α (1 + x)β ln(1 + x)|x − a|γ ln |x − a| ln(1 − x)(1 − x)α (1 + x)β |x − a|γ ln |x − a| ln(1 − x)(1 − x)α (1 + x)β ln(1 + x)|x − a|γ ln |x − a| (ln(1 − x)s )p (1 − x)α (1 + x)β (ln(1 + x)t )q |x − a|γ (ln |x − a|m )r

Parameters

α, β > −1 α, β > −1 α, β > −1 α, β > −1 α > −1, |a| < 1 α > −1, |a| < 1 α > −1, |a| < 1 α > −1, |a| < 1 α, β, γ > −1, |a| < 1 α, β, γ > −1, |a| < 1 α, β, γ > −1, |a| < 1 α, β, γ > −1, |a| < 1 α, β, γ > −1, |a| < 1 α, β, γ > −1, |a| < 1 α, β, γ > −1, |a| < 1 α, β, γ > −1, |a| < 1 α, β, γ > −1, |a| < 1, s, t , m ∈ R, p, q, r = 0, 1, 2, . . .

The singularities of algebraic or logarithmic type and possible high oscillations of the integrands make the integrals (1.1) very difficult to approximate accurately using standard methods, e.g., Gauss–Legendre rule. It should be also pointed out that many efficient methods, such as the Filon-type [8,9], Levin-type [10] and numerical steepest descent methods [11], cannot be applied directly to the integrals of the forms 3–18 in Table 1, since the nonoscillatory parts of the integrands are undefined at the endpoints or one point in the interior of the interval. Moreover, for singular weight functions G(x) of the forms 3–18 in Table 1, their moments may be awkward to compute. For example, for all α, β, γ > −1, |a| < 1, s, t , m ∈ R, p, q, r = 0, 1, 2, . . . , neither the moments



1

xk (ln(1 − x)s )p (1 − x)α (1 + x)β (ln(1 + x)t )q |x − a|γ (ln |x − a|m )r eiωx dx,

(1.2)

−1

nor the Chebyshev moments



1

Tk (x)(ln(1 − x)s )p (1 − x)α (1 + x)β (ln(1 + x)t )q |x − a|γ (ln |x − a|m )r eiωx dx,

(1.3)

−1

is often easily computed, where Tk (x) is the Chebyshev polynomial of the first kind of degree k. Unfortunately, since both the Filon-type method [8,9] and the Clenshaw–Curtis-type method [12] require that the above two moments (1.2)–(1.3) can be easily calculated, they are unavailable for the integrals of the forms 3–18 in Table 1. For the integrals (1.1), the method proposed in this paper first expands f into a series of Chebyshev polynomials, based on the ideas presented in [13]. Then, the factor eiωx is expanded into a series of Chebyshev polynomials and Bessel functions of the first kind. And, thanks to the properties of Chebyshev polynomials, we extend, in a sense, the ideas presented in [14], where the recurrence formula for the calculation of the modified moments



1

(1 − x)α (1 + x)β Tn (x)dx,

(1.4)

−1

was presented. It is worth noting that the modified moments (1.4) play an important role in the construction of quadrature 1 formulas for the integrals (1.1). Furthermore, all the required moments −1 Tk (x)G(x)eiωx dx, are expressed in terms of the simpler ones



1

Tk (x)G(x)dx,

(1.5)

−1

which can be computed efficiently. What is more important, the moments (1.5) satisfy recurrence relations that are stable in either forward or backward direction, which makes the whole algorithm quite simple. Hence, the presented method allows us to overcome some problems that usually appear in the case when the integrands are both singular and oscillating. Here, we would also like to mention several other papers related to the ideas given in this article. In [15,16], the algorithms 1 for computing integrals of the form −1 f (x)K (x)dx were proposed for a few simple singular/oscillating functions K . As in this paper, the algorithms mentioned in [15,16], are based on expanding the function f into a series of Chebyshev polynomials. The methods presented in [15,16] were constructed only for four simple functions K (x), but in [17,18] the idea was extended and generalized for a wide class of oscillatory or singular integrands. In [19,20], two different approaches based

H. Kang et al. / Journal of Computational and Applied Mathematics 242 (2013) 141–156

143

on Chebyshev series expansions were proposed for computing Cauchy principle value integrals of oscillatory functions. The ideas presented in [15–18,20] extend the Clenshaw–Curtis indefinite approach described in [13]. More precisely, a  function F is constructed (given in terms of Chebyshev polynomials), such that F (x) = f (x)K (x)dx. Then, one can simply

1

set −1 f (x)K (x)dx = F (1) − F (−1). It should be noted that we also partially use the similar indefinite approach in the construction of the devised method. 1 It should be also mentioned that for the integral −1 f (x)eiωx dx (the form 1 in Table 1) faster algorithms are already known, such as Clenshaw–Curtis-type method [14,12], Filon-type method [8,9], Levin method [21], Levin-type method [10], numerical steepest descent method [11] and other methods [22,15,17]. A fast algorithm for computing the integral 2 in Table 1 was also proposed in [23]. And the detailed algorithm for computing the integrals of the form 6 in Table 1 was presented in [17] as an example, where the complexity of that algorithm does not depend on ω. However, we do not remove the kernel functions of the forms 1, 2, 6 from the list given in Table 1 as it shows wide applications of the proposed method. of this paper is as follows. In Section 2 we present a method for the calculation of Fourier transformation  1 The outline iω x −1 f (x)e dx by expansions (2.1) and (2.4). In Sections 3 and 4, the proposed method is applied to compute the oscillatory integrals of the forms 2–17 in Table 1 with algebraic or logarithmic singularities at the end or interior points of the interval of integration. Finally, we finish in Section 5 by presenting some numerical examples to show the accuracy of the technique and the correctness of the theory. 2. On the calculation of Fourier transformation −1 f (x)eiωx dx

1

Let a sufficiently smooth f (x) be expanded as a series of Chebyshev polynomials of the form (see [24, Chapter 4], [25, Chapter 5] and [26]), f ( x) =

∞ 

′

ak Tk (x),

(2.1)

k=0

where the prime, as usual, indicates that the first term is to be halved, and ak =

2



1

f (x)Tk (x)

√

π

1 − x2

−1

dx.

Moreover, the first N coefficients ak can be computed efficiently by the FFT in O(N log N ) operations (see [24, Chapter 4] and [26]). In the sequel, we begin the exposition by presenting a series of expansions. From Abramowitz and Stegun [27, p. 361], it is shown that cos(z cos θ ) = J0 (z ) + 2

∞  (−1)−k J2k (z ) cos(2kθ );

(2.2)

k =1

sin(z cos θ ) = 2

∞  (−1)−k J2k+1 (z ) cos((2k + 1)θ ).

(2.3)

k=1

By using the substitution x = cos θ , x ∈ [−1, 1] and the trigonometric identity Tk (x) = cos kθ together with (2.2)–(2.3), we have by observation eiωx = cos(ω cos θ ) + i sin(ω cos θ )

= J0 (ω) + 2

∞ ∞   (−1)−j J2j (ω) cos(2jθ ) + 2i (−1)−j J2j+1 (ω) cos((2j + 1)θ ) j =1

=2

∞ 

j =1

i Jj (ω)Tj (x),

′ j

(2.4)

j =0

where Jj (ω) is the Bessel function of the first kind and of order j. Similarly to Clenshaw–Curtis indefinite approach described in [13], a function F (x) given in terms of Chebyshev polynomials is constructed, such that F (x) =



f (x)eiωx dx.

Now, substituting (2.1) and (2.4) into the indefinite integral F (x), we have F (x) = 2

∞  ∞  ′

j=0

k=0

′

ak ij Jj (ω)



Tj (x)Tk (x)dx.

144

H. Kang et al. / Journal of Computational and Applied Mathematics 242 (2013) 141–156

Owing to some properties of Chebyshev polynomials (see [25, Chapter 2]), Tn (1) = 1,

Tn (−1) = (−1)n , 1

Tk (x)Tj (x) =

2

(Tk+j (x) + T|k−j| (x)),

(2.5)

and



Tn (x)dx =

   1 Tn+1 (x) Tn−1 (x)   − ,   2 n+1 n−1

if n ≥ 2,

1

 T2 (x),   4 T1 (x),

if n = 1, if n = 0,

we then obtain by direct calculation 1



f (x)eiωx dx ≈ 2

I1 [ f ] = −1

N1 ∞  

′

j =0

′

ak ij Jj (ω)Ck,j ,

(2.6)

k=0

where 1



Tk (x)Tj (x)dx

Ck,j =

=

−1 0, 

1

−

(k + j)2 − 1

+

1

(k − j)2 − 1



if |k − j| is odd,

,

if |k − j| is even.

In (2.6), a sequence of Bessel function values Jj (ω)(j = 0, 1, . . . , N1 ), are required. Fortunately, the Bessel functions satisfy the three term recurrence relation, which should absolutely be used for computing the whole sequence of values Jj (ω). The recurrence relation Jn+1 (x) − 2n J (x) + Jn−1 (x) = 0, can be used to compute J0 (x), J1 (x), J2 (x), . . . , successively x n provided that n < x, otherwise severe accumulation of rounding errors will occur. Since, however, Jn (x) is a fast decreasing function of n when n > x, recurrence can always be carried out in the direction of decreasing n. Here, the two initial values J0 (ω), J1 (ω) for the three term recurrence relation can be obtained by truncating the asymptotic expansions after a suitable number of terms

 Jj (ω) ≈

2



πω

 K K   k A2k (j) k A2k+1 (j) (−1) (−1) cos x − sin x , ω2k ω2k+1 k=0 k=0

if j ≪ ω,

(2.7)

where A0 (ν) = 1,

(4ν 2 − 12 )(4ν 2 − 32 ) · · · (4ν 2 − (2k − 1)2 ) , k!8k  

Ak (ν) =

1

x=ω−

2

j+

1

4

k = 1, 2, 3, . . . ,

π.

For more details, one can refer to, e.g., [28] and [27, p. 385]. Obviously, when programming the proposed algorithm in a language like Matlab, we can use the built-in subroutine to calculate the values of Bessel functions. Especially, in the new versions of Matlab, the ‘besselj(j, ω)’ function is vectorized with respect to both arguments. In other words, we can compute the whole sequence J0 (ω), J1 (ω), . . . , JN1 (ω) with a single subroutine call. In fact, we still need to truncate the expansion (2.6) after an appropriate number of terms. When j is a little larger than ω, we may find that Jj (ω) can decrease to zero quite rapidly. Therefore, choosing N1 to be a little larger than ω, is a good strategy for obtaining the final approximations of such integrals. Moreover, we will provide an approximation to the proper value of N1 in Section 5. 3. Computation of the Fourier integrals with algebraic or logarithmic singularities at the end-points 3.1. Algebraic singularities at the endpoints In this subsection we first consider the problem of computing the integrals with algebraic singularities at the endpoints



1

I2 [ f ] =

(1 − x)α (1 + x)β f (x)eiωx dx,

−1

where α, β > −1.

(3.8)

H. Kang et al. / Journal of Computational and Applied Mathematics 242 (2013) 141–156

145

Likewise, substituting (2.1) and (2.4) into (3.8) yields

I2 [ f ] ≈ 2

N1 N  

ak ij Jj (ω)dk,j ,

(3.9)

(1 − x)α (1 + x)β Tk (x)Tj (x)dx,

(3.10)

′

j=0

′

k=0

where 1

 dk,j =

−1

which, using (2.5), can be written as dk,j =

1 2

[Mk[2+]j + M|[k2−] j| ].

(3.11)

Here the modified moments Mn[2] =

1



(1 − x)α (1 + x)β Tn (x)dx,

(3.12)

−1

can be computed efficiently by the three-term recurrence relation [14]

(α + β + n + 2)Mn[2+]1 + 2(α − β)Mn[2] + (α + β − n + 2)Mn[2−]1 = 0.

(3.13)

The main results, given by the four formulas (3.9)–(3.12), are common for all the considered integrals in Table 1. Therefore, in this paper, for each of the considered integrals, the similar results are not repeated over and over again. In [14], it is indicated that for (3.13) the forward recursion is perfectly numerically stable, except in two cases: 1 1 3 5

(a) α < β and α = − , , , , . . . ; 2 2 2 2 1 1 3 5

(b) α > β and β = − , , , , . . . . 2 2 2 2

Generally, by using the change of variable x = 1 − 2t and the definition and properties of the Beta function [27] B(a, b) =

1



xa−1 (1 − x)b−1 dx, 0

B(a, b + 1) =

b a+b

B(a, b),

[2]

[2]

the two starting values M0 and M1 may be expressed easily by [2]

M0



1

=

(1 − x)α (1 + x)β dx = 2α+β+1 B(α + 1, β + 1),

(3.14)

−1 [2]

M1



1

=

(1 − x)α (1 + x)β xdx =

−1

β −α [2] M , α+β +2 0

(3.15)

which can be computed explicitly since the Beta function B(a, b) can be implemented in Matlab by invoking the built-in function ‘beta(a, b)’. Remark 1. Using the representation of Tn (x) in terms of hypergeometric functions and Fasenmyer’s technique (see Rainville [29, Chapter 14]), the authors of [14] derived the three-term recurrence formula (3.13). By defining the integral (3.17), using properties of Tn (x) and integration by parts, a simpler proof is given as follows. Proof of (3.13). We have known the following elegant properties of Tn (x) [25]: 2xTn (x) = Tn+1 (x) + Tn−1 (x),

(1 − x )Tn (x) = 2

′

n 2

(Tn−1 (x) − Tn+1 (x)).

(3.16a) (3.16b)

146

H. Kang et al. / Journal of Computational and Applied Mathematics 242 (2013) 141–156

Let 1

 U =

(1 − x)α+1 (1 + x)β+1 Tn′ (x)dx.

(3.17)

−1

Using the identity (3.16b) and (3.12), Eq. (3.17) can be written in the form 1



(1 − x)α (1 + x)β (1 − x2 )Tn′ (x)dx

U = −1

=

 n  [2] [2] Mn−1 − Mn+1 . 2

(3.18)

On the other hand, integration by parts of Eq. (3.17) yields U = (α + 1)Mn[2] (α, β + 1) − (β + 1)Mn[2] (α + 1, β).

(3.19)

In fact, using the identity (3.16a), we obtain Mn[2] (α, β + 1) =



1

(1 − x)α (1 + x)β+1 Tn (x)dx

−1



1

=

(1 − x)α (1 + x)β Tn (x)dx +



−1

1

(1 − x)α (1 + x)β xTn (x)dx

−1

1

= Mn[2] + (Mn[2−]1 + Mn[2+]1 ).

(3.20)

2

Accordingly, it is easy to attain Mn[2] (α + 1, β) = Mn[2] −

1 2

(Mn[2−]1 + Mn[2+]1 ).

(3.21)

Substituting (3.20) and (3.21) into (3.19) together with (3.18), we complete the desired proof.



Remark 2. In [30], three methods for constructing recurrence relations for modified moments were proposed. Using the method I presented in [30], we can obtain lower order recurrence relations for modified moments Mn[c ] considered in the paper. Also, in the proof of (3.13), we use the same relations involving Chebyshev polynomials as used in the method given in [30]. There are similarities among their ideas. For more details one can refer to [30,29,14]. 3.2. Algebraic and logarithmic singularities at the endpoints The succeeding part is used to illustrate how the presented method can be employed to compute the integrals containing algebraic or logarithmic singularities at the endpoints



1

I3 [ f ] =

ln(1 − x)(1 − x)α (1 + x)β f (x)eiωx dx,

(3.22)

(1 − x)α (1 + x)β ln(1 + x)f (x)eiωx dx,

(3.23)

ln(1 − x)(1 − x)α (1 + x)β ln(1 + x)f (x)eiωx dx,

(3.24)

−1



1

I4 [ f ] = −1



1

I5 [ f ] = −1

where α, β > −1. Similarly, to compute the integrals (3.22)–(3.24), it is required to compute the modified moments Mn[c ] , c = 3, 4, 5, respectively. It is worth noting that

∂ [2] M = Mn[3] , ∂α n

∂ [2] M = Mn[4] , ∂β n

∂2 M [2] = Mn[5] . ∂α∂β n

Therefore, by differentiating the three-term recurrence relation (3.13) with respect to α or β , we find Mn[c ] , c = 3, 4, 5 satisfying the following three recurrence relations:

(α + β + n + 2)Mn[c+] 1 + 2(α − β)Mn[c ] + (α + β − n + 2)Mn[c−] 1 = rn[c ] , [c ]

where c is the number of the integral from Table 1 and the coefficients rn are gathered in Table 2:

(3.25)

H. Kang et al. / Journal of Computational and Applied Mathematics 242 (2013) 141–156

147

Table 2 The coefficients rn[c ] of the recurrence relations (3.25). c

rn[c ]

2 3 4 5

−Mn[2+]1 − 2Mn[2] − Mn[2−]1 −Mn[2+]1 + 2Mn[2] − Mn[2−]1 −Mn[3+]1 +2Mn[3] −Mn[3−]1 −Mn[4+]1 −2Mn[4] −Mn[4−]1

0

• Derivation of the initial moments: By using derivatives of the Beta function [31], d da d

B(a, b) = B(a, b)[ψ0 (a) − ψ0 (a + b)],

B(a, b) = B(a, b)[ψ0 (b) − ψ0 (a + b)], db 2 d B(a, b) = B(a, b){[ψ0 (a) − ψ0 (a + b)] × [ψ0 (b) − ψ0 (a + b)] − ψ1 (a + b)}, dadb where ψ0 (z ) is the Psi function or Digamma function which can be implemented in Matlab by calling ‘psi(z )’, and the Polygamma function (for definition, see [27, pp. 258–260]),

ψn (z ) =

dn+1 dz n+1

ln[Γ (z )] =

dn Γ ′ (z ) dz n

Γ (z )

=

dn dz n

ψ0 (z ),

can be implemented in Matlab by calling ‘psi(n, z )’, the required starting values can be attained by differentiating (3.14) and (3.15) with respect to α or β , as follows: [3]

M0

[3]

M1

[4]

M0

[4]

M1

[5]

M0

[5]

M1

= M0[2] [ln 2 + ψ0 (α + 1) − ψ0 (α + β + 2)], β −α 2(β + 1) [2] [3] M + M , =− (α + β + 2)2 0 α+β +2 0 = M0[2] [ln 2 + ψ0 (β + 1) − ψ0 (α + β + 2)], β −α 2(α + 1) [2] [4] M0 + M , = 2 (α + β + 2) α+β +2 0  [2] = M0 (ln 2)2 + [ψ0 (α + 1) + ψ0 (β + 1) − 2ψ0 (α + β + 2)] ln 2  + [ψ0 (α + 1) − ψ0 (α + β + 2)] × [ψ0 (β + 1) − ψ0 (α + β + 2)] − ψ1 (α + β + 2) ,   β −α β +1 [2] = 2M0 − [ln 2 + ψ0 (β + 1) − ψ0 (α + β + 2)] (α + β + 2)3 (α + β + 2)2 β −α 2(α + 1) [3] [4] + M + M . (α + β + 2)2 0 α+β +2 0

The three recurrence relations (3.25) for c = 3, 4, 5 are non-homogeneous recurrence relations, the homogeneous parts of which are just the three-term recurrence relation (3.13). Following the idea of the third section in [14], the remarks concerning the numerical stability of forward recursion for (3.13), remain valid here. It is then an easy matter to use forward recursion for the recurrence relations (3.25) for c = 2, 3, 4, 5 to compute Mn[c ] , c = 2, 3, 4, 5 efficiently. 4. Extensions to the Fourier integrals with algebraic or logarithmic singularities at the interior or end points 4.1. Singularity at only one interior point In this section we first investigate the calculation of the Fourier integrals with algebraic or logarithmic singularities at one interior point



1

I6 [ f ] =

|x − a|α f (x)eiωx dx,

(4.26)

|x − a|α sgn(x − a)f (x)eiωx dx,

(4.27)

−1



1

I7 [ f ] = −1

148

H. Kang et al. / Journal of Computational and Applied Mathematics 242 (2013) 141–156 Table 3 The coefficients rn[c ] of the recurrence relations (4.30).



1

I8 [ f ] =

c

rn[c ]

6 7 8 9

2[(−1)n (1 + a)α+1 − (1 − a)α+1 ] 2[(−1)n+1 (1 + a)α+1 − (1 − a)α+1 ] [6] [6] 2[(−1)n (a + 1)α+1 ln(1 + a)−(1 − a)α+1 ln(1 − a)]−(n − 1)Mn+1 +(n + 1)Mn−1 [7] [7] 2[(−1)n+1 (a + 1)α+1 ln(1 + a)−(1 − a)α+1 ln(1 − a)]−(n − 1)Mn+1 +(n + 1)Mn−1

|x − a|α ln |x − a|f (x)eiωx dx,

(4.28)

|x − a|α ln |x − a|sgn(x − a)f (x)eiωx dx,

(4.29)

−1



1

I9 [ f ] = −1

where α, β, γ > −1 and |a| < 1. Similarly, to compute the integrals (4.26)–(4.29), it is required to compute the modified moments Mn[c ] , c = 6, 7, 8, 9, which satisfy the recurrence relations

(n − 1)(n + α + 2)Mn[c+] 1 − 2a(n2 − 1)Mn[c ] + (n + 1)(n − α − 2)Mn[c−] 1 = rn[c ] ,

(4.30)

[c ]

where c is the number of the integral from Table 1 and the coefficients rn are shown in Table 3: • Derivation of the initial moments: By integration by parts, the starting values for the recurrence relations (4.30) for c = 6, 7, can be written as 1

[6]

M0 = [6]

M1

[7]

M0

[7]

M1

= = =

α+1 1

α+2 1

α+1 1

α+2

[(1 − a)α+1 + (1 + a)α+1 ], [(1 − a)α+2 − (1 + a)α+2 ] + aM0[6] , [(1 − a)α+1 − (1 + a)α+1 ], [(1 − a)α+2 + (1 + a)α+2 ] + aM0[7] .

By differentiating M0 , M1 represented by [6]

[6]

with respect to α , the starting values for the recurrence relation (4.30) for c = 8, can be

[8]

= G(α, a) + G(α, −a),

[8]

= G(α + 1, −a) − G(α + 1, a) + aM0[8] ,

M0 M1 where

G(α, a) =

(1 + a)α+1 (1 + a)α+1 ln(1 + a) − . α+1 (α + 1)2

By differentiating M0 , M1 with respect to α respectively, the starting values for the recurrence relation (4.30) for c = 9, can be expressed by [7]

[7]

[9]

= G(α, −a) − G(α, a),

[9]

= G(α + 1, −a) + G(α + 1, a) + aM0[9] .

M0 M1

The theoretical proof of the stability of the forward recursion in the case of the recurrence relation (4.30) is not known yet, and so we decide to verify the stability in an experimental way. This can be seen in Table 4 where we report absolute errors of forward recursion applied on the recurrence relation (4.30). For c = 6, Table 4 shows numerical stability of the recurrence relations (4.30) for different parameters. Analogously, for c = 7, 8, 9, numerical stability of the recurrence relations (4.30) remains valid. The similar numerical experiments are not repeated here. 4.2. Singularities at the endpoints and one point in the interior of interval The part is devoted to describing how the proposed method is used to deal with the computation of more complicated integrals with algebraic or logarithmic singularities at the interior or end points



1

I10 [f ] = −1

(1 − x)α (1 + x)β |x − a|γ f (x)eiωx dx,

(4.31)

H. Kang et al. / Journal of Computational and Applied Mathematics 242 (2013) 141–156

149

Table 4 Absolute errors for computing Mn[6] by forward recursion for the recurrence relations (4.30). a

α

n = 5000

n = 10000

n = 20000

n = 30000

0

− 31

2.0 × 10−15

1.4 × 10−15

2.0 × 10−15

4.3 × 10−15

− 12 0.99

4.6 × 10−15

2.2 × 10−15

2.5 × 10−15

1.9 × 10−15

−

1 3

1.8 × 10

−15

3.5 × 10

6.9 × 10

−15

7.7 × 10−15

−

1 2

5.3 × 10

−15

4.7 × 10

3.4 × 10

−15

8.3 × 10−15

−15 −15

Table 5 The coefficients rn[c ] of the recurrence relations (4.39). c

rn[c ]

10 11 12 13 14 15 16 17

0 [10] [10] 2Mn[10] − Mn+2 − Mn−2 [10] [10] [10] [10] 2(a − 1)(Mn+1 + Mn−1 ) + 2(2a − 1)Mn[10] − Mn+2 − Mn−2 [10] [10] [10] [10] [10] 2(a + 1)(Mn+1 + Mn−1 ) − 2(2a + 1)Mn − Mn+2 − Mn−2 [12] [12] [12] [12] [13] [13] [13] [13] 2(a + 1)(Mn−1 + Mn+1 ) − 2(2a + 1)Mn[12] − Mn+2 − Mn−2 + 2(a − 1)(Mn+1 + Mn−1 ) + 2(2a − 1)Mn[13] − Mn+2 − Mn−2 [11] [11] [11] [11] [13] [13] 2(a + 1)(Mn+1 + Mn−1 ) − 2(2a + 1)Mn[11] − Mn+2 − Mn−2 + 2Mn[13] − Mn+2 − Mn−2 [11] [11] [11] [11] [12] [12] 2(a − 1)(Mn+1 + Mn−1 ) + 2(2a − 1)Mn[11] − Mn+2 − Mn−2 + 2Mn[12] − Mn+2 − Mn−2 [16] [16] [16] [16] [15] [15] [15] [15] [14] [14] 2(a+1)(Mn−1 +Mn+1 )−2(2a+1)Mn[16] −Mn+2 −Mn−2 +2(a−1)(Mn+1 +Mn−1 )+2(2a−1)Mn[15] −Mn+2 −Mn−2 +2Mn[14] −Mn+2 −Mn−2



1

I11 [f ] =

(1 − x)α (1 + x)β |x − a|γ ln |x − a|f (x)eiωx dx,

(4.32)

ln(1 − x)(1 − x)α (1 + x)β |x − a|γ f (x)eiωx dx,

(4.33)

(1 − x)α (1 + x)β ln(1 + x)|x − a|γ f (x)eiωx dx,

(4.34)

ln(1 − x)(1 − x)α (1 + x)β ln(1 + x)|x − a|γ f (x)eiωx dx,

(4.35)

(1 − x)α (1 + x)β ln(1 + x)|x − a|γ ln |x − a|f (x)eiωx dx,

(4.36)

ln(1 − x)(1 − x)α (1 + x)β |x − a|γ ln |x − a|f (x)eiωx dx,

(4.37)

ln(1 − x)(1 − x)α (1 + x)β ln(1 + x)|x − a|γ ln |x − a|f (x)eiωx dx,

(4.38)

−1



1

I12 [f ] = −1



1

I13 [f ] = −1



1

I14 [f ] = −1



1

I15 [f ] = −1



1

I16 [f ] = −1



1

I17 [f ] = −1

where α, β, γ > −1 and |a| < 1. Similarly, to compute the integrals (4.31)–(4.38), it is required to compute the modified moments Mn[c ] , c 10, 11, . . . , 17, which satisfy the recurrence relations

=

(α + β + γ + n + 3)Mn[c+] 2 + 2[(1 − a)(α + 1) − (1 + a)(β + 1) − na]Mn[c+] 1 + 2[(1 − 2a)(α + 1) + (1 + 2a)(β + 1) − (γ + 1)]Mn[c ] + 2[(1 − a)(α + 1) − (1 + a)(β + 1) + na]Mn[c−] 1 + (α + β + γ − n + 3)Mn[c−] 2 = rn[c ] , [c ]

(4.39)

where c is the number of the integral from Table 1 and the coefficients rn are listed in Table 5: • Derivation of the initial moments: Because of the symmetry of the recurrence relation of the Chebyshev polynomials Tn (x), it is convenient to set T−n (x) = [10] Tn (x), n = 1, 2, . . . , and consequently M−n = Mn[10] . Therefore, when the starting values Mn[10] , n = 1, 2, are known, [10] Mn , n = −1, −2, can be obtained accordingly. Furthermore, numerical experiments show that all Mn[10] can be computed accurately by the forward recursion. First, we shall study how to compute the three starting values Mn[10] , n = 0, 1, 2. By using the integral representation of Gauss hypergeometric function 2 F1 (µ, ν; λ; z ) [27, p. 558] 2 F1 (µ, ν; λ; z ) =

1 B(ν, λ − ν)

 0

1

xν−1 (1 − x)λ−ν−1 (1 − zx)−µ dx,

150

H. Kang et al. / Journal of Computational and Applied Mathematics 242 (2013) 141–156

where ℜ(λ) > ℜ(ν), provided |z | ≤ 1, and the change of variables x = −1 + (1 + a)t and x = 1 + (a − 1)t, the starting [10] value M0 can be denoted by [10]

M0

1

 =

(1 − x)α (1 + x)β |x − a|γ dx

−1 a



α

β

(1 − x) (1 + x) (a − x) dx +

= −1 α

1



γ

(1 − x)α (1 + x)β (x − a)γ dx

a

= 2 (1 + a)

β+γ +1



1

t

β+1−1



β+γ +2−(β+1)−1

(1 − t )

1−

1+a

0

+ 2β (1 − a)α+γ +1



1

t α+1−1 (1 − t )α+γ +2−(α+1)−1

 1−

0

2

−(−α) t

1−a 2

dt

−(−β) t

dt

= H (α, β, γ , a) + H (β, α, γ , −a),

(4.40)

where H (α, β, γ , a) = 2α (1 + a)β+γ +1 B(β + 1, 1 + γ )2 F1 (−α, 1 + β; β + γ + 2; (1 + a)/2). Since ℜ(β + γ + 2) > ℜ(1 + β) and 0 < Substituting

1+a 2

(4.41)

< 1, the definition of (4.41) is valid.

T1 (x) = (1 + x) − 1, and T2 (x) = 2x2 − 1 = 1 − 2(1 + x)(1 − x), into, [10]

M1



1

=

(1 − x)α (1 + x)β |x − a|γ T1 (x)dx,

−1 [10]

M2



1

=

(1 − x)α (1 + x)β |x − a|γ T2 (x)dx,

−1

respectively, and following the idea of obtaining (4.40), we have by direct calculation [10]

M1

= M0[10] − H (α + 1, β, γ , a) − H (β, α + 1, γ , −a) = H (α, β + 1, γ , a) + H (β + 1, α, γ , −a) − M0[10] ,

[10]

M2

= M0[10] − 2[H (α + 1, β + 1, γ , a) + H (β + 1, α + 1, γ , −a)] = M0[10] + 2[H (α, β + 2, γ , a) + H (β + 2, α, γ , −a)] − 4[H (α, β + 1, γ , a) + H (β + 1, α, γ , −a)] = 2[H (α, β + 2, γ , a) + H (β + 2, α, γ , −a)] + 4[H (α + 1, β, γ , a) + H (β, α + 1, γ , −a)] − 7M0[10] = M0[10] + 2[H (α + 2, β, γ , a) + H (β, α + 2, γ , −a)] − 4[H (α + 1, β, γ , a) + H (β, α + 1, γ , −a)] = 2[H (α + 2, β, γ , a) + H (β, α + 2, γ , −a)] + 4[H (α, β + 1, γ , a) + H (β + 1, α, γ , −a)] − 7M0[10] .

In order to find the values of the initial moments M0 , M1 , M2 , c = 11, 12, . . . , 17, it is necessary to compute derivatives of the Gauss hypergeometric function 2 F1 with respect to its parameters. From Abramowitz and Stegun [27, p. 556], we know that the Gauss hypergeometric function 2 F1 (µ, ν; λ; z ) is defined for the circle of convergence of |z | < 1 by the power series [c ]

2 F1

(µ, ν; λ; z ) =

[c ]

[c ]

∞ Γ (λ)  Γ (µ + n)Γ (ν + n) z n , Γ (µ)Γ (ν) n=0 Γ (λ + n) n!

provided that λ is not 0, −1, −2, . . .. The behavior of this series on its circle of convergence is: (a) Divergence when ℜ(λ − µ − ν) ≤ −1; (b) Absolute convergence when ℜ(λ − µ − ν) > 0; (c) Conditional convergence when −1 < ℜ(λ − µ − ν) ≤ 0; the point z = 1 is excluded.

(4.42)

H. Kang et al. / Journal of Computational and Applied Mathematics 242 (2013) 141–156

151

The explicit differentiation formulas with respect to parameters µ, ν, λ, z are discussed in [27,32–35]. The first derivatives for only some special values of the parameters are presented in [22, p. 85] and references therein. If µ is a negative integer, the first derivatives with respect to ν or λ can be found in [33]. Recently, the derivatives of higher order, such as any n-th derivatives with respect to single parameter, have been studied in detail in [35]. The required two derivatives also appear in the website [34], as follows: d dµ

2 F1

(µ, ν; λ; z ) =

∞ zn Γ (λ)  Γ (µ + n)Γ (ν + n) [ψ0 (µ + n) − ψ0 (µ)] Γ (µ)Γ (ν) n=0 Γ (λ + n) n!

∞ Γ (λ)  Γ (µ + n)Γ (ν + n) zn ψ0 (µ + n) Γ (µ)Γ (ν) n=0 Γ (λ + n) n!   ∞ n−1 Γ (λ)  Γ (µ + n)Γ (ν + n) z n  1 , = Γ (µ)Γ (ν) n=0 Γ (λ + n) n! k=0 µ + k

= −ψ0 (µ)2 F1 (µ, ν; λ; z ) +

d dλ

2 F1

(µ, ν; λ; z ) = −

(4.43)

∞ Γ (λ)  Γ (µ + n)Γ (ν + n) zn [ψ0 (λ + n) − ψ0 (λ)] Γ (µ)Γ (ν) n=0 Γ (λ + n) n! ∞ Γ (λ)  zn Γ (µ + n)Γ (ν + n) ψ0 (λ + n) Γ (µ)Γ (ν) n=0 Γ (λ + n) n!   ∞ n−1  Γ (µ + n)Γ (ν + n) z n  1

= ψ0 (λ)2 F1 (µ, ν; λ; z ) − =−

Γ (λ) Γ (µ)Γ (ν)

n =0

Γ (λ + n)

n! k=0 λ + k

 from [27, p. 258], ψ0 (z + n) − ψ0 (z ) =

n−1  k=0

1 z+k

 .

(4.44)

Similarly, we have d2

2 F1 (µ, ν; λ; z ) =

dµdν d2 dµdλ

2 F1

∞ Γ (λ)  zn Γ (µ + n)Γ (ν + n) [ψ0 (µ + n) − ψ0 (µ)][ψ0 (ν + n) − ψ0 (ν)] Γ (µ)Γ (ν) n=0 Γ (λ + n) n!

(µ, ν; λ; z ) = −

∞ Γ (λ)  Γ (µ + n)Γ (ν + n) Γ (µ)Γ (ν) n=0 Γ (λ + n)

× [ψ0 (µ + n) − ψ0 (µ)][ψ0 (λ + n) − ψ0 (λ)] d2 dν dλ

2 F1

(µ, ν; λ; z ) = −

dµdν dλ

2 F1

(4.46)

n!

∞ Γ (λ)  Γ (µ + n)Γ (ν + n) [ψ0 (ν + n) − ψ0 (ν)] Γ (µ)Γ (ν) n=0 Γ (λ + n)

× [ψ0 (λ + n) − ψ0 (λ)] d3

zn

(4.45)

(µ, ν; λ; z ) = −

zn

(4.47)

n!

∞ Γ (λ)  Γ (µ + n)Γ (ν + n) Γ (µ)Γ (ν) n=0 Γ (λ + n)

× [ψ0 (µ + n) − ψ0 (µ)][ψ0 (ν + n) − ψ0 (ν)][ψ0 (λ + n) − ψ0 (λ)]

zn

n! (from[27, p. 258], Γ ′ (z ) = Γ (z )ψ0 (z )). (4.48) As µ and ν are symmetric in (4.42), the first derivative with respect to ν may then be obtained by exchanging µ and ν in (4.43). For other higher order derivatives needed here, one can refer to [34,35]. Since all these power series in (4.43)–(4.48) converge uniformly inside the disk |z | ≤ ρ < 1, the required derivative values can be obtained by truncating after a suitable number of terms and assigning values to parameters under the given conditions that z is not close to the endpoints −1 and 1 for all α, β, γ > −1. For numerical computations we can use only the disk |z | ≤ ρ < 1, with ρ depending on numerical requirements, such as precision and efficiency [36, p. 30]. Note that for 2 F1 (−α, 1 + β; β + γ + 2; (1 + a)/2), one has λ − µ − ν = α + γ + 1 > −1, 0 < |z | = (1 + a)/2 < 1 and λ = β + γ + 2 ̸= 0, −1, −2, . . . . Following (4.43)–(4.48), the derivatives of 2 F1 (−α, 1 + β; β + γ + 2; (1 + a)/2) with respect to α, β, γ can be obtained accordingly by truncating after a suitable number of terms when a is not close to the endpoints −1 and 1 for |a| < 1. Consequently, by using derivatives of the Beta function in Section 3.2 and differentiating [10] [10] [10] [c ] [c ] [c ] M0 , M1 , M2 with respect to α, β, γ , we can get the initial moments M0 , M1 , M2 , c = 11, 12, . . . , 17, for a not being close to −1 and 1. However, when a is close to the endpoints −1 and 1 for |a| < 1, we will develop an efficient and stable method to compute the initial moments in our future work.

152

H. Kang et al. / Journal of Computational and Applied Mathematics 242 (2013) 141–156 Table 6 Absolute errors for computing Mn[c ] (c = 10) by forward recursion for the recurrence relations (4.39). n = 1000

n = 10000

n = 20000

n = 40000

n = 100000

β = − 21 , γ = − 13

4.8 × 10−15

2.6 × 10−15

3.7 × 10−15

6.6 × 10−15

7.6 × 10−15

β = − 31 , γ = − 14

9.7 × 10−15

5.3 × 10−15

1.9 × 10−15

7.4 × 10−15

8.1 × 10−14

a = 0, α = − 14 ,

a = 0.5, α = − 12 ,

Table 7 Comparison (absolute errors and computation time in seconds) of the proposed method and n-point Gauss–Legendre rule for I2 (α = −1/3, β = −1/4, f (x) = 2

cos x) and I5 (α = −1/4, β = −1/3, f (x) = xe−x ). Absolute errors Time Absolute errors Time Absolute errors Time Absolute errors Time

E [2] (12, 267, 200) 4.9 × 10−15 0.415

EGL (3200, 200) 5.7 × 10−6 0.406

EGL (100000, 200) 6.0 × 10−8 239.578

E [2] (12, 477, 400) 4.0 × 10−15 0.722

EGL (5000, 400) 3.4 × 10−6 0.719

EGL (100000, 400) 6.0 × 10−8 232.000

E [5] (20, 477, 400) 4.4 × 10−15 1.246

EGL (7000, 400) 2.8 × 10−5 1.407

EGL (100000, 400) 8.4 × 10−7 240.656

E [5] (20, 897, 800) 4.5 × 10−15 2.337

EGL (7200, 800) 3.2 × 10−5 2.428

EGL (100000, 800) 8.5 × 10−7 245.578

[2]

[2]

[5]

[5]

[2]

[2]

[5]

[5]

For c = 10, Table 6 shows numerical stability of forward recursion for the recurrence relations (4.39) for different parameters. Likewise, for c = 11, 12, . . . , 17 and the cases of a not being close to −1 and 1, the similar numerical experiments can also verify numerical stability of forward recursion for the recurrence relations (4.39). Remark 3. Similarly, by constructing more complicated recurrence relations, the integrals −1 (ln(1 − x)s )p (1 − x)α (1 + x)β (ln(1 + x)t )q |x − a|γ (ln |x − a|m )r f (x)eiωx dx, α, β, γ > −1, s, t , m ∈ R, |a| < 1 and p, q, r = 0, 1, 2, . . . , can be also computed accurately by the proposed method. When the singularities become stronger (if p, q, r are large) or ω is very large, the presented method may be the only tool for accurate computation of such integrals.

1

5. Numerical examples In this section, we present the results of numerical experiments obtained using the proposed method. Our algorithm is compared to the n-point Gauss–Legendre rule. In order to conduct the experiments, we require the knowledge on the exact values of the integrals of the form (1.1). The values that we assume to be accurate are computed in Maple 13.02 using the 32 decimal digits precision arithmetic. The compared algorithms (the Gauss–Legendre rule and the presented one) are implemented in Matlab 7.0.1, taking advantage of its fast vectorized arithmetic operations. The experiments are performed on the desktop computer with AMD Athlon(tm) 64 X2 Dual Core Processor 4000+(2100 MHz) and 1 GB of RAM. [c ] By E [c ] (N , N1 , ω), EGL (n, ω), where c is the number of the integral from Table 1, we denote the absolute errors of the approximations to the integrals (1.1) computed by the proposed method and the n-point Gauss–Legendre rule, respectively. Here, we use the Matlab code written by Greg von Winckel in [37] to find the nodes and weights of the Gauss–Legendre rules. The process of determining the nodes and weights is included in the time measurements. Moreover, the order n of Gauss–Legendre rules is selected randomly, e.g. n = 100000 for not very large ω ≤ 10000. To judge the efficiency of the proposed method, computation time in seconds is also given below. Here, we provide the linear approximation, for example, N1 ≡ N1 (ω) = ⌊1.05ω + 57.5⌋.

(5.49)

This can guarantee that JN1 (ω) ≤ 10 for all ω > 0, and when ω is large, the optimal value of N1 is overestimated by only a few percent. In the following examples, we give approximations to the proper values of N1 by using (5.49). −15

Example 1. For the integrals Ic = −1 (1 − x)α (1 + x)β f (x)eiωx dx(c = 2), and Ic = −1 ln(1 − x)(1 − x)α (1 + x)β ln(1 + x)f (x)eiωx dx(c = 5), Table 7 shows the efficiency and accuracy of the proposed method, compared to the Gauss–Legendre rule.

1

1

Example 2. For Ic = −1 |x − a|α f (x)eiωx dx, (c = 6) and Ic = −1 |x − a|α ln |x − a|f (x)eiωx dx, (c = 8), Tables 8–9 show the efficiency and accuracy of the proposed method, compared to the Gauss–Legendre rule.

1

1

H. Kang et al. / Journal of Computational and Applied Mathematics 242 (2013) 141–156

153

Table 8 Comparison (absolute errors and computation time in seconds) of the proposed method and n-point Gauss–Legendre rule for I6 (α = −1/2, a = 1/2, f (x) = sin x) and I8 (α = −1/3, a = 1/2, f (x) = ex ). Absolute errors Times Absolute errors Time Absolute errors Time Absolute errors Time

E [6] (12, 897, 800) 7.7 × 10−15 1.462

EGL (7200, 800) 1.1 × 10−2 1.656 [6]

EGL (100000, 200) 2.1 × 10−3 246.812 [6]

E [6] (12, 1737, 1600) EGL (11000, 1600) 1.1 × 10−14 4.2 × 10−3 2.873 3.313

EGL (100000, 1600) 2.1 × 10−3 243.954

E [8] (12, 897, 800) 6.3 × 10−14 1.576

EGL (100000, 800) 8.0 × 10−3 292.141

[6]

EGL (7700, 800) 3.4 × 10−2 2.016 [8]

E [8] (12, 1737, 1600) EGL (11500, 1600) 8.1 × 10−14 3.0 × 10−2 2.477 4.375 [8]

[6]

[8]

EGL (100000, 1600) 8.0 × 10−3 291.375 [8]

Table 9 Comparison (absolute errors and computation time in seconds) of the proposed method and n-point Gauss–Legendre rule for I6 (α = −1/3, a = 0.99, f (x) = sin x) and I8 (α = −1/2, a = 0.99, f (x) = ex ). Absolute errors Time Absolute errors Time Absolute errors Time Absolute errors Time

E [6] (12, 897, 800) 2.1 × 10−14 1.365

EGL (7500, 800) 1.7 × 10−4 1.687 [6]

E [6] (12, 1737, 1600) EGL (11100, 1600) 3.4 × 10−14 4.5 × 10−4 2.431 3.391 [6]

E [8] (12, 897, 800) 7.2 × 10−13 1.412

EGL (7600, 800) 3.2 × 10−1 1.985 [8]

E [8] (12, 1737, 1600) EGL (10400, 1600) 1.2 × 10−12 1.5 × 10−1 2.490 3.516 [8]

EGL (100000, 200) 1.1 × 10−4 251.141 [6]

EGL (100000, 1600) 1.1 × 10−4 248.875 [6]

(100000, 800) 1.1 × 10−1 290.89 [8] GL

EGL (100000, 1600) 1.1 × 10−1 292.766 [8]

Table 10 Comparison (absolute errors and computation time in seconds) of the proposed method and n-point Gauss–Legendre rule for I10 (α = −1/3, β = −1/4, γ = −1/2, a = 0 and f (x) = cos x) and I17 (α = −1/3, β = −1/4, γ = −1/2, a = 0 and f (x) = cos x). Absolute errors Time Absolute errors Time Absolute errors Time Absolute errors Time

E [10] (12, 477, 400) 1.3 × 10−15 0.911

EGL (7500, 400) 1.1 × 10−2 1.703

EGL (100000, 400) 4.5 × 10−3 288.797

E [10] (12, 897, 800) 8.2 × 10−15 1.598

EGL (8800, 800) 1.5 × 10−2 2.750

EGL (100000, 800) 4.5 × 10−3 295.437

E [17] (12, 477, 400) 3.4 × 10−14 2.133

EGL (8000, 400) 4.4 × 10−2 2.216

EGL (100000, 400) 8.8 × 10−3 295.141

E [17] (12, 897, 800) 9.7 × 10−14 4.121

EGL (11800, 800) 3.8 × 10−2 4.665

EGL (100000, 800) 7.3 × 10−3 298.375

[10]

[10]

[17]

[17]

[10]

[10]

[17]

[17]

Example 3. For Ic = −1 (1 − x)α (1 + x)β |x − a|γ f (x)eiωx dx, (c = 10) and Ic = −1 ln(1 − x)(1 − x)α (1 + x)β ln(1 + x)|x − a|γ ln |x − a|f (x)eiωx dx, (c = 17), Tables 10, 11 and 13 show the efficiency and accuracy of the proposed method, compared to the Gauss–Legendre rule.

1

1

Table 12 shows that for I17 , when a = 0.99 is close to the endpoint 1, the proposed method cannot obtain high precision. This is because in this case z = (1 + a)/2 is close to the endpoint 1, and truncating (4.43)–(4.48), the obtained derivative values of the Gauss hypergeometric function 2 F1 with respect to its parameters are not very accurate, which makes the [17] [17] [17] values of the initial moments M0 , M1 , M2 unreliable. Similarly, for a close to −1 or 1, the truncation errors of the

154

H. Kang et al. / Journal of Computational and Applied Mathematics 242 (2013) 141–156 Table 11 Comparison (absolute errors and computation time in seconds) of the proposed method and n-point Gauss–Legendre rule for I10 (α = −1/3, β = −1/4, γ = −1/2, a = 0.99 and f (x) = cos x). Absolute errors Time Absolute errors Time

E [10] (12, 477, 400) EGL (7300, 400) EGL (100000, 400) 7.2 × 10−14 6.9 × 10−3 5.2 × 10−3 0.882 1.844 291.985 [10]

[10]

E [10] (12, 897, 800) EGL (9200, 800) EGL (100000, 800) 6.8 × 10−14 1.1 × 10−2 5.2 × 10−3 1.586 2.859 295.282 [10]

[10]

Table 12 Comparison (absolute errors and computation time in seconds) of the proposed method and n-point Gauss–Legendre rule for I17 (α = −1/3, β = −1/4, γ = −1/2, a = 0.99 and f (x) = cos x). Absolute errors Time Absolute errors Time

E [17] (12, 477, 400) EGL (8800, 400) 6.3 × 10−4 5.9 × 10−2 2.247 2.634 [17]

EGL (100000, 400) 6.4 × 10−3 291.985 [17]

E [17] (12, 897, 800) EGL (11000, 800) EGL (100000, 800) 7.5 × 10−4 2.3 × 10−2 8.1 × 10−3 3.278 3.675 298.373 [17]

[17]

Table 13 Comparison (absolute errors and computation time in seconds) of the proposed method and n-point Gauss–Legendre rule for I17 (α = −1/3, β = −1/4, γ = −1/2, a = 0.75 and f (x) = cos x). Absolute errors Time Absolute errors Time

E [17] (12, 477, 400) EGL (8800, 400) 3.9 × 10−13 3.9 × 10−2 2.597 2.634 [17]

EGL (100000, 400) 4.5 × 10−3 287.456 [17]

E [17] (12, 897, 800) EGL (11000, 800) EGL (100000, 800) 7.5 × 10−14 3.6 × 10−2 8.1 × 10−3 4.639 4.889 294.396 [17]

[17]

Table 14 Absolute errors and computation time in seconds of the proposed method for I3 (α = −1/2, β = −1/3), I9 (α = −1/2, a = 0) and I10 (α = −1/2, β = −1/3, γ =

−1/4, a = 0.99), with f (x) = 10e−(2x−1) (10(4x + 1)2 + 1)−1 and N = 360. 2

Absolute errors Time Absolute errors Time Absolute errors Time

E [3] (360, 477, 400) 3.4 × 10−15 21.231

E [3] (360, 897, 800) 6.7 × 10−15 39.126

E [9] (360, 582, 500) 8.3 × 10−15 23.786

E [9] (360, 1107, 1000) 4.5 × 10−15 46.145

E [10] (360, 477, 400) 3.2 × 10−15 21.637

E [10] (360, 897, 800) 9.1 × 10−15 38.472

power series from (4.43) to (4.48) are also dissatisfactory. In the future we will develop an efficient method to compute the [c ] [c ] [c ] initial moments M0 , M1 , M2 , c = 11, 12, . . . , 17, for a close to −1 or 1. In the above numerical Examples 1–3, for very simple f , we have compared the proposed method to the Gauss–Legendre rule. In the following we omit the comparison. Here, we consider a little more challenging function, such as, f (x) = 10e−(2x−1) (10(4x + 1)2 + 1)−1 , 2

(5.50)

which can be closely approximated by a polynomial of degree a few hundreds. Example 4. Let us consider I3 = −1 ln(1 − x)(1 − x)α (1 + x)β f (x)eiωx dx, I9 = −1 |x − a|α ln |x − a|sgn(x − a)f (x)eiωx dx, 1 I10 = −1 (1 − x)α (1 + x)β |x − a|γ f (x)eiωx dx for the function f given by (5.50), respectively (see Table 14).

1

1

Example 5. Let us consider I4 = −1 (1 − x)α (1 + x)β ln(1 + x)f (x)eiωx dx, I10 = −1 (1 − x)α (1 + x)β |x − a|γ f (x)eiωx dx for large ω, respectively (see Table 15).

1

1

H. Kang et al. / Journal of Computational and Applied Mathematics 242 (2013) 141–156

155

Table 15 Absolute errors and computation time in seconds of the proposed method for I10 (α = −1/4, β = −1/3, γ = −1/2, a = 0 and f (x) = cos x) and I4 (α = −1/3, β = −1/4 and f (x) = ex ). Absolute errors Time Absolute errors Time

E [10] (10, 5307, 5000) 7.8 × 10−15 21.539

E [10] (10, 10558, 10000) 6.3 × 10−15 39.459

E [4] (12, 5307, 5000) 2.5 × 10−15 23.984

E [4] (12, 10558, 10000) 8.1 × 10−15 46.637

Remark 4. Numerical examples show that the presented method can produce quite accurate approximations, which in turn verifies that the approximation value of the required modified moments can be of great precision. And Tables 4 and 6 show that forward recursion is numerically stable. It is well known that Mn[c ] can be calculated accurately by forward recursion only if it is a dominant solution of the recurrence relations (see Gautschi [28]). The recurrence relations for c = 6, 7, . . . , 17 are Poincaré difference equations (see Milne-Thomson [38, p. 548]). Then, following the ideas of [14] and references therein, the numerical stability analysis of forward recursion for the recurrence relations for c = 6, 7, . . . , 17 and all α, β, γ > −1 and |a| < 1 will be further discussed and be proved theoretically in our future work. Also, in the future, we will consider some target applications of the oscillatory and singular integrals (1.1), and test concrete numerical examples which occur in practice, e.g., Radon transform, crack problems and Helmholtz equations. 6. Conclusion Owing to a Chebyshev expansion approximating f , numerically stable recurrence relations for the modified moments and properties of Chebyshev polynomials of the first kind, the proposed method for oscillatory integrals with algebraic or logarithmic singularities at the end or interior points of the interval of integration, is constructed by using the expansion (2.4). We consider many different kernel functions G(x), which is a big advantage of the presented approach. The abovementioned numerical experiments show the accuracy of the proposed method. Acknowledgments The authors are deeply grateful to the anonymous referees for their useful comments and helpful suggestions for improvement of this paper. Moreover, they would like to thank editor, Dr. L. Wuytack for his hard and diligent work. The first author also thanks Dr. Xiaojie Wang for his suggestions on the English writing of the manuscript. References [1] H. Brunner, On the numerical solution of first-kind Volterra integral equations with highly oscillatory kernels, Isaac Newton Institute, HOP 13–17, September, 2010: highly oscillatory problems: from theory to applications. [2] H. Brunner, Open problems in the computational solution of Volterra functional equations with highly oscillatory kernels, Isaac Newton Institute, HOP 2007: effective computational methods for highly oscillatory solutions. [3] S.Z. Lu, A class of oscillatory singular integrals, Int. J. Appl. Math. Sci. 2 (2005) 47–64. [4] D.H. Phong, E.M. Stein, Singular integrals related to the Radon transform and boundary value problems, Proc. Nat. Acad. USA 80 (1983) 7697–7701. [5] M.A. Hamed, B. Cummins, A numerical integration formula for the solution of the singular integral equation for classical crack problems in plane and antiplane elasticity, J. King Saud Univ., Eng. Sci. (2) 3 (1991) 217–231. [6] V. Domínguez, I.G. Graham, V.P. Smyshlyaev, Stability and error estimates for Filon–Clenshaw–Curtis rules for highly-oscillatory integrals, IMA J. Numer. Anal. 31 (4) (2011) 1253–1280. [7] V. Domínguez, I.G. Graham, V.P. Smyshlyaev, A hybrid numerical-asymptotic boundary integral method for high-frequency acoustic scattering, Numer. Math. 106 (2007) 471–510. [8] A. Iserles, S.P. Nørsett, Efficient quadrature of highly oscillatory integrals using derivatives, Proc. Royal Soc. A. 461 (2005) 1383–1399. b [9] S. Xiang, Efficient Filon-type methods for a f (x)eiωg (x) dx, Numer. Math. 105 (2007) 633–658. [10] S. Olver, Moment-free numerical integration of highly oscillatory functions, IMA. J. Numer. Anal. 26 (2006) 213–227. [11] D. Huybrechs, S. Vandewalle, On the evaluation of highly oscillatory integrals by analytic continuation, SIAM J. Numer. Anal. 44 (2006) 1026–1048. [12] S. Xiang, X. Chen, H. Wang, Error bounds for approximation in Chebyshev points, Numer. Math. 116 (2010) 463–491. [13] C.W. Clenshaw, A.R. Curtis, A method for numerical integration on an automatic computer, Numer. Math. 2 (1960) 197–205. [14] R. Piessens, M. Branders, The evaluation and application of some modified moments, BIT Numer. Math. 13 (1973) 443–450. [15] T. Hasegawa, T. Torii, Indefinite integration of oscillatory functions by the Chebyshev series expansion, J. Comput. Appl. Math. 17 (1987) 21–29. [16] T. Hasegawa, T. Torii, Application of a modified FFT to product type integration, J. Comput. Appl. Math. 38 (1991) 157–168. [17] P. Keller, A method for indefinite integration of oscillatory and singular functions, Numer. Algorithms 46 (2007) 219–251. [18] P. Keller, P. Wozny, On the convergence of the method for indefinite integration of oscillatory and singular functions, Appl. Math. Comput. 216 (2010) 989–998. [19] H. Wang, S. Xiang, Uniform approximations to Cauchy principal value integrals of oscillatory functions, Appl. Math. Comput. 215 (2009) 1886–1894. [20] P. Keller, A practical algorithm for computing Cauchy principal value integrals of oscillatory functions, Appl. Math. Comput. 218 (2012) 4988–5001. [21] D. Levin, Procedures for computing one-and-two dimensional integrals of functions with rapid irregular oscillations, Math. Comp. 38 (1982) 531–538. [22] R. Chen, S. Xiang, Y. Zhou, A parameter method for computing highly oscillatory integrals, Comput. Math. Appl. 58 (2009) 1830–1837. [23] R. Piessens, M. Branders, On the computation of Fourier transforms of singular functions, J. Comput. Appl. Math. 43 (1992) 159–169. [24] G. Dahlquist, A. Björck, Numerical Methods in Scientific Computing, Vol. I, SIAM, Philadelphia, PA, 2007. [25] J.C. Mason, D.C. Handscomb, Chebyshev Polynomials, Chapman and Hall/CRC, 2003.

156 [26] [27] [28] [29] [30] [31] [32] [33] [34] [35] [36] [37] [38]

H. Kang et al. / Journal of Computational and Applied Mathematics 242 (2013) 141–156 L.N. Trefethen, Is Gauss quadrature better than Clenshaw–Curtis? SIAM Rev. 50 (2008) 67–87. M. Abramowitz, I.A. Stegun, Handbook of Mathematical Functions, National Bureau of Standards, Washington, DC, 1964. W. Gautschi, Computational aspects of three-term recurrence relations, SIAM Rev. 9 (1967) 24–82. E.D. Rainville, Special Functions, The Macmillan Company, New York, 1960. S. Lewanowicz, Construction of a recurrence relation for modified moments, J. Comput. Appl. Math. 5 (1979) 193–206. http://mathworld.wolfram.com/BetaFunction.html. Y.A. Brychkov, Handbook of Special Functions: Derivatives, Integrals, Series and other Formulas, CRC Press, 2008. J. Froehlich, Parameter derivatives of the Jacobi polynomials and Gaussian hypergeometric function, Integral Transforms Spec. Funct. 2 (1994) 252–266. http://functions.wolfram.com/HypergeometricFunctions/Hypergeometric2F1/20/01/. H. Kang, S. Xiang, Diffrentiation formulas of some hypergeometric functions with respect to parameters, Technical Report, 2011, Central South University. A. Gil, N. Temme, J. Segura, Numerical Methods for Special Functions, SIAM, 2008. http://www.mathworks.com/matlabcentral/fileexchange/4540-legendre-gauss-quadrature-weights-and-nodes/content/lgwt.m. L.M. Milne-Thomson, The Calculus of Finite Differences, AMS Chelsea Publishing, 2000.