# The Modified Trapezoidal Rule for Computing Hypersingular Integral on Interval

Authors: Jin Li; Xiuzhen Li
Journal: Journal of Applied Mathematics
 (2013)

---

## 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.


---

## Body

## 1. Introduction

Consider the following integral:(1)I(f,s):=⨎abf(x)(x-s)p+1dx,s∈(a,b),p=1,2,
where ⨎ab denotes a Hadamard finite-part integral (p=1 is called hypersingular integral and p=2 is called supersingular integral) and s 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:(2)⨎abf(x)(x-s)2dx=limε→0{∫as-εf(x)(x-s)2dx+∫s+εb‍f(x)(x-s)2dx-2f(s)ε},s∈(a,b).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 s is close to a mesh point. In order to make the mesh be selected in such a way that s 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 O(h|lnh|). 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 O(h) 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 Section2, 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

Leta=x0<x1<⋯<xn-1<xn=b be a uniform partition of the interval [a,b] with mesh size h=(b-a)/n and set
(3)x0i=xi-1+h6,i=1,2,…,n;
then we get the new partition:
(4)a=x00<x01<⋯<x0n<x0,n+1=b.We definefL(x), the modified trapezoidal interpolation for f(x), as
(5)fL(x)=x-x0ix0,i+1-x0if(x0,i+1)+x0,i+1-xx0,i+1-x0if(x0i),x∈[x0i,x0,i+1],0≤i≤n,
and a linear transformation
(6)x=x^0i(τ):=(τ+1)x0,i+1-x0i2+x0i,i=1,…,n-1,τ∈[-1,1],
from the reference element [-1,1] to the subinterval [x0i,x0,i+1]. For the two subintervals [a,x01] and [x0n,b] near the end of the interval [a,b], τ values in [2/3,1] and [-1,2/3], respectively.Replacingf(x) in (1) with fL(x) gives the new composite trapezoidal rule:
(7)In(f,s):=⨎abfL(x)(x-s)2dx=∑i=0n+1ωj(s)f(x0i)=I(f,s)-En(f),
where ωi(s) is the Cotes coefficients:
(8)ωi(s)=1-δi0x0,i-x0,i-1ln|x0i-sx0,i-1-s|-1-δi,n+1x0,i+1-x0iln|x0,i+1-sx0i-s|+δi0x00-s+δi,n+1s-x0,n+1,0≤i≤n+1, δij denotes the Kronecker delta, and En(f) denotes the error functional.Theorem 1.
Assumef(x)∈C1+α[a,b], α∈[0,1). For the trapezoidal rule In(f,s) defined in (7), there exists a positive constant C, independent of h and s, such that
(9)|En(f)|≤C(|lnh|+|lnγ(τ)|)hα,
where
(10)γ(τ)=min0≤i≤n+1|s-xi|h=1-|τ|2,τ∈(-1,1).Proof.
LetR(x)=f(x)-fL(x); then we have |R(x)|≤Ch1+α, as
(11)I(f,s)-In(f,s)=⨎abf(x)-fL(x)(x-s)2dx=⨎x0mx0,m+1f(x)-fL(x)(x-s)2dx+∑i=0,i≠mn∫x0ix0,i+1f(x)-fL(x)(x-s)2dx.
For the first part of (11), since R(x)∈C1+α[a,b], by Taylor expansion, we have
(12)|R(x)|≤Ch1+α-i,i=0,1,2.
By the definition of finite-part integral,(13)⨎abf(x)(x-s)2dx=(b-a)f(s)(s-a)(b-s)+f′(s)lnb-ss-a+∫abf(x)-f(s)-f′(s)(x-s)(x-s)2dx,
we have
(14)⨎x0mx0,m+1R(x)(x-s)2dx=hR(s)(s-x0m)(x0,m+1-s)+R′(s)lnx0,m+1-ss-x0m+∫x0mx0,m+1R(x)-R(s)-R′(s)(x-s)(x-s)2dx.
Now, we estimate the right hand side of (14) term by term. Since R(x0m)=0, we have
(15)|hR(s)(s-x0m)(x0,m+1-s)|=|h[R(s)-R(x0m)](s-x0m)(x0,m+1-s)|=|hRm′(ξm)(s-x0,m+1)|≤Chα,ξm∈(x0m,s),|R′(s)lnx0,m+1-ss-x0m|≤C|lnγ(τ)|hα,|∫x0mx0,m+1R(x)-R(s)-R′(s)(x-s)(x-s)2dx|≤Chα,ηm∈(x0m,x0,m+1).
For the second part of (11),
(16)|(∫ax01‍+∑i=1,i≠mn-1‍∫x0ix0,i+1‍+∫x0nb)R(x)(x-s)2dx|≤|(∫ax01+∑i=1,i≠m-1,m,m+1n-1∫x0ix0,i+1+∫x0nb∑i=1,i≠m-1,m,m+1n-1∫x0ix0,i+1)R(x)(x-s)2dx|+|∫x0,m-1x0,mR(x)(x-s)2dx+∫x0,m+1x0,m+2R(x)(x-s)2dx|≤|(∫ax01+∑i=1,i≠m-1,m,m+1n-1∫x0ix0,i+1+∫x0nb∑i=1,i≠m-1,m,m+1n-1∫x0ix0,i+1)R(x)(x-s)2dx|+|∫x0,m-1x0mR(x)-R(x0m)(x-s)2dx+∫x0,m+1x0,m+2R(x)-R(x0,m+1)(x-s)2dx|=|(∫ax01+∑i=1,i≠m-1,m,m+1n-1∫x0ix0,i+1+∫x0nb)R(x)(x-s)2dx|+|∫x0,m-1x0mR′(ξ0m)(x-x0m)(x-s)2dx+∫x0,m+1x0,m+2R′(ξ0,m+1)(x-x0,m+1)(x-s)2dx|≤Ch1+α(∫ax01+∑i=1,i≠m-1,m,m+1n-1∫x0ix0,i+1+∫x0nb)1(x-s)2dx+Chα(∫x0,m-1x0mdx|x-s|+∫x0,m+1x0,m+2dx|x-s|)≤C[|lnγ(τ)|+|lnh|]hα.
Combining (14), (15), and (16) leads to (9) and the proof is completed.Firstly, we set(17)ϕ1(x)={-12⨍-11τ2-1τ-xdτ,|x|<1,-12∫-11τ2-1τ-xdτ,|x|>1;
then we have
(18)ϕ1′(x)={-12⨎-11τ2-1(τ-x)2dτ,|x|<1,-12∫-11τ2-1(τ-x)2dτ,|x|>1.
By straight calculation, we get
(19)ϕ1′(x)=xlog1+x1-x-2.
Let J:=(-∞,-1)⋃(-1,1)⋃(1,+∞) and the operator W:C(J)→C(-1,1) be defined as
(20)Wf(τ):=f(τ)+∑i=1∞[f(2i+τ)+f(-2i+τ)],τ∈(-1,1).
Obviously, the operator W is linear operator. Then we set
(21)S1(τ)=Wϕ1′(τ)=ϕ1′(τ)+∑i=1∞[ϕ1′(2i+τ)+ϕ1′(-2i+τ)].Now we present our main results below. The proof will be given in the next section.Theorem 2.
Assumef(x)∈C3[a,b]. For the trapezoidal rule In(f,s) defined in (7), there exists a positive constant C, independent of h and s, such that
(22)En(f)=f′′(s)h2S1(τ)+ℛf(s),
where s=x0i+(1+τ)h/2, i=1,2,…,n, and
(23)|ℛf(s)|≤C(η(s)+|lnh|+|lnγ(τ)|)h2
(γ(τ) defined as (10)) and
(24)η(s)=max{1s-a,1b-s}.


## 3. Proof of Main Results

### 3.1. Preliminaries

In the following section,C denotes certain constant independent of h and s, and its value varies with places.Lemma 3.
Assume thatf(x)∈C3[a,b] and fL(x) are defined by (5); there holds
(25)f(x)-fL(x)=f′′(s)2(x-x0,i+1)(x-x0i)+ℛi1(x)+ℛi2(x),
where
(26)ℛi1(x)=(x-x0,i+1)(x-x0i)6(x0,i+1-x0i)[f(3)(ξ0i)(x-x0,i+1)2hhhhhhhhhhhhhhhhhh-f(3)(η0i)(x-x0i)2(x-x0,i+1)2],ℛi2(x)=f(3)(β0i)2(x-x0,i+1)(x-x0i)(x-s),ξ0i,η0i,β0i∈(x0i,x0,i+1), and
(27)|ℛi1(x)|≤Ch3.Proof.
Taking Taylor expansion forf(x0i), f(x0,i+1) at x, there holds
(28)f(x)-fL(x)=f′′(x)2!(x-x0,i+1)(x-x0i)+(x-x0,i+1)(x-x0i)6(x0,i+1-x0i)[f(3)(ξ0i)(x-x0,i+1)2-f(3)(η0i)(x-x0i)2(x-x0,i+1)2]:=f′′(s)2(x-x0,i+1)(x-x0i)+ℛi1(x)+ℛi2(x);
by applying Taylor expansion to f′′(x) in (28) at s, we have
(29)f′′(x)=f′′(s)+f(3)(ηi)(x-s).
Therefore, (25) can be obtained directly from (28) and (29).Lemma 4.
Assume thats∈(x0i,x0,i+1) and ci=2(s-x0i)/h-1,1≤i≤n-1; there holds
(30)ϕ1′(ci)={-12h⨎x0ix0,i+1(x-x0i)(x-x0,i+1)(x-s)2dx,i=m,-12h∫x0ix0,i+1(x-x0i)(x-x0,i+1)(x-s)2dx,i≠m.Proof.
According to (2) and the linear transformation (6) for i=m, we have
(31)⨎x0mx0,m+1(x-x0m)(x-x0,m+1)(x-s)2dx=hlimε→0{τ2-1(τ-cm)2dτ-2(τ2-1)ε(∫-1cm-(2ε/h)+∫cm+(2ε/h)1)×τ2-1(τ-cm)2dτ-2(τ2-1)ε}=h⨎-11τ2-1(τ-cm)2dτ=-2hϕ1′(cm),
where we have used x=x^0m(τ). The second identity can be similarly obtained.Lemma 5.
Supposeϕ′(x) and η(s) are defined by (18) and (24), respectively; then one has
(32)|∑i=m∞′ϕ1′(2i+τ)+∑i=n-m∞′ϕ1′(-2i+τ)|≤Cη(s)h,
where ∑′ denotes that the first interval is certain part of the reference element.Proof.
By straightforward calculation, we have(33)|ϕ1′(x)|≤C∫-11dτ|τ-x|2.
Noting that s=x0m+((τ+1)/2)h=a+(m+(τ/2))h-(h/3), we have 2(s-a)/h=3τ+6m-2 and
(34)|∑i=m∞′ϕ1′(2i+τ)|≤C[∫2/31‍dt|2i+τ-t|2+∑i=m+1∞∫-11‍dt|2i+τ-t|2]=C∫3τ+6m-2∞dxx2=C3τ+6m-2=Chs-a.
On the other hand, since b=a+nh, we have 6(b-s)/h=6(n-m)-(3τ-2) and
(35)|∑i=n-m∞′ϕ1′(τ-2i)|≤C[∫-12/3dt|2i-τ+t|2+∑i=n-m+1∞∫-11dt|2i-τ+t|2]=C∫6(n-m)-(3τ-2)∞dxx2=C6(n-m)-(3τ-2)=Chb-s.
Combining (34) and (35), we get (32).Set(36)Hm(x)=f(x)-fL(x)-f′′(s)2(x-x0,m+1)(x-x0m).Lemma 6.
Under the same assumptions of Theorem2, for Hm(x) in (36), there holds that
(37)|⨎x0mx0,m+1Hm(x)(x-s)2dx|≤C|lnγ(τ)|h2,
where γ(τ) is defined in (10).Proof.
Sincef(x)∈C3[a,b], by Taylor expansion, we have
(38)|Hm(i)(x)|≤Ch3-i,i=0,1,2.
By the definition of finite-part integral (13), we have
(39)⨎x0mx0,m+1Hm(x)(x-s)2dx=hHm(s)(s-x0m)(x0,m+1-s)+Hm′(s)lnx0,m+1-ss-x0m+∫x0mx0,m+1Hm(x)-Hm(s)-Hm′(s)(x-s)(x-s)2dx.
Now, we estimate the right hand side of (39) term by term. Since Hm(x0m)=0, we have
(40)|hHm(s)(s-x0m)(x0,m+1-s)|=|h[Hm(s)-Hm(x0m)](s-x0m)(x0,m+1-s)|=|hHm′(ξm)(s-x0,m+1)|≤Ch2,ξm∈(s,x0,m+1),(41)|Hm′(s)lnx0,m+1-ss-x0m|≤C|lnγ(τ)|h2,(42)|∫x0mx0,m+1Hm(x)-Hm(s)-Hm′(s)(x-s)(x-s)2dx|=|∫x0mx0,m+112Hm′′(ηm)dx|≤Ch2,ηm∈(x0m,x0,m+1).
Combining (40), (41), and (42) leads to (37) and the proof is completed.Proof of Theorem2.
Consider(43)(∫ax0m+∫x0,m+1b)f(x)(x-s)2dx-∑i=0,i≠mn∫0ix0,i+1fL(x)(x-s)2dx=∑i=0,i≠mn-1∫x0ix0,i+1f(x)-fL(x)(x-s)2dx=f′′(s)2∑i=1,i≠mn-1∫x0ix0,i+1(x-x0i)(x-x0,i+1)(x-s)2dx+∑i=0,i≠mn-1∫x0ix0,i+1ℛi1(x)(x-s)2dx+∑i=0,i≠mn-1∫x0ix0,i+1ℛi2(x)(x-s)2dx+∫ax01f(x)-fL(x)(x-s)2dx+∫x0nbf(x)-fL(x)(x-s)2dx.
By the definition ofEm(x), we have
(44)⨎x0mx0,m+1f(x)-fL(x)(x-s)2dx=⨎x0mx0,m+1ℋm(x)(x-s)2dx+f′′(s)2⨎x0mx0,m+1(x-x0i)(x-x0,i+1)(x-s)2dx.
Putting (43) and (44) together, we have
(45)⨎abf(x)-fL(x)(x-s)2dx=∑i=1n-1⨎x0ix0,i+1f(x)-fL(x)(x-s)2dx+∫ax01f(x)-fL(x)(x-s)2dx+∫x0nbf(x)-fL(x)(x-s)2dx=-f′′(s)2hS1(τ)+ℛf(s),
where
(46)ℛf(s)=ℛ1(s)+ℛ2(s)+ℛ3(s),ℛ1(s)=⨎x0mx0,m+1ℋm(x)(x-s)2dx,ℛ2(s)=(∫ax01+∑i=1,i≠mn-1∫x0ix0,i+1+∫x0nb)ℛi1(x)(x-s)2dx+(∫ax01+∑i=1,i≠mn-1∫x0ix0,i+1+∫x0nb)ℛi2(x)(x-s)2dx,ℛ3(s)=hf′′(s)2[∑i=m∞′ϕ1′(2i+τ)+∑i=n-m∞′ϕ1′(-2i+τ)].
For the first part of ℛ2(s),
(47)|(∫ax01+∑i=1,i≠mn-1∫x0ix0,i+1+∫x0nb)ℛi1(x)(x-s)2dx|≤|(∫ax01+∑i=1,i≠m-1,m,m+1n-1∫x0ix0,i+1+∫x0nb)ℛi1(x)(x-s)2dx|+|∫x0,m-1x0,mℛm-11(x)(x-s)2dx+∫x0,m+1x0,m+2ℛm+11(x)(x-s)2dx|≤|(∫ax01+∑i=1,i≠m-1,m,m+1n-1∫x0ix0,i+1+∫x0nb)ℛi1(x)(x-s)2dx|+|∫x0,m-1x0mℛm-11(x)-ℛm-11(x0m)(x-s)2dx+∫x0,m+1x0,m+2ℛm+11(x)-ℛm+11(x0,m+1)(x-s)2dx|=|(∫ax01+∑i=1,i≠m-1,m,m+1n-1∫x0ix0,i+1+∫x0nb)ℛi1(x)(x-s)2dx|+|∫x0,m-1x0mℛm-11′(ξ0m)(x-x0m)(x-s)2dx+∫x0,m+1x0,m+2ℛm+11′(ξ0,m+1)(x-x0,m+1)(x-s)2dx|≤Ch3(∫ax01+∑i=1,i≠m-1,m,m+1n-1∫x0ix0,i+1+∫x0nb)1(x-s)2dx+Ch2(∫x0,m-1x0mdx|x-s|+∫x0,m+1x0,m+2dx|x-s|)≤C[|lnγ(τ)|+|lnh|]h2,
where ξ0m∈(x0,m-1,x0m) and ξ0,m+1∈(x0,m,x0,m+1). We have also used the identity ℛm-11(x0m) = ℛm+11(x0,m+1)=0, |x-x0m|<|x-s|, and |x-x0,m+1|<|x-s|.
For the second part ofℛ2(s),
(48)|(∫ax01+∑i=1,i≠mn-1∫x0ix0,i+1+∫x0nb)ℛi2(x)(x-s)2dx|≤Ch2(∫ax01+∑i=1,i≠mn-1∫x0ix0,i+1+∫x0nb)1|x-s|dx≤C[|lnγ(τ)|+|lnh|]h2.
By Lemmas 5 and 6, we have
(49)|Rf(s)|≤|ℛ1(s)|+|ℛ2(s)|≤C[η(s)+|lnh|+|lnγ(τ)|]h2.
Then the proof is completed.


### 3.2. The Calculation ofS1(τ)

LetQn(x) be the function of the second kind associated with the Legendre polynomial Pn(x), defined by (cf. [24])
(50)Q0(x)=12ln|x+1x-1|,Q1(x)=xQ0(x)-1.
We also define
(51)W(f,τ):=f(τ)+∑i=0∞[f(2i+τ)+f(-2i+τ)],τ∈(-1,1).
Then, by the definition of W,
(52)W(Q0)(τ)=12ln1+τ1-τ+12∑i=1∞(ln2i+1+τ2i-1+τ+ln2i-1-τ2i+1-τ)=12limi→∞ln2i+1+τ2i+1-τ=0,W(xQ0′)(τ)=τ1-τ2-∑i=1∞(2i+τ(2i+τ)2-1+-2i+τ(-2i+τ)2-1)=12limn→∞∑k=-nk=n1k+(1/2)-(τ/2)=π2tanπτ2.
It follows that
(53)S(Q1′,τ)=W(Q0+xQ0′,τ)=π2tanτπ2,
which means
(54)S(Q1,τ)=∫‍π2tanτπ2τ=-lncosπτ2+C.What remains is to determine the constantC. By using the identities (cf. [24, Chapter 1, Section 1.2]),
(55)xcotx=1+∑k=1∞(-1)kB2k(2x)2k(2k)!,ln(2cosx)=-∑j=1∞1jcos(2jx),x∈(0,π),
where B2k denote the Bernoulli numbers, we have
(56)∑k=1∞(-1)kB2k(2x)2k+1(2k+1)!=2xln(sinx)+2[(ln2-1)x+∑j=1∞12j2sin(2jx)].
Setting x=π/2 gives
(57)∑k=1∞(-1)kB2k(π)2k+1(2k+1)!=ln2-1.
Then we have
(58)W(Q1;0)=-1+2∑i=1∞Q1(2i)=-1+2∑i=1∞∑k=1∞1(2k+1)(2i)2k=-1+∑k=1∞(-1)k+1B2k(π)2k+1(2k+1)!=-ln2,
where we have used the formulae (cf. [24, Chapter 1, Section 1.2])
(59)Q1(x)=∑k=1∞1(2k+1)x2k,|x|>1,∑i=1∞1i2k=(-1)k+122k-1(2k)!B2kπ2k.
Then we have
(60)S1(τ)=-ln2cosτπ2.By Theorem2, we get the following error expansion:
(61)En(f)=f′′(s)h2log2cosπτ2+O(h2).Let(62)T(hl)=I2l-1n0(f,s),l=1,2,…,
where n0 is the starting meshes, l is the refining numbers, and hl=(b-a)/2l-1n0. Then we have the following.Corollary 7.
Under the same assumption of Theorem2 and (62), there holds
(63)T(hl)-T(hl-1)=O(hl-12).


## 3.1. Preliminaries

In the following section,C denotes certain constant independent of h and s, and its value varies with places.Lemma 3.
Assume thatf(x)∈C3[a,b] and fL(x) are defined by (5); there holds
(25)f(x)-fL(x)=f′′(s)2(x-x0,i+1)(x-x0i)+ℛi1(x)+ℛi2(x),
where
(26)ℛi1(x)=(x-x0,i+1)(x-x0i)6(x0,i+1-x0i)[f(3)(ξ0i)(x-x0,i+1)2hhhhhhhhhhhhhhhhhh-f(3)(η0i)(x-x0i)2(x-x0,i+1)2],ℛi2(x)=f(3)(β0i)2(x-x0,i+1)(x-x0i)(x-s),ξ0i,η0i,β0i∈(x0i,x0,i+1), and
(27)|ℛi1(x)|≤Ch3.Proof.
Taking Taylor expansion forf(x0i), f(x0,i+1) at x, there holds
(28)f(x)-fL(x)=f′′(x)2!(x-x0,i+1)(x-x0i)+(x-x0,i+1)(x-x0i)6(x0,i+1-x0i)[f(3)(ξ0i)(x-x0,i+1)2-f(3)(η0i)(x-x0i)2(x-x0,i+1)2]:=f′′(s)2(x-x0,i+1)(x-x0i)+ℛi1(x)+ℛi2(x);
by applying Taylor expansion to f′′(x) in (28) at s, we have
(29)f′′(x)=f′′(s)+f(3)(ηi)(x-s).
Therefore, (25) can be obtained directly from (28) and (29).Lemma 4.
Assume thats∈(x0i,x0,i+1) and ci=2(s-x0i)/h-1,1≤i≤n-1; there holds
(30)ϕ1′(ci)={-12h⨎x0ix0,i+1(x-x0i)(x-x0,i+1)(x-s)2dx,i=m,-12h∫x0ix0,i+1(x-x0i)(x-x0,i+1)(x-s)2dx,i≠m.Proof.
According to (2) and the linear transformation (6) for i=m, we have
(31)⨎x0mx0,m+1(x-x0m)(x-x0,m+1)(x-s)2dx=hlimε→0{τ2-1(τ-cm)2dτ-2(τ2-1)ε(∫-1cm-(2ε/h)+∫cm+(2ε/h)1)×τ2-1(τ-cm)2dτ-2(τ2-1)ε}=h⨎-11τ2-1(τ-cm)2dτ=-2hϕ1′(cm),
where we have used x=x^0m(τ). The second identity can be similarly obtained.Lemma 5.
Supposeϕ′(x) and η(s) are defined by (18) and (24), respectively; then one has
(32)|∑i=m∞′ϕ1′(2i+τ)+∑i=n-m∞′ϕ1′(-2i+τ)|≤Cη(s)h,
where ∑′ denotes that the first interval is certain part of the reference element.Proof.
By straightforward calculation, we have(33)|ϕ1′(x)|≤C∫-11dτ|τ-x|2.
Noting that s=x0m+((τ+1)/2)h=a+(m+(τ/2))h-(h/3), we have 2(s-a)/h=3τ+6m-2 and
(34)|∑i=m∞′ϕ1′(2i+τ)|≤C[∫2/31‍dt|2i+τ-t|2+∑i=m+1∞∫-11‍dt|2i+τ-t|2]=C∫3τ+6m-2∞dxx2=C3τ+6m-2=Chs-a.
On the other hand, since b=a+nh, we have 6(b-s)/h=6(n-m)-(3τ-2) and
(35)|∑i=n-m∞′ϕ1′(τ-2i)|≤C[∫-12/3dt|2i-τ+t|2+∑i=n-m+1∞∫-11dt|2i-τ+t|2]=C∫6(n-m)-(3τ-2)∞dxx2=C6(n-m)-(3τ-2)=Chb-s.
Combining (34) and (35), we get (32).Set(36)Hm(x)=f(x)-fL(x)-f′′(s)2(x-x0,m+1)(x-x0m).Lemma 6.
Under the same assumptions of Theorem2, for Hm(x) in (36), there holds that
(37)|⨎x0mx0,m+1Hm(x)(x-s)2dx|≤C|lnγ(τ)|h2,
where γ(τ) is defined in (10).Proof.
Sincef(x)∈C3[a,b], by Taylor expansion, we have
(38)|Hm(i)(x)|≤Ch3-i,i=0,1,2.
By the definition of finite-part integral (13), we have
(39)⨎x0mx0,m+1Hm(x)(x-s)2dx=hHm(s)(s-x0m)(x0,m+1-s)+Hm′(s)lnx0,m+1-ss-x0m+∫x0mx0,m+1Hm(x)-Hm(s)-Hm′(s)(x-s)(x-s)2dx.
Now, we estimate the right hand side of (39) term by term. Since Hm(x0m)=0, we have
(40)|hHm(s)(s-x0m)(x0,m+1-s)|=|h[Hm(s)-Hm(x0m)](s-x0m)(x0,m+1-s)|=|hHm′(ξm)(s-x0,m+1)|≤Ch2,ξm∈(s,x0,m+1),(41)|Hm′(s)lnx0,m+1-ss-x0m|≤C|lnγ(τ)|h2,(42)|∫x0mx0,m+1Hm(x)-Hm(s)-Hm′(s)(x-s)(x-s)2dx|=|∫x0mx0,m+112Hm′′(ηm)dx|≤Ch2,ηm∈(x0m,x0,m+1).
Combining (40), (41), and (42) leads to (37) and the proof is completed.Proof of Theorem2.
Consider(43)(∫ax0m+∫x0,m+1b)f(x)(x-s)2dx-∑i=0,i≠mn∫0ix0,i+1fL(x)(x-s)2dx=∑i=0,i≠mn-1∫x0ix0,i+1f(x)-fL(x)(x-s)2dx=f′′(s)2∑i=1,i≠mn-1∫x0ix0,i+1(x-x0i)(x-x0,i+1)(x-s)2dx+∑i=0,i≠mn-1∫x0ix0,i+1ℛi1(x)(x-s)2dx+∑i=0,i≠mn-1∫x0ix0,i+1ℛi2(x)(x-s)2dx+∫ax01f(x)-fL(x)(x-s)2dx+∫x0nbf(x)-fL(x)(x-s)2dx.
By the definition ofEm(x), we have
(44)⨎x0mx0,m+1f(x)-fL(x)(x-s)2dx=⨎x0mx0,m+1ℋm(x)(x-s)2dx+f′′(s)2⨎x0mx0,m+1(x-x0i)(x-x0,i+1)(x-s)2dx.
Putting (43) and (44) together, we have
(45)⨎abf(x)-fL(x)(x-s)2dx=∑i=1n-1⨎x0ix0,i+1f(x)-fL(x)(x-s)2dx+∫ax01f(x)-fL(x)(x-s)2dx+∫x0nbf(x)-fL(x)(x-s)2dx=-f′′(s)2hS1(τ)+ℛf(s),
where
(46)ℛf(s)=ℛ1(s)+ℛ2(s)+ℛ3(s),ℛ1(s)=⨎x0mx0,m+1ℋm(x)(x-s)2dx,ℛ2(s)=(∫ax01+∑i=1,i≠mn-1∫x0ix0,i+1+∫x0nb)ℛi1(x)(x-s)2dx+(∫ax01+∑i=1,i≠mn-1∫x0ix0,i+1+∫x0nb)ℛi2(x)(x-s)2dx,ℛ3(s)=hf′′(s)2[∑i=m∞′ϕ1′(2i+τ)+∑i=n-m∞′ϕ1′(-2i+τ)].
For the first part of ℛ2(s),
(47)|(∫ax01+∑i=1,i≠mn-1∫x0ix0,i+1+∫x0nb)ℛi1(x)(x-s)2dx|≤|(∫ax01+∑i=1,i≠m-1,m,m+1n-1∫x0ix0,i+1+∫x0nb)ℛi1(x)(x-s)2dx|+|∫x0,m-1x0,mℛm-11(x)(x-s)2dx+∫x0,m+1x0,m+2ℛm+11(x)(x-s)2dx|≤|(∫ax01+∑i=1,i≠m-1,m,m+1n-1∫x0ix0,i+1+∫x0nb)ℛi1(x)(x-s)2dx|+|∫x0,m-1x0mℛm-11(x)-ℛm-11(x0m)(x-s)2dx+∫x0,m+1x0,m+2ℛm+11(x)-ℛm+11(x0,m+1)(x-s)2dx|=|(∫ax01+∑i=1,i≠m-1,m,m+1n-1∫x0ix0,i+1+∫x0nb)ℛi1(x)(x-s)2dx|+|∫x0,m-1x0mℛm-11′(ξ0m)(x-x0m)(x-s)2dx+∫x0,m+1x0,m+2ℛm+11′(ξ0,m+1)(x-x0,m+1)(x-s)2dx|≤Ch3(∫ax01+∑i=1,i≠m-1,m,m+1n-1∫x0ix0,i+1+∫x0nb)1(x-s)2dx+Ch2(∫x0,m-1x0mdx|x-s|+∫x0,m+1x0,m+2dx|x-s|)≤C[|lnγ(τ)|+|lnh|]h2,
where ξ0m∈(x0,m-1,x0m) and ξ0,m+1∈(x0,m,x0,m+1). We have also used the identity ℛm-11(x0m) = ℛm+11(x0,m+1)=0, |x-x0m|<|x-s|, and |x-x0,m+1|<|x-s|.
For the second part ofℛ2(s),
(48)|(∫ax01+∑i=1,i≠mn-1∫x0ix0,i+1+∫x0nb)ℛi2(x)(x-s)2dx|≤Ch2(∫ax01+∑i=1,i≠mn-1∫x0ix0,i+1+∫x0nb)1|x-s|dx≤C[|lnγ(τ)|+|lnh|]h2.
By Lemmas 5 and 6, we have
(49)|Rf(s)|≤|ℛ1(s)|+|ℛ2(s)|≤C[η(s)+|lnh|+|lnγ(τ)|]h2.
Then the proof is completed.


## 3.2. The Calculation ofS1(τ)

LetQn(x) be the function of the second kind associated with the Legendre polynomial Pn(x), defined by (cf. [24])
(50)Q0(x)=12ln|x+1x-1|,Q1(x)=xQ0(x)-1.
We also define
(51)W(f,τ):=f(τ)+∑i=0∞[f(2i+τ)+f(-2i+τ)],τ∈(-1,1).
Then, by the definition of W,
(52)W(Q0)(τ)=12ln1+τ1-τ+12∑i=1∞(ln2i+1+τ2i-1+τ+ln2i-1-τ2i+1-τ)=12limi→∞ln2i+1+τ2i+1-τ=0,W(xQ0′)(τ)=τ1-τ2-∑i=1∞(2i+τ(2i+τ)2-1+-2i+τ(-2i+τ)2-1)=12limn→∞∑k=-nk=n1k+(1/2)-(τ/2)=π2tanπτ2.
It follows that
(53)S(Q1′,τ)=W(Q0+xQ0′,τ)=π2tanτπ2,
which means
(54)S(Q1,τ)=∫‍π2tanτπ2τ=-lncosπτ2+C.What remains is to determine the constantC. By using the identities (cf. [24, Chapter 1, Section 1.2]),
(55)xcotx=1+∑k=1∞(-1)kB2k(2x)2k(2k)!,ln(2cosx)=-∑j=1∞1jcos(2jx),x∈(0,π),
where B2k denote the Bernoulli numbers, we have
(56)∑k=1∞(-1)kB2k(2x)2k+1(2k+1)!=2xln(sinx)+2[(ln2-1)x+∑j=1∞12j2sin(2jx)].
Setting x=π/2 gives
(57)∑k=1∞(-1)kB2k(π)2k+1(2k+1)!=ln2-1.
Then we have
(58)W(Q1;0)=-1+2∑i=1∞Q1(2i)=-1+2∑i=1∞∑k=1∞1(2k+1)(2i)2k=-1+∑k=1∞(-1)k+1B2k(π)2k+1(2k+1)!=-ln2,
where we have used the formulae (cf. [24, Chapter 1, Section 1.2])
(59)Q1(x)=∑k=1∞1(2k+1)x2k,|x|>1,∑i=1∞1i2k=(-1)k+122k-1(2k)!B2kπ2k.
Then we have
(60)S1(τ)=-ln2cosτπ2.By Theorem2, we get the following error expansion:
(61)En(f)=f′′(s)h2log2cosπτ2+O(h2).Let(62)T(hl)=I2l-1n0(f,s),l=1,2,…,
where n0 is the starting meshes, l is the refining numbers, and hl=(b-a)/2l-1n0. Then we have the following.Corollary 7.
Under the same assumption of Theorem2 and (62), there holds
(63)T(hl)-T(hl-1)=O(hl-12).


## 4. Numerical Examples

In this section, computational results are reported to confirm our theoretical analysis.Example 1.
Consider the hypersingular integral(64)⨎01x3(x-s)2dx=3s+32+1s-1+3s2log1-ss.We examine the dynamic points=x[0n/4]+(1+τ)h/2, in Table 1 show that when the local coordinate of singular point τ=±2/3, the quadrature reach the convergence rate of O(h2) as for the nonsupersingular point, there are no convergence rate which agree with our theorematically analysis. For the case of s=b-(τ+1)h/2, Table 2 shows that there is no superconvergence phenomenon because of the influence of η(s) which coincides with our theoretical analysis.Table 1
Errors of the mod-trapezoidal rule withs=x[0n/4]+(1+τ)h/2.
n
τ
=
0
τ
=
2
/
3
τ
=
-
2
/
3
τ
=
1
/
2
32
3.5392
e
-
002
-
7.1780
e
-
004
1.1040
e
-
003
1.7473
e
-
002
64
1.6961
e
-
002
-
1.9364
e
-
004
2.7065
e
-
004
8.4165
e
-
003
128
8.3004
e
-
003
-
5.0195
e
-
005
6.6991
e
-
005
4.1332
e
-
003
256
4.1057
e
-
003
-
1.2773
e
-
005
1.6664
e
-
005
2.0484
e
-
003
512
2.0417
e
-
003
-
3.2214
e
-
006
4.1554
e
-
006
1.0198
e
-
003
1024
1.0181
e
-
003
-
8.0885
e
-
007
1.0375
e
-
006
5.0878
e
-
004
h
α
1.0239
1.9587
2.0111
1.0204Table 2
Errors of the mod-trapezoidal rule withs=x0n-(1+τ)h/2.
n
τ
=
0
τ
=
2
/
3
τ
=
-
2
/
3
τ
=
1
/
2
32
2.5227
e
-
002
2.9756
e
-
002
1.2230
e
-
001
-
8.0605
e
-
002
64
1.1997
e
-
002
1.4808
e
-
002
6.0946
e
-
002
-
4.1289
e
-
002
128
5.8224
e
-
003
7.3643
e
-
003
3.0401
e
-
002
-
2.0913
e
-
002
256
2.8619
e
-
003
3.6669
e
-
003
1.5177
e
-
002
-
1.0529
e
-
002
512
1.4173
e
-
003
1.8283
e
-
003
7.5815
e
-
003
-
5.2839
e
-
003
1024
7.0490
e
-
004
9.1253
e
-
004
3.7886
e
-
003
-
2.6471
e
-
003
h
α
1.0323
1.0054
1.0025
0.9857Example 2.
Consider the hypersingular integral(65)⨎01x4+1(x-s)2dx=4s2+2s+43+s+1s(s-1)+4s3log1-ss.The numerical results show that the convergence rate reaches(h2) 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 O(h2) which agrees with our theoretical analysis.Table 3
Errors of the mod-trapezoidal rule.
n
s
=
1
/
4
h
α
s
=
1
/
2
h
α
s
=
5
/
8
h
α
32
-
2.1121
e
-
003
-
1.8532
e
-
003
-
1.0260
e
-
003
64
-
5.3080
e
-
004
1.9925
-
4.6835
e
-
004
1.9844
-
2.6489
e
-
004
1.9536
128
-
1.3304
e
-
004
1.9963
-
1.1770
e
-
004
1.9924
-
6.7234
e
-
005
1.9781
256
-
3.3301
e
-
005
1.9982
-
2.9502
e
-
005
1.9963
-
1.6933
e
-
005
1.9894
512
-
8.3306
e
-
006
1.9991
-
7.3850
e
-
006
1.9982
-
4.2487
e
-
006
1.9948
1024
-
2.0833
e
-
006
1.9996
-
1.8474
e
-
006
1.9991
-
1.0641
e
-
006
1.9974Table 4
A posteriori error estimate of the mod-trapezoidal rule.
n
s
=
1
/
4
h
α
s
=
1
/
2
h
α
s
=
5
/
8
h
α
32
64
-
1.5814
e
-
003
-
1.3849
e
-
003
-
7.6112
e
-
004
128
-
3.9776
e
-
004
1.9912
-
3.5065
e
-
004
1.9817
-
1.9765
e
-
004
1.9452
256
-
9.9737
e
-
005
1.9957
-
8.8202
e
-
005
1.9911
-
5.0301
e
-
005
1.9743
512
-
2.4971
e
-
005
1.9979
-
2.2117
e
-
005
1.9956
-
1.2684
e
-
005
1.9875
1024
-
6.2473
e
-
006
1.9989
-
5.5375
e
-
006
1.9978
-
3.1846
e
-
006
1.9939Example 3.
Now we consider an example of less regularity. Leta=-b=-1, s=0, and
(66)f(x)=ℱi(x):=x2+(2+sign(x))|x|2-i+0.5,i=-1,0,1.
Obviously, ℱi(x)∈C3-i+0.5[-1,1](i=-1,0,1). The exact value of the integral is
(67)ℐ2(ℱi(x),0)=14-4i3-2i.The numerical results are presented in Tables5 and 6. When the density function f(x) is smooth enough i=-1, the error bound is O(h2), and if the density function has less regularity i=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.
n
i
=
-
1
h
α
i
=
0
h
α
i
=
1
h
α
32
-
1.9220
e
-
002
-
5.6378
e
-
002
1.1169
64
-
4.9807
e
-
003
1.9481
-
2.1100
e
-
002
1.4179
7.8885
e
-
001
5.0167
e
-
001
128
-
1.2750
e
-
003
1.9659
-
7.7489
e
-
003
1.4451
5.5758
e
-
001
5.0058
e
-
001
256
-
3.2386
e
-
004
1.9770
-
2.8116
e
-
003
1.4626
3.9421
e
-
001
5.0021
e
-
001
512
-
8.1850
e
-
005
1.9843
-
1.0120
e
-
003
1.4742
2.7874
e
-
001
5.0007
e
-
001
1024
-
2.0617
e
-
005
1.9892
-
3.6228
e
-
004
1.4820
1.9709
e
-
001
5.0003
e
-
001Table 6
A posteriori error estimate of the mod-trapezoidal rule.
n
i
=
-
1
h
α
i
=
0
h
α
i
=
1
h
α
32
64
-
1.4239
e
-
002
-
3.5278
e
-
002
3.2804
e
-
001
128
-
3.7057
e
-
003
1.9420
-
1.3351
e
-
002
1.4019
2.3128
e
-
001
5.0427
e
-
001
256
-
9.5112
e
-
004
1.9621
-
4.9373
e
-
003
1.4351
1.6337
e
-
001
5.0150
e
-
001
512
-
2.4201
e
-
004
1.9746
-
1.7996
e
-
003
1.4560
1.1548
e
-
001
5.0053
e
-
001
1024
-
6.1233
e
-
005
1.9827
-
6.4972
e
-
004
1.4698
8.1643
e
-
002
5.0019
e
-
001



---

*Source: 736834-2013-12-10.xml*