Numerical methods for hypersingular integrals on the real line

Numerical methods for hypersingular integrals on the real line

Maria Carmela De Bonis^a·Donatella Occorsio^b

Communicated by G. Milovanovi´c

Abstract

In the present paper the authors propose two numerical methods to approximate Hadamard transforms of the type

H_p(f w_β,t) = Z

=

R

f(x)

(x−t)^p+1w_β(x)d x,

wherepis a nonnegative integer andw_β(x) =e^−|x|β,β >1, is a Freud weight. One of the procedures employed here is based on a simple tool like the “truncated” Gaussian rule conveniently modified to remove numerical cancellation and overflow phenomena. The second approach is a process of simultaneous approximation of the functions{Hk(f w_β,t)}^p_k=0. This strategy can be useful in the numerical treatment of hypersingular integral equations. The methods are shown to be numerically stable and convergent and some error estimates in suitable Zygmund-type spaces are proved. Numerical tests confirming the theoretical estimates are given. Comparisons of our methods among them and with other ones available in literature are shown.

1 Introduction

In the present paper we propose some global strategies to approximate finite part integrals (shortly FP integrals) Z

=

R

f(x)

(x−t)^p+1w_β(x)d x, t∈R,

wherepis a nonnegative integer,w_β(x) =e^−|x|β,β >1, is a Freud weight and the density function f is assumed in such a way that the integral exists in the Hadamard sense.

Integrals of this kind are of interest, for instance, in the solution of hypersingular BIE, which model many different Physical and Engineering problems (see[26]and the references therein,[8],[11],[18],[1],[31]).

The numerical treatment of FP integrals over finite ranges has been extensively treated[25],[16]and limiting ourselves to interpolatory rules based on the zeros of orthogonal polynomials, we cite among them[14],[29],[24]whenever the integrand contains Jacobi weights. For a comparison among different numerical approaches in[−1, 1]see also[32]. On the other hand the use of the tools introduced for finite ranges to integrals over infinite ranges can be inadequate. For instance, non linear transformation can produce slower convergent results, since the density function can be significantly worsened. Moreover, the analysis of the stability of the rules and the estimates of the errors in different functional spaces can present some additional difficulties. For instance, the appealing Gaussian rule looses its high performance as well as product integration rules based on the zeros of orthogonal polynomials, and by a truncation strategy introduced in[19](see also[7],[22]) the problem can be successfully overcome.

An additional hindrance depends on the weightw_β, since the coefficients of the three term recurrence relation of the corresponding orthogonal polynomials are unknown in the general case. However, it is possible to calculate them with high accuracy in extended arithmetic usingWolfram Mathematicaroutines[2](details are given in Section 5).

aDipartimento di Matematica ed Informatica, Università degli Studi della Basilicata, Viale dell’Ateneo Lucano 10, 85100 Potenza, Italy, GNCS member, email mariacarmela.debonis@unibas.it.

bDipartimento di Matematica ed Informatica, Università degli Studi della Basilicata, Viale dell’Ateneo Lucano 10, 85100 Potenza, Italy, GNCS member, email donatella.occorsio@unibas.it.

Following a very standard way, we start from the decomposition Hp(f w_β,t) :=

Z

=

R

f(x)

(x−t)^p+1w_β(x)d x

= Z

R

f(x)−Pp k=0

f^(k)(t) k! (x−t)^k

(x−t)^p+1 w_β(x)d x+

p

X

k=0

f^(k)(t) k!

Z

=

R

w_β(x) (x−t)^p+1−kd x

=: F_p(f,t) +

p

X

k=0

f^(k)(t)

k! H_p−k(w_β,t)

= 1

p!

d^p

d t^pF0(f,t) +

p

X

k=0

f^(k)(t)

k! H_p−k(w_β,t), (1)

where

F0(f,t) = Z

R

f(x)−f(t)

x−t w_β(x)d x, (2)

focusing the attention on the first right-hand integral, since the remaining ones can be efficiently computed by standard routines for many choices ofβ(see Section 5). In the first method, for any fixedt, we select a suitable subsequence of Gaussian rules based on Freud zeros, which leads to a stable and fast convergent procedure. By this way, the severe numerical cancellation arising whentis “close” to a Gaussian node is avoided, without changes of variables which double the samples of f and depend on the singularityt (see for instance[15],[25]). In addition, since we use “truncated” Gaussian rules[20], the number of function evaluations is drastically reduced and possible overflow phenomena are overcome. To complete the framework, fort

“large” (large in the sense we specify later) andf “sufficiently” smooth, we propose a shrewd use of the Gauss-Freud rule.

Since in some contexts the approximation of{Hk(f w_β,t)}^p_k=0for a “large” number oftand/or the uniform convergence of the approximation process are required, following an idea in[5](see also[28],[21],[24],[4]), we approximate{F^(k)₀ (f)}^p_k=0by means of the successive derivatives of a suitable Lagrange polynomial interpolatingF0(f)at Freud zeros. Since in the general case the samples ofF0(f)at the interpolation knots cannot be exactly computed, we approximate them by a “truncated” Gauss-Freud rule (see[20]). For a correct error estimate in weighted uniform spaces, we determine the class ofF0(f)depending on the Zygmund-type space f belongs to.

Moreover, whenever the derivatives{f^(k)(t)}^p_k=1are not available, this method can be performed avoiding the differentiation of the density function f. Indeed, in this case, we propose to approximate f^(k)(t)by the derivatives of a suitable Lagrange polynomial interpolating f, without requiring additional samples of the density functionf.

The paper is organized as follows. In Section 2 some notations and basic results needed to introduce the main results are collected. In Sections 3 and 4 we describe the numerical methods proposed to approximateHp(f w_β,t), getting some results about the stability and the rate of convergence. Section 5 contains some computational details useful in the implementation process and in Section 6 we present some numerical tests, by comparing both the procedures among them and with other methods available on the same topic. As we will show, the numerical results agree with the theoretical analysis and confirm the efficiency of the proposed procedures. Finally, in Section 7 we give the proofs of the main results.

2 Notations and preliminary results

In the sequelCwill denote any positive constant which can be different in different formulas. MoreoverC6=C(a,b, ..)will be used to say that the constantCis independent ofa,b, .... The notationA∼B, whereAandBare positive quantities depending on some parameters, will be used if and only if(A/B)^±1≤C, withCpositive constant independent of the above parameters.

Throughout the paperθ will denote a fixed real number belonging to the interval(0, 1)which can be different in different formulas. Moreover,Pmis the space of all algebraic polynomials of degree at mostm.

Let{p_m(w_β)}mbe the sequence of orthonormal polynomials w.r.t. the weightw_β(x) =e^−|x|β,β >1, with positive leading coefficients, i.e.

p_m(w_β,x) =γm(w_β)x^m+ terms of lower degree, γm(w_β)>0.

Denote byx_m,k,k=1, . . . ,_m

2

, the positive zeros ofp_m(w_β)and byx_m,−k=−x_m,kthe negative ones, settingx_m,0=0 form odd. The zeros ofp_m(w_β)satisfies

−am<x

m,−_m

2

<· · ·<xm,−1<· · ·<xm,1<xm,2<· · ·<x

m,_m

2

<a_m, wherea_m:=a_m(w_β)is the Mhaskar-Rachmanoff-Saff number w.r.t.w_β(in the sequel M-R-S number).

2.1 Function spaces and best approximation estimates Settingu(x) =e^−|x|

β

2 ,β >1, we define the function space C_u = n

f : f u∈C⁰(R), lim

x→±∞f(x)u(x) =0o ,

equipped with the norm

kfkC_u:=kf uk:=sup

x∈R|f(x)|u(x).

In the sequel we will writekf ukE=supx∈E|f(x)|u(x)for anyE⊂R. For smoother functions, we denote byW_s(u),s∈Z⁺, the Sobolev space

W_s(u) =

f ∈C_u:f^(s−1)∈AC(R):kf^(s)uk<∞ ,

whereAC(R)denotes the set of the functions which are absolutely continuous on every closed subset ofR, equipped with the norm

kfkW_s(u)=kf uk+kf^(s)uk.

Denoting by

E_m(f)u= inf

P∈Pmk(f−P)uk

the error of the best approximation inC_u, the following weaker version of the Jackson theorem holds E_m(f)u ≤ C

Z ^am_m

0

Ω^r(f,t)u

t d t, C6=C(m,f), (3)

wherea_m:=a_m(u)denotes the M-R-S number w.r.tuand Ω^r(f,t)u= sup

0<h≤tk(∆h^rf)ukI_rh,

is the main part of ther−th modulus of continuity (see[9]), withI_rh= [−Arh^−β−11 ,Arh^−β−11 ],A>0 a constant and

∆hf(x) =f

 x+h

2

‹

−f

 x−h

2

‹

, ∆_h^r=∆h(∆_h^r−1).

Sincea_m(u)∼a_m(w_β)∼m^β1, throughout the paper we employ the same symbola_mto denote both of them.

By means of the main part of ther−th modulus of continuity, the following Zygmund spaceZ_s(u), withs∈R⁺, can be defined (see[27])

Z_s(u) = §

f ∈C_u: sup

t>0

Ω^r(f,t)u

t^s <∞, r>s ª

, equipped with the norm

kfkZ_s(u)=kf uk+sup

t>0

Ω^r(f,t)u

t^s , r>s.

By (3) it is easy to deduce

E_m(f)u≤C a_m

m s

kfkW_s(u), ∀f ∈W_s(u), s∈Z⁺, C6=C(m,f), (4) E_m(f)u≤C

a_m m s

kfkZ_s(u), ∀f∈Z_s(u), s∈R⁺, C6=C(m,f). (5) 2.2 Truncated Lagrange interpolation at Freud zeros

For any fixed 0< θ <1, let

x_m,j:=x_m,j(m)=min

§

x_m,k:x_m,k≥a_mθ, k=1, 2, ..,

›m 2

žª

. (6)

Let L_m+2(w_β,g)be the Lagrange polynomial interpolating a given functiongat the zeros of p_m(w_β,x)(a²_m−x²)and denote by χjthe characteristic function of the interval(−x_m,j,x_m,j). The Lagrange polynomialL_m+2(w_β,g):=L_m+2(w_β,gχj), introduced in [20], can take the following expression

L_m+2(w_β,g,x) =

j

X

k=−j

l_m,k(x) a²_m−x²

a²_m−x²_m,kf(x_m,k) =:

j

X

k=−j

`m,k(x)f(x_m,k), (7)

wherel_m,k(x) = _p0 ^pm(wβ,x)

m(w_β,x_m,k)(x−x_m,k). The polynomialL_m+2(w_β,g)belongs toP_m+1^∗ ⊂Pm+1, with P^∗_m+1=

p∈Pm+1:p(x_m,k) =p(±a_m) =0, |k|>j , and the operatorL_m+2(w_β)projectsC_uontoP_m+^∗ ₁.

About the simultaneous approximation of a function and its successive derivatives, we state the following result.

Theorem 2.1. If f ∈Z_p+λ(u), p∈N,0< λ <1, we have k(f−L_m+2(w_β,f))^(k)uk ≤ CkfkZ_p+λ(u)logm

a_m m

p+λ−k

, k=1, . . . ,p, C6=C(m,f). In particular, if f ∈W_p+r(u), p∈N, r∈Z⁺,we get

k(f−L_m+2(w_β,f))^(k)uk ≤ CkfkW_p+r(u)logm a_m

m p+r−k

, k=1, . . . ,p, C6=C(m,f). (8)

2.3 Truncated Gauss-Freud rule

In[20](see also[7],[23]), it was introduced a truncated rule of the the following type Z

R

f(x)w_β(x)d x=X

|k|≤j

λm,kf(x_m,k) +ρm(f), (9)

whereλm,k,k=−_m

2

, . . . ,_m

2

, are the Christoffel numbers w.r.t. the weightw_β andρm(f)is the remainder term. We recall the following result useful in the next and proved in[20]in a more general context.

Theorem 2.2. For any f ∈C_u

|ρm(f)| ≤C E_M(f)u+e^−Amkf uk

, C6=C(m,f), A6=A(m,f), (10) where M=cm,0<c<1fixed.

3 The first method

Starting from the decomposition (1), lettbe fixed in the interval[−x_m,j,x_m,j], where the index jhas been defined in (6). Since there exist efficient routines to computeH_p−k(w_β,t)for some choices ofβ(see Section 5), we focus our attention to the term Fp(f,t)assumingf s.t. the integral exists as an improper Riemann integral. By using the truncated Gauss-Freud rule in (9) we have

Fp(f,t) =X

|k|≤j

λm,k

f(x_m,k)−Pp k=0

f^(k)(t)

k! (x_m,k−t)^p−k

(x_m,k−t)^p+1 +e_p,m(f,t) =:Fp,m(f,t) +e_p,m(f,t).

However, as observed also in[25], formulae of this kind are very appealing since its coefficients have a very simple form, but they can present severe numerical cancellation whentis “close” to one of the Gaussian nodes.

Thus, to overcome this problem we propose to select, for any fixedt, a proper subsequence of the truncated Gauss-Freud sequence{Fp,m(f,t)}m, in such a way that the distances|xm,k−t|,k=−j, . . . ,j, are always large enough.

Formsufficiently large (saym≥m₀∈N), there exists an indexds.t.x_m,d≤t<x_m,d+1. Since the zeros ofp_m(w_β)interlace those ofp_m+1(w_β), two cases are possible:

(1) x_m,d≤t<x_m+1,d+1, (2) x_m+1,d+1≤t<x_m,d+1. In the case (1), if

t∈

x_m,d, x_m+1,d+1+x_m,d 2

•

then we approximateF_p(f,t)withF_p,m+1(f,t)and we setm^∗=m+1, while if t∈

•x_m+1,d+1+x_m,d

2 , x_m+1,d+1

then we approximateFp(f,t)withFp,m(f,t)and we setm^∗=m. Similarly we proceed in the case (2).

By this way we have

Fp(f,t) := Fp,m∗(f,t) +e_p,m∗(f,t), whereep,m∗(f,t)is the remainder term and therefore

Hp(f w_β,t) := Fp,m∗(f,t) +

p

X

k=0

f^(k)(t)

k! Hp−k(w_β,t) +e_p,m∗(f,t). (11) In order to highlight the possible instability of the ordinary truncated Gaussian rule, we want to propose just an example.

Consider the integral

H1(f w_β),t) = Z

=

R

sin(x+5) (x−t)² e^−x2d x.

In Figure 1 we compare the absolute errors obtained by implementing the truncated Gauss-Freud rule{Fp,m(f,t)}m(red curve) and the modified Gauss-Freud rule{Fp,m^∗(f,t)}m(blue curve), for increasing values ofm.

As the graph shows, the highest errors are attained by using the ordinary truncated Gaussian rule and, for the samet, also for different values ofm. Since the functionf is smooth, this bad behavior depends on the closeness oftto some of the Gaussian abscissae. The shortcoming is successfully overcome by using the method in (11).

Next theorem deals with the stability and the convergence of the proposed quadrature rule.

Theorem 3.1. Let p∈N0,0< λ <1and0< θ <1fixed. For any fixed t∈(−θa_m,θa_m), if f ∈Z_p+λ(u)then

u(t)|Fp,m∗(f,t)| ≤ CkfkZ_p+λ(u)logm (12) and

u(t)|Fp(f,t)−F_p,m∗(f,t)| ≤ CkfkZ_p+λ(u)logm a_m

m λ

, (13)

where0<C6=C(m,f,t).

Figure 1:Comparison between the absolute errors of the rulesF_p,m(f,t)andF_p,m∗(f,t)

For smoother functions the following result holds.

Corollary 3.2. Let p∈N0, r∈N,0< λ <1and0< θ <1fixed. For any fixed t∈(−θa_m,θa_m), if f∈Z_p+r+λ(u)then u(t)|Fp(f,t)−Fp,m∗(f,t)| ≤ CkfkZ_p+r+λ(u)logm

a_m m

r+λ

, 0<C6=C(m,f,t), (14) and if f ∈W_p+r(u)then

u(t)|Fp(f,t)−F_p,m∗(f,t)| ≤ CkfkW_p+r(u)logm a_m

m r

, 0<C6=C(m,f,t). 3.1 The case t“large”

The method in (11) can be used whentstays between two Gaussian nodes, i.e. for|t|< θa_m, or equivalently,m> ^|t|_θβ

. Thus, the “larger”tis, the “larger”mwill be. For instance, for t=100,β=2 andθ =¹2, it should be chosenm>2500. In such cases the computation of weights and nodes in the Gaussian rule is too much expensive, and sometimes unfeasible. For all these reasons, we propose here an alternative procedure, which is essentially a shrewd application of the Gaussian rule again, in the sense we go to precise. For any fixedt, letmbe such that

|t|>x_m,j+1,

wherejis defined in (6). SettingG_t(x) =_(x−t)^f(x)p+1, we approximate the integral by them−th truncated Gauss-Freud rule, i.e.

Hp(f w_β,t) =X

|i|≤j

G_t(x_m,i)λm,i+R_m(G_t). (15)

Since f andG_t(f)for|t|>x_m,j+1 have the same smoothness, by (10) and (4), iff ∈W_r(u),r∈N, we get

|Rm(G_t)| ≤C

kG^(r)_t uka_m m

r

+e^−AmkGtuk

, (16)

where 0<C6=C(m,G_t)and 0<A6=A(m,Gt).

We observe that this procedure can be successfully applied whenf is “smooth” andtis “large”, say|t|>x_m,j+1, wherejis defined in (6). Moreover, smoother is the functionf and smaller is the value ofmto achieve the desired accuracy. Of course, the error bound in (16) holds for a fixedmand therefore the limit onmofR_m(F_t)has no meaning.

In conclusion, according to the position oftwe use (11) or (15), i.e.

H_p(f,t) =Ψp,m(f,t) +τp,m(f,t), (17) where

Ψp,m(t) =















Fp,m∗(f,t) +

p

X

k=0

f^(k)(t)

k! H_p−k(w_β,t), |t| ≤θa_m,

X

|i|≤j

G_t(x_m,i)λm,i+

p

X

k=0

f^(k)(t)

k! Hp−k(w_β,t), |t|>x_m,j+1, approximatingHp(f w_β,t)for a “wide” range oft.

4 The second method

LetF0(f)be defined in (2). For any fixedθ∈(0, 1), denoting by{x_m+1,k}^b

m+1 2 c

k=−b^m+1₂ cthe zeros ofp_m+1(w_β), let us define the index q:=q(m+1)by

x_m+1,q= min

k=1,2,..,b^m+1₂ c

x_m+1,k:x_m+1,k≥θa_m+1 , and set

F_0,m+1(f,t) =

q

X

i=−q

λm+1,i

f(x_m+1,i)−f(t)

x_m+1,i−t . (18)

By approximatingF₀(f)by the Lagrange polynomialL_m+2(w_β,F₀(f))defined in (7), in view of (1) Fp(f,t) = 1

p!L^(p)_m+2(w_β,F0(f),t) +ηp,m(f,t).

However, since the samplesF0(f,xm,k)are integrals containing the Freud weightw_β, we use formula (18) to evaluate them, i.e.

F_p(f,t) = 1

p!L_m+2^(p) (w_β,F_0,m+1(f),t) +ρp,m(f,t).

Therefore we set

Hp(f w_β,t) = 1

p!L^(p)_m+2(w_β,F_0,m+1(f),t) +

p

X

k=0

f^(k)(t)

k! Hp−k(w_β,t) +ρp,m(f,t)

=: Φp,m(f,t) +ρp,m(f,t), (19)

whereρp,m(f,t)is the remainder term.

We remark that in view of[27, Proposition 1], Gaussian knots and interpolation nodes are sufficiently far among them. This good “distance" prevents numerical cancellation when computingx_m+1,i−x_m,k. About the stability and the convergence of the procedure, we are able to prove the following results

Theorem 4.1. Let p∈N0and0< λ <1. Then, for any function f ∈Z_p+λ(u)

kΦp,m(f)uk ≤CkfkZ_p+λ(u)log²m, 0<C6=C(m,f). (20) Moreover, if f ∈Z_p+r+λ(u),r∈N, then

kρp,m(f)uk ≤CkfkZ_p+r+λ(u)log²m a_m

m r

, 0<C6=C(m,f). (21)

In particular, forf belonging to Sobolev spaces, we get

Corollary 4.2. Let p∈N0and r∈N. Then, for any function f ∈W_p+1(u)

kΦp,m(f)uk ≤CkfkW_p+1(u)log²m, 0<C6=C(m,f), and if f ∈W_p+r+1(u)then

kρp,m(f)uk ≤CkfkW_p+r+1(u)log²m a_m

m r

, 0<C6=C(m,f).

As we have announced, it is possible to implement the method without computing the derivatives{f^(k)(t)}_k=1^p . Indeed, we can approximatef^(k)(t)withL^(k)_m+2(w_β,f,t)fork=0, 1, 2, . . . ,p, whereL_m+2(w_β,f)is the truncated Lagrange polynomial defined in (7). By this way, reusing the same samples involved in the evaluation ofΦp,m(f,t)and taking into account thatH_p−k(w,t)can be computed with the desired accuracy, we get

Hp(f w_β,t) = Φp,m(f,t) +

p

X

k=0

p k

‹

L^(k)_m+2(w_β,f,t)Hp−k(w_β,t) +Tp,m(f,t), (22) where

Tp,m(f,t) = ρp,m(f,t) +

p

X

k=0

p k

‹

τk,m(f,t)H_p−k(w_β,t), with

τk,m(f,t) = (f(t)−L_m+2(w_β,f,t))^(k). The next theorem deals with the pointwise estimate of the errorT_p,m(f,t):

Theorem 4.3. Let p≥0. Then, for any function f ∈Z_p+r+λ(u)with0< λ <1, r≥1, we have

|Tp,m(f,t)u(t)| ≤ C

|t|^pkfkZ_p+r+λ(u)

a_m m

r

log²m, 0<C6=C(m,f). Remark1. For|t| ≥C,Cbeing an arbitrary positive constant, we get

|Tp,m(f,t)u(t)| ≤CkfkZ_p+r+λ(u)

a_m m

r

log²m, 0<C6=C(m,f).

5 Computational aspects

First of all, we point out that sincew_β(x) =e^−|x|β,β >1, is not a classical weight, the coefficients of the three-term recurrence relation of the corresponding orthonormal polynomials are unknown, in the general case. However, it is possible to calculate them, for instance, by the Chebyshev algorithm, which requires the moments

µk= Z

R

x^kw_β(x)d x, k=0, 1, . . . .

Even if the algorithm suffers of ill conditioning, we have successfully implemented it inWolfram Mathematica Language, in extended arithmetic with high accuracy by using the functionsaChebyshevAlgorithmandaGaussianNodesWeightsof the software packageOrthogonalPolynomialsby Cvetkovic and G.V. Milovanovic (see[2]).

Now we give some details about the computation of the Hadamard transforms of the weightw_β(x) =e^−|x|β. In order to use the relation

H_p−k(w_β,t):= 1 (p−k)!

d^p−k d t^p−k

Z

−

R

w_β(x)

x−t d x, (23)

we first deduce the expression ofH₀(w_β,t),k= 0, 1, . . . ,p, for some choices of the parameterβ >1. SinceH₀(w_β,t) =

−H0(w_β,−t), we assumet>0. Then, we have H0(w_β,t) =

Z

−

+∞

0

e^−xβ x−td x−

Z +∞

0

e^−xβ x+td x= 1

β

¨Z

−

+∞

0

e^−xx^1β−1 x^1β −t

d x− Z +∞

0

e^−xx^1β−1 x^1β +t

d x

« . In the caseβ=p/q,p,q∈N,p>q, we have forpeven

H0(w_β,t) =2q p

p

X2 j=1

t^2j−1 Z

−

+∞

0

x^{qp(p−2j+1)−1}

x^q−t^p e^−xd x (24)

and forpodd

H0(w_β,t) =q p

p

X

j=1

t^j−1

¨Z

−

+∞

0

x^{qp(p−j+1)−1}

x^q−t^p e^−xd x+ (−1)^j Z +∞

0

x^{qp(p−j+1)−1} x^q+t^p e^−xd x

«

. (25)

Integrals in (24), (25) can be found in a closed form for some values ofq. In particular forq=1, taking into account[30, p. 325 n. 16]

Z

−

+∞

0

e^−xx^α

x−t d x = −πt^αe^−tcot((1+α)π) +Γ(α)1F₁(1, 1−α,−t), α6=0 Z

−

+∞

0

e^−x

x−td x = −e^−tEi(t), and[10, p.217, n.(17)]

Z+∞

0

e^−xx^α

x+td x=Γ(1+α)t^αe^tΓ(−α,t),

where1F₁(a,b,x)is the Confluent Hypergeometric function, Ei(t)is the exponential integral function, andΓ(a,b)is the incomplete Gamma function, we get forpeven

H0(w_β,t) = 2 p

p

X2 j=1

t^2j−1 Z

−

+∞

0

e^−xx^1−2jp x−t^p d x

= −2π p

p 2−1

X

j=1

e^−tpcot



1+1−2j p

‹ π

‹ +Γ

1−2j p

‹

1F₁



1, 1−1−2j p ,−t^p

‹

−1

pe^−tpEi(t^p) and forpodd

H₀(w_β,t) = 1 p

p

X

j=1

t^j−1

¨Z

−

+∞

0

e^−xx^1−jp

x−t^p d x+ (−1)^j Z +∞

0

e^−xx^1−jp x+t^p d x

«

= −1

pe^−tpEi(t^p) +1 p

p

X

j=2

•

−πe^−tpcot



1+1−j p

‹ π‹

+t^j−1Á1−j p

‹

1F₁



1, 1−1−j p ,−t^p

‹˜

+ 1

p

p

X

j=1

(−1)^jΓ

 1+1−j

p

‹ e^tpΓ



−1−j p ,t^p

‹ .

Thus, taking into account[13, p. 1086, 9.213]

d

d tEi(t) =−d d t

Z+∞

−t

e^−x x d x=e^t

t , d

d t¹F₁(a,b;t) = a

b¹F₁(a+1,b+1,t),

by (23), the functions{Hp−k(w_β,t)}^p_k=0are completely determined. We point out that the computation of the functionEiand1F₁ have been performed by theWolfram MathematicaroutinesExpIntegralEiandHypergeometric1F1, respectively.

In the caseq=2 andq=4, the functionH₀(w_β)can be obtained by combining integrals in[13, 3.354.1 p. 359 and 3.358 p.

361].

Finally, we give some details about the computation of the derivatives of the fundamental Lagrange polynomials{`⁽ⁱ⁾m,k(t)}^j_k=−j, i=1, 2, . . . ,p. Setting

S_r,k(t) = X

|h| ≤_m

2

h6=k

r!

[(t−xm,h)(t²−a²_m)]^r+1, t6=±am, t6={xm,h}b^m2c

h=−b^m2c,

it is no hard to prove that

`⁽ⁱ⁾m,k(t) =

i−1

X

r=0

i−1 r

‹

(−1)^r`⁽ⁱ⁻m,k^1−r)(t)S_r,k(t), for−j≤k≤jandi=1, 2, . . . ,p.

6 Numerical Tests

In this section we propose some numerical experiments obtained by implementing the above introduced methods. To be more precise, we show the numerical results obtained approximatingHp(w_β,f)by the sequence{Ψp,m(f,t)}mof the first method (17) and by{Φp,m(f,t)}mof the second method (19). Since the exact values of the integrals are unknown, we will retain exact the values computed withm=1000, setting

τ¯p,m(f,t) =u(t)|Ψp,m(f,t)−Ψp,1000(f,t)|, ρ¯p,m(f,t) =u(t)|Φp,m(f,t)−Φp,1000(f,t)|, as weighted absolute errors.

The “truncation intervals”, depending on a fixedθ, have been empirically detected. In particular, denoting byepsthe machine precision, in the first method we set

j:=j1−j2

with

j₁:=













 max

k=0,...,b^m∗2c

(

k : λm∗,k

f(x_m∗,k)−Pp k=0

f^(k)(t)

k! (x_m∗,k−t)^p−k (x_m∗,k−t)^p+1

≥eps )

in (11)

max

k=0,...,b^m∗₂c

k : |Gt(x_m,k)|λm,k≥eps in (15)

j₂:=













 min

k=−b^m∗₂c,...,−1

(

k : λm^∗,k

f(x_m∗,k)−Pp k=0

f^(k)(t)

k! (x_m∗,k−t)^p−k (x_m∗,k−t)^p+1

≥eps )

in (11)

min

k=−b^m∗₂c,...,−1

k : |Gt(xm,k)|λm,k≥eps in (15)

and for the second method in (19), setting

¯

q_k:= max

i=0,...,b^m+12 c

i :

λm+1,i

f(x_m+1,i)−f(x_m,k) x_m+1,i−x_m,k

≥eps

,

ˆ

q_k:= min

i=−b^m+1₂ c,...,−1

i :

λm+1,i

f(x_m+1,i)−f(x_m,k) x_m+1,i−xm,k

≥eps

,

¯

q:= max

k=j2,...,j₁q¯_k, ˆq:= min

k=j2,...,j₁ˆq_k, ˆı₁:= max

k=0,...,b^m₂c

( k :

`^(p)m,k(t)

¯q_k

X

i=ˆq_k

λm+1,i

f(x_m+1,i)−f(x_m,k) x_m+1,i−xm,k

≥eps )

,

ˆı2:= min

k=−b^m₂c,...,−1

( k :

`^(p)m,k(t)

¯ q_k

X

i=ˆq_k

λm+1,i

f(x_m+1,i)−f(x_m,k) x_m+1,i−x_m,k

≥eps )

,

we let

ˆı:=ˆı₁−ˆı₂, q:=q¯−ˆq.

We point out that for the simultaneous approximation ofH0(f w_β,t), . . . ,Hp(f w_β,t)atsvalues oft, the first method requires j+ (p+1)ssamples of the density functionf, while the second needsq+ˆı+ (p+1)sevaluations off. In both the cases, the effort doesn’t depend ons.

In the numerical tests proposed below, we compare our results among them and with those achieved by a method discussed in[25, p. 14-15]for finite part integrals on bounded intervals. Thus we have properly adapted this procedure, in the sense we go to specify. Starting from the decomposition

H_p(f w_β,t) = F_p(f,t) +

p

X

k=0

f^(k)(t)

k! H_p−k(w_β,t),

by an idea in[15](see also[25]), to approximateFp(f,t), the interval(−∞,+∞)is broken up into(−∞,t)and(t,+∞)

F_p(f,t) = Z+∞

0

f(t−y)−Pp k=0

f^(k)(t) k! (−y)^k

(−y)^p+1 e^{−|t−y|β+y}e^−yd y+ Z+∞

0

f(y+t)−Pp k=0

f^(k)(t) k! (y)^k

(y)^p+1 e^−|y+t|β+ye^−yd y, and using them-th truncated Gauss-Laguerre rule w.r.t. the weighte^−y[19], we obtain

H_p(f w_β,t) =

k₁

X

i=1

f(t−ym,i)−Pp k=0

f^(k)(t) k! (−ym,i)^k

(−y_m,i)^p+1 e^{−|t−ym,i|β+ym,i}λ¯m,i+

k₂

X

i=1

f(ym,i+t)−Pp k=0

f^(k)(t) k! (ym,i)^k

(y_m,i)^p+1 e^{−|ym,i−t|β+ym,i}λ¯m,i

+

p

X

k=0

f^(k)(t)

k! Hp−k(w_β,t) +σp,m(f,t)

:= J_p,m(f,t) +σp,m(f,t), (26)

wherey_m,iand ¯λm,i,i=1, . . . ,m, are the knots and the weights of them−th Gauss-Laguerre rule and the parametersk₁andk₂ are empirically detected with criteria similar to those used for the rule (11). We will set

σˆp,m(f,t) =|Jp,m(f,t)−J_p,1000(f,t)|.

We point out that by implementing formula (26), for any value ofttwo changes of variables depending ontare required and none of the samples of the functionf can be saved. This means that the computation{Hi(f,t)}^p_i=0atsvalues oftrequires almosts(k₁+k₂+p+1)samples off, which is far greater than the number of samples required in our methods.

All the computations have been performed in double-machine precision (eps≈2.22044e−16) and in all the tables the symbol “-” means that the machine accuracy has been achieved.

Example6.1.

Hp(f w₂,t) = Z

R

sin ^x₂

cos(x−e)

(x−t)^p+1 e^−x2d x, p=0, 1, 2.

Since the function f(x) =sin ^x₂

cos(x−e)is very smooth we expect for a fast convergence. As the test shows, the sequence {Jp,m(f,t)}mappears slower convergent w.r.t. both the sequences{Ψp,m(f,t)}mand{Φp,m(f,t)}m. Indeed, looking at Tables 1-2, H0(f;t),H1(f;t),H2(f;t)are computed for different values oft, attaining the machine precision withm=30 in the first method and withm=40 in the second one. On the contrary, as Table 3 evidences, by the method (26) the machine precision is not always achieved form=160.

In Figure 2 we show the graphs of the approximating functionsΦp,40(f,t),p=0, 1, 2.

m j τ¯0,m(f,−3) τ¯0,m(f,−0.5) τ¯0,m(f, 4) τ¯0,m(f, 10)

10 8 5.11e−7 3.67e−6 1.07e−12 −

20 18 8.50e−16 1.55e−14 − −

30 27 − − − −

m j τ¯1,m(f,−3) τ¯1,m(f,−0.5) τ¯1,m(f, 4) τ¯1,m(f, 10)

10 8 9.71e−6 1.01e−6 5.87e−9 −

20 18 8.47e−16 7.90e−15 − −

30 27 − − − −

m j τ¯2,m(f,−3) τ¯2,m(f,−0.5) τ¯2,m(f, 4) τ¯2,m(f, 10)

10 8 5.85e−5 8.35e−7 1.51e−6 −

20 18 − 4.77e−15 − −

30 27 − − − −

Table 1:Example6.1: ¯τp,m(f,t),p=0, 1, 2, witht=−3,−0.5, 4, 10