A numerical method of computing Hadamard finite-part integrals with an integral power singularity at the endpoint on a half infinite interval
aa r X i v : . [ m a t h . NA ] O c t A numerical method of computing Hadamardfinite-part integrals with an integral powersingularity at the endpoint on a half infiniteinterval
Hidenori Ogata ∗ October 3, 2019
Abstract
In this paper, we propose a numerical method of computing Hadamardfinite-part integrals with an integral power singularity at the endpointon a half infinite interval, that is, a finite value assigned to a divergentintegral with an integral power singularity at the endpoint on a half infiniteinterval. In the proposed method, we express a desired finite-part integralusing a complex integral, and we obtain the integral by evaluating thecomplex integral by the DE formula. Theoretical error estimate and somenumerical examples show the effectiveness of the proposed method.
The integral Z ∞ x − f ( x )d x, where f ( x ) is an analytic function on the half infinite interval [0 , + ∞ ) such that f (0) = 0 and f ( x ) = O( x − α ) with α >
0, is divergent. However, we can assigna finite value to this divergent integral. In fact, we consider the integral Z ∞ ǫ x − f ( x )d x ∗ Department of Computer and Network Engineering, Graduate School of Informatics andEngineering, The University of Electro-Communications, 5 1-5-1 Chofugaoka, Chofu, Tokyo182-8585, Japan, (e-mail) [email protected] ǫ >
0, and, by integrating by part, we have Z ∞ ǫ x − f ( x )d x = Z ∞ ǫ (log x ) ′ f ( x )d x = (cid:20) log xf ( x ) (cid:21) ∞ ǫ − Z ∞ ǫ log xf ′ ( x )d x = − f ( ǫ ) log ǫ − Z ∞ ǫ log xf ′ ( x )d x = − f (0) log ǫ − Z ∞ log xf ′ ( x )d x + O(1) ( as ǫ ↓ . Then, the limit lim ǫ ↓ (cid:26)Z ∞ ǫ x − f ( x )d x + f (0) log ǫ (cid:27) exists and is finite. We call this limit the Hadamard finite-part (f.p.) integraland denote it by f . p . Z ∞ x − f ( x )d x. In general, we can define the f.p. integral I ( n ) [ f ] = f . p . Z ∞ x − n f ( x )d x ( n = 1 , , . . . ) , (1)where f ( x ) is an analytic function such that f (0) = 0 and f ( x ) = O( x n − α − )( x → + ∞ ) with α > . p . Z f ( x )( x − λ ) n d x ( 0 < λ < , n = 1 , , . . . ) , (2)many methods were proposed. Elliot and Paget proposed Gauss-type numericalintegration formulas for (2) [3, 9]. Bialecki proposed Sinc numerical integrationformulas for (2) [1, 2], where the trapezoidal formula is used together with vari-able transform technique as in the DE formula [11]. Ogata and et al. improvedthem and proposed a DE-type numerical integration formula for (2) [8].2he remainder of this paper is structured as follows. In Section 2, we definethe f.p. integral (1) and propose a numerical method of computing it. Inaddition, we show theoretical error estimate of the proposed method. In Section3, we show some numerical example which show the effectiveness of the proposedmethod. In Section 4, we give a summary of this paper. The f.p. integral (1) is defined by I ( n ) [ f ] = f . p . Z ∞ x − n f ( x )d x = lim ǫ ↓ (Z ∞ ǫ x − n f ( x )d x − n − X k =0 ǫ k +1 − n k !( n − − k ) f ( k ) (0) + log ǫ ( n − f ( n − (0) ) ( n = 1 , , . . . ) , (3)where f ( x ) is an analytic function on [0 , + ∞ ) such that f (0) = 0 and f ( x ) =O( x n − − α ) as x → + ∞ with α >
0, and the second term on the right-handside is zero if n = 1. We can show that it is well-defined as follows. In fact, for ǫ >
0, we can show by integrating by part Z ∞ ǫ x − n f ( x )d x = ǫ − n n − f ( ǫ ) + 1 n − Z ∞ ǫ x − n f ′ ( x )d x = ǫ − n n − f ( ǫ ) + ǫ − n ( n − n − f ′ ( ǫ ) + 1( n − n − Z ∞ ǫ x − n f ′′ ( x )d x = · · · = n − X k =0 ǫ k +1 − n ( n − n − · · · ( n − − k ) f ( k ) ( ǫ ) − log ǫ ( n − f ( n − ( ǫ )+ 1( n − Z ∞ ǫ log xf ( n ) ( x )d x = n − X k =0 ǫ k +1 − n ( n − n − · · · ( n − − k ) ( ∞ X l =0 f ( k + l ) (0) l ! ǫ l ) − log ǫ ( n − ∞ X k =0 f ( k ) (0) k ! ǫ k + 1( n − Z ∞ ǫ log xf ( n ) ( x )d x = n − X l =0 (cid:26) l X k =0 l − k )!( n − n − · · · ( n − − k ) | {z } ( ∗ ) (cid:27) ǫ l +1 − n f ( l ) (0) − log ǫ ( n − f ( n − (0) + O(1) ( as ǫ ↓ , ∗ ) = 1 l !( n −
1) + 1( l − n − n − · · · + 11!( n − n − · · · ( n − l + 1)( n − l )+ 1( n − n − · · · ( n − l + 1)( n − l )( n − l − l !( n −
1) + 1( l − n − n −
2) + · · · + 12!( n − n − · · · ( n − l + 2)( n − l + 1)+ 11!( n − n − · · · ( n − l + 2)( n − l + 1)( n − l − l !( n −
1) + 1( l − n − n −
2) + · · · + 13!( n − n − · · · ( n − l + 3)( n − l + 2)+ 12!( n − n − · · · ( n − l + 3)( n − l + 2)( n − l − · · · = 1 l !( n − l − . Then, we have Z ∞ ǫ x − n f ( x )d x = n − X l =0 ǫ l +1 − n l !( n − l − f ( l ) (0) − log ǫ ( n − f ( n − (0)+ O(1) ( as ǫ ↓ . 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 f ( z ) is analytic in a complex domain D , whichcontains the half infinite interval [0 , + ∞ ) in its interior. Then, the f.p. integral(3) is expressed as I ( n ) [ f ] = 12 π i I C z − n f ( z ) log( − z )d z, (4) where C is a complex integral path such that it encircles [0 , + ∞ ) in the positivesense and is contained in D . Proof of Theorem 1
From Cauchy’s integral theorem, we have12 π i I C z − n f ( z ) log( − z )d z = 12 π i (cid:18)Z Γ (+) ǫ + Z C ǫ + Z Γ ( − ) ǫ (cid:19) z − n f ( z ) log( − z )d z, ǫ >
0, where Γ (+) ǫ and C ǫ are complex integral paths respectively defined byΓ (+) ǫ = { x ∈ R + i0 | + ∞ > x ≧ ǫ } , Γ ( − ) ǫ = { x ∈ R − i0 | ǫ ≦ x < + ∞ } ,C ǫ = { ǫ e i θ | ≦ θ ≦ π } (see Figure 1), and the complex logarithmic function log z is the principal value,that is, the branch such that it takes a real value on the positive real axis. AsPSfrag replacements O C ǫ ǫ Γ (+) ǫ Γ ( − ) ǫ Figure 1: The complex integral paths Γ ( ± ) ǫ and C ǫ .to the integrals on Γ ( ± ) ǫ , we have12 π i (cid:18)Z Γ (+) ǫ + Z Γ ( − ) ǫ (cid:19) z − n f ( z ) log( − z )d z = 12 π i (cid:26) − Z ∞ ǫ x − n f ( x ) log( x e − i π )d x + Z ∞ ǫ x − n f ( x ) log( x e i π )d x (cid:27) = 12 π i (cid:26) − Z ∞ ǫ x − n f ( x )(log x − i π )d x + Z ∞ ǫ x − n f ( x )(log x + i π )d x (cid:27) = Z ∞ ǫ x − n f ( x )d x. As to the integral on C ǫ , we have12 π i Z C ǫ z − n f ( z ) log( − z )d z = 12 π i Z π ( ǫ e i θ ) − n f ( ǫ e i θ ) log( ǫ e i( θ − π ) )i ǫ e i θ d θ = ǫ − n π Z π e i(1 − n ) θ ( ∞ X k =0 f ( k ) (0) k ! ǫ k e i kθ ) { log ǫ + i( θ − π ) } d θ = 12 π ∞ X k =0 ǫ k − n +1 k ! log ǫf ( k ) (0) Z π e i( k − n +1) θ d θ + i2 π ∞ X k =0 ǫ k − n +1 k ! f ( k ) (0) Z π ( θ − π )e i( k − n +1) θ d θ, (5)5here we exchanged the order of the integral and the infinite sum since theinfinite series is uniformly convergent on 0 ≦ θ ≦ π . Since Z π e i( k − n +1) θ d θ = 2 πδ k,n − , Z π ( θ − π )e i( k − n +1) θ d θ = − π i( n − − k ) ( 0 ≦ k ≦ n − k = n − , we have(5) = log ǫ ( n − f ( n − (0) − n − X k =0 ǫ k − n +1 k !( n − − k ) f ( k ) (0) + O( ǫ ) ( ǫ ↓ . Summarizing the above calculations, we have12 π i I C z − n f ( z ) log( − z )d z = Z ǫ x − n f ( x )d x − n − X k =0 ǫ k − n +1 k !( n − − k ) f ( k ) (0)+ log ǫ ( n − f ( n − (0) + O( ǫ ) ( ǫ ↓ , and, taking the limit ǫ ↓
0, 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 DEformula [11], that is, Z ∞−∞ g ( u )d u = Z ∞−∞ g ( ψ DE ( v )) ψ ′ DE ( v )d v ≃ h N + X k = − N − g ( ψ DE ( kh )) ψ ′ DE ( kh ) , where h > u = ψ DE ( v ) is the DEtransform ψ DE ( v ) = ( sinh(sinh v ) ( g ( u ) = O( | u | − α − ) as u → ±∞ , α > v ( g ( u ) = O(exp( − c | u | )) as u → ±∞ , c > , and N ± is a positive integer such that the transformed integrand g ( ψ DE ( kh )) ψ ′ DE ( kh ) is sufficiently small at k = − N − , N + . We can take N ± small since g ( ψ DE ( v )) ψ ′ DE ( v ) decays double exponentially as v → ±∞ . Then,we have the approximation formula I ( n ) [ f ] ≃ I ( n ) h,N + ,N − [ f ] ≡ h π i N + X k = − N − ϕ ( ψ DE ( kh )) − n f ( ϕ ( ψ DE ( kh ))) log( − ϕ ( ψ DE ( kh ))) × ϕ ′ ( ψ DE ( kh )) ψ ′ DE ( kh ) , (6)where z = ϕ ( u ) , −∞ < u < + ∞ is a parameterization of the complex integralpath C .If f ( x ) is an analytic function on the real axis and C is an analytic curve, theproposed approximation (6) converges exponentially as shown in the followingtheorem. For the simplicity, we take N + = N − ≡ N ′ .6 heorem 2 We suppose that1. the parameterization function ϕ ( w ) of C is analytic in the strip D d = { w ∈ C | | Im w | < d } ( d > such that ϕ ( D d ) = { ϕ ( w ) | w ∈ D d } , is contained in C \ [0 , + ∞ ] ,2. N ( f, ϕ, ψ DE , D d ) ≡ lim ǫ → I ∂ D d ( ǫ ) | ϕ ( ψ DE ( w )) − n f ( ϕ ( ψ DE ( w ))) log( − ϕ ( ψ DE ( w ))) ψ ′ DE ( w ) | < ∞ , where D d ( ǫ ) ≡ { w ∈ C | | Re w | < /ǫ, | Im w | < d (1 − ǫ ) } . and3. there exist positive numbers C , c and c such that | f ( ϕ ( ψ DE ( v ))) | ≦ C exp( − c exp( c | v | )) ( ∀ v ∈ R ) . Then, we have the inequality | I ( n ) [ f ] − I ( n ) h,N ′ [ f ] | ≦ π N ( f, ϕ, ψ DE , D d ) exp( − πd/h )1 − exp( − πd/h )+ C ( f, ϕ, ψ DE , D d ) exp( − c exp( c N ′ h )) , (7) where I ( n ) h,N ′ [ f ] = I ( n ) h,N ′ ,N ′ [ f ] and C ( f, ϕ, ψ DE , D d ) is a positive number dependingon f ( z ) , ϕ , ψ DE and D d only. This theorem shows that the approximation formula (6) converges exponentiallyas the mesh h decreases and the number of sampling points 2 N ′ + 1 increases. Proof of Theorem 2
We have (cid:12)(cid:12)(cid:12)(cid:12) f . p . Z ∞ x − n f ( x )d x − I ( n ) h,N [ f ] (cid:12)(cid:12)(cid:12)(cid:12) ≦ (cid:12)(cid:12)(cid:12)(cid:12) f . p . Z ∞ x − n f ( x )d x − I ( n ) h [ f ] (cid:12)(cid:12)(cid:12)(cid:12) + h X | k | >N (cid:12)(cid:12) ϕ ( ψ DE ( kh )) − n f ( ϕ ( ψ DE )( kh )) ϕ ′ ( ψ DE ( kh )) ψ ′ DE ( kh ) (cid:12)(cid:12) , (8)where I ( n ) h [ f ] = lim N →∞ I ( n ) h,N [ f ]. For the first term on the right-hand side of(8), we have (cid:12)(cid:12)(cid:12)(cid:12) f . p . Z ∞ x − n f ( x )d x − I ( n ) h [ f ] (cid:12)(cid:12)(cid:12)(cid:12) ≦ π N ( f, ϕ, ψ DE , D d ) exp( − πd/h )1 − exp( − πd/h )7y Theorem 3.2.1 in [10]. For the second term on the right-side hand, we have | the second term | ≦ C h X | k | >N exp( − c exp( c kh )) ≦ C Z ∞ kh exp( − c exp( c x ))d x ≦ C Z ∞ kh exp( c x ) exp( − c exp( c x ))d x = 2 C c exp( − c exp( c N h )) . Then, we obtain (7).We remark here that we can reduce the number of sampling points by half ifthe integrand f ( x ) is real valued on the real axis. In fact, we have f ( z ) = f ( z )by the reflection principle, taking the integral path C symmetric with respectto the real axis, that is, ϕ ( − u ) = ϕ ( u ), which leads to ϕ ′ ( − u ) = − ϕ ′ ( u ) , and taking the DE transform ψ DE ( v ) to be an even function, we have I ( n ) [ f ] ≃ I ( n ) ′ h,N [ f ] ≡ h π Im (cid:8) ϕ ( ψ DE (0)) − n f ( ϕ ( ψ DE (0))) log( − ϕ ( ψ DE (0))) ϕ ′ ( ψ DE (0)) ψ ′ DE (0) (cid:9) + hπ Im (cid:26) N X k =1 ϕ ( ψ DE ( kh )) − n f ( ϕ ( ψ DE ( kh ))) log( − ϕ ( ψ DE ( kh ))) × ϕ ′ ( ψ DE ( kh )) ψ ′ DE ( kh ) (cid:27) . (9) In this section, we show some numerical examples which show the effectivenessof the proposed method.We computed the f.p. integrals(i) f . p . Z ∞ x − n x d x = ( ( π/ − m ( n = 2 m (even) )0 ( n = 2 m + 1 (odd) )(ii) f . p . Z ∞ x − n e − x d x = − γ ( n = 1 ) − γ ( n = 2 ) − γ ( n = 3 ) − + γ ( n = 4 ) (10)for n = 1 , , ,
4, where γ is Euler’s constant, by the formula (9). All the com-putations were performed using programs coded in C++ with double precisionworking. The complex integral path C in (4) is taken as C : z = ϕ ( u ) = u + 0 . π log (cid:18) u + 0 . − i( u + 0 . (cid:19) , + ∞ > u > −∞ N for given mesh h = 2 − , − , . . . by truncating the infinite sum at the k -th term such that hπ × | the k -th term | < ( − × | I ( n ) ′ h,N [ f ] | if I ( n ) ′ h,N [ f ] = 010 − otherwise . Figure 3 shows the relative errors of the proposed approximation formula (9) -3-2-1 0 1 2 3 -1 0 1 2 3 4 5
PSfrag replacements Re z I m z Figure 2: The complex integral path C . ε ( 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 conver-gence of the proposed formula as the number of sampling points N increases. PSfrag replacements Nε ( n ) N [ f ] PSfrag replacements Nε ( n ) N [ f ]integral (i) integral (ii)Figure 3: The errors of the proposed approximation formula (9) applied to thef.p. integrals (10). 9 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-infiniteinterval. In the proposed method, we express the desired f.p. integral using acomplex integral, and we obtain the f.p. integral by evaluating the complexintegral by the DE formula. Theoretical error estimate and some numericalexamples show the exponential convergence of the proposed method in the casethat the integrand is an analytic function.
References [1] B. Bialecki. A sinc-hunter quadrature rule for cauchy principal value inte-grals.
Math. Comput. , 55:665–681, 1990.[2] B. Bialecki. A sinc quadrature rule for hadamard finite-part integrals.
Numer. Math. , 57:263–269, 1990.[3] D. Elliot and D. F. Paget. Gauss type quadrature rules for cauchy principalvalue integrals.
Math. Comput. , 33:301–309, 1979.[4] R. Estrada and R. P. Kanwal. Regularization, pseudofunction, andhadamard finite part.
J. Math. Anal. Appl. , 141:195–207, 1989.[5] H. Ogata. A numerical method for computing hadamard finite-partintegrals with a non-integral power singularity at an endpoint, 2019.arXiv:1909.11398v1 [math.NA].[6] H. Ogata. A numerical method for hadamard finite-part integrals withan integral power singularity at an endpoint, 2019. arXiv:1909.08872v1[math.NA].[7] H. Ogata and H. Hirayama. Numerical integration based on hyperfunctiontheory.
J. Comput. Appl. Math. , 327:243–259, 2018.[8] H. Ogata, M. Sugihara, and M. Mori. De-type quadrature formulae forcauchy principal-value integrals and for hadamard finite-part itnegrals. In
Proceedings of the Second ISAAC Congress , volume 1, pages 357–366, 2000.[9] D. F. Paget. The numerical evaluation of hadamard finite-part integrals.
Numer. Math. , 36:447–453, 1981.[10] F. Stenger.
Numerical Methods Based on Sinc and Analytic Functions .Springer-Verlag, New York, 1993.[11] H. Takahasi and M. Mori. Double exponential formulas for numerical in-tegration.