- Email: [email protected]
On uniform approximations to hypersingular finite-part integrals
JID:YJMAA
AID:19947 /FLA
Doctopic: Applied Mathematics
[m3L; v1.166; Prn:11/11/2015; 10:00] P.1 (1-19)
J. Math. Anal. Appl. ••• (••••) •••–•••
Contents lists available at ScienceDirect
Journal of Mathematical Analysis and Applications www.elsevier.com/locate/jmaa
On uniform approximations to hypersingular finite-part integrals ✩ Shuhuang Xiang a , Chunhua Fang a , Zhenhua Xu b a
Department of Applied Mathematics and Statistics, Central South University, Changsha, Hunan 410083, PR China b College of Mathematics and Information Science, Zhengzhou University of Light Industry, Zhengzhou, Henan 450002, PR China
a r t i c l e
i n f o
Article history: Received 11 October 2014 Available online xxxx Submitted by W.L. Wendland Keywords: Finite-part integral Chebyshev approximation Uniform convergence Error bound Stability
a b s t r a c t In this paper, new uniform approximation schemes for computation of hypersingular finite-part integrals are studied. The methods are verified to be supremely qualified for oscillatory integrands. In addition, based on the results on Chebyshev approximations, the uniform convergence and computational error bounds are considered. Preliminary numerical results show the stability and efficiency of the schemes. © 2015 Published by Elsevier Inc.
1. Introduction Evaluation of the hypersingular finite-part integral ˆ1 Im (f ; c) = = −1
f (x) eiωx dx, (x − c)m+1
c ∈ (−1, 1),
(1.1)
can be transferred into the mth-derivative of the Cauchy principle value integral [14,43] ˆ1 1 dm f (x) iωx e dx, Im (f ; c) = − m! dcm x − c
(1.2)
−1
✩ This work is supported partly by NSF of China (No. 11371376), the Innovation-Driven Project and the Mathematics and Interdisciplinary Sciences Project of Central South University. E-mail addresses: [email protected] (S. Xiang), [email protected] (C. Fang), [email protected] (Z. Xu).
http://dx.doi.org/10.1016/j.jmaa.2015.11.002 0022-247X/© 2015 Published by Elsevier Inc.
JID:YJMAA 2
AID:19947 /FLA
Doctopic: Applied Mathematics
[m3L; v1.166; Prn:11/11/2015; 10:00] P.2 (1-19)
S. Xiang et al. / J. Math. Anal. Appl. ••• (••••) •••–•••
where f and its derivatives f (j) (j ≤ m−1) are assumed to be continuous on [−1, 1], f (m) is Hölder continuous [14,43], i2 = −1 and ω is a parameter. For m = 0, the integral is the usual Cauchy principle value integral or finite Hilbert transform. For m = 1, the integral is the Hadamard finite-part integral [28,43,49,50]. Integrals of the form (1.1) appear frequently in boundary element methods (BEMs) and other numerical computations [3,4,31,43,51], particularly in the area of applied mechanics [31,43]. The efficiency of BEMs often depends upon the efficiency of numerical evaluation of such hypersingular integrals. In the case ω = 0, the computation of integral (1.1) has been extensively studied, such as the Gauss-type method [14,16,17,23,32,33,41,42], the (composite) Newton–Cotes method [21,39,49–51], and others [12,23, 27,37]. A typical procedure for approximating (1.1) for ω = 0 starts by subtracting out the singularity at c as follows [14,16,17,23,41]:
Im (f ; c) =
ˆ1 f (x) − m j=0 −1
(x −
f (j) (c) j! (x m+1 c)
− c)j
ˆ m f (j) (c) dx = dx + , j! (x − c)m+1−j j=0 1
(1.3)
−1
followed by applying some ordinary quadrature rules, such as Gaussian [32,33], Fejér or Clenshaw–Curtis ˆ1 1−c dx = log and rule [9,10,37] to the first integral on the right-hand side in (1.3). While − x−c 1+c −1
ˆ1 = −1
ˆ1 dm−j 1 dx dx , = − (x − c)m+1−j (m − j)! dcm−j x − c
j = 1, . . . , m,
−1
can be computed exactly. Although these schemes are simple in general, there are severe numerical cancellations if one of the nodes happens to be close to c [28,29]. Then Gori and Santi [24,25] used the Gauss–Turán quadrature for the Gori–Michelli weight functions to avoid this problem. Milovanović and Spalević [41] generalized the results from [24,25] by using the Chakalov–Popoviciu quadrature. These developed methods can get high accuracy [41], but they involve computation of the higher derivatives of f at quadrature nodes and are dependent on c. More recently, Sun and Wu [49,50] developed (composite) Newton–Cotes rules for the computation of Hadamard finite-part integrals, and showed the superconvergence. However, these methods are possible to suffer from computational efforts or instability [16,53] as the nodes increase. Efficiently uniform approximation algorithms are proposed by Hasegawa and Torii [28,29] for the finite Hilbert transform and Hadamard finite-part integral ˆ1 f (x) dx, − x−c
ˆ1 =
−1
−1
f (x) dx, (x − c)2
respectively, based on Chebyshev approximations to f (x). Let
pN (x) =
N
N ak Tk (x)
k=0
be the Chebyshev approximate polynomial of f (x) at the Clenshaw–Curtis points cos(jπ/N ) (0 ≤ j ≤ N ), where Tk (x) is the Chebyshev polynomial of the first kind given by Tk (x) = cos(kθ) and x = cos(θ), and
JID:YJMAA
AID:19947 /FLA
Doctopic: Applied Mathematics
[m3L; v1.166; Prn:11/11/2015; 10:00] P.3 (1-19)
S. Xiang et al. / J. Math. Anal. Appl. ••• (••••) •••–•••
3
the double prime denotes a sum whose first and last terms are halved. The coefficients aN k can be computed by FFT [7,13,52]. Here we cite a Matlab code from Trefethen [52]. function I = FFT(f, N)
% (N + 1)-pt Fast Fourier Transform
x = cos(pi ∗ (0 : N) /N);
% Chebyshev points of the first kind
fx = feval(f, x)/(2 ∗ N);
% f evaluated at these points
g = real(fft(fx([1 : N + 1N : −1 : 2])));
% FFT
a = [2g(1); g(2 : N) + g(2 ∗ N : −1 : N + 2); 2g(N + 1)];
% Chebyshev coefficients
For the Hilbert transform in the case ω = 0, Hasegawa and Torii [28] introduced the approximation form ˆ1 QN 0 (f ; c)
1−c pN (x) − pN (c) dx + f (c) log , x−c 1+c
= −1
where (pN (x) − pN (c))/(x − c) can be rewritten in terms of Tk (x) as follows N −1 pN (x) − pN (c) N = dk Tk (x). x−c
(1.4)
k=0
Here the prime denotes the summation whose first term is halved. In addition, the coefficients dN k satisfy a three-term recurrence relation [28] N N N dN k−1 = 2ak + 2cdk − dk+1 , k = N − 1, . . . , 1
(1.5)
N N and the first two initial values are given by dN N = 0 and dN −1 = aN , which yields a new integration formula
N/2−1
QN 0 (f ; c)
=2
j=0
dN 1−c 2j + f (c) log 1 − 4j 2 1+c
(1.6)
with O(N log N ) operations. To extend to the Hadamard finite-part, Hasegawa and Torii [29] constructed a uniform scheme for the Hilbert transform and then applied to the Hadamard finite-part integral, by using the following identity [22] Tk+1 (x) − Tk+1 (c) = 2(x − c)
k
Tj (c)Uk−j (x) = 2(x − c)
j=0
k
Tj (x)Uk−j (c),
k ≥ 1,
(1.7)
j=0
sin((k + 1)θ) (x = cos(θ)) is the second kind of Chebyshev polynomial of order k. A uniform sin(θ) scheme [29] for the Hilbert transform I0 (f ; c) is given as follows where Uk (x) =
ˆ1 QN 0 (f ; c)
= −1
=4
1−c pN (x) − pN (c) dx + f (c) log x−c 1+c
N −1 k=0
AN k Tk (c) + f (c) log
1−c , 1+c
(1.8)
JID:YJMAA
AID:19947 /FLA
Doctopic: Applied Mathematics
[m3L; v1.166; Prn:11/11/2015; 10:00] P.4 (1-19)
S. Xiang et al. / J. Math. Anal. Appl. ••• (••••) •••–•••
4
where (N −k−1)/2
AN k
=
aN 2n+k+1 , 2n + 1
n=0
1 ≤ k ≤ N − 1,
which leads to a uniform scheme [27] for the Hadamard finite-part integral I1 (f ; c) QN 1 (f ; c) =
N −1 2f (c) d N 1−c Q0 (f ; c) = 4 − kAN . k Uk−1 (c) + f (c) log dc 1 + c 1 − c2
(1.9)
k=1
These two schemes (1.8) and (1.9) are stable and convergent independent of c as N trends to infinity, based on the complicated proofs in [27–29]. However they require O(N 2 ) operations. Moreover, the higher order finite-part is not considered in [27–29] and there are few results on this topic. In the case |ω| 1, the integrand is highly oscillatory and the computation of (1.1) will suffer from difficulty by the existed standard algorithms even though there are many algorithms on calculation of highly oscillatory integrals with weak or algebraic singularities or Hilbert transforms [8,11,18–20,26,30,34, 35,44,45,54–59]. In this paper, we shall be concerned with uniform approximation schemes for computation of hypersingular finite-part integral (1.1), particularly for |ω| 1. The overall computational complexity for the schemes is of O(N log N ). Without loss of generality, we assume ω ≥ 0. All the numerical tests in this paper are carried out by using Matlab R2012a. The rest of the paper is organized as follows. In Section 2, we present a uniform approximation scheme for the computation of hypersingular finite-part integrals Im (f ; c). In Section 3 we show the uniform convergence and present error bounds for these new schemes. In Section 4, we illustrate the effectiveness and accuracy by preliminary numerical examples. 2. Uniformly convergent schemes The second integral on the right-hand side in (1.3) with weight eiωx can be evaluated by ˆ m m 1 f (j) (c) eiωx dx = = j! (x − c)m+1−j m! j=0 j=0 1
−1
m j
f
(j)
ˆ1 iωx dm−j e dx (c) m−j − dc x−c −1
⎞
⎛
ˆ1 iωx 1 dm ⎝ e = dx⎠ f (c)− m! dcm x−c −1
=:
1 (m) g (c) m!
and then (1.1) has the form of ˆ1 Im (f ; c) =
e −1
iωx
f (x) −
m j=0
(x −
f (j) (c) j! (x m+1 c)
− c)j
dx +
1 (m) g (c), m!
(2.1)
where g(c) is represented by [36] as ˆ1 g(c) = f (c) −1
eiωx dx = eiωc (Ei(1, iω(1 + c)) − Ei(1, iω(c − 1)) + sgn(ω)πi), x−c
(2.2)
JID:YJMAA
AID:19947 /FLA
Doctopic: Applied Mathematics
[m3L; v1.166; Prn:11/11/2015; 10:00] P.5 (1-19)
S. Xiang et al. / J. Math. Anal. Appl. ••• (••••) •••–•••
and Ei(a, z) is the extended exponential integral Ei(a, z) = (z) > 0 or a > 1 and (z) ≥ 0) and dk Ei(a, iω(1 + z)) = (−iω)k Ei(a − k, iω(1 + z)), dz k
´∞ 1
5
t−a e−zt dt = z a−1 Γ(1 − a, z) (a ≥ 1 and
dk Ei(a, iω(z − 1)) = (−iω)k Ei(a − k, iω(z − 1)) dz k
(see [2]), which together yields that g (m) (c) can be exactly computed. Particularly, for ω = 0, g(c) = 1−c . f (c) log 1+c From (2.1) we define a uniform scheme for (1.1) by ˆ1 QN m (f ; c) =
e
p (x) − iωx N
m (x −
−1
(j)
pN (c) j! (x m+1 c)
j=0
− c)j
dx +
1 (m) g (c). m!
(2.3)
The computation of (2.3) can be done by using the second identity in (1.7) repeatedly following Hasegawa and Torii’s work [28,29]. For example, in the case ω = 0, using the second identity in (1.7), we have the ˆ1 f (x) following uniform approximation for the Hadamard finite-part integral = dx: (x − c)2 −1
ˆ1 QN 1 (f ; c)
= −1
=2
2f (c) 1−c pN (x) − pN (c) − pN (c)(x − c) − dx + f (c) log 2 (x − c) 1 + c 1 − c2
N
aN k
k−1
− pN (c) log N
aN k
k−1
N
aN k
=
+
N
−1
+2
k=3
⎞ 1 − c⎠ Tj (x) − Tj (c) dx + Tj (c) log x−c 1+c
2f (c) 1−c 1−c + f (c) log − 1+c 1 + c 1 − c2
Uk−1−j (c)
j−1
Uj−1−m (c)
m=0
k−1
Uk−1−j (c)Tj (c) log
1 + (−1)m 1−c − pN (c) log 2 1−m 1+c
2f (c) 1−c 1−c + f (c) log − 1+c 1 + c 1 − c2
aN k
k=3 N
ˆ1
j=2
k=3
4aN 2
Uk−1−j (c) ⎝
j=0
k=1
+2
j=0
k=1
=4
⎛
k−1
1−c 1−c Uk−2 (c)(4 + 4c log ) − (k − 1)Uk−1 (c) log 1+c 1+c ⎛
1−c + 2Uj−1 (c) + aN Uk−1−j (c) ⎝Tj (c) log k 1+c j=2
+ f (c) log
2f (c) 1−c − . 1 + c 1 − c2
[(j−1)/2]
=1
⎞ Uk−1−2 (c) ⎠ 1 − 42 (2.4)
m+2 In a similar way, QN ) operations. But this extension does suffer from big m (f ; c) can be achieved in O(N difficulty in computation and become useless for m ≥ 1. However, expression (2.3) is useful in the numerical analysis later.
JID:YJMAA
AID:19947 /FLA
Doctopic: Applied Mathematics
[m3L; v1.166; Prn:11/11/2015; 10:00] P.6 (1-19)
S. Xiang et al. / J. Math. Anal. Appl. ••• (••••) •••–•••
6
It is observed, from (1.2) and (2.1), that ˆ1 e
iωx
f (x) −
m j=0
(x −
−1
f (j) (c) j! (x c)m+1
− c)j
1 dm dx = m! dcm
ˆ1 eiωx −1
f (x) − f (c) dx. x−c
Similarly, we have ˆ1 e
m
p (x) − iωx N
(j)
pN (c) j! (x m+1 c)
j=0
(x −
−1
− c)j
1 dm dx = m! dcm
ˆ1 eiωx −1
pN (x) − pN (c) dx, x−c
which, together with (2.3), establishes a new approximate formula for m ≥ 1 QN m (f ; c)
1 dm = m! dcm
ˆ1 eiωx −1
1 (m) pN (x) − pN (c) dx + g (c). x−c m!
(2.5)
Furthermore, inserting (1.4) into (2.5) gives 1 dm QN m (f ; c) = m! dcm
N −1
N dk Mk (ω)
+
k=0
N −1 1 (m) 1 dm dN 1 (m) g (c) = g (c), Mk (ω) mk + m! m! dc m!
(2.6)
k=0
´1 where Mk (ω) = −1 Tk (x)eiωx dx is called modified moments, which can be efficiently computed by the following recursion [47] together with the Oliver’s algorithm [19,46] iωMk+2 (ω) + 2(k + 2)Mk+1 (ω) − 2iωMk (ω) − 2(k − 2)Mk−1 (ω) + iωMk−2 (ω) = 0 for ω = 0. The Matlab codes can be downloaded from [1]. If ω = 0, Mk (ω) = 0 for k being odd, while 2 Mk (ω) = 1−k 2 for k being even. Based on the above discussion, we outline our algorithm with estimation of computational complexity as follows for m ≥ 1: Algorithm 1. Computation of QN m (f ; c): 1. Compute the coefficients aN k (k = 0, 1, . . . , N ) by using FFT in O(N log N ) operations. 2. Compute Mk (ω) (k = 0, . . . , N ) in O(N ) operations by [1] (ω = 0); or Mk (0) = 0 for k being 2 odd, Mk (0) = 1−k 2 for k being even. dm N 3. Compute dcm dk (k = N, . . . , 0) with the following recurrence relation derived from (1.5) in O(N ) operations dm N dcm dN m
N d dcm dk
= 0,
..., m−1
4. Compute
m
= 0,
N N d d = 2 dc m−1 dk+1 + 2c dcm dk+1 −
1 dm m! dcm
dm N dcm dN −m−1
iωc
f (c)e
dm N dcm dk+2 ,
k = N − m − 2, . . . , 0.
(2.7)
[Ei(1, iω(1 + c)) − Ei(1, iω(c − 1)) + sgn(ω)πi] .
5. Calculate QN m (f, c) (2.6) in O(N ) operations. From Steps 1–5, we can see that this new scheme can be achieved in O(N log N ) operations for any m. Moreover, the new scheme has the N -order algebraic accuracy if we do not consider the rounding errors. However, this scheme is not stable when N or m increases. We will see that in Section 4.
JID:YJMAA
AID:19947 /FLA
Doctopic: Applied Mathematics
[m3L; v1.166; Prn:11/11/2015; 10:00] P.7 (1-19)
S. Xiang et al. / J. Math. Anal. Appl. ••• (••••) •••–•••
7
Notice that the linear system (2.7) can be represented by linear system Ax = b as ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝
⎞⎛
1 1 ..
.
1 1 −2c 1 1 −2c 1 1 −2c 1 .. .. .. . . . 1 −2c 1
⎟⎜ ⎟⎜ ⎟⎜ ⎟⎜ ⎟⎜ ⎟⎜ ⎟⎜ ⎟⎜ ⎟⎜ ⎟⎜ ⎟⎜ ⎟⎜ ⎟⎜ ⎟⎜ ⎟⎜ ⎟⎜ ⎠⎝
⎛
⎞ 0 ⎟ ⎟ ⎜ 0 ⎟ ⎟ ⎜ ⎜ ⎟ ⎟ ⎜ .. .. ⎟ ⎟ ⎜ . ⎟ . ⎟ ⎜ ⎟ ⎟ ⎜ dm N ⎟ 0 ⎟ d m ⎜ ⎟ N −m−1 dc ⎟ ⎜ dm−1 N ⎟ dm N ⎟ 2 d m−1 = d ⎜ ⎟. N −m−1 dc dcm N −m−2 ⎟ ⎜ ⎟ m−1 m ⎟ ⎜2 d N d N ⎟ dcm−1 dN −m−2 ⎟ dcm dN −m−3 ⎟ ⎜ ⎟ ⎜ dm−1 N dm N ⎟ ⎜ 2 dcm−1 dN −m−3 ⎟ ⎟ dcm dN −m−4 ⎟ ⎟ ⎟ ⎜ .. .. ⎜ ⎟ ⎠ ⎝ . . ⎠ m m−1 d N d N 2 dcm−1 d1 dcm d0 dm N dcm dN dm N dcm dN −1
⎞
(2.8)
The coefficient matrix is well-conditioned and then the system (2.7) can be efficiently computed by GMRES method based on (2.8) for m ≥ 1. Algorithm 2. Computation of QN m (f ; c): 1. Compute the coefficients aN k (k = 0, 1, . . . , N ) by using FFT with O(N log N ) operations. 2. Compute Mk (ω) (k = 0, . . . , N ) in O(N ) operations by [1] (ω = 0); or Mk (0) = 0 for k being 2 odd, Mk (0) = 1−k 2 for k being even. dm N 3. Compute the coefficients dcm dk (k = N, . . . , 0) by using the GMRES
for (2.8). 4. Compute
1 dm m! dcm
f (c)eiωc [Ei(1, iω(1 + c)) − Ei(1, iω(c − 1)) + sgn(ω)πi] .
5. Calculate QN m (f ; c) in O(N ) operations. Then Scheme (2.6) is uniformly convergent independent of c and ω as N trends to infinity which can be seen in Section 3. Furthermore, the accuracy increases when oscillation becomes faster. Remark 1. In particular, for ω = 0,
1−c 1 dm 1 dm iωc f (c)e f (c) log . [Ei(1, iω(1 + c)) − Ei(1, iω(c − 1)) + sgn(ω)πi] = m! dcm m! dcm 1+c Remark 2. To obtain the higher order for the convergence rate with respect to ω, here we also present a new modified scheme for (2.6) N (f ; c) = QN (f ; c) + 1 Q m m iω
−Pm (1) iω Pm (−1) −iω , e + e (1 − c)m+1 (−1 − c)m+1
(2.9)
m f (j) (c)−p(j) N (c) where Pm (x) = j=0 (x − c)j . In Section 3, we will prove the modified scheme’s accuracy is j! −2 O(ω ). In Section 4, numerical examples are presented to illustrate the effectiveness. Remark 3. For small values of N , Algorithm 1 and Algorithm 2 have nearly same accuracy (see Figs. 2, 5 and Tables 1–3 in Section 4). 3. Uniform convergence and computational error bounds In this section, we consider the uniform convergence and error bound of Scheme (2.6) independent of c and ω. Here we assume that dN k can be computed accurately, which can be done by Algorithm 1 for small values of N and m, or by Algorithm 2 based on the GMRES.
JID:YJMAA
AID:19947 /FLA
Doctopic: Applied Mathematics
[m3L; v1.166; Prn:11/11/2015; 10:00] P.8 (1-19)
S. Xiang et al. / J. Math. Anal. Appl. ••• (••••) •••–•••
8
From (2.1) and (2.2), it follows that ˆ1 Im (f ; c) −
QN m (f ; c)
=
e
iωx
F (x) −
m j=0
(x −
−1
F (j) (c) (x j! m+1 c)
− c)j
dx,
(3.1)
where F (x) = f (x) − pN (x). Thus, the error bound of Scheme (2.6) can be estimated by the remainder of the Taylor’s expansion of F (x) as Im (f ; c) − QN m (f ; c) ≤
2 (m+1) (m+1) − pN f . (m + 1)! ∞
(3.2)
In the following, we will focus on the estimate of the right-hand side of (3.2). Assumption 1. Suppose that f (x) has an absolutely continuous (k0 − 1)st derivative f (k0 −1) on [−1, 1] for ´x some k0 ≥ max{1, m} with f (k0 −1) (x) = f (k0 −1) (−1) + −1 h(y)dy, where h is absolutely integrable and of bounded variation Var(h) < ∞ on [−1, 1], and define Vk0 = inf
f (k0 −1) (x) = f (k0 −1) (−1) + ´ x h(y)dy for all x ∈ [−1, 1] with h being −1 Var(h) . absolutely integrable and of bounded variation
Let {bn } denote the exact spectral coefficients of f (x). Then the Chebyshev series for f is defined as [53]
f (x) =
∞
bj Tj (x),
bj =
j=0
2 π
ˆ1 −1
f (x)Tj (x) √ dx. 1 − x2
Lemma 3.1. (i) (Bernstein [5]) If f is analytic with |f (z)| ≤ M in the region Eρ bounded by the ellipse with foci {−1, 1} and major and minor semiaxis lengths summing to ρ > 1, then for each j ≥ 0, |bj | ≤
2M . ρj
(3.3)
(ii) (Trefethen [53]) Suppose that f (x) satisfies Assumption 1 for some k0 ≥ max{1, m}, then for each j ≥ k + 1, |bj | ≤
2Vk0 . πj(j − 1) · · · (j − k0 )
(3.4)
(iii) (Boyd [6, (4.55), (4.57)]) The coefficients of the interpolating polynomial pN (x) are related to those of f (x) via aN n = bn +
∞ (bn+2jN + b−n+2jN ),
n = 0, 1, . . . , N,
|f (x) − pN (x)| ≤ 2
j=1
∞
|bj |.
(3.5)
j=N +1
Moreover, the error of the interpolation polynomial pN can be estimated by the following formula [6,15] f (x) − pN (x) =
N ∞ b0 − aN (bj − aN bj Tj (x). 0 + j )Tj (x) + 2 j=1 j=N +1
JID:YJMAA
AID:19947 /FLA
Doctopic: Applied Mathematics
[m3L; v1.166; Prn:11/11/2015; 10:00] P.9 (1-19)
S. Xiang et al. / J. Math. Anal. Appl. ••• (••••) •••–•••
9
Then it implies that (m+1)
|f (m+1) (x) − pN
N
(x)| ≤
(m+1)
|bj − aN j | Tj
∞
∞ +
j=m+1
(m+1)
|bj | Tj
∞ .
(3.6)
j=N +1
In addition, for each n ≥ p, it follows from (A.12), (A.13) and (A.52) in [6, pp. 497, 503] that p−1 n2 − k2 dp n+p T (x)| = (±1) n ±1 p dx 2k + 1
dp (p) Tn (x) = n2p−1 (p − 1)!Cn−p (x), dxp
k=0
and (p)
(p)
|Cn−p (x)| ≤ |Cn−p (1)|,
x ∈ [−1, 1],
(p)
where Cn−p (x) is the Gegenbauer polynomial [6]. Then it establishes that p d dxp Tn
= ∞
p−1 k=0
n2p n2 − k2 ≤ . 2k + 1 (2p − 1)!!
(3.7)
This, together with (3.3), Lemma 3.1 and (4.55) in [6, p. 96], gives N ∞ j 2(m+1) j 2(m+1) (m+1) (m+1) bj − aN + (x) − pN (x) ≤ |b | f j j (2m + 1)!! (2m + 1)!! j=m+1 j=N +1
≤ N 2(m+1)
N j=0
∞
≤ N 2(m+1) ∞
|bj |j 2(m+1)
j=N +1
|bj | +
j=N +1
≤2
∞
|bj − aN j |+
∞
|bj |j 2(m+1) ([6, (4.55)])
j=N +1
|bj |j 2(m+1)
j=N +1
≤
∞ j=N +1
4Vk0 ≤ π
=
4Vk0 j 2m+2 ([6, Lemma 3.1]) πj(j − 1) · · · (j − 2m − 1)(j − k0 )k0 −2m−1
N +1 N − 2m
N +1 N − 2m
2m+1 ˆ∞
2m+1
N
1 dx (x − k0 )k0 −2m−1
2Vk0 , (k0 − 2m)π(N − k0 )k0 −2m
where we used the following inequalities for the proof of the last inequality in (3.8) j j j < < ··· < , j−1 j−2 N − 2m
j = N + 1, N + 2, . . . .
From the equations (3.2) and (3.8) we can deduce the following theorem directly.
(3.8)
JID:YJMAA
AID:19947 /FLA
Doctopic: Applied Mathematics
[m3L; v1.166; Prn:11/11/2015; 10:00] P.10 (1-19)
S. Xiang et al. / J. Math. Anal. Appl. ••• (••••) •••–•••
10
Theorem 3.1. Suppose that f (x) satisfies Assumption 1 for some k0 ≥ 2m + 1. Then Scheme (2.6) is uniformly convergent independent of c and satisfies that Im (f ; c) − QN m (f ; c) ≤
N +1 N − 2m
2m+1
4Vk0 . (m + 1)!(k0 − 2m)π(N − k0 )k0 −2m
(3.9)
In particular, if f is analytic with |f (z)| ≤ M in the region Eρ bounded by the ellipse with foci {−1, 1} and major and minor semiaxis lengths summing to ρ > 1, then for each N ≥ 0, based on the line integral for the error of the Chebyshev approximation polynomial pN [15, p. 390], an estimation is given by f (x) − pN (x) =
1 2πi
˛ Eρ
ωN +1 (x)f (z) dz, (z − x)ωN +1 (z)
where ωN +1 (x) = TN +1 (x) − TN −1 (x). By a similar proof to the Cauchy’s differentiation formula [38], (m+1) f (m+1) (x) − pN (x) can be represented by (m+1)
f (m+1) (x) − pN
(m+1)
Therefore, f (m+1) (x) − pN
(x) =
1 2πi
˛ Eρ
f (z) ωN +1 (z)
ωN +1 (x) z−x
(m+1) dz.
(x) can be estimated as follows
(m+1) (m+1) − pN f
∞
(m+1) 1 ˛ ωN +1 (x) f (z) = dz 2πi ωN +1 (z) z−x Eρ
.
(3.10)
∞
From (3.7), we see that p d ω dxp N +1
∞
p 2p d ≤ 2(N + 1) . = (T − T ) N −1 dxp N +1 (2p − 1)!! ∞
(3.11)
Moreover, for j = 0, 1, . . . , m + 1, we have (N + 1)2j+2 (N + 1)2(m+1) (N + 1)2j ≤ ≤ ··· ≤ , j!(2j − 1)!! (j + 1)!(2j + 1)!! (m + 1)!(2m + 1)!! where N ≥
√
2m2 + 3m + 1 − 1. Then for each fixed z ∈ Eρ we get
(m+1) ω N +1 (x) z−x
∞
m+1 m + 1 (m − j + 1)! (j) = ωN +1 (x) m−j+2 (z − x) j j=0 ≤ 2(m + 1)!
m+1 j=0
∞
2j
(N + 1) j!(2j − 1)!!ρm+2−j
≤
m+1 2(N + 1)2(m+1) 1 m+2−j (2m + 1)!! ρ j=0
≤
2(N + 1)2(m+1) . (2m + 1)!!(ρ − 1)
(3.12)
JID:YJMAA
AID:19947 /FLA
Doctopic: Applied Mathematics
[m3L; v1.166; Prn:11/11/2015; 10:00] P.11 (1-19)
S. Xiang et al. / J. Math. Anal. Appl. ••• (••••) •••–•••
11
Additionally, from (1.50) in [40] it follows that for z ∈ Eρ 1 N +1 (ρ − ρ−N −1 − ρN −1 + ρ−N +1 ) ≤ |ωN +1 (z)| = |TN +1 (z) − TN −1 (z)| , 2
(3.13)
which, together with (3.10), (3.12) and (3.13), derives that (m+1) (m+1) − pN f
∞
(m+1) 1 ˛ ωN +1 (x) f (z) = dz 2πi ωN +1 (z) z−x Eρ ≤
−1
(2m +
∞
2(m+1)
M (ρ + ρ )(N + 1) . − ρ−N −1 − ρN −1 + ρ−N +1 )(ρ − 1)
1)!!(ρN +1
Theorem 3.2. If f is analytic with |f (z)| ≤ M in the region Eρ bounded by the ellipse with foci {−1, 1} and √ major and minor semiaxis lengths summing to ρ > 1, then for N ≥ 2m2 + 3m + 1 − 1, the scheme (2.6) is uniformly convergent independent of c and satisfies that Im (f ; c) − QN m (f ; c) ≤
2M (ρ + ρ−1 )(N + 1)2(m+1) . (2m + 1)!!(m + 1)!(ρN +1 − ρ−N −1 − ρN −1 + ρ−N +1 )(ρ − 1)
(3.14)
It is worth noting that the accuracy of Scheme (2.6) increases when oscillation becomes faster. Theorem 3.3. Let f be a smooth function on [−1, 1]. Then
−Pm (1) iω 1 Pm (−1) N −iω e + e I(f ; c) − Qm (f ; c) ≤ iω (1 − c)m+1 (−1 − c)m+1 (m+2) (m+2) (m+3) − pN + 2 f (m+3) − pN f ∞ ∞ , + ω 2 (m + 1)! where Pm (x) =
m j=0
(j)
f (j) (c)−pN (c) (x j!
(3.15)
− c)j .
Before proving Theorem 3.3, we introduce the following lemmas. Lemma 3.2. (See van der Corput [48].) Suppose g(x) is real-valued and smooth in (a, b) and that |g (k) (x)| ≥ 1 for all x ∈ (a, b) for a fixed value of k. Then b ˆ eiωg(x) dx ≤ c(k)ω −1/k a
holds when (i) k ≥ 2, or (ii) k = 1 and g (x) is monotonic. The bound c(k) = 5 · 2k−1 − 2 is independent of g and ω. Lemma 3.3. (See [48, p. 334].) Under the assumptions on g(x) in Lemma 3.2, we can conclude that b ⎤ ⎡ ˆ ˆb eiωg(x) ϕ(x)dx ≤ c(k)ω −1/k ⎣|ϕ(b)| + |ϕ (x)|dx⎦ . a
a
JID:YJMAA
AID:19947 /FLA
Doctopic: Applied Mathematics
[m3L; v1.166; Prn:11/11/2015; 10:00] P.12 (1-19)
S. Xiang et al. / J. Math. Anal. Appl. ••• (••••) •••–•••
12
Proof of Theorem 3.3. Set Pm (x) =
m F (j) (c) j=0
j!
(x − c)j to be the Taylor polynomial of F (x) = f (x) −pN (x)
with degree deg(Pm ) ≤ m. Integrating by parts and from Lemma 3.3, we have ˆ1 F (x) − Pm (x) iωx 1 Im (f ; c) − QN = (f ; c) de m iω (x − c)m+1 −1 1 −Pm (1) iω Pm (−1) −iω ≤ ( e + e ) iω (1 − c)m+1 (−1 − c)m+1 ˆ1
1 F (x) − P (x) m iωx + e dx m+1 iω (x − c) −1 1 −Pm (1) iω Pm (−1) −iω ≤ ( e + e ) m+1 m+1 iω (1 − c) (−1 − c)
F (x) − Pm (x) 3 F (x) − P (x) m + 2 + 2 , ω (x − c)m+1 (x − c)m+1 ∞ ∞
(3.16)
F (m+1) (c) F (x) − Pm (x) F (x) − Pm (x) to define the value of = at x = c. We m+1 x→c (x − c) (m + 1)! (x − c)m+1
where we used the limit lim have
F (x) − Pm (x) (x − c)m+1
⎧ (x))(x − c) − (m + 1)(F (x) − Pm (x)) (F (x) − Pm ⎪ ⎪ , x = c ⎨ (x − c)m+2 = ⎪ F (m+2) (c) ⎪ ⎩ , x=c (m + 2)!
and
F (x) − Pm (x) (x − c)m+1
⎧ ⎨
=
⎩
F (x)−Pm (x))(x−c)2 −2(m+1)(F (x)−Pm (x))(x−c)+(m+1)(m+2)(F (x)−Pm (x)) , (x−c)m+3 (m+3)
2F (c) (m+3)! ,
x = c x = c.
Define G(x) = (F (x) − Pm (x))(x − c) − (m + 1)(F (x) − Pm (x)), it is easy to verify that
G(c) = G (c) = · · · = G(m) (c) = 0 and (m+2) (x))(x − c) = F (m+2) (x)(x − c). G(m+1) (x) = (F (m+2) (x) − Pm
Thus, the Taylor’s expansion of G(m+1) (x) at x = c has the form G(x) =
F (m+2) (ζ)(ζ − c) (x − c)m+1 , (m + 1)!
ζ = c + θ(x − c) (0 < θ < 1),
and for all x ∈ [−1, 1] ⎧ G(x) F (m+2) ∞ |ζ − c| ⎪ ⎪
F (x) − P (x) ⎪ ⎨ (x − c)m+2 ≤ (m + 1)! |x − c| , x = c F (m+2) ∞ m ≤ . = m+1 (m+2) ⎪ (x − c) (m + 1)! ∞ ⎪ F ⎪ , x=c ⎩ (m + 2)!
(3.17)
JID:YJMAA
AID:19947 /FLA
Doctopic: Applied Mathematics
[m3L; v1.166; Prn:11/11/2015; 10:00] P.13 (1-19)
S. Xiang et al. / J. Math. Anal. Appl. ••• (••••) •••–•••
13
Table 1 Absolute errors in 2n + 1-points approximation on Scheme (2.6) based on (2.7) and (2.8) to the Hadamard finite-part integral ˆ1 x6 = dx. (x − c)2 0
1 64
n
c=
c=
2 3 4
0.4336285 7.438494e-15 7.771561e-15
1 16
0.2110565 1.137979e-15 2.803313e-15
c = 0.5
c = 0.99999
0.0125 MP MP
0.7311613 MP MP
Similarly, let G(x) = (F (x) −Pm (x))(x −c)2 −2(m +1)(F (x) −Pm (x))(x −c) +(m +1)(m +2)(F (x) −Pm (x)), there holds
F (x) − Pm (x) (x − c)m+1
=
G(x) (c) = · · · = G (m) (c) = 0, G (m+1) (x) = F (m+3) (x)(x − c)2 , G(c) =G (x − c)m+3
and
F (x) − P (x) F (m+3) m ∞ . ≤ (x − c)m+1 (m + 1)!
(3.18)
A combination of (3.16)–(3.18) establishes the desired result. 2 From the proof of Theorem 3.3, we can deduce the following corollary immediately. Corollary 3.1. Let f be a smooth function on [−1, 1]. Then (m+2) (m+3) (m+2) − pN + 2 f (m+3) − pN f ∞ ∞ N (f ; c) ≤ . I(f ; c) − Q m 2 ω (m + 1)!
(3.19)
4. Numerical examples In this section, all the experiments are performed on the computer with 2.40 GHz laptop personal computer with 2 GB of RAM. We are using the R2012a version of the Matlab system. As stated in the following Tables 1–3, “MP” indicates that all the data have reached machine precision in Matlab. Of particular interest is the fact that for small values of N , Scheme (2.6) based on (2.7) and (2.8) have nearly same accuracy (see Tables 1–3 and Figs. 2 and 5). However, for large values of N , the difference becomes more significant as m increases. Example 1. We consider Scheme (2.6) based on (2.7) and (2.8) respectively for the finite-part integrals [49] ˆ1 = 0
ˆ1 1 x6 (t + 1)6 = dx = dt, (x − c)2 32 (t − 2c + 1)2
c ∈ (0, 1)
−1
(see Fig. 1 and Table 1). The exact solution [49] is 1 1−c 6 3 + c + 2c2 + 3c3 + 6c4 + + 6c5 log . 5 2 c−1 c Here, the GMRES for (2.8) converges at iteration 5 to a solution with relative residual less than 1e − 13.
JID:YJMAA
AID:19947 /FLA
Doctopic: Applied Mathematics
[m3L; v1.166; Prn:11/11/2015; 10:00] P.14 (1-19)
S. Xiang et al. / J. Math. Anal. Appl. ••• (••••) •••–•••
14
Fig. 1. The relative errors of Scheme (2.6) based on (2.7) and (2.8) for
Fig. 2. The relative errors of Scheme (2.6) based on (2.7) and (2.8) for respectively.
´ x6 =1 0 (x−c)2
dx with c = 0.549 and c = 0.99999, respectively.
´1 1−a2 = −1 (1−2ax+a2 )(x−c)3
dx with c = 0.549 and c = 0.99999,
Example 2. We illustrate the accuracy of Scheme (2.6) for computing the finite-part integrals [29] ˆ1 = −1
1 − a2 dx, (1 − 2ax + a2 )(x − c)3
ˆ1 = −1
1 − a2 dx, (1 − 2ax + a2 )(x − c)4
a = 0.8,
c ∈ (−1, 1)
1 dI1 (a, c) 1 d2 I1 (a, c) and , respectively, where 2 dc 6 dc2
1−c 1−a 2(1 − a2 ) 2a(1 − a2 ) log − 2 log − . I1 (a; c) = (1 − 2ac + a2 )2 1+c 1+a (1 − 2ac + a2 )(1 − c2 )
(see Figs. 2–4). The exact solutions are
Here, the GMRES for (2.8) converges at iteration 1 to a solution with relative residual less than 1e − 10 for N ≥ 200. Example 3. We demonstrate Scheme (2.6) for computing the Cauchy principal value integral I0 (f, c) = ´ 1 x sin(ωx) − e x−c dx with c = 0 (see Table 2). From maple 14, we see that the integral I0 (f, c) can be explicitly −1 expressed by the following formula
JID:YJMAA
AID:19947 /FLA
Doctopic: Applied Mathematics
[m3L; v1.166; Prn:11/11/2015; 10:00] P.15 (1-19)
S. Xiang et al. / J. Math. Anal. Appl. ••• (••••) •••–•••
Fig. 3. The relative errors of Scheme (2.6) based on (2.7) and (2.8) for respectively.
Fig. 4. The relative errors of Scheme (2.6) based on (2.7) and (2.8) for respectively.
15
´ 1−a2 =1 −1 (1−2ax+a2 )(x−c)3
dx with c = 0.549 and c = 0.99999,
´ 1−a2 =1 −1 (1−2ax+a2 )(x−c)4
dx with c = 0.549 and c = 0.99999,
Table 2
ˆ1 x sin(ωx) Absolute errors in 2n + 1-point approximation to the finite-part integral − e dx by Scheme (2.6). x −1
n
ω = 10
ω = 102
ω = 104
ω = 106
1 2 3 4
1.396989e-3 3.163167e-5 −1.418378e-9 MP
8.533057e-6 1.1859500e-7 7.234213e-13 MP
5.447682e-10 8.905765e-12 MP MP
6.217249e-14 8.881784e-16 MP MP
I0 (f, c) =
where E1 (z) =
´∞ z
i c(1+iω) e E1 (1 − c)(1 + iω) − E1 − (1 + c)(1 + iω) + ln(−1 − iω) 2 ! − ln(1 + iω) + e−c(−1+iω) E1 (1 + c)(−1 + iω) − E1 (1 − c)(1 − iω) ! + ln(1 − iω) − ln(−1 + iω) , −1 < c < 1,
e−t t dt
is the exponential integral.
(4.1)
JID:YJMAA
AID:19947 /FLA
16
Doctopic: Applied Mathematics
[m3L; v1.166; Prn:11/11/2015; 10:00] P.16 (1-19)
S. Xiang et al. / J. Math. Anal. Appl. ••• (••••) •••–•••
Fig. 5. The relative errors scaled by ω of Scheme (2.9) based on (2.7) and (2.8), respectively, for the computation of with ω = 108 .
Fig. 6. The absolute errors scaled by ω of Scheme (2.6) based on (2.7) for the computation of k = 1 : 2000 and N = 9.
´ =1 −1
´ =1 −1
ex sin(ωx) dx x2
ex sin(ωx) dx with ω = 1 + 10k, x2
´1 We also consider Scheme (2.6) and (2.9) for computing the finite-part integral I1 (f, 0) = =−1 ex sin(ωx) x2 dx d (see Figs. 5–7 and Table 3). The exact solution is dc I0 (f, c) c=0 , where I0 (f, c) is denoted by (4.1). Here, the GMRES for (2.8) converges at iteration 3 to a solution with relative residual less than 1e − 10 for N ≥ 200. As can be seen from Figs. 6–7, there is a good agreement between the theory results and the computations. 5. Conclusions We have presented efficient methods for the computation of hypersingular finite-part integrals. The convergence analysis for the corresponding methods is also given. The new uniform scheme allows calculating the Hadamard finite-part integral in O(N logN ) operations, while the existing algorithms require O(N 2 ) operations. Moreover, this method fills the gaps in fast computation of hypersingular integrals, especially for oscillatory integrands.
JID:YJMAA
AID:19947 /FLA
Doctopic: Applied Mathematics
[m3L; v1.166; Prn:11/11/2015; 10:00] P.17 (1-19)
S. Xiang et al. / J. Math. Anal. Appl. ••• (••••) •••–•••
Fig. 7. The absolute errors scaled by ω 2 of Scheme (2.9) based on (2.7) for computation of k = 1 : 2000 and N = 9.
´ =1 −1
17
ex sin(ωx) dx with ω = 1 + 10k, x2
Table 3
ˆ1 x sin(ωx) dx with Scheme (2.9). The absolute errors of 2n + 1-point approximation values to the finite-part integral = e x2 −1
n
ω = 102
ω = 104
ω = 106
ω = 108
1 2 3 4
1.458692e-6 8.270350e-9 7.00000e-15 MP
9.380800e-11 7.130000e-13 1.000000e-15 MP
1.000000e-14 MP MP MP
MP MP MP MP
Acknowledgments The authors are grateful to the anonymous referee for his constructive comments and helpful suggestions for improvement of this paper greatly. The authors also thank Dr. Guo He and Junjie Ma for their checking up every detail. References [1] http://www.unavarra.es/personal/victor_dominguez/clenshawcurtisrule/index.html. [2] M. Abramowitz, I.A. Stegun, Handbook of Mathematical Functions, National Bureau of Standards, Washington, DC, 1964. [3] M. Ainsworth, B. Guo, An additive Schwarz preconditioner for p-version boundary element approximation of the hypersingular operator in three dimensions, Numer. Math. 85 (2000) 343–366. [4] G. Bao, W. Sun, A fast algorithm for the electromagnetic scattering from a large cavity, SIAM J. Sci. Comput. 27 (2005) 553–574. [5] S.N. Bernstein, Sur l’ordre de la meilleure approximation des fonctions continues par les polynômes de degré donné, Mem. Cl. Sci. Acad. Roy. Belg. 4 (1912) 1–103. [6] J.H. Boyd, Chebyshev and Fourier Spectral Methods, Dover Publications, Inc., New York, 2000. [7] M. Branders, R. Piessens, An extension of Clenshaw–Curtis quadrature, J. Comput. Appl. Math. 1 (1975) 55–65. [8] M.R. Capobiancoa, G. Criscuoloa, On quadrature for Cauchy principal value integrals of oscillatory functions, J. Comput. Appl. Math. 156 (2003) 471–486. [9] M.M. Chawla, N. Jayarajan, Quadrature formulas for Cauchy principle value integrals, Computing 15 (1975) 347–355. [10] M.M. Chawla, S. Kumar, Convergence of quadrature for Cauchy principle value integrals, Computing 23 (1979) 67–72. [11] R. Chen, C. An, On evaluation of Bessel transforms with oscillatory and algebraic singular integrands, J. Comput. Appl. Math. 264 (2014) 71–81. [12] U.J. Choi, S.W. Kim, B.I. Yun, Improvement of the asymptotic behaviour of the Euler–Maclaurin formula for Cauchy principle value and Hadamard finite-part integrals, J. Comput. Appl. Math. 61 (2004) 496–513. [13] C.W. Clenshaw, A.R. Curtis, A method for numerocal integration on an automatic computer, Numer. Math. 2 (1960) 197–205.
JID:YJMAA 18
AID:19947 /FLA
Doctopic: Applied Mathematics
[m3L; v1.166; Prn:11/11/2015; 10:00] P.18 (1-19)
S. Xiang et al. / J. Math. Anal. Appl. ••• (••••) •••–•••
[14] G. Criscuolo, A new algorithm for Cauchy principle value and Hadamard finite-part integrals, J. Comput. Appl. Math. 78 (1997) 255–275. [15] G. Dahlquist, A. Björck, Numerical Methods in Scientific Computing, SIAM, Philadelphia, PA, 2007. [16] P.J. Davis, P. Rabinowitz, Methods of Numerical Integration, second edition, Academic Press, 1984. [17] K. Diethelm, A method for the practical evaluation of the Hilbert transform on the real line, J. Comput. Appl. Math. 112 (1999) 45–53. [18] V. Domínguez, Filon–Clenshaw–Curtis rules for a class of highly-oscillatory integrals with logarithmic singularities, J. Comput. Appl. Math. 261 (2014) 299–319. [19] V. Domínguez, I.G. Graham, V.P. Smyshlyaev, Stability and error estimates for Filon–Clenshaw–Curtis rules for highlyoscillatory integrals, IMA J. Numer. Anal. 31 (2011) 1253–1280. [20] V. Domínguez, I.G. Graham, M. Kim, Filon–Clenshaw–Curtis rules for highly-oscillatory integrals with algebraic singularities and stationary points, SIAM J. Numer. Anal. 51 (2013) 1542–1566. [21] Q. Du, Evaluations of certain hypersingular integrals on interval, Internat. J. Numer. Methods Engrg. 51 (2001) 1195–1210. [22] D. Elliott, Truncation errors in two Chebyshev series approximations, Math. Comp. 19 (1965) 234–248, 445–462. [23] D. Elliott, E. Venturino, Sigmoidal transformations and the Euler–Maclaurin expansion for evaluating Hadamard finite-part integrals, Numer. Math. 77 (1997) 453–465. [24] L. Gori, E. Santi, On the evaluation of Hilbert transforms by means of a particular class of Turán quadrature rules, Numer. Algorithms 10 (1995) 27–39. [25] L. Gori, E. Santi, Quadrature rules based on s-orthogonal polynomials for evaluating integrals with strong singularities, in: W. Gautschi, G.H. Golub, G. Opfer (Eds.), Applications and Computation of Orthogonal Polynomials, Oberwolfach, 1998, in: Internat. Ser. Numer. Math., vol. 131, Birkhäuser, Basel, 1999, pp. 109–119. [26] A.I. Hascelik, Suitable Gauss and Filon-type methods for oscillatory integrals with an algebraic singularity, Appl. Numer. Math. 59 (2009) 101–118. [27] T. Hasegawa, Uniform approximations to finite Hilbert transform and its derivatives, J. Comput. Appl. Math. 163 (2004) 127–138. [28] T. Hasegawa, T. Torii, An automatic quadrature for Cauchy principle value integrals, Math. Comp. 56 (1991) 741–754. [29] T. Hasegawa, T. Torii, Hilbert and Hadamard transforms by generalized Chebyshev expansion, J. Comput. Appl. Math. 51 (1994) 71–83. [30] G. He, S. Xiang, An improved algorithm for the evaluation of Cauchy principal value integrals of oscillatory functions and its application, J. Comput. Appl. Math. 280 (2015) 1–13. [31] C. Hui, D. Shia, Evaluations of hypersingular integrals using Gaussian quadrature, Internat. J. Numer. Methods Engrg. 44 (1999) 205–214. [32] N.I. Ioakimidid, Application of finite-part integrals to the singular integral equation of crack problems in plane and three-dimensional elasticity, Acta Mech. 45 (1982) 31–47. [33] N.I. Ioakimidid, On the uniform convergence of Gaussian quadrature rules for Cauchy principle value integrals and their derivatives, Math. Comp. 44 (1985) 191–198. [34] H. Kang, S. Xiang, On the calculation of highly oscillatory integrals with an algebraic singularity, Appl. Math. Comput. 217 (2010) 3890–3897. [35] H. Kang, S. Xiang, Efficient quadrature of highly oscillatory integrals with algebraic singularities, J. Comput. Appl. Math. 237 (2013) 576–588. [36] P. Keller, A practical algorithm for computing Cauchy principal value integrals of oscillatory functions, Appl. Math. Comput. 218 (2012) 4988–5001. [37] P. Kim, S. Ahn, U. Choi, On the Chebyshev quadrature rule of finite part integrals, J. Comput. Appl. Math. 180 (2005) 147–159. [38] A. Lars, Complex Analysis: An Introduction to the Theory of Analytic Functions of One Complex Variable, McGraw–Hill, New York, 1979. [39] P. Linz, On the approximate computation of certain strongly singular integrals, Computing 35 (1985) 345–353. [40] J.C. Mason, D.C. Handscomb, Chebyshev Polynomials, Chapman and Hall/CRC, New York, 2003. [41] G.V. Milovanović, M. Spalević, Quadrature rules with multiple nodes for evaluating integrals with strong singularities, J. Comput. Appl. Math. 189 (2006) 689–702. [42] G. Monegato, Numerical evaluation of hypersingular integrals, J. Comput. Appl. Math. 50 (1997) 9–31. [43] G. Monegato, Definitions, properties and applications of finite-part integrals, J. Comput. Appl. Math. 229 (2009) 425–439. [44] G.E. Okecha, Quadrature formulae for Cauchy principal value integrals of oscillatory kind, Math. Comp. 49 (1987) 259–268. [45] G.E. Okecha, Hermite interpolation and a method for evaluating cauchy principal value integrals of oscillatory kind, Kragujevac J. Math. 29 (2006) 91–98. [46] J. Oliver, The numerical solution of linear recurrence relations, Numer. Math. 11 (1968) 349–360. [47] R. Piessens, M. Branders, On the computation of Fourier transforms of singular functions, J. Comput. Appl. Math. 43 (1992) 159–169. [48] E. Stein, Harmonic Analysis: Real-Variable Methods, Orthogonality, and Oscillatory Integrals, Princeton University Press, Princeton, NJ, 1993. [49] W. Sun, J. Wu, Interpolatory quadrature rules for Hadamard finite-part integrals and their supperconvergence, IMA J. Numer. Anal. 28 (2008) 580–597. [50] W. Sun, J. Wu, The supperconvergence of Newton–Cotes rules for the Hadamard finite-part integral on an interval, Numer. Math. 109 (2008) 143–165. [51] W. Sun, N.G. Zamani, Adaptive mesh redistribution for the boundary element method in elastostatics, Comput. Struct. 36 (1990) 1081–1088. [52] L.N. Trefethen, Is Gauss quadrature better than Clenshaw–Curtis?, SIAM Rev. 50 (2008) 67–87.
JID:YJMAA
AID:19947 /FLA
Doctopic: Applied Mathematics
[m3L; v1.166; Prn:11/11/2015; 10:00] P.19 (1-19)
S. Xiang et al. / J. Math. Anal. Appl. ••• (••••) •••–•••
19
[53] L.N. Trefethen, Approximation Theory and Approximation Practice, SIAM, 2013. [54] H. Wang, S. Xiang, Uniform approximations to Cauchy principal value integrals of oscillatory functions, Appl. Math. Comput. 215 (2009) 1886–1894. [55] H. Wang, S. Xiang, On the evaluation of Cauchy principal value integrals of oscillatory functions, J. Comput. Appl. Math. 234 (2010) 95–100. [56] H. Wang, L. Zhang, D. Huybrechs, Asymptotic expansions and fast computation of oscillatory Hilbert transforms, Numer. Math. 123 (2013) 709–743. [57] S. Xiang, X. Chen, H. Wang, Error bounds in Chebyshev points, Numer. Math. 116 (2010) 463–491. [58] S. Xiang, G. He, Y. Cho, On error bounds of Filon–Clenshaw–Curtis quadrature for highly oscillatory integrals, Adv. Comput. Math. 41 (2015) 573–597. [59] Z. Xu, S. Xiang, G. He, Efficient evaluation of oscillatory Bessel Hilbert transforms, J. Comput. Appl. Math. 258 (2014) 57–66.
AID:19947 /FLA
Doctopic: Applied Mathematics
[m3L; v1.166; Prn:11/11/2015; 10:00] P.1 (1-19)
J. Math. Anal. Appl. ••• (••••) •••–•••
Contents lists available at ScienceDirect
Journal of Mathematical Analysis and Applications www.elsevier.com/locate/jmaa
On uniform approximations to hypersingular finite-part integrals ✩ Shuhuang Xiang a , Chunhua Fang a , Zhenhua Xu b a
Department of Applied Mathematics and Statistics, Central South University, Changsha, Hunan 410083, PR China b College of Mathematics and Information Science, Zhengzhou University of Light Industry, Zhengzhou, Henan 450002, PR China
a r t i c l e
i n f o
Article history: Received 11 October 2014 Available online xxxx Submitted by W.L. Wendland Keywords: Finite-part integral Chebyshev approximation Uniform convergence Error bound Stability
a b s t r a c t In this paper, new uniform approximation schemes for computation of hypersingular finite-part integrals are studied. The methods are verified to be supremely qualified for oscillatory integrands. In addition, based on the results on Chebyshev approximations, the uniform convergence and computational error bounds are considered. Preliminary numerical results show the stability and efficiency of the schemes. © 2015 Published by Elsevier Inc.
1. Introduction Evaluation of the hypersingular finite-part integral ˆ1 Im (f ; c) = = −1
f (x) eiωx dx, (x − c)m+1
c ∈ (−1, 1),
(1.1)
can be transferred into the mth-derivative of the Cauchy principle value integral [14,43] ˆ1 1 dm f (x) iωx e dx, Im (f ; c) = − m! dcm x − c
(1.2)
−1
✩ This work is supported partly by NSF of China (No. 11371376), the Innovation-Driven Project and the Mathematics and Interdisciplinary Sciences Project of Central South University. E-mail addresses: [email protected] (S. Xiang), [email protected] (C. Fang), [email protected] (Z. Xu).
http://dx.doi.org/10.1016/j.jmaa.2015.11.002 0022-247X/© 2015 Published by Elsevier Inc.
JID:YJMAA 2
AID:19947 /FLA
Doctopic: Applied Mathematics
[m3L; v1.166; Prn:11/11/2015; 10:00] P.2 (1-19)
S. Xiang et al. / J. Math. Anal. Appl. ••• (••••) •••–•••
where f and its derivatives f (j) (j ≤ m−1) are assumed to be continuous on [−1, 1], f (m) is Hölder continuous [14,43], i2 = −1 and ω is a parameter. For m = 0, the integral is the usual Cauchy principle value integral or finite Hilbert transform. For m = 1, the integral is the Hadamard finite-part integral [28,43,49,50]. Integrals of the form (1.1) appear frequently in boundary element methods (BEMs) and other numerical computations [3,4,31,43,51], particularly in the area of applied mechanics [31,43]. The efficiency of BEMs often depends upon the efficiency of numerical evaluation of such hypersingular integrals. In the case ω = 0, the computation of integral (1.1) has been extensively studied, such as the Gauss-type method [14,16,17,23,32,33,41,42], the (composite) Newton–Cotes method [21,39,49–51], and others [12,23, 27,37]. A typical procedure for approximating (1.1) for ω = 0 starts by subtracting out the singularity at c as follows [14,16,17,23,41]:
Im (f ; c) =
ˆ1 f (x) − m j=0 −1
(x −
f (j) (c) j! (x m+1 c)
− c)j
ˆ m f (j) (c) dx = dx + , j! (x − c)m+1−j j=0 1
(1.3)
−1
followed by applying some ordinary quadrature rules, such as Gaussian [32,33], Fejér or Clenshaw–Curtis ˆ1 1−c dx = log and rule [9,10,37] to the first integral on the right-hand side in (1.3). While − x−c 1+c −1
ˆ1 = −1
ˆ1 dm−j 1 dx dx , = − (x − c)m+1−j (m − j)! dcm−j x − c
j = 1, . . . , m,
−1
can be computed exactly. Although these schemes are simple in general, there are severe numerical cancellations if one of the nodes happens to be close to c [28,29]. Then Gori and Santi [24,25] used the Gauss–Turán quadrature for the Gori–Michelli weight functions to avoid this problem. Milovanović and Spalević [41] generalized the results from [24,25] by using the Chakalov–Popoviciu quadrature. These developed methods can get high accuracy [41], but they involve computation of the higher derivatives of f at quadrature nodes and are dependent on c. More recently, Sun and Wu [49,50] developed (composite) Newton–Cotes rules for the computation of Hadamard finite-part integrals, and showed the superconvergence. However, these methods are possible to suffer from computational efforts or instability [16,53] as the nodes increase. Efficiently uniform approximation algorithms are proposed by Hasegawa and Torii [28,29] for the finite Hilbert transform and Hadamard finite-part integral ˆ1 f (x) dx, − x−c
ˆ1 =
−1
−1
f (x) dx, (x − c)2
respectively, based on Chebyshev approximations to f (x). Let
pN (x) =
N
N ak Tk (x)
k=0
be the Chebyshev approximate polynomial of f (x) at the Clenshaw–Curtis points cos(jπ/N ) (0 ≤ j ≤ N ), where Tk (x) is the Chebyshev polynomial of the first kind given by Tk (x) = cos(kθ) and x = cos(θ), and
JID:YJMAA
AID:19947 /FLA
Doctopic: Applied Mathematics
[m3L; v1.166; Prn:11/11/2015; 10:00] P.3 (1-19)
S. Xiang et al. / J. Math. Anal. Appl. ••• (••••) •••–•••
3
the double prime denotes a sum whose first and last terms are halved. The coefficients aN k can be computed by FFT [7,13,52]. Here we cite a Matlab code from Trefethen [52]. function I = FFT(f, N)
% (N + 1)-pt Fast Fourier Transform
x = cos(pi ∗ (0 : N) /N);
% Chebyshev points of the first kind
fx = feval(f, x)/(2 ∗ N);
% f evaluated at these points
g = real(fft(fx([1 : N + 1N : −1 : 2])));
% FFT
a = [2g(1); g(2 : N) + g(2 ∗ N : −1 : N + 2); 2g(N + 1)];
% Chebyshev coefficients
For the Hilbert transform in the case ω = 0, Hasegawa and Torii [28] introduced the approximation form ˆ1 QN 0 (f ; c)
1−c pN (x) − pN (c) dx + f (c) log , x−c 1+c
= −1
where (pN (x) − pN (c))/(x − c) can be rewritten in terms of Tk (x) as follows N −1 pN (x) − pN (c) N = dk Tk (x). x−c
(1.4)
k=0
Here the prime denotes the summation whose first term is halved. In addition, the coefficients dN k satisfy a three-term recurrence relation [28] N N N dN k−1 = 2ak + 2cdk − dk+1 , k = N − 1, . . . , 1
(1.5)
N N and the first two initial values are given by dN N = 0 and dN −1 = aN , which yields a new integration formula
N/2−1
QN 0 (f ; c)
=2
j=0
dN 1−c 2j + f (c) log 1 − 4j 2 1+c
(1.6)
with O(N log N ) operations. To extend to the Hadamard finite-part, Hasegawa and Torii [29] constructed a uniform scheme for the Hilbert transform and then applied to the Hadamard finite-part integral, by using the following identity [22] Tk+1 (x) − Tk+1 (c) = 2(x − c)
k
Tj (c)Uk−j (x) = 2(x − c)
j=0
k
Tj (x)Uk−j (c),
k ≥ 1,
(1.7)
j=0
sin((k + 1)θ) (x = cos(θ)) is the second kind of Chebyshev polynomial of order k. A uniform sin(θ) scheme [29] for the Hilbert transform I0 (f ; c) is given as follows where Uk (x) =
ˆ1 QN 0 (f ; c)
= −1
=4
1−c pN (x) − pN (c) dx + f (c) log x−c 1+c
N −1 k=0
AN k Tk (c) + f (c) log
1−c , 1+c
(1.8)
JID:YJMAA
AID:19947 /FLA
Doctopic: Applied Mathematics
[m3L; v1.166; Prn:11/11/2015; 10:00] P.4 (1-19)
S. Xiang et al. / J. Math. Anal. Appl. ••• (••••) •••–•••
4
where (N −k−1)/2
AN k
=
aN 2n+k+1 , 2n + 1
n=0
1 ≤ k ≤ N − 1,
which leads to a uniform scheme [27] for the Hadamard finite-part integral I1 (f ; c) QN 1 (f ; c) =
N −1 2f (c) d N 1−c Q0 (f ; c) = 4 − kAN . k Uk−1 (c) + f (c) log dc 1 + c 1 − c2
(1.9)
k=1
These two schemes (1.8) and (1.9) are stable and convergent independent of c as N trends to infinity, based on the complicated proofs in [27–29]. However they require O(N 2 ) operations. Moreover, the higher order finite-part is not considered in [27–29] and there are few results on this topic. In the case |ω| 1, the integrand is highly oscillatory and the computation of (1.1) will suffer from difficulty by the existed standard algorithms even though there are many algorithms on calculation of highly oscillatory integrals with weak or algebraic singularities or Hilbert transforms [8,11,18–20,26,30,34, 35,44,45,54–59]. In this paper, we shall be concerned with uniform approximation schemes for computation of hypersingular finite-part integral (1.1), particularly for |ω| 1. The overall computational complexity for the schemes is of O(N log N ). Without loss of generality, we assume ω ≥ 0. All the numerical tests in this paper are carried out by using Matlab R2012a. The rest of the paper is organized as follows. In Section 2, we present a uniform approximation scheme for the computation of hypersingular finite-part integrals Im (f ; c). In Section 3 we show the uniform convergence and present error bounds for these new schemes. In Section 4, we illustrate the effectiveness and accuracy by preliminary numerical examples. 2. Uniformly convergent schemes The second integral on the right-hand side in (1.3) with weight eiωx can be evaluated by ˆ m m 1 f (j) (c) eiωx dx = = j! (x − c)m+1−j m! j=0 j=0 1
−1
m j
f
(j)
ˆ1 iωx dm−j e dx (c) m−j − dc x−c −1
⎞
⎛
ˆ1 iωx 1 dm ⎝ e = dx⎠ f (c)− m! dcm x−c −1
=:
1 (m) g (c) m!
and then (1.1) has the form of ˆ1 Im (f ; c) =
e −1
iωx
f (x) −
m j=0
(x −
f (j) (c) j! (x m+1 c)
− c)j
dx +
1 (m) g (c), m!
(2.1)
where g(c) is represented by [36] as ˆ1 g(c) = f (c) −1
eiωx dx = eiωc (Ei(1, iω(1 + c)) − Ei(1, iω(c − 1)) + sgn(ω)πi), x−c
(2.2)
JID:YJMAA
AID:19947 /FLA
Doctopic: Applied Mathematics
[m3L; v1.166; Prn:11/11/2015; 10:00] P.5 (1-19)
S. Xiang et al. / J. Math. Anal. Appl. ••• (••••) •••–•••
and Ei(a, z) is the extended exponential integral Ei(a, z) = (z) > 0 or a > 1 and (z) ≥ 0) and dk Ei(a, iω(1 + z)) = (−iω)k Ei(a − k, iω(1 + z)), dz k
´∞ 1
5
t−a e−zt dt = z a−1 Γ(1 − a, z) (a ≥ 1 and
dk Ei(a, iω(z − 1)) = (−iω)k Ei(a − k, iω(z − 1)) dz k
(see [2]), which together yields that g (m) (c) can be exactly computed. Particularly, for ω = 0, g(c) = 1−c . f (c) log 1+c From (2.1) we define a uniform scheme for (1.1) by ˆ1 QN m (f ; c) =
e
p (x) − iωx N
m (x −
−1
(j)
pN (c) j! (x m+1 c)
j=0
− c)j
dx +
1 (m) g (c). m!
(2.3)
The computation of (2.3) can be done by using the second identity in (1.7) repeatedly following Hasegawa and Torii’s work [28,29]. For example, in the case ω = 0, using the second identity in (1.7), we have the ˆ1 f (x) following uniform approximation for the Hadamard finite-part integral = dx: (x − c)2 −1
ˆ1 QN 1 (f ; c)
= −1
=2
2f (c) 1−c pN (x) − pN (c) − pN (c)(x − c) − dx + f (c) log 2 (x − c) 1 + c 1 − c2
N
aN k
k−1
− pN (c) log N
aN k
k−1
N
aN k
=
+
N
−1
+2
k=3
⎞ 1 − c⎠ Tj (x) − Tj (c) dx + Tj (c) log x−c 1+c
2f (c) 1−c 1−c + f (c) log − 1+c 1 + c 1 − c2
Uk−1−j (c)
j−1
Uj−1−m (c)
m=0
k−1
Uk−1−j (c)Tj (c) log
1 + (−1)m 1−c − pN (c) log 2 1−m 1+c
2f (c) 1−c 1−c + f (c) log − 1+c 1 + c 1 − c2
aN k
k=3 N
ˆ1
j=2
k=3
4aN 2
Uk−1−j (c) ⎝
j=0
k=1
+2
j=0
k=1
=4
⎛
k−1
1−c 1−c Uk−2 (c)(4 + 4c log ) − (k − 1)Uk−1 (c) log 1+c 1+c ⎛
1−c + 2Uj−1 (c) + aN Uk−1−j (c) ⎝Tj (c) log k 1+c j=2
+ f (c) log
2f (c) 1−c − . 1 + c 1 − c2
[(j−1)/2]
=1
⎞ Uk−1−2 (c) ⎠ 1 − 42 (2.4)
m+2 In a similar way, QN ) operations. But this extension does suffer from big m (f ; c) can be achieved in O(N difficulty in computation and become useless for m ≥ 1. However, expression (2.3) is useful in the numerical analysis later.
JID:YJMAA
AID:19947 /FLA
Doctopic: Applied Mathematics
[m3L; v1.166; Prn:11/11/2015; 10:00] P.6 (1-19)
S. Xiang et al. / J. Math. Anal. Appl. ••• (••••) •••–•••
6
It is observed, from (1.2) and (2.1), that ˆ1 e
iωx
f (x) −
m j=0
(x −
−1
f (j) (c) j! (x c)m+1
− c)j
1 dm dx = m! dcm
ˆ1 eiωx −1
f (x) − f (c) dx. x−c
Similarly, we have ˆ1 e
m
p (x) − iωx N
(j)
pN (c) j! (x m+1 c)
j=0
(x −
−1
− c)j
1 dm dx = m! dcm
ˆ1 eiωx −1
pN (x) − pN (c) dx, x−c
which, together with (2.3), establishes a new approximate formula for m ≥ 1 QN m (f ; c)
1 dm = m! dcm
ˆ1 eiωx −1
1 (m) pN (x) − pN (c) dx + g (c). x−c m!
(2.5)
Furthermore, inserting (1.4) into (2.5) gives 1 dm QN m (f ; c) = m! dcm
N −1
N dk Mk (ω)
+
k=0
N −1 1 (m) 1 dm dN 1 (m) g (c) = g (c), Mk (ω) mk + m! m! dc m!
(2.6)
k=0
´1 where Mk (ω) = −1 Tk (x)eiωx dx is called modified moments, which can be efficiently computed by the following recursion [47] together with the Oliver’s algorithm [19,46] iωMk+2 (ω) + 2(k + 2)Mk+1 (ω) − 2iωMk (ω) − 2(k − 2)Mk−1 (ω) + iωMk−2 (ω) = 0 for ω = 0. The Matlab codes can be downloaded from [1]. If ω = 0, Mk (ω) = 0 for k being odd, while 2 Mk (ω) = 1−k 2 for k being even. Based on the above discussion, we outline our algorithm with estimation of computational complexity as follows for m ≥ 1: Algorithm 1. Computation of QN m (f ; c): 1. Compute the coefficients aN k (k = 0, 1, . . . , N ) by using FFT in O(N log N ) operations. 2. Compute Mk (ω) (k = 0, . . . , N ) in O(N ) operations by [1] (ω = 0); or Mk (0) = 0 for k being 2 odd, Mk (0) = 1−k 2 for k being even. dm N 3. Compute dcm dk (k = N, . . . , 0) with the following recurrence relation derived from (1.5) in O(N ) operations dm N dcm dN m
N d dcm dk
= 0,
..., m−1
4. Compute
m
= 0,
N N d d = 2 dc m−1 dk+1 + 2c dcm dk+1 −
1 dm m! dcm
dm N dcm dN −m−1
iωc
f (c)e
dm N dcm dk+2 ,
k = N − m − 2, . . . , 0.
(2.7)
[Ei(1, iω(1 + c)) − Ei(1, iω(c − 1)) + sgn(ω)πi] .
5. Calculate QN m (f, c) (2.6) in O(N ) operations. From Steps 1–5, we can see that this new scheme can be achieved in O(N log N ) operations for any m. Moreover, the new scheme has the N -order algebraic accuracy if we do not consider the rounding errors. However, this scheme is not stable when N or m increases. We will see that in Section 4.
JID:YJMAA
AID:19947 /FLA
Doctopic: Applied Mathematics
[m3L; v1.166; Prn:11/11/2015; 10:00] P.7 (1-19)
S. Xiang et al. / J. Math. Anal. Appl. ••• (••••) •••–•••
7
Notice that the linear system (2.7) can be represented by linear system Ax = b as ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝
⎞⎛
1 1 ..
.
1 1 −2c 1 1 −2c 1 1 −2c 1 .. .. .. . . . 1 −2c 1
⎟⎜ ⎟⎜ ⎟⎜ ⎟⎜ ⎟⎜ ⎟⎜ ⎟⎜ ⎟⎜ ⎟⎜ ⎟⎜ ⎟⎜ ⎟⎜ ⎟⎜ ⎟⎜ ⎟⎜ ⎟⎜ ⎠⎝
⎛
⎞ 0 ⎟ ⎟ ⎜ 0 ⎟ ⎟ ⎜ ⎜ ⎟ ⎟ ⎜ .. .. ⎟ ⎟ ⎜ . ⎟ . ⎟ ⎜ ⎟ ⎟ ⎜ dm N ⎟ 0 ⎟ d m ⎜ ⎟ N −m−1 dc ⎟ ⎜ dm−1 N ⎟ dm N ⎟ 2 d m−1 = d ⎜ ⎟. N −m−1 dc dcm N −m−2 ⎟ ⎜ ⎟ m−1 m ⎟ ⎜2 d N d N ⎟ dcm−1 dN −m−2 ⎟ dcm dN −m−3 ⎟ ⎜ ⎟ ⎜ dm−1 N dm N ⎟ ⎜ 2 dcm−1 dN −m−3 ⎟ ⎟ dcm dN −m−4 ⎟ ⎟ ⎟ ⎜ .. .. ⎜ ⎟ ⎠ ⎝ . . ⎠ m m−1 d N d N 2 dcm−1 d1 dcm d0 dm N dcm dN dm N dcm dN −1
⎞
(2.8)
The coefficient matrix is well-conditioned and then the system (2.7) can be efficiently computed by GMRES method based on (2.8) for m ≥ 1. Algorithm 2. Computation of QN m (f ; c): 1. Compute the coefficients aN k (k = 0, 1, . . . , N ) by using FFT with O(N log N ) operations. 2. Compute Mk (ω) (k = 0, . . . , N ) in O(N ) operations by [1] (ω = 0); or Mk (0) = 0 for k being 2 odd, Mk (0) = 1−k 2 for k being even. dm N 3. Compute the coefficients dcm dk (k = N, . . . , 0) by using the GMRES
for (2.8). 4. Compute
1 dm m! dcm
f (c)eiωc [Ei(1, iω(1 + c)) − Ei(1, iω(c − 1)) + sgn(ω)πi] .
5. Calculate QN m (f ; c) in O(N ) operations. Then Scheme (2.6) is uniformly convergent independent of c and ω as N trends to infinity which can be seen in Section 3. Furthermore, the accuracy increases when oscillation becomes faster. Remark 1. In particular, for ω = 0,
1−c 1 dm 1 dm iωc f (c)e f (c) log . [Ei(1, iω(1 + c)) − Ei(1, iω(c − 1)) + sgn(ω)πi] = m! dcm m! dcm 1+c Remark 2. To obtain the higher order for the convergence rate with respect to ω, here we also present a new modified scheme for (2.6) N (f ; c) = QN (f ; c) + 1 Q m m iω
−Pm (1) iω Pm (−1) −iω , e + e (1 − c)m+1 (−1 − c)m+1
(2.9)
m f (j) (c)−p(j) N (c) where Pm (x) = j=0 (x − c)j . In Section 3, we will prove the modified scheme’s accuracy is j! −2 O(ω ). In Section 4, numerical examples are presented to illustrate the effectiveness. Remark 3. For small values of N , Algorithm 1 and Algorithm 2 have nearly same accuracy (see Figs. 2, 5 and Tables 1–3 in Section 4). 3. Uniform convergence and computational error bounds In this section, we consider the uniform convergence and error bound of Scheme (2.6) independent of c and ω. Here we assume that dN k can be computed accurately, which can be done by Algorithm 1 for small values of N and m, or by Algorithm 2 based on the GMRES.
JID:YJMAA
AID:19947 /FLA
Doctopic: Applied Mathematics
[m3L; v1.166; Prn:11/11/2015; 10:00] P.8 (1-19)
S. Xiang et al. / J. Math. Anal. Appl. ••• (••••) •••–•••
8
From (2.1) and (2.2), it follows that ˆ1 Im (f ; c) −
QN m (f ; c)
=
e
iωx
F (x) −
m j=0
(x −
−1
F (j) (c) (x j! m+1 c)
− c)j
dx,
(3.1)
where F (x) = f (x) − pN (x). Thus, the error bound of Scheme (2.6) can be estimated by the remainder of the Taylor’s expansion of F (x) as Im (f ; c) − QN m (f ; c) ≤
2 (m+1) (m+1) − pN f . (m + 1)! ∞
(3.2)
In the following, we will focus on the estimate of the right-hand side of (3.2). Assumption 1. Suppose that f (x) has an absolutely continuous (k0 − 1)st derivative f (k0 −1) on [−1, 1] for ´x some k0 ≥ max{1, m} with f (k0 −1) (x) = f (k0 −1) (−1) + −1 h(y)dy, where h is absolutely integrable and of bounded variation Var(h) < ∞ on [−1, 1], and define Vk0 = inf
f (k0 −1) (x) = f (k0 −1) (−1) + ´ x h(y)dy for all x ∈ [−1, 1] with h being −1 Var(h) . absolutely integrable and of bounded variation
Let {bn } denote the exact spectral coefficients of f (x). Then the Chebyshev series for f is defined as [53]
f (x) =
∞
bj Tj (x),
bj =
j=0
2 π
ˆ1 −1
f (x)Tj (x) √ dx. 1 − x2
Lemma 3.1. (i) (Bernstein [5]) If f is analytic with |f (z)| ≤ M in the region Eρ bounded by the ellipse with foci {−1, 1} and major and minor semiaxis lengths summing to ρ > 1, then for each j ≥ 0, |bj | ≤
2M . ρj
(3.3)
(ii) (Trefethen [53]) Suppose that f (x) satisfies Assumption 1 for some k0 ≥ max{1, m}, then for each j ≥ k + 1, |bj | ≤
2Vk0 . πj(j − 1) · · · (j − k0 )
(3.4)
(iii) (Boyd [6, (4.55), (4.57)]) The coefficients of the interpolating polynomial pN (x) are related to those of f (x) via aN n = bn +
∞ (bn+2jN + b−n+2jN ),
n = 0, 1, . . . , N,
|f (x) − pN (x)| ≤ 2
j=1
∞
|bj |.
(3.5)
j=N +1
Moreover, the error of the interpolation polynomial pN can be estimated by the following formula [6,15] f (x) − pN (x) =
N ∞ b0 − aN (bj − aN bj Tj (x). 0 + j )Tj (x) + 2 j=1 j=N +1
JID:YJMAA
AID:19947 /FLA
Doctopic: Applied Mathematics
[m3L; v1.166; Prn:11/11/2015; 10:00] P.9 (1-19)
S. Xiang et al. / J. Math. Anal. Appl. ••• (••••) •••–•••
9
Then it implies that (m+1)
|f (m+1) (x) − pN
N
(x)| ≤
(m+1)
|bj − aN j | Tj
∞
∞ +
j=m+1
(m+1)
|bj | Tj
∞ .
(3.6)
j=N +1
In addition, for each n ≥ p, it follows from (A.12), (A.13) and (A.52) in [6, pp. 497, 503] that p−1 n2 − k2 dp n+p T (x)| = (±1) n ±1 p dx 2k + 1
dp (p) Tn (x) = n2p−1 (p − 1)!Cn−p (x), dxp
k=0
and (p)
(p)
|Cn−p (x)| ≤ |Cn−p (1)|,
x ∈ [−1, 1],
(p)
where Cn−p (x) is the Gegenbauer polynomial [6]. Then it establishes that p d dxp Tn
= ∞
p−1 k=0
n2p n2 − k2 ≤ . 2k + 1 (2p − 1)!!
(3.7)
This, together with (3.3), Lemma 3.1 and (4.55) in [6, p. 96], gives N ∞ j 2(m+1) j 2(m+1) (m+1) (m+1) bj − aN + (x) − pN (x) ≤ |b | f j j (2m + 1)!! (2m + 1)!! j=m+1 j=N +1
≤ N 2(m+1)
N j=0
∞
≤ N 2(m+1) ∞
|bj |j 2(m+1)
j=N +1
|bj | +
j=N +1
≤2
∞
|bj − aN j |+
∞
|bj |j 2(m+1) ([6, (4.55)])
j=N +1
|bj |j 2(m+1)
j=N +1
≤
∞ j=N +1
4Vk0 ≤ π
=
4Vk0 j 2m+2 ([6, Lemma 3.1]) πj(j − 1) · · · (j − 2m − 1)(j − k0 )k0 −2m−1
N +1 N − 2m
N +1 N − 2m
2m+1 ˆ∞
2m+1
N
1 dx (x − k0 )k0 −2m−1
2Vk0 , (k0 − 2m)π(N − k0 )k0 −2m
where we used the following inequalities for the proof of the last inequality in (3.8) j j j < < ··· < , j−1 j−2 N − 2m
j = N + 1, N + 2, . . . .
From the equations (3.2) and (3.8) we can deduce the following theorem directly.
(3.8)
JID:YJMAA
AID:19947 /FLA
Doctopic: Applied Mathematics
[m3L; v1.166; Prn:11/11/2015; 10:00] P.10 (1-19)
S. Xiang et al. / J. Math. Anal. Appl. ••• (••••) •••–•••
10
Theorem 3.1. Suppose that f (x) satisfies Assumption 1 for some k0 ≥ 2m + 1. Then Scheme (2.6) is uniformly convergent independent of c and satisfies that Im (f ; c) − QN m (f ; c) ≤
N +1 N − 2m
2m+1
4Vk0 . (m + 1)!(k0 − 2m)π(N − k0 )k0 −2m
(3.9)
In particular, if f is analytic with |f (z)| ≤ M in the region Eρ bounded by the ellipse with foci {−1, 1} and major and minor semiaxis lengths summing to ρ > 1, then for each N ≥ 0, based on the line integral for the error of the Chebyshev approximation polynomial pN [15, p. 390], an estimation is given by f (x) − pN (x) =
1 2πi
˛ Eρ
ωN +1 (x)f (z) dz, (z − x)ωN +1 (z)
where ωN +1 (x) = TN +1 (x) − TN −1 (x). By a similar proof to the Cauchy’s differentiation formula [38], (m+1) f (m+1) (x) − pN (x) can be represented by (m+1)
f (m+1) (x) − pN
(m+1)
Therefore, f (m+1) (x) − pN
(x) =
1 2πi
˛ Eρ
f (z) ωN +1 (z)
ωN +1 (x) z−x
(m+1) dz.
(x) can be estimated as follows
(m+1) (m+1) − pN f
∞
(m+1) 1 ˛ ωN +1 (x) f (z) = dz 2πi ωN +1 (z) z−x Eρ
.
(3.10)
∞
From (3.7), we see that p d ω dxp N +1
∞
p 2p d ≤ 2(N + 1) . = (T − T ) N −1 dxp N +1 (2p − 1)!! ∞
(3.11)
Moreover, for j = 0, 1, . . . , m + 1, we have (N + 1)2j+2 (N + 1)2(m+1) (N + 1)2j ≤ ≤ ··· ≤ , j!(2j − 1)!! (j + 1)!(2j + 1)!! (m + 1)!(2m + 1)!! where N ≥
√
2m2 + 3m + 1 − 1. Then for each fixed z ∈ Eρ we get
(m+1) ω N +1 (x) z−x
∞
m+1 m + 1 (m − j + 1)! (j) = ωN +1 (x) m−j+2 (z − x) j j=0 ≤ 2(m + 1)!
m+1 j=0
∞
2j
(N + 1) j!(2j − 1)!!ρm+2−j
≤
m+1 2(N + 1)2(m+1) 1 m+2−j (2m + 1)!! ρ j=0
≤
2(N + 1)2(m+1) . (2m + 1)!!(ρ − 1)
(3.12)
JID:YJMAA
AID:19947 /FLA
Doctopic: Applied Mathematics
[m3L; v1.166; Prn:11/11/2015; 10:00] P.11 (1-19)
S. Xiang et al. / J. Math. Anal. Appl. ••• (••••) •••–•••
11
Additionally, from (1.50) in [40] it follows that for z ∈ Eρ 1 N +1 (ρ − ρ−N −1 − ρN −1 + ρ−N +1 ) ≤ |ωN +1 (z)| = |TN +1 (z) − TN −1 (z)| , 2
(3.13)
which, together with (3.10), (3.12) and (3.13), derives that (m+1) (m+1) − pN f
∞
(m+1) 1 ˛ ωN +1 (x) f (z) = dz 2πi ωN +1 (z) z−x Eρ ≤
−1
(2m +
∞
2(m+1)
M (ρ + ρ )(N + 1) . − ρ−N −1 − ρN −1 + ρ−N +1 )(ρ − 1)
1)!!(ρN +1
Theorem 3.2. If f is analytic with |f (z)| ≤ M in the region Eρ bounded by the ellipse with foci {−1, 1} and √ major and minor semiaxis lengths summing to ρ > 1, then for N ≥ 2m2 + 3m + 1 − 1, the scheme (2.6) is uniformly convergent independent of c and satisfies that Im (f ; c) − QN m (f ; c) ≤
2M (ρ + ρ−1 )(N + 1)2(m+1) . (2m + 1)!!(m + 1)!(ρN +1 − ρ−N −1 − ρN −1 + ρ−N +1 )(ρ − 1)
(3.14)
It is worth noting that the accuracy of Scheme (2.6) increases when oscillation becomes faster. Theorem 3.3. Let f be a smooth function on [−1, 1]. Then
−Pm (1) iω 1 Pm (−1) N −iω e + e I(f ; c) − Qm (f ; c) ≤ iω (1 − c)m+1 (−1 − c)m+1 (m+2) (m+2) (m+3) − pN + 2 f (m+3) − pN f ∞ ∞ , + ω 2 (m + 1)! where Pm (x) =
m j=0
(j)
f (j) (c)−pN (c) (x j!
(3.15)
− c)j .
Before proving Theorem 3.3, we introduce the following lemmas. Lemma 3.2. (See van der Corput [48].) Suppose g(x) is real-valued and smooth in (a, b) and that |g (k) (x)| ≥ 1 for all x ∈ (a, b) for a fixed value of k. Then b ˆ eiωg(x) dx ≤ c(k)ω −1/k a
holds when (i) k ≥ 2, or (ii) k = 1 and g (x) is monotonic. The bound c(k) = 5 · 2k−1 − 2 is independent of g and ω. Lemma 3.3. (See [48, p. 334].) Under the assumptions on g(x) in Lemma 3.2, we can conclude that b ⎤ ⎡ ˆ ˆb eiωg(x) ϕ(x)dx ≤ c(k)ω −1/k ⎣|ϕ(b)| + |ϕ (x)|dx⎦ . a
a
JID:YJMAA
AID:19947 /FLA
Doctopic: Applied Mathematics
[m3L; v1.166; Prn:11/11/2015; 10:00] P.12 (1-19)
S. Xiang et al. / J. Math. Anal. Appl. ••• (••••) •••–•••
12
Proof of Theorem 3.3. Set Pm (x) =
m F (j) (c) j=0
j!
(x − c)j to be the Taylor polynomial of F (x) = f (x) −pN (x)
with degree deg(Pm ) ≤ m. Integrating by parts and from Lemma 3.3, we have ˆ1 F (x) − Pm (x) iωx 1 Im (f ; c) − QN = (f ; c) de m iω (x − c)m+1 −1 1 −Pm (1) iω Pm (−1) −iω ≤ ( e + e ) iω (1 − c)m+1 (−1 − c)m+1 ˆ1
1 F (x) − P (x) m iωx + e dx m+1 iω (x − c) −1 1 −Pm (1) iω Pm (−1) −iω ≤ ( e + e ) m+1 m+1 iω (1 − c) (−1 − c)
F (x) − Pm (x) 3 F (x) − P (x) m + 2 + 2 , ω (x − c)m+1 (x − c)m+1 ∞ ∞
(3.16)
F (m+1) (c) F (x) − Pm (x) F (x) − Pm (x) to define the value of = at x = c. We m+1 x→c (x − c) (m + 1)! (x − c)m+1
where we used the limit lim have
F (x) − Pm (x) (x − c)m+1
⎧ (x))(x − c) − (m + 1)(F (x) − Pm (x)) (F (x) − Pm ⎪ ⎪ , x = c ⎨ (x − c)m+2 = ⎪ F (m+2) (c) ⎪ ⎩ , x=c (m + 2)!
and
F (x) − Pm (x) (x − c)m+1
⎧ ⎨
=
⎩
F (x)−Pm (x))(x−c)2 −2(m+1)(F (x)−Pm (x))(x−c)+(m+1)(m+2)(F (x)−Pm (x)) , (x−c)m+3 (m+3)
2F (c) (m+3)! ,
x = c x = c.
Define G(x) = (F (x) − Pm (x))(x − c) − (m + 1)(F (x) − Pm (x)), it is easy to verify that
G(c) = G (c) = · · · = G(m) (c) = 0 and (m+2) (x))(x − c) = F (m+2) (x)(x − c). G(m+1) (x) = (F (m+2) (x) − Pm
Thus, the Taylor’s expansion of G(m+1) (x) at x = c has the form G(x) =
F (m+2) (ζ)(ζ − c) (x − c)m+1 , (m + 1)!
ζ = c + θ(x − c) (0 < θ < 1),
and for all x ∈ [−1, 1] ⎧ G(x) F (m+2) ∞ |ζ − c| ⎪ ⎪
F (x) − P (x) ⎪ ⎨ (x − c)m+2 ≤ (m + 1)! |x − c| , x = c F (m+2) ∞ m ≤ . = m+1 (m+2) ⎪ (x − c) (m + 1)! ∞ ⎪ F ⎪ , x=c ⎩ (m + 2)!
(3.17)
JID:YJMAA
AID:19947 /FLA
Doctopic: Applied Mathematics
[m3L; v1.166; Prn:11/11/2015; 10:00] P.13 (1-19)
S. Xiang et al. / J. Math. Anal. Appl. ••• (••••) •••–•••
13
Table 1 Absolute errors in 2n + 1-points approximation on Scheme (2.6) based on (2.7) and (2.8) to the Hadamard finite-part integral ˆ1 x6 = dx. (x − c)2 0
1 64
n
c=
c=
2 3 4
0.4336285 7.438494e-15 7.771561e-15
1 16
0.2110565 1.137979e-15 2.803313e-15
c = 0.5
c = 0.99999
0.0125 MP MP
0.7311613 MP MP
Similarly, let G(x) = (F (x) −Pm (x))(x −c)2 −2(m +1)(F (x) −Pm (x))(x −c) +(m +1)(m +2)(F (x) −Pm (x)), there holds
F (x) − Pm (x) (x − c)m+1
=
G(x) (c) = · · · = G (m) (c) = 0, G (m+1) (x) = F (m+3) (x)(x − c)2 , G(c) =G (x − c)m+3
and
F (x) − P (x) F (m+3) m ∞ . ≤ (x − c)m+1 (m + 1)!
(3.18)
A combination of (3.16)–(3.18) establishes the desired result. 2 From the proof of Theorem 3.3, we can deduce the following corollary immediately. Corollary 3.1. Let f be a smooth function on [−1, 1]. Then (m+2) (m+3) (m+2) − pN + 2 f (m+3) − pN f ∞ ∞ N (f ; c) ≤ . I(f ; c) − Q m 2 ω (m + 1)!
(3.19)
4. Numerical examples In this section, all the experiments are performed on the computer with 2.40 GHz laptop personal computer with 2 GB of RAM. We are using the R2012a version of the Matlab system. As stated in the following Tables 1–3, “MP” indicates that all the data have reached machine precision in Matlab. Of particular interest is the fact that for small values of N , Scheme (2.6) based on (2.7) and (2.8) have nearly same accuracy (see Tables 1–3 and Figs. 2 and 5). However, for large values of N , the difference becomes more significant as m increases. Example 1. We consider Scheme (2.6) based on (2.7) and (2.8) respectively for the finite-part integrals [49] ˆ1 = 0
ˆ1 1 x6 (t + 1)6 = dx = dt, (x − c)2 32 (t − 2c + 1)2
c ∈ (0, 1)
−1
(see Fig. 1 and Table 1). The exact solution [49] is 1 1−c 6 3 + c + 2c2 + 3c3 + 6c4 + + 6c5 log . 5 2 c−1 c Here, the GMRES for (2.8) converges at iteration 5 to a solution with relative residual less than 1e − 13.
JID:YJMAA
AID:19947 /FLA
Doctopic: Applied Mathematics
[m3L; v1.166; Prn:11/11/2015; 10:00] P.14 (1-19)
S. Xiang et al. / J. Math. Anal. Appl. ••• (••••) •••–•••
14
Fig. 1. The relative errors of Scheme (2.6) based on (2.7) and (2.8) for
Fig. 2. The relative errors of Scheme (2.6) based on (2.7) and (2.8) for respectively.
´ x6 =1 0 (x−c)2
dx with c = 0.549 and c = 0.99999, respectively.
´1 1−a2 = −1 (1−2ax+a2 )(x−c)3
dx with c = 0.549 and c = 0.99999,
Example 2. We illustrate the accuracy of Scheme (2.6) for computing the finite-part integrals [29] ˆ1 = −1
1 − a2 dx, (1 − 2ax + a2 )(x − c)3
ˆ1 = −1
1 − a2 dx, (1 − 2ax + a2 )(x − c)4
a = 0.8,
c ∈ (−1, 1)
1 dI1 (a, c) 1 d2 I1 (a, c) and , respectively, where 2 dc 6 dc2
1−c 1−a 2(1 − a2 ) 2a(1 − a2 ) log − 2 log − . I1 (a; c) = (1 − 2ac + a2 )2 1+c 1+a (1 − 2ac + a2 )(1 − c2 )
(see Figs. 2–4). The exact solutions are
Here, the GMRES for (2.8) converges at iteration 1 to a solution with relative residual less than 1e − 10 for N ≥ 200. Example 3. We demonstrate Scheme (2.6) for computing the Cauchy principal value integral I0 (f, c) = ´ 1 x sin(ωx) − e x−c dx with c = 0 (see Table 2). From maple 14, we see that the integral I0 (f, c) can be explicitly −1 expressed by the following formula
JID:YJMAA
AID:19947 /FLA
Doctopic: Applied Mathematics
[m3L; v1.166; Prn:11/11/2015; 10:00] P.15 (1-19)
S. Xiang et al. / J. Math. Anal. Appl. ••• (••••) •••–•••
Fig. 3. The relative errors of Scheme (2.6) based on (2.7) and (2.8) for respectively.
Fig. 4. The relative errors of Scheme (2.6) based on (2.7) and (2.8) for respectively.
15
´ 1−a2 =1 −1 (1−2ax+a2 )(x−c)3
dx with c = 0.549 and c = 0.99999,
´ 1−a2 =1 −1 (1−2ax+a2 )(x−c)4
dx with c = 0.549 and c = 0.99999,
Table 2
ˆ1 x sin(ωx) Absolute errors in 2n + 1-point approximation to the finite-part integral − e dx by Scheme (2.6). x −1
n
ω = 10
ω = 102
ω = 104
ω = 106
1 2 3 4
1.396989e-3 3.163167e-5 −1.418378e-9 MP
8.533057e-6 1.1859500e-7 7.234213e-13 MP
5.447682e-10 8.905765e-12 MP MP
6.217249e-14 8.881784e-16 MP MP
I0 (f, c) =
where E1 (z) =
´∞ z
i c(1+iω) e E1 (1 − c)(1 + iω) − E1 − (1 + c)(1 + iω) + ln(−1 − iω) 2 ! − ln(1 + iω) + e−c(−1+iω) E1 (1 + c)(−1 + iω) − E1 (1 − c)(1 − iω) ! + ln(1 − iω) − ln(−1 + iω) , −1 < c < 1,
e−t t dt
is the exponential integral.
(4.1)
JID:YJMAA
AID:19947 /FLA
16
Doctopic: Applied Mathematics
[m3L; v1.166; Prn:11/11/2015; 10:00] P.16 (1-19)
S. Xiang et al. / J. Math. Anal. Appl. ••• (••••) •••–•••
Fig. 5. The relative errors scaled by ω of Scheme (2.9) based on (2.7) and (2.8), respectively, for the computation of with ω = 108 .
Fig. 6. The absolute errors scaled by ω of Scheme (2.6) based on (2.7) for the computation of k = 1 : 2000 and N = 9.
´ =1 −1
´ =1 −1
ex sin(ωx) dx x2
ex sin(ωx) dx with ω = 1 + 10k, x2
´1 We also consider Scheme (2.6) and (2.9) for computing the finite-part integral I1 (f, 0) = =−1 ex sin(ωx) x2 dx d (see Figs. 5–7 and Table 3). The exact solution is dc I0 (f, c) c=0 , where I0 (f, c) is denoted by (4.1). Here, the GMRES for (2.8) converges at iteration 3 to a solution with relative residual less than 1e − 10 for N ≥ 200. As can be seen from Figs. 6–7, there is a good agreement between the theory results and the computations. 5. Conclusions We have presented efficient methods for the computation of hypersingular finite-part integrals. The convergence analysis for the corresponding methods is also given. The new uniform scheme allows calculating the Hadamard finite-part integral in O(N logN ) operations, while the existing algorithms require O(N 2 ) operations. Moreover, this method fills the gaps in fast computation of hypersingular integrals, especially for oscillatory integrands.
JID:YJMAA
AID:19947 /FLA
Doctopic: Applied Mathematics
[m3L; v1.166; Prn:11/11/2015; 10:00] P.17 (1-19)
S. Xiang et al. / J. Math. Anal. Appl. ••• (••••) •••–•••
Fig. 7. The absolute errors scaled by ω 2 of Scheme (2.9) based on (2.7) for computation of k = 1 : 2000 and N = 9.
´ =1 −1
17
ex sin(ωx) dx with ω = 1 + 10k, x2
Table 3
ˆ1 x sin(ωx) dx with Scheme (2.9). The absolute errors of 2n + 1-point approximation values to the finite-part integral = e x2 −1
n
ω = 102
ω = 104
ω = 106
ω = 108
1 2 3 4
1.458692e-6 8.270350e-9 7.00000e-15 MP
9.380800e-11 7.130000e-13 1.000000e-15 MP
1.000000e-14 MP MP MP
MP MP MP MP
Acknowledgments The authors are grateful to the anonymous referee for his constructive comments and helpful suggestions for improvement of this paper greatly. The authors also thank Dr. Guo He and Junjie Ma for their checking up every detail. References [1] http://www.unavarra.es/personal/victor_dominguez/clenshawcurtisrule/index.html. [2] M. Abramowitz, I.A. Stegun, Handbook of Mathematical Functions, National Bureau of Standards, Washington, DC, 1964. [3] M. Ainsworth, B. Guo, An additive Schwarz preconditioner for p-version boundary element approximation of the hypersingular operator in three dimensions, Numer. Math. 85 (2000) 343–366. [4] G. Bao, W. Sun, A fast algorithm for the electromagnetic scattering from a large cavity, SIAM J. Sci. Comput. 27 (2005) 553–574. [5] S.N. Bernstein, Sur l’ordre de la meilleure approximation des fonctions continues par les polynômes de degré donné, Mem. Cl. Sci. Acad. Roy. Belg. 4 (1912) 1–103. [6] J.H. Boyd, Chebyshev and Fourier Spectral Methods, Dover Publications, Inc., New York, 2000. [7] M. Branders, R. Piessens, An extension of Clenshaw–Curtis quadrature, J. Comput. Appl. Math. 1 (1975) 55–65. [8] M.R. Capobiancoa, G. Criscuoloa, On quadrature for Cauchy principal value integrals of oscillatory functions, J. Comput. Appl. Math. 156 (2003) 471–486. [9] M.M. Chawla, N. Jayarajan, Quadrature formulas for Cauchy principle value integrals, Computing 15 (1975) 347–355. [10] M.M. Chawla, S. Kumar, Convergence of quadrature for Cauchy principle value integrals, Computing 23 (1979) 67–72. [11] R. Chen, C. An, On evaluation of Bessel transforms with oscillatory and algebraic singular integrands, J. Comput. Appl. Math. 264 (2014) 71–81. [12] U.J. Choi, S.W. Kim, B.I. Yun, Improvement of the asymptotic behaviour of the Euler–Maclaurin formula for Cauchy principle value and Hadamard finite-part integrals, J. Comput. Appl. Math. 61 (2004) 496–513. [13] C.W. Clenshaw, A.R. Curtis, A method for numerocal integration on an automatic computer, Numer. Math. 2 (1960) 197–205.
JID:YJMAA 18
AID:19947 /FLA
Doctopic: Applied Mathematics
[m3L; v1.166; Prn:11/11/2015; 10:00] P.18 (1-19)
S. Xiang et al. / J. Math. Anal. Appl. ••• (••••) •••–•••
[14] G. Criscuolo, A new algorithm for Cauchy principle value and Hadamard finite-part integrals, J. Comput. Appl. Math. 78 (1997) 255–275. [15] G. Dahlquist, A. Björck, Numerical Methods in Scientific Computing, SIAM, Philadelphia, PA, 2007. [16] P.J. Davis, P. Rabinowitz, Methods of Numerical Integration, second edition, Academic Press, 1984. [17] K. Diethelm, A method for the practical evaluation of the Hilbert transform on the real line, J. Comput. Appl. Math. 112 (1999) 45–53. [18] V. Domínguez, Filon–Clenshaw–Curtis rules for a class of highly-oscillatory integrals with logarithmic singularities, J. Comput. Appl. Math. 261 (2014) 299–319. [19] V. Domínguez, I.G. Graham, V.P. Smyshlyaev, Stability and error estimates for Filon–Clenshaw–Curtis rules for highlyoscillatory integrals, IMA J. Numer. Anal. 31 (2011) 1253–1280. [20] V. Domínguez, I.G. Graham, M. Kim, Filon–Clenshaw–Curtis rules for highly-oscillatory integrals with algebraic singularities and stationary points, SIAM J. Numer. Anal. 51 (2013) 1542–1566. [21] Q. Du, Evaluations of certain hypersingular integrals on interval, Internat. J. Numer. Methods Engrg. 51 (2001) 1195–1210. [22] D. Elliott, Truncation errors in two Chebyshev series approximations, Math. Comp. 19 (1965) 234–248, 445–462. [23] D. Elliott, E. Venturino, Sigmoidal transformations and the Euler–Maclaurin expansion for evaluating Hadamard finite-part integrals, Numer. Math. 77 (1997) 453–465. [24] L. Gori, E. Santi, On the evaluation of Hilbert transforms by means of a particular class of Turán quadrature rules, Numer. Algorithms 10 (1995) 27–39. [25] L. Gori, E. Santi, Quadrature rules based on s-orthogonal polynomials for evaluating integrals with strong singularities, in: W. Gautschi, G.H. Golub, G. Opfer (Eds.), Applications and Computation of Orthogonal Polynomials, Oberwolfach, 1998, in: Internat. Ser. Numer. Math., vol. 131, Birkhäuser, Basel, 1999, pp. 109–119. [26] A.I. Hascelik, Suitable Gauss and Filon-type methods for oscillatory integrals with an algebraic singularity, Appl. Numer. Math. 59 (2009) 101–118. [27] T. Hasegawa, Uniform approximations to finite Hilbert transform and its derivatives, J. Comput. Appl. Math. 163 (2004) 127–138. [28] T. Hasegawa, T. Torii, An automatic quadrature for Cauchy principle value integrals, Math. Comp. 56 (1991) 741–754. [29] T. Hasegawa, T. Torii, Hilbert and Hadamard transforms by generalized Chebyshev expansion, J. Comput. Appl. Math. 51 (1994) 71–83. [30] G. He, S. Xiang, An improved algorithm for the evaluation of Cauchy principal value integrals of oscillatory functions and its application, J. Comput. Appl. Math. 280 (2015) 1–13. [31] C. Hui, D. Shia, Evaluations of hypersingular integrals using Gaussian quadrature, Internat. J. Numer. Methods Engrg. 44 (1999) 205–214. [32] N.I. Ioakimidid, Application of finite-part integrals to the singular integral equation of crack problems in plane and three-dimensional elasticity, Acta Mech. 45 (1982) 31–47. [33] N.I. Ioakimidid, On the uniform convergence of Gaussian quadrature rules for Cauchy principle value integrals and their derivatives, Math. Comp. 44 (1985) 191–198. [34] H. Kang, S. Xiang, On the calculation of highly oscillatory integrals with an algebraic singularity, Appl. Math. Comput. 217 (2010) 3890–3897. [35] H. Kang, S. Xiang, Efficient quadrature of highly oscillatory integrals with algebraic singularities, J. Comput. Appl. Math. 237 (2013) 576–588. [36] P. Keller, A practical algorithm for computing Cauchy principal value integrals of oscillatory functions, Appl. Math. Comput. 218 (2012) 4988–5001. [37] P. Kim, S. Ahn, U. Choi, On the Chebyshev quadrature rule of finite part integrals, J. Comput. Appl. Math. 180 (2005) 147–159. [38] A. Lars, Complex Analysis: An Introduction to the Theory of Analytic Functions of One Complex Variable, McGraw–Hill, New York, 1979. [39] P. Linz, On the approximate computation of certain strongly singular integrals, Computing 35 (1985) 345–353. [40] J.C. Mason, D.C. Handscomb, Chebyshev Polynomials, Chapman and Hall/CRC, New York, 2003. [41] G.V. Milovanović, M. Spalević, Quadrature rules with multiple nodes for evaluating integrals with strong singularities, J. Comput. Appl. Math. 189 (2006) 689–702. [42] G. Monegato, Numerical evaluation of hypersingular integrals, J. Comput. Appl. Math. 50 (1997) 9–31. [43] G. Monegato, Definitions, properties and applications of finite-part integrals, J. Comput. Appl. Math. 229 (2009) 425–439. [44] G.E. Okecha, Quadrature formulae for Cauchy principal value integrals of oscillatory kind, Math. Comp. 49 (1987) 259–268. [45] G.E. Okecha, Hermite interpolation and a method for evaluating cauchy principal value integrals of oscillatory kind, Kragujevac J. Math. 29 (2006) 91–98. [46] J. Oliver, The numerical solution of linear recurrence relations, Numer. Math. 11 (1968) 349–360. [47] R. Piessens, M. Branders, On the computation of Fourier transforms of singular functions, J. Comput. Appl. Math. 43 (1992) 159–169. [48] E. Stein, Harmonic Analysis: Real-Variable Methods, Orthogonality, and Oscillatory Integrals, Princeton University Press, Princeton, NJ, 1993. [49] W. Sun, J. Wu, Interpolatory quadrature rules for Hadamard finite-part integrals and their supperconvergence, IMA J. Numer. Anal. 28 (2008) 580–597. [50] W. Sun, J. Wu, The supperconvergence of Newton–Cotes rules for the Hadamard finite-part integral on an interval, Numer. Math. 109 (2008) 143–165. [51] W. Sun, N.G. Zamani, Adaptive mesh redistribution for the boundary element method in elastostatics, Comput. Struct. 36 (1990) 1081–1088. [52] L.N. Trefethen, Is Gauss quadrature better than Clenshaw–Curtis?, SIAM Rev. 50 (2008) 67–87.
JID:YJMAA
AID:19947 /FLA
Doctopic: Applied Mathematics
[m3L; v1.166; Prn:11/11/2015; 10:00] P.19 (1-19)
S. Xiang et al. / J. Math. Anal. Appl. ••• (••••) •••–•••
19
[53] L.N. Trefethen, Approximation Theory and Approximation Practice, SIAM, 2013. [54] H. Wang, S. Xiang, Uniform approximations to Cauchy principal value integrals of oscillatory functions, Appl. Math. Comput. 215 (2009) 1886–1894. [55] H. Wang, S. Xiang, On the evaluation of Cauchy principal value integrals of oscillatory functions, J. Comput. Appl. Math. 234 (2010) 95–100. [56] H. Wang, L. Zhang, D. Huybrechs, Asymptotic expansions and fast computation of oscillatory Hilbert transforms, Numer. Math. 123 (2013) 709–743. [57] S. Xiang, X. Chen, H. Wang, Error bounds in Chebyshev points, Numer. Math. 116 (2010) 463–491. [58] S. Xiang, G. He, Y. Cho, On error bounds of Filon–Clenshaw–Curtis quadrature for highly oscillatory integrals, Adv. Comput. Math. 41 (2015) 573–597. [59] Z. Xu, S. Xiang, G. He, Efficient evaluation of oscillatory Bessel Hilbert transforms, J. Comput. Appl. Math. 258 (2014) 57–66.