# 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 finite-part integral with a non-integral power singularity at an endpoint, that is, a finite part of a divergent integral 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 formula. Theoretical error estimate and some numerical examples show the effectiveness of the proposed method.

09/19/2019

### 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 fini...
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/02/2019

### A numerical method of computing Hadamard finite-part integrals with an integral power singularity at the endpoint on a half infinite interval

In this paper, we propose a numerical method of computing Hadamard finit...
01/12/2018

### Bayesian Quadrature for Multiple Related Integrals

Bayesian probabilistic numerical methods are a set of tools providing po...
12/02/2020

### A Fast Numerical solution of the quark's Dyson-Schwinger equation with Ball-Chiu vertex

In this paper, we present two feasible and efficient methods to numerica...
12/03/2020

### Computing the matrix fractional power based on the double exponential formula

Two quadrature-based algorithms for computing the matrix fractional powe...
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...

## 1 Introduction

The integral

 ∫10xα−2f(x)dx(0<α<1),

where is an analytic function on the closed interval , is divergent. However, we can assign a finite value to this divergent integral as follows. For , using integration by part, we have

 ∫1ϵxα−2f(x)dx= 1α−1∫1ϵ(xα−1)′f(x)dx = 1α−1[xα−1f(x)]1ϵ−∫1ϵxα−1f′(x)dx = ϵα−1f(ϵ)−f(1)1−α−∫1ϵxα−1f′(x)dx = ϵα−1f(0)1−α+(terms finite as ϵ↓0),

and the limit

 limϵ↓0{∫1ϵxα−2f(x)dx−ϵα−1f(0)1−α}

is finite. We call this limit an Hadamard finite-part (f.p.) integral and denote it by

 f.p.∫10xα−2f(x)dx.

Similarly, we can define a f.p. integral

 f.p.∫10xα−1−nf(x)dx (1)

for , and an analytic function on [6].

In this paper, we propose a numerical method of computing f.p. integrals (1). In the proposed method, we express the f.p. integral using a complex loop integral, and we obtain the desired f.p. integral by evaluating the complex integral by the trapezoidal formula with equal mesh. Theoretical error estimate and numerical examples will show that the approximation by the proposed method converges exponentially as the number of sampling points increases.

Previous works related to this paper are as follows. The author and Hirayama proposed a numerical integration method for ordinary integrals related to hyperfunction theory [8], where a desired integral is expressed using a complex loop integral, and it is obtained by evaluating the complex integral by the trapezoidal formula with equal mesh. The author proposed a numerical method of computing f.p. integrals with an integral order singularity [10]. For Cauchy principal value integrals and Hadamard finite-part integrals with a singularity in the interior of the integral interval

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

many numerical methods were proposed. Elliot and Paget proposed a Gauss-type numerical integration formulas for f.p. integrals (2) [5, 11]. Bialecki proposed Sinc numerical integration formulas for f.p. integrals [3, 2], where the trapezoidal formula with the variable transform technique are used as in the DE formula for ordinary integrals [12]. The author et al. improved them and proposed a DE-type numerical integration formulas for f.p. integrals (2) [9].

The remainder of this paper is structured as follows. In Section 2, we define the f.p. integrals and propose a numerical method of computing them. Then, we give a theorem on error estimate of the proposed method. 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.

The Hadamard finite-part integral is defined by

 (3)

We can see that it is well-defined using integration by part. In fact, repeating integration by part, we have

 ∫1ϵxα−n−1f(x)dx = ϵα−nn−αf(ϵ)+ϵα−n+1(n−α)(n−1−α)f′(ϵ)+ϵα−n+2(n−α)(n−1−α)(n−2−α)f′′(ϵ) +⋯+ϵα−1(n−α)(n−1−α)(n−2−α)⋯(1−α)f(n−1)(0) +1(n−α)(n−1−α)(n−2−α)⋯(1−α)∫1ϵxα−1f(n)(x)dx +(terms finite as n↓0, which will be denoted by% ⋯'' below) = ϵα−nn−αn−1∑k=0ϵkk!f(k)(0)+ϵα−n+1(n−α)(n−1−α)n−2∑k=0ϵkk!f(k)(0) +ϵα−n+2(n−α)(n−1−α)(n−2−α)n−3∑k=0ϵkk!f(k)(0) +⋯+ϵα−1(n−α)(n−1−α)(n−2−α)⋯(1−α)f(n−1)(0)+⋯ = ϵα−nn−αf(0)+ϵα−n+1n−α(1+1n−1−α)f′(0) +ϵα−n+2n−α{12!+1n−1−α+1(n−1−α)(n−2−α)}f′′(0) +ϵα−n+3n−α{13!+12!(n−1−α)+1(n−1−α)(n−2−α) +1(n−1−α)(n−2−α)(n−3−α)}f′′′(0) +⋯ +ϵα−1n−α{1(n−1)!(n−1−α)+1(n−2)!(n−1α)(n−2−α) +1(n−3)!(n−1−α)(n−2−α)(n−3−α) +⋯+1(n−1−α)(n−2−α)⋯(1−α)}f(n−1)(0) +⋯ = ϵα−nn−αf(0)+ϵα−n+1n−1−αf′(0)+ϵα−n+22!(n−2−α)f′′(0)+⋯+ϵα−1(n−1)!(1−α)f(n−1)(0) +⋯.

If the integrand is an analytic function on the closed interval , the f.p. integral (3) is expressed using a complex loop integral as in the following theorem.

###### Theorem 1

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

 f.p.∫10xα−1−nf(x)dx=12πi∮Cz−nf(z)Ψα(z)dz+n−1∑k=0f(k)(0)k!(α−n+k), (4) where Ψα(z)=α−1z−1F(α,1;α+1;z−1), (5)

and is a closed complex integral path contained in and encircling the interval in the positive sense.

### Proof of Theorem 1

From Cauchy’s integral theorem, the complex integral of the first term on the right-hand side of (4) is modified into

 12πi∮Cz−nf(z)Ψα(z)dz=(∫C(0)ϵ+∫Γ(−)ϵ+∫C(1)ϵ+∫Γ(+)ϵ)z−nf(z)Ψα(z)dz,

where , , and are complex paths respectively defined by

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

with (see Figure 1).

From the formula 15.3.7 in [1], we have

 Ψα(z)=−πsinπα(−z)α−1−1α−1F(1−α,1;2−α;z)(|arg(−z)|<π).

Then, as to the integrals on , we have

 12πi(∫Γ(+)ϵ+∫Γ(1)ϵ)z−nf(z)Ψα(z)dz = −12isinπα{−∫1−ϵϵx−nf(x)(−(x+i0))α−1dx+∫1−ϵϵx−nf(x)(−(x−i0))α−1dx} = −12isinπα(−e−iπ(α−1)+eiπ(α−1))∫1−ϵϵxα−n−1f(x)dx = ∫(1+ϵ)−1ϵxα−n−1f(x)dx,

where we remark that is a single-valued analytic function on the interval . As to the integral on , we have

 12πi∫C(0)ϵz−nf(z)Ψα(z)dz = −12isinπα∫C(0)ϵz−nf(z)(−z)αdz−12πi(α−1)∫C(0)ϵz−nf(z)F(1−α,1;2−α;z)dz.

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

 (the first term)= (−1)n+12isinπα∫C(0)ϵf(z)(−z)α−n−1dz = (−1)n+12isinπα∫2π0f(ϵeiθ){ϵei(θ−π)}α−n−1iϵeiθdθ = ϵα−ne−iπα2sinπα∫2π0f(ϵeiθ)ei(α−n)θdθ = = −n−1∑k=0ϵα−n+kk!(n−k−α)f(k)(0)+O(ϵα),

and the second term is written as

 (the second term) = −12πi(α−1)∫C(0)ϵz−n{∞∑k=0f(k)(0)k!zk}{∞∑l=01−α1−α+lzl}dz = −∑k+l=n−1f(k)(0)k!(α−l−1)=−n−1∑k=0f(k)(0)k!(α−n+k).

As to the integral on , from the formula 15.3.10 in [1], we have

 Ψα(z)=z−1∞∑k=0(α)kk!{ψ(k+1)−ψ(k+α)−log(1−z−1)}(1−z−1)k, where ψ(z) is the Digamma function: ψ(z)=Γ′(z)/Γ(z) and (α)0=1,(α)k=α(α+1)(α+2)⋯(α+k−1)(k=1,2,…),

and, then, the integral on is of . Summarizing the above calculations, we have

 12πi∮Cz−nf(z)Ψα(z)dz= ∫1ϵxα−n−1f(x)dx−n−1∑k=0ϵα−n+kk!(n−k−α)f(k)(0) −n−1∑k=0f(k)(0)k!(α−n+k)+O(ϵlogϵ),

and, taking the limit , we have (4).

We can obtain the desired f.p. integral (3) by evaluating the complex integral in (4) on the closed integral path , which is parameterized by , , by the trapezoidal formula with equal mesh as follows.

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

The hypergeometric function in the definition of in (5) is easily evaluated using the continued fraction expansion (see §12.5 in [7]). If is an analytic curve, the complex loop integral is an integral of an analytic periodic function on an interval of one period, to which the trapezoidal formula with equal mesh is very effective, and, then, the approximation formula (6) is very accurate. In fact, applying the theorem in §4.6.5 in [4] to the approximation of the complex integral in (6), we have the following theorem on error estimate of the proposed approximation.

###### Theorem 2

We suppose that

• the parameterization function of is analytic in the strip domain

 Dd={z∈C||Imz|0),
• the domain

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

is contained in , and

• the function is analytic in .

Then, we have for arbitrary

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

This theorem says that the proposed approximation (6) converges exponentially as the number of sampling points increases if is an analytic periodic function and is an analytic curve.

We remark here that, if is real valued on the real axis, we can reduce the number of sampling points by half. In fact, in this case, we have

 f(¯¯¯z)=¯¯¯¯¯¯¯¯¯¯f(z)

from the reflection principle, and, taking the integral path to be symmetric with respect to the real axis, that is,

 φ(−u)=¯¯¯¯¯¯¯¯¯¯¯φ(u),φ′(−u)=−¯¯¯¯¯¯¯¯¯¯¯¯φ′(u),

we have

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

## 3 Numerical examples

In this section, we show some numerical examples which show the effectiveness of the proposed method. We computed the integrals

 (i) f.p.∫10xα−n−1exdx=F(α−n;α+1−n;1)α−n, (10) (ii)

with by the approximation formula (9). 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≦up(ρ>0)

with for the integral (i) and for the integral (ii). Figure 2 show the relative errors of the proposed method applied to the integrals (i) and (ii) as functions of the number of sampling points . From these figures, the errors of the proposed method decays exponentially as increases, and the decay rate of the error does not depend much on . Table 1 shows the decay rates of the errors of the proposed method applied to the f.p. integrals (i) and (ii).

## 4 Summary

In this paper, we proposed a numerical method of computing Hadamard finite part integrals with a non-integral power singularity on an endpoint. In the proposed method, we express the desired f.p. integral using a complex loop integral, and obtain the f.p. integral by evaluating the complex integral by the trapezoidal formula with equal mesh. Theoretical error estimate and some numerical examples showed the exponential convergence of the proposed method.

We can obtain similarly f.p. integrals on an infinite interval. This will be reported in a future paper.

## References

• [1] M. Abramowitz and I. A. S. (eds.) (1965) Handbook of mathematical functions with formulas, graphs and mathematical tables. Dover, New York. Cited by: §2.
• [2] B. Bialecki (1990) A sinc quadrature rule for hadamard finite-part integrals. Numer. Math. 57, pp. 263–269. External Links: Document Cited by: §1.
• [3] 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.
• [4] P. J. Davis and P. Rabinowitz (1984) Methods of numerical integration, second ed.. Academic Press, San Diego. Cited by: §2.
• [5] 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.
• [6] 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.
• [7] P. Henrici (1977) Applied and computational complex analysis. Vol. 2, John Wiley & Sons, New York. Cited by: §2.
• [8] H. Ogata and H. Hirayama (2018) Numerical integration based on hyperfunction theory. J. Comput. Appl. Math. 327, pp. 243–259. Cited by: §1.
• [9] 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.
• [10] H. Ogata (2019) A numerical method for hadamard finite-part integrals with an integral power singularity at an endpoint. Note: arXiv:1909.08872v1 [math.NA] Cited by: §1.
• [11] D. F. Paget (1981) The numerical evaluation of hadamard finite-part integrals. Numer. Math. 36, pp. 447–453. External Links: Document Cited by: §1.
• [12] 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.