Numerical quadrature : a comparison of available procedures at the Computing Centre of the Eindhoven University of Technology

at the Computing Centre of the Eindhoven University of

Technology

Citation for published version (APA):

Geurts, A. J., & Haan, den, P. J. (1981). Numerical quadrature : a comparison of available procedures at the Computing Centre of the Eindhoven University of Technology. (EUT report. WSK, Dept. of Mathematics and Computing Science; Vol. 81-WSK-04). Technische Hogeschool Eindhoven.

Document status and date: Published: 01/01/1981

Document Version:

Publisher’s PDF, also known as Version of Record (includes final page, issue and volume numbers)

Please check the document version of this publication:

• A submitted manuscript is the version of the article upon submission and before peer-review. There can be important differences between the submitted version and the official published version of record. People interested in the research are advised to contact the author for the final version of the publication, or visit the DOI to the publisher's website.

• The final author version and the galley proof are versions of the publication after peer review.

• The final published version features the final layout of the paper including the volume, issue and page numbers.

Link to publication

General rights

Copyright and moral rights for the publications made accessible in the public portal are retained by the authors and/or other copyright owners and it is a condition of accessing publications that users recognise and abide by the legal requirements associated with these rights. • Users may download and print one copy of any publication from the public portal for the purpose of private study or research. • You may not further distribute the material or use it for any profit-making activity or commercial gain

• You may freely distribute the URL identifying the publication in the public portal.

If the publication is distributed under the terms of Article 25fa of the Dutch Copyright Act, indicated by the “Taverne” license above, please follow below link for the End User Agreement:

www.tue.nl/taverne Take down policy

If you believe that this document breaches copyright please contact us at: openaccess@tue.nl

providing details and we will investigate your claim.

ONDERAFDELING DER WISKUNDE EN INFORMATICA DEPARTMENT OF MATHEMATICS AND COMPUTING SCIENCE

N~i1[erical Quadrature

A comparison of available procedures at the

Computing Centre of the Eindhoven University of Technology by

A.J. Geurts p.J. den Haan

T.H.-Report 81-WSK-04 Augustus 1981

Computing Centre of the Eindhoven University of Technology. The pro-cedures are compared wit:i respect to the criteria: reliability and efficiency.

The results may be of help for a user, dealing with an integral that must be computed numerically, to select a procedure from the available

software.

1. Introduction

Comparison of a number of individual items of numerical software, such as programs or procedures/routines, all of which are intended to solve the same group of probLems, may be done for several reasons, for in-stance:

- to show that, or under which conditions, a new developed procedure is superior to a number of already existing ones;

- to select a subset out of a set of procedures, e.g. in constructing a program library;

- to supply information with the help of which someone (a user) can select for his problem an appropriate procedure from a collection of available procedures (a program library).

As is well known, it is very difficult to set up a comparative test such that the results are of practical use; the major difficulties being:

simplification; the appreciation of a piece of software depends on a large number of aspects. But only a few of them, often simplified, are involved in the test, and these may even not be the relevant ones;

extrapolation; a comparison is performed on a set of test problems, the conclusions are of value only if they hold outside this set; interpretation; a comparison yields a lot of data. These data must be compressed into a readable form without loosing too much infor-mation. The information from which conclusions have to be drawn may be indistinct, a concl~sion may be very personal.

In this report we describe a comparative test performed on the set of procedures for automatic quadrature available in Algol for general use on the B7700 at the Computing Centre of the Eindhoven University of Technology. The:'goal ef the test is the third above mentioned, that is, to provide information for the users of the computing centre which may help them to make their choice out of the software on hand.

Numerical quadrature belongs to the class of so-called "unreliable arithmetic algorithms", cf. Lyness [6J. This means that, in principle, the algorithm does not solve the problem, but rather furnishe~ an answer that might be an approximation of the required solution. In ac-tual practice, the situation with numerical quadrature is not as bad as suggested; there are procedures which quite often give an approxi-mation within the required accuracy. Yet there is a fundamental

uncer-tainty about the obtained result, unless extra information, eg., the magnitude of some derivative, available. "So the question is not whether the program is right or wrong, but rather how well it performs" Rice [10, p.5J.

The choice of a procedure for a particular problem made by a user is subjective and possibly based on both quantitative criteria (eg., CPU-time, storage requirement, complexity, failure probability) and non-quantitative criteria (eg., readability of the documentation, users convenience, confidence in the designer). A comparison test is necessarily restricted to quantitative criteria, although non-quan-titative criteria might be more pertinent to the user.

The test presented in this report treats the following two criteria: - reliability, measured by the fraction of successes within a test

set of integrals.

efficiency, measured ty the number of function evaluations ~n the computation of an integral.

With respect to reliability we register not only whether a result is correct or not, but also the indication given by the procedure (if any)

about the correctness of the result. If this indication is adequate, it should be taken into account by the user, if not, the user should not be troubled with it. We give our conclusion on the adequacy of the

considered indications as well.

With respect to efficiency we agree that the number of function eval-uations is only a substitute, since it is the total CPU-time needed for the evaluation of the integral that counts. Furthermore, the use of the former to measure the latter can be misleading (Lyness [6J, Robinson [llJ). On the other hand, it is very difficult, if not im-possible, to measure the real CPU-time in a computing system with time-sharing facilities, and measuring time by allowing the test pro-gram the sole and uninterrupted use of the computer (Robinson [11]) may bt misleading as well.

Another practical argument for taking the number of function evaluations as a measure of efficiency is that it allows one to spread the compu-tations of the test in time, since it is to some extent insensitive to modifications of the computing system, as long as the arithmetic re-mains unchanged. Even the addition of a new procedure afterwards can

easily be done.

An important tool for a test is the set of problems on which the test is performed. Test prob:ems may be distinguished in "real life" and "artificial" problems, the former generally hard to obtain and often expensive in computer time, the latter probably too specific and con-sequently unrealistic. In any case the test set will depend on the goal of the test. We think that for our goal "real life" problems would be ideal, i.e., perfect, but not likely to be achieved. So we have chosen for sets of artificial problems. And we hope they are not too far from real life.

The test method used belongs to the category of the so-called battery experiments, the pros and cons of which are discussed several times in literature, cf. Lyness [7J, [8J, Robinson [IIJ. A major difficulty with the battery experiment is the abundance of numbers the reader confronted with, cf. Kahaner [5J, Robinson [11J. We use the idea of a pairwise comparison matrix, cf. Saaty [12J, to compress the information with respect to efficiency.

In section 2 a description is given of the procedures considered. Tn section 3 we describe how the population of test integrals is built up. The frame-work of the test method is explained in section 4. In section 5 and 6 we give numQrical results with our conclusions.

Robinson [llJcontains a test set of integrals which is comparable with the set used for our test. We used this set to check the validity of our results. This check is reported in section 7.

2. Description of the procedures

The set of procedures to be compared in this report consists of all ALGOL procedures, available for general use on the B7700 at the Com-puting Centre of the Findhoven University of Technology, that claim

to compute the integ .. :al

b

I :=

f

f(x)dx

a

for finite a and b and a "general" integrand f, with a specified accu-racy.

In this section a description of these procedures is given, consisting of a global indication of the method used and furthermore answering the following questions:

- is the algorithm adaptive *);

- does the procedure ask for an absolute or a relative tolerance. In the latter case, is the tolerance relative to

b

f

f(x)dxl a b or to

I

If(x)ldx; a

- does the procedure indicate whether the required accuracy has been reached, and if it does, in which way;

- what is the maximum number of function evaluations allowed by the procedure.

*'A procedure is called adaptive if the sequence of points for which the integrand is evaluated, depends on the function values encountered.

The B7700 NUMERALS, Numerical Analysis Program Library contains the following three procedures for automatic quadrature. For a complete description and references, see [14].

ROMBERGINT.

Computes the integral by means of Romberg-extrapolation with polynomials, combined with interval bisection. The algorithm is non-adaptive. The

b

procedure asks for a relative tolerance with respect to

f

If(x)ldx.

a

There is no indication on whether the required accuracy has been reached. At most 15 subdivisions are applied, so the maximum number of function

1 . . 215 1 32769

eva uat10ns 1S + = •

BULIRSCHINT.

Computes the integral by means of Romberg-extrapolation with rational functions and subdivisions according to the so-called Bulirsch sequence (h := (b - a)/N, with N = 2,3,4,6, ••• ). The algorithm is non-adaptive. The tolerance is the same as in ROMBERGINT. Whether the desired accu-racy has been attained, might be concluded from the parameter ERR of the procedure. This parameter contains, in case of convergence, a con-servative error bound and with no convergence the value 10 (but this is not mentioned in the documentation). The maximum number of

subdivi-sions is IS, so the maximum number of function evaluations is 3 x 27 + I

=

385.

SIMPSONINT.

The procedure is based on the adaptive Simpson method of McKeeman, but employs interval bisection in place of trisection. The tolerance is the

same as ~n ROMBERGINT. There is no indication on whether the required accuracy has been reached. Since the depth of recursion is maximal 24J the number of function evaluations is bounded, but may be very

large (less than 226 +

~).

The local library of the university, called THE LIBRARY. contains four b

procedures for the computation of

J

f(x)dx with a prescribed accuracy.

a

For a complete description and references, see [16J.

INTEGRAL TRAPEX.

The procedure is based on Romberg-extrapolation with rational functions and interval subdivisions according to the Bulirsch sequence. So the algorithm is similar to BULIRSCHINT. The procedure asks for an absolute and a relative tolerance,

AE

and RE respectively, and the convergence criterion then is

b

lOll ~

AE

+ RE x

J

f(x)dxl·

a

A Boolean parameter reveals whether the required accuracy has been reached. The procedure asks for the maximum number of interval sub-divisions, say M. The maximum number of function evaluations then is ZM/2+1 + 1 if M is even and 3 x 2(M-l)/2 + 1 if M is odd.

INTEGRAL.

An adaptive procedure based on Simpson's rule with varying step size. An absolute and a relative tolerance AE and RE, respectively, are

required. The convergence criterion in the n-th step then is

x

n

lOIn'

5: AE + RE x

J

If(x)ldx.

X_{n -1}

A minimal steplength HMIN must be supplied. The procedure returns the number of times an interval of length HMIN is encountered which does not satisfy the convergence criterion. This number is used as an indi-cation with respect to the attained accuracy. On account of HMIN the number of function evaluations is bounded, but may be very large.

QADRAT •.

An adaptive procedure based on Gaussian quadrature formulae of high order and a recursive interval bisection strategy. Tolerances and error bound are the same as those of INTEGRAL.

The procedure uses an internal minimal steplength HMIN which is pro-portional to RE, viz. HMIN

=

32 x RE x (b - a). Similar to INTEGRAL the procedure returns the number of subintervals of length HMIN for which the required accuracy is not reached. So with respect to the accuracy indication and the number of function evaluations, see INTE-GRAL.

INTEGRML.

An adaptive procedure based on a 5-point Gauss-Legendre quadrature formula and a varying step size strategy. The procedure asks for an

*)

absolute tolerance. The procedure returns a realistic approximation to the error in the computed value, which indicates whether the de-sired accuracy has been obtained. The procedure uses a fixed internal minimal steplength FJ:.llN = 10-4 x Ib - al. The maximum number of

func-tion evaluafunc-tions therefore is approx. 105•

The NAG LIBRARY Mark 7 (ALGOL 60) contains the following four proce-dures for automatic numerical quadrature. For a complete description and references, see [15J.

D01AAA.

Computes the integral by means of the Clenshaw-Curtis method with the O'Hara-Smith error estimate. This method is based on finite Chebyshev series approximations, and is non-adaptive. The procedure requires a tolerance relative to the value of the integral. A parameter IFAIL in-dicates whether the accuracy has been reached. The maximum number of function evaluations is M + 1 where M. the maximum number of subintervals. which should be a power of 2 not greater than 128, must be supplied.

Remark. The weights and abscissae used in the Chebyshev series are cal-culated before the evaluation of the integral. Subsequent calls to the procedure bypass this calculation. Therefore the first call of the procedure is much more expensive than subsequent calls.

*) This is a quantity that is generally of the same order of magnitude. It needs not to be an upper bound, and it may occasionally be a bad approximation, especially when the accuracy is not reached.

DOIABA.

The procedure is based on Romberg-extrapolation with polynomials and interval bisection,cf. ROMBERGINT. For tolerance and accuracy indica-tion, see DOlAAA.

The maximum number of function evaluations is M + 1 analogous to DOlAAA. In this case M may not exceed 1024.

DOIACA.

Computes the integral by means of a method described by Patterson. This method uses a family of quadrature formulae of increasing order over the whole interval. so it is non-adaptive. The procedure requires an absolute and a relative (to the integral) accuracy and will termi-nate whenever either accuracy has been reached.

The accuracy indication is realized by the parameter IFAIL. The pro-cedure supplies a realistic approximation to the error in the computed integral as well. The maximum number of function evaluations in 255.

Remark. The procedure makes no use of the values of the integrand at the endpoints of the interval.

DOIAGA.

This procedure evaluates the integral by the Clenshaw-Curtis method and a double-adaptive interval subdivision strategy. An absolute tole-rance must be supplied. The accuracy indication is realized by the parameter lEAIL, furthermore a realistic approximation to the error in

the computed value is returned.

The procedure uses a minimum steplength HMIN

=

Ib - al 2-m, where m must be supplied by the user. So the number of function evaluations

is bounded, but may bp very large.

The description given are summarized in the following table.

Method adap- re re acc. error max

tive ae f(x) If(x)1 indic. estim. fu.ev.

ROMB ERG INT Romberg + 32769

BULIRSCHINT Romberg-Bulirsch + + + 385

SIMPSONINT Simpson + + < 00

INT. TRAP EX Romberg-Bulirsch + + + 2 !11+ 1 + 1

INTEGRAL Simpson + + + + < 00

QUADRAT Gauss + + + + < 00

INTEGRAAL Gauss + + + + < 00

DOlAAA Clenshaw-Curtis + + M+I

DOIABA Romberg + + M+I

DOIACA Patterson + + + + 255

3. Description of the population of test functions

As already has been pointed out in the introduction, we have used "artificial" problems. We gathered these problems from test sets of others, from textbooks dnd tables, etc. We have chosen functions which are finite in the ivterval of integration, easy to evaluate and of which the integral is exactly known, or at least to machine precision.

Like others (Robinson [11J, Lyness [7J) we believe that a population of all sorts of integrable functions is too broad and too various to be used as a population for the comparison of quadrature procedures. So we distinguished the following groups and performed comparisons on each group separately (cf. Robinson [lIJ).

I. Smooth or well-behaved functions.

2. Functions with one or more sharp peaks within the interval. 3. Functions with a singularity in (one of) the endpoints. 4. Functions wi th a singularity wi thin the open interval. S. Rapidly oscillating functions.

During the computations we encountered a number of functions for which nearly all the procedures failed. These were unbounded functions (that is why the functions in all five groups are finite) and a number of functions we have gathered in a separate class of miscellaneous

func-tions. We called class

6. Pest functions.

In order to enlarge the number of functions in a simple way and to avoid to some extent partiality in the choice of test functions we have

used the following transformation, due to O'Hara and Smith, cf [9J,

(3.1 )

:= (qb - a)x - (g - I)ab

y (q - I)x - ag + b

This transformation ~ps, for any positive value of the parameter g,

the interval [a,bJ onto itself, and if the function g(q,x) is defined by (3.2) then a g(q,x) b .- fey)

.:!l

dx 2 q(b - a)

«q -

l)x - aq + b)2 b

J

f(x)dx =

J

g(q,x)dx a f (~b - a)x - (q - l)ab) \ q - 1) x - aq + b '

For the effect of the transformation on the integrand, see Van der Vorst [l3J.

Remark. From the detai.led results of the tests it follows that integrals corresponding to the same basic function and different values of q

usually behave similar, but there are cases with considerable differences. So in our opinion the use of the transformation (3.1) makes sense.

b

Since the relative tolerance is related either to I

f

f(x)dxl or to

b a

f

If(x)ldx we used only positive functions (if necessary, we have made

a

that if f is positive, then g is positive as well. Another argument for the choice of positive functions is that the integral of a posi-tive function is well-conditioned with respect to perturbations in the function values.

A'list of all test functions, arranged 1n the various groups, 1S given in Appendix A.

4. Framework of the test

In order to obtain a fair comparison the parameters of the va~ious

procedures must be chosen with care. The specific choices we have made for that reason are the following:

- All integrands are positive functions.

If possible, then a relative tolerance £ is chosen; if not, then the absolute tolerance £ fbf(x)dx is used.

a

- Some procedures ask for the maximum number m of interval subdivisions, others ask for a minimal steplength h . • Sometimes the documentation

m1n

gives restrictions, sometimes recommendations, for the values of these parame ters • documentation recommends/demands we took INT. TRAPEX m::; 10 20 INTEGRAL h min ~ 10-5 Ib - al 10-6 Ib - al DOlAAA 8 ~ m 128 128 DOIABA 4 ::; m ::; 1024 1024 2-m Ib - al 2

DOIAGA h min

=

m

= -

log(£)

A test set consists of a group of Nintegrals, and one tolerance. In a test run all integrals are evaluated with tolerance E by all (eleven) procedures. This yields the following 11 x N test results:

{I, ~I, IND, NFE}.

I the computed approximation of the integral I

t,1 the relative error, viz. &1 := 1([ - 1)/11 IND the accuracy indication, default value is true NFE the number of function evaluations.

From these results we first derive a reliability score. We distinguish the following four cases, coded from 0 to 3:

code accuracy obtained

o

2

3

The reliabili ty score (in

%)

is table. Example. Group 4, £

=

J 0-9 , N Method 3 SIMPSONINT 4 INT. TRAPEX 5 INTEGRAL 6 QADRAT 7 INTEGRAAL I 1 DOIAGA true true false false registered 200. 0 95 I 1 3 45 51 46 18 41 68 indication true false false true in a so-called reliability 2 3 5 82 4 37 18 3 41 4 27

From a reliability table we may derive conclusions about the reliability of the procedure in solving the quadrature problem (columns (0 + 1) versus columns (2 + 3»; and about the reliability of the accuracy

indicator (columns· (0 + 2) versus columns (1 + 3».

With respect to effLciency we dislike the presentation of a bulk of numbers, from which everyone can, and indeed has to, make his own com-parison, cf. Kahaner [5J, Robinson [IIJ. We would rather express

the relative efficiency of one procedure with respect to another by a single value. To this end we proceed as follows.

Let P

k and P 9, be two procedures to be compared. Let J := {Ii

11

s i s N} be a set of integrals correctly computed by both P_k and P~. Let Fmi be the number of function evaluations used by P , m = k,9" for the

m

computation of Ii' and finally let Qk9,i := F_ki/F_2i, l s i s N. Then the relative efficiency of P_k with respect to P

2 may be expressed by the geometric mean

(4.1)

Now let vm' m k,!(', be defined by

N 1

(4.2) v _m :=

(i~l

F .

\N"

mL

)

then

(4.3) _{Ek2 vk/v9,'}

So we may assign to each procedure the corresponding value (4.2) as an efficiency number, of which only the relative magnitude counts.

Now if the procedures P , m = 1,2, ••• ,M are to be compared and we use

m

one and the same test set J, then for any pair :k and P£ the above statements hold. This implies, in particular, that we may use the num-bers v , m _m = 1,2, ••• ,M as an efficiency score to compare the procedures mutually.

Taking one and the same problem set J for all procedures would quite probably imply a considerable loss of information, since this set, being the set of integrals correctly computed by all procedures, will often be much smaller than any of the sets correctly computed by a pair of

procedures. This loss of information may be avoided in the following way.

Let J be a given set of N problems. Let Jk~ be the subset of J consis-ting of all integrals correctly computed by both P_k and P₂ and let Nkt be the number of problems in J

k£* Let Ek£ be defined by

(4.4)

where the product is taken over all i with Ii E J

k£. The numbers Ekt,

I ~ k,~ ~ M constitute a positive matrix

(4.5)

In Appendix B we show that the positive eigenvector, w say, correspon-ding to the eigenvalue p(E) m~y be used as a replacement for the vector v defined by (4.2), and that

(4.6) fl : = LP~( E~):.-...-....,...;;.:;M M-]

called consistency, is a workable measure for the quality of this replacement.

As a result, the efficiency score of a test set consists of the eigen-vector w corresponding to the eigenvalue p(E) of the matrix E, and the number ~. We record the results in a so-called efficiency table. The lower triangle ~ ~ < k ~ M contains the fractions 100 x Nkt/N,

the upper triangle ~ k < £ ~ M contains the matrix elements E kt, which describes the entire matrix, since Ek£ ==

E~~

for all k.L The

( M )

third last column contains the mean values of the rows, i.e., \

I

Ek£

1M

t= I for the k-th row, as an initial quess for the eigenvector w. The second last column contains w. These two columns are normalized such that the minimum value equals 1. The difference between the mean value sum and w is also an indication for the consistency. The last column contains 2log w. The value of p is given separately.

Example:

EFFICIENCY TABLE

INPUT IDENTIFICATION ! GI~OUl;:' 4

EPS=@'-9 j.94 NUMBER OF SELECTED INTEGRALS:

METHOD 3 4 5 7

F<EL • EFF • NUMBI::F:S

MEAN EIG. LOG.

ROW VEC. SCAt

SUM

---_._--3: SIMPSONINT 0.2 0.5 :L.O 0.5 1.6 1.6 1.6 0.69

4:

INT.TRAPEX 14 _{..: .. +} '') .! \J 7.5 2.9 8.5 8.5 8.4 3.06 5: INTEGRAL 46 13 2.4 1.0 2.9 3.1 3.1 1.61 6: OArtFMT 98 14 46 0.4 1.5 1.5 1.4 0.46 7: INTEGF~AAL 60 :1.3 40 60 2. <J 3.1 3.0 1. 60 11: D01AGA 70 11 :35 70 45 1.0 1.0 0.00 CONSISTENCY: 6.163E··-03

In establishing the reliability and efficiency we assumed the computed integral to be correct if the actual error is less than the given tolerance 8. Many users are less rigorous, that is they rather accept

a result if the error is not much larger than the tolerance. To get an idea of the infl~~nce a less stringent acceptance criterion may have on the results of the test we investigated how often the error in the computed value is only a little larger than E. To this end we

have considered, for every problem group and every tolerance, the set of integrals which are evaluated with an error greater than E, and

we determined the percentage of integrals in this set with an error less than 2E. These percentages are given 1n the following table.

group e

=

10-3 25 10-6 19 10-9 ₄₆ 2 6 3 5 3 42 18 7 4 28 9 5 5 all 9 18 8 9 5 7

The overall percentage amounts to 10%. So we see that errors larger than the tolerance are mostly considerably larger, and the results of the test will hardly be influenced by a small enlargement of the accep-tance of the computed result.

5. Reliability and efficiency, numerical results

In this section we give the numerical results with respect to the cri,teria reliability and efficiency, We remind that there are five

groups of integrals, namely:

1. Smooth or well-behaved functions.

2. Functions with one or more sharp peaks within the interval. 3. Functions with a sLngularity in (one of) the endpoints. 4. Functions with a singularity within the open interval. 5. Rapidly oscillating functions.

For every group we have approximately 20 basic integrals; these in-tegrals are listed in Appendix A. From every basic integral we have derived 10 integrals by means of transformation (3.1) with the follow-ing values of q : 0.33, 0.50, 0.75, 0.78, 1.00, 1.05, 1.25, 1.50, 1.63 and 2.14. So every problem group contains approximately 200 in-tegrals.

-3 -6 -9

We have used three tolerances: £

=

10 ,10 and 10 ,the last one being about the smallest possible, since the machine-precision of the Burroughs B7700 equals

~

x 8-12

~

0.7 x 10-11.

On each of the pages 23-27 the results concerning reliability

and efficiency of one group are given. The upper half of a page con-tains the reliability tables of all three tolerances. For convenience of the reader we repeat that the columns of a table contain the following scores, in percentages:

Column 0 integrals correctly computed, agreed by the indication used; Column integral correctly computed, disagreed by the indication used; Column 2 integral incorrectly computed, failure indicated by the

pro-cedure;

Column 3 integral incorrectly computed, failure not indicated by the procedure.

The lower half of the page contains histograms of the last columns of

2

the efficiency tables, i.e. log w, with the methods arranged in de-creasing order of efficiency. If for some method the percentage of correctly computed inte;rals is too low (actually less than 5%), then this method is excluded from the histogram. We repeat that for a pair-wise comparison only the difference of the corresponding values of

2

log w counts: a difference of I cm in height of two bars in a histo-gram corresponds to a factor

1:2

in efficiency, in favour of the method corresponding to the lower bar.

With these results it LS easy to compare the performances of all

methods together with respect to one problem group.

On the next pages we discuss the performances of the methods on account of these results. The discussion is far from exhaustive; the recommen-dation of some procedures for each group, on account of both reliability and efficiency, is rather personal.

Finally we give the complete efficiency tables of every test set of integrals.

METHOD 1: 2: 3: 4: 5: 6: 7: 8: 9: 10: 11: ROMBERGINT BULIRSCHINT SIMPSONINT INT.TRAPEX INTEGRAL QADRAT INTEGRAAL D01AAA DOl ABA n01ACA D01AGA E= -3 10 0 100 49 100 99 100 100 100 100 99 100 100 2810411693715 ~ ~ Reliability tables .~ -e-3 e-6 1 2 3 0 1 2 100 27 24 10 48 42 100 1 100 100 100 94 100 :I. 100 100 100 Efficiency histogra~s 2861110479315 e-9 3 0 1 2 3 l.OO 7 59 34 100 100 97 3 100 6 9"'" ,J 5 100 99 1 100 100 £:= 10- 9 r- - 4 - 3 - 2 ,...""- _ 1 8 6 11 2 10 4 7 9 1 3 5

METHOD 1: 2: 3: 4: 5: 6: 7: 8: 9: 10: 11: ROMBERGINT BULIRSCHINT SIMF'SONINT INT.TRAPEX INTEGRAL QADRAT INTEGRAAL D01AAA D01ABA D01ACA D01AGA E;= -3 10

-386115107491 0 98 1 99 54

leO

23 66 13 49 28 73 (!-3 1 2 3 0 1 2 80 1 90 8 6 1 100 2 35 9 48 3 100 19 47 11 100 12 22 78 20 1 3 83 7 2 4 46 1 17 21 21 51 12 17 11 16 99 Reliability tables £= -6 10

- r-- ..-.--~ r--

,..

-l -l 6 8 1 0 3 7 5 4 9 1 Efficiency histograms (!-6 2 3 20 94 48 1 2 86 5 62 '71 1 0 78 96 30 34 87 46 6 7 84 s= -9 10

r

-1 1 14 4 8 13 8

-(!-9 2 3 22 99 4 56 8 10 c' >J 4~S 3 93 1 87 8~j 16 - - 5 -4 _ 3 -2 - 1

METHOD 1 : 2: 3: 4: 5: 6: 7: 8: 9: 10: 11: ROMBERGINT BULIRSCHINT SIMPSONINT INT.TRAPEX INTEGRAL QADRAT INTEGRAAL D01AAA D01ABA D01ACA D01AGA E;:: -3 10 0 100 14 :1.00 '72 94 67 100 98 98 100 99 10 8 6 11 2 4- 3 9 1 7 5 Reliability tables ~-'3 ~--6 1 2 3 0 1 2 95 9 73 4 1 1 97 1.00 28 55 29 6 46 40 9 7 26 44 56 70 29 2 FJ6 5 l 1 42 10 48 100 1 97 E,fficiency histograms

Q

r - .r -

,'

,

-I

8 10113 6 7 5 4- 9 1 ~-9 3 0 1 2 5 C::'") >J"-1 1 99 95 16 lO 2 86

...

_{>J 2~~~ 7 46} ~56 44 1 12 56 3") .:.. 9 34 1.4

c''''

>J.,. 7 2 9l 71 5 24 3 98 s= 10-9

Q

_,....

, - '----

,

-

--810116735491 ~~ 48 5 2 '")I:!' ~ .. .J 2 4 _ 3 - 2 _ 1 - 0

Reliability tables

e-3 e-6 e-9

METHOD 0 1 2 3 0 1 2 3 0 1 2 3 1: ROMBERGINT 97 3 94 6 50 50 2: BULIRSCHINT 12 9 75 4 2 2 96 1 1 98 3: SIMPSONINT 100 100 95 0::' .J 4: INT.TRAPEX 73 27 33 10 38 19 11 3 82 4 5: INTEGRAL 96 4 70 14 16 45 37 1.8 6: QADRAT 48 12 1 39 43 53 2 2 51 46 3 7: INTEGRAAL 89 11 46 43 6 5 18 41 4:1, 8: D01AAA 34 66 2 8 90 1 87 12 9: D01ABA 93 1 6 34 20 45 1 5 7 88 10: D01ACA 82 1 4 13 15 7 76 2 1 99 11: D01AGA 73 1 26 67 2 3 28 68 1 4 27 Efficiency histograms E= 10-6

o

r-- 4

-r - _ 3

r

- r-r--- 2 r-r - r--- 1 r -

r-[

I

_{- 0} 8362J1Jt91D175 3J16571D941

n

6 3 7 5 9 4 1

METHOD 0 1: ROMBERGINT 96 2: BULIRSCHINT 7 3: SIMPSONINT 87 4: INT.TRAPEX 79 5: INTEGRAL 96 6: QADRAT 66 7: INTEGRAAL 90 8: D01AAA 73 9: D01ABA 82 10: D01ACA 100 11: [l01AGA 91 r--_ r - " ....-,...-

r-

,...

-8 2 10 6 3 4- Jl 7 9 5 1 Reliability tables (!-3 (!-6 1 2 3 0 1 4 87 15 75 3 4 26 13 94 21 83 1 2 2 87 16 4 14 87 13 10 84 7 2 25 80 2 6 12 67 10 88 8 7 2 9:3 Efficiency histograms E: ::: 10- 6 ,..-r- ...- ....- -

r-,...

r-

.--I"

2 70 14 , 13 8 15 23 3 7 (!-9 3 0 1 2 3 13 81 19 7 36 57 6 86 14 2 71 5 20 4 48 39 85 10 5 1 74 5 20 1 3 72 2 21 5 42 15 43 1 76 1 22 1 95 1 4

Q

r " - ,..-_ 3 ,..- ₂

.--

...--

,...--1

-r

_-

o

8 10 6 2 11 7 4- 3 9 10 5 8 10 2 6 11 7 4- 9 1 2 5

Discussion General remarks.

I. For the sake of limitation of the CPU-time, procedures are not allowed to use more than 35000 function evaluations per integral. If this maximum has been reached, then the computation of the inte-gral is aborted and no results are registered. Consequently, the

Sum of the reliability scores may be less than 100%, see for in--9

stance: group 2, ~ = 10 ,INTEGRAL or QADRAT.

2. The procedure BULIRSCHINT is useless because of its bad reliability score for all groups. Therefore BULIRSCHINT will be left out of the discussion in the sequel.

3. The behaviour of QADRAT fore= 10-3 is remarkably bad in all groups except group 1. It turns out that, if the required accuracy has not been achieved, then the procedure has used only a few function values. This suggests that the termination criterion is not ade-quate for low accuracy.

4. The first call of DOlAAA requires a lot of extra computations (see page ]0). Therefore, in the discussion of the individual groups this procedure is not recommended in the cases when its efficiency score belongs to the best.

5. If the eigenvector w has a large component, corresponding to one procedure (e.g. ROMBERGINT), then the consistency may improve con-siderably if this procedure is omitted from the efficiency table. On the other hand, this omission hardly influences the other compo-nents and has therefore no effect on the assessment of the procedures.

2

Example. The consistency and the components of log w 1n group 3 with E = 10-9, with ROMBERGINT (see pag40, last column) and with-out ROMBERGINT are, respectively,

consistency I 3 4 5 6 7 8 9 10 I 1

6.51

10-2 6.61 2.26 3.31 3.24 0.98 1. 91 0.00 4.01 0.02 0.90

1.13

10-2

----

2.23 2.97 3.18 1.00 1.90 0.00 3.50 0.03 0.89

Remark::: on the individual groups.

Group I: Smooth or well-behaved functions. I. All procedures behave very reliably.

The procedures based upon Romberg-extrapolation with polynomials or Simpson's rule are significantly less efficient than the other ones.

2. Recommended procedures. QADRAT, DOIACA, DOIAGA.

Group 2: Functions with one or more sharp peaks within the interval. I. SIMPSONINT is very reliable, followed by ROMBERGINT and DOIAGA.

QPDRAT and DOIAGA are the most efficient procedures.

The non-adaptive procedures are evidently not suited for this type of integrals.

2. Recommended procedures.

DOIAGA; SIMPSONINT if reliability is highly wanted; QADRAT for high accuracy only.

Group 3: Functions with a singularity in (one of) the endpoints.

1. SIMPSONINT and DOIAGA are highly reliable, DOIACA is only for -9

E

=

10 less reliable, yet very faithful, i.e., indicates properly that the result is incorrect.

The procedures based upon Romberg-extrapolation are not suited for this type of integrals.

The most efficient procedure 1S DOIACA. Remind that this procedure

makes no use of the endpoints of the interval. 2. Recommended procedures.

SIHPSONINT, DOIACA, DOIAGA; QADRAT for high accuracy only. Group 4: Functions with a singularity within the open interval.

I. SIHPSONINT is the only procedure with high reliability for all tolerances. The rather bad reliability of DOIAGA for all tolerances is remarkable. This bad performance turns out to be neither re-stricted to some basic problems of the group, nor to specific values of q in the transformation (3.1). The incorrectly computed integrals with indication IFAIL

=

0 (column (3) of the reliability table) are more or less homogeneously distributed over the group.

The non-adaptive procedures are not suited for this type of il1te-grals.

SIMPSONINT, QADRAT and DOIAGA are significantly more efficient than the other procedures.

2. Recommended procedures.

5: Rapidly oscillating functions. ----~~

I. The procedures ROMBERGINT, SIMPSONINT, INTEGRAAL, DOIACA and DOIAGA are reliable for all tolerances.

The procedures based upon Romberg-extrapolation with polynomials or Simpson's rule are the least efficient. The most efficient of the reliable procedures is DOIACA.

2. Recommended procedures.

REL .. EFF.NUMBERS

MEAN EIG. LOG.

ROW VEC. SCA.

METHOD 1 2 3 4 5 6 7 8 9 10 11 SUM

---1: ROMBERGINT 2.1 1.3 1.9 0.9 1.8 1.0 2.1 1.4 2.0 1.9 2.:' 2.2 1.13 2: BULIRSCHINT 76 0.5 0.9 0.4 0.8 0.5 0.9 0.7 0.9 0.9 1.0 1.0 0.00 3: SIMPSONINT 100 76 1.5 0.7 1.4 0.8 1.7 1.1 1.6 1.5 1.8 1.7 0.81 4: INT.TRAPEX 99 75 99 0.5 1.0 0.5 1.1 0.7 1.0 1.0 1.1 1.1 0.20 IoU 5: INTEGRAL 100 76 100 99 2.0 1.1 2.3 1.5 2.1 2.0 2.4 2.4 1.25 N 6: QADRAT 100 75 100 98 100 0.6 1.2 0.8 1.1 1.0 1.2 1.2 0.29 7:: INTEGRAAL 100 76 100 99 100 100 2.1 1.4 1.9 1.8 2.1 2.1 1.10 8: D01AAA 100 76 100 99 100 100 100 0.7 0.9 0.9 1..0 1.0 0.04 9: [l01ABA 99 76 99 97 99 98 99 99 1.4 1.3 1.6 1.6 0.64 10: D01ACA 100 76 100 99 100 100 100 100 99 1.0 1.1 1.1 0.16 11: J)olAGA 100 76 100 99 100 100 100 100 99 100 1.2 1.2 0.22 CONSISTENCY: 1.827£-04

REL.EFF.NUMBERS

MEAN EIG. LOG.

ROW VEC. SCA ..

METHOD 1 2 3 4 5 6 7 8 9 10 11 SUM

---1: ROMBERGINT 3.0 1.0 2.1 0.7 3.2 1.9 3.3 1 .. 1 2.7 2 .. 9 3 .. 3 3.2 1 .. 69 2: BULIRSCHINT 58 0.3 On8 0 .. 2 0 .. 9 0.6 0 .. 8 0.4 0.8 0.9 1 .. 0 1.0 0 .. 00 3: SIMF'SONINT 100 58 2.1 0.6 3.1 1.8 3.2 1.1 2.6 2 .. 9 3.3 3 .. 2 1.68 4: INT.TRAF'EX 100 58 100 0 .. 3 1 .. 5 0.9 1.5 0.5 1.2 1.4 1.5 1.S 0.S8 w 5: INTEGRAL 100 S8 100 100 4 .. 9 2 .. 9 5.0 1.7 4.1 4.5 5.1 5.0 2 .. 33 w 6: QADRAT 100 S8 100 99 100 0.6 1.0 0.3 0.8 0.9 1.0 1.0 0.03 7: INTEGRAAL 94 55 94 93 94 93 1.7 0.6 1.4 1.6 1.7 1 .. 7 0.78 8: D01AAA 100 58 100 100 100 100 94 0.3 0 .. 8 0.9 1.0 1.0 0.02 9: [l01ABA 100 58 100 100 100 100 94 100 2.4 206 3.0 2 .. 9 l .. SS 10: [l01ACA 100 58 100 100 100 100 94 100 100 1.1 1.3 1.2 0.30 11: [l01AGA 100 58 100 100 100 100 94 100 100 100 1.1 1 .. 1 0.16 CONSISTENCY: 6.321E-04

METHOD 1 1: ROMBERGINT 3: SIMPSONINT 97 4: TNT. TRAPEX 56 5: INTEGRAL 98 6: QADHAT 4'') "-"7n .'

.

I NTEGRr.IAL. '7 ' l b B: [to:l. F1At':1 14 9: [101ABA 53 10: D01ACA 49 U. = DO:l.AGA 73 CONSISTENCY: 2.977E-02 "1 'J 27.0 '36 99 4') ,-77 14 53 49 73 4 ~:j 1.(114.0

o

u ~.~ 0 .. 5 3.3 :::;6 :35 42 44 'i'B 6 14 46 ~;3 36 49 48 74 6 B.5

LO

4 .. 5 2.0 T7 11 4r) "-:37 41 ('y _:LO ""

'"

8 1:1. 5.4 14.0 :I. • :,'$ ;5.7 16.::;

o

.

~1 _{.;.. l'l2 0.1 0.4 0 .. 9} 1.4 3.1 O. '7 1..5 4.7 0 .. 5 2.9 O. :.3 0.8 L? 0.2 1 • :3 0.1 0.4 0 .. 8 4.0 0.6 1.7 3 .. 9 13 0.1 0.6 0.6 47 11 2.6 5.6 4':' •• J 1.4 41 2.0 62 14 49 46 REL.EFF.NUMBERS

MEAN EIG. LOG.

ROW VEC. SCAD

SUM 1<J.~~ 17.:\. 4.10 :1..1 LO 0.00 t:- 1::-oJ tt • .J 505 2.47 2 .. 3 2.0 1.02 1.2 L2 () .. 22 4.6 4 .. :1. 2 u O~j 1.0 1..:1. 0.10 7 _{. J} t:' 7.7 2.94 2.5 2.6 1.3(1 1..4 1 • :;,~ 0.29 UJ LI'I

MEJiN EIG. LOG.

ROW VEC. SCAR

METHOD 1 3 4 ~) ₆ 7 8 <1 10 11 SUM

---1: FWMBERGINT LL<1 ,f") P) 7.0 30.4 12.5 19.9 :L7 (/.4 30 .. 6 2~j .. 6 24. :5 4.60 ~ • ., 4,... 3: SIMPSONINT f.1O

o

.

'7 _{..; 0.5} '"j -y &.. . . . :)

o •

<y 1.7 0.3 1.3 -) A- ... '"'1 ' 2ltl '1 ,,;;., .. A .... '') 1.11 4: TNT. TRAPEX 51 51 1.6 7.1 3.0 7.0 0 .. 8 3.2 6. <1 6.B 7.3 2.86 5: INTEGRAL 80 100 51 4.1 1.7 4.4 0.6 r) , ..:,."0 4.3 4.

:3

4.3 2.0<1 6: QADRAT 80 100 51 100 0.4 1.2 () • 1

o

.

'7 , _{1.0 1.1 1..0 0.05 w}

'"

'?~ INTEGf-\:AAL. 80 98 51 98 98 2.9 0.3 l.6 2.6 :::: h C) ...., ,::-.-:.~ n ... J :1. .. 33 8: D01AAA 9 9 9 9 (1 9 0.1 0.9 1.0 1.0 1.0 0.05 <1: D01ABA 38 38 38 38 38 :58 8 5.6 7.5 8.9 9.3 3 .. 22 10: DOIACA :!~3 29 '') 1 ';'.,t. 29 29 29 9 21 1.6 1.5 1. '7 0.76 11.: D01AG~1 79 99 :'51 99 99 97 9 3D 29 :1,.0 1..0 0.00 CONSISTENCY: 1.716E-02

...

1:~El. EFF "NUMBER~3

MEAN

EIG. U)(3.

ROW VEC. SCA.

METHOD , 3 4 J::* 6 "1 ₈ 9 10 11 SUM .L ,.} I

---1: ROMBERGINT 6.3 ~L 1 1.4 39.4 12.4 58n8 2.2 24.1 42.1 39.9 32.9 5.04 3: SIMP~3()NINT 79

o •

7' 0.4 6" 1 2.8 4.0 1.2 4.8 '7" :2 6.0 6.2 ?63 4 = I NT . TI;:APEX 44 44 0.4 7.1 3.4 LLO 0.9 6" ;2 7.4 8.6

B.l

3.01 I::' I NT EG f-t: A I... 28 3"'" :~~6 15.5 7.6 10.0 3.7 11.9 14.7 :1.4.6 :1.6.7 4.06 -.. J :: . J 6: QADF~AT 79 91 44 35 0.5

On?

0.2 0.8 :l..2

LO

LO

o .0:.5 w -.J "7 II INTEGRA(-1L. 51 54 42 30 ~54 2 .. 0 0.4 1.7 2., :~ ;? .. :~~ r~ r) 1.16 l • _{.... " 4} 8: DO:lAAA 7 7 I". 7 7 7 0.1 1.0

LO

:I. • 1 1.1 O.OS -..}

9: DOl ABA 13 13 13 1. ~~ 13 12 :1. 4.3 ~:i .. 0 5.9 6.n 2.77

10: DOl.ACA 14 15 12 14 1:5 le:-. ..1 ...•

"

'7 1.l 1.2 1."2 0.32

11: D01AGA 76 86 44 3~5 El6 54 /' :1.3 15 1.0 :1..0 0.00

MEAN I::.IG. LOG.

ROW VET. SCAN

METHOD 1 2 "

..

""J

'*

r.:-' .•. } 6 "? 8 9 1.0 1.1 SUM

---1.: ROMBEI:;:GINT 2.6 _2"7 2.9 On :'i :3 • f:l 1.0 4. "7 1.4 5.0 3.8 4.4 4.6 ~? • 1 (}> 2: BULIF<SCHHH "-Z _{L....J 0.6 1.0 0.1 1.0 0.4 L l 0.5 1.1 0.9 1.2 1.3 0.41} 3: SIMF'BONINT 100 23 1.3

o ')

.

~" L6 0.4 1.8

o

• ..J 0:.:- 1.9 1.4 1.8 1.8 0.88 4: INT. TRAF'EX 72 19 72 0.1 1..2 0.3 1 • :~

o "'"

• ..J 1..3 1.l. 1.4 1.4 0.50 5: INTEGRAL 100 23 100

.,

I ... r::. , "7.9 1.8 9.5 2.7 9 _...J r.:- _{"7.2 9.5 9 .. ~:!} . . 'ill:" V .. A,' •• """ W ():) 6: WHiRAT 74 21 74 5B ""1 , _{0.3 1.0 0.4 1..1 0.9 1.1} i ' i _0.23 I I "+ .1. q 14M 7: I NTEGf-<AAL 100 2] 100 72 100 74 ~5 .. 0 1.5 5112 3.9 4.6 4.7 2.24 8: D01AAA 98 -'"Z _{98 72 98 71 98 0.3 1.0 0.8 1.0 1.0 0.04} ... ..1 9: D01ABA ('19 ₂₃ 99 71 99 72 99 97 3.5 2.6 3.:1. 3.2 1.68 10: D01ACA 100 ... _~ ) ... ..

"

} 100 ?:~ 100 74 100 9B 99 0 .. 8 :1. .. 0 1..0 0 .. 00 11: DO:l.AGA 99 23 ("19 71 99 74 99 9"7 SiB 99 L~~ 1.3 0.36 CONSISTENCY: 5.442E-03

REL.EFF.NUMBERS

MEAN EIG. LOG.

ROW VEC. SCAM

METHO[I 1 3 4 5 6 7 8 9 10 11 SUM ---~--- --- 1: ROMBERGINT 21.0 3.4 6 .. 5 16.2 7 .. 3 43.5 1.6 44 .. 0 24.2 46.1 43_0 5.43 3: SIMPSONINT 95 0.4 0 .. 3 0.8 0.3 2 .. 4 0.2 2.0 1.2 2.4 2.4 1 .. 29 4: INT.TRAPEX 55 55 1.0 2.5 1.1 8.1 0.5 7.1 4 .. 2 7.8 8.0 3 .. 01 5:: INTEGRAL 85 86 53 2.3 1.0 7.2 0 .. 6 6.2 3.6 7.1 7.3 2.87 w 6: QADRAT 94 100 55 85 0.4 3.0 0.2 2.6 1..5 3.0 3.0 1.61 1.0 7: INTEGRAAL 94 99 54 85 99 6.9 0.5 5.8 3.4 6.8 7.0 2.80 8: It01AAA 85 86 55 77 85 85 0.1 0.9 0.5 1.0 1.0 0.00 9: [l01ABA 52 52 46 51 51 51 51 12.3 7.5 13.9 15.2 3.93 10:: [101ACA 95 100 55 86 100 99 86 52 0.6 1.1 1.1 0.19 11: D01AGA 92 97 55 83 97 96 84 51 97 2.0 2.0 0.97 CONSISTENCY: 1.056E-02

REL.EFF.NUMBERS

MEAN E1G. LOG.

ROlJ VEC. SCAM

METHOD 1 3 4 5 6 7 8 9 10 11 SUM

---1: ROMBERGINT 29.4 3.4 9.5 67.4 34.4120.3 1.7131.6 70.4 118.1 97.9 6.61 3: SIMPSONINT 52 0.7 0.5 1.8 1.2 4.3 0.7 3.6 2.2 4.0 4.8 ~.26 4: INT .. TRAPEX 12 12 0.9 4.8 2.0 8.5 0.5 9.5 4 .. 1 8.3 9.9 3.31 5: INTEGRAL 29 29 11 4.5 2.2 9.0 0.9 10.1 4.8 9.0 9.4 3.24 .c::--6: QADRAT 52 95 3.2 29 0.6 1.9 0.1 1.8 1.2 1.9 2.0 0.98 0 7: INTEGRAAl 51 68 11 28 68 3.7 0.4 3.2 1.9 3 .. 5 3.8 1.91 8; D01AAA 48 48 12 29 48 47 0.1 1.1 0.6 1.0 1.0 0.00 9: D01ABA 9 9 8 8 9 9 9 12.1 5.5 11.6 16.1 4.01 10: D01ACA 52 76 12 29 76 68 48 9 0 .. 6 1.0 1.0 0.02 11: D01AGA 51 93 12 28 98 67 47 9 75 1.7 1.9 0.90 CONSISTENCY: 6.515E-02

REL.EFF.NUMBERS

MEAN EIG. LOG.

ROY VEC. SCA.

METHOD 1 2 3 4 5 6 7 8 9 10 11 SUM

---1: ROMBERGINT 1.9 3.5 1.4 0.8 2.7 0.9 3.4 1.4 1.4 2 .. 2 3N~ _3.5 1.79 * 2: BULIRSCHINT 21 0.9 0.8 0.3 0.9 0.3 1.2 0.7 0.4 0.7 1.3 1.3 0.41 3: SIMPSONINT 96 21 0.4 0.2 1..0 0.3 1.5 0.4 0.4 0.7 1.2 1.2 0.22 4: I NT • TRAF'EX 71 19 73 0.5 1.6 0.6 2.8 1.0 0.9 1.4 2.3 2.3 1.22 ~ 5: INTEGRAL 97 21 100 73 3.7 1.1 4.8 1.7 1.8 2.9 4.7 4.6 2.20

-6: QADRAT 57 17 59 43 60 0.3 1.1 0.5 0.5 0.8 1.3 1.3 0.37 7: INTEGRAAL 86 21 89 65 89 56 4.1 1.5 1.6 2 .. 6 4.1 4.0 2.02 8: D01AAA 33 11 34 26 34 25 31 0.4 0.5 0.7 1.0 1.0 0.00 9: I101ABA 93 21 93 71 93 56 84 32 1.0 1.6 2.5 2.5 1.31 10; [I01ACA 80 20 82 60 83 52 74 28 79 1.6 2.6 2.6 1.36 11: D01AGA 70 17 73 55 73 47 66 28 68 60. 1.6 1.6 0.68 CONSISTENCY: 6.630E-03

REL.EFF.NUMBERS

MEAN EIG. LOG.

ROW VEC. SCA.

METHOD 1 3 4 5 6 7 9 10 11 SUM

---1: ROHBERGINT 26.4 1.5 6.8 17.1 5.3 1 .. 8 1.6 18.1 22.1 18.4 4.20 3: SIHf'SONINT 94 0.1 0.3 0.6 0.2 0 .. 2 0.3 0.8 1.0 1.0 0.00 4: INT.TRAF'EX 43 43 2.5 4.4 1.5 1.2 1.1 5.6 7.0 7.4 2.89 5: INTEGRAL 80 84 43 1.8 0.6 0.5 0.8 2.2 2.9 2.9 1.56 .j::-6: QADRAT 91 96 42 81 0.3 0.3 0.5 1.2 1.6 1.6 0.67 N 7: INTEGRAAL 87 88 38 75 86 0.8 1.1 3.9 5.0 4.8 2.27 9: D01ABA 54 54 37 54 52 49 1.2 4.6 5.8 6.3 2.65 10: D01ACA 22 22 16 22 22 19 22 3.1 3.9 5.1 2.34 11: D01AGA 65 68 31 59 66 61 40 16 1.3 1.3 0.34 CONSISTENCY: 3.390E-02

REL .. EFF .. NUMBERS MEAN EIG. LOG ..

ROW VEC. seA.

METHOD 1 3 4 5 6 7 9 11 SUM

---1: ROMBERGINT 33.3 2.0 15.1 42.0 17.9 1.8 46.5 48.5 39 .. 4 5.30 3: SIMPSON!NT 51 0.2 0.5 1 .. 0 0.5 0.4 1.6 1.6 1.7 0.78 4: INT.TRAPEX 14 14 2.6 7.5 2.9 1.0 8.5 8.5 9.4 3.23 5: INTEGRAL 35 46 13 2.4 1.0 0.6 2.9 3.1 3.1 1.63 .j:-. 6: QADRAT 51 98 14 46 0.4 0.2 1.5 1 .. 5 1.3 0.40 w 7: INTEGRAAL 43 60 13 40 60 0 .. 4 2 .. 9 3.0 2 .. 9 1.52 9: D01ABA 12 12 9 11 12 11 5 .. 7 6 .. 0 8 .. 1 3 .. 01 11: D01AGA 39 70 11 35 70 45 10 1 .. 0 1.0 0.00 CONSISTENCY: 5.600E-02

METHOD

1 2 '7 ' • .1 4 !=) 6 7 8 9 :1.0 1.1 t:

RDMJ.3ERGINT

2.9 6.5 2.2 2.0 4.0 3.6 6.8 1.2 7.6 3.0 2:

BULIRSCHINT

10

0.4

0.9

0.2 0.6 0.4

1.0

0.4 0.6 0.3 3:

SIMF'SDNINT

81 10 4 =

INT.

TRAPEX 67 9 ~S6 5:

INTEGRAL

92 10 84

6: OADRAT

57 9 46 7:

INTEGRAAL

85 9 76 8:

D01AAA

4') ₈ ~n 9: DOtABA 73 9 59

10:

D01ACA 94 10 f:\4 11:

D01AGA

'77 9 67

CONSISTENCY: 3.481[-02

0.6 0.3 1.4 0.5 2.4 0.3 1.2 0.7 0.6 1.6 1.2 2.9 0.5 2.0 1.2 66 3.9 1.7 6.8 1.0 3.8 2.1 49 60 37 61 6tl 63 c·.!:' ;;;J "I 87 4:1 71 96 77 0.5 1.7 0.3 1.0 0.5 51 3.9 0.5 2.2 1.1 35 i~'r ... ,.J 57 ~:;4 40 64 89 72 0.2 0.6 0.3 41 3.9 2.0 42 73 0.5 41 67 79

REL..EFF.NUMBERS

MEAN

EIG.

LOG.

ROW

VEC.

SCAR

SliM

8,,3 " '") .I, ... 2 n ~3 :;.~., 9 6.4 1 .7 3 .. 5 1.0 5,,6 1.7 3,,4 9·,0 3.16 1" 5 0" ~;f.l 2 .. 1 1.10 3.2 L6'7 6.1 2.60 1.7 0 .. 79 3 .. 4 1.15 LO 0 .. 00 5"9 2 .. ~:j5 1.6 0.71 3.2 1. bB ~ ~ [

REL.FFF.NUMBERS

MEPIN EIG. LOG.

ROW

VFC.

SCAD METHOD 1 2 3 4 5 6 7 8 9

to

11 SUM

---1: ROMBERGINT 4.7 3.1 2.9 L 1 10 .. 9 6.2 10.5 1«1 13.3 6.5 14.0 13 .. 8 3.79 2: BULIRSCHINT 29 0.3 0.7 0.1 :1..0 0.6 1.7 0 .. 2 1.4 0»

el

t ,.8 :1..8 0.88 3: SIMPSONINT 79 28 1 ..,

.

,

0.4 3.2 2.0 6.2 0.9 4.1 '") '") ..:.." ,-;,.. 5" <y 5.7 '") 4· . . . . .... '::-1 J J, 4: INT.TRAPEX 64 27 ~j'7 0.2

..,

.... " 4

"

1.3 3.7 0 .. 4 2.9 1.3 :3 It 6 3.5 1-82 ~ V1 5: INTEGRAL 85 29 '78 64 9.'7 5.6 :1.4.9 2.0 :1.2.2 5.8 l.~j", 4 14,,8 3.89 6: QADRAT 87 29 93 64 85 0.6 1.6 0.2 1.2 0.6 1.6 :1..6 0.67 7: INTEGRAAL 87 29 84 64 85 91 2.9 0.4 2.2 1.1 2 .. 9 2 .. 8 1.48 8: DOl.AAA 63 29 58 59 63 64 63 0.1 0.7 0.4

LO

:l .0 0 .. 00

9: n01ABA ei4 ~~8 47

t::'"

,.J4 54 ~.:i4 54 ~~;2 5.8 2~5 7.5 7 .. 7 2.94

10: DOlf:ICA 86 29 86 64 85 94 88 64 54 0.5 1.3 :1..3 0.34

11.: DOlAGA 8'7 29 82 64 P"" .Id 9() 138 1.;3 54 B8

"

~_ II . .r

.,

2,,6 1.:W

ME?IN ETG. LOG. RDW

VEe.

SCAn METHOD '1 _{... 2 3 4 5 6 7 8 9 10 11 SUM}

---1: RClMBERGINT 6 .. 0 1.2 3.9 0.3 :1.2.1 "7.6 16.7 1.4 22.6 10 .. 0 23·.7 ,.,,o) t::" .. ·_A .. ft ...J 4.49 2: BlILIRSCHINT 33 0.1 0.7 0.1 0.8 0.5 2.3 0.2 1.6 0.7 2.3 :~ 11 4 1.26 3= SIMPSONINT /9 ~~2 5.3 0.4 9.9 6.1 25.3 2 .. 4 18.0 8.1 25n5 23.7 4. ~.;i6 4: INT.TRAPFX :::;4 30 :;3 0.1 1.8 1.1 4.4 0.4 3.3 1.4 4.4 4 .. 3 2 .. 11.

.,..

0\ 5: INTEGRAL 28 1.8 ~:.~8 27 19.5 11.4 46.9 4.1 3~~.7 1,6 .. 2 48·.3 49.5 5.6:.3 6: QADRAT 81 33 81 54 28 0.6 2.6 0.3 1.8 0.8 2.6 ''l> 110', . . . . "" J 1 .. 30 7: INTEGRAAL 79 ~52 79 S2 27 81 4.4 0.4 :5.0 1.;3 4 • :3 4.1 ;.~. 02 8: [l01AAA 50 33 50 47 2B 50 50 0.1 0.7 0.3 1.0 LO 0.00 9: nOl.ABA 37 29 37 :36 24 ~57 ~57 37 7.4 3.3 10 .. 6 1.0.9 :3.45 10: DOl.ACA 76 :53 '76 :'51 2B 76 '7/.) 50 37 0.5 1.4 1.3 0.43 11: D01AGA 80 32 81 !:;3 27 94 ~:W !':jO ~f,7 '76 3 .. ':;~ 3,,0 1 n :57 CONSISTENCY: 1.316E-02

6. Relevance of accuracy indicators

In this section we discuss the accuracy indicators we have used in the test. Some of them, like IFAIL, are described as such in the docu-mentation, other One?, like }1 in QADRAT or ERROR in INTEGRAAL, have been adopted as such by us.

In order to investigate the relevance of a certain parameter as an accuracy indicator we have established from the reliability tables (cf. pp. 23-27) for every group and tolerance the following numbers: The percentage of cases in which the parameter has given the wrong indication, i.e. column (1) + column (3);

The percentage of cases in which the indication "accuracy not obtained" is wrong, i.e. column (I)/(column (1) + column (2).

These numbers are given in two separate tables, for every considered indicator.

The parameter M of QADRAT

The procedure QADRAT uses an internally defined minimal step length HMIN. It is possible that during the computation of the integral a subinterval is encountered with length less than HMIN, for which the computed approximation of the integral does not satisfy the accuracy condition used. If so, this approximation is nevertheless accepted. The number of such subintervals is recorded in the output parameter M. We have regarded M being zero or not an indication as to whether the procedure has succeeded in computing the integral with the required

group e: = 10-3 10-6 1:)-9

o

o

o

2 30

o

4 3 33 56

44

4 51 55 46 5 38 15 20

table 1: M gives the wrong indication

Over all five problem groups and all three tolerances together the percentage of cases in which M gives the wrong indication turns out to be 27%. group E: =

o

o

o

2 29

o

44

3 100 100 100

4

92 96 94 5 47 100 100

table 2: M

¥

0 is the wrong indication

Over all five groups and all three tolerances together the percentage of problems with M

¥

0 in which the error in the computed integral is less than the tolerance turns out to be 81%.

From these results we conclude that the parameter M of QADRAT is irrelevant as an accuracy indicator.

Remark. The parameter M of INTEGRAL with the same meaning is

some-what more relevant, but in our opinion equally not usable. Particularly, the performance in the difficult problem groups 3 and 4 is rather bad.

Upon exit the parameter ERROR of INTEGRivu" contains, according to the

documentation, a realistic approximation (see footnote p. 10) of the

error in the computed integral. This value being less than the

toler-ance might be used as a criterion of acceptance of the computed

in-tegral. group £ 10-3 10-6 10-9

o

6 5 2 34 20 11 3

o

30 56 4 II 48 41 .5 10 6 6

table 3: ERROR induces the wrong indication

The overall percentage of cases In which the actual error and the

value of ERROR lie on different sides of the tolerance turns out to

be 18%.

group 2 3 q. 5

'" =10- 3 0 H'O 0 0 0

0 91 100 88 40

10-9 _{0 J6 64 50 22}

table if! ection 011 account of ERROR ~s wrong

The overall percentage of problems with ERROR> € in which the actual

error is less than € turns out to be 60%.

From these results we conclude that the value of ERROR being less than

the percentages of the cases in which the error in the computed value is less than twice the value of ERROR (table 5) and the percentages of the cases in which the ratio of the two lies between! and 2 (table 6).

~.I.0Up 2 3 4

e: = 10-3 60 60 99 63

10-6 _{54 72 99 83}

10-9 _{66 86 100 96}

table 5: actual error < 2 value of ERROR

group 1 2 3 4

e: = 10-3 19 10 I 21

10-6 _{35 26 I 8}

10-9 _{40 21 6 16}

table 6:

!

< actual error <

2

value of ERROR 5 74 97 99 5 20 26 13

From these last two tables we see that the error estimate supplied by the procedure is conservative rather than realistic.

The parameters IFAIL and ACC of DOIACA

Upon exit the procedure DOIACA indicates via the parameter IFAIL whether the computed integral is acceptable. The tables 7 and 8 reveal the

significance of this indication.

group 2 3 4 5

e:

=

10- 3 0 21 0 14 II

10-6 _{0 17 0 9 23}

10-9 _{0 8 5 0 15}

table 7: IFAIL gives the wrong indication

group 2 3 4 5

e:

=

10-3 0 29 0 20 100

10-6 _{0 19 0 8} 81

10-9 _{0 9 17 0} 37

table 8: IFAIL

=

I is the wrong indication

The overall percentage of cases with IFAIL

= ]

(accuracy not obtained) in which the error is less than the tolerance turns out to be 13%. From these two tables we conclude that IFAIL is a significant indi-cator, with a possible exception of group 5.

T!he parameter ACC of DOIACA has the same meaning as ERROR of INTEGRAAL. For this parameter the tables corresponding to the tables 5 and 6 are:

group 2 3 4 5

e:

=

10-3 100 86 100 79 90

10-6 _{89 86 100 94 79}

10-9 _{6 86 91 95 70}

table 9: actual error < 2 value of ACe

group 2 3 4 5

I::

=

10-3 0 10 31 ]6

10-6 _{6 9 5 19 22}

10-9 ₂₂ II 52 20 27

table ]0: ~ < actual error < 2 value of ACC

From these two tables we see that the value of ACC contains in most cases (convergent or not) an upper bound for the actual error, but not a good estimate.

The parameter IFAIL and ERROR of DOIAGA

The parameter IFAIL of DOIAGA has the same meaning as the same parameter of DOIACA. The table corresponding to table 7 is

group 2 3 4 5

£ == 10-3 0 16 26 5

10-6

_a

_{3 30 0}

10-9 _{0 16 2} 28 5

table 11: IFAIL gives the wrong indication

From this table and the fact that IFAIL = I and yet the integral is

-6 -9

computed correctly, occurs only twice (group 4, £

=

10 ,10 ) we conclude that IFAIL ~s a significant indicator.

Upon exit the parameter ERROR contains, according to the documentation, an estimate of the absolute error in the computed integral. For this parameter the tables corresponding to the tables 5 and 6 are:

group 2 3 4 5

c: 10-3 lOa 73 90 48 79

10-6 _{100 99 97 47 90}

10-9 _{100 100 99 58 96}

table 12: actual error value of ERROR < 2

group IS

=

10-3 5 10-6 10-9 ₂ table 13:

!

< 2 3 4 7 17 15 10 17 14 1 1 12 actual error value of ERROR 5 2 3 < 2

From these two tables we see that ERROR contains nearly always an upper bound for the actual error, as the documentation states.

The results with respect to group 4 for both IFAIL and ERROR are re-markably bad (cf. Remarks on group 4, page 30).

7. The Robinson test set

In [11J Robinson gives a test set of integrals divided in groups that are more or less comparable with the groups on page 13, namely:

A - Well-behaved functions;

B - Functions with one or more integrable singularities within the integration interval [a,bJ;

C - Functions with one or more sharp peaks 1n [a,bJ, or with a singu-larity "near" [a,bJ;

D - Rapidly oscillating functions;

E - Step functions; functions with a number of finite discontinuities in the function or its first derivative within [a,bJ.

A test performed with Robinson's set may give an idea of the validity of the results and the ensuing conclusions in section 5. In this section

we will discuss the confrontation of our results with the results ob-tained by test runs with the types A,B,C and D of Robinson's set.

Remark. In contrast to our test set the functions 1n Robinson's set are not made positive, but the integrals as given in [IIJ are computed.

Type A

This set of integrals agrees with group 1. The reliability score· agrees practically with that of group 1, though most of the pro-cedures do not have, 100% correctly computed integrals.

The differences in the efficiency are small. The corresponding entries of the log. scale columns of the efficiency tables differ about 0.2 - 0.3, except INTEGRAL and DOIAGA where the difference is about 0.6. If the function A2I: (x - 1)(x - 2) ••• (x - n)/n:, which might be considered oscillatory rather than smooth particularly for the larger values of n, is omitted, then the difference for INTEGRAL and DOIAGA reduces to the level of the other procedures.

Type B

The integrals of type B may be divided into the following three groups. i) Unbounded functions: BI, B2, B8, BII, B17-22. These functions

be-long to none of the five groups.

All procedures, except SIMPSONINT sometimes, fail in computing the correct value of these integrals, in agreement with our earlier