Документ взят из кэша поисковой машины. Адрес оригинального документа : http://num-meth.srcc.msu.ru/english/zhurnal/tom_2009/pdf/v10r145.pdf Дата изменения: Wed Nov 25 16:42:35 2009 Дата индексирования: Mon Oct 1 23:28:00 2012 Кодировка:

Numerical Methods and Programming, 2009, Vol. 10 (http://num-meth.srcc.msu.ru)

391

PRINCIPLES OF THE BAYESIAN AUTOMATIC ADAPTIVE QUADRATURE Gh. Adam , S. Adam
1 1

The addition of Bayesian inference to the automatic adaptive quadrature scheme developed in QUADPACK (R. Piessens, E. deDoncker-Kapenga, C.W. Uberhuber, and D.K. Kahaner, QUADPACK, a quadrature package for automatic integration, Springer, Berlin, 1983) signi cantly increases the output reliability under the lack of a priori knowledge on the behavior of the function over the integration domain. In the present paper, we review the progress obtained until now along these lines.

Keywords: numerical integration, automatic adaptive quadrature, integrand pro le analysis, reliability, Bayesian inference. 1. Problem statement. The evaluation of (proper or improper) one-dimensional Riemann integrals by automatic adaptive quadrature 1{3] has been implemented and is available in the most authoritative program libraries (e.g., SLATEC, IMSL, and NAG). To alleviate severe code failures in cases of practical interest 4{6], separate codes, each being able to solve a given speci c class of integrals, have been implemented 7, 8]. It is then user's responsibility to choose an appropriate code from a library for the problem of interest. Such an approach remains, however, useless in the case of parametric integrals arising in various physical models (see, e.g., 9{13]). Since the variation of the model parameters results in the occurrence of integrals falling in di erent classes, a priori decisions concerning the assignment of the right code cannot be taken. We are thus left with the trial and error approach, with an unacceptably high rate of failure and frustration. If the reliability of the local quadrature rule output pairs (q; e) is explicitly questioned via the use of post validation consistency criteria 14, 15], then the overwhelming fraction of the spurious (q; e) pairs is ruled out, with the consequence that the class conscious decisions 2] of the automatic adaptive quadrature algorithms get signi cantly improved. In the present paper, we discuss a set of necessary consistency criteria allowing the identi cation of the spurious (q; e) outputs prior to the activation of the local quadrature rules. This is an instance of Bayesian inference 16] which, by elimination of the guaranteed spurious outputs (q; e), increases the chance of obtaining meaningful (q; e) pairs under complete lack of a priori know ledge on the integrand function. 2. De nitions and notations. 2.1. The integral. We consider the (proper or improper) one-dimensional Riemann integral
I I a; b]f =
Zb

where the integrand function f : a; b] ! R is assumed to be continuous almost everywhere on a; b] such that (1) exists and is nite. If the integrand is factorized as a product g(x)f (x), where the weight function g(x) absorbs an analytically integrable di cult part of the integrand (e.g., an endpoint singular or oscillatory function), then the following considerations are equally valid for this integrand function f (x). 2.2. Local quadrature rules. Given ; ] a; b], a local quadrature rule produces an approximate solution of I ; ]f as a couple fq; eg, where q Q ; ]f denotes a quadrature sum approximation of I ; ]f , while e > 0 denotes a probabilistic bound of the error (estimate of the error) associated to q. If e > jeq j, where eq = I ; ]f q is the actual error associated to q, then the couple fq; eg is reliable, otherwise it is unreliable. In the rst case, the decisions of a class conscious automatic adaptive algorithm are meaningful, while in the second a wrong decision branch may be chosen and the numerical solution fails. 2n X We assume that q is a (2n + 1)-knot interpolatory quadrature sum q q2n Q2n ; ]f = wi f (xi ) at Laboratory of Information Technologies, Joint Institute for Nuclear Research, 141980, Dubna, Russia and Horia Hulubei National Institute for Physics and Nuclear Engineering (IFIN-HH), 407 Atomistilor, Magurele, Bucharest, 077125, Romania; e-mail: adamg@jinr.ru, e-mail: adams@jinr.ru c Research Computing Center, Moscow State University
1

a

f (x) dx;

(1)

i

=0


392

Numerical Methods and Programming, 2009, Vol. 10 (http://num-meth.srcc.msu.ru)

the 2n + 1 abscissas (called the quadrature knots ) inside ; ]: 6 x0 < x1 < : : : < x2n 6 : (2) An automatic adaptive quadrature algorithm may operate with one or more local quadrature rules. In what follows, we assume that these may be either Gauss{Kronrod (GK) or Clenshaw{Curtis (CC) quadrature rules 1{3]. Both of them result in symmetric quadrature sums, with the interpolation abscissas (2) given by xi = c + hyi ; 1 c = 2 ( + ); h = 1 ( ) ; i = 0; 1; : : : ; 2n, where 2

0 = yn < yn+1 < yn+2 < : : : < y2n 6 1; yn i = yn+i ; i = 1; : : : ; n; (3) denote the reduced quadrature knots yi de ned over 1; 1]. 2.3. Standard automatic adaptive quadrature. The automatic adaptive quadrature was systematically developed in QUADPACK 1], the de facto standard of one-dimensional globally adaptive numerical integration. A globally adaptive quadrature algorithm involves the following two fundamental steps (QUADPACK 1], p. 60). (i) Initialization: the number of integration subranges is set to N = 1 and the local quadrature rule (q; e) is used to solve the given integral over the whole integration domain to yield the initial global pair (Q1 = q, E1 = e > 0). Using Q1 and the input accuracy speci cations, an initial estimate 1 of the acceptable tolerance associated to the initial solution is computed. (ii) Error decrease by subrange subdivision: if EN > N (i.e., the global error estimate EN > 0 associated to the composite quadrature sum approximation QN exceeds the tolerance level N ), then the local quadrature errors are decreased by subrange bisection (hence, N is increased to N + 1) and the global quantities QN ; EN ; N are updated until the error tolerance level is achieved (EN < N ). 2.4. The integrand pro le. The set of all currently computed values of the integrand over the current subrange ; ] a; b] de nes the integrand pro le over ; ]. Within the Bayesian automatic adaptive quadrature, the integrand pro le comes from two kinds of abscissas: inherited from the ancestor subranges and the local quadrature knots (3) inside ( ; ). While not currently needed for the computation of the local (q; e) pair, the inherited abscissas and integrand values provide valuable enhancement of the quality of the integrand behavior analysis over ; ]. The role of the local quadrature knots (2) or (3) is twofold. First, they serve to the derivation of hints on the integrand conditioning over ; ]. Second, under ful llment of all the well-conditioning criteria by the integrand function, they serve to the derivation of the outputs (q; e) for the local quadrature rules. The union of the inherited and currently computed abscissas over ; ] de nes the ne discretization set of abscissas over ; ]. This may be de ned either in terms of xl < xl+1 < : : : < xr ; the absolute abscissa values, or in terms of yl < yl+1 < : : : < yr ; the reduced abscissa values. Both sets of abscissas are uniquely de ned for any arbitrary subrange ; ] a; b]. In what follows, the distance between two reduced abscissas will be of interest: (4) j k = yk yj ; j; k 2 l; l + 1; : : : ; rg: behavior around a certain quadrature knot x terms of the already available integrand pro neighborhoods. The left ne discretization neighborhood Fl (xk ) = xk The left coarse discretization neighborhood of

2.5. Discrete neighborhoods over the integrand pro le sampling. The analysis of the integrand k 2 ; ] asks for the use of neighborhoods which are de ned in
le sampling. This analysis is local and involves several kinds of (5) of xk is given by 2 ; xk 1; xk ; xk +1 \ xl ; : : : ; xr : xk is given by

lateral neighborhood is complete provided it contains exactly four points. Similar de nition holds for the right lateral neighborhood Lr (xk ).

Cl (xk ) = xk 2; xk \ xl ; : : : ; xr : (6) Similar de nitions hold for the right discretization neighborhoods Fr (xk ) and Cr (xk ), respectively. The left lateral neighborhood of xk is given by Ll (xk ) = xk 3; xk 2; xk 1; xk \ xl ; : : : ; xr . The left


Numerical Methods and Programming, 2009, Vol. 10 (http://num-meth.srcc.msu.ru)

393

The inner extremal point xk of the integrand pro le is isolated to the left provided Ll (xk ) is complete and the sequence fk 3 ; fk 2; fk 1; fk is monotonic. An inner extremal point xk is isolated provided it is isolated both to the left and to the right. A monotonic subset of sequential integrand pro le values fl ; fl+1 ; : : : ; fl+q de nes a monotonicity subrange inside ; ] provided its length q + 1 > 6. 3. Integrand features sub ject to Bayesian inferences. There are integrand features which result in conspicuously unreliable local quadrature rule (q; e) outputs if properly questioned and identi ed: | severe precision loss due to cancellation by subtraction; | occurrence of a range of variation of a monotonic integrand which exceeds the worst case bound inferred from the polynomial set spanning the interpolatory quadrature sum; | occurrence of a nite jump with nite lateral derivatives, immersed into a monotonicity subrange of the integrand; | same as previous, but turning point; | integrand oscillations at a rate of variation beyond the current quadrature knot set resolving power; | isolated irregular integrand extremum. The speci c decisions following from the identi cation of one or another of the above-mentioned cases depends on the diagnostic. One of the following three continuations is possible: (i) stop immediately the computation and return the appropriate error ag (when there is no hope to improve the output for the present problem formulation); (ii) proceed immediately to symmetric subrange bisection (when it is expected that the re nement of the discretization into subranges will result into a better resolved integrand pro le); (iii) proceed immediately to the solution of a number of auxiliary problems. If the occurrence of an inner isolated o ending point xs was inferred, then resolve its location inside ( ; ) to machine accuracy. Proceed then to the splitting ; ] = ( ; xs) (xs ; ) : The abscissa xs will be locked from now on at subrange boundaries within the subrange subdivision process of a; b]. If xs = a or xs = b, then solve one lateral boundary layer problem 17, 18] at x+ or xs , respectively, in order to determine the nature of s the integrand behavior at xs as well as appropriate integrand lateral limits. If xs 2 (a; b), then solve two lateral boundary layer problems at xs and x+ , respectively. s The solutions of the auxiliary problems de ne the further continuation of the algorithm. If xs is an essential singular point (i.e., it associates a singularity of f (xs ) together with in nitely many oscillations of f (x) in its neighborhood (like, e.g., sin (1=x) at x = 0+ )), then further continuation is useless. The computation is stopped immediately and the appropriate error ag is returned. If xs corresponds either to a nite jump or a turning point with nite lateral derivative, then its contribution to the original Riemann integral is nil. The local quadrature outputs (q; e) become insensitive to the occurrence of the nearby o ending locked endpoint. If there is a lateral singularity at xs in the integrand and/or its rst order derivative, then the local quadrature outputs (q; e) remain sensitive to the occurrence of the nearby o ending isolated singularity. Moreover, slow convergence under further symmetric subrange bisection occurs. However, convergence acceleration is possible by the use of extrapolation techniques. Therefore, a ag explicitly pointing to the allowance of the activation of a convergence acceleration procedure is set. This discussion points to the need of three pointers for the integrand behavior characterization over a subrange, corresponding to the left end, the right end, and the subrange interior respectively. The pointer ipinn characterizing the interior of a subrange carries the output analysis information given in Table 1. The pointer ipend characterizing a subrange end carries the information given in Table 2. 4. Order of integrand computation at quadrature knots over subranges. Within the standard automatic adaptive quadrature, all the integrand values asked by the local quadrature rule pair (q; e) are computed, irrespective of the meaningfulness of the (q; e) output or not. Within the Bayesian automatic adaptive quadrature, the computation of the integrand values and the analysis of the integrand behavior, prior to the activation of the local quadrature rule, are done in separate distinct procedures. A noticeable decrease of the number of integrand evaluations may be obtained provided the


394

Numerical Methods and Programming, 2009, Vol. 10 (http://num-meth.srcc.msu.ru)

Table 1 ipinn 0 -2 -1 2 ipend 0 1 2 3 All consistency criteria passed. Local quadrature rule activation allowed. Insu ciently resolved integrand pro le. Do immediate symmetric bisection. Localized integrand di culty. Solve auxiliary and/or boundary layer problems. Heavy cancellation by subtraction. Stop all computations and issue error ag. Table 2 Free, presumably regular, subrange end. Locked singularity in function and/or rst order derivative. Allow extrapolation. Locked zero measure subrange end. Extrapolation denied. Essential singularity. Stop all computations and issue error ag.

analysis is done as soon as possible after the computation of a new integrand value. The computation/analysis process may be optimized as a two-step interlacing involving the computation of an appropriate subset of integrand values, followed by the check of corresponding consistency criteria. The minimal subset of integral values requested by the analysis comprises the subrange ends and the inner local quadrature knots inside ( ; ) (a; b). In order to make the following discussion independent of the open (GK) or closed (CC) character of the local quadrature sum of interest, we consider the union of the quadrature knots (3) with the reduced endpoint abscissas: y < y +1 < y2n+ = y0 < y1 < : : : < y2n f 1; 1g. The integrand values f = f ( ) and f = f ( ) at the endpoints y = 1 and y2n+ = +1 are either inherited from the parent subrange or computed during the root initialization of the binary subrange tree. A consequence of the use of local quadrature sums of interpolatory type is the characteristic symmetric distribution of the quadrature knots inside every subrange ( ; ) (a; b): sparser towards the subrange center (such that the norm of distribution (4) is given by max j j kj = yn+1 = yn 1 ) and denser towards the subrange ends (with the outermost two quadrature knots lying signi cantly nearer to each other and to the corresponding subrange end as compared to the remaining ones). The average inter-knot distance
l l r l r

= (2n + r

l

)

1

(7)

provides a convenient threshold for the separation from each other of the sparse and dense knot regions respectively inside the left and right subrange halves. The center yn = 0 itself (or, in absolute units, = ( + )=2) plays a special role within the symmetric bisection since the pairs (f ; f ) and (f ; f ) provide endpoint inheritance for the descendent subranges. This discussion points to a computation/analysis two stage process which is to be done for seven distinct groups of inner quadrature knots: the center of ; ]; the left and right pairs of quadrature knots lying nearest to the endpoints and respectively inside ; ]; the left-half and right-half subsets covering the remaining dense quadrature knot distributions; ibid., for sparse quadrature knot distributions. Under the inheritance of previously computed integrand values over the ancestor ranges, the analysis process is completed by the operation of merging the inherited and just computed sequences with the purpose of enhancing the reliability of the Bayesian inferences. 5. Severe precision loss due to cancellation by subtraction. Following the idea rst developed in QUADPACK (1], p. 71), we formulate the following criterion which ends quickly the computation of a vanishing integral value under non-vanishing integrand. Criterion C1. If the integrand analysis at the initialization step of the automatic adaptive quadrature 2n+ 2n+ X X returns the result f (xk ) < 100 "0 f (xk ) , where "0 is the epsilon with respect to addition, then a
r r

roundo error ag is set and the computation is stopped.
l l

k

=

k

=


Numerical Methods and Programming, 2009, Vol. 10 (http://num-meth.srcc.msu.ru)

395

derivation of such a bound is Fact 1. Let XM = max subrange ; ] a; b]. The oating point degree of precision equal to d is an e K=
8 >d > >h < >h > > :

6. Upper bound to the global range of variation of a monotonic integrand. The key to the
provided by j j; j j denote the endpoint maximum absolute value of the current quadrature
of precision, K of an interpolatory quadrature sum showing an algebraic degree

ective value de ned as follows:
i i

ln "0 = ln XM

i i

XM 2 xm ; xM ]; xm = "1=d ; xM = xm1, 0 XM < xm , XM > xM ,
d X k
=0

ln "0 = ln XM i where a] denotes the integer part of a > 0.

Proof. The result follows from the consideration of the fundamental polynomial pd (x) =

xK and

derivation of those conditions under which al l the monomials of the fundamental power series 1; x; : : : ; xK bring signi cant bits to the computed values of pd (x) at various arguments x. From pd (x) it follows that 2K is a worst case upper bound for the range of variation of the integrand values over the triplet ff ; f ; f g. We may therefore formulate Criterion C2. Let ; ] a; b], f = f ( ), f = f ( ), f = f ( ), = ( + )=2; = jf f j; = jf f j. If f ; f ; f de ne a monotonic sequence, then the infringement of the condition 2K min ; > max ; points to the need of exiting the computation/analysis process and to proceed to the immediate symmetric bisection of ; ]. (i) ; ] a; b] and a sampling f0 ; f1 ; f2 around the endpoint x0 (x0 = or x0 = ) of ; ] over the abscissas set x0; x1; x2 , where x1 and x2 denote the two abscissas lying nearest to x0 within the merged set of currently generated and inherited abscissas, 0 e (ii) the three estimates f00(k) (k = 1; 2; 3) of f0 = f 0 (x0) following from the samplings S1 = f0 ; f1 and f e S2 = f0 ; f2 as a rst order divided di erence f00(k) = dk0 = xk f0 = fk f0 , k = 1; 2, or from the h k0 n k x0 o e e e sampling S3 = f0 ; f1; f2 as f00(3) = d10 + (d10 d20) 10;20, 10;20 = 10 , then the set f00(1); fe00(2); f00(3) is 20 taken for being consistent provided
e fe00(1) f00
(3)

7. Inferences from integrand slope approximation at a subrange endpoint. Criterion C3 (the endpoint slope consistency criterion ). Given:

6 fe00(2) fe00(3) 6 max f00(1) ; f00(2) ; f00(3) ; 1 ;

(8)

where is an empirical value set to = 1=3. If the consistency criterion (8) is infringed, then: (a) under jhj = j 2 j > 1 or a nonmonotonic sequence S3 , immediate symmetric bisection of ; ] is recommended, since it is expected that an insu ciently resolved integrand pro le by the set (3) of the local quadrature knots will result over ; ]; (b) the solution of a boundary layer problem at x0 inside ; ] is asked otherwise. 8. Check of Nyquist threshold for oscillatory integrands. The reconstruction of periodic signals (17], Chap. 12) shows that the structural details recovered by analysis cannot be ner than the norm of the discretization sampling. The Nyquist theorem established in this context has two straightforward implications in the Bayesian automatic adaptive quadrature with respect to the faithful representation of the integrand function structure by the integrand pro le at the set of the local quadrature knots (3). Criterion C4 (Nyquist local ). The faithful representation of a non-monotonic integrand variation by the pro le derived at the local abscissa set (3) asks for a lower bound of the distance between two successive extrema not smaller than 3=4 , with given by (7) and 3=4 being an empirical factor. Criterion C5 (Nyquist global ). The integrand pro le derived at the local abscissa set (3) is faithful if the number of counted oscillations inside it does not exceed the Nyquist threshold 2=yn+1.


396

Numerical Methods and Programming, 2009, Vol. 10 (http://num-meth.srcc.msu.ru)

The infringement of any of these two criteria points to an insu ciently resolved integrand pro le and this asks for immediate symmetric subrange bisection. 9. Integrand behavior at isolated inner extrema. If x0 is an extremal point of a smooth function f (x), then f 0 (x0) = 0 and there is a nite neighborhood V (x0) of x0 inside which f 00(x0 ) 6= 0. These fundamental properties result into the following Bayesian hint: Criterion C6. The integrand f (x) is smooth at the isolated extremal knot xk provided 0 0 (I) fF (xk ) < fC (xk ) ; where fF ; C g is a symbolic notation for the pairs Fl (xk ); Cl (xk ) and respectively Fr (xk ); Cr (xk ) , Eqs. (5) and (6); (II) the curvature of f (x) at xk 1; xk ; xk+1 keeps constant sign irrespective of the manifold Ll (xk ), Fl (xk ), Fr (xk ), or Lr (xk ) over which it is computed from third degree interpolatory polynomials. The infringement of any of conditions (I) and (II) points to an irregular extremum xk . Scale invariance of the diagnostic of an irregular extremum is checked under symmetric bisection until it is detected at a half-width jhj < 1, therefrom the analysis follows the general pattern described at point (iii) of Section 3. 10. Checking integrand regularity over monotonicity subranges. The analysis done in this case is k based 15] on the comparison of the rst order divided di erences dk+1;k = fh+1 fk and is formalized in the k+1;k following Bayesian hint: Criterion C7. A quadrature knot xk is assumed to belong to a neighborhood (xk 1; xk+1) inside which the integrand f (x) is continuous provided the rst order divided di erences satisfy min jdk;k 1j; jdk+1;kj > max jdk;k 1j; jdk+1;kj ; (9) where is an empirically de ned threshold. Under infringement of (9), the analysis follows the general pattern (iii), Section 3. If Criterion C7 is ful lled, then the analysis is re ned using Criterion C8. If the integrand pro le is monotonic and the curvature keeps constant sign everywhere inside ( ; ), then the activation of the local quadrature rules is accepted irrespective of the subrange width. Criterion C9. If the integrand pro le is monotonic and the pattern of the sign of the curvature over a sequence of three consecutive intervals is either + + or + , then a turning point with nite lateral derivatives has to be resolved or disproved. 11. Mesoscopic analysis of the boundary layer 17, 18]. Assume that f (x) is a continuous twice di erentiable function over a; b] and let xr 2 a; b] denote a reference argument value. Then there exists a non-vanishing neighborhood V (xr ) a; b] of xr inside which a linear Taylor series expansion holds true within a prede ned threshold 0 < " 1. Numerical check of the continuity of f (x) is done from a sampling of its computed values,
n

fi = f l f (xi ) i = 0; 1; : : : ; m ;

o

over a set of machine number arguments Sm (xr ) = xi 2 V (xr ) i = 0; 1; : : : ; m , m > 3, chosen such that f l(xr ) 2 Sm (xr ), where f l( ) denotes the oating point representation of 2 R . Let f (xi ) j i = 0; 1; : : : ; m denote the set of actual values of f (xi ) over Sm (xr ). In general, due to the round-o , f (xi ) f l f (xi ) 6= 0, hence the best information on the smoothness properties of f (x) at xr following from the set xi ; fi is obtained from the scrutiny of the properties of a second degree polynomial least squares t to the oating point data. It is convenient to perform the scale transformation xi = x0 + i hr , i = 0; 1; : : : ; m, i 2 Z where hr , denotes the distance from xr to its nearest machine number inside a; b]. This leads to the second degree tting polynomial y2 (xi ) = 0 + 1hr p1( i ) + 2 h2p2 ( i ) , spanned by the orthogonal basis polynomials pk ( i ), r k = 0; 1; 2, of norms Nk , respectively. Under negligible 2, the rst order derivative of f (x) at xr is given by m X 0 f 0 (xr ) y2 (xr ) = 1 = N1 1 p1 ( i )fi . i=0 The smallest sampling Sm (xr ) suitable for a least squares analysis providing insight on the smoothness properties of f (x) at xr = a and xr = b respectively consists of four distinct abscissas (i.e., m = 3). We choose them such that the set fx0; x1; x2g de nes a uniform mesoscopic mesh 0 = 0, 1 = p, 2 = 2p, 3 = q, q 6= f0; p; 2pg. Then the validity of a linear Taylor expansion around the reference abscissa xr is found to hold true within prescribed accuracy " provided the minimal sampling yields scale invariant approximations of the rst order derivative f 0 (x). Details and implementation are reported in 17, 18].


Numerical Methods and Programming, 2009, Vol. 10 (http://num-meth.srcc.msu.ru)

397

12. The priority queue associated to the binary subrange tree. Within the standard automatic adaptive quadrature, the magnitude of the local quadrature error e provides the simple key pointing to the subrange to be bisected next. Within the Bayesian automatic adaptive quadrature, consistent subrange handling is secured by a composite priority queue key. The magnitude of the local quadrature error is the primary priority queue key which secures the storage of the subrange showing the largest local error at root. For a subrange in unde ned state, the conventional value e = of low, where of low is a value near to the machine over ow threshold. The depth of the terminal nodes in the binary subrange tree provides the secondary priority queue key. If a smaller depth is detected in the depth list, then the corresponding subrange is moved at the root with two important consequences: (1) systematic elimination of spurious well-conditioning diagnostics over large subranges; (2) consistent extrapolation procedure activation. The subranges the local quadrature errors of which fall below a signi cance threshold are ruled out from the priority queue. This keeps the priority queue length to a minimal value. 13. Conclusions. A review of the principles of the Bayesian automatic adaptive quadrature has been given. Special emphasis was put on the description of the integrand properties which result in conspicuously unreliable local quadrature rule (q; e) outputs prior to their e ective computation. Consistency criteria enabling Bayesian hints are formulated. Detailed considerations of such criteria in separate publications pointed to the substantial increase of the reliability of the automatic globally adaptive quadrature under their use. Acknowledgments. The authors acknowledge the partial nancial support from pro ject CEX05{D11{ 68/11.10.2005, the Hulubei{Meshcheryakov Programme, JINR Orders 324/28.05.09, p. 30, and 328/29.05.09, p. 8, and the LIT-JINR theme 05{6{1060/2005{2010. References
1. Piessens R., deDoncker-Kapenga E., Uberhuber C.W., Kahaner D.K. QUADPACK, a subroutine package for automatic integration. Berlin: Springer Verlag, 1983. 2. Davis P.J., Rabinowitz P. Methods of Numerical Integration, Second edition. Orlando (Fla), USA: Academic Press, 1984. 3. Krommer A.R., Ueberhuber C.W. Computational Integration. Philadelphia: SIAM, 1998. 4. Lyness J.N. When not to use an automatic quadrature routine // SIAM Rev. 1983. 25. 63{87. 5. Adam Gh., Nobile A. Product integration rules at Clenshaw{Curtis and related points: a robust implementation // IMA J. Numerical Analysis. 1991. 11. 271{296. 6. Adam Gh. Case studies in the numerical solution of oscillatory integrals // Romanian J. Phys. 1993. 38. 527{538. 7. NAG Ltd., NAG Fortran Library Manual | Mark 17, Oxford, UK, 1996. 8. Visual Numerics Inc., IMSL MATH/LIBRARY: User's Manual | Version 3.0, Houston, TX, 1994. 9. Plakida N.M., Anton L., Adam S., Adam Gh. Exchange and spin- uctuation mechanisms of superconductivity in cuprates // Zh. Exp. Teor. Fiz. 2003. 124. 367{378; English transl.: J. Exp. Theor. Phys. 2003. 97. 331{342. 10. Marchetti P.A., Zhao-Bin Su, Lu Yu U(1) SU(2) Chern{Simons gauge theory of underdoped cuprate superconductors // Phys. Rev. B 1998. 58. 5808{5824. 11. Marchetti P.A., Jian-Hui Dai, Zhao-Bin Su, Lu Yu. Gauge eld theory of transport and magnetic relaxation in underdoped cuprates // J. Phys.: Cond. Matter 2000. 12. L329{L336. 12. Chu C.W. High Temperature Superconducting Materials: Present Status, Future Challenges, and One Recent Example | the Superconducting Ferromagnet // Physica C 2000. 341{348. 25{30. arstoiu F., Florescu A., Greiner W. Role of the higher static deformation of fragments 13. Sandulescu A., Misicu S., C^ in the cold binary ssion of 252 Cf // Phys. Rev. C 1998. 57. 2321{2328. 14. Adam Gh., Adam S. Increasing reliability of Gauss{Kronrod quadrature by Eratosthenes' sieve method // Computer Phys. Commun. 2001. 135. 261{277. 15. Adam Gh., Adam S. Reliability Conditions in Quadrature Algorithms // Computer Phys. Commun. 2003. 154. 49{64. 16. Justice J.H. Maximum Entropy and Bayesian Methods in Applied Statistics. Cambridge: Cambridge University Press, 1986. E.T. Jaynes. Bayesian methods: An introductory tutorial. 17. Adam Gh., Adam S., Tifrea A., Neacsu A. Resolving thin boundary layers in numerical quadrature // Romanian Reports in Physics. 2006. 58, N 2. 107{122. 18. Adam Gh., Adam S. The boundary layer problem in Bayesian adaptive quadrature // Physics of Particles and Nuclei Letters. 2008. 5, N 3. 269{273. 19. Press W.H., Teukolsky S.A., Vetterling W.T., Flannery B.P. Numerical Recipes in Fortran 77 | The Art of Scienti c Computing. Second Edition. Cambridge: Cambridge University Press, UK, 2001.

Received November 5, 2009