# 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 finite-part integrals with an integral power singularity at the endpoint on a half infinite interval, that is, a finite value assigned to a divergent integral with an 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 integral by evaluating the complex integral by the DE formula. Theoretical error estimate and some numerical examples show the effectiveness of the proposed method.

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/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...
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...
12/19/2020

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

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

### Matroid-Based TSP Rounding for Half-Integral Solutions

We show how to round any half-integral solution to the subtour-eliminati...
07/04/2021

### Interval probability density functions constructed from a generalization of the Moore and Yang integral

Moore and Yang defined an integral notion, based on an extension of Riem...

## 1 Introduction

The integral

 ∫∞0x−1f(x)dx,

where is an analytic function on the half infinite interval such that and with , is divergent. However, we can assign a finite value to this divergent integral. In fact, we consider the integral

 ∫∞ϵx−1f(x)dx

with , and, by integrating by part, we have

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

Then, the limit

 limϵ↓0{∫∞ϵx−1f(x)dx+f(0)logϵ}

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

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

In general, we can define the f.p. integral

 I(n)[f]=f.p.∫∞0x−nf(x)dx(n=1,2,…), (1)

where is an analytic function such that and with [4].

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

Previous works related to this paper are as follows. The author and Hirayama proposed a numerical method of computing ordinary integrals related to hyperfunction theory [5], a theory of generalized functions based on complex function theory. The author proposed numerical methods for computing Hadamard finite-part integrals with a singularity at an endpoint on a finite interval [8, 7]. In these methods, we express a desired integral using a complex integral, we obtain the integral by evaluating the complex integral by conventional numerical integration formulas. For Cauchy principal-value integrals or Hadamard finite-part integrals on a finite interval with a singularity in the interior of the integral interval

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

many methods were proposed. Elliot and Paget proposed Gauss-type numerical integration formulas for (2) [3, 9]. Bialecki proposed Sinc numerical integration formulas for (2) [2, 1], where the trapezoidal formula is used together with variable transform technique as in the DE formula [11]. Ogata and et al. improved them and proposed a DE-type numerical integration formula for (2) [6].

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 theoretical error estimate of the proposed method. In Section 3, we show some numerical example which show the effectiveness of the proposed method. In Section 4, we give a summary of this paper.

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

The f.p. integral (1) is defined by

 (3)

where is an analytic function on such that and as with , and the second term on the right-hand side is zero if . We can show that it is well-defined as follows. In fact, for , we can show by integrating by part

 ∫∞ϵx−nf(x)dx = ϵ1−nn−1f(ϵ)+1n−1∫∞ϵx1−nf′(x)dx = ϵ1−nn−1f(ϵ)+ϵ2−n(n−1)(n−2)f′(ϵ)+1(n−1)(n−2)∫∞ϵx2−nf′′(x)dx = ⋯ = n−2∑k=0ϵk+1−n(n−1)(n−2)⋯(n−1−k)f(k)(ϵ)−logϵ(n−1)!f(n−1)(ϵ) +1(n−1)!∫∞ϵlogxf(n)(x)dx = n−2∑k=0ϵk+1−n(n−1)(n−2)⋯(n−1−k){∞∑l=0f(k+l)(0)l!ϵl} −logϵ(n−1)!∞∑k=0f(k)(0)k!ϵk+1(n−1)!∫∞ϵlogxf(n)(x)dx = n−2∑l=0{l∑k=01(l−k)!(n−1)(n−2)⋯(n−1−k)(∗)}ϵl+1−nf(l)(0) −logϵ(n−1)!f(n−1)(0)+O(1)(as ϵ↓0),

and

 (∗)= 1l!(n−1)+1(l−1)!(n−1)(n−2) +⋯+11!(n−1)(n−2)⋯(n−l+1)(n−l) +1(n−1)(n−2)⋯(n−l+1)(n−l)(n−l−1) = 1l!(n−1)+1(l−1)!(n−1)(n−2)+⋯ +12!(n−1)(n−2)⋯(n−l+2)(n−l+1) +11!(n−1)(n−2)⋯(n−l+2)(n−l+1)(n−l−1) = 1l!(n−1)+1(l−1)!(n−1)(n−2)+⋯ +13!(n−1)(n−2)⋯(n−l+3)(n−l+2) +12!(n−1)(n−2)⋯(n−l+3)(n−l+2)(n−l−1) = ⋯ = 1l!(n−l−1).

Then, we have

 ∫∞ϵx−nf(x)dx= n−2∑l=0ϵl+1−nl!(n−l−1)f(l)(0)−logϵ(n−1)!f(n−1)(0) +O(1)(as ϵ↓0).

Therefore, the limit in (3) exists and is finite.

The f.p. integral (3) is expressed using a complex integral.

###### Theorem 1

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

 I(n)[f]=12πi∮Cz−nf(z)log(−z)dz, (4)

where is a complex integral path such that it encircles in the positive sense and is contained in .

### Proof of Theorem 1

From Cauchy’s integral theorem, we have

 12πi∮Cz−nf(z)log(−z)dz=12πi(∫Γ(+)ϵ+∫Cϵ+∫Γ(−)ϵ)z−nf(z)log(−z)dz,

for , where and are complex integral paths respectively defined by

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

(see Figure 1), and the complex logarithmic function is the principal value, that is, the branch such that it takes a real value on the positive real axis.

As to the integrals on , we have

 12πi(∫Γ(+)ϵ+∫Γ(−)ϵ)z−nf(z)log(−z)dz = 12πi{−∫∞ϵx−nf(x)log(xe−iπ)dx+∫∞ϵx−nf(x)log(xeiπ)dx} = 12πi{−∫∞ϵx−nf(x)(logx−iπ)dx+∫∞ϵx−nf(x)(logx+iπ)dx} = ∫∞ϵx−nf(x)dx.

As to the integral on , we have

 12πi∫Cϵz−nf(z)log(−z)dz = 12πi∫2π0(ϵeiθ)−nf(ϵeiθ)log(ϵei(θ−π))iϵeiθdθ = = 12π∞∑k=0ϵk−n+1k!logϵf(k)(0)∫2π0ei(k−n+1)θdθ +i2π∞∑k=0ϵk−n+1k!f(k)(0)∫2π0(θ−π)ei(k−n+1)θdθ, (5)

where we exchanged the order of the integral and the infinite sum since the infinite series is uniformly convergent on . Since

 ∫2π0ei(k−n+1)θdθ=2πδk,n−1, ∫2π0(θ−π)ei(k−n+1)θdθ=⎧⎪⎨⎪⎩−2πi(n−1−k)(0≦k≦n−2)0(k=n−1),

we have

 (???)=logϵ(n−1)!f(n−1)(0)−n−2∑k=0ϵk−n+1k!(n−1−k)f(k)(0)+O(ϵ)(ϵ↓0).

Summarizing the above calculations, we have

 12πi∮Cz−nf(z)log(−z)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)+O(ϵ)(ϵ↓0),

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

We obtain the f.p. integral by evaluating the complex integral on the right-hand side of (4) by a conventional numerical integration formula such as the DE formula [11], that is,

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

where is the mesh of the trapezoidal formula, is the DE transform

 ψDE(v)={sinh(sinhv)(g(u)=O(|u|−α−1)asu→±∞,α>0)sinhv(g(u)=O(exp(−c|u|))asu→±∞,c>0),

and is a positive integer such that the transformed integrand
is sufficiently small at . We can take small since decays double exponentially as . Then, we have the approximation formula

 I(n)[f]≃ I(n)h,N+,N−[f] ≡ h2πiN+∑k=−N−φ(ψDE(kh))−nf(φ(ψDE(kh)))log(−φ(ψDE(kh))) ×φ′(ψDE(kh))ψ′DE(kh), (6)

where is a parameterization of the complex integral path .

If is an analytic function on the real axis and is an analytic curve, the proposed approximation (6) converges exponentially as shown in the following theorem. For the simplicity, we take .

###### Theorem 2

We suppose that

1. the parameterization function of is analytic in the strip

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

such that

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

is contained in ,

2.  N(f,φ,ψDE,Dd) ≡ limϵ→0∮∂Dd(ϵ)|φ(ψDE(w))−nf(φ(ψDE(w)))log(−φ(ψDE(w)))ψ′DE(w)| < ∞,

where

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

and

3. there exist positive numbers , and such that

 |f(φ(ψDE(v)))|≦C0exp(−c1exp(c2|v|))(∀v∈R).

Then, we have the inequality

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

where and is a positive number depending on , , and only.

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

### Proof of Theorem 2

We have

 ∣∣∣f.p.∫∞0x−nf(x)dx−I(n)h,N[f]∣∣∣ ≦ ∣∣∣f.p.∫∞0x−nf(x)dx−I(n)h[f]∣∣∣ +h∑|k|>N∣∣φ(ψDE(kh))−nf(φ(ψDE)(kh))φ′(ψDE(kh))ψ′DE(kh)∣∣, (8)

where . For the first term on the right-hand side of (8), we have

 ∣∣∣f.p.∫∞0x−nf(x)dx−I(n)h[f]∣∣∣≦12πN(f,φ,ψDE,Dd)exp(−2πd/h)1−exp(−2πd/h)

by Theorem 3.2.1 in [10]. For the second term on the right-side hand, we have

 |the second term|≦ C0h∑|k|>Nexp(−c1exp(c2kh)) ≦ 2C0∫∞khexp(−c1exp(c2x))dx ≦ 2C0∫∞khexp(c2x)exp(−c1exp(c2x))dx = 2C0c2exp(−c1exp(c2Nh)).

Then, we obtain (7).

We remark here that we can reduce the number of sampling points by half if the integrand is real valued on the real axis. In fact, we have by the reflection principle, taking the integral path symmetric with respect to the real axis, that is, , which leads to

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

and taking the DE transform to be an even function, we have

 I(n)[f]≃I(n)′h,N[f] ≡ h2πIm{φ(ψDE(0))−nf(φ(ψDE(0)))log(−φ(ψDE(0)))φ′(ψDE(0))ψ′DE(0)} +hπIm{N∑k=1φ(ψDE(kh))−nf(φ(ψDE(kh)))log(−φ(ψDE(kh))) ×φ′(ψDE(kh))ψ′DE(kh)}. (9)

## 3 Numerical examples

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

We computed the f.p. integrals

 (i) (10) (ii) f.p.∫∞0x−ne−xdx=⎧⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪⎩−γ(n=1)−1+γ(n=2)34−12γ(n=3)−1136+16γ(n=4)

for , where is Euler’s constant, by the formula (9). All the computations were performed using programs coded in C++ with double precision working. The complex integral path in (4) is taken as

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

(see Figure 2). We took the number of sampling points for given mesh by truncating the infinite sum at the -th term such that

 hπ×|the k-th term|<{10−15×|I(n)′h,N[f]|ifI(n)′h,N[f]≠010−15otherwise.

Figure 3 shows the relative errors of the proposed approximation formula (9)

 ε(n)N[f]=⎧⎨⎩|I(n)h,N[f]−I(n)[f]|/|I(n)[f]|(I(n)[f]≠0)|I(n)h,N[f]|(otherwise)

applied to the f.p. integrals (10). These figures shows the exponential convergence of the proposed formula as the number of sampling points increases.

## 4 Summary

In this paper, we proposed a numerical integration formula for Hadamard finite-part integrals with an integral power singularity at the endpoint on a half-infinite interval. In the proposed method, we express the desired f.p. integral using 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 the exponential convergence of the proposed method in the case that the integrand is an analytic function.

## 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] 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.
• [4] 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.
• [5] H. Ogata and H. Hirayama (2018) Numerical integration based on hyperfunction theory. J. Comput. Appl. Math. 327, pp. 243–259. Cited by: §1.
• [6] 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.
• [7] 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.
• [8] 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.
• [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.