The Modified Trapezoidal Rule for Computing Hypersingular Integral on Interval

Abstract

The modified trapezoidal rule for the computation of hypersingular integrals in boundary element methods is discussed. When the special function of the error functional equals zero, the convergence rate is one order higher than the general case. A new quadrature rule is presented and the asymptotic expansion of error function is obtained. Based on the error expansion, not only do we obtain a high order of accuracy, but also a posteriori error estimate is conveniently derived. Some numerical results are also reported to confirm the theoretical results and show the efficiency of the algorithms.

1. Introduction

Consider the following integral: 𝐼𝑓,𝑠∶=^𝑏_𝑎𝑓(𝑥)(𝑥−𝑠)^𝑝+1𝑑𝑥,𝑠∈(𝑎,𝑏),𝑝=1,2, (1) where ⨎^𝑏_𝑎 denotes a Hadamard finite-part integral (𝑝=1 is called hypersingular integral and 𝑝=2 is called supersingular integral) and 𝑠 is the singular point. The formulation of these classes of boundary value problems in terms of hypersingular integral equations has drawn lots of interest. Many scientific and engineering problems, such as acoustics, electromagnetic scattering, and fracture mechanics, can be reduced to boundary integral equations with hypersingular kernels. There exist several definitions, equivalent mathematically, for such kind of integrals in some literatures.

We mention the following one: ^𝑏_𝑎𝑓(𝑥)(𝑥−𝑠)²𝑑𝑥=lim_𝜀→0𝑎^𝑠−𝜀𝑓(𝑥)(𝑥−𝑠)²𝑑𝑥+^𝑏_𝑠+𝜀𝑓(𝑥)(𝑥−𝑠)²(𝑠)𝑑𝑥−2𝑓𝜀,𝑠∈(𝑎,𝑏). (2)

Accurate calculation of boundary element methods (BEM) arising in boundary integral equations has been a subject of intensive research in recent years. The hypersingular integrals have certain properties different from regular and weak singular integrals. One of the major problems arising from boundary element method, for solving such integral equations, is how to evaluate the hypersingular integrals on the interval or on the circle efficiently.

Hypersingular integral must be considered in Hadamard finite-part sense. Numerous works [1–18] have been devoted towards developing efficient quadrature formulas. In 1983, the series expansion of hypersingular integral kernel on circle was firstly suggested by Yu [19]. He solved the harmonic and biharmonic natural boundary integral equations successfully. The Newton-Cotes methods to compute the hypersingular integral on interval were firstly studied by Linz [20] with generalized trapezoidal and Simpson rules which fail altogether when the singular point 𝑠 is close to a mesh point. In order to make the mesh be selected in such a way that 𝑠 falls near the center of a subinterval, two shorter subintervals at the end of the interval were allowed. Then Yu [21] gave new quadrature formulae to compute the case of singular point coinciding with the mesh point which presented that the error estimate is 𝑂(ℎ|lnℎ|). In 1999, Wu and Yu [22] presented simple, easy to be implemented methods not affected by the location of singular point with calculation of double. In recent years, the case of singular point coincided with the mesh point, and Wu et al. [23] presented a modified trapezoidal rule and proved the 𝑂(ℎ) convergence rate.

In this paper, for the case of singular point coinciding with the mesh point a new quadrature rule is introduced. Based on the expansion of the error functional, the error estimate is presented and a posteriori error estimate is given. Then not only do we obtain a high order of accuracy, but also a posteriori error estimate is conveniently derived.

The rest of this paper is organized as follows. In Section 2, after introducing some basic formulas of the general (composite) trapezoidal rule and notations, we present our main result. In Section 3, the corresponding theoretical analysis is given. Finally, several numerical examples are given to validate our analysis.

2. Main Result

Let 𝑎=𝑥₀<𝑥₁<⋯<𝑥_𝑛−1<𝑥_𝑛=𝑏 be a uniform partition of the interval [𝑎,𝑏] with mesh size ℎ=(𝑏−𝑎)/𝑛 and set 𝑥_0𝑖=𝑥_𝑖−1+ℎ6,𝑖=1,2,…,𝑛; (3) then we get the new partition: 𝑎=𝑥₀₀<𝑥₀₁<⋯<𝑥_0𝑛<𝑥_0,𝑛+1=𝑏. (4)

We define 𝑓_𝐿(𝑥), the modified trapezoidal interpolation for 𝑓(𝑥), as 𝑓_𝐿(𝑥)=𝑥−𝑥_0𝑖𝑥_0,𝑖+1−𝑥_0𝑖𝑓𝑥_0,𝑖+1+𝑥_0,𝑖+1−𝑥𝑥_0,𝑖+1−𝑥_0𝑖𝑓𝑥_0𝑖,𝑥∈𝑥_0𝑖,𝑥_0,𝑖+1,0≤𝑖≤𝑛, (5) and a linear transformation 𝑥=𝑥_0𝑖(𝜏)()𝑥∶=𝜏+1^0,𝑖+1−𝑥_0𝑖2+𝑥_0𝑖,𝑖=1,…,𝑛−1,𝜏∈[−1,1], (6) from the reference element [−1,1] to the subinterval [𝑥_0𝑖,𝑥_0,𝑖+1]. For the two subintervals [𝑎,𝑥₀₁] and [𝑥_0𝑛,𝑏] near the end of the interval [𝑎,𝑏], 𝜏 values in [2/3,1] and [−1,2/3], respectively.

Replacing 𝑓(𝑥) in (1) with 𝑓_𝐿(𝑥) gives the new composite trapezoidal rule: 𝐼_𝑛𝑓,𝑠∶=^𝑏_𝑎𝑓_𝐿(𝑥)(𝑥−𝑠)²𝑑𝑥=^𝑛+1_𝑖=0𝜔_𝑗(𝑠)𝑓𝑥_0𝑖=𝐼𝑓,𝑠−𝐸_𝑛𝑓, (7) where 𝜔_𝑖(𝑠) is the Cotes coefficients: 𝜔_𝑖(𝑠)=1−𝛿_𝑖0𝑥_0,𝑖−𝑥_0,𝑖−1𝑥ln^0𝑖−𝑠𝑥_0,𝑖−1−−𝑠1−𝛿_𝑖,𝑛+1𝑥_0,𝑖+1−𝑥_0𝑖𝑥ln^0,𝑖+1−𝑠𝑥_0𝑖+𝛿−𝑠𝑥^𝑖0₀₀+𝛿−𝑠𝑠−𝑥^𝑖,𝑛+1_0,𝑛+1, (8)0≤𝑖≤𝑛+1, 𝛿_𝑖𝑗 denotes the Kronecker delta, and 𝐸_𝑛(𝑓) denotes the error functional.

Theorem 1. Assume 𝑓(𝑥)∈𝐶^1+𝛼[𝑎,𝑏], 𝛼∈[0,1). For the trapezoidal rule 𝐼_𝑛(𝑓,𝑠) defined in (7), there exists a positive constant 𝐶, independent of ℎ and 𝑠, such that 𝐸_𝑛𝑓≤𝐶|lnℎ|+ln𝛾(𝜏)ℎ^𝛼, (9) where 𝛾(𝜏)=min_{0≤𝑖≤𝑛+1}𝑠−𝑥_𝑖ℎ=1−|𝜏|2,𝜏∈(−1,1). (10)

Proof. Let 𝑅(𝑥)=𝑓(𝑥)−𝑓_𝐿(𝑥); then we have |𝑅(𝑥)|≤𝐶ℎ^1+𝛼, as 𝐼𝑓,𝑠−𝐼_𝑛𝑓,𝑠=^𝑏_𝑎𝑓(𝑥)−𝑓_𝐿(𝑥)(𝑥−𝑠)²=𝑑𝑥_{𝑥0,𝑚+1𝑥0𝑚}𝑓(𝑥)−𝑓_𝐿(𝑥)(𝑥−𝑠)²+𝑑𝑥_{𝑛𝑖=0,𝑖≠𝑚𝑥0,𝑖+1𝑥0𝑖}𝑓(𝑥)−𝑓_𝐿(𝑥)(𝑥−𝑠)²𝑑𝑥. (11) For the first part of (11), since 𝑅(𝑥)∈𝐶^1+𝛼[𝑎,𝑏], by Taylor expansion, we have |𝑅(𝑥)|≤𝐶ℎ^{1+𝛼−𝑖},𝑖=0,1,2. (12)
By the definition of finite-part integral, ^𝑏_𝑎𝑓(𝑥)(𝑥−𝑠)²()𝑓(𝑠)𝑑𝑥=𝑏−𝑎(𝑠−𝑎)(𝑏−𝑠)+𝑓^′(𝑠)ln𝑏−𝑠𝑠−𝑎+^𝑏_𝑎𝑓(𝑥)−𝑓(𝑠)−𝑓^′(𝑠)(𝑥−𝑠)(𝑥−𝑠)²𝑑𝑥, (13) we have ^𝑥0,𝑚+1_𝑥0𝑚𝑅(𝑥)(𝑥−𝑠)²=(𝑠)𝑑𝑥ℎ𝑅𝑠−𝑥_0𝑚𝑥_0,𝑚+1−𝑠+𝑅^′(𝑠)𝑥ln^0,𝑚+1−𝑠𝑠−𝑥_0𝑚+^𝑥0,𝑚+1_𝑥0𝑚𝑅(𝑥)−𝑅(𝑠)−𝑅^′(𝑠)(𝑥−𝑠)(𝑥−𝑠)²𝑑𝑥. (14) Now, we estimate the right hand side of (14) term by term. Since 𝑅(𝑥_0𝑚)=0, we have ℎ𝑅(𝑠)𝑠−𝑥_0𝑚𝑥_0,𝑚+1=ℎ𝑅(𝑠)𝑥−𝑠−𝑅_0𝑚𝑠−𝑥_0𝑚𝑥_0,𝑚+1=−𝑠ℎ𝑅^′_𝑚𝜉_𝑚𝑠−𝑥_0,𝑚+1≤𝐶ℎ^𝛼,𝜉_𝑚∈𝑥_0𝑚,𝑅,𝑠_′(𝑠)𝑥ln^0,𝑚+1−𝑠𝑠−𝑥_0𝑚≤𝐶ln𝛾(𝜏)ℎ^𝛼,_{𝑥0,𝑚+1𝑥0𝑚}𝑅(𝑥)−𝑅(𝑠)−𝑅^′(𝑠)(𝑥−𝑠)(𝑥−𝑠)²𝑑𝑥≤𝐶ℎ^𝛼,𝜂_𝑚∈𝑥_0𝑚,𝑥_0,𝑚+1. (15) For the second part of (11), ^𝑥01_𝑎+^𝑛−1_{𝑖=1,𝑖≠𝑚}^𝑥0,𝑖+1_𝑥0𝑖+^𝑏_𝑥0𝑛𝑅(𝑥)(𝑥−𝑠)²≤𝑑𝑥_𝑥01𝑎+^𝑛−1_{𝑖=1,𝑖≠𝑚−1,𝑚,𝑚+1}^𝑥0,𝑖+1_𝑥0𝑖+^𝑏_𝑥0𝑛𝑅(𝑥)(𝑥−𝑠)²+𝑑𝑥_{𝑥0,𝑚𝑥0,𝑚−1}𝑅(𝑥)(𝑥−𝑠)²𝑑𝑥+^𝑥0,𝑚+2_𝑥0,𝑚+1𝑅(𝑥)(𝑥−𝑠)²≤𝑑𝑥_𝑥01𝑎+^𝑛−1_{𝑖=1,𝑖≠𝑚−1,𝑚,𝑚+1}^𝑥0,𝑖+1_𝑥0𝑖+^𝑏_𝑥0𝑛𝑅(𝑥)(𝑥−𝑠)²+𝑑𝑥_{𝑥0𝑚𝑥0,𝑚−1}𝑅(𝑥)−𝑅𝑥_0𝑚(𝑥−𝑠)²+𝑑𝑥_{𝑥0,𝑚+2𝑥0,𝑚+1}𝑅(𝑥)−𝑅𝑥_0,𝑚+1(𝑥−𝑠)²=𝑑𝑥_𝑥01𝑎+^𝑛−1_{𝑖=1,𝑖≠𝑚−1,𝑚,𝑚+1}^𝑥0,𝑖+1_𝑥0𝑖+^𝑏_𝑥0𝑛𝑅(𝑥)(𝑥−𝑠)²+𝑑𝑥_{𝑥0𝑚𝑥0,𝑚−1}𝑅^′𝜉_0𝑚𝑥−𝑥_0𝑚(𝑥−𝑠)²𝑑𝑥+^𝑥0,𝑚+2_𝑥0,𝑚+1𝑅^′𝜉_0,𝑚+1𝑥−𝑥_0,𝑚+1(𝑥−𝑠)²𝑑𝑥≤𝐶ℎ^1+𝛼𝑥01_𝑎+^𝑛−1_{𝑖=1,𝑖≠𝑚−1,𝑚,𝑚+1}^𝑥0,𝑖+1_𝑥0𝑖+^𝑏_𝑥0𝑛1(𝑥−𝑠)²𝑑𝑥+𝐶ℎ^{𝛼𝑥0𝑚}_{𝑥0,𝑚−1}𝑑𝑥||+𝑥−𝑠^𝑥0,𝑚+2_𝑥0,𝑚+1𝑑𝑥||(𝜏+|ℎ𝑥−𝑠≤𝐶ln𝛾)lnℎ|^𝛼. (16)
Combining (14), (15), and (16) leads to (9) and the proof is completed.

Firstly, we set 𝜙₁(𝑥)=−12⨍¹₋₁𝜏²−1|𝑥|−1𝜏−𝑥𝑑𝜏,<1,2∫¹₋₁𝜏²−1|𝑥|𝜏−𝑥𝑑𝜏,>1; (17) then we have 𝜙^′₁(𝑥)=−12⨎¹₋₁𝜏²−1(𝜏−𝑥)²|𝑥|−1𝑑𝜏,<1,2∫¹₋₁𝜏²−1(𝜏−𝑥)²𝑑𝜏,|𝑥|>1. (18) By straight calculation, we get 𝜙^′₁(𝑥)=𝑥log1+𝑥1−𝑥−2. (19) Let 𝐽∶=(−∞,−1)⋃(−1,1)⋃(1,+∞) and the operator 𝑊∶𝐶(𝐽)→𝐶(−1,1) be defined as 𝑊𝑓(𝜏)∶=𝑓(𝜏)+^∞_𝑖=1𝑓()(,().2𝑖+𝜏+𝑓−2𝑖+𝜏)𝜏∈−1,1 (20) Obviously, the operator 𝑊 is linear operator. Then we set 𝑆₁(𝜏)=𝑊𝜙^′₁(𝜏)=𝜙^′₁(𝜏)+^∞_𝑖=1𝜙^′₁(2𝑖+𝜏)+𝜙^′₁(−2𝑖+𝜏). (21)

Now we present our main results below. The proof will be given in the next section.

Theorem 2. Assume 𝑓(𝑥)∈𝐶³[𝑎,𝑏]. For the trapezoidal rule 𝐼_𝑛(𝑓,𝑠) defined in (7), there exists a positive constant 𝐶, independent of ℎ and 𝑠, such that 𝐸_𝑛𝑓=𝑓^′′(𝑠)ℎ2𝑆₁(𝜏)+ℛ_𝑓(𝑠), (22) where 𝑠=𝑥_0𝑖+(1+𝜏)ℎ/2, 𝑖=1,2,…,𝑛, and ℛ_𝑓(𝑠)≤𝐶𝜂(𝑠)+|lnℎ|+ln𝛾(𝜏)ℎ² (23) (𝛾(𝜏) defined as (10)) and 𝜂(𝑠)1=max,1𝑠−𝑎.𝑏−𝑠 (24)

3. Proof of Main Results

3.1. Preliminaries

In the following section, 𝐶 denotes certain constant independent of ℎ and 𝑠, and its value varies with places.

Lemma 3. Assume that 𝑓(𝑥)∈𝐶³[𝑎,𝑏] and 𝑓_𝐿(𝑥) are defined by (5); there holds 𝑓(𝑥)−𝑓_𝐿(𝑥)=𝑓^′′(𝑠)2𝑥−𝑥_0,𝑖+1𝑥−𝑥_0𝑖+ℛ¹_𝑖(𝑥)+ℛ²_𝑖(𝑥), (25) where ℛ¹_𝑖(𝑥)=𝑥−𝑥_0,𝑖+1𝑥−𝑥_0𝑖6𝑥_0,𝑖+1−𝑥_0𝑖𝑓⁽³⁾𝜉_0𝑖𝑥−𝑥_0,𝑖+1²−𝑓⁽³⁾𝜂_0𝑖𝑥−𝑥_0𝑖²,ℛ²_𝑖(𝑥)=𝑓⁽³⁾𝛽_0𝑖2𝑥−𝑥_0,𝑖+1𝑥−𝑥_0𝑖(𝑥−𝑠), (26)𝜉_0𝑖,𝜂_0𝑖,𝛽_0𝑖∈(𝑥_0𝑖,𝑥_0,𝑖+1), and ℛ¹_𝑖(𝑥)≤𝐶ℎ³. (27)

Proof. Taking Taylor expansion for 𝑓(𝑥_0𝑖), 𝑓(𝑥_0,𝑖+1) at 𝑥, there holds 𝑓(𝑥)−𝑓_𝐿(𝑥)=𝑓^′′(𝑥)2!𝑥−𝑥_0,𝑖+1𝑥−𝑥_0𝑖+𝑥−𝑥_0,𝑖+1𝑥−𝑥_0𝑖6𝑥_0,𝑖+1−𝑥_0𝑖𝑓⁽³⁾𝜉_0𝑖𝑥−𝑥_0,𝑖+1²−𝑓⁽³⁾𝜂_0𝑖𝑥−𝑥_0𝑖²𝑓∶=^′′(𝑠)2𝑥−𝑥_0,𝑖+1𝑥−𝑥_0𝑖+ℛ¹_𝑖(𝑥)+ℛ²_𝑖(𝑥); (28) by applying Taylor expansion to 𝑓^′′(𝑥) in (28) at 𝑠, we have 𝑓^′′(𝑥)=𝑓^′′(𝑠)+𝑓⁽³⁾𝜂_𝑖(𝑥−𝑠). (29) Therefore, (25) can be obtained directly from (28) and (29).

Lemma 4. Assume that 𝑠∈(𝑥_0𝑖,𝑥_0,𝑖+1) and 𝑐_𝑖=2(𝑠−𝑥_0𝑖)/ℎ−1,1≤𝑖≤𝑛−1; there holds 𝜙^′₁𝑐_𝑖=−1⨎2ℎ^{𝑥0,𝑖+1𝑥0𝑖}𝑥−𝑥_0𝑖𝑥−𝑥_0,𝑖+1(𝑥−𝑠)²−1𝑑𝑥,𝑖=𝑚,∫2ℎ^{𝑥0,𝑖+1𝑥0𝑖}𝑥−𝑥_0𝑖𝑥−𝑥_0,𝑖+1(𝑥−𝑠)²𝑑𝑥,𝑖≠𝑚. (30)

Proof. According to (2) and the linear transformation (6) for 𝑖=𝑚, we have ^𝑥0,𝑚+1_𝑥0𝑚𝑥−𝑥_0𝑚𝑥−𝑥_0,𝑚+1(𝑥−𝑠)²𝑑𝑥=ℎlim_𝜀→0^{𝑐𝑚−(2𝜀/ℎ)}₋₁+¹_{𝑐𝑚+()2𝜀/ℎ}×𝜏²−1𝜏−𝑐_𝑚²2𝜏𝑑𝜏−²−1𝜀=ℎ¹₋₁𝜏²−1𝜏−𝑐_𝑚²𝑑𝜏=−2ℎ𝜙^′₁𝑐_𝑚, (31) where we have used 𝑥=𝑥_0𝑚(𝜏). The second identity can be similarly obtained.

Lemma 5. Suppose 𝜙^′(𝑥) and 𝜂(𝑠) are defined by (18) and (24), respectively; then one has ^∞′_𝑖=𝑚𝜙^′₁(2𝑖+𝜏)+^∞′_{𝑖=𝑛−𝑚}𝜙^′₁(−2𝑖+𝜏)≤𝐶𝜂(𝑠)ℎ, (32) where ∑^′ denotes that the first interval is certain part of the reference element.

Proof. By straightforward calculation, we have 𝜙^′₁(𝑥)≤𝐶¹₋₁𝑑𝜏|𝜏−𝑥|². (33) Noting that 𝑠=𝑥_0𝑚+((𝜏+1)/2)ℎ=𝑎+(𝑚+(𝜏/2))ℎ−(ℎ/3), we have 2(𝑠−𝑎)/ℎ=3𝜏+6𝑚−2 and ^∞′_𝑖=𝑚𝜙^′₁(2𝑖+𝜏)≤𝐶¹_2/3𝑑𝑡|2𝑖+𝜏−𝑡|²+^∞_𝑖=𝑚+1¹₋₁𝑑𝑡|2𝑖+𝜏−𝑡|²=𝐶^∞_{3𝜏+6𝑚−2}𝑑𝑥𝑥²=𝐶=3𝜏+6𝑚−2𝐶ℎ.𝑠−𝑎 (34) On the other hand, since 𝑏=𝑎+𝑛ℎ, we have 6(𝑏−𝑠)/ℎ=6(𝑛−𝑚)−(3𝜏−2) and ^∞′_{𝑖=𝑛−𝑚}𝜙^′₁(𝜏−2𝑖)≤𝐶^2/3₋₁𝑑𝑡|2𝑖−𝜏+𝑡|²+^∞_{𝑖=𝑛−𝑚+1}¹₋₁𝑑𝑡|2𝑖−𝜏+𝑡|²=𝐶^∞_{6()−()𝑛−𝑚3𝜏−2}𝑑𝑥𝑥²=𝐶6()−()=𝑛−𝑚3𝜏−2𝐶ℎ.𝑏−𝑠 (35) Combining (34) and (35), we get (32).

Set 𝐻_𝑚(𝑥)=𝑓(𝑥)−𝑓_𝐿(𝑥)−𝑓^′′(𝑠)2𝑥−𝑥_0,𝑚+1𝑥−𝑥_0𝑚. (36)

Lemma 6. Under the same assumptions of Theorem 2, for 𝐻_𝑚(𝑥) in (36), there holds that ^𝑥0,𝑚+1_𝑥0𝑚𝐻_𝑚(𝑥)(𝑥−𝑠)²𝑑𝑥≤𝐶ln𝛾(𝜏)ℎ², (37) where 𝛾(𝜏) is defined in (10).

Proof. Since 𝑓(𝑥)∈𝐶³[𝑎,𝑏], by Taylor expansion, we have 𝐻_𝑚^(𝑖)(𝑥)≤𝐶ℎ^3−𝑖,𝑖=0,1,2. (38)
By the definition of finite-part integral (13), we have ^𝑥0,𝑚+1_𝑥0𝑚𝐻_𝑚(𝑥)(𝑥−𝑠)²=𝑑𝑥ℎ𝐻_𝑚(𝑠)𝑠−𝑥_0𝑚𝑥_0,𝑚+1−𝑠+𝐻^′_𝑚(𝑠)𝑥ln^0,𝑚+1−𝑠𝑠−𝑥_0𝑚+^𝑥0,𝑚+1_𝑥0𝑚𝐻_𝑚(𝑥)−𝐻_𝑚(𝑠)−𝐻^′_𝑚(𝑠)(𝑥−𝑠)(𝑥−𝑠)²𝑑𝑥. (39) Now, we estimate the right hand side of (39) term by term. Since 𝐻_𝑚(𝑥_0𝑚)=0, we have ℎ𝐻_𝑚(𝑠)𝑠−𝑥_0𝑚𝑥_0,𝑚+1=ℎ𝐻−𝑠^𝑚(𝑠)−𝐻_𝑚𝑥_0𝑚𝑠−𝑥_0𝑚𝑥_0,𝑚+1=−𝑠ℎ𝐻^′_𝑚𝜉_𝑚𝑠−𝑥_0,𝑚+1≤𝐶ℎ²,𝜉_𝑚∈𝑠,𝑥_0,𝑚+1,𝐻^′_𝑚(𝑠)𝑥ln^0,𝑚+1−𝑠𝑠−𝑥_0𝑚≤𝐶ln𝛾(𝜏)ℎ²,_{𝑥0,𝑚+1𝑥0𝑚}𝐻_𝑚(𝑥)−𝐻_𝑚(𝑠)−𝐻^′_𝑚(𝑠)(𝑥−𝑠)(𝑥−𝑠)²=𝑑𝑥_{𝑥0,𝑚+1𝑥0𝑚}12𝐻_𝑚^′′𝜂_𝑚𝑑𝑥≤𝐶ℎ²,𝜂_𝑚∈𝑥_0𝑚,𝑥_0,𝑚+1. (40)(41)(42)
Combining (40), (41), and (42) leads to (37) and the proof is completed.

Proof of Theorem 2. Consider ^𝑥0𝑚_𝑎+^𝑏_𝑥0,𝑚+1𝑓(𝑥)(𝑥−𝑠)²𝑑𝑥−^𝑛_{𝑖=0,𝑖≠𝑚}^𝑥0,𝑖+1_0𝑖𝑓_𝐿(𝑥)(𝑥−𝑠)²=𝑑𝑥_{𝑛−1𝑖=0,𝑖≠𝑚𝑥0,𝑖+1𝑥0𝑖}𝑓(𝑥)−𝑓_𝐿(𝑥)(𝑥−𝑠)²=𝑓𝑑𝑥_′′(𝑠)2^𝑛−1_{𝑖=1,𝑖≠𝑚}^𝑥0,𝑖+1_𝑥0𝑖𝑥−𝑥_0𝑖𝑥−𝑥_0,𝑖+1(𝑥−𝑠)²+𝑑𝑥_{𝑛−1𝑖=0,𝑖≠𝑚𝑥0,𝑖+1𝑥0𝑖}ℛ¹_𝑖(𝑥)(𝑥−𝑠)²𝑑𝑥+^𝑛−1_{𝑖=0,𝑖≠𝑚}^𝑥0,𝑖+1_𝑥0𝑖ℛ²_𝑖(𝑥)(𝑥−𝑠)²+𝑑𝑥_𝑥01𝑎𝑓(𝑥)−𝑓_𝐿(𝑥)(𝑥−𝑠)²𝑑𝑥+^𝑏_𝑥0𝑛𝑓(𝑥)−𝑓_𝐿(𝑥)(𝑥−𝑠)²𝑑𝑥. (43)
By the definition of 𝐸_𝑚(𝑥), we have ^𝑥0,𝑚+1_𝑥0𝑚𝑓(𝑥)−𝑓_𝐿(𝑥)(𝑥−𝑠)²𝑑𝑥=^𝑥0,𝑚+1_𝑥0𝑚ℋ_𝑚(𝑥)(𝑥−𝑠)²+𝑓𝑑𝑥_′′(𝑠)2^𝑥0,𝑚+1_𝑥0𝑚𝑥−𝑥_0𝑖𝑥−𝑥_0,𝑖+1(𝑥−𝑠)²𝑑𝑥. (44)
Putting (43) and (44) together, we have ^𝑏_𝑎𝑓(𝑥)−𝑓_𝐿(𝑥)(𝑥−𝑠)²𝑑𝑥=^𝑛−1_𝑖=1^𝑥0,𝑖+1_𝑥0𝑖𝑓(𝑥)−𝑓_𝐿(𝑥)(𝑥−𝑠)²+𝑑𝑥_𝑥01𝑎𝑓(𝑥)−𝑓_𝐿(𝑥)(𝑥−𝑠)²+𝑑𝑥_{𝑏𝑥0𝑛}𝑓(𝑥)−𝑓_𝐿(𝑥)(𝑥−𝑠)²𝑓𝑑𝑥=−^′′(𝑠)2ℎ𝑆₁(𝜏)+ℛ_𝑓(𝑠), (45) where ℛ_𝑓(𝑠)=ℛ¹(𝑠)+ℛ²(𝑠)+ℛ³(𝑠),ℛ¹(𝑠)=^𝑥0,𝑚+1_𝑥0𝑚ℋ_𝑚(𝑥)(𝑥−𝑠)²ℛ𝑑𝑥,₂(𝑠)=^𝑥01_𝑎+^𝑛−1_{𝑖=1,𝑖≠𝑚}^𝑥0,𝑖+1_𝑥0𝑖+^𝑏_𝑥0𝑛ℛ¹_𝑖(𝑥)(𝑥−𝑠)²+𝑑𝑥_𝑥01𝑎+^𝑛−1_{𝑖=1,𝑖≠𝑚}^𝑥0,𝑖+1_𝑥0𝑖+^𝑏_𝑥0𝑛ℛ²_𝑖(𝑥)(𝑥−𝑠)²ℛ𝑑𝑥,₃(𝑠)𝑓=ℎ^′′(𝑠)2^∞′_𝑖=𝑚𝜙^′₁(2𝑖+𝜏)+^∞′_{𝑖=𝑛−𝑚}𝜙^′₁(−2𝑖+𝜏). (46) For the first part of ℛ²(𝑠), ^𝑥01_𝑎+^𝑛−1_{𝑖=1,𝑖≠𝑚}^𝑥0,𝑖+1_𝑥0𝑖+^𝑏_𝑥0𝑛ℛ¹_𝑖(𝑥)(𝑥−𝑠)²≤𝑑𝑥_𝑥01𝑎+^𝑛−1_{𝑖=1,𝑖≠𝑚−1,𝑚,𝑚+1}^𝑥0,𝑖+1_𝑥0𝑖+^𝑏_𝑥0𝑛ℛ¹_𝑖(𝑥)(𝑥−𝑠)²+𝑑𝑥_{𝑥0,𝑚𝑥0,𝑚−1}ℛ¹_𝑚−1(𝑥)(𝑥−𝑠)²𝑑𝑥+^𝑥0,𝑚+2_𝑥0,𝑚+1ℛ¹_𝑚+1(𝑥)(𝑥−𝑠)²≤𝑑𝑥_𝑥01𝑎+^𝑛−1_{𝑖=1,𝑖≠𝑚−1,𝑚,𝑚+1}^𝑥0,𝑖+1_𝑥0𝑖+^𝑏_𝑥0𝑛ℛ¹_𝑖(𝑥)(𝑥−𝑠)²+𝑑𝑥_{𝑥0𝑚𝑥0,𝑚−1}ℛ¹_𝑚−1(𝑥)−ℛ¹_𝑚−1𝑥_0𝑚(𝑥−𝑠)²+𝑑𝑥_{𝑥0,𝑚+2𝑥0,𝑚+1}ℛ¹_𝑚+1(𝑥)−ℛ¹_𝑚+1𝑥_0,𝑚+1(𝑥−𝑠)²=𝑑𝑥_𝑥01𝑎+^𝑛−1_{𝑖=1,𝑖≠𝑚−1,𝑚,𝑚+1}^𝑥0,𝑖+1_𝑥0𝑖+^𝑏_𝑥0𝑛ℛ¹_𝑖(𝑥)(𝑥−𝑠)²+𝑑𝑥_{𝑥0𝑚𝑥0,𝑚−1}ℛ^1′_𝑚−1𝜉_0𝑚𝑥−𝑥_0𝑚(𝑥−𝑠)²+𝑑𝑥_{𝑥0,𝑚+2𝑥0,𝑚+1}ℛ^1′_𝑚+1𝜉_0,𝑚+1𝑥−𝑥_0,𝑚+1(𝑥−𝑠)²𝑑𝑥≤𝐶ℎ^3𝑥01_𝑎+^𝑛−1_{𝑖=1,𝑖≠𝑚−1,𝑚,𝑚+1}^𝑥0,𝑖+1_𝑥0𝑖+^𝑏_𝑥0𝑛1(𝑥−𝑠)²𝑑𝑥+𝐶ℎ^2𝑥0𝑚_{𝑥0,𝑚−1}𝑑𝑥||+𝑥−𝑠^𝑥0,𝑚+2_𝑥0,𝑚+1𝑑𝑥||(𝜏+|ℎ𝑥−𝑠≤𝐶ln𝛾)lnℎ|², (47) where 𝜉_0𝑚∈(𝑥_0,𝑚−1,𝑥_0𝑚) and 𝜉_0,𝑚+1∈(𝑥_0,𝑚,𝑥_0,𝑚+1). We have also used the identity ℛ¹_𝑚−1(𝑥_0𝑚) = ℛ¹_𝑚+1(𝑥_0,𝑚+1)=0, |𝑥−𝑥_0𝑚|<|𝑥−𝑠|, and |𝑥−𝑥_0,𝑚+1|<|𝑥−𝑠|.
For the second part of ℛ²(𝑠), ^𝑥01_𝑎+^𝑛−1_{𝑖=1,𝑖≠𝑚}^𝑥0,𝑖+1_𝑥0𝑖+^𝑏_𝑥0𝑛ℛ²_𝑖(𝑥)(𝑥−𝑠)²𝑑𝑥≤𝐶ℎ^2𝑥01_𝑎+^𝑛−1_{𝑖=1,𝑖≠𝑚}^𝑥0,𝑖+1_𝑥0𝑖+^𝑏_𝑥0𝑛1||(𝜏+|ℎ𝑥−𝑠𝑑𝑥≤𝐶ln𝛾)lnℎ|². (48) By Lemmas 5 and 6, we have 𝑅_𝑓(𝑠)≤ℛ¹(𝑠)+ℛ²(𝑠𝜂(𝑠)+||+(𝜏ℎ)≤𝐶lnℎln𝛾)². (49) Then the proof is completed.

3.2. The Calculation of 𝑆₁(𝜏)

Let 𝑄_𝑛(𝑥) be the function of the second kind associated with the Legendre polynomial 𝑃_𝑛(𝑥), defined by (cf. [24]) 𝑄₀(𝑥)=12ln𝑥+1𝑥−1,𝑄₁(𝑥)=𝑥𝑄₀(𝑥)−1. (50) We also define 𝑊𝑓,𝜏∶=𝑓(𝜏)+^∞_𝑖=0𝑓()(,().2𝑖+𝜏+𝑓−2𝑖+𝜏)𝜏∈−1,1 (51) Then, by the definition of 𝑊, 𝑊𝑄₀(𝜏)=12ln1+𝜏+11−𝜏2^∞_𝑖=1ln2𝑖+1+𝜏2𝑖−1+𝜏+ln2𝑖−1−𝜏=12𝑖+1−𝜏2lim_𝑖→∞ln2𝑖+1+𝜏𝑊2𝑖+1−𝜏=0,𝑥𝑄^′₀(𝜏)=𝜏1−𝜏²−^∞_𝑖=12𝑖+𝜏(2𝑖+𝜏)²+−1−2𝑖+𝜏(−2𝑖+𝜏)²=1−12lim_𝑛→∞^𝑘=𝑛_𝑘=−𝑛1()−()=𝜋𝑘+1/2𝜏/22tan𝜋𝜏2. (52) It follows that 𝑆𝑄^′₁,𝜏=𝑊𝑄₀+𝑥𝑄^′₀=𝜋,𝜏2tan𝜏𝜋2, (53) which means 𝑆𝑄₁=𝜋,𝜏2tan𝜏𝜋2𝜏=−lncos𝜋𝜏2+𝐶. (54)

What remains is to determine the constant 𝐶. By using the identities (cf. [24, Chapter 1, Section 1.2]), 𝑥cot𝑥=1+^∞_𝑘=1(−1)^𝑘𝐵_2𝑘(2𝑥)^2𝑘()!,()2𝑘ln2cos𝑥=−^∞_𝑗=11𝑗cos2𝑗𝑥,𝑥∈(0,𝜋), (55) where 𝐵_2𝑘 denote the Bernoulli numbers, we have ^∞_𝑘=1(−1)^𝑘𝐵_2𝑘(2𝑥)^2𝑘+1()!()()2𝑘+1=2𝑥lnsin𝑥+2ln2−1𝑥+^∞_𝑗=112𝑗²sin2𝑗𝑥. (56) Setting 𝑥=𝜋/2 gives ^∞_𝑘=1(−1)^𝑘𝐵_2𝑘(𝜋)^2𝑘+1(2𝑘+1)!=ln2−1. (57) Then we have 𝑊𝑄₁;0=−1+2^∞_𝑖=1𝑄₁(2𝑖)=−1+2^∞∞_𝑖=1𝑘=11(2𝑘+1)(2𝑖)^2𝑘=−1+^∞_𝑘=1(−1)^𝑘+1𝐵_2𝑘(𝜋)^2𝑘+1(2𝑘+1)!=−ln2, (58) where we have used the formulae (cf. [24, Chapter 1, Section 1.2]) 𝑄₁(𝑥)=^∞_𝑘=11(2𝑘+1)𝑥^2𝑘,|𝑥|>1,_∞𝑖=11𝑖^2𝑘=(−1)^𝑘+12^2𝑘−1()!𝐵2𝑘^2𝑘𝜋^2𝑘. (59) Then we have 𝑆₁(𝜏)=−ln2cos𝜏𝜋2. (60)

By Theorem 2, we get the following error expansion: 𝐸_𝑛𝑓=𝑓^′′(𝑠)ℎ2log2cos𝜋𝜏2+𝑂ℎ². (61)

Let 𝑇ℎ_𝑙=𝐼₂^𝑙−1_𝑛0𝑓,𝑠,𝑙=1,2,…, (62) where 𝑛₀ is the starting meshes, 𝑙 is the refining numbers, and ℎ_𝑙=(𝑏−𝑎)/2^𝑙−1𝑛₀. Then we have the following.

Corollary 7. Under the same assumption of Theorem 2 and (62), there holds 𝑇ℎ_𝑙−𝑇ℎ_𝑙−1=𝑂ℎ²_𝑙−1. (63)

4. Numerical Examples

In this section, computational results are reported to confirm our theoretical analysis.

Example 1. Consider the hypersingular integral ¹₀𝑥³(𝑥−𝑠)²3𝑑𝑥=3𝑠+2+1𝑠−1+3𝑠²log1−𝑠𝑠. (64)

We examine the dynamic point 𝑠=𝑥_[0𝑛/4]+(1+𝜏)ℎ/2, in Table 1 show that when the local coordinate of singular point 𝜏=±2/3, the quadrature reach the convergence rate of 𝑂(ℎ²) as for the nonsupersingular point, there are no convergence rate which agree with our theorematically analysis. For the case of 𝑠=𝑏−(𝜏+1)ℎ/2, Table 2 shows that there is no superconvergence phenomenon because of the influence of 𝜂(𝑠) which coincides with our theoretical analysis.

Table 1 

Errors of the mod-trapezoidal rule with 𝑠=𝑥_[0𝑛/4]+(1+𝜏)ℎ/2.

Table 2 

Errors of the mod-trapezoidal rule with 𝑠=𝑥_0𝑛−(1+𝜏)ℎ/2.

Example 2. Consider the hypersingular integral ¹₀𝑥⁴+1(𝑥−𝑠)²𝑑𝑥=4𝑠²4+2𝑠+3+𝑠+1𝑠(𝑠−1)+4𝑠³log1−𝑠𝑠. (65)

The numerical results show that the convergence rate reaches (ℎ²) when the singular point coincides with the mesh point in Table 3. In Table 4, a posteriori error estimate is presented and the convergence rate is also 𝑂(ℎ²) which agrees with our theoretical analysis.

Table 3 

Errors of the mod-trapezoidal rule.

Table 4 

A posteriori error estimate of the mod-trapezoidal rule.

Example 3. Now we consider an example of less regularity. Let 𝑎=−𝑏=−1, 𝑠=0, and 𝑓(𝑥)=ℱ_𝑖(𝑥)∶=𝑥²+2+sign(𝑥)|𝑥|^2−𝑖+0.5,𝑖=−1,0,1. (66) Obviously, ℱ_𝑖(𝑥)∈𝐶^3−𝑖+0.5[−1,1](𝑖=−1,0,1). The exact value of the integral is ℐ₂ℱ_𝑖(𝑥),0=14−4𝑖.3−2𝑖 (67)

The numerical results are presented in Tables 5 and 6. When the density function 𝑓(𝑥) is smooth enough 𝑖=−1, the error bound is 𝑂(ℎ²), and if the density function has less regularity 𝑖=0,1, there is no hyperconvergence phenomenon, which means the regularity of density function cannot be reduced.

Table 5 

Errors of the mod-trapezoidal rule.

Table 6 

A posteriori error estimate of the mod-trapezoidal rule.

Acknowledgments

The work of Jin Li was supported by the National Natural Science Foundation of China (nos. 11101247, 11201209, and 11101246), China Postdoctoral Science Foundation Fund Project (no. 2013M540541), the Shandong Provincial Natural Science Foundation of China (no. ZR2011AQ020), and a Project of Shandong Province Higher Educational Science and Technology Program (no. J11LE8).