- Email: [email protected]
Nested quadrature rules for Cauchy principal value integrals
Journal of Computational North-Holland
and Applied
Mathematics
25 (1989) 251-266
251
Nested quadrature rules for Cauchy principal value integrals * Annamaria
PALAMARA
Dipartimento
di Matematica,
ORSI
Politecnico di Torino, 10129 Torino, Italy
Received 30 October 1987 Revised 22 November 1988 Abstract:
We consider
quadrature
bf----dx, (xl
formulas
for the numerical
evaluation,
with error estimate,
of integrals
of the form
- co
f u x-A
and we discuss them, pointing out the computational compare our numerical results with the ones of another
aspects. Then we present an automatic published integrator for Cauchy principal
Keywords:
rules, algorithms.
Cauchy principal
value integrals,
quadrature
integrator and we value integrals.
1. Introduction
In this paper we examine the problem of the construction of an automatic routine for the numerical evaluation, with error estimate, of Cauchy principal value integrals of the form V(x) f .x-A
dx y
-
cx
+CQ.
0)
Quadrature rules of interpolatory type, explicitly constructed to evaluate (l), can be of two different kinds, according to whether among the nodes we also include the point X or not. By interpolating the function f(x) on n distinct points { t,}l= 1 with a polynomial of degree n - 1, we obtain a formula of the type f
;isdx
= i
w,(A)f(t,)
+ R;‘(f)
(4
i=l
with degree of exactness n - 1, i.e. with R’,l’( f) = 0 whenever f(x) is any polynomial of degree n - 1. If in addition to the points { t, } :=i we also interpolate at the point t, = A, with h # t,, we obtain a formula of the type f
;,$dx
=A(X)f(h)
+ 5
w,(A)f(t,)
+ ~‘,~‘(f)
r=l
with degree of exactness at least n. * Work sponsored 0377-0427/89/$3.50
by the Minister0
della Pubblica
Istruzione
0 1989, Elsevier Science Publishers
of Italy.
B.V. (North-Holland)
252
A. Palamara
Orsi / Cauchy principal
value
If X = fj, for some i, we obtain a rule that involves the derivative f’(h) While the degree of (2) in general does not exceed n - 1, (3) can become fact, if d 3 n - 1 is the degree of exactness of the interpolatory rule
[4]. a “Gaussian”
I’ J’(x) dx = 2 v,F(t,)
(4)
i=l
“-1
rule; in
where the nodes { t, }ycl coincide with those of (3), then the degree of (3) is d + 1; hence when (4) is Gaussian we have d + 1 = 2n. In general, the computation of the coefficients of a formula of type (2) is not a simple matter and if it is not done carefully it may cause numerical cancellation when h is close to one of the t,‘s; furthermore it is fairly onerous. When the associated rule (4) is Gaussian, stable algorithm for the evaluation of the w;(A)‘s has been given (see [4]). In the computation of the approximation given by (3) we have numerical cancellation when h approaches one of the nodes (see, for example, [5,10]). Based on the above observations, to reduce the cancellation phenomenon Monegato [4] proposes a particular quadrature formula of type (3) to approximate (l), in which also the nodes depend on X. In this paper we assume this formula as basic integration rule for the evaluation of integrals of type (1) and we construct a suitable corresponding extension to estimate the error efficiently.
2. Basic quadrature formula We wish to compute the Cauchy principal value integral (l), where f(x) is such that the integral (1) exists. Monegato [4] suggests to split the interval (a, b) into two subintervals so that X coincides with the midpoint of one of them. Assuming, for example, that h > :( a + b), we set c = 2A - b and we obtain
‘ftX) fa
x-x
dx=
“dX)
J x-x
dx+
a
(likewise in the case X < $( a + b)). Because of the symmetry with respect can be reduced to the form
1 -dx. g(x)
f -1
‘dX) f (’ X-A
dx
to h of the interval
(5)
(c, b), the integral
between
c and b
(6)
x
For the evaluation of (6) Monegato [4] proposes the use of 2m-points Gauss-Legendre rule. Indeed, in this way X = 0 does not coincide with any node and, moreover, is a zero of each Legendre function of the second kind qZm(z), so that (see [4]) the 2m-points Gauss-Legendre rule
J-1’ F(x)
dx = E H,F(<,) I=1
A. Palamara
applied
directly
l f -1
to the function
g(x)
-dx x
g( x)/x
Orsi / Cauchy principal value
253
gives
z ii CH!/5,) ’ [ g(5i) -
(8)
EC-ti>]
i=l
where the remainder term vanishes whenever g(x) is a polynomial of degree 4m. For the evaluation of the integral on (a, c) Monegato [4] proposes a rule of type (3) and, in order to estimate (1) with the highest degree of exactness, he chooses as the nodes { t, }:=i in (3) those of formula (7) (with n = 2~2) so that it also has degree of exactness equalling 4~2. He thus deduces the following quadrature formula of type (3) constructed on 4m + 1 nodes dependent on A, with degree of exactness 4m: f
Q’=$+dx=A”,“f(X)
+ F
[&‘f(x,)
+A;~‘+~)]
i=l
where xi=+(C-a)&++(a+C),
z, = $(b - c)& + :(b + c),
Fly= :(c-a)Hi/(xi-h),
A(‘) = 0
A12’ = f&/5, >
ln
I5 I- ‘c”/p. i=l
This formula is, therefore, sufficiently stable and its coefficients modification of those of the underlying rule chosen. In the case X < i( a + b), we obtain an analogous formula.
are obtained
by a simple
3. Extended quadrature formula In order to obtain a realistic and efficient estimate of the error made in the evaluation of integral (l), we intend to compare the estimates given by the rule (9) and by a corresponding extended rule, with an optimal degree of exactness quadrature rule, which uses all the function values previously computed in (9). Let us consider the Kronrod extension to rule (7), with degree of exactness 6m + 1 (see, for example, [ 31) ’ F(x) J -1
dx = E aiF
+ ‘F&-(yJ
r=l
(10)
i=l
where the nodes { y, } fZl+’ are the zeros of the (2m + l)th-degree defined by the orthogonality relations 1
J
_l~2m(4%z+l
(x)x”
dx=O,
k=O,
l,...,
2m
where P2m( x) is the (2m)th-degree Legendre polynomial. In (5) to evaluate the integral on (a, c) we use a Kronrod with Gauss-Legendre nodes and we obtain
JQcsdx
= @‘f(X)
+ ‘c” Bl”)f(x,) i=l
+ ‘F1 i=l
Stieltjes polynomial
E2m+l(~)
(11)
extension
c,c’)f(t;)
of the rule of type (3)
(14
254
A. Palamara Orsi / Cauchy principal value
where
x,=~(c-a)~,+:(a+C), B(l)
1
=
I
2
B(l) = ln 0
t, = :(c - a)y, + &2 + c),
dc- 4 x,-x
c(l)=
1
NC-4
’
2
ti-h
’
’
C-X
1 I a-A
i=l
1=7
This formula has degree of exactness 6m + 2. In order to extend the quadrature formula concerning the integral on (b, c), we need to construct, if it exists, an optimal extension to formula (8). Rabinowitz [9], taking into consideration the Kronrod extensions to Gaussian integration rules for the evaluation of Cauchy principal value integrals of the form (l), introduces a new extension, for the case A = 0, which has the advantage of not involving the evaluation of the derivative of the integrand function. In fact he extends the 2m-points Gaussian rule adding 2m + 2 points, rather than 2m + 1 as in the Kronrod case; in this way he avoids the point x = h = 0. The extension of the formula (8) he constructs is of the type
with degree of exactness 6m + 3, where the nodes { JJ }F!+,’ are the zeros, whenever they exist, of the (2m + 2)th-degree Stieltjes polynomial i 2m+2(~) which satisfies the orthogonality conditions 1
(x)xk=O, J_lP,mbP2m+2
k=O,l,...,
04)
2m+l.
The coefficients of the polynomial ,!?2m+2(~) are determined, orthogonality conditions (14) which characterize it. The polynomial E *m+*(x)= X&m+lb) for any constant
up to a constant
by the
05)
- a
a is such that
1
(x)x” J-1 p2mb)E2m+2 Moreover
factor,
dx=O,
k=O,
l,...,
2m-1.
(16)
if we define
&n+z(x)= xE2m+lb)-
(17)
a”
where 1 a”=
/
h7(4L+lbb2m+1
Wj-;lf’d4x2n’
dx,
-1
this polynomial satisfies (14) and Monegato [2] proved that it is the unique one (up to a constant factor) satisfying it. We can construct (13) when the polynomial I&,,+2 (x) has 2m + 2 zeros which are real, distinct, located in the interior of (- 1, 1) and separated by the zeros of P2m(~); in this case for
A. Palamara Orsi / Cauchy principal value
the weights of (13) we can easily find the following
t=
tia,E2m+*(Sr)
-
t,E2m+1(5i)
CiH,
_~
,
i=
expressions:
1, 2,...,2m,
22”bn~,
s;= k,*,P2,(j,)[E,,+,(y,)
255
+y;EE;,+,(y,)] ’
I=1923.-.T2m+27
(18)
where the leading coefficient of E,,+,(x) is 22” and the leading coefficient of Pz,,,(x) is k2*, = [l . 3. 5 . . _. *(4m - 1)]/(2m)!. However the new quadrature rule (13) does not exist for any m; in the examined case of weight function equal to 1, we need to make m odd so that I? 2n,+2(~) can have 2m + 2 real and distinct zeros in (- 1, 1). Furthermore Rabinowitz [9] notices that, in this situation, these zeros exist, but, due of the decreasing of the smallest positive zero of E”2,,,+2(~) (which we denote with Jr) while m increases, the term of (13) T = (s,/Jl)[g( j+) - g( -y,)] can lead to the loss of significant figures. Nevertheless, from the numerical experiments we have performed, we notice that, if one has numerical cancellation in the quantity g( jl) - g( -jr), this is compensated by the fact that, in general, even if jj, is very small, the value of T is however negligible with regard to the value of the total sum and, therefore, the “erroneous” significant digits do not contribute to the formation of the final value. Anyway, in order to reduce the eventual numerical instability of (13) we have used a new polynomial of type (15) with m odd, in which the parameter a = a(m) is such that the polynomial still has 2m + 2 real and distinct zeros all contained in ( - 1, 1) (i.e., xTE2,,,+ 1(XT ) < a -c 0, where x,* 1s the smallest positive abscissa of the minimum points of xE2,,+,( x)), but the distance of his smallest positive zero from the origin is greater than J, (3, is the smallest positive zero of the polynomial E”2m+2(~) with the same m). Among the values of a for which these conditions are verified we have chosen the value ~7= 0.4x:E2,+,(x;)
-=Ia”,
because, based on numerical experiments done with m = 3(2)25, we have obtained the amplification factor K(a) with a variable, defined (in the case a = a”) as
the result that
(19) has a minimum when a = ii. Let us indicate this new polynomial ~2m+zb>
=
XE2m+lb)
by
(20)
-a
and its zeros by { j, } fJJc*. We can derive the new formula
(21) (where the symbols have a clear meaning) (16) is satisfied.
with degree of exactness
6m + 1, since in this case only
256
A. Palamara
Therefore f
we obtain
the following
y&lx =B(yf(A)
+
Orsi / Cauchy principal
global rule, in the case A > :( a + b):
‘c” lpf(xJ
+
yf
2y Ci’l’f(t,) r=l
i=l
+
ualue
B!2)f(Z,) +
2~2c;2)f(ui)
(22)
i=l
i=l
where u, = :(b - c)y, + :(b + c),
z, = :(b - c)& + :(b + c), B!2’ = 7,/& )
c,‘2’ = S,/jJl
and with the previous meaning of the other symbols. Formula (22) is constructed on 8m + 4 nodes and has degree of exactness We obtain an analogous formula in the case X < i( a + b).
6m + 1.
4. Numerical results For the values m = 3(2)25, we have determined the polynomial Elm+i( x) which satisfies (11) using the method described in [3]. In order to compute the above-mentioned polynomial, the polynomial P2m( X) and their first derivatives we have used Clenshaw’s algorithm applied to the corresponding finite expansions in Chebyshev polynomials. As initial approximations for computing the roots of E”2m+2(~) and E,,+,(x) by Newton’s method we have taken the midpoints between the zeros and the minimum points of E,,+,(x) evaluated previously. Tables 1 and 2 show the differences between formulas (13) and (21).
Table 1 Behaviour of the smallest positive zeros 9, and J, of the polynomials E”Zm+2(x) = x&,+,(x)- ii, respectively xEzm+ 1(x) m
d
a
Jl
d and E,,+,(x)
Yl
3
- 1.5D-2
4.2D-2
- lSD-1
5 7
- 5.7D-3 - 2.5D-3
2.OD-2 1.2D-2
- 9.5D-2 - 7.OD-2
1.4D-1 8.6D-2 6.3D-2
- i.4D-4
i.8D-3
- i.OD-2
i.8D-2
2;
Table 2 Amplification factors associated with formulas (8), (13), (21) N K
14 7.0
Z i?
62 9.9
78 10.4
102 10.9
10.1
30 8.5 15.6
20.3
22.4
25.1
6.5
8.0
9.5
10.0
10.5
=
257
A. Palamara Orsi / Cauchy principal value
In Table 1 we list the values of the constant a” which appears in the expression (17) of the Stieltjes polynomial &m+2(~) = xE~~+~(x) - a” whose zeros J, are nodes of the rule (13) and correspondingly the values of Z of the new polynomial (20) El,,,+ 2( x) = xEzm+ ,( x) - Z whose zeros y, are nodes of the r_ule (21). The latter polynomial E2,,,+2 is obtained by translation of the former ,!?Zm+2: acting on the constant a” we move the smallest positive zero away from the origin and we reduce the numerical cancellation in the corresponding quadrature rule. In Table 2, K, I? and K indicate the amplification factors associated with formulas (8) (13) and (21) respectively and N is the total number of nodes of the formula chosen. Then we have compared numerically the rule of type (9) with the corresponding one of type (22) (that is with m having the same value) applying them to integrals of type (1). The examples considered give us a measurement of the better approximation furnished by (22) with respect to that of (9) and show that using (22) we can obtain a reasonable estimate of the error in (9). In Table 3 we give abscissas and weights for a representative case (m = 3) of the new formulas (9) and (22). They are given for the interval (- 1, 1); because of symmetry only the positive abscissas and their corresponding weights are given. In Table 3 the symbols have the following meaning: 5, = positive nodes H, = weights of the (Y; = weights of the corresponding y, = positive zeros p, = weights of the corresponding r, = weights of the j, = positive zeros s, = weights of the
Table
of the 6-point Gauss-Legendre rule, 6-point Gauss-Legendre rule, 13-point Kronrod extension of the 6-point Gauss-Legendre to the abscissas Ei, of the 7th-degree Stieltjes polynomial ET(x), 13-point Kronrod extension of the 6-point Gauss-Legendre to the abscissa 0.0 and to the abscissas t,, 13-point new formula (21), corresponding to the abscissas of the 8th-degree new polynomial E,(x) = xE,(x) - ii, 13-point new formula (21), corresponding to the abscissas
3
Abscissas
and weights for the new formulas
(9) and (22), in the case m = 3
5,
H,
a,
c
0.2386191860831969
0.4679139345726911
0.2337708641169944
0.03756616183158088
0.6612093864662646
0.3607615730481386
0.1810719943231376
0.2092717906299459
0.9324695142031520
0.1713244923791704
0.08369444044690674
0.06260430386803677
Y,
P,
Y,
sz
0.0
0.2410725801734648
0.1375457357181660
0.2596716622719610
0.4631182124753046
0.2132096522719623
0.4264557065135209
0.2542968856122309
0.8213733408650279
0.1373206046344469
0.8362466005176528
0.1368677681526986
0.9887032026126789
0.03039615411981977
0.9831454575784491
0.03972142763354597
rule,
rule, <,, J,.
A. Palamara Orsi / Cauchy principal value
258
Table 4 Relative errors in the 13-point and 2%point rules ((9) and (22), respectively) for the integral H,(X) = f’ ,cos(5x +2)/(x - A) dx A
H,(X)
13 points
28 points
13 points with respect to 28 points
0.01 0.05 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 0.95 0.99
-
2.8D-5 1.7D-5 8.8D-6 2.3D-6 3.4D-7 5.1D-8 5.6D-9 1.7D-8 3.5D-7 5.7D-6 3.8D-5 2.7D-5 2.1D-5
3.2D-12 1.4D-12 4SD-13 4.5D-14 2.7D-15 4.5D-16 1.7D-16 1.9D-16 3.8D-16 4.6D-15 7.1D-14 8.1D-14 9.1D-14
2.8D-5 1.7D-5 8.8D-6 2.3D-6 3.4D-7 5.1D-8 5.6D-9 1.7D-8 3.5D-7 5.7D-6 3.8D-5 2.7D-5 2.1D-5
2.74829889791444 2.39910773684172 1.82773804017123 0.376974897586441 1.18225431509743 2.47184332651756 3.17940827525908 3.13424923411046 2.34722742476011 1.00041044109658 - 0.642524886674366 - 1.57953528558164 - 2.92679456816624
Table 5 Relative errors in the (4m + l)-point with m = 3, 5, 25, for the integral
and (8~ +4)-point H,(h)
= f’,dg/(x
rules ((9) and (22), respectively), - X) dx
A
H,(X)
13 points
28 points
13 points with respect to 28 points
21 points
44 points
21 points with respect to 44 points
101 points
204 points
101 points with respect to 204 points
0.01 0.05 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 0.95 0.99
- 0.031415926536 -0.157079632679 -0.314159265359 -0.628318530718 - 0.942477796077 - 1.256637061436 - 1.570796326794 - 1.884955592154 -2.199114857513 - 2.513274122872 - 2.827433388231 - 2.984513020910 - 3.110176727054
2.9D-2 8.7D-3 4.3D-3 1.9D-3 l.lD-3 6.9D-4 4.4D-4 2.6D-4 1.3D-4 2.3D-5 6.2D-5 5.4D-5 4.9D-3
6.4D-3 1.4D-3 6.7D-4 3.1D-4 1.9D-4 1.3D-4 8.8D-5 6.2D-5 4.3D-5 2.7D-5 1.2D-5 3.7D-6 1.4D-5
2.3D-2 7.4D-3 3.6D-3 1.6D-3 9.2D-4 5.7D-4 3.5D-4 2.OD-4 8.7D-5 3.6D-6 7.4D-5 5.OD-5 4.9D-3
9.7D-3 2.2D-3 l.OD-3 4.6D-4 2.7D-4 1.7D-4 l.lD-4 6.4D-5 3.2D-5 6.4D-6 1.7D-5 2.8D-5 5.7D-4
1.5D-3 3.OD-4 1.4D-4 6.6D-5 4.OD-5 2.7D-5 1.9D-5 1.3D-5 8.9D-6 5.4D-6 2.3D-6 5.3D-7 3.3D-6
8.3D-3 1.9D-3 8.9D-4 4.OD-4 2.3D-4 1.4D-4 8.7D-5 5.OD-5 2.3D-5 9.4D-6 1.9D-5 2.8D-5 5.7D-4
l.OD-4 2.OD-5 9.4D-6 4.2D-6 2.4D-6 lSD-6 9.6D-7 5.8D-7 3.OD-7 6.3D-8 1.5D-7 2.7D-7 4.1D-7
9.6D-6 1.9D-6 9.OD-7 4.1D-7 2.5D-7 1.6D-7 l.lD-7 7.5D-8 4.8D-8 2.6D-8 6.4D-9 4.8D-9 1.8D-8
9.2D-5 l.SD-5 8.5D-6 3.8D-6 2.2D-6 1.3D-6 8.5D-7 5.OD-7 2.5D-7 3.7D-8 1.6D-7 2.6D-7 3.9D-7
we have
obtained
In Tables 4 and 5 we report some numerical results above-mentioned rules to the following two integrals:
by applying
the
dx = --TX. In Tables 4 and 5 we report the relative accuracy obtained using the rule (9) with 4m + 1 nodes, with respect to the exact value, the relative accuracy obtained using the rule (22) with
A. Palamara Orsi / Cauchy principal value
259
8m + 4 nodes, with respect to the exact value and the relative accuracy obtained using the rule (9) with respect to the value obtained using the rule (22). We assume the last error as error estimate of integral approximation obtained using the rule (9); we can also assume it as an overestimation of the error of integral approximation obtained using the rule (22). In Table 4 we can see that the value of the integral supplied by (22) is already accurate to nearly machine precision with 28 nodes; this is difficult to realize if one only considers the estimation of the error which is reported in the last column. On the contrary, in Table 5 we can see that the extended formula is not much better than the basic formula. This different behaviour is due to the different degree of smoothness of the two integrands; indeed in the first case the integrand function is much more smooth and easier to be approximated by a polynomial than in the second case.
5. Adaptive routine For the evaluation of (1) for a single value of A E (a, b), supposing the function f(x) has a low degree of smoothness in the interval [a, b], it is convenient to apply an adaptive strategy analogous to that used by Piessens in [6] to evaluate the integral labf( x) dx and in [7] to evaluate (1). Indeed an adaptive strategy locates a few nodes in the s&intervals in which the integrand is regular and many nodes in the neighborhood of difficult spots of the integrand. If, however, we want to evaluate (1) for a set of values of A and for smooth functions f(x), then the use of Gauss-Kronrod rules for Cauchy principal value integrals will be more efficient, as Rabinowitz [8] pointed out. Let us describe the automatic globally adaptive integration scheme used to obtain an approximation y, for some preassigned relative tolerance tol, to the integral (1). Given the interval I = (a, b), the function f and the parameter A, let
G*,(I; f> K 4,+1(~; f) G&,+i(I; f) K&+2( 1; f) G4*m+i(I; f) K,*,*,,(I; f)
= 2m-point Gauss-Legendre rule applied to the function f( x)/( x - A) with x E [a, b] (with degree of exactness 4m - l), = K ronrod extension to the GZm( I; f) integration rule, = rule of type (3) where the nodes { ti } fr, coincide with those of the Zm-point Gauss-Legendre rule (with degree of exactness 4m), = Kronrod extension to the G&+i( I; f), = rule (9) (with degrees of exactness 4m), = rule (22).
Algorithm Step 1. Input: {a, b, A, f, m, to/}, where a < X -c b. Step 2. Let 1, = (a, b) Set ier = 0 (ier: error flag) Calculate est( I,) = K,*,*,,( 1,; f), Step 3.
err(h) = I G4*m*+lUl; f> - KiLL(4; If err( I,) < tel. 1e.st(Il) 1, then y = e.st(ll) and exit.
f> I
260
Step 4.
Step 5. Step 6.
Step 7.
Step 8.
Step 9. SteplO.
Step 11. Step 12.
Step 13. Step 14. Step 15.
Step 16.
A. Palamara Orsi / Cauchy principal value
Divide the interval 1, into two subintervals I1 and I*, eliminating the previous I, and adding a new interval I,, and set k = 1, n = 2. The subdivision of the initial interval Ii must be made as follows: if h @ I1 then Ii is bisected; if a < h < (a + b)/2 then Ii is subdivided at the point (A + b)/2; otherwise subdivision occurs at the point (A + a)/2. For I = k and I= n do Let ((Y, p) denote the open interval I,. If h E I,, then calculate est( I,) = &*m*++(I,; f), err (IA = I G4*m+1U~; f) - GZ+dC f> I and go to Step 9 of the algorithm. Else if X>/3 and IX--/31 -C Ip-c~l or if h < (Y and A G I, but it is near I,), then calculate est( II) = K4*m+2(I,; f), err(h) = I G’k+d4; f> - KEm+2U1; f> and go to Step 9 of the algorithm. Else (i.e. if A 6 I, and it is far from II) calculate est( I!) = K,, +1( I,; f),
I A - a I < I j3 - a I (i.e. if
I
err(h) = I %M; f> - &m+I(4; f) I. End of Z-loop. ) the initial interval (a, b) is divided into n subintervals Ij, At step n (n=2,3,... j = 1, 2,. .-, n, for each of which the estimates est( I,), err( I,) are available, where these values are carried over from the previous steps or have been computed for the subintervals just generated. If the maximum number of integrand evaluations allowed has been achieved, then set ier = 1 and go to Step 14 of the algorithm. Else if the length of some interval I, is getting too small, then set ier = 2 and go to Step 14 of the algorithm. Else if tel. IC’Jclest( I,) I < E (6: relative machine accuracy), then if C/“,,err(Ij) G tol, then go to Step 14 of the algorithm. Else if E3=,err(Ij) G tol- IC’J=,est(Ij) 1, then calculate y = CT= ,est( I,) and exit. Else split up the interval I, on which err( I,) = maxi <, i ,,err( I,) into two subintervals using the same strategy described at Step 4 of the algorithm; append a new interval I n+l to the list of subintervals and replace the old Ik by a new one. Set n = n + 1 and go back to Step 5 of the algorithm.
We remark that the strategy indicated at Step 4 of the algorithm is used to avoid a subdivision occurring at the abscissa A of the non-integrable singularity or too close to the abscissa X; in this latter case, in fact, some subtractive cancellation in the computation of the total integral estimate would happen.
A. Palamara Orsi / Cauchy principal value Table 6 Integral H,(h) = f’IeX/(x - A) dx; number of function evaluations required Quadpack at requested relative accuracy to1 and respective really achieved relative
x
0.01 0.05 0.1
0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 0.95 0.99
Present method
Quadpack
to1= 10-12
to1=10-‘2
261
by the present accuracies
NE,
AERR,
NE,
AERR 2
28 28 28 28 28 28 28 28 28 28 28 28 28
4.2D-16 4.3D-16 5.6D-16 2.4D-16 4.1D-16 3.4D-16 0.0 l.lD-15 2.4D-16 6.6D-16 5.8D-16 4.5D-16 5.OD-16
25 25 25 25 25 25 25 25 25 25 25 25 25
8.4D-16 8.6D-16 5.6D-16 2.4D-16 5.5D-16 1.7D-16 1.2D-16 1.9D-15 l.lD-15 9.2D-16 9.2D-16 1.5D-16 6.7D-16
method
and
by
Concerning Step 8 of the algorithm, we remark that when A is far from the integration interval, we do not use the formula G -&+l(I; f) to avoid the computation of A( A) which is much more expensive than the test regarding the nearness or not of A to the integration interval. In Tables 6-10 some numerical results obtained by applying the routine which implements the previous algorithm and the analogous routine QA WC-QA WCE of Piessens [7] to the following integrals are shown: H,(h)=fl~dx,
NV
= f,
H,(h) = H,(A)
-l
f_l,“‘,IT2 dx,
= f_l, “‘,‘_I;
X-h
We asked the routines
-l
’ dx,
1 1x - 0.65 I’.’
to integrate
-l
-l
O
each test function
to relative
accuracies
to1 = 10-3, 10-6, 10-9, lo-12. The results indicate that the two routines have about the same average reliability (if we define a program to be reliable when the result obtained achieves the requested tolerance), but the first one appears more efficient; in fact the number of integrand evaluations required to calculate a Cauchy principal value integral with a given accuracy is significantly inferior for functions f(x)
28
28
28
28
28
28
28
28
28
28
28
28
0.05
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
0.95
0.99
9.1Dp14
8.lD-14
7.5D-14
6.OD-15
0.0
0.0
4.2D-16
3.6D-16
153
1.2D-16
96 138 112 112
2.7D-15
2.4Dp15
8.4D-16
3.OD-16
70
70
70
70
28
28 70
112
4.2D-16 70
112
3.6D-16
0.0
111
2.1Dp15
137
154
1.9Dp16
1.3Dp15
153
NE,
0.0
AERR,
0.0
28
28
28
70
2.1D-15
70
4.4D-14
70
4.5D-13
1.4D-12
70
NE,
3.2D-12
AERR,
NE,
28
tol=lO~”
number
of function
1.5Dp16
7.OD-16
3.6Dp15
2.7D-15
0.0
1.4D-16
1.4Dp16
3.6D-16
0.0
1.3Dp15
0.0
1.9D-16
1.6D-16
AERR,
25 25
5.9Dp16 0.0 3.6Dp16
324 222
25
7.OD-16
190
25
25
3.5D-15
190 1.5Dp16
25
2.2Dp15
164 190
25
25
2.8D-16 0.0
180
25
205
205
2.8Dp16
25
0.0
299
298
25
25
3.7D-16
289
25
3.2D-16
0.0
5.6D-16
4.OD-15
5.5D-15
1.7D-15
1.8D-15
7.OD-16
1.6D-15
l.lD-15
lSD-15
2.1Dp15
7.4Dp16
1.5D-15
method
25
65
65
65
65
65
65
65
65
65
65
65
65
NE,
0.0
3.1Dp8
7.4Dp8
4.3Dp8
1.5D-8
8.6D-9
5.4D-9
3.OD-9
1.6Dp9
2.6D-8
8.9D-9
7.7Dp9
7.3D-9
AERR,
tol= lo-”
by the present
AERR,
NE,
AERR,
NE, 247
Quadpack
required
to/=10-’
evaluations
t0l=10~”
relative accuracies
tol=lO~’
+ 2)/( x - X)dx;
really achieved
tol= 10-3
Present method
0.01
x
rol and respective
= fL ,cos(5x
H,(X)
Integral
accuracy
Table 7
375
435
555
405
315
325
285
445
595
835
445
415
385
NE,
2.7D-10
1.3Dp10
3.OD-10
1.8D-11
1.6D-11
1.2D-9
4.5Dp10
1.9D-11
4.9Dp10
7.5D-10
6.1D-10
7.1Dp10
l.lD-9
AERR,
to/=lo-’
and by Quadpack
6565
8245
11695
11795
7475
6165
6225
8925
15055
15045
12745
10825
9325
NE,
tol=lo~‘*
relative
3.8D-12
1.9Dp12
4.1Dp12
4.4D-12
1.5D-12
3.OD-12
7.9D-13
4.ODp13
9.ODp13
4.2D-12
1.9D-12
3.7D-12
4.4D-12
AERR,
at requested
%-
x
3. R G’ k
:
9 *’ \ 0 !S E
1.9Dp4
1.3Dp4
8.8D-5
6.2D-5
4.3D-5
2.7D-5
1.2Dp5
3.7D-6 2.3Dp6
28
28
28
28
28
28
28 112
0.4
0.5
0.6
0.7
0.8
0.9
0.95 0.99
2.5Dp4
70
28
0.2
0.3
1.7Dp8
9.5D-9
1.4Dp8 l.lD-8
513
555 613
1.9Dp8
5.5Dp9
497
455
439
1.6D-8
9.7Dp9
439
1.6D-8
6.3D-9
449
423
440
2.2D-9
475
385 395 405
1429 1606
5.ODp12 l.lDp11
893 951
8.5D-15 5.ODp15
345
5.3Dp16 4.lD-15
1345
335
1387
295
9.1D-15
1.2D-11
1370
1.3D-11
819 6.6Dp12
1354
1.3D-11
860
325 195
4.8Dp15 2.1Dp15
7.6D-15
315
315
375
405
495
877
1380
l.lD-11 1.3D-11
846 819
2.4D-15
8.8Dp15
7.1Dp15
l.lD-15
1.3D-14
861
1417 1416 1407
4.3Dp12 4.4Dp12
829
1.5D-12
881 896
1641 1562
1.9D-11
1.3Dp5
1771
1.8Dp13
908
70
7.8Dp10
985
0.1
528
1.8D-8
l.lD-6
605
2.6D-5
164
X15 865
l.lD-4
775
765
725
685
715
715
705
735
765
825
915
NE,
1.3Dp7 3.7D37
2.7D-7
3.8D-X
8.6Dp8
1.5D--3
9.2Dp8
1.4D-7
2.4D-7
1.9Dp7
1.2Dp8
2.7Dp7
method
1.3Dp9 1.4Dp9 1.4D-9
1515 1565
4.4Dp9 3.XDp9
l.OD-9
l.ODp9
6.7Dp10
1.3Dp9
5.4D-10
1505
1465
1425
1455
1355
4.8Dp9
5.5D-9
4.5D39
2.8D-9
4.OD-9
1545
4.ODp11 5.2Dp9
4.3D-11 1445
1685 1535
4.3D-11
1925
1.3Dp10
1.5Dp10 3.8D-11
2735
AERR,
y
NE,
toI = 10
and by Quadpack
9.7Dp11
9.4D-11
X.ODp11
l.lDp10
AERR,
rol = lo-”
by the present
AERR,
241
NE,
0.05
AERR,
0.01
NE,
AERR,
NE,
AERR,
NE, NE,
tol=lO~’
rol=lo~‘*
tot-10-h
required
rol=lO~~ tol=lO~’
evaluations
Quadpack
AERR,
of function
relative accuracies
- X)d:. 1. number
really achieved
Present method
x
to/ and respective
= f’,js/(~
H,(X)
Integral
accuracy
Table 8
4.2Dp12 4.4Dp12 4.9Dp12
11515 11645 11725
4.ODp12
10895 11535
1.3Dp12 3.4D-12
11885
l.lDp12
14465
2.7D-12
13615 10915
1.5D-14 9.7Dp14
15055
3.5Dp13 2.7D-13
15055
3.7Dp13
15055 15055
AERR
lo-‘*
relative
NE,
tol=
at requested
898
1.4D-4
2.8D-16
3.6D-3
2.7D-5
278
195
153
112
28
138
138
138
28
70
70
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
0.95
0.99
1.4D-5
0.0
3.8Dp4
1.9D-5
l.lD-16
1.7D-13
7.7D-5
70
242
216
242
242
138
28
242
366
195
450
0.0
2.213-8
8.4D-7
4.1D-7
7.8D-8
l.lDp16
1.7Dp13
3.1D-7
1.3D-8
2.8D-16
8.6Dp8
70
0.0
9.1D-11
7.7Dp10
346 372
1.4Dp9
398
164
2.4Dp13 0.0
502 96
2.4D-15
6.7Dp13
1.3Dp14
l.lD-16
2.1D-16
6.7Dp14
5.3Dp14
4.2Dp16
8.1D-14
1.6DpI 3
1.2D-13
AERR,
evaluations
580
606
580
l.lD-16 6.1D-11
96
528
782
273
840
2.1Dp16
7.6D-11
6.5D-11
2.8Dp16
424
138
70
398
548
195
606
8.4D-11
1194
1.3Dp11
7.2D-8
2.2D-8
716
5.4Dp12
534
700
908
l.lDp4
320
6.5D-4
570
0.01
NE,
AERR,
NE,
AERR,
NE,
NE,
AERR,
toi=10-’
to1 =10p6
rol=lO~”
1x1)/(x-A)d x; number of function really achieved relative accuracies
rol=lO-’
Present method
0.05
x
Table 9 Integral H3(X)=fl,(x4+ accuracy to/ and respective
65
65
65
225
225
135
65
195
285
65
335
445
605
NE,
2.6D-6
9.9Dp6
1.5D-4
4.5Dp6
l.ODp7
8.9D-9
l.lD-3
8.9Dp8
4.1Dp6
l.lDp2
6.7D-6
1.3D-6
1.4D-6
method
225
315
345
405
345
135
345
345
435
465
485
595
695
NE,
4.ODp6
2.2Dp10
9.OD-9
5.1 D-8
8.6Dp9
8.9Dp9
4.6Dp9
4.50-9
2.OD-9
2.8D-9
l.OD-8
5.2Dp9
6.2Dp8
AERR,
tol= 10-h
by the present
AERR,
tot =10-s
Quadpack
required
495
555
555
765
615
195
555
555
915
885
1025
1345
1985
NE,
4.ODp6
7.5D-IO
6.3D-10
7.7Dp10
9.7Dp10
1.2D-9
5.3D-10
1.6D-9
1.7D-10
8.2Dp10
7.2Dp10
7.8D-10
9.6D-10
AERR,
to/ = lomq
and by Quadpack
1905
2445
3045
5475
3195
2325
2565
2415
7335
7245
10175
14425
15125
NE,
relative
4.ODp6
4.9Dp12
5.9D-12
7.OD-12
5.2Dp12
4.8Dp12
4.3Dp12
3.1D-12
1.6Dp12
2.ODp12
2.5D-12
4.2Dp12
l.OD-11
AERR,
to1=10-‘2
at requested
+
ti iii-
t
$.
2 3.
\ 0 e %
i? a’ 0 a __
2 F
6.3D-4
5.1Dp4
3.1D-4
8.3Dp5
1.2D-3 2.4D-7
2.ODp16 5.3D-5
l.lD-6
2.2Dp4
l.OD-7
4.7D-4
69
69
28
112
69 277
2X 278
221
112
138
70
0.2
0.3
0.4
0.5 0.6
0.65 0.7
0.X
0.9
0.95
0.99
190
216
242
429
2X 408
256 319
257
121
147
199
450 39X
6.3D-7
4.ODp8
466
808
3.2D-9
5.9D-7
28 725
2.OD-16 1.6D-7
413
2.4Dp7 495 620
448
1.6Dp7 l.OD-8
449
6.4Dp9
422
449
4.5D-9
1.7D-7
6.9D-7
277
0.1
174
3.7Dp7
l.lDp7
96
0.05
173
1.5D-5
69
0.01
rol = lo-
‘*
l.ODp9
6.7Dp12
2.7Dp10
9.5D-12
1.5D-15 l.OD-13 1.3Dp12 1.5Dp13
846 700 752
2.2Dp13
0.0
1593
96
7.4D-14 6.1D-15
X90 1000 1354
l.lD-13
808
2.ODp16
9.6D-16 4.1D-15
2.1Dp13
734 843
1.4Dp12
880
5.XDp15
605
AERR,
6X3
NE,
2.6D-11
2.5D-11 X.4Dp12
1.6D-10
6.ODp12
X.2Dp12
4.5D-10
7.3D-12
3.3D-11
AERR,
195 25
225
355 215
385
145
65 235
155 65
185
155
25
NE,
7.1Dp7 4.5Dp4
2.2D-7
2.3Dp7
4.6Dp8
4.OD-3 9.9D-4
1.6Dp6
1.15Dp3
255
285
345
615
545
1145
365
395
355
285
255
275
245
155
NE,
to/-l0
1.2D-8
2.3D-8
1.6D-8
7.3D-7
2.1D-8
6.1Dp7
1.5D-6
4.6Dp10
6.2Dp9
1.4Dp9
1.3Dp8
2.7Dp8
2.1Dp8
1.3Dp6
AERR,
oh
by the present method
4.ODp6
1.9D-9
2.1D-6
7.7Dp7
l.lDp4
AERR,
tol=lO_’
ro/= lo-” NE,
AERR,
AERR,
ro/=lO_”
NE,
to/=lOK’
Quadpack
required
Present method
NE,
A
x. number of function evaluations
really achieved
relative accuracies
IO/ and respective
accuracy
10
Integral H4(X)=fbI_~-0.651’-‘/(x-X)d
Table
10
7.XD-10 5.3D-10
555
6.ODp 10
3.4Dp9
2.OD-9
6.6D--10
7.6Dp10
7.7D-
685
745
4505
1405
19x5
915
735
1.3Dp9
5.6Dp10 725
5.2Dp 10 715
X.3DpIO
l.ODp9
2.XD- 11
675
5X5
5X5
2.2Dp10 1.3D-12 1.9Dp12 2.lDp13
15075 9885 5585 2995
l.lDp11
l.OD-12 15105
7735
2.XDp12 4.3Dp12
76X5 9035
1.9D-12 3.7Dp12
8525 7825
2.ODp12
3435 7755
3.3Dp12 2.6Dp12
2565
3.1D-12
AERR,
relative
1935
NE,
425
to/=10-‘2
AERR;
NE,
at requested
tol= 10-9
and by Quadpack
A. Palamara Orsi / Cauchy principal value
266
having a low degree of smoothness and for high accuracies. Only for well-behaved integrands and low accuracies the two routines, generally, are competitive. In particular, we can notice the advantage of the former routine over the latter, also for low accuracies, in the case of the integral H4( X) with X = 0.65 in which the integrand is not singular (see Table 10). Moreover, for our test integrals, with the present method the actual error in the computed values is almost always at least one order of magnitude smaller than the requested error; only in a few cases (see Table 9: toZ= 10e3, h = 0.9; tol= 10p9, h = 0.8. Table 10: tol= 10p3, X = 0.5; to1 = 10-9, A = 0.99; to1 = lo-‘*, X = 0.05 and X = 0.95) it is a little higher (by a factor less or equal to 3.6) than the requested tolerance. Also with Quadpack quite often the actual error is smaller than the requested error, though the anomalous cases are more than a few. In Tables 6-10, NE, and NE, denote the total number of function evaluations while AERR, and AERR, denote the actual relative error in the computed values of the integrals, using our routine with m = 3 and QA WC-QA WCE routine respectively. All computations reported on this paper were carried out in double precision on the PDP-11/23+ computer which has about 16-digit accuracy in double precision.
Acknowledgment I gratefully
acknowledge
G. Monegato
for his constant
guidance
during
this work.
References [l] G. Monegato, A note on extended gaussian quadrature rules, Math. Cump. 30 (1976) 812-817. [2] G. Monegato, On polynomials orthogonal with respect to particular variable-signed weight functions, ZAMP 31 (1980) 549-555. [3] G. Monegato, Stieltjes polynomials and related quadrature rules, SIAM Rev 24 (1982) 137-158. [4] G. Monegato, The numerical evaluation of one-dimensional Cauchy principal value integrals, Computing 29 (1982) 337-354. [5] D.F. Paget and D. Elliott, An algorithm for the numerical evaluation of certain Cauchy principal value integrals, Numer. Math. 19 (1972) 373-385. [6] R. Piessens, An algorithm for automatic integration, Angew. Inform. 9 (1973) 399-401. [7] R. Piessens, E. de Doncker-Kapenga, C.W. ijberhuber and D.K. Kahaner, QUADPACK: a Subroutine Package for Automatic Integration (Springer, Berlin, 1983). [8] P. Rabinowitz, A stable Gauss-Kronrod algorithm for Cauchy principal-value integrals, Comp. & Maths. Appl. 12 B (5/6) (1986) 1249-1254. [9] P. Rabinowitz, Gauss-Kronrod integration rules for Cauchy principal value integrals, Math. Comp. 41 (1983) 63-78. [lo] I.H. Sloan, The numerical evaluation of principal value integrals, J. Comput. Phys. 3 (1968) 332-333.
and Applied
Mathematics
25 (1989) 251-266
251
Nested quadrature rules for Cauchy principal value integrals * Annamaria
PALAMARA
Dipartimento
di Matematica,
ORSI
Politecnico di Torino, 10129 Torino, Italy
Received 30 October 1987 Revised 22 November 1988 Abstract:
We consider
quadrature
bf----dx, (xl
formulas
for the numerical
evaluation,
with error estimate,
of integrals
of the form
- co
f u x-A
and we discuss them, pointing out the computational compare our numerical results with the ones of another
aspects. Then we present an automatic published integrator for Cauchy principal
Keywords:
rules, algorithms.
Cauchy principal
value integrals,
quadrature
integrator and we value integrals.
1. Introduction
In this paper we examine the problem of the construction of an automatic routine for the numerical evaluation, with error estimate, of Cauchy principal value integrals of the form V(x) f .x-A
dx y
-
cx
+CQ.
0)
Quadrature rules of interpolatory type, explicitly constructed to evaluate (l), can be of two different kinds, according to whether among the nodes we also include the point X or not. By interpolating the function f(x) on n distinct points { t,}l= 1 with a polynomial of degree n - 1, we obtain a formula of the type f
;isdx
= i
w,(A)f(t,)
+ R;‘(f)
(4
i=l
with degree of exactness n - 1, i.e. with R’,l’( f) = 0 whenever f(x) is any polynomial of degree n - 1. If in addition to the points { t, } :=i we also interpolate at the point t, = A, with h # t,, we obtain a formula of the type f
;,$dx
=A(X)f(h)
+ 5
w,(A)f(t,)
+ ~‘,~‘(f)
r=l
with degree of exactness at least n. * Work sponsored 0377-0427/89/$3.50
by the Minister0
della Pubblica
Istruzione
0 1989, Elsevier Science Publishers
of Italy.
B.V. (North-Holland)
252
A. Palamara
Orsi / Cauchy principal
value
If X = fj, for some i, we obtain a rule that involves the derivative f’(h) While the degree of (2) in general does not exceed n - 1, (3) can become fact, if d 3 n - 1 is the degree of exactness of the interpolatory rule
[4]. a “Gaussian”
I’ J’(x) dx = 2 v,F(t,)
(4)
i=l
“-1
rule; in
where the nodes { t, }ycl coincide with those of (3), then the degree of (3) is d + 1; hence when (4) is Gaussian we have d + 1 = 2n. In general, the computation of the coefficients of a formula of type (2) is not a simple matter and if it is not done carefully it may cause numerical cancellation when h is close to one of the t,‘s; furthermore it is fairly onerous. When the associated rule (4) is Gaussian, stable algorithm for the evaluation of the w;(A)‘s has been given (see [4]). In the computation of the approximation given by (3) we have numerical cancellation when h approaches one of the nodes (see, for example, [5,10]). Based on the above observations, to reduce the cancellation phenomenon Monegato [4] proposes a particular quadrature formula of type (3) to approximate (l), in which also the nodes depend on X. In this paper we assume this formula as basic integration rule for the evaluation of integrals of type (1) and we construct a suitable corresponding extension to estimate the error efficiently.
2. Basic quadrature formula We wish to compute the Cauchy principal value integral (l), where f(x) is such that the integral (1) exists. Monegato [4] suggests to split the interval (a, b) into two subintervals so that X coincides with the midpoint of one of them. Assuming, for example, that h > :( a + b), we set c = 2A - b and we obtain
‘ftX) fa
x-x
dx=
“dX)
J x-x
dx+
a
(likewise in the case X < $( a + b)). Because of the symmetry with respect can be reduced to the form
1 -dx. g(x)
f -1
‘dX) f (’ X-A
dx
to h of the interval
(5)
(c, b), the integral
between
c and b
(6)
x
For the evaluation of (6) Monegato [4] proposes the use of 2m-points Gauss-Legendre rule. Indeed, in this way X = 0 does not coincide with any node and, moreover, is a zero of each Legendre function of the second kind qZm(z), so that (see [4]) the 2m-points Gauss-Legendre rule
J-1’ F(x)
dx = E H,F(<,) I=1
A. Palamara
applied
directly
l f -1
to the function
g(x)
-dx x
g( x)/x
Orsi / Cauchy principal value
253
gives
z ii CH!/5,) ’ [ g(5i) -
(8)
EC-ti>]
i=l
where the remainder term vanishes whenever g(x) is a polynomial of degree 4m. For the evaluation of the integral on (a, c) Monegato [4] proposes a rule of type (3) and, in order to estimate (1) with the highest degree of exactness, he chooses as the nodes { t, }:=i in (3) those of formula (7) (with n = 2~2) so that it also has degree of exactness equalling 4~2. He thus deduces the following quadrature formula of type (3) constructed on 4m + 1 nodes dependent on A, with degree of exactness 4m: f
Q’=$+dx=A”,“f(X)
+ F
[&‘f(x,)
+A;~‘+~)]
i=l
where xi=+(C-a)&++(a+C),
z, = $(b - c)& + :(b + c),
Fly= :(c-a)Hi/(xi-h),
A(‘) = 0
A12’ = f&/5, >
ln
I5 I- ‘c”/p. i=l
This formula is, therefore, sufficiently stable and its coefficients modification of those of the underlying rule chosen. In the case X < i( a + b), we obtain an analogous formula.
are obtained
by a simple
3. Extended quadrature formula In order to obtain a realistic and efficient estimate of the error made in the evaluation of integral (l), we intend to compare the estimates given by the rule (9) and by a corresponding extended rule, with an optimal degree of exactness quadrature rule, which uses all the function values previously computed in (9). Let us consider the Kronrod extension to rule (7), with degree of exactness 6m + 1 (see, for example, [ 31) ’ F(x) J -1
dx = E aiF
+ ‘F&-(yJ
r=l
(10)
i=l
where the nodes { y, } fZl+’ are the zeros of the (2m + l)th-degree defined by the orthogonality relations 1
J
_l~2m(4%z+l
(x)x”
dx=O,
k=O,
l,...,
2m
where P2m( x) is the (2m)th-degree Legendre polynomial. In (5) to evaluate the integral on (a, c) we use a Kronrod with Gauss-Legendre nodes and we obtain
JQcsdx
= @‘f(X)
+ ‘c” Bl”)f(x,) i=l
+ ‘F1 i=l
Stieltjes polynomial
E2m+l(~)
(11)
extension
c,c’)f(t;)
of the rule of type (3)
(14
254
A. Palamara Orsi / Cauchy principal value
where
x,=~(c-a)~,+:(a+C), B(l)
1
=
I
2
B(l) = ln 0
t, = :(c - a)y, + &2 + c),
dc- 4 x,-x
c(l)=
1
NC-4
’
2
ti-h
’
’
C-X
1 I a-A
i=l
1=7
This formula has degree of exactness 6m + 2. In order to extend the quadrature formula concerning the integral on (b, c), we need to construct, if it exists, an optimal extension to formula (8). Rabinowitz [9], taking into consideration the Kronrod extensions to Gaussian integration rules for the evaluation of Cauchy principal value integrals of the form (l), introduces a new extension, for the case A = 0, which has the advantage of not involving the evaluation of the derivative of the integrand function. In fact he extends the 2m-points Gaussian rule adding 2m + 2 points, rather than 2m + 1 as in the Kronrod case; in this way he avoids the point x = h = 0. The extension of the formula (8) he constructs is of the type
with degree of exactness 6m + 3, where the nodes { JJ }F!+,’ are the zeros, whenever they exist, of the (2m + 2)th-degree Stieltjes polynomial i 2m+2(~) which satisfies the orthogonality conditions 1
(x)xk=O, J_lP,mbP2m+2
k=O,l,...,
04)
2m+l.
The coefficients of the polynomial ,!?2m+2(~) are determined, orthogonality conditions (14) which characterize it. The polynomial E *m+*(x)= X&m+lb) for any constant
up to a constant
by the
05)
- a
a is such that
1
(x)x” J-1 p2mb)E2m+2 Moreover
factor,
dx=O,
k=O,
l,...,
2m-1.
(16)
if we define
&n+z(x)= xE2m+lb)-
(17)
a”
where 1 a”=
/
h7(4L+lbb2m+1
Wj-;lf’d4x2n’
dx,
-1
this polynomial satisfies (14) and Monegato [2] proved that it is the unique one (up to a constant factor) satisfying it. We can construct (13) when the polynomial I&,,+2 (x) has 2m + 2 zeros which are real, distinct, located in the interior of (- 1, 1) and separated by the zeros of P2m(~); in this case for
A. Palamara Orsi / Cauchy principal value
the weights of (13) we can easily find the following
t=
tia,E2m+*(Sr)
-
t,E2m+1(5i)
CiH,
_~
,
i=
expressions:
1, 2,...,2m,
22”bn~,
s;= k,*,P2,(j,)[E,,+,(y,)
255
+y;EE;,+,(y,)] ’
I=1923.-.T2m+27
(18)
where the leading coefficient of E,,+,(x) is 22” and the leading coefficient of Pz,,,(x) is k2*, = [l . 3. 5 . . _. *(4m - 1)]/(2m)!. However the new quadrature rule (13) does not exist for any m; in the examined case of weight function equal to 1, we need to make m odd so that I? 2n,+2(~) can have 2m + 2 real and distinct zeros in (- 1, 1). Furthermore Rabinowitz [9] notices that, in this situation, these zeros exist, but, due of the decreasing of the smallest positive zero of E”2,,,+2(~) (which we denote with Jr) while m increases, the term of (13) T = (s,/Jl)[g( j+) - g( -y,)] can lead to the loss of significant figures. Nevertheless, from the numerical experiments we have performed, we notice that, if one has numerical cancellation in the quantity g( jl) - g( -jr), this is compensated by the fact that, in general, even if jj, is very small, the value of T is however negligible with regard to the value of the total sum and, therefore, the “erroneous” significant digits do not contribute to the formation of the final value. Anyway, in order to reduce the eventual numerical instability of (13) we have used a new polynomial of type (15) with m odd, in which the parameter a = a(m) is such that the polynomial still has 2m + 2 real and distinct zeros all contained in ( - 1, 1) (i.e., xTE2,,,+ 1(XT ) < a -c 0, where x,* 1s the smallest positive abscissa of the minimum points of xE2,,+,( x)), but the distance of his smallest positive zero from the origin is greater than J, (3, is the smallest positive zero of the polynomial E”2m+2(~) with the same m). Among the values of a for which these conditions are verified we have chosen the value ~7= 0.4x:E2,+,(x;)
-=Ia”,
because, based on numerical experiments done with m = 3(2)25, we have obtained the amplification factor K(a) with a variable, defined (in the case a = a”) as
the result that
(19) has a minimum when a = ii. Let us indicate this new polynomial ~2m+zb>
=
XE2m+lb)
by
(20)
-a
and its zeros by { j, } fJJc*. We can derive the new formula
(21) (where the symbols have a clear meaning) (16) is satisfied.
with degree of exactness
6m + 1, since in this case only
256
A. Palamara
Therefore f
we obtain
the following
y&lx =B(yf(A)
+
Orsi / Cauchy principal
global rule, in the case A > :( a + b):
‘c” lpf(xJ
+
yf
2y Ci’l’f(t,) r=l
i=l
+
ualue
B!2)f(Z,) +
2~2c;2)f(ui)
(22)
i=l
i=l
where u, = :(b - c)y, + :(b + c),
z, = :(b - c)& + :(b + c), B!2’ = 7,/& )
c,‘2’ = S,/jJl
and with the previous meaning of the other symbols. Formula (22) is constructed on 8m + 4 nodes and has degree of exactness We obtain an analogous formula in the case X < i( a + b).
6m + 1.
4. Numerical results For the values m = 3(2)25, we have determined the polynomial Elm+i( x) which satisfies (11) using the method described in [3]. In order to compute the above-mentioned polynomial, the polynomial P2m( X) and their first derivatives we have used Clenshaw’s algorithm applied to the corresponding finite expansions in Chebyshev polynomials. As initial approximations for computing the roots of E”2m+2(~) and E,,+,(x) by Newton’s method we have taken the midpoints between the zeros and the minimum points of E,,+,(x) evaluated previously. Tables 1 and 2 show the differences between formulas (13) and (21).
Table 1 Behaviour of the smallest positive zeros 9, and J, of the polynomials E”Zm+2(x) = x&,+,(x)- ii, respectively xEzm+ 1(x) m
d
a
Jl
d and E,,+,(x)
Yl
3
- 1.5D-2
4.2D-2
- lSD-1
5 7
- 5.7D-3 - 2.5D-3
2.OD-2 1.2D-2
- 9.5D-2 - 7.OD-2
1.4D-1 8.6D-2 6.3D-2
- i.4D-4
i.8D-3
- i.OD-2
i.8D-2
2;
Table 2 Amplification factors associated with formulas (8), (13), (21) N K
14 7.0
Z i?
62 9.9
78 10.4
102 10.9
10.1
30 8.5 15.6
20.3
22.4
25.1
6.5
8.0
9.5
10.0
10.5
=
257
A. Palamara Orsi / Cauchy principal value
In Table 1 we list the values of the constant a” which appears in the expression (17) of the Stieltjes polynomial &m+2(~) = xE~~+~(x) - a” whose zeros J, are nodes of the rule (13) and correspondingly the values of Z of the new polynomial (20) El,,,+ 2( x) = xEzm+ ,( x) - Z whose zeros y, are nodes of the r_ule (21). The latter polynomial E2,,,+2 is obtained by translation of the former ,!?Zm+2: acting on the constant a” we move the smallest positive zero away from the origin and we reduce the numerical cancellation in the corresponding quadrature rule. In Table 2, K, I? and K indicate the amplification factors associated with formulas (8) (13) and (21) respectively and N is the total number of nodes of the formula chosen. Then we have compared numerically the rule of type (9) with the corresponding one of type (22) (that is with m having the same value) applying them to integrals of type (1). The examples considered give us a measurement of the better approximation furnished by (22) with respect to that of (9) and show that using (22) we can obtain a reasonable estimate of the error in (9). In Table 3 we give abscissas and weights for a representative case (m = 3) of the new formulas (9) and (22). They are given for the interval (- 1, 1); because of symmetry only the positive abscissas and their corresponding weights are given. In Table 3 the symbols have the following meaning: 5, = positive nodes H, = weights of the (Y; = weights of the corresponding y, = positive zeros p, = weights of the corresponding r, = weights of the j, = positive zeros s, = weights of the
Table
of the 6-point Gauss-Legendre rule, 6-point Gauss-Legendre rule, 13-point Kronrod extension of the 6-point Gauss-Legendre to the abscissas Ei, of the 7th-degree Stieltjes polynomial ET(x), 13-point Kronrod extension of the 6-point Gauss-Legendre to the abscissa 0.0 and to the abscissas t,, 13-point new formula (21), corresponding to the abscissas of the 8th-degree new polynomial E,(x) = xE,(x) - ii, 13-point new formula (21), corresponding to the abscissas
3
Abscissas
and weights for the new formulas
(9) and (22), in the case m = 3
5,
H,
a,
c
0.2386191860831969
0.4679139345726911
0.2337708641169944
0.03756616183158088
0.6612093864662646
0.3607615730481386
0.1810719943231376
0.2092717906299459
0.9324695142031520
0.1713244923791704
0.08369444044690674
0.06260430386803677
Y,
P,
Y,
sz
0.0
0.2410725801734648
0.1375457357181660
0.2596716622719610
0.4631182124753046
0.2132096522719623
0.4264557065135209
0.2542968856122309
0.8213733408650279
0.1373206046344469
0.8362466005176528
0.1368677681526986
0.9887032026126789
0.03039615411981977
0.9831454575784491
0.03972142763354597
rule,
rule, <,, J,.
A. Palamara Orsi / Cauchy principal value
258
Table 4 Relative errors in the 13-point and 2%point rules ((9) and (22), respectively) for the integral H,(X) = f’ ,cos(5x +2)/(x - A) dx A
H,(X)
13 points
28 points
13 points with respect to 28 points
0.01 0.05 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 0.95 0.99
-
2.8D-5 1.7D-5 8.8D-6 2.3D-6 3.4D-7 5.1D-8 5.6D-9 1.7D-8 3.5D-7 5.7D-6 3.8D-5 2.7D-5 2.1D-5
3.2D-12 1.4D-12 4SD-13 4.5D-14 2.7D-15 4.5D-16 1.7D-16 1.9D-16 3.8D-16 4.6D-15 7.1D-14 8.1D-14 9.1D-14
2.8D-5 1.7D-5 8.8D-6 2.3D-6 3.4D-7 5.1D-8 5.6D-9 1.7D-8 3.5D-7 5.7D-6 3.8D-5 2.7D-5 2.1D-5
2.74829889791444 2.39910773684172 1.82773804017123 0.376974897586441 1.18225431509743 2.47184332651756 3.17940827525908 3.13424923411046 2.34722742476011 1.00041044109658 - 0.642524886674366 - 1.57953528558164 - 2.92679456816624
Table 5 Relative errors in the (4m + l)-point with m = 3, 5, 25, for the integral
and (8~ +4)-point H,(h)
= f’,dg/(x
rules ((9) and (22), respectively), - X) dx
A
H,(X)
13 points
28 points
13 points with respect to 28 points
21 points
44 points
21 points with respect to 44 points
101 points
204 points
101 points with respect to 204 points
0.01 0.05 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 0.95 0.99
- 0.031415926536 -0.157079632679 -0.314159265359 -0.628318530718 - 0.942477796077 - 1.256637061436 - 1.570796326794 - 1.884955592154 -2.199114857513 - 2.513274122872 - 2.827433388231 - 2.984513020910 - 3.110176727054
2.9D-2 8.7D-3 4.3D-3 1.9D-3 l.lD-3 6.9D-4 4.4D-4 2.6D-4 1.3D-4 2.3D-5 6.2D-5 5.4D-5 4.9D-3
6.4D-3 1.4D-3 6.7D-4 3.1D-4 1.9D-4 1.3D-4 8.8D-5 6.2D-5 4.3D-5 2.7D-5 1.2D-5 3.7D-6 1.4D-5
2.3D-2 7.4D-3 3.6D-3 1.6D-3 9.2D-4 5.7D-4 3.5D-4 2.OD-4 8.7D-5 3.6D-6 7.4D-5 5.OD-5 4.9D-3
9.7D-3 2.2D-3 l.OD-3 4.6D-4 2.7D-4 1.7D-4 l.lD-4 6.4D-5 3.2D-5 6.4D-6 1.7D-5 2.8D-5 5.7D-4
1.5D-3 3.OD-4 1.4D-4 6.6D-5 4.OD-5 2.7D-5 1.9D-5 1.3D-5 8.9D-6 5.4D-6 2.3D-6 5.3D-7 3.3D-6
8.3D-3 1.9D-3 8.9D-4 4.OD-4 2.3D-4 1.4D-4 8.7D-5 5.OD-5 2.3D-5 9.4D-6 1.9D-5 2.8D-5 5.7D-4
l.OD-4 2.OD-5 9.4D-6 4.2D-6 2.4D-6 lSD-6 9.6D-7 5.8D-7 3.OD-7 6.3D-8 1.5D-7 2.7D-7 4.1D-7
9.6D-6 1.9D-6 9.OD-7 4.1D-7 2.5D-7 1.6D-7 l.lD-7 7.5D-8 4.8D-8 2.6D-8 6.4D-9 4.8D-9 1.8D-8
9.2D-5 l.SD-5 8.5D-6 3.8D-6 2.2D-6 1.3D-6 8.5D-7 5.OD-7 2.5D-7 3.7D-8 1.6D-7 2.6D-7 3.9D-7
we have
obtained
In Tables 4 and 5 we report some numerical results above-mentioned rules to the following two integrals:
by applying
the
dx = --TX. In Tables 4 and 5 we report the relative accuracy obtained using the rule (9) with 4m + 1 nodes, with respect to the exact value, the relative accuracy obtained using the rule (22) with
A. Palamara Orsi / Cauchy principal value
259
8m + 4 nodes, with respect to the exact value and the relative accuracy obtained using the rule (9) with respect to the value obtained using the rule (22). We assume the last error as error estimate of integral approximation obtained using the rule (9); we can also assume it as an overestimation of the error of integral approximation obtained using the rule (22). In Table 4 we can see that the value of the integral supplied by (22) is already accurate to nearly machine precision with 28 nodes; this is difficult to realize if one only considers the estimation of the error which is reported in the last column. On the contrary, in Table 5 we can see that the extended formula is not much better than the basic formula. This different behaviour is due to the different degree of smoothness of the two integrands; indeed in the first case the integrand function is much more smooth and easier to be approximated by a polynomial than in the second case.
5. Adaptive routine For the evaluation of (1) for a single value of A E (a, b), supposing the function f(x) has a low degree of smoothness in the interval [a, b], it is convenient to apply an adaptive strategy analogous to that used by Piessens in [6] to evaluate the integral labf( x) dx and in [7] to evaluate (1). Indeed an adaptive strategy locates a few nodes in the s&intervals in which the integrand is regular and many nodes in the neighborhood of difficult spots of the integrand. If, however, we want to evaluate (1) for a set of values of A and for smooth functions f(x), then the use of Gauss-Kronrod rules for Cauchy principal value integrals will be more efficient, as Rabinowitz [8] pointed out. Let us describe the automatic globally adaptive integration scheme used to obtain an approximation y, for some preassigned relative tolerance tol, to the integral (1). Given the interval I = (a, b), the function f and the parameter A, let
G*,(I; f> K 4,+1(~; f) G&,+i(I; f) K&+2( 1; f) G4*m+i(I; f) K,*,*,,(I; f)
= 2m-point Gauss-Legendre rule applied to the function f( x)/( x - A) with x E [a, b] (with degree of exactness 4m - l), = K ronrod extension to the GZm( I; f) integration rule, = rule of type (3) where the nodes { ti } fr, coincide with those of the Zm-point Gauss-Legendre rule (with degree of exactness 4m), = Kronrod extension to the G&+i( I; f), = rule (9) (with degrees of exactness 4m), = rule (22).
Algorithm Step 1. Input: {a, b, A, f, m, to/}, where a < X -c b. Step 2. Let 1, = (a, b) Set ier = 0 (ier: error flag) Calculate est( I,) = K,*,*,,( 1,; f), Step 3.
err(h) = I G4*m*+lUl; f> - KiLL(4; If err( I,) < tel. 1e.st(Il) 1, then y = e.st(ll) and exit.
f> I
260
Step 4.
Step 5. Step 6.
Step 7.
Step 8.
Step 9. SteplO.
Step 11. Step 12.
Step 13. Step 14. Step 15.
Step 16.
A. Palamara Orsi / Cauchy principal value
Divide the interval 1, into two subintervals I1 and I*, eliminating the previous I, and adding a new interval I,, and set k = 1, n = 2. The subdivision of the initial interval Ii must be made as follows: if h @ I1 then Ii is bisected; if a < h < (a + b)/2 then Ii is subdivided at the point (A + b)/2; otherwise subdivision occurs at the point (A + a)/2. For I = k and I= n do Let ((Y, p) denote the open interval I,. If h E I,, then calculate est( I,) = &*m*++(I,; f), err (IA = I G4*m+1U~; f) - GZ+dC f> I and go to Step 9 of the algorithm. Else if X>/3 and IX--/31 -C Ip-c~l or if h < (Y and A G I, but it is near I,), then calculate est( II) = K4*m+2(I,; f), err(h) = I G’k+d4; f> - KEm+2U1; f> and go to Step 9 of the algorithm. Else (i.e. if A 6 I, and it is far from II) calculate est( I!) = K,, +1( I,; f),
I A - a I < I j3 - a I (i.e. if
I
err(h) = I %M; f> - &m+I(4; f) I. End of Z-loop. ) the initial interval (a, b) is divided into n subintervals Ij, At step n (n=2,3,... j = 1, 2,. .-, n, for each of which the estimates est( I,), err( I,) are available, where these values are carried over from the previous steps or have been computed for the subintervals just generated. If the maximum number of integrand evaluations allowed has been achieved, then set ier = 1 and go to Step 14 of the algorithm. Else if the length of some interval I, is getting too small, then set ier = 2 and go to Step 14 of the algorithm. Else if tel. IC’Jclest( I,) I < E (6: relative machine accuracy), then if C/“,,err(Ij) G tol, then go to Step 14 of the algorithm. Else if E3=,err(Ij) G tol- IC’J=,est(Ij) 1, then calculate y = CT= ,est( I,) and exit. Else split up the interval I, on which err( I,) = maxi <, i ,,err( I,) into two subintervals using the same strategy described at Step 4 of the algorithm; append a new interval I n+l to the list of subintervals and replace the old Ik by a new one. Set n = n + 1 and go back to Step 5 of the algorithm.
We remark that the strategy indicated at Step 4 of the algorithm is used to avoid a subdivision occurring at the abscissa A of the non-integrable singularity or too close to the abscissa X; in this latter case, in fact, some subtractive cancellation in the computation of the total integral estimate would happen.
A. Palamara Orsi / Cauchy principal value Table 6 Integral H,(h) = f’IeX/(x - A) dx; number of function evaluations required Quadpack at requested relative accuracy to1 and respective really achieved relative
x
0.01 0.05 0.1
0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 0.95 0.99
Present method
Quadpack
to1= 10-12
to1=10-‘2
261
by the present accuracies
NE,
AERR,
NE,
AERR 2
28 28 28 28 28 28 28 28 28 28 28 28 28
4.2D-16 4.3D-16 5.6D-16 2.4D-16 4.1D-16 3.4D-16 0.0 l.lD-15 2.4D-16 6.6D-16 5.8D-16 4.5D-16 5.OD-16
25 25 25 25 25 25 25 25 25 25 25 25 25
8.4D-16 8.6D-16 5.6D-16 2.4D-16 5.5D-16 1.7D-16 1.2D-16 1.9D-15 l.lD-15 9.2D-16 9.2D-16 1.5D-16 6.7D-16
method
and
by
Concerning Step 8 of the algorithm, we remark that when A is far from the integration interval, we do not use the formula G -&+l(I; f) to avoid the computation of A( A) which is much more expensive than the test regarding the nearness or not of A to the integration interval. In Tables 6-10 some numerical results obtained by applying the routine which implements the previous algorithm and the analogous routine QA WC-QA WCE of Piessens [7] to the following integrals are shown: H,(h)=fl~dx,
NV
= f,
H,(h) = H,(A)
-l
f_l,“‘,IT2 dx,
= f_l, “‘,‘_I;
X-h
We asked the routines
-l
’ dx,
1 1x - 0.65 I’.’
to integrate
-l
-l
O
each test function
to relative
accuracies
to1 = 10-3, 10-6, 10-9, lo-12. The results indicate that the two routines have about the same average reliability (if we define a program to be reliable when the result obtained achieves the requested tolerance), but the first one appears more efficient; in fact the number of integrand evaluations required to calculate a Cauchy principal value integral with a given accuracy is significantly inferior for functions f(x)
28
28
28
28
28
28
28
28
28
28
28
28
0.05
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
0.95
0.99
9.1Dp14
8.lD-14
7.5D-14
6.OD-15
0.0
0.0
4.2D-16
3.6D-16
153
1.2D-16
96 138 112 112
2.7D-15
2.4Dp15
8.4D-16
3.OD-16
70
70
70
70
28
28 70
112
4.2D-16 70
112
3.6D-16
0.0
111
2.1Dp15
137
154
1.9Dp16
1.3Dp15
153
NE,
0.0
AERR,
0.0
28
28
28
70
2.1D-15
70
4.4D-14
70
4.5D-13
1.4D-12
70
NE,
3.2D-12
AERR,
NE,
28
tol=lO~”
number
of function
1.5Dp16
7.OD-16
3.6Dp15
2.7D-15
0.0
1.4D-16
1.4Dp16
3.6D-16
0.0
1.3Dp15
0.0
1.9D-16
1.6D-16
AERR,
25 25
5.9Dp16 0.0 3.6Dp16
324 222
25
7.OD-16
190
25
25
3.5D-15
190 1.5Dp16
25
2.2Dp15
164 190
25
25
2.8D-16 0.0
180
25
205
205
2.8Dp16
25
0.0
299
298
25
25
3.7D-16
289
25
3.2D-16
0.0
5.6D-16
4.OD-15
5.5D-15
1.7D-15
1.8D-15
7.OD-16
1.6D-15
l.lD-15
lSD-15
2.1Dp15
7.4Dp16
1.5D-15
method
25
65
65
65
65
65
65
65
65
65
65
65
65
NE,
0.0
3.1Dp8
7.4Dp8
4.3Dp8
1.5D-8
8.6D-9
5.4D-9
3.OD-9
1.6Dp9
2.6D-8
8.9D-9
7.7Dp9
7.3D-9
AERR,
tol= lo-”
by the present
AERR,
NE,
AERR,
NE, 247
Quadpack
required
to/=10-’
evaluations
t0l=10~”
relative accuracies
tol=lO~’
+ 2)/( x - X)dx;
really achieved
tol= 10-3
Present method
0.01
x
rol and respective
= fL ,cos(5x
H,(X)
Integral
accuracy
Table 7
375
435
555
405
315
325
285
445
595
835
445
415
385
NE,
2.7D-10
1.3Dp10
3.OD-10
1.8D-11
1.6D-11
1.2D-9
4.5Dp10
1.9D-11
4.9Dp10
7.5D-10
6.1D-10
7.1Dp10
l.lD-9
AERR,
to/=lo-’
and by Quadpack
6565
8245
11695
11795
7475
6165
6225
8925
15055
15045
12745
10825
9325
NE,
tol=lo~‘*
relative
3.8D-12
1.9Dp12
4.1Dp12
4.4D-12
1.5D-12
3.OD-12
7.9D-13
4.ODp13
9.ODp13
4.2D-12
1.9D-12
3.7D-12
4.4D-12
AERR,
at requested
%-
x
3. R G’ k
:
9 *’ \ 0 !S E
1.9Dp4
1.3Dp4
8.8D-5
6.2D-5
4.3D-5
2.7D-5
1.2Dp5
3.7D-6 2.3Dp6
28
28
28
28
28
28
28 112
0.4
0.5
0.6
0.7
0.8
0.9
0.95 0.99
2.5Dp4
70
28
0.2
0.3
1.7Dp8
9.5D-9
1.4Dp8 l.lD-8
513
555 613
1.9Dp8
5.5Dp9
497
455
439
1.6D-8
9.7Dp9
439
1.6D-8
6.3D-9
449
423
440
2.2D-9
475
385 395 405
1429 1606
5.ODp12 l.lDp11
893 951
8.5D-15 5.ODp15
345
5.3Dp16 4.lD-15
1345
335
1387
295
9.1D-15
1.2D-11
1370
1.3D-11
819 6.6Dp12
1354
1.3D-11
860
325 195
4.8Dp15 2.1Dp15
7.6D-15
315
315
375
405
495
877
1380
l.lD-11 1.3D-11
846 819
2.4D-15
8.8Dp15
7.1Dp15
l.lD-15
1.3D-14
861
1417 1416 1407
4.3Dp12 4.4Dp12
829
1.5D-12
881 896
1641 1562
1.9D-11
1.3Dp5
1771
1.8Dp13
908
70
7.8Dp10
985
0.1
528
1.8D-8
l.lD-6
605
2.6D-5
164
X15 865
l.lD-4
775
765
725
685
715
715
705
735
765
825
915
NE,
1.3Dp7 3.7D37
2.7D-7
3.8D-X
8.6Dp8
1.5D--3
9.2Dp8
1.4D-7
2.4D-7
1.9Dp7
1.2Dp8
2.7Dp7
method
1.3Dp9 1.4Dp9 1.4D-9
1515 1565
4.4Dp9 3.XDp9
l.OD-9
l.ODp9
6.7Dp10
1.3Dp9
5.4D-10
1505
1465
1425
1455
1355
4.8Dp9
5.5D-9
4.5D39
2.8D-9
4.OD-9
1545
4.ODp11 5.2Dp9
4.3D-11 1445
1685 1535
4.3D-11
1925
1.3Dp10
1.5Dp10 3.8D-11
2735
AERR,
y
NE,
toI = 10
and by Quadpack
9.7Dp11
9.4D-11
X.ODp11
l.lDp10
AERR,
rol = lo-”
by the present
AERR,
241
NE,
0.05
AERR,
0.01
NE,
AERR,
NE,
AERR,
NE, NE,
tol=lO~’
rol=lo~‘*
tot-10-h
required
rol=lO~~ tol=lO~’
evaluations
Quadpack
AERR,
of function
relative accuracies
- X)d:. 1. number
really achieved
Present method
x
to/ and respective
= f’,js/(~
H,(X)
Integral
accuracy
Table 8
4.2Dp12 4.4Dp12 4.9Dp12
11515 11645 11725
4.ODp12
10895 11535
1.3Dp12 3.4D-12
11885
l.lDp12
14465
2.7D-12
13615 10915
1.5D-14 9.7Dp14
15055
3.5Dp13 2.7D-13
15055
3.7Dp13
15055 15055
AERR
lo-‘*
relative
NE,
tol=
at requested
898
1.4D-4
2.8D-16
3.6D-3
2.7D-5
278
195
153
112
28
138
138
138
28
70
70
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
0.95
0.99
1.4D-5
0.0
3.8Dp4
1.9D-5
l.lD-16
1.7D-13
7.7D-5
70
242
216
242
242
138
28
242
366
195
450
0.0
2.213-8
8.4D-7
4.1D-7
7.8D-8
l.lDp16
1.7Dp13
3.1D-7
1.3D-8
2.8D-16
8.6Dp8
70
0.0
9.1D-11
7.7Dp10
346 372
1.4Dp9
398
164
2.4Dp13 0.0
502 96
2.4D-15
6.7Dp13
1.3Dp14
l.lD-16
2.1D-16
6.7Dp14
5.3Dp14
4.2Dp16
8.1D-14
1.6DpI 3
1.2D-13
AERR,
evaluations
580
606
580
l.lD-16 6.1D-11
96
528
782
273
840
2.1Dp16
7.6D-11
6.5D-11
2.8Dp16
424
138
70
398
548
195
606
8.4D-11
1194
1.3Dp11
7.2D-8
2.2D-8
716
5.4Dp12
534
700
908
l.lDp4
320
6.5D-4
570
0.01
NE,
AERR,
NE,
AERR,
NE,
NE,
AERR,
toi=10-’
to1 =10p6
rol=lO~”
1x1)/(x-A)d x; number of function really achieved relative accuracies
rol=lO-’
Present method
0.05
x
Table 9 Integral H3(X)=fl,(x4+ accuracy to/ and respective
65
65
65
225
225
135
65
195
285
65
335
445
605
NE,
2.6D-6
9.9Dp6
1.5D-4
4.5Dp6
l.ODp7
8.9D-9
l.lD-3
8.9Dp8
4.1Dp6
l.lDp2
6.7D-6
1.3D-6
1.4D-6
method
225
315
345
405
345
135
345
345
435
465
485
595
695
NE,
4.ODp6
2.2Dp10
9.OD-9
5.1 D-8
8.6Dp9
8.9Dp9
4.6Dp9
4.50-9
2.OD-9
2.8D-9
l.OD-8
5.2Dp9
6.2Dp8
AERR,
tol= 10-h
by the present
AERR,
tot =10-s
Quadpack
required
495
555
555
765
615
195
555
555
915
885
1025
1345
1985
NE,
4.ODp6
7.5D-IO
6.3D-10
7.7Dp10
9.7Dp10
1.2D-9
5.3D-10
1.6D-9
1.7D-10
8.2Dp10
7.2Dp10
7.8D-10
9.6D-10
AERR,
to/ = lomq
and by Quadpack
1905
2445
3045
5475
3195
2325
2565
2415
7335
7245
10175
14425
15125
NE,
relative
4.ODp6
4.9Dp12
5.9D-12
7.OD-12
5.2Dp12
4.8Dp12
4.3Dp12
3.1D-12
1.6Dp12
2.ODp12
2.5D-12
4.2Dp12
l.OD-11
AERR,
to1=10-‘2
at requested
+
ti iii-
t
$.
2 3.
\ 0 e %
i? a’ 0 a __
2 F
6.3D-4
5.1Dp4
3.1D-4
8.3Dp5
1.2D-3 2.4D-7
2.ODp16 5.3D-5
l.lD-6
2.2Dp4
l.OD-7
4.7D-4
69
69
28
112
69 277
2X 278
221
112
138
70
0.2
0.3
0.4
0.5 0.6
0.65 0.7
0.X
0.9
0.95
0.99
190
216
242
429
2X 408
256 319
257
121
147
199
450 39X
6.3D-7
4.ODp8
466
808
3.2D-9
5.9D-7
28 725
2.OD-16 1.6D-7
413
2.4Dp7 495 620
448
1.6Dp7 l.OD-8
449
6.4Dp9
422
449
4.5D-9
1.7D-7
6.9D-7
277
0.1
174
3.7Dp7
l.lDp7
96
0.05
173
1.5D-5
69
0.01
rol = lo-
‘*
l.ODp9
6.7Dp12
2.7Dp10
9.5D-12
1.5D-15 l.OD-13 1.3Dp12 1.5Dp13
846 700 752
2.2Dp13
0.0
1593
96
7.4D-14 6.1D-15
X90 1000 1354
l.lD-13
808
2.ODp16
9.6D-16 4.1D-15
2.1Dp13
734 843
1.4Dp12
880
5.XDp15
605
AERR,
6X3
NE,
2.6D-11
2.5D-11 X.4Dp12
1.6D-10
6.ODp12
X.2Dp12
4.5D-10
7.3D-12
3.3D-11
AERR,
195 25
225
355 215
385
145
65 235
155 65
185
155
25
NE,
7.1Dp7 4.5Dp4
2.2D-7
2.3Dp7
4.6Dp8
4.OD-3 9.9D-4
1.6Dp6
1.15Dp3
255
285
345
615
545
1145
365
395
355
285
255
275
245
155
NE,
to/-l0
1.2D-8
2.3D-8
1.6D-8
7.3D-7
2.1D-8
6.1Dp7
1.5D-6
4.6Dp10
6.2Dp9
1.4Dp9
1.3Dp8
2.7Dp8
2.1Dp8
1.3Dp6
AERR,
oh
by the present method
4.ODp6
1.9D-9
2.1D-6
7.7Dp7
l.lDp4
AERR,
tol=lO_’
ro/= lo-” NE,
AERR,
AERR,
ro/=lO_”
NE,
to/=lOK’
Quadpack
required
Present method
NE,
A
x. number of function evaluations
really achieved
relative accuracies
IO/ and respective
accuracy
10
Integral H4(X)=fbI_~-0.651’-‘/(x-X)d
Table
10
7.XD-10 5.3D-10
555
6.ODp 10
3.4Dp9
2.OD-9
6.6D--10
7.6Dp10
7.7D-
685
745
4505
1405
19x5
915
735
1.3Dp9
5.6Dp10 725
5.2Dp 10 715
X.3DpIO
l.ODp9
2.XD- 11
675
5X5
5X5
2.2Dp10 1.3D-12 1.9Dp12 2.lDp13
15075 9885 5585 2995
l.lDp11
l.OD-12 15105
7735
2.XDp12 4.3Dp12
76X5 9035
1.9D-12 3.7Dp12
8525 7825
2.ODp12
3435 7755
3.3Dp12 2.6Dp12
2565
3.1D-12
AERR,
relative
1935
NE,
425
to/=10-‘2
AERR;
NE,
at requested
tol= 10-9
and by Quadpack
A. Palamara Orsi / Cauchy principal value
266
having a low degree of smoothness and for high accuracies. Only for well-behaved integrands and low accuracies the two routines, generally, are competitive. In particular, we can notice the advantage of the former routine over the latter, also for low accuracies, in the case of the integral H4( X) with X = 0.65 in which the integrand is not singular (see Table 10). Moreover, for our test integrals, with the present method the actual error in the computed values is almost always at least one order of magnitude smaller than the requested error; only in a few cases (see Table 9: toZ= 10e3, h = 0.9; tol= 10p9, h = 0.8. Table 10: tol= 10p3, X = 0.5; to1 = 10-9, A = 0.99; to1 = lo-‘*, X = 0.05 and X = 0.95) it is a little higher (by a factor less or equal to 3.6) than the requested tolerance. Also with Quadpack quite often the actual error is smaller than the requested error, though the anomalous cases are more than a few. In Tables 6-10, NE, and NE, denote the total number of function evaluations while AERR, and AERR, denote the actual relative error in the computed values of the integrals, using our routine with m = 3 and QA WC-QA WCE routine respectively. All computations reported on this paper were carried out in double precision on the PDP-11/23+ computer which has about 16-digit accuracy in double precision.
Acknowledgment I gratefully
acknowledge
G. Monegato
for his constant
guidance
during
this work.
References [l] G. Monegato, A note on extended gaussian quadrature rules, Math. Cump. 30 (1976) 812-817. [2] G. Monegato, On polynomials orthogonal with respect to particular variable-signed weight functions, ZAMP 31 (1980) 549-555. [3] G. Monegato, Stieltjes polynomials and related quadrature rules, SIAM Rev 24 (1982) 137-158. [4] G. Monegato, The numerical evaluation of one-dimensional Cauchy principal value integrals, Computing 29 (1982) 337-354. [5] D.F. Paget and D. Elliott, An algorithm for the numerical evaluation of certain Cauchy principal value integrals, Numer. Math. 19 (1972) 373-385. [6] R. Piessens, An algorithm for automatic integration, Angew. Inform. 9 (1973) 399-401. [7] R. Piessens, E. de Doncker-Kapenga, C.W. ijberhuber and D.K. Kahaner, QUADPACK: a Subroutine Package for Automatic Integration (Springer, Berlin, 1983). [8] P. Rabinowitz, A stable Gauss-Kronrod algorithm for Cauchy principal-value integrals, Comp. & Maths. Appl. 12 B (5/6) (1986) 1249-1254. [9] P. Rabinowitz, Gauss-Kronrod integration rules for Cauchy principal value integrals, Math. Comp. 41 (1983) 63-78. [lo] I.H. Sloan, The numerical evaluation of principal value integrals, J. Comput. Phys. 3 (1968) 332-333.