# A numerical method for Hadamard finite-part integrals with an integral power singularity at an endpoint

In this paper, we propose a numerical method for computing Hadamard finite-part integrals with an integral-power singularity at an endpoint, the part of the divergent integral which is finite as a limiting procedure. In the proposed method, we express the desired finite-part integral using a complex loop integral, and obtain the finite-part integral by evaluating the complex integral by the trapezoidal rule. Theoretical error estimate and some numerical examples show the effectiveness of the proposed method.

## Authors

• 8 publications
09/25/2019

### Numerical method for computing Hadamard finite-part integrals with a non-integral power singularity at an endpoint

In this paper, we propose a numerical method of computing a Hadamard fin...
10/08/2019

### Numerical method of computing Hadamard finite-part integrals with a non-integral power singularity at the endpoint over a half infinite interval

In this paper, we propose a numerical method of computing an Hadamard fi...
09/11/2019

### A numerical method of computing oscillatory integral related to hyperfunction theory

In this paper, we propose a numerical method of computing an integral wh...
10/21/2021

### A boundary integral method for 3D nonuniform dielectric waveguide problems via the windowed Green function

This paper proposes an efficient boundary-integral based "windowed Green...
10/25/2021

### Cubature Method for Stochastic Volterra Integral Equations

In this paper, we introduce the cubature formulas for Stochastic Volterr...
03/23/2022

### Efficient extraction of resonant states in systems with defects

We introduce a new numerical method to compute resonances induced by loc...
11/24/2021

### A topology optimisation of acoustic devices based on the frequency response estimation with the Padé approximation

We propose a topology optimisation of acoustic devices that work in a ce...
##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

The integral

 ∫10x−1f(x)dx,

where is an analytic function on the closed interval such that , is divergent. However, if is analytic on the closed interval , for such that , we have the following using integral by part.

 ∫1ϵx−1f(x)dx= ∫1ϵ(logx)′f(x)dx = [logxf(x)]1ϵ−∫1ϵlogxf′(x)dx = −f(ϵ)logϵ−∫1ϵlogxf′(x)dx = −f(0)logϵ+(terms finite as ϵ↓0).

Therefore, we can define the so-called Hadamard finite-part (f.p.) integral by

 f.p.∫10x−1f(x)dx=limϵ↓0{∫1ϵx−1f(x)+f(0)logϵ}

[5]. Similarly, we can define the f.p. integrals

 f.p.∫10x−nf(x)dx(n=1,2,…). (1)

We propose a numerical method for computing the f.p. integrals (1). In the proposed method, we express the f.p. integral by a complex loop integral, and we obtain the f.p. integral by evaluating the complex integral by the trapezoidal formula with equal mesh.

Previous works related to this paper are as follows. Ogata and Hirayama proposed a numerical integration method based on hyperfunction theory, a theory of generalized functions based on complex function theory, where they obtain ordinary integrals by expressing them as complex integrals and evaluating them by the conventional numerical integral formulas [6]. For Cauchy principal value integrals and Hadamard f.p. integrals with a singularity inside the integral interval

 f.p.∫10f(x)(x−λ)ndx(n=1,2,…), (2)

many approximation methods have been proposed. Elliot and Paget proposed a Gauss type numerical integration formula for Cauchy principal value integrals (2) with [4], and Paget proposed a Gauss type formula for Hadamard finite-part integrals (2) with [8]. Bialecki proposed approximation formulas for (2) based on the Sinc method [2, 1], that is, methods using the trapezoidal formula together with variable transforms as in the DE formula [9]. The author et al. improved them and proposed a DE-type numerical integration formula for Cauchy principal-value integrals and Hadamard finite-part integrals with an integral power singularity inside the integral interval [7].

The remainder of this paper is structured as follows. In Section 2, we define the f.p. integrals (1) and show the expression of them by complex loop integral. Then, we give an approximation formula for the desired f.p. integral. In Section 3, we show some numerical examples which show the effectiveness of the proposed method. In Section 4, we give a summary of this paper.

## 2 Hadamard finite-part integrals and its approximation

We define the Hadamard finite-part integrals (1) by

 f.p.∫10x−nf(x)dx=limϵ↓0{∫1ϵx−nf(x)−n−2∑k=0ϵk+1−nk!(n−1−k)f(k)(0)+logϵ(n−1)!f(n−1)(0)}(n=1,2,…), (3)

where the integrand is analytic on the closed interval , and the second term on the right-hand side is zero if . We can show that it is well-defined using integral by part as follows.

 ∫1ϵx−nf(x)dx = = −1n−1{[x−(n−1)f(x)]1ϵ−∫1ϵx−(n−1)f′(x)dx} = −f(1)n−1+ϵ1−nn−1f(ϵ)+1n−1∫1ϵx−(n−1)f′(x)dx = ϵ1−nn−1n−2∑k=0ϵkk!f(k)(0)−1(n−1)(n−2)∫1ϵ(x−(n−2))′f′(x)dx +(terms finite as ϵ↓0, which is denoted % by ⋯'' below) = ϵ1−nn−1n−2∑k=0ϵkk!f(k)(0)+ϵ2−n(n−1)(n−2)n−3∑k=0ϵkk!f(k+1)(0) −1(n−1)(n−2)(n−3)∫1ϵ(x−(n−3))′f′′(x)dx+⋯ = ⋯ = ϵ1−nn−1n−2∑k=0ϵkk!f(k)(0)+ϵ2−n(n−1)(n−2)n−3∑k=0ϵkk!f(k+1)(0) +ϵ3−n(n−1)(n−2)(n−3)n−4∑k=0ϵkk!f(k+2)(0)+⋯+ϵ−1(n−1)!f(n−2)(0)−logϵ(n−1)!f(n−1)(0) −1(n−1)!∫1ϵlogxf(n)(x)dx+⋯ = ϵ1−nn−1n−2∑k=0ϵkk!f(k)(0)+ϵ2−n(n−1)(n−2)n−2∑k=1ϵk−1(k−1)!f(k)(0) +ϵ3−n(n−1)(n−2)(n−3)n−2∑k=2ϵk−2(k−2)!f(k)(0)+⋯+ϵ−1(n−1)!f(n−2)(0)−logϵ(n−1)!f(n−1)(0) −1(n−1)!∫1ϵlogxf(n)(x)dx+⋯ = ϵ1−nn−1f(0)+ϵ2−nf′(0){1n−1+1(n−1)(n−2)} +ϵ1−nn−1n−2∑k=2ϵkk!f(k)(0)+ϵ2−n(n−1)(n−2)n−2∑k=2ϵk−1(k−1)!f(k)(0) +ϵ3−n(n−1)(n−2)(n−3)n−2∑k=2ϵk−2(k−2)!f(k)(0)+⋯+ϵ−1(n−1)!f(n−2)(0)−logϵ(n−1)!f(n−1)(0) −1(n−1)!∫1ϵlogxf(n)(x)dx+⋯ = ϵ1−nn−1f(0)+ϵ2−nn−2f′(0)+ϵ3−nn−1f′′(0){12!+1n−2+1(n−2)(n−3)} +ϵ1−nn−1n−2∑k=3ϵkk!f(k)(0)+ϵ2−n(n−1)(n−2)n−2∑k=3ϵk−1(k−1)!f(k)(0) +ϵ3−n(n−1)(n−2)(n−3)n−2∑k=3ϵk−2(k−2)!f(k)(0)+⋯+ϵ−1(n−1)!f(n−2)(0)−logϵ(n−1)!f(n−1)(0) +⋯ = ϵ1−nn−1f(0)+ϵ2−nn−2f′(0)+ϵ3−n2!(n−3)f′′(0) +ϵ4−nn−1f′′′(0){13!+12!(n−2)+1(n−2)(n−3)+1(n−2)(n−3)(n−4)1(n−2)(n−4)12!(n−4)} +ϵ1−nn−1n−2∑k=4ϵkk!f(k)(0)+ϵ2−n(n−1)(n−2)n−2∑k=4ϵk−1(k−1)!f(k)(0) +ϵ3−n(n−1)(n−2)(n−3)n−2∑k=4ϵk−2(k−2)!f(k)(0) +ϵ4−n(n−1)(n−2)(n−3)(n−4)n−2∑k=4ϵk−3(k−3)!f(k)(0)+⋯+ϵ−1(n−1)!f(n−2)(0) −logϵ(n−1)!f(n−1)(0)+⋯ = ϵ1−nn−1f(0)+ϵ2−nn−2f′(0)+ϵ3−n2!(n−3)f′′(0)+ϵ4−n3!(n−4)f′′′(0) +ϵ1−nn−1n−2∑k=4ϵkk!f(k)(0)+ϵ2−n(n−1)(n−2)n−2∑k=4ϵk−1(k−1)!f(k)(0) +ϵ3−n(n−1)(n−2)(n−3)n−2∑k=4ϵk−2(k−2)!f(k)(0) +ϵ4−n(n−1)(n−2)(n−3)(n−4)n−2∑k=4ϵk−3(k−3)!f(k)(0)+⋯+ϵ−1(n−1)!f(n−2)(0) −logϵ(n−1)!f(n−1)(0)+⋯ = ϵ1−nn−1f(0)+ϵ2−nn−2f′(0)+ϵ3−n2!(n−3)f′′(0)+ϵ4−n3!(n−4)f′′′(0)+⋯+ϵ−1(n−2)!f(n−2)(0) −logϵ(n−1)!f(n−1)(0)+⋯.

The f.p. integral is expressed using a complex loop integral as in the following theorem.

###### Theorem 1

We suppose that is analytic in a complex domain containing the closed interval in its interior. Then, the f.p. integral (3) is expressed as

 f.p.∫10x−nf(x)dx=12πi∮Cz−nf(z)log(zz−1)dz−n−2∑k=0f(k)(0)k!(n−1−k)(n=1,2,…), (4)

where is a closed complex integral path in encircling the interval in the positive sense, and the second term on the right-hand side of (4) is zero if .

### Proof of Theorem 1

Using Cauchy’s integral theorem, we have

 12πi∫Cz−nf(z)log(zz−1)dz=12πi(∫C(0)ϵ+∫C(1)ϵ+∫Γ(+)ϵ+∫Γ(−)ϵ)z−nf(z)log(zz−1)dz, (5)

where the integral paths , , and are respectively

 C(0)ϵ= {ϵeiθ|0≦θ≦π}, C(1)ϵ= {1+ϵeiθ|0≦θ≦π}, Γ(+)ϵ= {x∈R|1−ϵ≧x≧ϵ}, Γ(−)ϵ= {x∈R|ϵ≦x≦1−ϵ}

with small (see Figure 1).

From Table 1, we have

 12πi(∫Γ(+)ϵ+∫Γ(−)ϵ)z−nf(z)log(zz−1)dz = −12πi∫1−ϵϵx−nf(x){log(x1−x−iπ)}dx+12πi∫1−ϵϵx−nf(x){log(x1−x+iπ)}dx = ∫1−ϵϵx−nf(x)dx=∫1ϵx−nf(x)dx+O(ϵ).

The integral on is written as

 12πi∫C(0)ϵz−nf(z)log(zz−1)dz = 12πi∫2π0ϵ−ne−inθf(ϵeiθ)log(−ϵeiθ1−ϵeiθ)iϵeiθdθ = ϵ1−n2π∫2π0f(ϵeiθ)log(ϵei(θ−π))e−i(n−1)θdθ−ϵ1−n2π∫2π0f(ϵeiθ)log(1−ϵeiθ)e−i(n−1)θdθ.

The first integral on the right-hand side is written as

 ϵ1−n2π∫2π0f(ϵeiθ)log(ϵei(θ−π))e−i(n−1)θdθ = ϵ1−n2π∫2π0{∞∑k=0ϵkk!f(k)(0)eikθ}{logϵ+i(θ−π)}e−i(n−1)θdθ = ϵ1−n2πlogϵ∞∑k=0∫2π0ei(k−n+1)θdθ+i(n−1)ϵ1−n2π∞∑k=0ϵkk!f(k)(0)∫2π0(θ−π)ei(k−n+1)θdθ = logϵ(n−1)!f(n−1)(0)−n−2∑k=0ϵk−n+1k!(n−1−k)f(k)(0)+O(ϵ),

where we exchanged the order of the integral and the infinite summation on the second equality since the infinite sum is uniformly convergent on . Similarly, the second integral is written as

 ϵ1−n2π∫2π0f(ϵeiθ)log(1−ϵeiθ)e−i(n−1)θdθ=−n−2∑k=0f(k)(0)k!(n−1−k).

Then, we have

 12πi∫C(0)ϵz−nf(z)log(zz−1)dz=−n−2∑k=0ϵk−n+1k!(n−1−k)f(k)(0)+logϵ(n−1)!f(n−1)(0)+n−2∑k=0f(k)(0)k!(n−1−k)+O(ϵ).

As to the integral on , we have

 12πi∫C(1)ϵz−nf(z)log(zz−1)dz=O(ϵlogϵ)

since is analytic near . Summarizing the above calculations, we have

 12πi∮Cz−nf(z)log(zz−1)dz=∫1ϵx−nf(x)dx−n−2∑k=0ϵk−n+1k!(n−1+k)f(k)(0)+logϵ(n−1)!f(n−1)(0)+n−2∑k=0f(k)(0)k!(n−1+k)+O(ϵlogϵ).

Taking the limit , we have (4).

The complex integral in (4) is the integral of an analytic function over an interval of the length of one period, and it is accurately approximated by the trapezoidal formula with equal mesh. Using a parameterization of the closed integral path

 C:z=φ(u),0≦u≦up,

where is a periodic function of period , we obtain the following approximation formula for the f.p. integral.

 f.p.∫10x−nf(x)dx≃I(n)N[f]≡h2πiN−1∑k=0φ(kh)−nf(φ(kh))log(φ(kh)φ(kh)−1)φ′(kh)−n−2∑k=0f(k)(0)k!(n−1−k)(h=upN), (6)

where the second term on the right-hand side is zero if .

We remark here that, if the integrand is real valued on and the integral path is symmetric with respect to the real axis, we can reduce the number of sampling points by half. In fact, in this case, we have due to the reflection principle, and , . Then, we have

 f.p.∫10x−nf(x)dx≃I′(n)N[f] ≡ h2πIm{φ(0)−nf(φ(0))log(φ(0)1−φ(0))φ′(0) +φ(up2)−nf(φ(up2))log(φ(up/2)1−φ(up/2))φ′(up2)} +hπIm{N−1∑k=1φ(kh)−nf(φ(kh))log(φ(kh)φ(kh)−1)φ′(kh)} −n−2∑k=0f(k)(0)k!(n−1−k)(h=up2N), (7)

Applying the theorem in §4.6.5 in [3] to the approximation of the complex integral by the trapezoidal formula in (6), we have the following theorem on the error estimate of the approximation formula (6).

###### Theorem 2

We suppose that

• the strip domain

 Dd={w∈C||Imw|0)

is contained in ,

• the parameterization function of is analytic in , and

• the integrand is analytic in

 φ(Dd)={φ(w)|w∈Dd}.

Then, we have the following inequality for arbitrary .

 ∣∣∣f.p.∫10x−nf(x)dx−I(n)N[f]∣∣∣≦dπN(f,n,d′)exp(−2πd′N/up)1−exp(−2πd′N/up), (8) where N(f,n,d′)=max|Imw|=d′∣∣ ∣∣φ(w)−nf(φ(w))log(φ(w)1−φ(w)).∣∣ ∣∣ (9)

This theorem says that the approximation (6) converges exponentially as increases if the integrand function is analytic on and the integral path is an analytic curve.

## 3 Numerical examples

We computed the integrals

 (1) f.p.∫10x−nexdx=∞∑k=0(k≠n−1)1k!(k−n+1), (2) f.p.∫10x−n1+xdx=(−1)n{log2+n−1∑l=1(−1)ll},

where the second term on the right-hand side of the integral (2) is zero if , by the proposed method. All the computations were performed using programs coded in C++ with double precision working. The complex integral path was taken as the ellipse

 C:z=12+14(ρ+1ρ)cosu+14(ρ−1ρ)sinu,0≦u≦2π(ρ>1),

where the parameter is taken as for the integral (1) and for the integral (2). Figure 2 shows the relative errors of the approximation formula (7) applied to the integrals (1) and (2) as functions of the number of sampling points . From these figures, the errors decay exponentially as increases. Table 2 shows the decay rates of the errors of the proposed method.

## 4 Summary

We proposed a numerical method for Hadamard finite-part integrals with an integral order power singularity at an endpoint over a finite interval. In the proposed method, we express the desired f.p. integral using a complex loop integral, and we obtain the f.p. integral by evaluating the complex integral by the trapezoidal formula. Theoretical error estimate and numerical examples show the exponential convergence of the proposed method.

We can also give approximation methods for f.p. integrals with a non-integral power singularity and f.p. integrals over a half-infinite interval in a way similar to this paper. They will be reported in other papers.

## References

• [1] B. Bialecki (1990) A sinc quadrature rule for hadamard finite-part integrals. Numer. Math. 57, pp. 263–269. External Links: Document Cited by: §1.
• [2] B. Bialecki (1990) A sinc-hunter quadrature rule for cauchy principal value integrals. Math. Comput. 55, pp. 665–681. External Links: Document Cited by: §1.
• [3] P. J. Davis and P. Rabinowitz (1984) Methods of numerical integration, second ed.. Academic Press, San Diego. Cited by: §2.
• [4] D. Elliot and D. F. Paget (1979) Gauss type quadrature rules for cauchy principal value integrals. Math. Comput. 33, pp. 301–309. External Links: Document Cited by: §1.
• [5] R. Estrada and R. P. Kanwal (1989) Regularization, pseudofunction, and hadamard finite part. J. Math. Anal. Appl. 141, pp. 195–207. External Links: Document Cited by: §1.
• [6] H. Ogata and H. Hirayama (2018) Numerical integration based on hyperfunction theory. J. Comput. Appl. Math. 327, pp. 243–259. Cited by: §1.
• [7] H. Ogata, M. Sugihara, and M. Mori (2000) DE-type quadrature formulae for cauchy principal-value integrals and for hadamard finite-part itnegrals. In Proceedings of the Second ISAAC Congress, Vol. 1, pp. 357–366. Cited by: §1.
• [8] D. F. Paget (1981) The numerical evaluation of hadamard finite-part integrals. Numer. Math. 36, pp. 447–453. External Links: Document Cited by: §1.
• [9] H. Takahasi and M. Mori (1978) Double exponential formulas for numerical integration. Publ. Res. Inst. Math. Sci., Kyoto Univ. 339, pp. 721–741. Cited by: §1.