JID:APNUM AID:4602 /FLA [m3G; v1.338] P .1 (1-9)
Applied Numerical Mathematics ••• ( •••• ) ••• – •••
Contents lists available at ScienceDirect
Applied Nu m er i ca l Mathematics
journal homepage: www .elsevier .com/locate/apnum
Computing int egr als with an exponential wei gh t on the r eal
axis in floating point arithmetic
Te r e s a L a u d a d i o a , ∗ , Nicola Mastr onardi a , ∗ , Donat ella Occorsio b , c
a Istituto per le Applicazioni del Calcolo “M. Picone”, Consiglio Nazionale delle Ricerche, Bari, Italy
b Department of Mathematics, Computer Science and Economics, University of Basilicata, Potenza, Italy
c Istituto per le Applicazioni del Calcolo “M. Picone”, Consiglio Nazionale delle Ricerche, Nap ol i, Italy
a r t i c l e i n f o a b s t r a c t
Article history :
Rec e ive d 12 December 2022
Rec e ive d in revi se d form 28 May 2023
Accep ted 30 May 2023
Ava i l a b l e online xxxx
K eywor ds:
Gaussian quad ratu re rules
Golub and We l s c h alg orithm
Singular Value Decomposition
Product integration rules
The aim of this work is to propose a fast and r eliable alg orithm for computing integr als of
the type
∞
−∞
f ( x ) e − x 2 − 1
x 2 dx ,
where f ( x ) is a sufficientl y smoo th function, in floating point arithmetic. The alg orithm
is based on a product integr ation rule, whose rate of conver gence depends onl y on the
reg ul ar i ty of f , since the coefficients of the rule are “exactl y” comput ed by means of
suitable r ecurrence rel at i on s her e derived. We pro v e st ability and conv ergence in the space
of locally continuous functions on R eq uipped with wei g hte d uniform norm. By ex te ns ive
numerical tests , the accuracy of the pr oposed product rule is compar ed with that of the
Gauss–Hermite quad ra tu re formula w. r. t . the function f ( x ) e − 1
x 2 . The numerical res ul ts
confirm the effecti v eness of the method, supporting the pr o ven theor etical es timates.
© 2023 IMA CS. Published by Elsevier B.V. All rights res er ve d.
1. Introduction
The present paper deals with the computation of int egr als
∞
−∞
f ( x ) ˜
w ( x ) dx , (1)
where ˜
w ( x ) = e − x 2 − 1
x 2 is a combination of the Pollaczek–type wei g ht e − 1
x 2 and the Hermite wei gh t e − x 2 , and f is a
“smooth” function.
A str aightforwar d approac h to compute an appr o ximation of ( 1 )i s to consider the Gauss qua d ra t u r e rule ( GQR ) with
wei g h t function w ( x ) = e − x 2 and integr and function e − 1
x 2 f ( x ) . U nfortunatel y , the presence of the f actor e − 1
x 2 as part of the
* Corresponding auth ors .
E-mail addresses: teresa.laudadio@cnr .it (T . Laudadio), nicola.mastr onardi@cnr .it (N. Mastronar di), [email protected] (D. Occorsio).
https://doi.org/10. 1016/j.apnum.2023.05.025
0168-9274/ © 2023 IMACS. Published by Elsevier B.V. All rights re se rve d.
JID:APNUM AID:4602 /FLA [m3G; v1.338] P .2 (1-9)
T . Laudadio, N. Mastronardi and D. Occorsio Applied Numerical Mathematics ••• ( •••• ) ••• – •••
integr and function stron gly re du ce s the conver gence rate of GQR . On the ot he r hand, the computation of integr als ( 1 )i s
of inter est in numerical me thods for integral eq uations, which in turn are model for man y applications arising in applied
sciences. As an exa m pl e , integr al eq uations with the same wei g ht function ˜
w , but ext e nd e d to the positi ve real semiaxis,
model some mathematical financial problems [ 7 ]. Appr o ximation pr operties rela te d to Pollaczek–Laguerre wei gh t s of the
type W ( x ) = x γ e − x β − 1
x α on the positiv e semiaxis, ha ve been exte n s ive ly stu die d in [ 11 , 12 , 15 , 16 ]. In particular in [ 5 , 6 , 11 ]
the auth or s proposed a Gaussian qua d ra t u r e rule for appro ximating int egrals on the positi ve semiaxis with a Pollaczek–
Laguerre wei g ht and only with the use of the ext e nd e d precision they were able to over come the numerical instability in
computing the corresponding zer os and wei g ht s [ 6 , p. 230]. Such inst ability occurs only in com puting the zeros and the
coefficients of the Gaussian rule for a non standar d wei g h t, in applying the Cheby shev me thod of moments, and does not
affect the s tability of the Ga ussian rule and the rela ted methods considered ther e. Her e, we deal with the Pollaczek–Hermit e
wei g h t ˜
w ( x ) = e − x 2 − 1
x 2 proposing a fast and reliable algorithm for the efficient computation of int egr als on R of type ( 1 )
without using ex te nd e d precision. The product int egration rule we propose is based on the appro ximation of f by means of
a “truncated” Lagrang e pol ynomial interpolating f , and hence on the exa c t computation of the rule coefficients, by means
of suitable recurrence relat io ns here deri v ed. In view of the truncation, the rule requi res a red uc e d number of function
evaluatio ns. Moreov er , fo r functions f belo ngi ng to suitable spaces of locall y continuous functions, endow ed with wei g h te d
uniform norm, we prov e the rule is st able and conv ergent. Since we will pro ve that the qu a d ra t u re error has the same
beha vio r of the bes t polynomial appro ximation of f , the more regular f is, the faster the conv erg ence of the product rule
is.
The paper is organized as follow s. In Section 2 the main featu res of the Hermite polynomials ar e described. Section 3 is
devo ted to function spaces and to Lagrange interpolation processes associat ed to the considered product rule. In Section 4
the construction of the pr oposed product rule is described.
A number of numerical ex am p l es , confirming the stability properties of the pr oposed method, is pr ovided in Section 5 ,
follow e d by the proof s of the theorems stated in Section 6 . The paper ends with the concluding rem ark s .
2. Computing the zeros of orthog onal polynomials
Let H ( x ) = k x + − 1
j = 0 c j x j , = 0 , 1 , ... , be the seq uence of orthonormal Hermite pol ynomials with resp ec t to the
wei g h t function w ( x ) = e − x 2 in the interv al ( −∞ , ∞ ) , i.e.,
∞
−∞
H i ( x ) H j ( x ) w ( x ) dx = δ ij , with δ ij = 1 , if i = j ,
0 , if i = j ,
where k , c j ∈ R , j = 0 , 1 , ..., − 1, and k > 0. P olynomials H ( x ) satisfy the following three–t erm r ecurrence rela ti on [ 20 ]
⎧
⎪
⎨
⎪
⎩
H − 1 ( x ) = 0 ,
H 0 ( x ) = k 0 = 4
1
π ,
β + 1 H + 1 ( x ) = x H ( x ) − β H − 1 ( x ), ≥ 0 ,
(2)
where
β 0 = 0 ,β
= k
k − 1 =
2 , = 1 , 2 ,.... (3)
Using ( 2 ), we can write the n -step r ecurrence rela ti on
J h ( x ) = x h ( x ) − β n H n ( x ) e n ,
where
J =
⎡
⎢
⎢
⎢
⎢
⎢
⎢
⎣
0 β 1
β 1 0 β 2
β 2
. . . . . .
. . . 0 β n − 1
β n − 1 0
⎤
⎥
⎥
⎥
⎥
⎥
⎥
⎦
, h ( x ) = ⎡
⎢
⎢
⎢
⎢
⎢
⎣
H 0 ( x )
H 1 ( x )
.
.
.
H n − 2 ( x )
H n − 1 ( x )
⎤
⎥
⎥
⎥
⎥
⎥
⎦
. (4)
The matrix J is called Jacobi matrix [ 2 ]. The following theorem was pro ved in [ 2 , Th. 1 .31].
2
JID:APNUM AID:4602 /FLA [m3G; v1.338] P .3 (1-9)
T . Laudadio, N. Mastronardi and D. Occorsio Applied Numerical Mathematics ••• ( •••• ) ••• – •••
Theorem 2. 1 . Let J = QX Q T be the spectr al decomposition of J, where X ∈ R n × n is a diagonal matrix and
x := diag ( X ) = [ x 1 ,..., x n ] T , x 1 < x 2 < ··· < x n ,
and Q ∈ R n × n is an orthogonal matrix, i.e., Q T Q = I n . Then H n ( x i ) = 0 , i = 1 , ..., n, and Q = VD , with
V = [ h ( x 1 ), h ( x 2 ) ,..., h ( x n ) ] , D = ⎡
⎢
⎢
⎢
⎣
ˆ
w 1
ˆ
w 2
. . .
ˆ
w n
⎤
⎥
⎥
⎥
⎦ , (5)
and
ˆ
w i = 1
h ( x i ) 2 = 1
n − 1
= 0 H 2
( x i )
.
The eigen v alues x i , i = 1 , 2 , ..., n , ar e the nodes of GQR and the corresponding wei g ht s w i are defined as (see [ 20 ])
w i := w i ( x i ) = ˆ
w 2
i = 1
n − 1
= 0 H 2
( x i ) , i = 1 ,..., n . (6)
As shown in [ 21 ], the wei g ht s can be also obtaine d by the first row of Q as
w i = μ 0 q 2
1 , i , i = 1 ,..., n , (7)
where μ 0 = R e − x 2 dx . The Golub and We l s c h algorithm [ 3 ], rely in g on a modification of the QR –method [ 1 ], yields the
nodes and the wei g h ts of GQR by computing only the eigen v alues of the Jacobi matrix and the first ro w of the associated
eigen v ector matrix.
A different way to compute the wei g ht s of GQR in a stable manner has been described in [ 8 ].
3. Fun c ti o n spaces and Lagrang e interpolation pr ocesses
Fro m now on, C will deno te an y positiv e constant ha ving dif fer ent meanings at different occurr ences and the writing
C = C ( a , b , ..) means that C does no t depend on a , b , ...
For a fix ed δ ∈ R + let w δ ( x ) = √ w ( x )( 1 +| x | ) δ . Let us introduce the space C w δ defined as
C w δ = f ∈ C 0 ( R ) : lim
x →±∞ f ( x ) w δ ( x ) = 0 ,
eq uipped with the norm
f w δ = fw
δ ∞ = sup
x ∈ R | f ( x ) w δ ( x ) | .
Denoting by
E n ( f ) w δ := inf
P n ∈ P n ( f − P n ) w δ ∞ , (8)
the error of bes t pol ynomial appro ximation of f ∈ C w δ , it can be pr o ved that
lim
n →∞ E n ( f ) w δ = 0 ⇐⇒ f ∈ C w δ ,
[ 9 ] (see, e.g. [ 17 , p. 275–276]). Note that functions in C w δ can gro w e xponentially as | x | →∞ . For smoother functions, let
us consider the following Sobolev subspaces of C w δ
W r ( w δ ) = f ∈ C w δ : f ( r − 1 ) ∈ AC ( R ), f ( r ) w δ ∞ < ∞ , r ∈ N , r ≥ 1 ,
where AC ( R ) denot es the set of all functions whic h are absolut ely continuous on ev er y closed subset of R . Such spaces are
eq uipped with the norm
f W r ( w δ ) = fw
δ ∞ + f ( r ) w δ ∞ .
For f ∈ W r ( w δ ) , we re ca ll that the following es timate holds [ 18 , 19 , 13 ]
3
JID:APNUM AID:4602 /FLA [m3G; v1.338] P .4 (1-9)
T . Laudadio, N. Mastronardi and D. Occorsio Applied Numerical Mathematics ••• ( •••• ) ••• – •••
E n ( f ) w δ ≤ C f W r ( w δ )
( √ n ) r , C = C ( n , f ). (9)
Let { x k } n
k = 1 be the zeros of H n and, for a fixed 0 <θ < 1, le t j = j ( n ) be the index of the zer o of H n such that
x j = min
k = 1 , 2 ,..., n { x k : x k >θ √ 2 n } . (10)
Moreo v er , let P n − 1 the subspace of P n − 1 defined as
P n − 1 := { P ∈ P n − 1 : P ( x i ) = 0 , k > j , and k < n + 1 − j } .
Setting
P ∈ P n − 1 the polynomial defined by
( f −
P ) w δ ∞ = inf
P ∈ P n − 1 ( f − P ) w δ ∞ =:
E n − 1 ( f ) w δ ,
in [ 13 , Lemma 3. 1] (in a more gene ra l conte xt) the following es timate was pro ved
E n − 1 ( f ) w δ ≤ C E N ( f ) w δ + e − B n fw
δ ∞ , C = C ( n , f ), (11)
with N = θ
1 + θ 2 n
2 ∼ n , where ξ deno tes the larg est int eger no t exce e d in g ξ ∈ R + , and B ∈ R + , B = B ( n , f ) .
Let L 1 ( R ) be the space of all measur able functions f , eq uippe d with the norm
f L 1 ( R ) =
R
| f ( x ) | dx ,
and let L ∞ ( R ) =: C 0 .
3. 1 . T runcated Lagr ange interpolating polynomial
For a gi ven function f , with j = j ( n ) , defined in ( 10 ), the pol ynomial
L n − 1 ( w , f , x ) =
j
k = n − j + 1
n , k ( x ) f ( x k ), n , k ( x ) = H n ( x )
H
n ( x k )( x − x k ) , (12)
is the truncated Lagrange polynomial interpolating f [ 13 , 14 ]. L n − 1 ( w , f ) ∈ P n − 1 , int erpolates f at the zer os { x k } j
k = n + 1 − j ∈
( − θ √ 2 n , θ √ 2 n ) and v anishes at {{ x k } { k > j }∪{ k < n + 1 − j } } . Differentl y fr om the ordinary Lagrang e polynomial int erpolating f on
all the zeros of H n ( x ) , L n − 1 ( w , f ) is no t a project or in P n − 1 , but it pr eserves pol ynomials bel on gin g to the subspace P n − 1 ,
i.e., L n − 1 ( w , Q ) = Q , Q ∈ P n − 1 .
4. The Product rule
In this section we construct the product rule of interpolation type for computing an appr o ximation of the integral ( 1 ).
Hence, rep la ci ng f by L n − 1 ( w , f ) , we obtain the product rule
∞
−∞
˜
w ( x ) L n − 1 ( w , f , x ) dx =
j
k = n + 1 − j
C k f ( x k ) (13)
:=
j
k = n + 1 − j
f ( x k ) ⎡
⎣ w k
n − 1
= 0
H ( x k )
∞
−∞
˜
w ( x ) H ( x ) dx ⎤
⎦
=
j
k = n + 1 − j
f ( x k ) w k
n − 1
= 0
H ( x k ) M ,
where, for = 0 , 1 , ... ,
M =
∞
−∞
H ( x ) ˜
w ( x ) dx , (14)
are the so called modified moments .
4
JID:APNUM AID:4602 /FLA [m3G; v1.338] P .5 (1-9)
T . Laudadio, N. Mastronardi and D. Occorsio Applied Numerical Mathematics ••• ( •••• ) ••• – •••
4. 1 . Computing the Modified moments M
In this subsection we construct the recurrence rela ti on s for computing the modified moments ( 14 ). Defining
N :=
∞
−∞
H ( x )
x 2 ˜
w ( x ) dx , = 0 , 1 ,...,
the following lemma holds.
Lemma 4. 1 . The sequences M and N , = 0 , 1 , ... satisfy the follo wing recurr ence relations
⎧
⎪
⎨
⎪
⎩
M 0 = 4
√ π
e 2 ,
M = 2 −
√ ( − 1 ) M − 2 + 2
√ ( − 1 ) N − 2 , = 2 , 4 , 6 ,...,
M = 0 , = 1 , 3 , 5 ,...,
(15)
⎧
⎪
⎪
⎪
⎪
⎪
⎨
⎪
⎪
⎪
⎪
⎪
⎩
N 0 = 4
√ π
e 2 ,
N + 2 = 2
√ ( + 1 )( + 2 ) M −
+ 1
+ 2 + + 1
+ 2 N
− − 1
+ 1
+ 2 N − 2 , = 2 , 4 , 6 ,...,
N = 0 , = 1 , 3 , 5 ,....
(16)
4.2. St ability and error estimates
By means of the product rule ( 13 ), it is intr oduced the error
e n ( f ) :=
R
[ f ( x ) − L n − 1 ( w , f , x ) ] ˜
w ( x ) dx , (17)
i.e.
R
f ( x ) ˜
w ( x ) dx =
j
k = n + 1 − j
C k f ( x k ) + e n ( f ), (18)
bein g the rule exa c t fo r P n + 1 ∈ P n + 1 .
Now we state the stability and the con vergence of the formula in C w δ , providing err or estimates in the spaces W r ( w δ ) .
Theorem 4.2. Under the assumption δ> 1
2 , the rule ( 13 ) is stable in C w δ , i.e.,
sup
n
j
i = 1
| C i |
w δ ( x i ) < ∞ , (19)
and
| e n ( f ) |≤ C E N ( f ) w δ + e − B n fw
δ ∞ , (20)
where N = θ
1 + θ 2
n ∼ n, and 0 < B = B ( n , f ), C = C ( n , f ) .
By combining ( 20 ) and ( 9 ), the following err or estimat e follow s.
Corollary 4.3. For any f ∈ W r ( w δ ) , with δ> 1
2 , the following err or estimate holds
| e n ( f ) |≤ C f W r ( w δ )
( √ n ) r , C = C ( n , f ). (21)
5
JID:APNUM AID:4602 /FLA [m3G; v1.338] P .6 (1-9)
T . Laudadio, N. Mastronardi and D. Occorsio Applied Numerical Mathematics ••• ( •••• ) ••• – •••
Ta b l e 1
Values of the appro ximations of the integral ( 22 )o b t a i n e d TPR and GQR , with f ( x ) = cos ( x ) , and Mval =
8 . 945397612471845 × 10 − 2 .
n 2 j − n TPR | TPR − Mval
Mval | GQR | GQR − Mval
Mval |
88 8 . 945098794037276 × 10 − 2 3 . 34 × 10 − 5 6 . 949412400327926 × 10 − 2 2 . 23 × 10 − 1
16 16 8 . 945397611011636 × 10 − 2 1 . 63 × 10 − 10 9 . 723435717459904 × 10 − 2 8 . 70 × 10 − 2
32 32 8 . 945397612471852 × 10 − 2 6 . 21 × 10 − 16 8 . 918619331793731 × 10 − 2 2 . 99 × 10 − 3
64 48 8 . 945397612471877 × 10 − 2 3 . 41 × 10 − 15 8 . 936154227660445 × 10 − 2 1 . 03 × 10 − 3
128 68 8 . 945397612471848 × 10 − 2 1 . 55 × 10 − 16 8 . 945089722880344 × 10 − 2 3 . 44 × 10 − 5
256 90 8 . 945397612471870 × 10 − 2 2 . 64 × 10 − 15 8 . 945673065110885 × 10 − 2 3 . 08 × 10 − 5
512 124 8 . 945397612471838 × 10 − 2 9 . 31 × 10 − 16 8 . 945401734394182 × 10 − 2 4 . 61 × 10 − 7
1024 166 8 . 945397612471874 × 10 − 2 3 . 10 × 10 − 15 8 . 945397746086610 × 10 − 2 1 . 49 × 10 − 8
2048 242 8 . 945397612471845 × 10 − 2 0 . 00 × 10 0 8 . 945397612351721 × 10 − 2 1 . 34 × 10 − 11
4096 404 8 . 945397612471892 × 10 − 2 6 . 05 × 10 − 15 8 . 945397612473283 × 10 − 2 1 . 61 × 10 − 13
Ta b l e 2
Values of the appro ximations of the integral ( 22 )o b t a i n e d appl ying TPR and GQR , with f ( x ) = arctan ( 1 + x
4 ) , and Mval =
5 . 427697244322335 × 10 − 2 .
n 2 j − n TPR | TPR − Mval
Mval | GQR | GQR − Mval
Mval |
88 5 . 427686729587454 × 10 − 2 1 . 94 × 10 − 6 5 . 051927492275927 × 10 − 2 6 . 92 × 10 − 2
16 16 5 . 427697240051552 × 10 − 2 7 . 87 × 10 − 10 5 . 644085331675656 × 10 − 2 3 . 99 × 10 − 2
32 32 5 . 427697244322757 × 10 − 2 7 . 77 × 10 − 14 5 . 414250683007876 × 10 − 2 2 . 48 × 10 − 3
64 48 5 . 427697244322347 × 10 − 2 2 . 17 × 10 − 15 5 . 426208992962754 × 10 − 2 2 . 74 × 10 − 4
128 68 5 . 427697244322358 × 10 − 2 4 . 22 × 10 − 15 5 . 427544256103390 × 10 − 2 2 . 82 × 10 − 5
256 90 5 . 427697244322337 × 10 − 2 2 . 56 × 10 − 16 5 . 427767485244951 × 10 − 2 1 . 29 × 10 − 5
512 124 5 . 427697244322327 × 10 − 2 1 . 53 × 10 − 15 5 . 427698372784009 × 10 − 2 2 . 08 × 10 − 7
1024 166 5 . 427697244322347 × 10 − 2 2 . 17 × 10 − 15 5 . 427697278329974 × 10 − 2 6 . 27 × 10 − 9
2048 242 5 . 427697244322321 × 10 − 2 2 . 68 × 10 − 15 5 . 427697244286701 × 10 − 2 6 . 57 × 10 − 12
4096 407 5 . 427697244322308 × 10 − 2 4 . 99 × 10 − 15 5 . 427697244322690 × 10 − 2 6 . 53 × 10 − 14
5. Num e r i c a l Examples
In this section we st ate the numerical re su lts obtained by appro ximating integrals
I ( f ) =
∞
−∞
e − x 2 − 1
x 2 f ( x ) dx , (22)
for functions f belon gin g to different spaces of functions. We compare the resu lts achiev ed by the proposed product rule
( TPR ) with those obtained by the truncated Gauss–Hermit e qua d r a tu r e rule ( GQR ) [ 10 ], considering e − x 2 as the Hermite
wei g h t and g ( x ) = e − 1
x 2 f ( x ) as function. In the tables n denotes the number of points in volv ed in TPR and GQR , re sp ec t ively .
For each test, the appro ximations of ( 22 ) obtaine d by the two methods, TPR and GQR , for n = 2 i , i = 3 , 4 , ..., 12, are
computed in Matlab R2022a with machine pr ecision ∼ 2 . 22 × 10 − 16 and compared with the val u e of the integral
computed by the function NIntegrate of Mathematica , with a wo rk in g precision of 500 digits (denot ed by Mval ),
considered as the exa c t val u e of the integral. The com puted appro ximations and the corresponding rel at ive errors are
rep or te d in the associated tables.
It can be no ticed that the ra ti o ( 2 j − n )/ n b ecomes smaller and smaller as n incr eases. Therefor e, the pr oposed truncated
product rule has a computational cos t much smaller than the non truncat ed one for large n .
Example 5. 1 . In this ex am p le , f ( x ) = cos ( x ) in ( 22 ) and Mval = 8 . 945397612471845 × 10 − 2 . The re su lts are displa y ed in
T able 1 . They sho w that the rela tive errors of TPR ar e close to the machine precision with n ≈ 32 , since cos ( x ) is a ver y
regular function. About the slo w con verg ence of GQR , we observe that, even though g ( x ) := e − 1
x 2 cos ( x ) ∈ W r ( w δ ) for an y
r ≥ 1, the seminorms of g appearing in the err or estimates becomes larger and larger as r increases, by slowing do wn the
speed of conv erg ence of GQR . For instance, g W r ( w δ ) ∼ 5 × 10 3 fo r r = 5, and g W r ( w δ ) ∼ 6 × 10 11 for r = 10 .
Example 5.2. Her e, f ( x ) = ar ctan 1 + x
4 and Mval = 5 . 427697244322335 × 10 − 2 . The res ul ts are displa y ed in T able 2 . Simi-
larl y to Example 5. 1 , the re la tive errors of TPR are close to the machine precision with n ≈ 32 , bein g arctan 1 + x
4 a regular
function, while the conv erg ence of GQR is almost linear since it is applied to h ( x ) := e − 1
x 2 arctan 1 + x
4 ∈ W r ( w δ ) for an y
r ≥ 1, and whose seminorms become larger and larger as r incr eases, stron gly affecting the speed of converg ence. For
instance, h W r ( w δ ) ∼ 2 × 10 11 , for r = 10.
6
JID:APNUM AID:4602 /FLA [m3G; v1.338] P .7 (1-9)
T . Laudadio, N. Mastronardi and D. Occorsio Applied Numerical Mathematics ••• ( •••• ) ••• – •••
Ta b l e 3
Values of the appro ximations of the integral ( 22 )o b t a i n e d appl ying TPR and GQR , with f ( x ) =| cos ( x ) | 5 / 4 , and Mval =
9 . 105558869283804 × 10 − 2 .
n 2 j − n TPR | TPR − Mval
Mval | GQR | GQR − Mval
Mval |
88 9 . 485235885806660 × 10 − 2 4 . 17 × 10 − 2 7 . 301487850488275 × 10 − 2 1 . 98 × 10 − 1
16 16 9 . 074650084283779 × 10 − 2 3 . 39 × 10 − 3 9 . 863496922603453 × 10 − 2 8 . 32 × 10 − 2
32 32 9 . 182616784940598 × 10 − 2 8 . 46 × 10 − 3 9 . 162629410875198 × 10 − 2 6 . 27 × 10 − 3
64 48 9 . 079166478725698 × 10 − 2 2 . 90 × 10 − 3 9 . 070205616732352 × 10 − 2 3 . 88 × 10 − 3
128 68 9 . 121947847053928 × 10 − 2 1 . 80 × 10 − 3 9 . 121708479806487 × 10 − 2 1 . 77 × 10 − 3
256 90 9 . 104913976690782 × 10 − 2 7 . 08 × 10 − 5 9 . 105186176595942 × 10 − 2 4 . 09 × 10 − 5
512 124 9 . 109013869493021 × 10 − 2 3 . 79 × 10 − 4 9 . 109017834749777 × 10 − 2 3 . 80 × 10 − 4
1024 166 9 . 104375386675621 × 10 − 2 1 . 30 × 10 − 4 9 . 104375504983463 × 10 − 2 1 . 30 × 10 − 4
2048 240 9 . 106285825658801 × 10 − 2 7 . 98 × 10 − 5 9 . 106285824155236 × 10 − 2 7 . 98 × 10 − 5
4096 404 9 . 105641735268770 × 10 − 2 9 . 10 × 10 − 6 9 . 105630563127487 × 10 − 2 7 . 87 × 10 − 6
Example 5.3. In this test f ( x ) =| cos ( x ) | 5 / 4 and Mval = 9 . 105558869283804 × 10 − 02 .
The res ul ts are display e d in T able 3 . The function | cos ( x ) | 5 / 4 ∈ W 1 ( w δ ) and the spee d of con ver gence of bo th the rules
TPR and GQR is ver y slow .
6. Proof s
In this section, the proofs of Lemma 4. 1 and Theorem 4.2 ar e outlined.
Proof of Lemma 4. 1 . Ta k i n g into account the three–t erm r ecurrence relation ( 2 ), we have
β M =
∞
−∞
x H − 1 ( x ) e − x 2 − 1
x 2 dx − β − 1 M − 2 , ≥ 2 . (23)
Let us no w comput e the integr al in ( 23 ). Since [ 20 , p. 106],
e − x 2 H − 1 ( x ) =− 1
√ 2 ( − 1 )
d
dx e − x 2 H − 2 ( x ) , (24)
integr ating by parts ( 23 )w e obtain:
∞
−∞
x H − 1 ( x ) ˜
w ( x ) dx = 1
√ 2 ( − 1 )
∞
−∞
H − 2 ( x ) ˜
w ( x ) dx
+ 2
− 1
∞
−∞
H − 2 ( x )
x 2 ˜
w ( x ) dx
= 1
√ 2 ( − 1 ) M − 2 + 2
− 1 N − 2 .
Therefor e, the recurr ence rela ti on ( 15 )f o r M holds, where the coefficients β are gi ven in ( 3 ).
To compute the seq uence N , = 0 , 1 , 2 , ... , we need the following lemma.
Lemma 6. 1 . The sequence of even (odd) Hermite polynomials H k ( x ), k = 0 , 2 , 4 , ... ( k = 1 , 3 , 5 , ... ) satisfies the follo wing three–
recurr ence rel at io n
β k + 1 β k + 2 H k + 2 ( x ) = x 2 − (β 2
k + β 2
k + 1 ) H k ( x ) − β k − 1 β k H k − 2 ( x ),
Proof. Le t us consider ( 2 )e v a l u a t e d for k = k + 1 and k = k − 1, resp ec tive ly ,
H k + 1 ( x ) = β k + 2 H k + 2 ( x ) + β k + 1 H k ( x )
x (25)
and
H k − 1 ( x ) = β k H k ( x ) + β k − 1 H k − 2 ( x )
x . (26)
7
JID:APNUM AID:4602 /FLA [m3G; v1.338] P .8 (1-9)
T . Laudadio, N. Mastronardi and D. Occorsio Applied Numerical Mathematics ••• ( •••• ) ••• – •••
Rep l a ci n g ( 25 ) and ( 26 )i n t o ( 2 )f o r = k , we then obtain:
β k + 1
β k + 2 H k + 2 ( x ) + β k + 1 H k ( x )
x + β k
β k H k ( x ) + β k − 1 H k − 2 ( x )
x = x H k ( x ),
and thus
β k + 1 β k + 2 H k + 2 ( x ) = x 2 − (β 2
k + β 2
k + 1 ) H k ( x ) − β k − 1 β k H k − 2 ( x ), (27)
i.e., we obtain a three–term recurr ence rela ti on for the odd and even polynomials ext ra ct e d from the seq uence H ( x ) ,
= 0 , 1 , 2 , ... .
Therefor e, the recurr ence rela ti on ( 16 )i s obtained dividing ( 27 )b y x 2 and considering the wei g hte d integral, with wei gh t
˜
w , on both sides.
The va lu e s of the modified moments M 0 and N 0 can be retr ie ve d from [ 4 , p. 337] and [ 4 , p. 369], res pe ct ively .
Observe that M = 0 , ev en, since ˜
w is an even function and H is an odd function for odd.
Proof of Theor em 4.2 . First of all we pro ve that for any f ∈ C w δ ,
R
| L n − 1 ( w , f , x ) | ˜
w ( x ) dx ≤ C fw
δ ∞ , C = C ( n , f ). (28)
Since the fact or e − 1
x 2 ∈ C ∞ ( R ) satisfies 0 ≤ e − 1
x 2 < 1, we hav e
R
| L n − 1 ( w , f , x ) | e − x 2 − 1
x 2 dx ≤ ⎛
⎝
R
w ( x ) dx ⎞
⎠
1
2 ⎛
⎝
R
( L n − 1 ( w , f , x )) 2 w ( x ) dx ⎞
⎠
1
2
.
By using the n -th Gauss-Hermite rule whic h is exa c t fo r the polynomial ( L n − 1 ( w , f )) 2 ∈ P 2 n − 2 , we get
⎛
⎝
R
( L n − 1 ( w , f , x )) 2 w ( x ) dx ⎞
⎠
1
2
= ! n
k = 1
w k [ f ( x k ) ] 2 " 1
2
×
! n
k = 1
w k
( 1 +| x k | ) 2 δ [ f ( x k ) w δ ( x k ) ] 2 " 1
2
≤ C fw
δ ∞ ! n
k = 1
w k
( 1 +| x k | ) 2 δ w ( x k ) " 1
2
≤ C fw
δ ∞ ⎛
⎝
R
dx
( 1 +| x | ) 2 δ ⎞
⎠
1
2
≤ C fw
δ ∞ ,
taking into account the assump tion δ> 1
2 .
We omit the proof of ( 19 )s i n c e it can be easily deduced by standard arguments by ( 28 ). In order to prov e ( 20 ), let
P ∈ P n − 1 , as defined in ( 11 ).
Then
| e n ( f ) | ≤
R # # ( f ( x ) −
P ( x )) # # ˜
w ( x ) dx
+
R # # L n + 1 ( w , f −
P , x ) # # ˜
w ( x ) dx
≤ ( f −
P ) w δ ∞
R
w ρ ( x )
w δ ( x ) dx + C ( f −
P ) w δ ∞
≤ C ˜
E n − 1 ( f ) w δ .
The thesis follow s by ( 11 ).
8
JID:APNUM AID:4602 /FLA [m3G; v1.338] P .9 (1-9)
T . Laudadio, N. Mastronardi and D. Occorsio Applied Numerical Mathematics ••• ( •••• ) ••• – •••
7. Conclusions
The construction of a fast and reliable product rule for computing int egrals in volving functions of kind e − x 2 − 1
x 2 f ( x ) on
the whole rea l line in floating point arithmetic is considered.
It is shown that the con vergence pr operties of the proposed method onl y rely on the reg ul ar it y properties of f ( x ) rat he r
than those of f ( x ) e − 1
x 2 .
The numerical experiment s confirm the effecti veness of the pr oposed approac h.
Declaration of compe ting inter est
The auth or s declare that they ha ve no conflict of interes t.
Data av ailability
Data sharing not applicable to this article as no datasets wer e ge ne rate d or analyzed during the current stud y .
Ac kno w ledgements
This wor k wa s partially supported by Gruppo Nazionale Calcolo Scientifico (GN CS) of Istitut o N azionale di Alta Matem-
atica (INdD AM). The au tho rs thank the anon ymous Re f er ee for his/her comments and sugges tions which helped to impro v e
the paper . Third autho r is partially supported by INdAM GNCS project 2022 “Me t odi e softwar e per la modellistica int egrale
multiv ariata”.
Re fe re n c es
[1] H. Bowdler , R.S. Martin, C. Rei n sch , J.H. Wilkinson, The QR and QL algorithms for symmetric matrices, Num er. Math. 11 (1968) 293–306.
[2] W. Gautschi, Orthogonal Po lynom ial s: Computation and Appro ximation, Oxford Un ive r si ty Press, Oxford, UK, 2004.
[3] G.H. Golub, J.H. We l s c h , Calculation of Gauss qu adra tu re rules, Math. Comput. 23 (106) (1969) 221–230.
[4] I.S. Gradshteyn, I.M. Ryzhik , Ta b l e of Integr als, Series, and Pr oducts, seventh edition, Acad em ic Press, Boston, 2007 .
[5] P. Junghanns, G. Mastr oianni, I. Notara ngel o, On Nyström and pr oduct integr ation methods for Fre d ho l m integral equations, in: J. Dic k, F. Y. Kuo, H.
Wo ´ zniakow ski (Eds.), Cont empor ary Computational Mathematics -a Celebration of the 80th Birthday of Ian Sloan, Spring er Int ernational Publishing,
2018, pp. 645–673.
[6] P. Junghanns, G. Mastroianni, I. Notarange lo , Wei g h t e d Polynom ial Appro ximation and Nu me r ic a l Methods for Integral Equations, Pa thways in Mathe-
matics, Birkäuser–Springer , Cham, 2021 .
[7] M. Ki jima, Sto ch as ti c Processes with Applications to Finance, 2nd edn., Chapman & Hall/CR C Financial Mathematics Series, CRC Press, Boca Raton, 2013.
[8] T . Laudadio, N. Mastro nardi, P. Van Dooren, Computing Gaussian qua dra tu re rules with high rel at ive accuracy , Num er. Alg orithms 92 (2023) 767–793,
https://doi .org /10 .1007 /s11075 -022 -01297 -9 .
[9] D.S. Lubinsky , Which we i gh ts on R admit Jackson theor ems?, Isr . J. Math. 155 (2006) 253–280.
[10] G. Mastroianni, G.V. Milo v ano vi ´
c, Interpolation Processes – Basic Theory and Applications, Spring er Monographs in Mathematics, Spring er , Berlin,
German y , 2008.
[11] G. Mastr oianni, G.V. Milo vano vi ´
c, I. Notarange lo , Gaussian q uad rat ure rules with an exponential weig h t on the real semiaxis, IMA J. Num e r . Anal. 34
(2014) 1654–1685.
[12] G. Mastroianni, G.V. Milovano vi ´
c, I. Notara ngel o, A Ny ström method for a class of Fre dh o l m integr al equations on the real semiaxis, Calcolo 54 (2017)
567–585.
[13] G. Mastroianni, I. Notara ngel o, A Lagrange-type projector on the re al line, Math. Comput. 79 (269) (2010) 327–352.
[14] G. Mastroianni, I. Notara ngel o, Some Fourier-type operators for functions on un vounded interv al, Act a Math. Hung. 127 (4) (2010) 347–375.
[15] G. Mastroianni, I. Notara ngel o, Poly nom ial appro ximation with an exponential we ig ht on the real semiaxis, Ac ta Math. Hung. 142 (1) (2014) 167–198.
[16] G. Mastroianni, I. Notara ngel o, Lagrange interpolation at Pollaczek–Laguerre zeros on the real semiaxis, J. Approx. Theor y 245 (1) (2019) 83–100.
[17] G. Mastroianni, M.G. Russo, Fourier sums in wei gh te d spaces of functions. A survey , Jaen J. Appro x. 1( 2 ) (2009) 257–291 .
[18] G. Mastroianni, J. Szabados, Polyn omi al appro ximation on infinite intervals with weights having inner zeros, Act a Math. Hung. 96 (3) (2002) 221–258.
[19] G. Mastroianni, J. Szabados, Dir ect and conv erse polynomial appr o ximation theorems on the re al line with weights having zer os, in: Fro nt i er s in
Interpolation and Appro ximation, De dicated to the Memory of A. Sharma, in: Pure Appl. Math. (Boca Rat on), vol. 282, Chapman & Hall/CR C, Boca
Raton FL, 2007 , pp. 287–306.
[20] G. Szegö, Orthogonal Poly nom ial s, 4th ed., American Mathematical Society , Pro vidence, RI, 1975.
[21] H.S. Wilf, Mathematics for the Phy sical Sciences, Wiley , New Yo r k , 1962.
9
