# 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 finite-part integral, a finite value assigned to a divergent integral, with a non-integral power singularity at the endpoint on a half infinite interval. In the proposed method, we express a desired finite part integral using a complex integral, and we obtain the finite part integral by evaluating the complex integral by the DE formula. Theoretical error estimate and some numerical examples show the exponential convergence of the proposed method.

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...
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...
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...
12/18/2018

### Evaluating the squared-exponential covariance function in Gaussian processes with integral observations

This paper deals with the evaluation of double line integrals of the squ...
11/03/2019

### Dimension-free path-integral molecular dynamics without preconditioning

Convergence with respect to imaginary-time discretization (i.e., the num...
12/19/2020

### The Complex-Scaled Half-Space Matching Method

The Half-Space Matching (HSM) method has recently been developed as a ne...

## 1 Introduction

The integral

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

where is an analytic function on such that and as , is divergent. However, we can assign a finite value to this divergent integral. In fact, for , we have by integration by part

 ∫∞ϵxα−2dx= −11−α∫∞ϵ(xα−1)′f(x)dx = −11−α{[xα−1f(x)]∞ϵ−∫∞ϵxα−1f′(x)dx} = ϵα−1f(ϵ)1−α+11−α∫∞ϵxα−1f′(x)dx = ϵα−1f(0)1−α+O(1)% asϵ↓0,

and the limit

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

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

 f.p.∫∞0xα−2f(x)dx.

More generally, we can define the f.p. integral

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

for , and an analytic function on such that and as [2].

In this paper, we propose a numerical method of computing a f.p. integral (1). In the proposed method, we express the desired f.p. integral using a complex integral, and obtain the f.p. integral by evaluating the complex integral by the DE formula [11]. Theoretical error estimate and some numerical examples show that the proposed approximation formula converges exponentially as the number of sampling points increases.

Previous studies related to this paper are as follows. The author and Hirayama proposed a numerical method of computing ordinary integrals based on hyperfunction theory, a theory of generalized functions based on complex function theory [4]. In their method, we obtain the desired integral by evaluating a complex integral using a conventional numerical integration formula as in the method for computing f.p. integrals proposed in this paper. The author proposed numerical methods of computing a f.p. integral with a singularity at an endpoint on a finite interval [6, 7] and a f.p. integral with an integral power singularity at the endpoint on a half infinite interval [8]. Also in these methods, we obtain a desired integral using a complex integral, and obtain the integral by evaluating the complex integral by a conventional numerical integration formula. For the computation of a Cauchy principal value integral or a f.p. integral on a finite interval with a singularity in the interior of the integral interval

 f.p.∫baf(x)(x−λ)ndx(−∞

many numerical methods were proposed. Elliot and Paget proposed Gauss-type numerical integration formulas for (2) [1, 9]. Bialecki proposed Sinc numerical integration formula of computing (2), where the trapezoidal formula together with a variable transform technique are used as in the DE formula [11]. The author et al. improved these methods and proposed a DE-type numerical integration formula of computing (2) [5].

The remainder of this paper is structured as follows. In Section 2, we define the f.p. integral (1) and propose a numerical method of computing it. In addition, we show a theoretical error estimate which shows the exponential convergence 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.

## 2 Hadamard finite-part integral and a numerical method

Let , , and be an analytic function on such that and as . The Hadamard finite-part integral (1) is defined by

 I(n,α)[f]= f.p.∫∞0xα−1−nf(x)dx = limϵ↓0{∫∞ϵxα−1−nf(x)dx−n−1∑k=0ϵα−n+kk!(n−α−k)f(k)(0)}. (3)

We can show that (3) is well-defined as follows. In fact, by integration by part, we have for

 ∫∞ϵxα−1−nf(x)dx = −1n−α∫∞ϵ(xα−n)′f(x)dx = −1n−α{[xα−nf(x)]∞ϵ−∫∞ϵxα−nf′(x)dx} = ϵα−nf(ϵ)n−α+1n−α∫∞ϵxα−nf′(x)dx = ϵα−nf(ϵ)n−α+ϵα−n+1f′(ϵ)(n−α)(n−α−1)+1(n−α)(n−α−1)∫∞ϵxα−n+1f′′(x)dx = ⋯ = n−1∑k=0ϵα−n+kf(k)(ϵ)(n−α)(n−α−1)⋯(n−α−k) +1(n−α)(n−α−1)⋯(1−α)∫∞ϵxα−1f(n)(x)dx = n−1∑k=0ϵα−n+k(n−α)(n−α−1)⋯(n−α−k)∞∑l=0ϵll!f(l+k)(0) +1(n−α)(n−α−1)⋯(1−α)∫∞ϵxα−1f(n)(x)dx = n−1∑k=0n−k−1∑l=0ϵα−n+k+lf(l+k)(0)l!(n−α)(n−α−1)⋯(n−α−k)+O(1) = n−1∑m=0{∑k+l=m1l!(n−α)(n−α−1)⋯(n−α−k)}ϵα−n+mf(m)(0)+O(1) = n−1∑m=0{m∑k=01(m−k)!(n−α)(n−α−1)⋯(n−α−k)}(∗)ϵα−n+mf(m)(0) +O(1)as ϵ↓0,

and

 (∗)= 1m!(n−α)+1(m−1)!(n−α)(n−α−1)+⋯ +11!(n−α)⋯(n−α−m+2)(n−α−m+1) +1(n−α)⋯(n−α−m+2)(n−α−m+1)(n−α+m) = 1m!(n−α)+1(m−1)!(n−α)(n−α−1)+⋯ +12!(n−α)⋯(n−α−m+3)(n−α−m+2) +11!(n−α)⋯(n−α−m+3)(n−α−m+2)(n−α−m) = ⋯=1m!(n−α−m).

Then, we have

 ∫∞ϵxα−1−nf(x)dx=n−1∑m=0ϵα−n+mm!(n−α−m)f(m)(0)+O(1)as ϵ↓0,

and the limit of (3) exists and is finite.

As shown in the following theorem, a f.p. integral (3) is expressed using a complex integral, which is the bases of our numerical method.

###### Theorem 1

We suppose that is analytic in a domain containing the positive real axis in its interior. Then, we have

 I(n,α)[f]=(−1)n+12isinπα∮C(−z)α−1−nf(z)dz, (4)

where is a complex integral path such that it encircles the positive real axis in the positive sense and it is contained in , and is the principal value, that is, the branch such that it takes a real value on the positive real axis222The complex integral on the right-hand side of (4) coincides with the integral of as a hyperfunction [3]. .

### Proof of Theorem 1

By Cauchy’s integral theorem, we have

 ∮C(−z)α−1−nf(z)dz=(∫Cϵ+∫Γ(+)ϵ+∫Γ(−)ϵ)(−z)α−1−nf(z)dz,

where , and and are the complex integral paths respectively given by

 Γ(+)ϵ= {x+i0|+∞>x≧ϵ}, Γ(−)ϵ= {x−i0|ϵ≦x<+∞}, Cϵ= {ϵeiθ|0≦θ≦2π}

(see Figure 1).

Regarding the integrals on , we have

 (∫Γ(+)ϵ+∫Γ(−)ϵ)(−z)α−1−nf(z)dz = = −∫∞ϵ(xe−iπ)α−1−nf(x)dx+∫∞ϵ(xeiπ)α−1−nf(x)dx = 2i(−1)n+1sinπα∫∞ϵxα−1−nf(x)dx.

Regarding the integral on , we have

 ∫Cϵ(−z)α−1−nf(z)dz = ∫2π0(ϵei(θ−π))α−1−nf(ϵeiθ)iϵeiθdθ = = i(−1)n+1ϵα−ne−iπα∫2π0ei(α−n)θ∞∑k=0ϵkk!f(k)(0)eikθdθ = i(−1)n+1e−iπα∞∑k=0ϵα−n+kk!f(k)(0)∫2π0ei(α−n+k)θdθ = −2i(−1)n+1sinπαn−1∑k=0ϵα−n+kk!(n−α−k)f(k)(0)+o(1)as ϵ↓0,

where we exchanged the order of the integration and the infinite sum on the fourth equality since the infinite sum is uniformly convergent on . Summarizing the above calculations, we have

 ∮C(−z)α−1−nf(z)dz = 2i(−1)n+1sinπα{∫∞ϵxα−1−nf(x)dx−n−1∑k=0ϵα−n+kk!(n−α−k)f(k)(0)} +o(1)as ϵ↓0,

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

We obtain an approximation formula of computing the f.p. integral (3) by evaluating the complex integral in (4) by the DE formula [11]

 ∫∞−∞g(u)du= ∫∞−∞g(ψDE(v))ψ′DE(v)dv ≃ hN+∑k=−N−g(ψDE(kh))ψ′DE(kh), (5)

where is the DE transform

 ψDE(v)={sinh(sinhv)if g(u)=O(|u|−1−α)as u→±∞ (α>0)sinhvif g(u)=exp(−c|u|)as u→±∞ (c>0). (6)

We can take the positive integers small for a given mesh since the transformed integrand decays double exponentially as . Taking a parameterization of the integral path

 C:z=φ(u),−∞

and evaluating the complex integral of (4) by the DE formula, we obtain the approximation formula of the f.p. integral

 I(n,α)[f]= (−1)n+12isinπα∫∞−∞(−φ(u))α−1−nf(φ(u))φ′(u)du ≃ I(n,α)h,N+,N−[f] = (−1)n+1h2isinπαN+∑k=−N−(−φ(ψDE(kh)))α−1−nf(φ(ψDE(kh))) ×φ′(ψDE(kh))ψ′DE(kh). (7)

A theoretical error estimate of the approximation (7) is given in the following theorem, where are taken to be for the simplicity.

###### Theorem 2

We suppose that

• is an analytic function in the strip

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

and the domain

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

is contained in ,

• satisfies

 N(n,α)(f,φ,ψD,Dd) ≡ < ∞,

where

 Dd(ϵ)={w∈C||Rew|<1/ϵ,|Imw|

and

• there exists positive numbers and such that

 ∣∣(−φ(ψDE(v)))α−1−nf(φ(ψDE(u)))φ′(ψDE(v))ψ′DE(v)∣∣≦Cexp(−c1exp(c2|v|))(∀v∈R).

Then, we have the inequality

 |I(n,α)[f]−I(n,α)h,N[f]| ≦ 12πN(n,α)(f,φ,ψDE,Dd)exp(−2πd/h)1−exp(−2πd/h) +C(n,α)(f,φ,ψDE,Dd)exp(−c1exp(c2Nh)), (8)

where

 I(n,α)h,N[f]=I(n,α)h,N,N[f]

and is a positive number depending on and only.

This theorem shows that the proposed approximation (7) converges exponentially as the mesh decreases and the number of sampling points increases.

### Proof of Theorem 2

We have

 |I(n,α)[f]−I(n,α)h,N[f]| ≦ |I(n,α)[f]−I(n,α)h[f]| +∣∣ ∣∣h2πi∑|k|>N(−φ(ψDE(kh)))α−1−nf(φ(ψDE(kh)))φ′(ψDE(kh))ψDE(kh)∣∣ ∣∣, (9)

where . Regarding the first term on the right-hand side of (9), from (4) and Theorem 3.2.1 of [10], we have

 |the first term|≦12πN(f,φ,ψDE,Dd)exp(−2πd/h)1−exp(−2πd/h).

Regarding the second term, we have

 |the second term|≦ Ch∑|k|>Nexp(−c1exp(c2kh)) ≦ 2C∫∞Nhexp(−c1exp(c2x))dx ≦ 2C∫∞Nhexp(c2x)exp(−c1exp(c2x))dx = 2Cc1c2exp(−c1exp(c2Nh)).

Therefore, we have (8).

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. Then, taking the integral path to be symmetric with respect to the real axis, that is,

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

which leads to , we have

 I(n,α)[f]≃ I(n,α)h,N[f] = hsinπαIm{12(−φ(ψDE(0)))α−1−nf(φ(ψDE(0)))φ′(ψDE(0))ψ′DE(0) +N∑k=1(−φ(ψDE(kh)))α−1−nf(φ(ψDE(kh)))φ′(ψDE(kh))ψ′DE(kh)} (10)

## 3 Numerical examples

We computed the f.p. integrals

 (i) f.p.∫∞0xα−1−n1+x2dx={(−1)m(π/2)/sin(πα/2)n=2m (even)% (−1)m+1(π/2)/cos(πα/2)n=2m+1 (odd), (11) (ii) f.p.∫∞0xα−1−ne−xdx=(−1)nΓ(α)(1−α)(2−α)⋯(n−α)

with by the proposed approximation formula (10). We performed all the computations using programs coded in C++ with double precision working. We took the complex integral path as

 C:w=φ(u)=u+0.5iiπlog(1+i(u+0.5i)1−i(u+0.5i)),+∞>u>−∞

(see Figure 2). We decided the number of sampling points for given mesh by truncating the infinite sum of the right-hand side of (10) at the -th term satisfying

 hsinπα|the transformed integrand of the k-th term|<10−15×∣∣I(n,α)h,N[f]∣∣.

Figure 3 shows the relative errors of the proposed method (10) applied to the f.p. integrals (11). From these figures, the proposed formula converges exponentially as the number of sampling points increases.

## 4 Summary

In this paper, we proposed a numerical method of computing a f.p. integral with a non-integral power singularity at the endpoint over a half infinite interval. In the proposed method, we express the desired f.p. integral by a complex integral, and we obtain the f.p. integral by evaluating the complex integral by the DE formula. Theoretical error estimate and some numerical examples show that the proposed approximation converges exponentially as the mesh of the DE formula decreases and the number of sampling points increases for an analytic integrand.

The complex integral which expresses a desired f.p. integral and gives the basis of the proposed method coincides with the definition of the integral of a hyperfunction, a generalized function given by an analytic function. In hyperfunction theory [3], a hyperfunction is described by an analytic function called a defining function, and its integral is defined by an complex integral involving the defining function. In addition, in hyperfunction theory, we can deal with an ordinary integral and a f.p. integral in a unified way. We found the proposed method, that is, the method of computing a f.p. integral by evaluating a complex integral, from this viewpoint in hyperfunction theory. Therefore, we expect that hyperfunction theory is applicable to many numerical computations in science and engineering.

## References

• [1] 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.
• [2] 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.
• [3] U. Graf (2010) Introduction to hyperfunctions and their integral transforms — an applied and computational approach. Birkha̋user, Basel. Cited by: §4, footnote 2.
• [4] H. Ogata and H. Hirayama (2018) Numerical integration based on hyperfunction theory. J. Comput. Appl. Math. 327, pp. 243–259. Cited by: §1.
• [5] 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.
• [6] H. Ogata (2019) A numerical method for computing hadamard finite-part integrals with a non-integral power singularity at an endpoint. Note: arXiv:1909.11398v1 [math.NA] Cited by: §1.
• [7] 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.
• [8] H. Ogata (2019) A numerical method of computing hadamard finite-part integrals with an integer power singularity at the endpoint on a half infinite interval. Note: arXiv:1910.00807v1 [math.NA] Cited by: §1.
• [9] D. F. Paget (1981) The numerical evaluation of hadamard finite-part integrals. Numer. Math. 36, pp. 447–453. External Links: Document Cited by: §1.
• [10] F. Stenger (1993) Numerical methods based on sinc and analytic functions. Springer-Verlag, New York. Cited by: §2.
• [11] 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, §1, §2.