NUMERICAL METHODS IN MULTIPLE INTEGRATION

Dissertation for the Degree of Ph. D.
MICHIGAN STATE UNIVERSITY
WAYNE EUGENE HOOVER
1977

LIBRARY

Michigan State
University

 

This is to certify that the

thesis entitled

NUMERICAL METHODS IN MULTIPLE INTEGRATION
presented by

Wayne Eugene Hoover

has been accepted towards fulfillment
of the requirements for

PI; . D. degree in Mathematics

 

n fl “
w A $11 W
Major professor

I (J.S. Frame)

 

 

Dag Octdber 29, 1976

 

 

 

. —-

., ,_.~_
(‘sq

I MAY 2T3 2002

 

ABSTRACT

NUMERICAL METHODS IN MULTIPLE INTEGRATION

By

Wayne Eugene Hoover

To approximate the definite integral
1:" bl
[(1):] f(x1"”’xn)dxl”'dxn
an a‘
over the n-rectangle,
n
R = n [41, bi] ,
1' =1

conventional multidimensional quadrature formulas employ a weighted sum of function values

m
Q“) = Z ij(x,'1:""xjn)-
i=1

Since very little is known concerning formulas which make use of partial derivative data, the objective
of this investigation is to construct formulas involving not only the traditional weighted sum of function
values but also partial derivative correction terms with weights of equal magnitude and alternate signs at
the corners or at the midpoints of the sides of the domain of integration, R, so that when the rule is
compounded or repeated, the weights cancel except on the boundary.

For a single integral, the derivative correction terms are evaluated only at the end points of the
interval of integration. In higher dimensions, the situation is somewhat more complicated since as the
dimension increases the boundary becomes more complex. Indeed, in higher dimensions, most of the
volume of the n-rectande lies near the boundary. This is accounted for by the construction of multi-

dirnendonal integration formulas with boundary partial derivative correction terms, the number of which

increases as the dimension increases.

Wayne Eugene Hoover

The Euler-Maclaurin Summation formula is used to obtain new integration formulas including a
derivative corrected midpoint rule and a derivative corrected Romberg quadrature. Several new open
formulas with Euler-Maclaurin type asymptotic expansions are presented.

The identification and utilization of the inclusion property, the persistence of form property, and
the equal weight-altemate sign property coupled with the method of undetermined coefficients provide
the basis for the derivation of a number of new multidimensional quadrature formulas.

These new formulas are compared with conventional rules such as Gauss’ and Simpson’s rules
and the numerical results show the derivative corrected formulas to be more efficient and economical

than conventional integration rules.

NUMERICAL METHODS IN MULTIPLE INTEGRATION

Wayne Eugene Hoover

A DISSERTATION

Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of

DOCTOR OF PHILOSOPHY

Department of Mathematics

1977

(:> Copyfightby
WAYNEEUGENEHOOVER

1977

to
my wife Diane
and

daughter Susan

ACKNOWLEDGMENTS

It is a great pleasure to express my sincere appreciation to Professor J. Sutherland Frame, Michigan
State University, East Lansing, for suggesting the problem. His invaluable guidance, expert advice, stim-
ulating discussions, and useful insights made this work possible.

[would also like to thank Professors Bang-Yen Chen, Edward A. Nordhaus, and Mary Winter for
serving on my guidance committee.

For imparting to me a generous measure of their enthusiasm for and wide knowledge of the field
of numerical analysis, it is a pleasure to thank Professors E. Ward Cheney, J. Sutherland Frame, Gfinter
Meinardus, and Gerald D. Taylor.

The financial support for this work was provided by the U.S. Naval Air Test Center (NATC),
Patuxent River, Md. For this I am indebted to Mr. James C. Raley, Jr., NATC Staff, Mr. Samuel C. Brown,
Computer Services Directorate, and Mr. Theodore W. White, Employee Development Division.

I wish to thank Mr. Howard 0. Norfolk and Mr. J. William Rymer, Technical Support Directorate,
for making available the computing facilities of the Real Time Telemetry Processing System. I also wish
to thank Mr’. Larry E. McFarling for providing the excellent 3-D plots, and Mr. Durwood Murray for
capable programming assistance.

NATC also provided for the composition and reproduction of this dissertation. The typesetting
was done by Engineering Services and Publications, Inc., Gaithersburg, MD, and the printing by the
Technical Information Department, NATC.

Finally, I thank my wife, Diane, and daughter, Susan, for their confidence, patience, and under-

standing displayed throughout the duration of this investigation.

iv

TABLE OF CONTENTS

LIST OF TABLES
LIST OF FIGURES
KEY TO SYMBOLS
INTRODUCTION

FUNCTIONS OF ONE VARIABLE

2.1 Bernoulli Polynomials and Numbers

2.2 The Euler-Maclaurin Summation Formula
2.3 Error Estimates

2.4 A Numerical Example

APPLICATIONS OF THE EULER-MACLAURIN SUMMATION FORMULA

3.1 Quadrature Formulas with Asymptotic Expansions
3.2 Error Estimates

3.3 Newton-Cotes Formulas

3.4 The Midpoint Rule and Some Open Formulas

3.5 Romberg Quadrature

3.6 Derivative Corrected Romberg Quadrature

3.7 A Numerical Example

FUNCTIONS OF TWO VARIABLES

4.1 The Euler-Maclaurin Summation Formula
4.2 Error Estimates

4.3 Additional Error Estimates

4.4 Sharper Error Estimates

4.5 A Numerical Example

APPLICATIONS OF THE 2-DIMENSIONAL EULER-MACLAURIN
SUMMATION FORMULA

5.1 Cubature Formulas with Asymptotic Expansions

5.2 Some Cubature Formulas Obtained from the Euler-Maclaurin
Summation Formula

5.3 The Midpoint Rule and Various Other Formulas

5.4 A Numerical Example

Page

viii

M

NMO‘M

l4

14
21
22
22
31
33
34

38

38
42
45
45
48

54
54

57
92
94

MULTIDIMENSIONAL QUADRATURE FORMULAS WITH PARTIAL
DERIVATIVE CORRECTION TERMS

6.1
6.2
6.3
6.4

The Equal Weight-Alternate Sign Property

Derivation of MINTOV

Comparison of Several Multidimensional Quadrature Formulas of Precision Five
Construction of 47 New Cubature Formulas with Partial Derivative Correction
Terms and Error Estimates

NUMERICAL RESULTS

7.1

7.2

Double Integrals and the DC-MQUAD Algorithm

7.1.1 The Application of MINTOV

7.1.2 MINT OV vs JPL’s MQUAD

7.1.3 Error Estimates

7.1.4 A Rational Function with a Singularity Approaching the Domain
of Integration

7.1.5 The Comparison of 45 New Cubature Formulas with 12
Conventional Rules

Triple Integrals

7.2.1 A Function with Vanishing Mixed Higher-Order Partial Derivatives
7.2.2 MINTOV vs JPL’s MQUAD
7.2.3 An Example Illustrating the Persistence of Form

CONCLUSIONS AND RECOMMENDATIONS

LIST OF REFERENCES

General References

Bibliography

Page

96

96
99
106

110

118

118

121
126
127

136

137
141

142
146
147

150

152

157

Table

2.1.1
2.4.1
3.1.1
3.1.2
3.1.3
3.1.4
3.4.1
3.4.2
3.4.3
3.4.4
3.4.5
3.7.1
3.7.2
3.7.3
3.7.4
3.7.5
3.7.6
4.5.1
4.5.2

4.5.3
4.5.4
4.5.5
4.5.6
4.5.7
5.4.1
6.1.1
6.1.2
6.1.3
6.3.1
6.3.2

6.3.3
6.3.4

LIST OF TABLES

Bernoulli Numbers

Computation of n

Derivation of Closed Quadrature Formulas

Quadrature Weights and Derivative Correction Terms

Principal Errors and General Terms in the Asymptotic Expansions
Leading Terms in the Asymptotic Expansions

The Midpoint, DC Midpoint, and Simpson’s Rules Applied to I36 — = 1n 2
Derivation of Open Quadrature Formulas

Chradrature Weights and Derivative Correction Terms

Principal Errors and General Terms in the Asymptotic Expansions
Leading Terms' in the Asymptotic Expansions

Romberg Quadrature T-table for ID 1% sin (1rx)1rdx - 1

Derivative Corrected Romberg Quadrature (3 — l) Cl-table
Derivative Corrected Romberg Quadrature (s - 2) C 2 -table
Romberg Error for fol V2 sin (1rx)1rdx = 1

Derivative Corrected Romberg Error (3 = l)

Derivative Corrected Romberg Error (s = 2)

Partial Derivative Correction Sums @(h, h; 23, 23)

Euler-Maclaurin Summation Formula: [(1) = 0.523 248 144

z T(h,h) - <I>(h, h; 23, 23)

ErrorR(h, h; 23, 23) =I(f) - T(h,h) + <I>(h,h; 23, 23)

Estimates for R(h,h; 23, 23) using (4.3.4)

Estimates for R(h,h; 23, 23) using (4.3.5)

Estimates for R(h,h; 23, 23) using (4.3.6)

|(Error Estimate 4.3.k)/Error|, k = 4,5,6

Several Cubature Rules Applied to [01101 exydxdy = 1.317 902‘151 454 4
Elements

Generators

Generator Values

Number of Function Evaluations for Several Fifth-Order Formulas
MINTOV: Number of Function Evaluations Required for n Subdivisions
in d-Dimensions

Maximum Number of Subdivisions n such that nfe < 106
Maximum Usable Dimension d Assuming n > 4 and nfe < 10“

Page

l3
l7
I8
19
20
24
27
28
29
30
35
35
35
36
36
36
51

51
51
52
52
52
53
94
97
98
98
108

108
109
109

Table

6.4.1
7.1.1

7.1.2.1

7.1.3.1

7.1.3.2

7.1.3.3

7.1.3.4

7.1.4.1

7.1.5.1

7.2.1

7.2.1.1

7.2.2.1

7.2.3.1

Cubature Formulas
Nfe for Various Cubature Formulas Compounded nm Times. These are
Fifth-Order Formulas Except for Simpson’s Rule Which is Third-Order

1011010 +x2y2)-1dxdy = 0.915 965 594 177 30

MINTOV vs MQUAD for the Integral 112" f 2"

Relative Error Requested = 10'“, a = 1(1)10.
Error Estimates for Approximating (if: m dxdy =
6.859 942 640 334 65 when h = k = 2 (n = m =1).
Error Estimates for Approximating [id \/3—+x—+-y dxdy =
6.859 942 640 334 65 when h = k = 1/5 (n = m =10).
{:13 m dxdy = 6.859 942 640 334 65. Guaranteed MINTOV Error =
10‘“, a =1(1)12. (Grid Size )1 #5 k.)
[if-i m dxdy = 6.859 942 640 334 65. Guaranteed MINTOV Error =
10’“, a =1(1)12. (Grid Size h = k.)
Application of Various Cubature Formulas to I: fl [4(w + 2 + x +y)]’1dxdy
withh=k= l/50 (n=m = 100).
The Comparison of 45 New Orbature Formulas with 12 Conventional Rules for the
Double Integral fol I01 me" + 1) sin (1ry)dxdy = e/1r
Nfe for Several Multiple Quadrature Formulas Repeated n ln2n3 Times
flzflzflzln(xyz)dxdydz = 1.158 883 083 359 67
MINTOV vs MQUAD for the Integral f-ZI: 1:3: [:2 cos (x) cos (y) cos (z)dxdydz =
8. Relative Errors Requested = 10““, a = 1(1)10.

"/2 "/2 "/2 -l—+!- sin (x) sin (y) sin (2) e'wdxdydz =

f0 f0 f0 xyz
1.531670 226 93, w2 =x2 +y2 + 22.

viii

1 (xy)-1dxdy = 0.550 471 023 504 079.

Page
1 12

119
124

127

131

132

133

134

137

139

142

143

146

148

Figure

2.3.1
4.1.1
4.1.2
4.5.1
5.2.1
5.2.2
5.2.3
5.2.4
5.2.5
5.2.6
5.2.7
5.2.8
5.2.9
5.2.10
5.2.11
5.2.12
5.2.13
5.2.14
5.2.15
5.2.16
5.2.17
5.2.18
5.2.19
5.2.20
5.3.1
5.3.2
5.3.3
5.3.4
5.4.1
6.2.1
6.2.2

LIST OF FIGURES

Values of h and 3 for the Best Error Estimate
Trapezoidal Weights

Partial Derivative Weight Assignments

Graph ofz = (x +y + 1)‘1 on [0,1]2

Trapezoidal Rule

DC Trapezoidal Rule

Simpson’s Rule

Component Trapezoidal Sums for Simpson’s Rule
DC Simpson’s Rule

Simpson’s Second Rule

Component Trapezoidal Sums for Simpson’s Second Rule
DC Simpson’s Second Rule

52 Point Rule

DC 52 Rule

Boole’s Rule

Component Trapezoidal Sums for Boole’s Rule
DC Boole’s Rule

Weddle’s Rule

DC Weddle’s Rule

Newton-Cotes’ 72 Rule

Component Trapezoidal Sums for Newton-Cotes’ 72 Rule
DC Newton-Cotes’ 72 Rule

Romberg’s 92-Point Rule

Component Trapezoidal Sums for Romberg’s 92 Rule
Squire’s Rule

Ewing’s Rule

Tyler’s Rule

Miller’s Rule

Graph ofz = e” on [0,1]2

Sign Arrangement for DI (f)

Sign Arrangement for 0120)

ix

Page

10
39
40
50
58
59

60
61
62
63

65
66
67
68
69
70
71
72
73
79
81
82
93
93
93
93
95
100
100

Figure

6.4.1

7.1.1.1
7.1.1.2
7.1.1.3

7.1.2.1
7.1.2.2
7.1.3.1
7.1.3.2
7.1.4.1
7.1.5.1
7.2.1.1
7.2.1.2
7.2.2.1

7.2.3.1

7.2.3.2

Node Arrangements
Graph ofz = (1 +x2y2)'1 on [0,1] 2

f(x.y) = (1 +x2y2)’l
Error Curves in Approximating fol f 01 (l + x2y2)'ldxdy

Performance Comparison in Approximating [12'1 f12.1(xy)_1dxdy
Performance Comparison in Approximating [12.1 f12'1(xy)'1dxdy
Graph ofz =W on [-1,1]2

Error Curves in Approximating f. i 1:: m dxdy

Graph ofz = %(2 +x +y)‘1 on (-1, 1]2

Graph ofz = %(e" + 1) sin (ny) on [0,1]2

Graph ofz = ln(xy) on [1,2] 2

Error Curves in Approximating flz 1'12 {12 1n (xyz)dxdydz

Graph of z = cos (3:) cos (y) on [-7r/2, 1r/2] 2

+5277:
’0’

Graph ofz = 1 sin (x) sin (y) e‘ x + y

on [—7r/2, 1r/2] 2

Error Curves

Page

116
122
123
125

128
129
130
135
136
138
143
144
146

147
149

Symbol

Ac"
al"'“s+l
A 231.1. ..., 253-1(h)
a'
”i
301)
3'0!)
B

19a(t)

i1":
Iris

KEY TO SYMBOLS

A (a1 ’al)(a2, a1) -~-(as,a1)(62562) ... (as,as) (h,k)

(2511'1'2512‘1)(2521'1’2522'1)'"(23n1‘1’23n2'1)
Cubature element
Cubature element
Cubature element
Cubature element
Cubature element
Cubature element

Quadrature formula based on the Trapezoidal rule

Lower limit of integration

Lower limit of integration, i = 1(1)N
Boole’s rule

DC Boole’s rule

a—th Bernoulli number

Bernoulli polynomial of degree a
Upper limit of integration

Upper limit of integration, i= l(l)N
a-th modified Bernoulli number
Midpoint rule

DC Midpoint rule

Midpoint rule

Function space

Function space

Function space

DC Trapezoidal Sum

(m + l, k+ 1) entry in the DC Romberg C‘-table

5" I 5,- 2 8' 8
A1112Ia11 a}: + “1'2 all

Page

54
97
97
97
97
97
97

14

102
l7
17

102

22
23
92

38
42
33

34

55

13,0)

Fhk("“)

Vector (cl , - ° ' , c”) = vertex ofH or R
Real constant

i-th component of c = #hi, (1,, or bi

Derivative correction term, 2 01(c) 91.1 f (c)
C
Derivative correction term, 262 ol-k(c) 911.1 Q) 11‘ f (c)

Derivative correction term, 91.1 [f (v(b,)) - flv(ai))]

Derivative correction term, 9 1191.1 [f(v (“1" 0,,» - f (v(ai, b k»
- f(v(b,. a,» + f(v(b,-. bk))l

Derivative correction term, fl“)(b) - f1“)(a)

Partial derivative correction term

Derivative Corrected

Dimension, number of variables

Partial derivative correction term

a-th partial derivative off with reSpect to the i-th variable
Truncation error

hzaD 2 a-l

Partial derivative correction term

m-l n-l
ha+1ke+1 Z Z f“5(x,-+th,y,-+uk)
i=0 i=0
n-l
Moving average h 2 f (x,- + th)
i=0
m-l n-l
Z f(x,- 1’ ’th’j ‘1' l4k)
i=0 i=0

Moving average hk

(b - a)(d - e)[h"kl'lir"{2 + flk’Mfliz]
Cubature element
Cubature element
Cubature element
Cubature element
Cubature element
Cubature element
Cubature element
Cubature element

real valued function

Page

99
14
99

99

99

104

40
16
106
39

96
101

41

41

38
111

116
116
110
110
110
111
110
111

Symbol

f0!)

f9?
fe

:-

”145‘ .3‘ 3‘

1(f)

1(f)

[If]

1(f)

I“

2

a3

aBIII7
Mjk...1

a—th derivative off, Olaf (x)

Partial derivative off, W39? f (x, y)
Function evaluation(s)

Cubature generator
Step size, (b - a)/n

Hypervolume ofH,2Nh1x hz x ”'x hN

Hypervolume of subregion,h1 x hz x - ' - x hN
Limit of integration

length of subinterval, (bi - ai)/n]-

Set of positive integers

Definite integral

ff(x)dx

f f bftx.y)dxdy

bN b1
J‘ f(xl"”’xN)dxl’.H’dx/V
”N ‘1

Definite integral
index

x/I

index

Step size, (d-c)/m
Definite integral
Real constant

max If“)(x)l

a<x<b

max lafix. |
(x,y)eRf. ( y)

max WI ”-959.“ x
xeR ’ k If( )I '

Number of partitions of [c, (1]
Vector (ml, ' - ' ,mN)

Page

39
13
97

103

103

96
103
11

38

\J

12
38
38
41
11

42

101

38
102

nfe

QU)

Qafi
a1 mas-+1
02,1-l,...,2,s-1(h)

R

R

R(h, a)

R(h.k; 23,29
IV!)

S (h)
S '(h)

T(h)

T(h, k)

Midpoint, (”i + bj)/2, j = 1(1)N
Dimension, number of variables
Newton Cotes’ 7-p0int rule

DC Newton Cotes’ 7-point rule

Number of partitions of la, b]
Number of subdivisions of R

Step size, 0’1" i)/hi’ j = 1(l)d
Number of function evaluations
Multidimensional quadrature formula

Partial derivative correction terms
Quadrature formula based on the Midpoint rule

Rectangle [a, b] x [c, d]
N

N-rectangle, I] [ai, bi]
i=1

Remainder integral, F1110)Ba(t)dt

Remainder integral

:21]!1+ R (h, a)
Simpson’s rule
DC Simpson’s rule
I!
Trapezoidal sum, h Z f (xi)
1' =0

m n
Trapezoidal sum, hk Z" Z "flxi' yf)
i=0 i=0

Trapezoidal sum, T(ah)

Trapezoidal sum, Ail-[T(aih, a7-k) + T(aih, a,k)]
Trapezoidal sum

(m + 1, k + 1) entry in Romberg T-table
Simpson’s Second rule

DC Simpson’s Second rule

(a1 +h1(r'1-0),'°',a~ +hN (iN -0))
Weddle’s rule

DC Weddle’s rule

xiv

Page

102
96
17
17
42

107
103
11
101
55

25

38

102

41
ll

17
17

39

17

54
32
32
17
17
103
17

17

Symbol

v(x,-)

V(Xj,xk)

0j(C)
o/k(c)

<1>(h, 3)

our, k; 23, 23)
49,,(1)

4,00

5201, k; 23, 23)
0.74

1.48

Page

(a1 +ilhl, -",xi , aj+1+ii+1hj+1""'“1v +l'NhN) 104
(“1 +1'1’11' ""xj' ”7+1 +i,-,1h,-,1, 'xk' ak+1 + ik+lhk+l' 104
“N + ilvhlv)
Romberg’s 9-Point rule 17
DC Romberg’s 9-point rule 17
wlxwzx-“wa 102
Length of j-th interval, bi - aj 102
Quadrature node, a + ih, i = 0(1)n 7
Quadrature node, c + jk, j = O(1)m 38
index 6
index 14
Positive integer 14
Area of rectangle (b - a)(d - c) 42
172.15 2a-l 3
{1 if 5,,- = 0 54
25,-, otherwise
is an element of 10
Positive integer 54
{ Vz if i= j 54
1 otherwise
Cartesian product 1
Pi, 3.14159265358979 ~--' 12
Summation 1
{-1 if c]- = —h]- or aj 99
+1 otherwise 102
o,(c)ok(c) 99
S
Derivative correction sum, Z A2a_1 8
a =1
Partial derivative correction sum 41
Bernoulli polynomial of degree a 39
Bernoulli polynomial of degree a. 39
Truncation error 56
n2/6\/5— 42
n2/3\/5_ 42

XV

Symbol

2.17
2.66
2.95
3.00

5.31

5.93

7.17

10.61

11.17

r4/45
«flux/E
2n2/3\/§
31r4l40\/g

1r4\/_6—/45

_4_1r3_[_1_+ n2 + 71/12]
3V5 2 12¢? 1.02/21»2

 

317—2}... "2 ,1
3x/521N5— 9

27r4\/6—/45
——"4‘/6-[1 +6”2 + —_h/6 ]
45 1 -(h/21r)32

xvi

Page

42
23

46

44

48

47

I. INTRODUCTION
We wish to approximate multidimensional integrals of the form

ij b1
[(1') =f ... [(xl,...’xN)dxl...de (11)
“N 01

over the N-rectangle

N
R = II [c,-.12.]
i=1

where a,- and b, are real numbers.
The traditional methods of multidimensional numerical integration employ a weighted sum of m

function values

Qtf) = : Wif(x,-1."’,xuv)- (1-2)
i=

The w, are called weights and the (in , - - - , xi”) are called nodes. The difference
E(f) = 10') - Q0) (1.3)
is the truncation error (or error).

Let pk = pk(xl, ‘ ' ' , x”) be a polynomial of degree k in N variables. We say that the multi-
dimensional quadrature rule or formula Q(f) is of order k or has degree of precision k if for any pk,
50),) = 0, but 5(Pk+1) =/= 0 for at least one polynomial pk“. ‘

Since it is not uncommOn for numerical procedures to make use of partial derivatives, e.g. in
optimization techniques, it is surprising that except for work by Tanimoto [50] 1, Obreschkoff [38] ,
and Ionescue [23] , very little is known concerning the use of partial derivatives in nonproduct

multidimensional quadrature rules.

 

1The numbers in brackets refer to entries in the Bibliography.

1

2

Therefore, the objective of this investigation is to construct a number of new multidimensional
quadrature formulas using first- and mixed second-order partial derivatives of the integrand in addition
to function values of the integrand.

It is shown that the use of partial derivatives of the integrand evaluated on the boundary of R
increases both the efficiency and accuracy of composite multidimensional integration formulas.

The accuracy and efficiency are achieved by the proper combination of the following three
properties.

The first is the “inclusion property”. It is well known that m-point Gaussian integration rules for
the N-rectangle, R, are extremely accurate. However, when R is subdivided into 3 subrectangles or cells
and the m-point Gauss rule is applied to each cell, the total number of nodes is ms since the nodes are
interior to each cell.

A more efficient procedure is to employ an integration rule in which some nodes coincide with
the boundary of the domain of integration. Then when the domain is subdivided, these nodes are
included in more than one cell. Thus the total number of nodes is considerably less than the sum of
their numbers in each cell. We call this the “inclusion property”.

The second is known as “persistence of form”. Briefly, this means that for many functions,
it requires approximately the same if not less computer time to evaluate the partial derivative at a
point where the function is being evaluated as it does to evaluate the function at another point.

Finally, efficiency and accuracy result from applying the “equal weight-alternate sign property”.
Essentially this means that in the composite formulation of a rule, the weights of the partial derivative
correction terms cancel at interior points and consequently, the partials need be evaluated only on the
boundary of R. This results in a substantial increase in efficiency.

Numerical multiple integration is currently receiving considerable attention by numerical analysts.
The first book on the subject, written by Stroud [49] , appeared only recently. An excellent introduc-
tion may be found in Davis and Rabinowitz [13]. Comprehensive bibliographies are given by

Fritsch [l9] , Stroud [49] , and De Doncker and Piessens [14] .

3

Brief surveys of the literature may be found in Ahlin [2] , Hammer [21] , and Hammer and Wymore
[22]. A history of the Euler-Maclaurin Summation formula is given by Barnes [4] .

Squire [47] devotes an entire chapter to derivative corrected quadrature formulas. Stroud [49]
states only one cubature formula, C2:2-l, which uses partial derivatives of the integrand. Lanczos [28]
and Davis and Rabinowitz [l3] discuss quadrature rules using derivative data. Tanimoto [50] was one
of the first to consider cubature rules with partial derivative correction terms. Burnside [l l] pub-
lished the first nonproduct fifth-order cubature rule. Tyler [51] gave the first derivation of the 8-point
Burnside formula. Finally, Price [41] gives some interesting examples.

In this dissertation, Chapter 2 provides a background for the 1-dimensional Euler-Maclaurin
Summation formula (Euler [15] , Maclaurin [33] ). Several error estimates and an algorithm due to
Frame [18] are stated.

It is shown in Chapter 3 that the Euler-Maclaurin Summation formula may be used to obtain
asymptotic expansions for at least 5 of the Newton-Cotes’ quadrature formulas, the midpoint formula,
and several new open formulas including some with end derivative correction terms. The third-order
derivative corrected midpoint rule is shown to be more efficient than the classic Simpson’s rule [46].
After reviewing Romberg quadrature, we define a new technique called derivative corrected Romberg
quadrature.

The 2-dimensiona1 Euler-Maclaurin Summation formula is stated in Chapter 4. Also, several error
estimates are given.

In Chapter 5, the double Euler-Maclaurin Summation formula is used to obtain asymptotic
expansions for a variety of formulas. New asymptotic expansions are given for Squire’s [48] ,

Ewing’s [15] , Tyler’s [51] and Miller’s [35] rules.

The method of undetermined coefficients, the inclusion pr0perty, the persistence of form
property, and the equal weight-alternate sign property are employed in Chapter 6 to construct 47 new
derivative corrected cubature rules of orders 1, 3, 5, and 7. These formulas are generalizations of the
midpoint, trapezoidal, Squire’s [48] , Ewing’s [l6], Tyler’s [51] , Miller’s [35] , Simpson’s, and
Albrecht, Collatz [3] and Meister’s [34] rules. One of these, caUed MINTOV (for Multiple lNTegration,

Order 5), may be considered a generalization of Lanczos’ [28] result which Lanczos calls Simpson’s

4
rule with end corrections. Error bounds are included for the 47 new cubature rules. As is shown in the
case of MINTOV, these formulas may be generalized to multidimensional derivative corrected quadrature
rules.

In Chapter 7, we compare the first-, third-, and fifth-order formulas of Chapter 6 with existing
formulas of comparable degrees of precision. The numerical results indicate that partial derivative correc-
tion terms substantially enhance the efficiency of composite numerical integration formulas.

Conclusions and recommendations for further study are given in Chapter 8.

Lyness [29] states than an important objective of those concerned with the formulation of
numerical integration techniques is “to determine which rule requires the fewest function evaluations to
obtain a result of particular accuracy or degree”. Thus, we will pay particular attention to minimizing
the number of points at which the function or partial derivatives are to be evaluated. This will require
the generalization of the equal weight—alternate sign property to higher dimensions in order to obtain
efficient composite integration formulas with partial derivative correction terms. In cases where
derivative information is easily obtained, we will show the superiority of the derivative corrected
formulas over the traditional rules which use only function information, provided that the first- and/or
mixed second-order partial derivatives are easily evaluated.

We conclude this introduction by recalling Oliver’s [40] classic statement: “Now for any given
formula or algorithm a pathological problem can always be devised for which an arbitrary small accuracy
cannot be attained; we can therefore never argue the universality of any particular method, and we do
not attempt this.” Indeed, it is sufficient to take for the integrand, f = Mpk, where p is a polynomial
having zeros which coincide with the nodes of the integration formula, [C > 2, and M is a sufficiently

large positive constant.

2. FUNCTIONS OF ONE VARIABLE
2.1 BERNOULLI POLYNOMIALS AND NUMBERS
The Bernoulli polynomial Ba(t) of degree a has the form

3.0) = 1;) (wk)!

 

and the following properties:

Bo(t) = I

1
t =t.._
Bl() 2

Ba-1(f) = 3&0)

-..-1 bk _ lifa=l
Ba(l) — 30(0) — g0 (a—k)! - 0 if 0 >1.

(2.1.1)

(2.1.2)

The last property in (2.1.2) may be used to determine the coefficients ba. The coefficients Ba in

the expansion

a co
-3. _a-1 2a l_-__.x.£ i=2:
1 2 + :1( 1) Bax /(2a). 2 + 2 coth<2) baxa
a:

are the Bernoulli numbers and are related to the 1),, by the relation

61

(2a)! '

 

62,, = (-1)"'1 a > 0.

(2.1.3)

(2.1.4)

It can be shown that b2a+1 = O for a. > 0. Also,b2,ll =Ba(0) =Ba(l) for a > 1. For reference we

list the first 10 Bernoulli numbers in Table 2.1.1.

Table 2.1.1 Bernoulli Numbers

 

 

 

 

 

 

 

 

 

 

 

 

 

a. [ll 2 3 4 5 6 7 8 9 10
B 1 1 1 1 5 691 7 3617 43867 174611
a 6 30' If E E 2730 6 510 798 330

 

 

Adams [1] lists the first 62 Bernoulli numbers. The Unpublished Math. Tables Repository has the

most extensive list of Bernoulli numbers containing the first 836.

It is shown in Knopp [26] that

Ba(2fl)2a " I

 

, = Z - (2.1.5)
2(2a).k=1 k2a
This shows that the 3,, increase rapidly for a. > 5. Indeed
_BSLL > W ... a. (2.1.6)

Ba 2172

asa—o 00.

Frame [18] gives several interesting results concerning the Bernoulli polynomials and the Bernoulli

numbers.
Lemma 2.1.1
1 Ba
32 ___ (2.1.7)
J; a (t)dt —(2a)!
Lemma 2.1.2

Ba(2")2a ,4
(2a)! < :5 (2.1.8)

Comparing Lemmas 2.1.1 and 2.1.2 we see that

1
f B}(r)dr < 2.17(21r)'2‘1 (2.1.9)
0

where 2.17 z n4/45.

Lemma 2.1.3

1
Ibaal < 3b}, (2.1.10)
Next we will state the Euler-Maclaurin Summation formula.

2.2 THE EULER-MACLAURIN SUMMATION FORMULA
Let 025 [a,b] be the set of all real-valued functions defined on the finite closed interval [a,b]
with the property that the derivatives f (“)(x), a < 23, are continuous on [a, b] . Let 0,, denote the

following boundary derivative correction terms:

D, = f(“’(b) - f‘“)(a),ae 1’. (2.2.1)
Partition the interval [a,b] into n equal parts each of width h = (b - a)/n. Let the points of sub-
division be denoted by x, = a + ih,i = 0(l)n, where x0 =0 and x,, = b.

We now define the a-th remainder integral
1
R(h,a) = f Fig“) (t)Ba(t)dt (2.2-2)
0

in terms of the Bernoulli polynomials Ba(t) and the a-th derivative of the moving average

n—l
Fh(t) = a): f(x,+ th). (2.2.3)
i=0

The celebrated Euler-Maclaurin Summation formula expresses a sum of values of f (3:) evaluated

at the equally spaced points x i in terms of the definite integral

b
M) = I f (x)dx (2.2.4)

and a series consisting of constant multiples of the derivative correction terms 02 (1.1.

Theorem 2.2.1 (Euler [15], Maclaurin [33])

Forfe C74 [a,b]

b n .. S , 1
I f(x)dx = hZ f(X,-) - Z h2“62a02a_1 + f F(2‘)(t)82‘(t)dt. (2.2.5)
“ i=0 a=l 0

Here the double prime signifies that the first and last terms in the sum are assigned weights 1/2

and the remaining terms receive weights 1.

The proof of Theorem 2.2.1 follows by integrating (2.2.2) by parts and noting that

I: Z "f(x.)
i=0

é—[Fh(0) + Fh(1)] = 1(f) + R(h,l), (2.2-6)

T(h)

l b
I F),(t)dt = I f(x)dx, (2.2.7)
0 a

and
R(h,a-1) = bahaDa_l —R(h,a). (2.2.8)
Introducing the notation

_ 2
E2.1-1 - ’1 “020-1 (2.2.9)

A2.1-1 "’ b2aE2a-1

3
Z A2a--1 ’
a=1

we may express the Euler-Maclaurin Summation formula with remainder in the more compact form:

901.8)

10‘) = T(h) - <I>(h,3) + R(h,23). (2.2.10)

For concreteness we write the first several terms of <I>(h,s) in (2.2.10):

b n h El
I ftxldx = h f(a+ih) - 3 [f(a)+f(b)l - T2—
(1 '-0

'-
(2.2.1 1)
E3 £5 E7 139

+ —— .. + — —— +
720 30240 1209600 47900160

 

+ R(h,23).

From this we see that the Euler-Maclaurin Summation formula is a generalization of the
Trapezoidal Rule, T(h). It relates a sum of equally spaced function values and an integral and states
explicitly the boundary derivative correction terms which allow one to be converted into the other.

Consequently, it may be used for summation as well as for numerical quadrature. Also, it has

provided the basis for some useful results in the theory of asymptotic expansions. (See Oliver [39] .)

2.3 ERROR ESTIMATES

Effective application of the Euler-Maclaurin Summation formula in the approximation of sums or
integrals depends on a close estimate of the remainder integral, R(h,3), and a judicious selection of h
and 3 which will achieve the required accuracy with a minimum of computational effort.

For many functions, for example rational functions, the remainder integrals R(h,s) for a given h
may first become smaller, but then grow without bound as 3 increases to infinity.

The problem is to predict the h and s which will yield the desired accuracy without compromising

computational economy.

Applying the Cauchy-Schwarz inequality to (2.2.8) and applying Lemmas 2.1.1 and 2.1.2, Frame

[18] obtained the following upper bounds for the quadrature error in the Euler-Maclaurin Summation

formula.
Theorem 2.3.1
h S
lR(h,s)I < l.48(b—a)MS (5;) , s > 1 (2.3.1)
where
M, = Max lf(‘)(x)|, (2.3.2)
a<x<b

1.48 as n2/3\/5T

For appropriate 11 and 3, a sharper error estimate is the following.
Theorem 2.3.2

h 28+2
IR(h,23+l)| < 2.66(b-a)M2,,2 <27) , s > 1 (23.3)

where 2.66 approximates 1r4/15\/E
We observe that we now have in fact three error estimates. To see this, substitute b2,” = 0
in (2.2.14) to obtain

R0623) = -R(h.23+ 1). (23.4)

Then using (2.3.4) in (2.3.1) we may write the estimates as follows:

23
lR(h.2s)I < 1.48(b-a)M2, (5"?) (2.3.5)
I! 28+I
lR(h,23+ l)| < L48(b-a)M2,,l (3;) (23.6)
h 23+2
lR(h,23+ l)! < 2.66(b -a)M2,.,2 (27-) (2,3,7)

The selection of the apprOpriate estimate depends on the step size, h, the number of derivative

correction terms, 3, and the modulus of certain derivatives of f (3:) over the domain of integration.

10

For example, for rational integrands, the sharpest error estimate may depend on h and s as

roughly indicated in Figure 2.3.1.

 

 

h S l 2 3. 4
l
1 / 2
1/ 3 (2.3.7) (2 .3 .6) (2.3 .5)
1/4

 

 

 

 

Figure 2.3.1 Values of h and s for the Best Error
Estimate
Of considerable interest is the question for which functions f and for which values of h and a
does R(h,a) converge to zero. We are assuming f e C “[a,b] ; hence f is Riemann integrable on [a,b]

and we have

lim R(h,a) = 0. (23.8)
-hfl

On the other hand, many functions f e C,°°[a,b] , e.g. rational functions, satisfy

limR(h,a) = 9°. (23.9)

a»-

However, if we assume f e C'[a,b] and if there are constants M and c such that

Ma = max |f(“)(x)|<Mc", (2.3.10)
a<x<b

then for O < .h < 21r/c we have

lim R(h,a) = 0. (23.11)

a .. a
Examples of functions for which (2.3.11) apply are products of functions such as exponential
functions ek", trigonometric functions sin (10:), cos (10:), and Bessel functions J", (10:).

Finally, in the case f (x) e C"[a,b] is a rational function, f (x) is a sum of partial fractions of the

11

form ai/(x, - x) and f (“)(x) is a sum of terms each of the form a! ai/(xj - x)“+l. Moreover, there
exist constants M and c such that Ma. satisfies a weaker inequality

Ma. < a! MC“. (2.3.12)
Thereforeffor 0 < h < 21r/c, Frame [18] has shown that

(1-3.0.)
r(h) = min R(h,a) < (b-a)M(hc)"’51o llhc

a e] + (2 .3 .13)
and suggests the following algorithm for selecting a priori the step size, h, and s, the number of deriva-
tive correction terms, for use in the Euler-Maclaurin Summation formula:

(i) Given f (x), a, b, and an error requirement E: > 0, calculate M and c according to

Ma < a! Mo“. (2.3.14)

(ii) Choose h sufficiently small so that the error satisfies

r(h) < (b -a)M(hc)'V110(l' H39'75) < e- (2.3.15)
(iii) Determine 3 6 1+ by
3 = [n/hc] , (2.3.16)
the greatest integer < 1r/hc.
(iv) Apply the EuleroMaclaurin Summation formula to approximate a definite integral or a sum
using the above values of h and s. 0
Thus, even though inn. R(h,a) = 0°, that is, the asymptotic series diverges, it may be useful in
certain applications.
The problem of a priori estimating the h and s which not only guarantee a stated error require-
ment but also result in maximum computational efficiency by minimizing the number of function

and/or derivative evaluations (nfe) is not completely solved. It may be stated as follows.

Given f(x) 6 C“[a,b] , and e: > 0, find n = (b - a)/h and s < a/2 which minimize

nfe 23 + n + 1 (2.3.17)

subject to

h 2(S+I)
lR(h,23+ l)I < 2.66(b"a) ‘5; M2$+2 < 5. (2.3.18)

12

2.4 A NUMERICAL EXAMPLE
Consider the definite integral

‘ M):

= 11. (2.4.1)
0 1+):2

 

1(f) =

Expand the integral into partial fractions and differentiate a times to obtain

f‘“’(x) = a! 2110-10“ -(-i-x)’l'“]- (242)
Since the maximum occurs at x = 0, we have
Ma < 4a! (2.43)

Now set M = 4 and c = 1 in (2.3.15) to obtain

_1_o_
r(h) < 4h"*510I1 11”). (2.4.4)

Then using (2.4.4) and (2.3.16) we compute the values in Table 2.4.1.

The values of D2a_1 are calculated as follows.

02.1-1 f(2“'l’(l) - f‘2“'”(0) (2.4.5)

(2a - 1)!22’asin (arr/2).

Hence, employing the Euler-Maclaurin Summation formula, we have

n=hin 4

 

(s+l)/2
- Z (—1)“h4a'2b4a-2(4a - 3)! 23‘2“ + R(h, 23)

a=l

 

n 2 6 10 14
~ 4 h h
h Z: + £6 ' 5%4' 1056 ' 384
M 1 + (1702
» (2.4.6)
1118 43 867 _h22 854 513
1 838 592 1 554 432

1126 8 553103 _ 1:30 8 615 841276 005
319 488 3 519 774 720

 

s-l

+.--+ (-1) 2 rushes—nu“ + R(h,23), s = 1,3,5,-~

12
2.4 A NUMERICAL EXAMPLE
Consider the definite integral

l4d):

01+):2

 

1(f) =

=17.

Expand the integral into partial fractions and differentiate a times to obtain

f(“)(x) = a! 2i[(i-x)'1'“ —(-1—x)’1'a].
Since the maximum occurs at x = O, we have
Ma < 40.!

Now set M = 4 and c = l in (2.3.15) to obtain

(22)
1'0!) < 4h'V‘10 “ "

Then using (2.4.4) and (2.3.16) we compute the values in Table 2.4.1.

The values of DZa-l are calculated as follows.

DZa-l fan-”(1) _ f(2a-l)(0)

(2a — l)!22'“sin (arr/2).

Hence, employing the Euler-Maclaurin Summation formula, we have

 

(rm/2
- Z (-1)ah4a'2b4a_2(4a - 3)! 23'2‘1 +R(h,23)

a=l

" .. 2 6 10 14
hZ 4 + L - _h_ + ——h 6 -——2
i=0 [+0702 6 504 105 84

 

+ h18 43 867 _h22 854 513
1 838 592 1 554 432

 

h26 8 553103 _ h30 8 615 841276 005
319 488 3 519 774 720

 

s-l

+---+ (-1) 2 h2‘b23(23-1)!22“ + R(h,23), 3 =1,3,5,---.

(2.4.1)

(2 .4 .2)

(2.43)

(2.4.4)

(2.4.5)

(2 .4 .6)

13

Computation to 14 decimals with h = 1/5 and s = 15 in (2.4.6) gives

1 1 100 100 100 100
- —,17 -— +—+——+—+——+
" R(s ) 5(2 26 29 34 39 1)

5-2 5-6 5-10 544+.” (2.4.7)

+__——+——-
6 504 1056 384

3.141 592 653 590 07.

The actual error is 2.8 x 10‘13 while the error from Table 2.4.1 is 2.1 x 10‘”.
Finally we note in Table 2.4.1 that 50 decimals of 1r are guaranteed by taking h = 1/19 and s = 59
in (2.4.6). This requires 80 function evaluations (fe). It can be shown from (2.3.7) that 50 decimals

of 11 are also guaranteed by taking h = 1/29 and s = 25. The computational cost is 55 fe.

Table 2.4.1 Computation of 1r

 

 

h 3 Upper bound for r(h) rife

l 3 7.50-2* 6
1/2 6 993-5 9
1/4 12 986-10 17
1/5 15 2.07-12 22
1/8 25 1.7220 35
1/10 31 675-26 43
1/19 59 265-50 80

 

 

 

 

 

 

I'This means 7.50 x 10'2.

3. APPLICATIONS OF THE EULER-MACLAURIN SUMMATION FORMULA
3.1 QUADRATURE FORMULAS WITH ASYMPTOTIC EXPANSIONS
Sheppard [45] , Becker [6] , and Frame [18] suggested a class of quadrature formulas which are
obtained by taking a weighted average of T(h), T(2h), T(3h), - ' - with weights selected to eliminate one
or more of the terms containing DI, D3, D5, - ' - . This is accomplished as follows.

Let a1, a2, - - - , 03+; be factors of n = (b - a)/h. As before we write 5213-1 = hmDZB_1 and

a‘azlllas+l

A213-1 = bzeEze-l- Dem“ by Azo,-l,2o,-l.---.2a,

(h) the approximation to the integral
-1

b
[(f) =f f(x)dx based on the weighted average of T(alh), T(azh), - - ° , T(asflh) which eliminates the
a

terms A2514, 1' = 1(1)13, but possibly contains some derivative correction terms involving

A2714" ' ' . ’ A27t-l'

Without loss of generality, suppose that

1<al<a2<"'<as+1<n
1<Bl <Bz<"'<Bs (3.1.1)

1<71<72<-~-<7,,

We wish to find constants W1: 1' = l(l)s + 1 such that

3+1 t
a! toaas

+1 _
A231-1,"°,255-I(h) - Z WiT(a,h) + Z ch27i-l
i=1

1=l
(3.1.2)
= 1(f)+ C0A27_1+ "'

where 7 is the smallest positive integer distinct from the Bi and 71" and c], are certain constants. The

constants w, which are found by solving the linear system

14

15

 

 

 

 

 

 

* 1 1 1 T ”w, T "17
2‘3 251 25!
al I a2 . o o as+l wz = O
(3.13)
25 2B 26
Lal S a2 3 . . . as+i WS+I 0
_. L. .1 — _I
determine the asymptotic expansion for the resulting quadrature formula‘
a a 1 8+1
1 3+ =
A261_1’,.'2‘3_1(h) Z W ”am
1 3+1
Z A2 71,42;wa (3.1.4)
=1
- 3+1
= 1(f)+Z A2k-1ZWiai2k
k=v
k¢Bi
k=¥=7j
The principal error term is defined to be
3+1
7127,, Z rig-6,27 (3.1.5)
. i=1

A quadrature formula of degree at least 23+l may be constructed by taking a]. = j and 5k - k in

(3 1.3). In this case, it is easily seen that the following matrix has rank 3+1 and hence a unique

solution exists.

 

— ‘T
l l l ' ' ' 1
1 22 32 . . . (3.1.02
4 4 . . . . '-
1 2 3 (s+1)4 = [12"1’l,,- (3.1.6)
[_1 223 323 (3+1)23

 

The Vandermonde determinant associated with this matrix is nonzero

 

1.3.1110!ch

16
Thus there are countably many quadrature formulas in this class. Becker [6] conjectured but did not
prove this result.

The principal error term is

3+1

Azs-i Z W14” (3.1.7)
i=1

and the asymptotic expansion is

3+1

, ’..., 1 _ _
Al} (h) - Zwioh)
i=1

(3.1.8)
- 3+1
= 1(f) + Z AZk'l Z Wiizs _
k=s i=1

We hasten to point out that this may not be the best quadrature formula. A more efficient
quadrature rule may be obtained by selecting n = (b — a)/h to have as many factors as possible, 1 = 0.1 <
a; < - - - < as” = n. This may be seen by comparing the principal error terms in the 7-point
formulas A}? and .4336, which are h6DS [840 and 3hsD7/2800, respectively.

For reference we list in Tables 3.1.1 to 3.1.4 several formulas constructed by this method. The
entries in the tables may be understood by comparing them with the derivative corrected (DC)

Simpson’s Rule where h = (b - a)/2n and f,- = a + ih.

Agz [16T(h) - T(2h) - Ell/15

115% + 16fl + l4f2 +---+7r2,,1-;’§h21f'(b)- /'(a)1

6
1105 11807 111009

- — + - +
10) 9450 75 600 712 800

 

(3.1.9)

b ..
f f (x)dx + Z hz‘b2,DZS-1(l6 — 4‘)/15.

a 3=3

The principal error is the first term in the asymptotic expansion, 41605 [9450.

l7

. 51 u _m c5 23:. u u... 23 38:8 fission 88:2. 09.

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

mm. T 0mm 3m 81 3: SN 3%: 3; mm? is waneom on
T «m 3% T 25 v 2w N 33 am? a: Manson
W. T ms. «,8 T 0% 3 E 52 em? 4.: 3898302 8
T N: 581 ca _ mwlw 32 em? a...» .maovceaoz
1%. N R1 EN mm ES mm... 22.0; on
a or 2 pm. 3; mm? 583
1m _ ow- v8 _ Mme Em a”? 28m on
_ ON. 3 mm 33 em? 28m
..wlm T 03 WK v? “58% on
T 2 1mm L< “58$
mum T ; .o_lw 3.: m? 2.88 ”cease on
T a .w E: m "a. 258m yaoaaa
.mlw T 2 min 33 NT Saea on
T v m :5 «T 82.5%
W1 . _ 3: 2 1282:: on
_ _ 3: _< 388%;

_m a. a... a my a. Q co swam Ea? 23. 0.52

 

.8358... 22530 U320 do 3:230 36 0:3

 

l8

.:d>aleu.~v:a.5+auc

 

 

 

 

 

 

 

 

 

mflm Sn 2 w: R on 8 a: 3. 3n 8 02 mm a...“ 3 a: 3 3n 3 an 2 Ego: $3 an? ...—$.53on on
.5 v3 «8 ~ «2 v8 _ one «8 a N2 «.8 _ SN .mmlmm :35 ”NM? ...; wanes.
W1 ...: n mg m o: b m8 o So 22 2o 0 a: h 9: m Bub £2 em? #2 3895.32 on
T S e: 2 E 2 SN T bw. $2 9%? ...: 3258.32
We. 2 E 8 _ 2 8 S o... 5 mm 3? mm? 0633 on
d N n _ o ~ m _ $1 3; «MT. 0383
WI 5 3:. 2n 2.. in :~ 1mm 3:. em? 225 on
n 3 NM 2 S h Iwmw SE em? 28m
W 8 as ”2 § m2 3 1mm. «T .59: on
m o w w m m 1mm... c? EBA—1m
1W1 2 8 Z a 2 m“. g a a? 28% £935 on
_ N m m _ Wm 35 m? 28% 3.5355
WIN: 5 z o— N. Mal— }t. m NM< coma—Em Un—
~ N v g ..ml gm «7. .8355
.m. .m _ .m. a 3: _< 32835 on
1w _ w a 3: da. 12885
f .3. ma he a. m a m m c c. we awe"; :33. 23. 252

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

£83. concur—cu 339609 «:8 8.303 ouaflugo N. fin 03$.

 

 

 

 

19

 

 

 

Tamas u Tedd Ea Tan—3: u Tam.
as $2 $3 1 €503 + $33“ a 1 3. SN. 33 m: N? tame? E. a an? ...E 9353. on
5 Q Fee 1 $09.” 133.2 1 08: 1&4 25 35mm: 2:3 mm? .....o 32.53.
8o N: :8 1 €va 1.0.38. m 1 as e: If So 3.23:1 32 an? .5. .3898302 8
94: $8 1 Ama N: 1.13% 1 9a: Taq 8a a}: 32 em? .5. 3,582.52
2Q 5% + $32 1 22 3.3 8m 5}? 2; mm? 2803 on
o: 7% 1.30 1 m: an 9.x}. 3; a”? 583
ai 7.2 1.38 1 e8: Tfiq RN 81m? 3:. ..m? 28m 8
3 $2 138 1 E. Esq aim: Em .M ”< 28m
3Q $2 1 £2 Tao. m8 21m? 3.? 321m 8
22mm: 1 o: qu 3mm .7. .583
9: _ma 1 a. Taq cc :6? E: ...? 258m ...oaea on
a re 1 a Tfiq 8mm 35 n T 2.8% «hogan
m: we 1 o: Tfiq 02 «ET 3% a? Bag on
Q E. 1 a Tm~< 25mm :3 a? Saga
3.3 oEmT 3: ~.... 3282:: on
if N: _m 3:. 2 3282.5
=23.“qu “589:; 5 ESP 3250 8am 35055 305$. 23— 252

 

 

 

 

 

 

.3085me ozoEEAuc. 05 5 «Econ. 3350 was Eobm 22:35:— m._.m 035.

20

Table 3.1.4 [leading Terms in the Asymptotic Expansions“

 

 

 

 

 

 

 

 

 

 

 

Name Rule Abbrev E1 E3 E5 E7 89
T 'dal Al T(h) i- -'-1- 1 '1 l
rapezm 12 720 30 240 1 209 600 47 900 160
. l I -1 1 ‘1 1
DC Trapezo1dal A T ”1) 72—0 30 240 1 209 SW 47 960 136
. 12 l —1 1 -l7
Sunpson A1 SW 136 1121 "12766 W 395
unpson 3 9756 73366 7121176
, 13 1 -1 13 -41
Simpson 3 Second A1 Up!) §6 336 9 233 112
DC 8' ’ Seco 11 A13 U’(h) ‘3 3 ‘13
unpson s n 3 11 20 H 800 W
. 14 1 -l7 13 -4369
S-Pomt A, ‘4‘5‘ T1 0 W 993 7
_ 14 -8 17 -13
DC S-Pomt A3 16 065 80 325 151 470
Boole A13 30*! 975' 966 W9 4 0
DC Boole A35 B (h) 99 225 93 S33
‘ 123 _1- ‘1 7
Waddle A13 WM 840 2 400 6‘73 3 o
123 1 '3 1
DC Weddle A35 V W 137 200 129 360
, 1236 ._3- '5
Newton-Cotes 7-Pt. A13 5 NM} 2 800 3 3%
. 1236 1 ‘3
DC Newton-Cotes 7-Pt. A357 N (h) 154 000
Romberg 9-Pt. A135 W01) 4 725 13 7”
1248 I -512
DC Romberg 9-Pt. A357 W (h) 7 952 175

 

 

 

 

 

 

 

 

 

 

, 2s
‘Ezs-l ” D25-1.

21

USpensky [52] uses the expansion of a function in terms of Bernoulli polynomials (see Krylov
[27]) to derive asymptotic expansions for the trapezoidal, Simpson’s, Simpson’s Second, Boole’s, and
Newton-Cotes’ 7-Point rules. Except for these Newton-Cotes’ rules, Weddle’s rule, the DC trapezoidal
rule, and the DC Simpson’s rule given by Tanimoto [50] , the asymptotic expansions of this section are
believed to be new; the derivation is based on the Euler-Maclaurin Summation formula.

Becker [6] constructed Simpson’s, Boole’s, Newton-Cotes’ 7-Point, Weddle’s, Romberg’s 9-Point,
and several other rules, but did not construct the correSponding derivative corrected rules as we have

done.

3.2 ERROR ESTIMATES

The quadrature error in (3.1.4) may be estimated from the principal error term by

s+l
a ...a
|I(f) — ABI‘H,BSS+1(}1)I 2». (b-a)|bz,IhZ7M2,Z 1111,1637. (3.2.1)
i=1
Here we have replaced 027.1 in (3.1.5) with (b - a)M27. Also, as before,
7:min[1+-{311...5B51713H'97t}]' (3.2.2)
For example, let us compare Simpson’s formula S(h) with the DC Simpson’s formula S '(h).
The error in Simpson’s formula is less than (b - a)h4M4/180, whereas, the error in the DC

Simpson’s formula is less than (b - a)h6M6/9450. Now if

21121116
< 1
105 M4

 

(3.2.3)

then this additional accuracy of fifth-order vs. third-order is gained at the one-time minimal cost

of computing the correction term -hz[ f ’(b) - f '(a)] [15.

22

3.3 NEWTON-COTES FORMULAS

Examination of Tables 3.1.] to 3.1.4 reveals that we have obtained asymptotic expansions for
five of the Newton-Cotes’ quadrature formulas: the trapezoidal rule, T(h), Simpson’s rule, S(h), Simp-
son’s Second rule, U(h), Boole’s rule, B(h), and the Newton-Cotes’ 7-Point rule, N(Iz).

Examination of

A1501) [26 T(h) - T(Sh)]/25

h
5.69% + 52f1 + 52f2 + 52f3 + 52f4 + 42f5 +---+ 21f,,] (3.3.1)

10”) + 59911403/18 000 - - --
reveals that it is impossible using the technique described in Section 3.1 to obtain the 6-point Newton-
Cotes’ rule. Therefore, this class of quadrature formulas is not merely a restatement of the Newton-
Cotes’ rules.

As previously noted, Uspensky [52] shows how another technique may be used to obtain an

asymptotic expansion for any Newton-Cotes’ quadrature formula.

3.4 THE MIDPOINT RULE AND SOME OPEN FORMULAS
An asymptotic expansion for the midpoint or centroid formula, C(h), may be derived from the

Euler-Maclaurin Summation formula by writing

C(h) 2701/2) - T(h)

11 Z f(a +171/2)
i=0

b (3.4.1)
= f f(x)dx - 191/24 + 7E3/5760

-3155/967 680 +---+ (2145- 05234112,

+ R(h,2s)

23

where h = (b - a)/n and the absolute error lR(h, 23)l = lR(h, 23+ 1) l is estimated by

1
IR(h.2s+1)|= I] IZFlié’Wr) - FAM’lemuwr
O

 

(3.4.2)
4
- 7r . . h 2 2
Of course, (2.3.5) and (2.3.6) may also be employed to obtain additional estimates for R(h, 2.5).
As before,
M = max I (2" x I 3.4.3
25 a<x<b f ( ) ( )
and
DZS-l = f(23-l)(b) _ f(23-l)(a)
521.1 = hZSDZH (3.4.4)
AZs-l = b2sE2s-l-
From (3.4.1) we immediately obtain the derivative corrected (DC) midpoint formula
11
C'(h) = h Z f(a+ih/2) + 151/24
i= 0
b (3 4 5)
=f f(x)dx + 753/5760 — 31E5/967 680 ' '
a
+ + (2143— 1)h25b25.023_1 + 1201,21).
The error R(h, 2s) may be estimated by
IR(h,2s)l < 3(b-a)M4(h/21r)4. (3.4.6)

Here 3 is an approximation for (3114 [40% ),

The principal error term of the DC midpoint quadrature formula is 7h4D3/5760. This compares
favorably with the Simpson’s principal error term, h4D3 / 180.

Squire [47] observed that Simpson’s rule is “probably the most widely used integration formula.”

Thus it is of interest to compare the third-order DC midpoint and Simpson’s quadrature formulas.

24

Clearly the DC midpoint rule is easier to apply than Simpson’s rule since the weights in the
former are unity whereas there are three different weights for (composite) Simpson’s formula.

If n represents the number of times the basic or holistic quadrature rule is applied, then tra-
ditionally, for Simpson’s rule, it = (b - a)/2n while for the midpoint rule, h = (b - a)/n. Thus one is
tempted to compare the midpoint and Simpson results for the same step size h. Rather one should
compare results when the holistic rule is applied an equivalent number of times, 11. (Note that this
is done in Chapters 6 and 7.) Milne [36] also discusses this.

When this method of comparison is used, then one sees that for 11 applications, the DC midpoint
principal error is exactly 3.5 times the Simpson principal error. However, the DC midpoint rule uses
only n+2 function evaluations (counting a derivative evaluation as one function evaluation), while
Simpson’s rule requires 2n + 1 function evaluations.

Therefore, for the same order of magnitude error, the Derivative Corrected midpoint formula
(3.4.4) requires approximately half the number of function evaluations as Simpson’s mile. This

represents a considerable savings in computer time.

Table 3.4.1 The Midpoint, DC Midpoint, and Simpson’s Rules
6 (1):
A 1' dt — = ln 2
pp 1e 0 f; x

 

Midpoint DC Midpoint Simpson

 

Error nfe Error nfe Error nfe

 

3 3.39-3* 3 -7.97-5 5 -2.26-5 7
6 8.63-4 6 —5.20-6 8 -1.48-6 13
12 2.17-4 12 —3.28-7 l4 -9.4l-8 25
24 5.42-5 24 -2.08-8 26 -6.20-9 49

 

 

 

 

 

 

 

 

 

‘This means 3.39 x 10‘3.

25

As can be seen in the Table 3.4.1, for n = 6 the DC midpoint rule produces an error of 5.20 X 10'6
using only 8 function evaluations while Simpson’s rule requires 13 function evaluations to produce an
error of 1.48 X 10’6.

Thus, in situations where Simpson’s rule is considered adequate for the approximation of a
definite integral, the DC midpoint rule should be considered as a viable alternative to Simpson’s rule.

Next we use the midpoint rule (3.4.1) to derive several Open formulas based on the methods of
Section 3.2. A numerical integration formula is defined to be Open if all of the nodes are interior to
the domain of integration.

Let a1, ° ' ' , 5+1 be factors of n = (b - a)/h. Denote by QZaalltiaffl

2‘3 _1 (h) the approximation
S

b
to [(f) = f f (x)dx based on the weighted average C (alh), ' ° ' , C(as+1h) which eliminates the
0

terms containing AZfix-l’ . . . , A2534, but possibly contains 1 derivative correction terms
involving A271_1, ' ' ' , A271-1.

Assuming that (3.1.1) holds, we seek constants W1, 1' = 1(l)s+ 1, such that

3+1 1

“1"'as+1 _
Q251_1,...,233-1(h) — Z w.C(a,-h) + 962,].-.
' l

"1 1': (3.4.7)

= [(f) + c0(21'27-1)A27_1 + .-

Here 7 is the smallest positive integer distinct from the 31' and 7i and c}- are certain constants. Now
using (3.1.3) to find the numbers wi, we obtain an asymptotic expansion for the resulting quadrature

formula

3+1 1 5+1
at"'as+1 __ 1-27. 27‘
Q21.-1.---,21a,-1(") .- Z w,C(a,-h) - Z (2 1-1)A271_IZwiai 1
i=1 i=1 i=1
(3.4.8)
b 6- 3+1
=f f(x)dx + Z (21'2k-1)A2k_12w,.a}"
a k=7 (=1
119213,.

k¢7i

26
having the principal error
3+1
(2“27 - 1)Az,_1 Z 111,637 (3.4.9)
1': 1
Several open and derivative corrected open quadrature rules obtained by this technique are
presented in Tables 3.4.2 to 3.4.5.
For reference we state the open quadrature formula ng which is the analog of the derivative
corrected Simpson’s formula.
Q? = 5,502)

[16C(h) - C(2h) - 51/21/15

$312111 +1112) - fax +1) + 210 + 311/2) + ~-

2
+ 210 -3172) - 101-11> + 2101-1112)] £5110)- 1(a)] (3.4.10)

3111605 127 11807

= + .-
IU) 302400 9676800+

 

b
= f f(x)cbc + 11211222,-101-21-11(16-411/15.
0 s=3

Ico’as+l

The open quadrature formula Q:;:;-1 213 _1(h) based on the asymptotic expansion of the
’ , s

midpoint rule is the analog of the closed formula A; 6114a“ ’2 35401) which is based on the Euler.
Maclaurin Summation formula.

In particular, the open formulas S0, U0. Bo, V0, N0, and W0 are the analogs of Simpson’s,
Simpson’s Second, Boole’s, Weddle’s, Newton-Cotes’ 7-point, and Romberg’s 9-point formulas,
respectively.

Finally‘we note the following relationships:

30111) = 14c<10~ «2101/3

3001) [1650(h) - 50(211)]/15 (3.4.11)

111001) [648001) - 30011)] /63.

27

 

 

 

 

 

 

 

..a «an _m Es :36 new.
as _ 1 1 2.... oz 6 an
lowl _ 6mm 4% _N 3; SN _ E .3 SN .0

_1 S can _1 68 e Mmmm $63 an “o
9% _1 m3. 2: m1 68 3. bolmb 3? ammo
1 1 E 6 n2
_ N: Sm £33 _ E 2 £20
mb 1 .mlvm o 2
m N 2 EN . _ As .> 26
1 .E. o 2
_ o 2 _ E > n20
1mm E o 2“
N 2 ow. 2.8 _ _ E .m z _o
1 o m—
_ 8 3 mm 2: m 4:0
w 1 .8 6 ...
pm _ 5 .._.1 3 .2 m _o
— 0 1M1 5: 0: M —0
pm _1 2 Mm Eon mo
1 1m. 6 _
~ v _ it m N _0
Hum _ _ xtfb ~O
_ _ \50 _O
2 a e v n N _ ab
m u u o u o o .8 .22.... 55.2 63.

 

 

 

 

 

 

 

 

 

 

5836.55 2335.30 :30 15 20:92.28 N66 632.

 

28

q

 

 

 

 

 

 

 

 

2.3:. 1 a 1 N 223 NE 1. N 1...
NN: NNN NNN NNN N2 NNN NNN NNN _1 NNN NNN NNN N2 NNN NNN NNN NNNSN E63 NNNo
NN NN N1 NN NN N1 NN NN N1 NN NN N1 NN NN . NNN.
NNN o N:
NE NNN1 N3 N4 N2 NNN1 N: T N: NNN1 NNN N4 N2 NNN1 N: 1m“ E 3 N620
NW1 NNN N 8N T 8o N NNN N NNN T NNN N 8° N SN _1 NNN N BNE sz ommmo
2N NNT NNN SN 2:1 SN NNN NNT 2N old E62 N20
N N:
.NIN. 1 1 1 MR 6 NN
N NNN N_ a N: N_ N: a N. N: N E.> N20
.ol_ 6 N.
N T N N T N N T N N E > NN_o
B m... o a
N NNN 3.1 NNN N NNN 91 NNN 5. EN NNNO
1 1 1va o N.
N: S 2 _ S 2 N. s. E N 320
B 9.1 1 N
N NN NN NN N E. a 20
N o _
N N N N E a 20
.8 1 m. o N
_ N N N N E.N So
1 11N1 o N
N _ N N E N 20
R _ .N Eb o
_ N
_ N Eu .0
. _
Nm 2.“ z. NNN N: 2. SN N. N. N. NN m. a N. NN : .6 8.83 .653. 23.

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

.3th 5:60:00 o>zm>ton Ea 3:363 232330 min 635.

 

29

. ~1m~m £0 .11 Tmmd b5 _ImND an: n TmNWe

 

 

 

 

 

 

NNN SENS 1 AmoceNN + NNNVSN _N 1 a: NNN. :1NN1_NV _1NN< N: 3318: EN; Nmmo
NNN N53 1 NNNQN 1. N33: 1 N8 3 21.1.33 TNNq 8N NN\NmNN_1 E63 Nwmo
o8 NENNN 1 AmeNNv 1. ANNE: N 1 NNN N: :1NN1_NV7NN< SNNNN NNENNN. E...z ammo
SNENN 1 NNeN: + NNENNN 1 NNN : 2121.8 TNN< SNNNQNENT Eoz ammo
NNNENeN + ANENN 1 NNN. :1NN12NV _1NN< 8... EN ENENN E...> . Nmmo

22% 1. ANN: 1 N: 21.1.39 _1NN< NNN NN\NN.:N1 E6> Nmo

$352 + NNNEN 1 No : :1NN1_NV_1NN< 8N N: N\NmNN_ Em. ammo

NENE .1 N38 1 S. :1NN1_NV_1NN< ON. ENE? Eom cmo

82% 1 _N_ :1NN1_NV_1NN< 8N Emma E...: Nmo

N236 1 N. :1NN1_NV _1NN< o§NmN1 E6: Lo

22.1... 1 N_:_1NN1_NV_1NN< 84 NNN}: ENN Nmo

szv 1 a :1NNLNV _1NN< o3 :NmN1 EON Nflo

:1NN12NV TNN< 8N N\NmN Eb _o

:1NN1NNV £3 E _m1 Eu Ho

coufimxm 633.333 5 ENE _Socoo Norm 390:2 >2nn< 23—

 

 

 

 

..NcomNcauxm ozoafixfix 2: 5 NENNF 38:00 23 team .2635 vim 03¢.

 

30

Table 3.4.5 Leading Terms in the Asymptotic Expansions“

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Rule Abbrev E 1 E3 E5 E7 E9
Q1 C111) 1 7 —31 127 -511
24 "5 76—0 967 680 "154—828 800 24—524 881 900
5—760 ‘6—9 7 68—0 "154 82 8800" '2'4'_5'2_4 8'8‘1‘9'00'“

12 -7 31 —127 1241

01 56”” ”—“1 440 4—8 38'4“ 1 843 200 —175"1'7'—7 7'28"
12 1 31 -127 511

03 sow 302 400" '9 6‘76 800' 364 953 600'
13 ;7_ 31 -1651 2993

Q1 ”6”" 640 10 75 2’ 2 45 7' 600 ’19 464192

013 U, (11) 93 -381 6643
3 o —358 400 W34 4 4—32 5_37 600
124 —31 127 —8687

Q13 80”” 15120 115 200" ‘18' 24—7 68—0“
124 1 127 -73

Q35 311/") 3 175 200 _3 421 440—
123 -31 127 -3577

013 V6”) 26 880 307 200‘ '32 44""‘0‘3‘2‘0
123 . 381 -73

Q35 Vo/h/ —17"5_6"1'60_0 *9 461 760
1236 —381 365

Q135 New 358 400 270 336
1236 1 219

0357 NJ") '1'1"2"64‘0‘00
1248 -127 1241

Q135 wow 37 800 _171 072'
ms 1 _L_

0357 wow 581644800

.. ZS
'Ezs-i ”h D23-1.

 

31

For the derivative corrected formulas we have
56 (11) = [16C'(h) - c'(211).]/15
8601) = [64S6(h) - Sd(2h)]/63 (3.4.12)
196(k) = [256Bd(h) - Bd(2h)]/255.

This leads to the observation that the derivative corrected Romberg quadrature to be defined
in Section 3.6 may be based on the asymptotic expansion for the midpoint rule, (3.4.1), in place of
the Euler-Maclaurin Summation formula. In this case the first 4 columns of the “open Romberg”

table are given by the open quadrature rules C (h),SO(h),BO(h), and W001), respectively.

3.5 ROMBERG QUADRATURE

The Euler-Maclaurin Summation formula is a tool of strategic theoretical importance. Indeed,
Romberg [44] proposed a new class of quadrature formulas for a finite closed interval [a,b] based on
Richardson’s extrapolation technique [24, 43] applied to the Euler-Maclaurin formula. It involves re-
peated halvings of the integration interval and successive elimination of higher order terms in the Euler-
Maclaurin expansion.

An extensive discussion of the theory is given by Bauer, Rutishauser, and Stiefel [5].

Romberg concluded that Richardson’s deferred approach to the limit would improve the accuracy

of the trapezoidal rule.

T(h) = h[l—f0 + f, + + ifn] = h Z"f(a+ih). (3.5.1)
2 2 i=0

For n = (b - a)/h even, he obtained Simpson’s formula by writing

S(h) = [4 T(h) - T(2h)]/3. (3.5.2)
Next for n a multiple of 4 he obtained the Newton-Cotes’ formula of order 6, Boole’s Rule:

301) = [16S(h) — 3(211)]/15. (3.5.3)

In the next step Romberg obtained the formula

“’01) = [643(k) - 8(2h)]/63
_ 4h
- 33;}? [217f0 +1024f1+ 352f2 +1024f3 (3.5.4)

+ 436f4+1024f5 +352f6+ 1024f7+434f8+ +217fn]

32
which is not a Newton~Cotes’ quadrature rule. Thus Romberg’s method is not a reformulation of the
NewtomCotes’ formulas. Apparently, W(h) was first derived by Sheppard [4S] and later rediscovered
by Becker [6].
Now let
kl!
To. = hZ f(a+ih) (3.5.5)
1' =0
be trapezoidal sums where h = (b - a)/2k. Recalling the Euler-Maclaurin formula
8
T(h) = 1(f) + Z cahZ“ — R(h,2s) (3.5.6)
a = l

where ca = b2aD2a-1’ it is easy to understand the definition

ka = [4me-l,k+1_ m-1,k]/[4m -1]. (3.5.7)
In fact
_ 20’21’,,,,2m b-d
ka - Al,2,-~,m (T), m > 0. (3..58)

From this the Romberg T-table is constructed:

(3 .5 .9)

We have already seen that the first three columns are the trapezoidal, Simpson, and Boole’s values,
respectively.

Bauer, Rutishauser, and Stiefel [5] show for any f e C2m+2[a,b] there is a f 6 [0,1] such that

IT [W < (b -a)82m.2 If‘2m*2’(r)l (3 510)
"'k (2m+2)! 2<2k-m><m+1> ' ' '

 

Combining this result with lemma 2.1.2 we obtain a new and much more convenient error estimate for

any entry in the Romberg table.

33
Theorem 3.5.1

If f e C2"'*2[a, b] then the error for any entry in the Romberg T-table may be estimated by

(b -U)M2m+2

T - I <
I mk UN 45 "2(m-1)2(m+1)(2+2k-m) (3'5'11)

 

where

M2,“; = aff§b|f(2’"+2)(x)|- (3.5.12)

3.6 DERIVATIVE CORRECTED ROMBERG QUADRATURE
Lanczos [28] shows how the addition of only one derivative correction term can significantly
improve the accuracy of a quadrature formula at the minimal expense of a small increase in computa-
tional effort. He also states the derivative corrected (DC) trapezoidal and Simpson’s Rules.
We note that the DC trapezoidal rule generates the DC Simpson’s rule.
S'(h) = [16T'(h) - T'(2h)]/15 (3.6.1)
which in turn generates the DC Boole’s rule
B'(h) = [64 S’(h) - S'(2h)]/63. (3.6.2)
Next, the DC Boole’s rule generates the DC 9-point Romberg formula
W'(h) = [256B'(h) - B’(2h)]/255 (3.6.3)
which is distinct from any DC Newton-Cotes’ rule.
This suggests a new quadrature scheme, the derivative corrected Romberg quadrature. Let
h = (b -a)/2" (3.6.4)

and

2" 3
C3,, — h Z: f(a+ih) — Z AZH
”0 “=1 (3.6.5)

TOk - (1’01”?)

be DC trapezoidal sums. Define

34

Crsnk = “much-mu ‘ qr-1,k)/(4m H '1) (3'6-6)

and construct the DC Romberg C3-table

3
C00
5'
C01 Cil
S 3
C32 C12 C22 (3.6.7)

5' s s 3
C03 C13 C23 C33

The first column is the Euler-Maclaurin formula with h = (b - a)/2" and the m-th column, m > 0, is

given by the quadrature formula

0 I .. _
cf", = .42 '2 " '2'" (L1). (3.6.8)

3+1,,+2,...,s+m 2):
For the case s '= 1, the first three columns of the Citable are the DC trapezoidal, DC Simpson’s,
and DC Boole’s rules, respectively.
In general, comparing (3.5.8) with (3.6.8), we see that the m-th column of the C‘-table is the

derivative corrected quadrature rule corresponding to the m-th column of the Romberg T-table.

3.7 A NUMERICAL EXAMPLE
To illustrate, we employ the Romberg and derivative corrected Romberg quadrature formulas to

estimate the integral

1
f l—sin(1rx)1rdx = 1. (3.7.1)
0 2
We note that
h = 2"‘

DZa-l = (_1)a1r2a

2k-1
C3,, = «21* Z sin(jn2"‘) (3”)
Fl

8

<I>(h,s) = Z(-1)an2“2'2kab20.

a=l

35

Table 3.7.1 Romberg Quadrature T-table for L: V2 sin(1rx)1rdx = l

 

(Newton-Cotes

k 0 (Trapezoidal) 1 (Simpson) 2 (Boole) 7-Point)

 

0 0.900 000 000 000
1 0.285 398 163 397 1.047 197 551 20

2 0.948 059 448 969 1.002 279 877 49 0.999 885 365 912
3 0.987 115 800 973 1.000 134 584 97 0.999 991 565 473 1.000 00; 774 99

 

 

 

 

 

 

 

Table 3.7.2 Derivative Corrected Romberg Quadrature (s = 1) Cl-table

 

 

m (Corrected (Corrected (Corrected (Corrected ,
k Trapezoidal) 1 Simpson) Boole) 3 Newton-Cotes
7-Point)
0 0.822 467 033 424
1 0991 014 921 753 1.00; 251 447 64
2 0.999 463 638 558 1.000 026 886 34 0.999 99_1_ 575 848
3 0.999 9_6_6 848 370 1.000 000 895 69 0.999 999 9:15 204 1.000 000 008 14

 

 

 

 

 

 

 

Table 3.7.3 Derivative Corrected Romberg Quadrature (s = 2) CZ-tablc

 

 

m (Corrected 1 (Corrected (Corrected (Corrected ,
k Trapezoidal) Simpson) 2 Boole) 3 Newton-Cotes
7-Pomt)
0 0.987 757 437 638
1 0.999 470 572 017 1.000132 685 26
2 0.999 99; 116 699 1.000 000 895 19 0.999 999 876 401
3 0.999 999 878 254 1.000 000 001 45 0.999 999 999 999 1.000 000 000 08

 

 

 

 

 

 

 

36
Table 3.7.4 Romberg Error for fol 16 sin (nx)1rdx = 1

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

m 3
k 0 (Trapezoidal) 1 (Simpson) 2 (Boole) (Newton-Cotes
7 Pornt)
0 1.00+0
1 2.15-1 -4.72-2
2 S.19-2 -2.28-3 7.15-4
3 1.29-2 -1.35-4 8.44-6 -2.78-6
Table 3.7.5 Derivative Corrected Romberg Error (3 = l)
m (Corrected (Corrected (Corrected (Corrected ,
k Trapezoidal) 1 Simpson) 2 Boole) 3 Newton-Cotes
7-Point)
0 1.78-1
1 899-3 -2.25-3
2 536-4 -2.69—5 8.42-6
3 332-5 -3.96-7 2.48-8 -8.14-9
Table 3.7.6 Derivative Corrected Romberg Error (s = 2)
m (Corrected 1 (Corrected (Corrected (Corrected .
1‘ Trapezoidal) Simpson) Boole) 3 Newton-Cotes
7-P01nt)
0 422-2
1 5.29—4 —1.33-4
2 788-6 -3 .95-7 1.24-7
3 122-7 -1 .45-9 909-” -2.98-11

 

 

 

 

 

 

‘This means 4.22 x 10‘2.

 

 

 

37

The values in the C l-'I'ab1e 3.7.2 show a marked improvement over those in the Romberg T—Table

3.7.1. For the calculation of the Cl-table, it should be emphasized that

D.” = f’(l) - f'(0) (3.7.3)
is computed only once, in fact before the calculation of the‘Cl-table commences. Thus the first deriva-
tive of the integrand is evaluated only at the two end points of the interval of integration.

The advantage of the derivative corrected Romberg quadrature over the classical Romberg
quadrature is its increased accuracy and efficiency. Indeed, it can be shown that the m-th column of the
CS-table is given by a quadrature formula of order h2lm +1“ 2‘ as compared with an h2(m +1) order for
the m-th column of the Romberg T-table (m=0,l ,2,‘ ' -).

This improved accuracy is gained by the minimal cost of evaluating the derivative correction terms,
Dza-li once. Moreover, the derivative corrected Romberg quadrature is appealing because the trape-
zoidal sums are closely related to Riemann sums, and the weights are easy to program. As noted on
page 31, an extrapolation procedure may be applied to an asymptotic expansion for the midpoint

rule in place of the trapezoidal rule.

4. FUNCTIONS OF TWO VARIABLES
4.1 THE EULER-MACLAURIN SUMMATION FORMULA
Let C” [R] denote the set of functions f (x, y) of two real variables where the partial derivatives
(see 4.1.5) f “5(x, y), 0 < a, 5 < 23, exist and are continuous on the rectangle R = [c, b] X [c,d] and

let f (x, y) e C23 [R]. We wish to estimate the integral
(1 b
Itf) = f f f (x.y)dxdy. (4.1.1)
c a

Divide [a, b] into n equal parts each of width in = (b - a)/n by points x,- = a + ih setting x0 = a and
x" = b. Similarly divide [c,d] into m equal parts each of width k = (d - c)/m by points yl- = c + jk
where yo = c andym = d.

Now for 0 < t, u < 1, average the values of (b - a)(d - c) f in each subrectangle
Rf], = [x,,xi+ 1 ] X D’i’YJ’HI at points at a distance of t from the left side and a distance of u from

the bottom by writing

m-l n -1

Fhk(t,u) = th Z f(x,- + my]. + uk). (4.1.2)
[=0 i=0

Then the double integral of this moving average over the unit square is exactly the integral I (f )

over the rectangle R:

1 1
J J Fhk(t,u)dtdu
O 0

m-l n-l

l l
= Z Z] J f(x,- + my]. + uk)hkdtdu
0 0

M i=0 (4.1.3)

""1 "'1J~H(i+l)kJ~a+(i+l)h

a+ih

f (x .y)dxdy

38

39

The average T(h,k) is the (composite) trapezoidal rule:

T(h,k) = %[F(0,0) + F(1,0) + F(O,l) + F(1,1)]

(4.1 .4)

m n n H
= ME 2 f(x,-.m-
i=0 i=0

The double primes on the summation signs indicate weights are to be assigned as indicated in

 

Figure 4.1.1.
14 V2 ‘26 14
V; 1 1 )‘z
16 1 1 95
$4 ‘26 V2 %

 

 

 

Figure 4.1.1 Trapezoidal Weights

The mid-value FhkOé, V2) is the composite centroid or midpoint rule and has been investigated
by Good and Gaskins [20].

Before proceeding, we define some notation which will simplify the writing of the Euler-
Maclaurin Summation formula. Let ¢a(t) and ih,,(u) represent the a-th Bernoulli polynomials in the
variables t and 11 respectively, a > 1.

For 1 <a,[3<2s denote

f“"(x.y) = 9f®ff(x,y) W5)
and
m
dao = Z "[f“0(b,yi) -fao(a,y,.)] (41 6)
i=0 . .

Z" [10504.40 - f°“(xi.c)l
i= 0

don

d... = f“"(b.d) - Mme) - f“"(a.d) + f“"(a.c).

40

The double primes signify trapezoidal weights, that is, the first and last terms in each sum are to be

assigned weights V: and the remaining terms are assigned weights 1. The weight assignment is illus-

trated in Figure 4.1.2.

-V2

-1

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

55 l 1 $6 -1
)6
l
1
)6
415 -1 -l 375 1
duo “'06

Figure 4.1.2 Partial Derivative Weight Assignments

and

Next for a,[3 = 1(2)Zs - l we define

[F“0(1,l) + F“0(1,0) — F°0(O,1) - F“0(0,0)]

tolu—

DaO =
1
Doc = 31" 0"(1,1) - F0305) + Foam) - F°B(o.0)1

0,, = F“"(l,l) - F“B(l,0) - F““(O,l) + F“"(0,0)

1
Ia,0 = if [F“°(t. 1) + F“°(t,0)]¢a(t)dt
0
= [F°B(1,u) + F03(o,u)] wp(u)du
1
[2:43 =f [F2515(t,1) - F2"”(r,0)] ¢23(t)dt
O
1
1,1,2, =f [Fa’23(l,u) - Fa’2‘(0,u)] (112$(u)du
0

l 1
la,“ =J’ J‘ F“a(t,u)¢a(t)¢a(u)dtdu.
0 0

-1

dad

(4.1.7)

(4.1.8)

41

Now since
m—I n-l
Fa5(t,u) = 11““chl Z Z f“‘3(x,-+th.y,-+uk)
i=0 i=0
we see that
.. 1
Duo — ha+ kdaO

- +1
- +1 +1
I)“, - h“ k5 dag.

Finally, for unequal positive integers a, B we define

an = an = DaO + DOa
Eafi = EBa = Dali + DBa
Eaa = Duo.

and
LaO = L0 = ’0 +100.
Lafi = Lfla = 121-,3 + [a,23

L=l

aa aa'

Theorem 4. I. 1 (EULER-MA CLA URIN S UMMA TION FORMULA )

Letfe C2‘[R]. Then
’0) = T(h.k) + <I>(h.k;2s,ZS) + R(h,k;25.25)

where

S a.
¢(h.k;28s281 = Zb2a|:-E2a-1,O * ZbZBEZa-l.28-l

a-1

6—1

and

S
R(hkfltzs) = L2s,o ”Z bzaLzs,2a-1 + L2$.25'

a-l

(4.1.9)

(4.1.10)

(4.1.11)

(4.1.12)

(4.1.13)

(4.1.14)

(4.1.15)

Lyness and MCHUgh [32] give an n-dimensional formulation of the Euler-Maclaurin Summation

formula.

It may be the case for some functions that over the rectangle R, higher partial derivatives exist

42

with respect to the second variable than exist with respect to the first. In this connection we observe

that the Euler-Maclaurin Summation formula may be generalized as follows.
Let C 28’2’ [R] denote the set of all functions f (x, y) where the partial derivatives f “5,

0 < a ‘\< s, 0 < B < r exist and are continuous on the rectangle R. Then

Theorem 4.1.2
Iffe C231” [R] and 0 < s < r, then
[(f) = T(h,k) + d>(h,k;2s.2r) +R(h,k;23,2r)

where

S a .|
¢(h.k;23.2r) = Z b2a[-E2a-l,0 + Z bZBEZa—l,2£3-ll

a=1 6:1 .5

r r s
+ Z bza'L‘Do,2a-1 + Zb26026-1,2a-1]
a=s+1 (i=1

and

S
R(h.k;28.2r) = 123.0 + [0,2r " Z b2a(12s,2a-1 + IZa-l.2r) + 125,25"

a=1

4.2 ERROR ESTIMATES

Let

Mae = max lf“5(x.y)l
' (XJHR

hzskaMzs a + haksta 2 a < 2s

S

N23, 0. =
10.74(hk/21r)2‘M23,23 a = 23

0.74 z n2/6\/§
1.48 z n2/3x/5—
2.95 as 2112/37?

7 = (b-aXd-c).

(4.1.16)

(4.1.17)

(4.1.18)

(4.2.1)

We now give an estimate for the remainder term in the Euler-Maclaurin Summation formula.

43
Theorem 4. 2. I

Letfe C25 [R] and

8
P2: = lbi'st,0 + Z Ib20|N23,2a-1 + N2s,2s- (4-2-2)
a=1
Then
|R(h,k;2s,2s)| g 2.95(b -a)(d —c)P2,(2n)-2S. (423)
Proof:

For0<a¢6<23wehave
m-l n-l
IF“B(t,u)I g ha+1k5+1 Z Z lf“5(xi+th,yj+uk)l g WWW”. (4.2.4)
i=0 i=0
Take absolute values in (4.1.8) and apply (2.3.1) to obtain the inequality
11,61 < 2.95 7Mak‘3MaB(21r)'2‘, a = 28 or B = 2s (4.25)

where A = 1/2 if a. or B = 0 and A = 1 otherwise. We estimate [233$ as follows.

112,3,1 < 1.4827(hk)25M28’23(21r)'43- (4.2.6)

Finally from (4.1.12) we obtain the result
3
IR(h,k,2s,2s)l g les,0' + Z lb2aL2s.2a-ll + les,2sl

a=1

S
< 2.957(21r)'23[|b1|N23’0 + Z IbzalN2,,za-1 + st,2sl (4.2.7)

a=1

= 2957134211)“ 0

Theorem 4.2.2
If there exists a positive constant C such that for 0 < a Q 23, 0 < k < h < 7r/2, M2”, < C, and
Ma,28 < C, then

lR(h, k;2s, 2s)1 < 5.93(b — a)(d - c)C(h/2rr)2‘. (4.2.8)

44
The obvious result holds if 0 < h < k < 1r/2. Here 5.93 is a bound for the expression

ffiL+ «2 + h/12 4
3x5 2 IN? 1-(11/217)2 (.2-9)

 

The proof of this and several succeeding theorems are similar to the proof of Theorem 4.2.3 and

therefore are omitted.

Theorem 4.2.3
If there exists a positive constant C such that for 0 < a < 28, 0 < k < h < 1r/c,

M < C25”, and Ma 2s < C23“, then

25',a
hC ZS
lR(h.k;23,28)| < 7.17 (0 -a)(d - ”(317) (4.2.10)
Proof:
Apply the bounds on the partial derivatives and the well-known result
2(21r)'2“
b < ——
| 2a| 1_ 21-2“ (4.2.11)

to (4.2.7) to obtain

4 2 hC 23 3 2 2a-1 2 25
lR(h.k;2s,2S)| < -—"—7(—) 14 Z 1&92— + .1— £9
3V5 _ 3(21r) a 12¢; Zn

2 23F ]_hC2 25' 2 hC 23
€51.7(E> 1—+<E>———( M) + n <—) (4.2.12)

 

 

 

 

3J5 2w 2 12 l-(hC/21r)2 1% 211
4"2 (hCY‘ 1 n2 710/12
g — 7 —- — + + _—
3\/§ 2" 2 12x/5 1-(hC/2rr)2
Finally since hC < 1r,
fifiF... "2 + hC/12 J
3x/5 2 12\/5_ 1-(hC/21r)2
(4.2.13)
47r2[1 112 1r]
<——-—+ +——<7.17.6
3x/3' 2 12¢? 9

45
4.3 ADDITIONAL ERROR ESTIMATES

For 0 < a, (i < 23 we note that
12w = ‘123+1.13’ [a,23 = _la,23+1’ [23,25 = I2s+1,2s+1-

Hence we have the inequalities:

”2.3-+1431 g 2.95 7KhZS+lkfiM2s+l'fi(2n)-(23+1)

Haas-+11 < 2.95 7Ahak23+lMaJS+1(2n)-(Zs+l)

”23.231 < 1-432 ”r(hk)2’+1M2s+1.2s+1(21r)'(4”2)

where A = 1/2 if a or B is zero and A = 1 otherwise. Then the analog of Theorem 4.2.1 is

Theorem 4.3.1
lR(h,k;25.23)| = lR(h.k;25+1.23+1)|

< 2.95(b — a)(d - 29102,, 1 (210423+ 1)

where
S
P2s+l = |b1|N28+1,0 + Z |b2a|N25+l,2a-l + N2s+1.2s+1-
a=1
_ n2 hk 2"
N28+l,2s+l _EV‘; E; M28+l,25+1'
Theorem4.3.2

Theorem 4.2.2 is true if 25 is replaced by 7.5 + 1.

Theorem 4.3.3

Theorem 4.2.3 holds if 23 is replaced by 28 + l.

4.4 SHARPER ERROR ESTIMATES

(4.3.1)

(4.3.2)

(4.3.3)

(4.3.4)

Due to the asymptotic nature of the Euler-Maclaurin series, the following error estimates will

often provide closer error estimates than those given in the two previous sections.

46
Theorem 4.4.1

lR(h,k;2s+l, 23+1)| < 5.31 (b — a)(d _ c)p23+2(2n)-(2s+2)

where
S
P2s+2 = 'bIIN2s+2,o + Z lb2a|N23+2,2a-l + N2s+2.2s+2-
a=1
_ hk 2‘ M2s+1,2s+1
N2s+2.2s+2 ‘ '2; ——\/'6— .
and
5.31 z 114\/6_/45.
Proof:

Recall from (2.1.2), (2.1.7), (2.1.8), and (2.1.10) the following properties:

1
J; ¢§(t)dt = (-1)"“b23 < 2.17(2n)'2‘3

1 2

2.17 z 114/45.

(4.4.1)

(4.4.2)

(4.4.3)

For convenience we let i signify + if B = 0 and - if B > 0. Also, let 1' signify - if {3 = 0 and +

if B > 0. Then integrating by parts and applying (4.4.3) we find

1
425,3 = -fO[F2"B(t,1) i F25’5(1,0)]¢és+1(t)dr

' szs+102sp + 128+l,{3

123+] .13

‘ iI’2s+21)2s+1,11 - [28+2,B

1
if [F2’+2"3(t.1) i F23+2’3(110)11ib2s+2 ‘ ¢2s+2(’)] ‘1’-
0

(4.4.4)

47

Apply the Cauchy-Schwarz inequality and use (4.4.3) to obtain

1123-0-14” < 27h28+2kBM25+2,fi(b225-+2 - b4s+4)%

(4.4.5)
< 5.317h23+2k3M25+2,B(21r)'(25+ 2).
Similarly
”111,2: = Ia.2s+1
(4.4.6)
11,135.11 < 5.317hak23+2Ma,25+2(21r)‘(2‘+2).
Moreover
[23,23 = [23+1,25+l ‘ (4.4.7)
Applying these results to (4.1.15) we have
R(h,k;28,23) = —R(h,k;28+1,2s+1)
8 (4.4.8)
= “L2s+1.0 + ZbZaL23+l,Za-l ‘ L25+l,2$+1'
a=1
Finally
IR(h,k;28, 2s)| = lR(h,k; 2s+l, 2s+l)l
< 5.317(2n)“2‘*2’llb1|N2.+2.o
5 (4.4.9)
+ Z Ib2aW2s+2,2a-1 + N25+2,2s+2]' D
a=1
Theorem 4.4.2

If there exists a constant C such that for 0 < a. < 28, 0 < k < h < 71, M23+2’a < C, and

M ’25....2 < Cthen

a

IR (h, k; 23, 23)| < ll.17(b- a)(d- c)C(h/2n)2s+2. (4.4.10)

Here i 1.17 approximates

(4.4.11)

__,,4\/g[1 + 6"/2 + —————h/6 ]
45 1 - (h/27r)2

48

Theorem 4.4.3

If there is a constant C such that for 0< a < 25, 0 < k< h < 1r/C, M23” a < C25”, and

M 25.4.2 < CZS+a, then

a!

lR(h.k;23,2s)| = lR(h.k;25+l,28+1)|

25+2
(1061(b_ a)(d- c)(h_£) [LL/8+ hC/l2 1

12 1 4126721024 (4'4‘12)

where

10.61 z 2114 \/6‘/45. (4.4.13)
Proof:

|R(h,k; 2s+1, 2s+1)|

 

<47r4\/6- 2" 3., a hC l
(b- aXd- 0(2), ) 22C2 2[lb1 1"“;lb2a'0702 1+(27r)2 2776‘]

 

1r4 25+2 5
<2 4\/6_(b_ a)(d- c)(h_C) [i+ z_(2 fl)-2a(hc)2a-l +\_/_6]

 

 

2 12 (4.4.14)
a=1
2145/— (b a)(d )(hc)2‘+2 [6 +~/6_ + «(hey -(hC/2n)2‘]
21r 12 6 211/1—(hC/2‘rr)2 °
4.5 A NUMERICAL EXAMPLE
The Euler-Maclaurin Summation Formula (4.1.13) applied to the integral
1 l dxdy
[(f) = J —— = 1n (27/16) = 0.523 248144 (4.5.1)
0 0 x +y + 1

results in the formulation

[(f) - R(h,h; 23, 25)

<1>(h,h;2s,2s)
"Znin 1
h2 _—
'+' +
i=0 i=0 0 ”h l

-Z=:b2ah2“[2h(2a-l)! Z {(1 +511)?“ - (2+Bh)'2“}

1 6= -—0

(4.5 .2)

a-l

+ Z 2b26h2"(2a + 20 -2)! {41‘(a*‘” — 314W“) - 1}
=1

13

+ 152,,112“(4a.-2)![41'2a - 3"4a - 1}]

49

where

h=—1—, n= 1,2,
r1

= (-1)a+fi(a+11)! .
(x +y + 11W”

 

fa‘3(x,y) (4.5.3)

Applying the error estimates (4.2.3), (4.3.3), and (4.4.1) we find

h 25 S 2 1 ’1 ZS ]
. __ a- __ 1 l
lR(h.h,2s,2s)| < 295(2fl) [(291 + 22:] Ibzalh (2s+2a 1)- + “(fin (43).

(4.5.4)
’1 ZS+I S
lR(h,h;2.s,2s)l < 2.94;) [(23+1)! + ZZlbzalhza'l(23+2a)!
a=1
h 25,, (4.5.5)
h 23+2 S
|R(h,h;Zs+ 1,2s+1)| < 5.31(2—-> [(2s+2)! + 2Z162a|(2s+2a+1)1 (4.5.6)
" a=1
2s
+ 5%???) (4s+2)!]

The results are presented in Tables 4.5.1-4.5.7. The partial derivative correction sums are given in
Table 4.5.1. Table 4.5.2 gives the results of the Euler-Maclaurin Summation formula for several values
of h and s; the associated errors are listed in Table 4.5.3.

Tables 4.5 .4 through 4.5.6 present the results of applying the error estimates (4.5.4) through
(4.5.6), respectively. Finally, in Table 4.5.7, we give the absolute value of the ratio (Error Estimate
4.5.k)/(Actua1 Error), k = 4,5,6. The results indicate that the choice of the error estimate depends
not only on the integrand f but also on the values of h and s.

For values of s < l/h, (4.5.6) provides the sharpest error estimate while for s > l/h, (4.5.4) should
be used. For 8 z l/h, (4.5.5) provides the best estimate of the truncation error.

In practice, one would fix the value of s and let )1 decrease to zero. In this case, error estimate

(4.5.6) should be applied.

50

O
90

4
4.444444444444444.

O

4
e 0
044444

O
O

O

O

O

O O 4

9

000
9000

O
O
006

e
0004444444444

....
4
4
4

O

4

O

4
o

04

4

4
4

4
4
4

4

4

4
4

44

4

4
4

44

4
4

4
4

4
4
4
4

4
4
4

4
4
4

4
4
4

4
4
4

4
4
4

 

4
4
4

4
4
4

 

4
4
4

4
4
4

4
4
4

4
4

4

4
4
44

4.

4

:4

4

Y

 

4
4
4
4

4
4
4

 

 

4
4
4
4

4
4
4
4

 

 

..
.
.
4
..
4

4
4
4

 

4
44

 

 

 

 

 

 

 

 

 

 

= (x +y + 1)"or110.112

Figure 4.5.1 Graph of z

51

Table 4.5.1 Partial Derivative Correction Sums <i>(h,h; 2s, 2s)

 

 

 

 

 

 

 

 

 

s
h l 2 3 4 5
l .063 143 004 .058 775 925 .058 555 694 .018 202 954 -2.l47 609 63
1/2 .014 501 993 .014 231 639 .014 255 757 .014 250 457 .014 251 156
1/3 .006 301 135 .006 247 312 .006 249 366 .006 249 187 .006 249 216
1/4 .003 513 700 .003 496 638 .003 497 003 .003 496 985
1/5 002 239 387 .002 232 394 .002 232 490
l/lO .005 566 435 .005 562 062

 

Table 4.5.2 Euler-Maclaurin Summation Formula: I( f ) = 0.523 248 144 z T(h,h) - <1>(h,h; 23, 2s)

 

 

 

 

 

 

 

 

 

 

 

 

 

 

h s 0 1 2 3 4 5
1 .583 333 333 .529 190 329 .524_ 557 409 .525 777 640 .585 130 379 3,730 942 96
1/2 527 500 000 .52; 998 007 .523 268 361 .523 241843 .523 24g 543 .523 248 844
1/3 .52_9_ 497 354 .523 126 230 .523 259 043 .523 247 289 .523 248 197 .523 248 138
1/4 .529 745 130 .523 281430 .523 2485192 .523 248 1_2_7 .523 248 14_5_
1/5 .52_5_ 480 630 .523 24; 243 .523 248 _2_36 .523 248 149
1/10 .523 804 351 .523 247 108 .523 248 148
Table 4.5.3 Error R(h,h;2s, 28) = [(f) — T(h,h) + <I>(h,h;29,28)
S
h 0 1 2 3 4 5
1 -6.01-2* 3.06-3 —1.31—3 -1.53-3 -4.19-2 -2.21+0
1/2 -1.43-2 2.50-4 _2.02-5 330-6 -l.40—6 4.004
1/3 —6 .25—3 5.19-5 -1 .90—6 155-7 —232-8 520-9
1/4 -3.50—3 1.67-5 -3.49-7 1.65-8 -1.40-9
1/5 -2.23—3 6.90-6 —9.24-8 330-9
1/10 -5.56-4 4.35-7 -1.80—9

 

 

 

 

 

 

 

‘-6.01-2 means —6.01 x 10-2.

 

 

 

 

52

Table 4.5.4 Estimates for R(h,h; 2s, 23) using (4.5.4)

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

s
h I 2 3 4 5
l 2.57-1 1.45-1 ' 5.23—1 1.1 1+1 7.50+2
1/2 4.87-2 436-3 107-3 6.57-4 1.34-3
1/3 1.97-2 732-4 693-5 1 .26-5 3.99-6
1/4 1.06—2 2.16-4 1.11-5 1.07-6 1.67-7
l/S 6.61-3 8.50-5 2.75-6 1.66-7 1.61-8
1/10 1.57-3 4.91—6 3.85-8 563—10 132—11
Table 4.5.5 Estimates for R(h,h; 2s, 23) using (4.5.5)
h s 1 2 3 4 5
1 1 .44-1 1 .88-1 1 .67+0 6.77+1 7.26+3
1/2 123-2 190-3 7.10—4 7.10—4 2.63-3
1/3 3.26-3 2.04-4 2.74-5 6.5 5—6 2.73-6
1/4 1 .3 0-3 4 .46-5 3.21-6 4.00-7 7.73-8
1/5 647-4 139-5 631-7 491-8 585-9
1/10 759-5 397-7 436-9 8.19-11 2.35-12
Table 4.5.6 Estimates for R(h,h; 2s, 2s) using (4.5.6)
h S 1 2 3 4 5
1 1.75-1 3.03-1 2.62+0 9.83+1 1.00+4
1/2 7.63-3 1.78—3 9.07-4 1.11—3 4.31—3
1/3 1.32-3 1.24-4 2.23—5 6.81-6 3.53-6
1/4 392-4 200-5 192-6 301-7 7 .04-8
1/5 154-4 495-6 299-7 291-8 4.17-9
1/10 8.87-6 6.94-8 1.01-9 238-11 822-13

 

 

 

 

 

 

 

 

 

 

53

Table 4.5.7 I(Error Estimate 4.3.k)/Errorl, k = 4,5,6

 

 

 

h 1 2 3 4 5

1 84 41 57 i1 144 231 8Q]0921712 2_6_51616 2346 38932844575
1/2 195 49 21 216 94 88 324 _2_1_§ 275 192 507 793 £13 3757 6157
1/3 380 63 _25 385107 65 447 177 £4 543 _2_§_2_ 294 767 5g 679
1/4 635 78 .22 619128 _51 673 195 M 764 286 2.15 835 387 222.
1/5 958 94 gg 920150 21 833 191 2;

1/10 3609174 30 2728221 12

 

 

 

 

 

 

 

 

5. APPLICATIONS OF THE 2-DIMENSIONAL EULER-MACLAURIN
SUMMATION FORMULA

5.1 CUBATURE FORMULAS WITH ASYMPTOTIC EXPANSIONS

The generalization of the quadrature formulas in Section 3.1 to functions of two variables involves
expressing the Euler-Maclaurin Summation formula (4.1.15) in terms of a rectangular grid and Taking
appropriate weighted trapezoidal sums for various grid sizes in order to eliminate the desired number
of terms in the asymptotic error expansion. The result is a cubature formula of the required degree of
precision. The details of the technique are as follows.

Let I <al<a2 < - ' -<a3+l be factors ofn = (b —-a)/h andm =(d-c)/k. For 1 <1,

i < s + 1, define the n + l = (s+l)(s+2)/2 trapezoidal sums

where
1/2 if i =j
" { 1 otherwise ( )

and T(h,k) is defined by (4.1.4). For convenience, we define E5” = 530 and [111,3 = 505.
Denote by

A _ A(a1.al1012.911)":(as.a1)(a2.a2)'“(as.as)

— (2511'19 2512'1)(21321-l. 2522‘l)-"(23m -l, 23n2-1)(h1k) (5.1.3)

the approximation to the double integral

d b
1(f)=] ] f(x.y)dxdy (544)
C 0

based on the weighted average of the n sums Tara," 1 < 1', j < s + l which eliminates the 1) terms
involving 5251141 231,24, 1 < 1' < 17, but possibly has some loweroorder partial derivative correction

terms, 52711-1, 27i2-1’ l g i < I.

Let

8,]. = { " (5.1.5)

2B,)- otherwise

and

54

55

iii: 5' 6. 5- 6-
C. . = >\ ' a 111a 122 + ”I 122
1112 1112 j, ['2 “1'2 air ' (5.1.6)

We wish to find 77 + 1 constants xi}- such that

n+1
A = Z xiiTal-aj = Hf) + CbzubQOzu-i,2v-1 (5-1'7)
i>i=l

where C is some constant, u and v are appropriate nonnegative integers, and Q is given by

n+1
Q2u-l,2v-I = Z xrrEzu-1,2v-1(arh»07'k)- (5-1-8)
1>i=1

If we write “’1 = x“, wz = x21, ° ' ' , wnfl = xs+1’s+l, then the 17+] constants w,- may possibly

be found by solving the linear system

 

f2)“ 27‘21 2’81 2‘22 2282 ”my: FWIT P1-
11 11 11 11 11 11
611 C21 Csl C22 Csz Cs+1,s+1 “’2 0
22 22 22 22 22 22 _
C11 C21 C81 C22 ’ ' ' C32 ' ' ' Cs+l,s+1 W3 ' 0 ' (5.1-9)
1111 1m
_CH C2212 ' ' ' C21" C22 ' ' ' C372" ' ' ° Cs+l.s+1_ _wrfil‘ _Oj

 

 

 

 

 

If a solution to (5.1.9) exists, then for some integer a we obtain an asymptotic expansion for the

cubature formula, A:

s+1 t
A = Z xi/Taiai ’ szmbhzz 92711-142712:l
i>j=l 1=1

(5.1.10)

[(0 + Q(h,k; 00).

A“ .1)(2.l)(2.2)(4.l)(4.2)(4.4)

Note that in some cases, for example, (1,0)(3,0)(3,1)(3,3)(5,o)

(h ,k), correction terms of the

form E a0 may not be eliminated.

Recalling (4.1.15), the definition of R(h, k; 2s, 2s), we write

Raiai = Ail-[R(aih,ajk;o,o) + R(afh,a,-k; (7,0)], (5.1.11)

56

Then the truncation error 52 in (5.1.10) is

S
902.1640) =. Z x.,R.,., (5.1.12)
i>f=l

and may be estimated by the methods of Chapter 4.

To illustrate the technique we will derive the partial derivative corrected (DC) Simpson rule:

. _ (1.1)(2.1)(2,2)
S(h,k)—A(3’O)(3,l) (h,k). (5.1.13)

Settings= 1, 17:2, (a1,a1)=(1,1), (a2,a1)=(2,1), (aQ,a2)=(2,2), (811,612)=(2,0),
(1321,22)=(2.1), (711.712): (1.0). and(721,722)= (1,1) we obtain the system
.— -[ " " T

1 18 32 w2 = 0 (5.1.14)

 

 

 

 

 

 

1 20 64Lw3 0

having the unique solution

(wl, W2, W3) = (256/225, -16/225,1/225). (5.1.15)
Thus we obtain (5.2.4), an asymptotic expansion for the DC Simpson’s rule. The error is given by
Q(h,k;2s, 23) = [256R11-16R21+R22]/225. (5.1.16)

This cubature rule is illustrated in Figure 5.2.5 and the trapezoidal sums are shown in Figure
5.2.4.

In Section 5.2 we give without comment several cubature formulas with asymptotic expansions,
error terms, and appropriate diagrams. These formulas are the 2-dimensional generalizations of the
quadrature formulas of Section 3.1. Some of the results, e.g., theDC Weddle’s rule, are believed to be
new. Sheppard [45] obtained the double Simpson’s and Weddle’s rules. A number of additional
cubature rules were obtained but are not presented here.

Using another technique, Tanimoto [50] derived what we call the DC Simpson’s rule, (5.2.4).

However, his paper has several errors. Moreover, he does not give an expression for the error.

57
Finally we note that

(l .l)(2.1)(2.2)(4.1)(4.2)(4.4)
(l ,0)(l.1)(3.0)(3,l)(5.1)

results in the linear system

~_—

"1 2 1 2 2 1 w,
1 6 8 20 48 64 w2
1 8 16 32 128 256 W3

1 18 32 260 576 1024 w,

l 20 64 272 1280 4096 “’5

 

 

 

1 68 256 4112 17408 65536 W6
.1 _1

 

in which the coefficient matrix is singular.1 In this case there are infinitely many solutions.

5.2 SOME CUBATURE FORMULAS OBTAINED FROM THE
EULER-MACLAURIN SUMMATION FORMULA

Euler-Maclaurin Summation Formula: Trapezoidal Rule

m n
T(h,k) _=_ th Z f(a+ih,c+jk)

i=0 i=0
d b E10 E11
= , dxd +— _—
fcfaflx”) y 12 144

_flg_531 _ E33

720 8640 518 400

E50 551 553

 

 

 

 

+ + — -
30 240 362 880 21 772 800 914 457 600

E70 E71 E73

 

 

_ _ + _ _________.
1209 600 14 515 200 870 912 000 36 578 304 000

(5.1.17)

(5.1.18)

(5.2.1)

(Equation (5 .2.1) continues)

1 . . . . . .
If the 1 on the right Side of (5.1.18) 15 changed to 2025, then a solution 18 given by (4096, -1280, 400, 64,
‘20. I). This system proved difficult to evaluate numerically. A modified version of the Fortran program recommended

by Forsythe and Molcr [17] was used.

E77

58

E90 E91 E93

 

" 1 463 032 160 000

E 95

 

+ —
1 448 500 838 400

 

+ + ——-———-— -
47 900 I60 574 801 920 34 488 115 200

E97 E99

 

 

57 940 033 536 000 - 2 294 425 328 025 000

S
- + b2s [En-1,0 + szeEzs-ma-I] +R(h,k§2‘12‘)'

[i=1
(5.2.1)

 

)4

5.4

 

%

th

)4

 

 

Figure 5.2.1 Trapezoidal Rule

DC T rapezoidal Rule

E10 E11
T(hJC) = T(h,k) -7 + m
E30 531

IU)

 

_ _ E33
720 8640 518 400

+ (5.2.2)

S
+ b25[E2S'1,0 + Zb23E25-1926-1] + R(h,k; 23, 23)‘

{i=1

59

 

 

 

 

 

 

 

 

 

 

-hkfxy -kfy -kfy hkfxy
144 12 12 144
hfx 41)“,
E V: V: 12
x hk

hf, -hfx
_1'2— V“ y‘ 12
hkfxy kfy kfy 41ka
144 7 1_2‘ 144

Figure 5.2.2 DC Trapezoidal Rule

Simpson ’s Rule

S(h,k) =

(1.1)(2.1)(2.2)
A(1.0)(1.1) 0””

[16 T11 -4 T21 '1‘ T221/9

(5.2.3)
Q30 533 Qso E53 I555

[(f) - — - + - .. + . . .
720 32 400 30 240 272160 2286144

 

 

s
l
+ b25[Q2s-1,0 + 32 b23[16’4(4s +4”) + 4S+BIE23-1,25-1] + 5201,16; 23. 23)-
B=l

60

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

1 4 1
l6
4 4 X ’31
9
1 4 I
Figure 5.2.3 Simpson’s Rule
V: 16 ’24 l I
1
‘A % 1
V: V2 V: 1 l
T“ = T(h,k) T21 = T(2h,k)+T(h,2k)
1 l
1
T22 = T (211.210

Figure 5.2.4 Component Trapezoidal Sums for Simpson’s Rule

61

DC Simpson 's Rule

(3.0)(3.1)
=[256T 16T +T 1/225 g9+i
950 E51 E55

 

= + - -
IU) 30 240 ' 141 750 89 302 500

S
l
+ b23[Q23_1,0 + 2752 1’23 [256 -16(4‘ + 43) + 4S+B]E23,1,2B_1]
a=1

 

 

 

 

 

 

+ Q(h,k; 23, 2s).
—hkfxy -7kfy -l6kfy ~7kfy
7hfx 49 l 12 49
1 12 256 1 12
l6hfx
7h)“, 49 112 49
hkfxy 7kfy 16kfy 7kfy

Figure 5.2.5 DC Simpson’s Rule

hkfxy

-7hfx

-16th

-7th

—hkfxy

.44..
225

(5 .2.4)

62

Simpson's Second Rule

 

 

 

 

 

 

 

 

 

 

(LOXIJ)
Q E Q E E
=IU)- 30__ 33+ 50_ 53_ 55 +...
720 6400 30 240 26 880 112 896
(5.2.5)
1 S
+ b25[Q25-1,0 + g; 2528181 '9(9‘ *9”) + 93+6]EZS-1,213—l]
a=1
+ Q(h,k;2s,2$).
1 3 3 l
9 9
3 3
9hk
(,4
9 9
3 3
I 3 3 1

Figure 5.2.6 Simpson’s Second Rule

1/4 1/2 1/2

 

 

 

 

 

 

 

 

I 1

1/2
1 1

1/2
1/4 1/2 1/2
T11 = T(h,k)

3/2

3/2

1/4

1/2

1/2

1/4

3/2

3/2

3/2

3/2

3/2 3/2

 

 

 

 

 

 

 

 

3/2 3/2

T31 = T(3h,k)+T(h,3k)

 

3/2

 

 

 

 

 

 

 

3/2

T33 = T(3h,3k)

Figure 5.2.7 Component Trapezoidal Sums for Simpson’s Second Rule

3/2

3/2

3/2

3/2

DC Simpson '8 Second Rule

Q10
A(1’1)(3’1)(3’3)(h,k) = [65617111 -31131+ T33]/6400 _ __ +

(3.0)(3.1)

4111:ny

26hfx

54h f,

54h f,

26th

411kny

12

950 9E 51 9555

 

[(f)

+ - — -————+
30 240 44800 1254 400

95”
1600

 

S
l
13:1

+ 9$+61 E25._1’2B_1] + S2(h,k,2S,25).

 

 

 

 

 

 

 

 

—26kfy —54kfy -54kfy -26kfy
169 351 351 169
351 729 729 351
351 729 729 351
169 351 351 169

26kfy 54kfy 54kfy 26kfy

Figure 5.2.8 DC Simpson’s Second Rule

41:1:ny

-26hfx

—54th

-54hfx

-26th

-4h kfxy

(5.2.6)

6—400

65

5 2 Point Rule

(1.1)(4.1)(4.4) _
Au’om’” (h,k) — [256T11-16T41 + 7441/225

Q3O E33 17553 289555 +
" 720 ' 2025 -‘ 85 050 ' 3 572100

 

(5.2.7)
1 S
(i=1

+ 16s+6152,_1,25-1] + S2(h,k;25,2s).

 

 

 

 

 

 

 

 

 

 

9/8 3 3 3 9/8
8 8 8
3 3
8 8 8 32hk
3 x _—
225
8 8 8
3 3
9/8 3 3 3 9/8

Figure 5.2.9 52 Point Rule

66

 

 

 

 

 

 

 

 

 

 

 

 

DC 52 Rule
0 165
(1.1)(4.1)(4.4) _ ~10 11
A(3’0)(3’1) (h,k) - [65 536T11- 256 T4] + T44]/65 025 -—12 + _2601
Q50 32E51 645 5
= [(f) + _ _ ____.5___ +
30 240 819 315 258 084 225
(5.2.8)
S
+ b Q + -——1—- b [65 536 256 (168 +163
2s 2s-1,0 65 025 Z 213 - )
+ 163+B]E2s_1,23_1] + Q(h,k;28,2S).
25 315 315 25
— ? hkfxy - 64 kfy -10kfy -10kfy -10kfy — H kfy .8. hkfxy
315 315
—- hf 3969 63 63 63 3969 - — hfx
64 x —— — 64
128 128
63 128 128 128 63
10th -10th
x 512hk
65 025
63 128 128 128 63
10th -10hfx
63 128 128 128 63
10th -10th
315 3969 3969 315
_ h __ __ _—
64 fx 128 63 63 63 128 64 hfx
25 315 315 25
—8—hkfxy a kfy lOkf), lOkfy IOkfy 64— kfy --8- hkfxy

Figure 5.2.10 DC 52 Rule

67
Boole '8 Rule

__. (l,1)(2.1)(2.2)(4,1)(4,2)(4,4)
R(h’k) A(1.0)(l.1)(3.0)(3,1)(3.3) (h’k)

[4096 T11 - 1280 721 + 400 T22 + 64 T4, - 20 T42 + T44] /2025

[256511-16512 +S22]/225

+ Q50 4E55 Q70 E75 E77 +
30 240 ' 893 025 ' 1 209 600 ' 425 250 ' 810000

 

 

 

[(f)

S
+ b2$[Q2s-1,0 + $52 bZB [4096 ‘ 1280913 + 46) + 400W”)
5:1

h

+ 64(163+16B)-20(425+5+4S+23) + 16S+5]E25-1,23;1] + 52(h.k;2s,231.
(5.2.9)

 

 

 

 

 

 

 

 

 

 

49 224 84 224 49
1024 384 1024
224 224
84 384 144 384 W,
2025
1024 384 1024
224 7 224
49 224 84 224 49

Figure 5.2.11 Boole’s Rule

V;

1A

‘A

V;

 

 

 

 

 

 

 

 

 

 

V;

 

 

 

 

 

 

 

 

 

 

% 1/2 ‘2’2
1 1 1
l 1 l
1 1 I

v. 14 V2
T11 = T(h,k)

l 2 1

2
2 4 2
2
1 2 1

T21 = T(2h,k) + T(h,2k)

2

 

 

 

 

 

 

 

 

 

 

2

T22 = T(2h,2k)

V1

68

 

 

 

 

 

 

 

 

 

2 2 2
T41 = T(4h,k)+ T(h,4k)

4

 

 

 

 

 

 

 

 

 

4
T42 = 7(4h,2k) + T(2h,4k)

 

 

 

 

 

 

 

 

 

T 44 = T(4h,4k)

Figure 5.2.12 Component Trapezoidal Sums for Boole’s Rule

 

 

 

69

DC Boole’s Rule

A“ .l)(2.1)(2.2)(4.1)(4.2)(4.4)

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Figure 5.2.13 DC Boole’s Rule

3'0““ ._. (3.0)(3.1)(5.0)(3,3)(5.1) 0"")
= [1048 576711 —81 920721 +64007‘22 + 102471,,- 80T42
Q10 1615‘11
+ T 893 025 -— + —.——
441/ 12 3969
= ' .. ' + '
[4096s11 64s21 $22]/3969
Q70 167571 16577 (5.2.10)
= [(f) - —————- - - +
1 209 600 6 251 175 9 845 600 625
l S
+ b + —— b 1 048 576
23[Q2$-1,0 893 025 g1 2131
— 81920(4s +4?) + 6400(4W) + 1024(16s + 165)
— 80(423+5 + 43*”) + 16W] 52,,1,2,,_1]+ 52(h,k;2s,2s) s = 4,5,---
225
_2—2§ kfxy ”gig/cf). -960kfy -810kfy ~960kfy —E:—;9kfy T-hkfxy
6_5i0hf 47 089 6 944 5 859 6 944 47 089 _6510hf
16 x 16 16 ’1? x
6944 16 384 13 824 16 384 6944
960th —960hfx
5859 13 824 11 664 13 824 5859
810th -810th 896:2;
6944 16 384 13 824 16 384 6944
960th -960hfx
6510“ 47 089 47089 _6510U
16 x 16 6 944 5 859 6 944 16 16 'x
225 6510 , , 6510 -225
—4-hkfxy —l—6~kfy 960ka 810ka 960kfy To, Thkfxy

70

Weddle ’8 Rule

(1.1)(2.1)(2.2)(3.1)(3,2)(3,3)
Au,0)(1,1)(3.0)(3,1)(3.3) 0"")

= [225T11-90T21+36T22 +15T31-6T32 + T33l/100

Q50 Q70 E75 577
[(f) + - _ _ _ __._ + . . .
30 240 l 209 600 2 016 000 5 760 000

 

(5.2.11)

S
+ b2.[Q25-1,0 + 1002b28l225 -90(4s +4?) + 36(4S+B) + 15(93 +96)
a=1

— 6(9‘43 + 4599) + 9s+9152,,1,2,,_1] + S2(h,k;2s,23).

 

 

 

 

 

 

 

 

 

 

 

 

 

 

1 5 1 6 1 5 1
25 5 30 5 25
5 5
5 1 6 1 5
1 1
30 6 36 6 30 9hk
6 100
5 1 6 1 5
1 1
25 s 30 5 25
5 5
1 5 ‘ 1 6 1 5 1

Figure 5.2.14 Weddle’s Rule

71

DC Weddle’s Rule

(1 .1)(2,1)(2.2)(3,1)(3.2)(3,3)
A(3.0)(3.1)(3.3)(5,0)(5,1) 0"")

 

 

Q10 9511
= [72 900 Tn - 7290 721+ 729 T22 + 540 T31 - 54 T32 + 4 T33]/60 025 "17— + 240T
Q70 45571 9577
= [(f) _ - _ + . . .
1 209 600 6 722 800 18 823 840 000 (52.12)

S
+ 62, [023-1 ,0 + 6—01025 Z 62, [72 900 - 7290(46 + 46) + 729(4s+6) + 540(96 + 96)
a=1

- 54(9646 + 4696) + 95113119234334] + 52(h,k; 2s, 2s).

 

 

 

 

 

 

 

 

 

 

 

 

 

 

25711:ny -185kfy 450M, -360kfy 4601(7), ~360kfy 45074,, 4851.], 25hkfxy
185th 1369 3 330 2 664 3 404 2 664 3 330 1369 485th
3330 8 100 6 480 8 280 6 480 8 100 3330
45014,, 45072;,
2664 6 480 5 184 6 624 5 184 6 480 2664
360hfx 660th
X 911k
3404 8 280 6 624 8 464 6 624 8 280 3404 60025
460hfx -460th
2664 6 480 5 184 6 624 5 184 6 480 2664
360th 860th
3330 8 100 6 480 8 280 6 480 8 100 3330
450hfx 450117,,
1369 3 330 2 664 3 404 2 664 3 330 1369
185th 485th
251:1:ny 1851(7), 450kfy 360kfy 460w), 3601(7), 450kfy 185kfy 4511ka

Figure 5.2.15 DC Weddle’s Rule

 

Newton-Cores’ 72 Rule

N(ll,k)

=A

(l,0)(l.l)(3.0)(3.l)(3,3)(5,0)(5.1)(5,3)(5.5)

(1,1)(2.1)(2.2)(3.1)(3.2)(3.3)(6.1)(6.2)(6.3)(6,6)(h k)

— 1296 761+ 567 762 -112 763 + 7661/705 600

+ b25|:Q2$-1,O +

+ 145 152(95 + 96) - 63 504(9646 + 4696) + 12 544(967‘6) - 1296(366 + 366)

Q70

9577

Q90

E97

 

 

 

 

10)-

l
705 600

_ + _ _ +
1 209 600 7 840 000 47 900 160 689 920 13 660 416

[1679 616 T11 - 734 832 721+ 321489 722 +145152 T31 - 63 504 732 +12 544 T33

(5.2.13)

S
Z62, [1 679 616 - 734 832(4‘ + 46) + 321(4S+6)
a=1

+ 567(36646 + 46366) — “2069‘3 + 96366) + 365+31E2H 213-1] + S2(h,k;2s, 2s).

1681

8856

1107

11152

1107

8856

1681

8856

1107

11152

1107

8856

 

46

656

5 832

58 752

5 832

46 656

 

832

729

7 344

729

5 832

 

58

752

7 344

73 984

7 344

58 752

 

832

729

7 344

729

5 832

 

46

656

5 832

58 752

5 832

46 656

 

 

 

 

 

 

 

 

 

8856

1107

11152

Figure 5.2.16 Newton-Cotes’72 Rule

1107

8856

1681

8856

1107

11152

1107

8856

1681

hk
19 600

73

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

56 ‘6 56. ’72 1/2 92 ‘24:
l 1 1 l I
1/2 '72
1 l 1 1 1
1/2 %
1 1 1 1 1
‘22 52
l 1 1 1 l
15 %
1 l 1 1 1
ti 14
’74 ‘A ‘A 1/z '72 $2 ’74
T“ = T(h,k)
l 1 2 l 2 1 l
2 2
l l
2 4 2 4 2
2 2
2 2
l 1
2 4 2 4 2
2 2
2 2
l 1
1 1 2 1 2 1 1

72, = T(2h,k)+T(h.2k)

Figure 5.2.17 Component Trapezoidal Sums for Newton-Cotes’ 72 Rule

74

 

 

 

 

 

 

 

 

 

 

 

 

 

 

1 2 2
4 4
2
4 4
2
1 2 2

7“22 = 7(2):, 2k)

 

 

 

 

 

 

 

 

 

 

 

 

 

 

3/2 3/2 3/2 3 3/2 3/2
3
3/2
3
3/2
3 3 6 3 3
3
3
3/2
3
3/2
3/2 3/2 3/2 3 3/2 3/2

T31 = T(3h,k) ‘1’ T(h,3k)

Figure 5.2.17 (cont’d.)

3/2

3/2

3/2

3/2

3/2

3/2

75

 

 

 

 

 

 

 

 

 

 

 

 

 

 

3 3 3 3 3
6

3 3

6 6

3 3
6

3 3

3 3 3 3 3

T32 = T(3h, 2k) + T(Zh,3k)

 

 

 

 

 

 

 

 

 

 

 

 

 

 

9/4 9/2 9/4
9

9/2 9/2

9/4 9/2 9/4

733 = T(3h,3k)

Figure 5.2.17 (cont’d.)

76

 

 

 

 

 

 

 

 

 

 

 

 

 

 

3 3 3 3 3
761 = T(6h,k)+T(h,6k)

 

 

 

 

 

 

 

 

 

 

 

 

 

 

6 6 6
6
6
6 6 6

762 = T(6h,2k)+T(2h,6k)

Figure 5.2.17 (cont’d.)

77

 

 

 

 

 

 

 

 

 

 

 

 

 

 

9 9 9
9 9
9 9 9

T63 = T(Oh,3k) '1' T(3h,6k)

 

 

 

 

 

 

 

 

 

 

 

 

 

 

766 = T(6h,6k)

Figure 5.2.17 (cont’d.)

78
Dc Newton-Cotes’ 72 Rule

(1.l)(2,1)(2,2)(3.l)(3.2)(3.3)(6.1)(6.2)(6.3)(6.6)(}2 k)
(3.0)(3.l)(3.3)(5 ,0)(5,1)(5 ,3)(5.5)(7,0)(7 ,3) ’

N'(h,k) = A
= [2176 782 336 711 - 238 085 568 721
+ 26 040 609 722 + 20 901 888 7,1 — 2286144732 + 200 704 T33

— 46 656761 + 5103 762 — 448T63 + 7661/1 764 000000

Q10 + 9511
12 2500
Q90 9591 9E99

(5.2.14)

 

 

= 1U) + - - +
47 900 160 7 700 000 23 716 000 000

1
+ __.—__.—
b2‘[Q2"160 1 764 000 000

+

S
Z 62,, [2 176 782 336 — 238 085 568(46 + 46)
a=1

+

26 040 609(46+6) + 20 901 888(96 + 96) — 2 286 144(9646 + 4696)

+

200 704(96+6) — 46 656(366 + 366) + 5103(36646 + 46366)

448(36‘9" ’6 9’36") ’6 36”“1525-126-1]

+

Q(h,k;23,2s).

9

7

coo coo av
63

 

x

85.88 2-

$88 82-

$88 68.

$82 68.

U$88 68-

982 88.

$88 88-

U$88 82-

.3888 :

.588 82

.588 88

.382 m8

988 68

.582 “8

use 1 688.8252 8. 2.88 28m

.988 8m

9888 8:

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

838 88 $8 88 2.: 88 28 88 :2 88 88 .838
8888 6:6 88 88 8:. 88 886 88 8:. 6:6 88 8838
83:2 88 8:. 88 88 88 83. 88 88 88 8:. 85:2
88 28 88 88 88 88 88 88 88 88 88 88 8828
8322 88... 8:. 88 88 88 88 88 88 88 8:. 88:2
8888 83 $8 88 8:. 88 88 88 8:. 8; 88 883.8
838 88 3.8 88 :2 88 $8 83 :2 88 88 8.38

$8882. $8888- $8288- $8868- $8288. .988 88.

8.388% :

R388 82

$88 88

4.882 m8

$88 68

982 88

$88 68

$88 82

$88 82. 55.88 :-

80

Romberg's 9 2-Point Rule

A“ ,1)(2,1)(2,2)(4,l )(4.2)(4,4)(8.l)(8,2)(8,4)(8.8)(h k)
(l ,0)(l ,1)(3,0)(3,1)(3,3)(5,0)(5,1)(5.3)(5.5) ’

[16777 216TH - 5505 024T21 + 1806336T22

+ 3440647‘41 -112896T42 + 7056T44

— 4096T81 + 1344 T82 — 84T84 + T88]/8 037 225
= [40961211 — 64821 + 8221/3969

Q7O 256E77 Q90 2176E97 18 4961599
.. + _ _ ...
1 209 600 22 325 625 47 900160 88 409 475 350 101 521

 

 

 

10)-

(5.2.15)

4.

S
1
b25[Q25-l,0 + WEI)” [16 777 216 — 5505(4S +44)

+

1806 336mm) + 344 064(16‘ + 165) - 112 896(16‘4fi + 45163)

+

7056065”) - 4096(643 + 643) + 1344(64545 +4564‘3)

84(64316‘i +16564‘3) + 645+B]Ezs_1,23,1]

+

”(hvk; 25) 25)‘

81

gm :6 m

923—

x

amok. v

QONN NN

vwme s

mONN mm

chm

womm NN

vwmo N.

wONN NN

one» v

22% “Eomdm MwSQEoM

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

wONN NN vane N. momm mm 23. a wo- mm vwmc N. momN NN
2.3 X: $56 on 053 X: wove 3 2.3 3: 336 cm 2.3 2.:
wvvo on 3mm 2 wvvo on N53.” 2 3.3 cm coon N. 3.3 on
2.3 1: 3.3 em 22 2: $3 3 2.3 X: wvvo on 2.3 X:
vove 3. 23 w— vocw 3. oaoo m— wove 3 2.3 2 93% 3.
3.3 X: 3.3 on 2.3 X: wove 3. 3.3 X: 3.3 on 2.3 X:
3.3 on voam 2 wvvo on 2% n. wvvo cm 33 S 38.0 on
2.3 3: avg on 2.3 2: vote 3. 2.3 we. wvwo on cam vo—
woNN NN vane N. womm NN «_ov a .woNN NN 336 n momm mm

2.3 285

one. v

womm «m

gmo h

wONN mm

«Sea

womm mm

vwmo 5

menu NN

omen v

82

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

‘54: 1/’z ‘A 95 1/z lé IA %

l 1 l l 1 1 l
‘A

1 l l l 1 1 1
‘24

l 1 l l l l 1
‘A

l 1 1 l 1 l l
56

l 1 1 l l 1 1
‘5

l 1 l l 1 l 1
1/2

l 1 l l 1 1 l
1A
‘4 5’2 ‘5 V2 % 1A 'A Vz

Tll = T(h,k)

Figure 5.2.20 Component Trapezoidal Sums for Romberg’s 92 Rule

54

83

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

1 1 2 1 2 1
2 2
1
2 4 2 4 2
2
2 2
1
2 4 2 4 2
2
2 2
1
2 4 2 4 2
2
2 2
1
1 1 2 1 2 1

T21 = T(2h,k) + T(h,2k)

Figure 5.2.20 (cont’d.)

84

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

l 2 2
4 4
2
4 4
2
4 4
2
1 2 2

T22 = T(2h,2k)

Figure 5.2.20 (c0nt’d.)

85

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

2 2 2 2 4 2 2 2
4
2
4
2
4
2
4 4 4 8 4 4 4
4
4
2
4
2
4
2
2 2 2 2 4 2 2 2

T41 = T(4h,k) + T(h,4k)

Figure 5.2.20 (cont’d.)

h)

86

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

4 4 8 4
8
4
8 l6 8
8
8
4
4 4 8 4

T42 = T(4h,2k) + T(2h,4k)

Figure 5.2.20 (cont’d.)

87

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

4 8
16

8

4 8

T44 = T(4h, 4k)

Figure 5.2.20 (cont’d.)

88

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

4 4 4 4 4 4

T81 = T(8h,k)+T(h,8k)

Figure 5.2.20 (c0nt’d.)

89

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

8 8 8
8
8
8
8 8 8

T82 = T(8h,2k)+T(2h,8k)

Figure 5.2.20 (cont’d.)

90

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

16 16 16
16 16
16 16 16

T84 = T(8h,4k) + 71(4h,8k)

Figure 5.2.20 (cont’d.)

16

91

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

16

Figure 5.2.20 (cont’d.)

 

T33 = T(8h , 8k)

16

16

92

DC Romberg's 92—P0im Rule

(l.l)(2.1)(2.2)(4.1)(4.2)(4,4)(8,l)(8,2)(8,4)(8,8)(h k)

A13,O)(3.11133115411154)15.3)(5.5117.01<7,1)

[68 719 476 736 Tll - 5 637 144 576 T21 + 462 422 016 T22

+ 88 080 384 T41 - 7 225 344 T42 + 112 896 T44 - 262 144 T81

Q10 125441?“

[65 536/31'l — 256Bé1 + 32'21/65 025

 

 

Q90 8 192591 262 144E99
= I + — - +
(f) 47 900 160 2 027 804 625 63 237 087 230 625
(5.2.16)
5'
+ b _1___
2 Q - + b [68 719476736
5[ 23 1.0 (240 975)2 g1 26

S 637 144 576(45 + 43) + 462 422 016(45+5)

...

88 080 384(16‘ +163) — 7 225 344(16‘43 + 45163)

+

112 896065”) - 262 144(645 + 643)

+

21 504(64543 +4364“) — 336(64516‘3 + 16564“) + 64””1523-126-1]

+

Q(h,k;2s,2s).

5.3 THE MIDPOINT RULE AND VARIOUS OTHER FORMULAS
The double Euler-Maclaurin Summation formula (5.2.1) may be used to derive an asymptotic

expansion for the midpoint or centroid formula C(lr,k) by writing

[1 k h k
4T(—2—,3) - 2[T(3,k) + T( ,9] + T(h,k)

hki ifla + 171/2, 6' + jk/2)

i=1 i=1

d b
fff(x.y)dxdy +

C(Iz,k)

(5.3.1)

93

1f the first- and mixed second-order partial derivative correction terms in (5.3.1) are transposed,

a third-order derivative corrected midpoint formula analogous to (3 .45) is obtained. The resulting DC
midpoint formula is the same as formula EX183S of Table 6.4.1 and consequently the details are not
given here.

The technique described in Section 5.1 may be applied to (5.3.1) to obtain a number of Open
cubature rules including generalizations of the Open quadrature rules of Section 3.4. Because of space
limitations, the results will be omitted. However, we will indicate how to obtain asymptotic expansions
for several nonproduct cubature formulas.

Figures 5.3.1 to 5.3.4 illustrate Squire’s [48], Ewing’s [16] , Tyler’s [51] and Miller’s [35]

cubature rules.

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

1 l 1
hk hk
1 0 0 1 x 0 _
T 8 12
1 1 1
Figure 5.3.1 Squire’s Rule Figure 5.3.2 Ewing’s Rule
1 -1 4 -1
hk hk
1 o o 11 1 x — 0 (1 x __
2 6 4 4 12
l -l 4 -1

Figure 5.3.3 Tyler’s Rule Figure 5.3 .4 Miller’s Rule

94

Asymptotic expansions for these rules may be obtained from the following:
Squire = -T(h,k) + [T(h,k/2) + T(h/2,k)]

Ewing = T(h,k) -—:—[T(h,k/2) + T(h/2,k)] +2—T(h/2,k/2)

1 4 (5.3.2)
Tyler = - ? T(h,k) + —3— T(h/2,k/2)
, 5 4
Miller = —?T(h,k) +?[T(h,k/2) + T(h/2,k)].

These expansions are believed to be new. Another approach to these 4 cubature formulas is via the

Taylor series; this is done in Chapter 6.

5.4 A NUMERICAL EXAMPLE

The results of applying the midpoint and DC midpoint rules to the double integral

1 l
J- J. exydxdy 3 1.317 9021514544 (5.4.1)
0 0

for grid sizes ofh = k = 1/5 and h = k =1/10 are shown in Table 5.4.1.

Table 5.4.1 Several Cubature Rules Applied to fo’fo’ exydxdy =
1.317 902151454 4

 

 

 

 

h=k=1/5 h=k=ll10
Rule
nfe“ Error nfe Error
Midpoint 25 165-31 100 4.16-4
DC Midpoint 49 -l.22-6 144 -7.62-8

 

 

 

 

 

 

 

‘Number of function evaluations
mus means 1.65X10‘3

The DC Simpson’s rule shows a similar improvement over Simpson’s formula. Additional numerical

results are given in Chapter 7.

95

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

1111

 

 

 

 

 

 

 

 

 

 

 

 

 

 

MHLLHJWM‘

Figure 5.4.1 Graph ofz = ex" on [0,1] 2

 

 

 

 

 

 

 

 

 

 

 

 

 

 

6 MULTIDIMENSIONAL QUADRATU RE FORMULAS WITH
PARTIAL DERIVATIVE CORRECTION TERMS

6.1 THE EQUAL WEIGHT-ALTERNATE SIGN PROPERTY
Let f(xl , - ° - ,xN) be a real valued function defined on the symmetrically placed N-dimensional

rectangle

2:
‘=2

II
p—

[411311,], N 2 2.
I

We wish to estimate the multiple integral

hN h,
1(‘f)=‘JI ... f(xl,...,xN)dxl...dx1v (6.1.1)

by a multidimensional quadrature formula with partial derivative correction terms.
Partial derivatives are denoted by
’1’? . . .(jftjiaf(xl’...’xN) (6.1.2)
where ’1'.“ denotes the a-th partial derivative of f with respect to the i-th variable.

Before proceeding we make several observations for the Specific case N = 2.

The study of the Euler-Maclaurin Summation formula for a function of two variables led to an
investigation to search for additional cubature formulas involving first- and perhaps second-order partial
derivative correction terms with weights of equal magnitude and alternate signs at the four corners or at
the midpoints of the sides of the rectangular domain of integration so that when the rule was com-
pounded or repeated, the weights would cancel except on the boundary. We call this the “equal-
weight, alternate-sign property.”

The objective is to select nodes and weights which result in efficient cubature formulas of high
precision. In this connection we observe that for a composite formula of given degree of precision, a
rule in which most of the nodes coincide with the boundary will be more efficient than one in which
almost all of the nodes lie in the interior of the domain of integration. Moreover, additional economy
can be gained by requiring the equal-weight, alternate-sign property.

Because of the endless variety of possible combinations, it was decided to limit the study to the

six “cubature elements” listed in Table 6.1.1.

96

97

Table 6.1.1 Elements

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Name Cubature Element Diagram
Centroid A (f) = 4h h f(O 0)
Value 0 l 2 ’ 0
Corner
Sum ACU) = 4’11h2[f(h1:h2)+f('h1rhz)+f(’h19'h2)+f(h1"h2)l
Mid oint T
Sump AmU) = 4h1h2[f(h1:0)+f(0.h2)+f(‘h1-0)+f(0,'h2)l " - ‘
CDzrnerr Ac:(f) = 4h12hzlfx011’h2) 'fx(‘h1:h2)'fx('h1v’h2)+fx(h1v‘h2)] —+ ++
rivative
2
Correction + hlhz Uy(h1'h2)+fy('hirh2) -fy(-h1,-h2) -fy(h1,-h2)] __ +
$415616: A,;,(f) = 4h 12712140116)- fx(-h1.0)] :
rivative _4 0+
Correction + hlhz2 [6409112) ‘fym' 4'2“ :
Comer ,. 2 2 - +
Mixed Ac U) = 4hikz [fxy(h12h2) "fxy(‘h11h2)
Derivative
- - _ _ + _
Correction + fxy(h1' hz) fxy(h1. ’12)]

 

 

 

 

 

The selection of the derivative correction cubature elements AC', A5,, and AC" is based on the

following considerations. For a and B nonnegative integers, define the “cubature generators” G,-
listed in Table 6.1.2.

Now suppose f(x1,x2) can be expanded in a Taylor Series about the point (0,0) as far as may
be required. If the series converges on R. then the 01' assume the values indicated in Table 6.1.3.

Examination of Table 6.1.3 reveals that odd/even constraints must be placed upon a and B to
avoid nonzero generators Gr" which are components of the partial derivative correction elements A C ,
Ag], and Ac". These cubature elements are then candidates for' inclusion in cubature formulas.

Now in order to construct multiple integration formulas with partial derivative correction terms,
it is necessary to generalize the cubature elements in Table 6.1.1. Consideration of the geometry of the
N-rectangle suggests how to proceed for arbitrary N > 1. We are now ready to derive MINTOV (an

acronym for “Multiple lNTegration, Order 5”).

98

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Table 6.1.2 Generators
Name Cubature Generator Diagram
_ +
a, 914‘; [f(hlJiz) - f(—h1.h2) +r(-h1.-h2) 4011.412)!
.1.. __
G 6 . a. - +
2 v» ,u 1 0011.112) -f(-h1.h2) —f(-h1.-hz)+f(h1,-h2)1 _ 1
+ +
03 9’39? [f(h1.h2)+f('h1.h2) ““411. 412) 'f(h1. 412)]
0. 93911f(h1.01—f(-h1.011 -» 1+
1'
G, 930110682) - f(O. 412)]

 

 

 

 

 

 

 

Table 6.1.3 Generator Values

 

 

 

 

 

 

 

 

 

61,5 a Odd a Odd (1 Even a. Even
Generator [3 Odd [3 Even 6 Odd B Even
GI Nonzero 0 0 0
02 0 Nonzero 0 0
G3 0 0 Nonzero 0
G4 0 Nonzero 0 0
05 0 0 Nonzero 0

 

 

 

99

6.2 DERIVATION OF MINTOV

Denote the vertices of the N-rectangle
N
R = [I [-h,-,h,-]
i =1
by c = (cl, 1 ° - , cN) where cf = —h]- or hf, and the volume by
N
h = 1.11 Zhj.
’3

Define the sign functionals

( ) {-1 if Ci = 41"
a c =
I +1 otherwise
(6.2.1)
ojk(c) = oi(c)ok(c)
and the first. and second-order partial derivative correction terms
0,0) = 2 0,-(c)-"l,-'f(c)
C
(6.2.2)

Dij) = §0jk(c)€/;9,éf(c)

where the sums are over the 2” corner points of R. Figures 6.2.1 and 6.2.2 illustrate two of the
partial derivative correction terms for the 4-cube. Moreover, careful inspection will reveal the sign

arrangements for D1 (f) and Dlz( f ) for the 2- and 3-cubes.

100

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

+

Figure 6.2.1 Sign Arrangement for D1(f)

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Figure 6.2.2 Sign Arrangement for 012(f)

101
Since [0") and the “cubature elements"
A00") = hf(O)
4.0) = 11ch {(c)

.v
A'c(f) = thzl-DIU) (6.2.3)

1' =1
A: (f) = h Z hjthij)
j<k
vanish for functions which are odd in any variable, we may approximate [(f) by the linear
combination
00) = MAO + >‘2Ac + A3142 + 44A; (624)
4 2 x2

which is exact for the even functions 1, xi x1. and x1 x2 . By symmetry, Q(f) will then also be exact

for all polynomials of degree at most 5.

Now let x = (x1. ' ' ' , xN) and
Mf‘fi'”- max 197 v“v.“f(x)l. (6.2.5)
1km] xeR k I

Theorem 6.2.1

Iff(x] , ~ - ° ,xN) has continuous partial derivatives of the first six orders on R, then

[(7) = [8A0 + (74,. — 4' -A'c' /3)/2N1/15 +E(f) (6.2.6)

C

is a multidimensional quadrature formula with degree of precision 5. The truncation error. E(fl), is

bounded by

N

|E(f)l< 94—-5—-hOthjM6 +35 2h; IrkM 1k

1k 1
j¢k
(6.2.7)

+ 280 Z h] zhzth?“”
j<k<l

102

Proof:

Applying Taylor’s theorem for N-variables to (6.2.4) and equating coefficients of similar terms

we obtain the linear system

 

 

 

 

 

 

F T T 7
1 21V 0 01 811 1
o 2‘"1 2” 0 1.2 %
1 l = 1 (6.2.8)
... N-3 _ zN-l A __.
0 3 2 3 O 3 120
l
2N-2 2N 2N _
0 l. A4 36
L. A _l l. .1
having the unique solution
(8102.83.81) = (8/15, 7/2N15, -1/2N15, -1/2N45). (6.2.9)
The bound on the truncation error is a consequence of Taylor’s theorem and is a straightforward
calculation. Observe that the last term in the error appears only for N > 2. D
Next we transform the variables to obtain a formula for an N-fold integral over an arbitrary
N-rectangle,
N
R = H [aybl J
i=1
N
Letwj=bj—aj, w= H w', m=(ml, - ' ' ,mN),andc=(c1,--',cN)where
i=1
m,- =-;(a,- + bi) and Ci = of or bi. Define
-1 if c- = a-
I I
o. c =
’( ) { +1 otherwise
(6.2.10)

Oik (c) = ol-(c)ok(c).

Corollary 6. 2. 1
N

For any N—dimensional hyperrectangle, R = H [lily bj] .
i=1

103

bN b,
‘J. ...I f(xl,...,xN)dxl...de
UN 01

 

 

 

8w 7W
= —- m + ——
Sf( ) 2(“Slim-f0) (6.2.11)
N
. + E
2N+115 ,ZilDi w 2N2“ Z W ’wk 0km (f)
where
W 6
M? + 35 w w
IE(f)I < 604800 Z w,- :21? W
k
,2 (6.2.12)
+ 280 Z w fwkwmfif
j<k<l
Proof:

Make a linear change of variables in (6.2.7) to obtain

bN b1
(6.2.13)
hN hl w
1 “W w :1
= . . . _— .——.—- ,.._d .. .d .
aNaJ-h J_hlf(2h1xl.'”, ZI’NXN) h x1 xN

Finally, we obtain the composite formulation of (6.2.11). Partition each interval [a-, bi] into n]-

subintervals each of length h]- = wf/nl- and write

N
= “If".
i=1

In order to condense notation we define
11(6) = (a1 + h1(i1-6), ' . . ,aN +hN(iN -0)) (6.2.14)

and

 

Furl}

104
1509.): (01+i1h1'. " ,xj,aj+1+ij+lhj+l, " ',aN+1NhN)
”(xi'xk) = (01 +"1”1""1"/*”j+1 ”Alb/41"” (6.2.15)

xk20k+1 +1.k+lhk+1"” .aN + l.-vh.v)°
Furthemiore
0,176» = v,’1f(v(b,)) - f(v(a,»1
Dig/(v1) = (4' v,'1r(v(a,-. 0,,» - f(v(a,-. bk» (6.2.16)
- f(v(bj9 ak)) + f(V(b]-, bk)” '

We now state our main result.

Corollary 6. 2.2 (MINTOV)

lff(x1. - - ' , xN) has continuous partial derivatives of the first six orders on the N-rectangle

N
R = H [ai'bi] then
i=1

bN b1
la. . “W.

"a.

N
=%Z Z] f(u(1/2)) + ,hlSZ Z f(u(0))

a=1i = a: ,- =0 (6.2.17)

Z Z Z h,- 0,06»

N+l
-2 15 i=1 a= 1i a=0
aatj

 

 

[V N

Z Z ZO' h,h,.D,-,,(f(v))+£(f)

/<k a=1ia=
a4=/,k

 

2A+245

is a composite multidimensional quadrature formula with first- and second-order partial derivative
correction terms having degree of precision five. The primes signify the weight 7 is to be assigned if

the node is common to 7 subregions. The truncation error, EU), is bounded by

«a
“r

t.-

105

 

EN: 222

IE(f)l< 604W800 12:1th6 + 35 114]}!ij + 280 2 h2 'hkhl 2M1“ (6118)
j,k= l j<k<f
j¢k

Proof:

The proof follows by repeated use of (6.2.1 1). 0

Note that the last term in (62.18) appears only for N>2. The 2-dimensional formulation is given
in (6.4.4).

The name MINTOV (Multiple lNTegration, Order V) is given to (6.2.14).

With appropriate interpretation of the corner and centroid nodes and the partial derivative
correction terms, MINTOV is a nonproduct N-dimensional generalization of the composite Simpson’s

formula with end corrections (Lanczos, [28] :

b
I f(x)dx = ——~Zf(a+h(i -1,» + ——Z f(a+ih)

(6.2.19)

 

-—lf '(b)- f(a)] + 116/68)

604 800

As before, the prime on the summation signifies that the weight 7 is to be assigned in case the
node is common to 7 subintervals. Also, :1 = (b - a)/h and § is some point in [(1, b] . The number of
function evaluations is 2n + 3.

It may also be stated that MINTOV is composite Ewing’s formula [16, 49] with partial
derivative correction terms. In the next section, a two-dimensional composite formulation of Ewing’s
formula is given as D0503.

The MINTOV truncation error estimate (6.2.18) is primarily of theoretical importance. It is
useful for comparing MINTOV with other fifth-order multiple integration formulas. In practice it is
usually not possible or at least not feasible to bound the required sixth-order partial derivatives. For

example, it would be tedious to calculate and estimate the requisite sixth-order partials for the function

w sin (x) sin (v) sin (2) e'w

 

f (x. y, 2) =

3 (6.2.20
on. say, the cube H [0,1r/2] where w2 =x2 +y2 + 22.
i=1

106

6.3 COMPARISON OF SEVERAL MULTIDIMENSIONAL QUADRATURE
FORMULAS OF PRECISION FIVE

The attempt to compare various quadrature routines leads to the discovery of the absence of
generally accepted standards of benchmarking techniques. Lyness and Kaganove [30] state that “any
individual who constructs a routine can find some problem for which it is more efficient than an
existing available routine, and with this evidence, arrange for its inclusion in the local subroutine
library. Existing routines are not removed because there are other problems for which they are more
efficient than the new routine.”

They add that numerical experiments must play a basic role in the comparison testing and then
distinguish between “battery experiments” and the “performance profile evaluation technique.”

The battery experiment method requires the use of many different integrand functions, limits
of integration, tolerances, and different quadrature routines. Two notable investigations were con-
ducted by Caselleto, Pickett and Rice [12] at Purdue University, and by Kahaner [25] at Los Alamos
Scientific laboratory.

Kahaner tested 11 quadrature routines on 21 integrands using 3 error tolerances and eventually
selected 3 routines based on average reliability and average speed.

Lyness and Kaganove [30] object to the technique used in these investigations on the basis of the
difficulty of interpreting the results and of extracting definite conclusions, and because integrand
functions which are “close” to one another were not used. Consequently, they recommend the use
of the more popular performance profile evaluation technique which is described in detail in [31].

We will apply neither the battery experiment nor the performance profile technique to assess the
merits of MINTOV. Instead, we will make general comparisons and conduct several numerical
experiments using some definite integrals of our selection and several found in the literature.

Now we observe that the number of function evaluations, nfe, required for the d-dimensional
MINTOV (6.2.14) is given by

d

d d a a
nfe = [I n,- + “(ni+l)+ZZn(ni+l)+4Z I] (ni+1). (6.3.1)
i=1

i=1 i=1i=1 j<k i=1
i¢j i¢j,k

107

Furthermore, in the case "i = (b,- - ai)/hi = n for all i,

nfe = nd + (n+ 1)d + 2d(n+ 1)d‘1 + 2d(d—1)(n+ 1)d'2
6.3.2
d (H d ( )
: ——‘ .. ' 2 i
2n + L [2(d 1) + ”(1.)".
i=0
Hereafter we refer to n as the “number of subdivisions” of R. Note that R is partitioned into nd sub-

regions. In (6.3.1) a partial derivative evaluation is counted the same as a function evaluation. For
some integrands it may be necessary to weight the partial derivative evaluations. However, in many
cases, because of persistence of form, the present enumeration technique will suffice for our purpose,
namely, to size MINTOV and compare it with several well-known formulas.

Indeed, for functions such as ln(xyz) the first order partials require only about 40% of the time to
evaluate the given function whereas functions similar to cos(x)cos(y)cos(z), exp (-xyz),

(l + x + y + z)'4, and even (6.2.18) require not more than 5% additional time to evaluate the first-
and second-order partials than the original functions.

We will compare MINTOV with the fifth-order composite multidimensional quadrature formulas
Lyness, Gauss, and Boole, which refer to the composite formulations of C" :5-5 as listed in Stroud [49]
and which is due to Mustard, Lyness and Blatt [37], the composite product Gauss, and the fifth-order
composite product Newton-Cotes formulas, respectively.

MINTOV and Lyness are nonproduct formulas and as such might be expected to be more
efficient than Gauss and Boole, as indeed they are.

The number of function evaluations for each rule is given in Table 6.3.1.

It can be shown that for n > d, MINTOV nfe < C" : 5-5 nfe; equality holds only for n = d = 2.
Also, for any n > 1 and any d, MINTOV nfe < Gauss rife. Moreover, for all n,d, MINTOV nfe <
Boole nfe. For example, for 8 partitions in 4 Space, MINTOV, Lyness, Gauss, and Boole require
18 433, 43 425, 331 776, and l 185 921 function evaluations, respectively.

Table 6.3.2 lists the number of function evaluations required by MINTOV for various subdivisions

n and dimensions d.

Table 6.3.1

108

Number of Function Evaluations for Several
Fifth-Order Formulas

 

 

 

Formula nfe
MINTOV nd + (n + l)d'2[(n + 1)2 + 2d(n +d)]
Lyness (n + DJ + (2d + l)nd
Gauss (3n)d
Boole (4n + 1)d

 

 

 

Table 6.3.2 MINTOV: Number of Function Evaluations Required for n Subdivisions
in d-Dimensions

 

 

 

 

 

 

 

 

 

 

 

 

 

 

"(1 1 2 3 4 5 6 7 8 9 10
1 5 17 57 177 513 I 409 3 713 9 473 23 553 57 345
2 7 29 125 529 2165 8 569 32 933 123 457 453 221 l 634 713
4 ll 65 399 2 481 15 399 94 721 575 759 3 456161 20 496 519
8 19 185 1835 18433 186587 1985 945 19 280 411
16 35 617 10 947 195 297 3 500 163
32 67 2 249 75 635 2 548 129
64 131 8 585 562 899
128 259 33 545 4 345 235
256 515 132 617
512 1027 527 369

 

MINTOV may be sized by computing the maximum number.of subdivisions n such that the

number of function evaluations is less than some preassigned limit, say 106. This is done in Table

6.3.3 and the results indicate that compared to the techniques listed, MINTOV is substantially

superior for dimensions l-S.

 

109

Table 6.3.3 Maximum Number of Subdivisions n such that nfe < 106

 

 

Rule d 1 2 3 4 5 6 7 8 9 1o
MINTOV 499 998 705 77 24 ~12 6 4 3 2 1
Lyness — 408 49 17 9 6 4 3 3 2
Gauss 333 333 333 33 10 s 3 2 1 1 1
Boole 249 999 249 24 7 3 2 1 1 — —

 

 

 

 

Finally, considering the speed of today’s fourth-generation computers (cycle times are measured
in nanoseconds), it is not unreasonable to expect a composite multidimensional quadrature formula to
perform at least n = 4 subdivisions using. say, 106 function calls. With these admittedly rough but
realistic guidelines, it is possible to estimate the maximum usable dimension for the formulas under
consideration. The results in Table 6.3 .4 indicate that the useful dimensional range for MINTOV is
1-7. Later we will show that MINTOV is particularly efficient and accurate in dimensions 2 and 3.

In particular, from Table 6.3.4 we see that if there is a requirement that a multiple quadrature
formula employ at least n = 4 partitions (i.e.,4d subregions), and not more than 106 function evaluations,

then the maximum usable dimensions are 5 for Gauss and 7 for MINTOV.

Table 6.3 .4 Maximum Usable Dimension 0' Assuming n 9 4
and nfe < 10“

 

 

"m" nfe
Rule 102 103 104 105 106
MINTOV 2 3 4 6 7
Lyness — 3 4 6 7
Gauss l 2 3 4 5
Boole 1 2 3 4 4

 

 

 

 

110

6.4 CONSTRUCTION OF 47 NEW CUBATURE FORMULAS WITH PARTIAL

DERIVATIVE CORRECTION TERMS AND ERROR ESTIMATES

In this section, for concreteness. the discussion will be limited to 2-dimensional quadrature or

“cubature” formulas. We will simplify the notation whenever possible.

We wish to approximate the integral

d b
1(f) =ff f(x.,V)dxd.V

(6.4.1)

over the rectangle R = [a, b] X [c, d] by a cubature formula Q(f) which contains partial derivative

correction terms.

Partition R into nm subrectangles each of size hk where h = (b -a)/n and k = (d~c)/m. Define

the following cubature elements:

likZZf(a+h(i-‘/z), c+k(j-'/z))

F0 =
j=li=l
FV = hk Z): f(a+1'h,c+jk)
i=0 i=0
FM = lsz [f(a+h(i-'/2).c) + f(a+h(i-‘/2).d)]
i=1
+ hk ma. c+ko°- a» + f(b.c+k(i-%))l
j=1
m n-l
+ thZ f(a+ih,c+k(i-‘/2))
j=l i=1
"'1 ——.
+ th Lf(a+h(i-‘/z),c+jk)
j=l i=1
FV] = hzk Z '[ fx(b. c+ fk)- fx(a.c+ 12)]
i=0

+ hk2 Z ' [fy(a +171, d) -fy(a +171, c)]
i=0

(6.4.2)

(Equation (6 .4.2) continues)

where

111

FM] = 1.sz mo, c+k(i - 14)) -fx(a. c+ko°- v2»)
j=1

+ hkz 2.161" + h(i - ma) -f,(a +h(i - ‘ri),c)]

i=1

FVll = h3k2[ fxyo, c)- fx},(b, c)+ fxy(b. .1)- fx},(a, d)]

g=vvom
g=vyww
fxy = ”1” } f(x,,v).

M11]2 = maxR l'l’z’j'lf(x.y)|
(x.)')e

Fij = (b-a)(d-c) [hikiM'iff hikifli{iz] = F/‘i.

mon to a subrectangles.

As in the case of MINTOV, the method of undetermined coefficients is used to construct a

(6.4.2)

(6.4.3)

(6.4.4)

The primes on the summations in (6.4.2) signify that weight a is to be assigned if the node is com-

variety of new cubature formulas. The results are compiled in Table 6.4.1. Of the 88 combinations

considered, 52 had nonvanishing determinants. Some of the unsuccessful attempts are indicated

by zeros.

The significance of the entries in Table 6.4.1 may be understood by observing that formula

DCSC5 is the 2-dimensional formulation of MINTOV (6.2.17):

fd bf(x )dxd ——8-F0+—7—FV i—FVi —1-FV11+F(f)
c]; 'y y 15 60 120 720 ‘ '

|E(f)| < (F60 + 35 F42)/604 800.

The number of function evaluations is

nfe = 2(nm) + 3(n+m) + 9.

(6.4.5)

(6.4.6)

(6.4.7)

Table 6.4.1 Cubature Formulas

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Elements [Error Bound nfc
No. Name
170 W FM 131 [MI W" 120 140 1522 160 142 nm n+m c
Midpoint 1
1 [£0101 1 271 1 0 0
— 131140 0 0 I 0 4
1 -7 1
2 EM143 1 24 m m I 2 0
l -7 -5
.’ _ __ __ 7
3 1.T183 1 48 5760 576 1 - 4
4 15111335 1 —'— J— :7— 1 2 4
’ 24 576 5760
.. . _L :5_ __'7
5 L(1(,3s 1 4a 576 5760 1 2 8
. ~ _L _5_ _‘7
6 1'5”” i 288 144 5760 l 4 4
— EHIGO 0 0 0 0 1 4 8
Trapezoidal
1 -1
7 T0401 4 1—2 1 1 1
— T1440 0 0 0 1 1 5
1 -1 1 -1
8 TM443 Z 12 7—2—0 3 1 3 l
1 -1 1 1
10 T\'483§ i "—1 i —'— 1 3 5
' " 4 12 72 720
11 T("4C3S l i —'— —'— 1 3 9
' 4 24 144 720
.. . . 1 —1 —1 1
12 T6403. 3- E 36 7—2—0- 1 5 5
_ T114011 0 0 0 0 1 5 9
Squire
1 -1
13 110401 Z 4—8 2 1 0
-— 311-440 0 0 ‘0 2 1 4
1 -1 1 1
l4 ““443 I W 11520 5‘76 2 3 0
1 -1 1 1
_ —- — 7
‘5 ”T433 4 96 11520 144 “ 3 4
1 —1 1 1
.' ‘ _ _ _ ‘)
l6 ”"4835 4 4a 576 11520 ‘ 3 4
, . l -1 1 1
‘7 M(.4(3S 4 '9—6 F 11530 2 3 3
18 MS4(‘3‘ 1— _l_ :1 I , 5 4
'3 4 2311 36 11520 '
— 11111400 0 0 0 0 2 5 8
[using
19 110503 3 i —"— LL 2 1 1
* 3 I2 2880 288
20 015435 3— i " Si 2 1 5
' .. 3 12 288 2880
8 1 1 -1
21 DMS43A 3 3—6- -3_6 m 2 3 1
22 D~15433 3 '7— ;1 ll— 2 3 1
' ‘ 15 60 60 180
4 5 -1 1
S —. ... __
23 DT.83A 9 36 72 4320 2 3 5
8 7 -1 -1
4 _ _
2 [”5833 15 60 120 720 2 3 5

 

 

 

 

 

 

 

 

 

 

 

113

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Table 6.4.1 (cont‘d.)
Elements l-lrror Bound nfe
No Name
to W I-M 17V] 1w 1=v11 140 1:22 1‘60 142 nm n+m c
____ _ _,._
8 7 -1 -1 l 7
25 0x585 ~1—5 m 66 “1+1—0 604800 69120 3 3 3
MINTOV
, 8 7 -1 -1 1 1
36 ”3‘33 F E 120 720 604800 17280 2 3 9
8 7 -1 1 1 1
_ __ __ _ __ __ 2 5 5
27 055“ 15 60 90 180 604800 23040
054 8 7 _7 1 1 1
28 15115055 73 60‘ 3710 R“ 240 6’——"04800 3 3 9
Tyler
29 x0503 i —' -1 l— 3 l 0
-~ 3 6 2880 576
30 x1-‘543s _l— i —' " 3 1 4
3 6 576 2880
1 7 -1 1
.4 3 ._ — — ——
31 x154.r 15 30 60 576 3 3 0
4 5 1 -l7
32 XTS83A 3 E _2_8 34560 3 3 4
1 7 -1 I7
33 XT583B '1—5' 3—0 '20 2880 3 3 4
1 7 -I 1 1 '7
34 xx585 E :5 60 57—6 604800138240 3 3 4
‘ ‘ 1 ‘7 -l 17 1 ‘13
35 ’“5‘5 E 30— 12?) 28% 60480013824013 3 8
. 1 7 1 -17 1 -1
3" X55“ E 37) 288 7271 __604800 30720 3 5 4
. 1 7 7 -13 —1 1
37 x11505s T5 301720 360 32—0 604800 3 S 8
Miller 1
. -1 1 , -1 1
38 00303 T2 3] 2880 1‘47 3 3 l
-1 1 i I -|
-43. — —- ——
39 0'3 3 12 3' 14 7880 13 3 3
40 OMB43A i 3’ i I 3 4 '
36 9] 36 4320
-1 4 1 -1 1
M84313 — —- _—
3' O 60 15] 60 360 3 4 l
, _ -1 4 1-1 I
2 883 — — ——
4 or r 60 15 120 144 3 4 3
-1 4 -1 1 I I
43 0“” a r5 5 no 604800 ram ‘ 4 5
, -1 4 -1 1 I -l
4" 003” 6—0 B 120 T4? ,604800 8640 3 3 9
, -1 4 I -5 l '1
45 058(5 60 T5 17?) 180 604800 23040 3 6 3
. —1 4 5 -2 —1 l 1
46 0119055 E 15 327) F 240 1604800 3 6 9
_l 1
Simpson
47 309035 4— '— i —" 4 7 1
9 36 9 2880 '
— 51940 0 0 0 0 I 4 2 5
8 1 8 -1 . l -'
48 531945 4—5— 33 4-5- 35 1604800 69120 ‘3 3 l
4 17 2 -1 ' l l
4 _. _. _ ._ _—
9 ”985 9 180 45 120 1604800 34560 4 4 3
2 7 7 -1 —1 1
50 SX985$ -9- W 3; K5 1735 604800 3 3 5
16 13 4 -1 1 1
5 4 4S _ ~— _ __ ______ -—-——'—
1 st9(-s 45 180 45 120 720 604300 3 3 9
4 1 2 -1 -1 l
52 559 55 ~ -— —— ——- —
C 15 20 15 360 90 604800 4 6 3
4 511900 0 0 0 o 0 0 4 6 9

 

 

 

 

 

 

 

 

 

 

114

The first two symbols of the names assigned to the fonnulas listed in Table 6.4.] were chosen
somewhat arbitrarily, whereas, a digit was selected for the third symbol to represent the number of
function evaluations required for the holistic or basic rule. that is for n = m = 1.

The fourth symbol represents the number of partial derivative evaluations required for the holistic
rule. Here C= 12 and G = 16. The fifth digit is the order or degree of precision.

The presence of a sixth symbol signifies that the method of undetermined coefficients led to more
equations than unknowns. The symbols A and B are used to indicate two successful combinations: whereas
an S or T indicates that only one combination was successful.

For completeness, formulas l, 7. and 47, which are the composite formulations of the midpoint,
trapezoidal, and Simpson’s rules, respectively, have been included. Good and Gaskins [20] recently
investigated the midpoint or centroid mle. The holistic formulations of formulas 13, 19, 29, and 38
were investigated by Squire [48] , Ewing [16], Tyler [51], and Miller [35] respectively. Thus except
for formulas 1, 7, 13, 19, 29, 38 and 47 the remaining 45 cubature formulas are new. Formula
X0503 was apparently discovered independently by Bickley [7] and Tyler [51]. We will follow Stroud
[49] and call it Tyler’s rule.

It is interesting to compare DF543S with Simpson’s rule, 809033, because both are third-order
formulas and both have the same error bounds; however, the former requires half as many function
evaluations as Simpson’s rule. Moreover, DF543S requires only 4 second-order mixed partial
derivatives evaluated at the vertices of the rectangle R. There may be applications where DF543S
should be considered as a viable alternative to Simpson’s Rule.

DF543S is a combined trapezoidal-midpoint or Ewing rule [16] with boundary correction terms.
This was the first cubature formula discovered which requires mixed second-order partial derivatives
evaluated only at the comers of the rectangular domain of integration R. Observe that the two formulas,
XFS43S and OF893S, which share this property with DF543S also have the same error bound as

Simpson’s rule. They are also more efficient than Simpson’s rule.

115

For smooth functions, DHSGSS will likely provide the greatest accuracy with the most economy.
For brevity, we shall refer to DHSGSS simply as “CSA” (Corrected Sth-order Approximation).

Except for DF543S, XF543S, and OF843S, the derivative corrected formulas listed in Table 6.4.1
exhibit the same property as the formulas of Chapter 4: namely, the partial derivative correction terms
are evaluated on the boundary as well as the vertices of R. Nevertheless, the cubature rules with
boundary correction terms are more efficient than conventional formulas.

Since the nodes do not form a lattice, it follows that most of the cubature formulas in Table
6.4.] are not product formulas.

Since the DC Simpson quadrature rule, (6.2.19), is constructed using the method of undetermined
coefficients, it was conjectured that

SH9G0 = MFO + AZFV + A3FM + A4FV1 + ASFMI + A6FV11

would be the product formulation of (6.2.19). However, the coefficient matrix of the resulting

 

 

 

 

 

 

system
"1 4 4 0 0 0 7 PM F 1 7
0 4/2! 2/2! 4 2 0 82 1/31
0 4/4! 2/4! 4/3! 2/3! 0 713 1/51
= (6.4.8)
0 4/2!2! 0 8/2! 0 4 8,, 1/313!
0 4/6! 2/6! 4/5! 2/5! 0 7.5 1/71
_0 4/4!2! 0 1/2 0 4/3!J J6- L0513!]

was singular and the attempt was unsuccessful.

A second and different derivation of the DC Simpson quadrature rule (6.2.19), employs the Euler-
Maclaurin Summation formula. This was done in (3.1.9). (Recall that prior to Chapter 6, h denoted
the distance between nodes, whereas here h = (b - a)/n and k = (d - c)/m).

Similarly, the double Euler-Maclaurin Summation formula may be used to construct the DC

Simpson cubature rule (5.2.7), which is a product formulation of (3.1.9).

116

The formulas in Table 6.4.1 are first-, third-, and fifth-order rules. It is somewhat disappointing

that SH9G¢ did not produce a seventh-order formula. However, by selecting different nodes, several

seventh-order derivative-corrected rules can be constructed. To this end, define the elements

2m n
Z Zf(a - h/2 +171. c + k/4 +jk/2)

FA =
i=1 i=1
m 2n
+ ZZf(a+h/4 + ih/2,c—k/2 + jk)
i=1 i=1
2m—I 2n-l
FB = Z Z f(a+h/4+ih/2,c+k/4 + jk/2)
i=0 i=0

with nodes arranged as shown in Figure 6.4.1.

 

 

 

 

 

 

 

 

O O O
O O
o o 0
FA FB
Figure 6.4.1 Node Arrangements
Thus we have the seventh-order formulas
fdfbfl )dd 44F0+13 FV+10FM
x. x = -— -— 1 —
C a y y 105 420 189
256 l l
+ —FA - -—FV1 + —FV11+E(f)

945 504 3360

|E(f)| < (F80 + 28 F62 + 364F44)/l 625 702 400

and

(6.4.9)

(6.4.10)

(6.4.11)

117

fdfbflx )dxdv = 1m + fl—FV + lé—FM
c a 'y ' 63 378) 315
(6412)
+13—8-FB - —1—FV1 + —l—FV11+E(f).
945 504 1440
16ml < (F80 + 112F62 + 854F44)/1625 702 400. (6.4.13)

Here (6.4.10) is a derivative corrected Tyler [51] cubature rule and (6.4.12) is a derivative corrected
Albrecht, Collatz [3] and Meister [34] rule. In both cases
nfe = 8(nm) + 4(n + m) + 9 (6.4.14)
which is substantially better than the 167m fe required by the seventh-order composite product
Gauss formula.
Finally. we note that as in the case of MINTOV. the formulas of this section can be generalized
to multidimensional quadrature rules. Indeed, with appropriate generalizations of F0, FV, etc., for

N23 it can be shown that

b . hi
IV .
f j(-x19333.-x‘~')datl ...d'Y/V
01‘! ° 01
8(14-5N) . 33 4
=—_——- +——."— I _ .11
135 27113533327“
(6.4.15)
1 PM 3 FM] ] FV11+F(f)
T70 '13s 21’540 ’ '

N
W 6,
[E(_f)|<b———04 3002,33 11f (6.4.16)
I:

7. NUMERICAL RESULTS

7.1 DOUBLE INTEGRALS AND THE DC—MQUAD ALGORITHM

The 2-dimensional formulation of MINTOV is an accurate and efficient cubature rule. As pre-
viously noted. it requires 201177) + 301 + m) + 9 function evaluations.

Let the rectangle R = [0, b] X [c, d] be partitioned into nm subrectangles each of size
hk = (b — a)(d — c)/nm and let Q(f) denote the approximation provided by a cubature rule.

For ease of reference the two-dimensional formulation of MINTOV is given below.

d b M
Li f(x.y)dxd.v=—:Zn lf(a+h(i-%)C+k(i- I5))

j= i=
7hk
ZEZZWH” ”1") (7.1.1)
0i= -0 i= -0 '
h2k"'
—1—ZOZ [f (b c+rk)- f (a c+/k)]
2
-’;——:O :5) [f, (a+ih, d)— f,(a+ih, c)]
h 3k2
_ 7__20 [fxy(a, c)— fxy(b,c) + fxy(b,d) -fxy(a,d)] + Etf).
b- d-
IEW < $655)- [016/116+ k6M6) + 350141.3er42 +11 1411424)] (7.1.2)

Mustard, Lyness, and Blatt [37] prOposed a 9-point, degree 5 cubature formula which when
compounded nm times requires 6nm + (n + m) + 1 function evaluations, fe. As before, we refer to
this as Lyness. This is the most efficient, non-derivative-corrected composite fifth-order cubature

formula given in the literature.

The Radon [42], Albrecht, Collatz [3] 7-point, degree 5 composite formula requires 7nm fe.

For brevity we will call it Radon’s formula.

118

119
The 9-point, degree 5 composite product Gauss cubature fomiula requires 911m fe, and the

25-p0int, degree 5 composite product Boole’s rule when compounded mn times requires

(411+ 1X4m+ 1) fe.

Tanimoto's [50] corrected Simpson rule requires 4nm + 6(n + m) + 9 fe. Tanimoto’s rule is

the only fifth-order derivative corrected cubature formula which we found in the literature.

Table 7.1.1 shows the number of function evaluations these formulas require for various

subdivisions. Clearly, in this respect. MINTOV and CSA are superior to Tanimoto, Lyness, Radon,

Gauss, Boole, and even Simpson‘s rule.

MINTOV and CSA (formulas 26 and 28, respectively, in Table 6.4.1) have been tested on a

variety of integrands with many different grid sizes and have produced excellent results with respect

to accuracy, computational efficiency, and economy.

The Lyness rule also performed well on the same examples; however. it required up to 3

times the number of function evaluations and was therefore much more expensive than MINTOV or

CSA. Tanimoto’s rule required nearly twice as many function evaluations as either MINTOV or CSA.

We now propose a new technique for cubature with error estimates called the DC-MQUAD

(Derivative Corrected—Multiple Quadrature) Algorithm.

 

 

 

 

 

Table 7.1.1 Nfe for Various Cubature Formulas Compounded nm Times. These are Fifth-Order
Formulas Except for Simpson’s Rule Which is Third-Order.
MINTOV CSA Simpson Tanimoto Lyness Radon Gauss Boole
n=m 2n2+6n+9 2n2+10n+9 4n2+4n+1 4112+1211+9 6112+ 2n+l 7n2 91:2 (4n+1)2
l 17 21 9 25 9 7 9 25
2 29 37 25 49 29 28 36 81
4 65 81 81 121 105 112 144 289
8 185 217 289 361 401 448 S76 1 089
16 617 681 1 089 1 225 1 569 1 792 2304 4 225
32 2 249 2 377 4 225 4 489 6 209 7168 9 216 16 641
64 8 585 8 841 16 641 17 161 24 705 28 672 36 864 66 049
100 20 609 21 009 40 401 41 209 60 50! 70 000 90 000 160 801

 

 

 

 

 

 

 

 

 

 

 

 

 

 

The DC-MQUAD Algwithm

]. For a given step size, compute the 4 cubature elements F0. FV, FVI, and FV11. (For
2-dimensions this requires 2nm + 3(n + m) + 9 function evaluations.)

2. Apply the appropriate weights given in Table 6.4.1 to these cubature elements, and compute
the first-order approximations, trapezoidal (T0401) and midpoint (E0101); the third-order
approximations, Ewing (D0503), DF543S, and DT583A; and the fifth-order approximation MINTOV
(DCSC5) (see Table 6.4.1). In the case of apparent convergence, accept MINTOV as the approximation
to [U ).

3. Estimate the actual and relative errors for the first five cubature rules by computing
MINTOV-Q(f) and [MINTOV-Q(f)] /MINTOV, respectively. Use MINTOV (hz, k2)—

MINTOV (hl , 1‘1) and [MINTOV (112. k2)-MINTOV (ht , k1)] /MINTOV (I12, k2) to estimate the
actual and relative errors for MINTOV. Here hz < hl and k2 < k1. Stop when the error estimate
meets the preselected requirement.

4. Decrease the grid size and repeat steps 1-3.

For a sufficiently smooth function. say one having continuous partials of the first six orders, this
technique provides a close error estimate which is not only practical but is also computable.

The convergence of the DC-MQUAD algorithm is guaranteed by the following theorem.

Theorem 7.1. I

If f(x, y) has continuous sixth-order partial derivatives over R = [0, b] X [c, d], then

= 0. (7.1.3)

lim d b
th + k~’—o0 MINTOV -f I f (£1061me
C 0

 

Hmf'
Since the sixth-order partial derivatives are continuous over a compact region R, they are

bounded on R. Then clearly. the truncation error is bounded by

(b - a)(d - c)

E <
I (f) I 604 800

[(126M? + 165113 )+ 350141121111; + 1121611111331” (714)

and goes to zero as V h2 + k2 goes to zero. D

121

Similar theorems can be proved for each of the formulas listed in Table 6.4.1.
Since the cubature formulas in Table 6.4.1 generalize to n-dimensions, it follows that the

DC-MQUAD algorithm can be generalized for the approximation of multiple integrals of the form

n h!
[(f) =1 J‘ f(x1,°'°,x,,)dx1~“dxn. (7.1.5)
at 0,

Finally we note that if the six cubature elements F0, FV, FM, FVl , FMl , and FVll are com-
puted, by using only 4nm + 6(n +m) + 9 function evaluations, all 52 approximations to the integral
[(f) furnished by the cubature fomiulas of Table 6.4.1 may be obtained. In this connection, see
Section 7.1.5 for numerical results.

Now we give the results of some numerical experiments.

7.1.1 THE APPLICATION OF MINTOV

Consider the approximation of

f J01d_x__dyz=L"/2J'O/2 dxdy _Il tan'1(x)dx
1+); 2(1+COS(JC)C05(y)) 0 X (7.11.1)

= 0.915 965 594177 30.

 

Takingf(x,y) = (1 +x2y2)‘l and h = k = 1/2 in (7.1.1) we obtain
~ 8 1 256 256 256

_ —— — —— + —— +
I“) 15 4l257 2(265) 337]

+ i —[1 + 2(1)+1/2 + 212 1111+ 212x4/51+4<16/ 17)]

 

60 4
l 1 . . .2

 

 

169281536 + 17213 + 57
344 270 775 40 800 24 000

0.49] 710445 + 0.421 887 255 + 0.002 375 000

0.915 972 700 0.
The error is 0.000 007 106 8. Here fx = £03211 + x2)32)’2,f‘, = -2x2)'(1 +9112)”, and

- _ , 2 ,2 ' 7 a -7 . . . . .
fxy ‘ 4200‘ .1 '1)“ + x3") “. The first-order partial derivatives ft and fr vanish on the axes while

the second-order mixed partial derivative fxv vanishes at the corners as well as on the axes.

 

 

 

 

 

 

 

 

 

 

 

 

ILL
/

 

 

 

 

 

 

 

 

 

LLL

 

 

 

K

122

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

N

/

J

2
ll
0,
2)-l on 1
2y
+ x
= (1
h of 2
P
Gra
l
1.
1.
7.
re
Figu

 

 

123

 

 

 

 

 

 

fxy=0 fy=0 f)’=-8/25 f,=-l/2 fxy=0
fx=0 1 4/5 3/2 fx=-1/2
256 256
' 265 ' 337
1 16/17 4/5
fx= 0 f, =-8/25
256 256
.— .-—.-
257 265
fx=0 1 l 1 fx=0
fx :0 fyzo fy=0 fy=0 fxy=0

Figure 7.1.1.2 f(x,y) = (1+x2y2)"

Figure 7.1.1.2 shows the function and partial derivative values used in the calculation (7.1.1.2).
Note that the partial derivative correction terms are evaluated only on the boundary and not at
interior points.

In Table 7.1.2 MINTOV is compared with several other cubature formulas.

The examples were run on the U.S. Naval Air Test Center’s Real Time Telemetry Processing
System Xerox Sigma 9 computer using the CPR version C008 operating system and Fortran IV in
double precision with 15 significant decimal digits.

Figure 7.1.1.3 shows a plot of the grid size 11 = k vs. the logarithm of the absolute error,
log |I(f) - Q(f)l. The theoretical results indicate and the numerical results confirm that in terms
of accuracy and economy, the derivative corrected cubature formulas provide approximations

superior to comparable conventional rules.

124

Table 7.1.2 10‘ f0'(l +x3yZ)-1dxdy = 0.915 965 594 177 30

 

 

 

 

 

 

 

h=k=l/5(n=m=5) h=k=1/10(n=m=10)

Rule nfe 382:]; Error nfe 121;]; Error Ordet
Trapezoidal 3 6 .006 1 .90-3 * 121 .023 4 .76-4 1
Midpoint 25 .005 —9.52—4 100 .020 -2.38-4 1
EM143 45 .015 -l.11-6 140 .039 —6.97-8 3
Ewing 61 .012 -3 .44-7 221 .043 -2 .04-8 3
DF543S 65 .014 -3.44—7 225 .045 -2.04-—8 3
Simpson 121 .023 -3. 16-7 441 .086 -1 .99-8 3
MINTOV 89 .026 -2.20—8 269 .065 -3.39-10 5
CSA 109 .035 431-10 309 .084 861—12 5
SC9CSS 149 .037 5 .46-10 489 .108 904-12 5
Tanimoto 169 .046 591-10 529 .128 922-12 5
Lyness 161 .028 566-9 621 .108 870-1 1 5
Radon 175 .030 -1 .84-9 700 .1 16 —2.81—ll 5
Gauss 225 .037 1.78-10 900 .145 283-12 5
Boole 441 .078 -1 .85-10 1681 .297 -2.77-12 5

 

 

 

 

 

 

 

 

 

 

 

 

 

*By 1.90-3 we mean 1.90 X 10'3
fDegree of precision

125

 

 

dxdy

1 ( )2 = 0.91596559417730
+ xy

ff

0

 

DCSA

O MINTOV

* FUNCTION EVALUATIONS

l l L l I I I l I I l I l

 

101
10° ~
.8
11 "
.2
I.“
§ _.
U)
Q
C
(9
10°1 _
to"
10‘“
Figure 7.1.1.3

10"2 1o“° to“ 10“5 to" 10" 10°

LOG (ABSOLUTE ERROR)

Error Curves in Approximating f1 f l (1 + x2y2)'1dxdy
o o

 

7.1.2 MINTOV vs. JPL’s MQUAD

Bunton, Diethelm, and Haigler [8] (hereafter referred to as Bunton) of Jet Propulsion

laboratory proposed the following example.
2.1 2.1
f f (xy)"dxdy = 1n2(2.1) -- 0.550471 023 504079. (7.1.2.1)
l l

MQUADl is a Modified Romberg Multiple quadrature routine in the JPL program library. Since the
single precision version computes only the first three columns of the Romberg Table, it has degree of
precision 5.

Bunton states that MQUAD was designed to “. . . satisfy all the needs of a large scientific and
engineering computer community. The routine is to be the ‘library standard,’ and its use should be
greater than the total of all other special purpose quadrature routines. The standard should be one

9

against which others could be measured. . . .’

Since both MQUAD and MINTOV are fifth-order methods, it seemed reasonable to compare them.
The relative error tolerance requested for MQUAD was 10“. The grid size for MINTOV was decreased

until the relative error as estimated by

|(M1NTOV (h/2, k/2) — MINTOV (h, k))/MINTOV (Ir/2, k/2)l
was smaller than 10'4. Under these conditions MQUAD required 441 function evaluations while
MINTOV used only 65 function evaluations.

The results of MINTOV and MQUAD for a variety of relative error tolerances are presented in
Table 7.1.2.1. Here MINTOV required from 1/2 to 1/26 the number of function evaluations required
by the IPL routine MQUAD. Of greater significance is that MINTOV required 2% to 28% the actual
computer time required by MQUAD. On the average, MINTOV is better than MQUAD by an order
of magnitude.

Figures 7.1.2.1 and 7.1.2.2 show the number of function evaluations vs. absolute error for 10

different cubature formulas. For this example, SC9C5S, Tanimoto, Lyness, and Radon exhibit

 

1I am grateful to Mr. Wiley R. Bunton of JPL for sending a I-‘ortran listing of MQUAD.

127

Table 7.1.2.1 MINTOV vs MQUAD for the Integral 1111112" (x_v)'1dxdy =
0.550 471 023 504 079. Relative Error Requested = 10'“, a = 1(1)10.

 

 

 

 

 

 

Relative MINTOV MQU AD

Error . . . .

Reeeeeee .r. 121‘? “£113? .r. (Tsar R812?
10'1 17 .002 141-3 441 .157 -1 .38-7
10'2 17 .002 141—3 441 .157 -1.38-7
10'3 29 .003 326-5 441 .157 -1 .38-7
10‘4 65 .010 591—7 441 .157 -l.38-7
10'5 185 .030 967-9 625 .217 -3.51-8
10‘6 617 .100 153—10 1 089 .370 -4.93-9
10'7 617 .100 153-10 2 025 .676 —3.05-10
10'8 2249 .365 243—12 4 761 1.567 -1 .19-11
10'9 8585 1.395 9.4144 15 129 4.958 —3.92—13
10'10 8585 1.395 9.41-14 42 849 13.965 -3.27-l4

 

 

 

 

 

 

 

 

 

 

 

*Relative error = [1(f)-Q(f)I/l(f)

approximately the same efficiency. That is, for a given number of function evaluations, they produce
approximately the same absolute error.

The graphs indicate that the formulas XSSCS, CSA, and MINTOV are the most efficient cubature
rules for this double integral, and MQUAD is the least efficient.

Here the Lyness truncation error is estimated by

(b - a)(d — c)

IE I <
(f) 12 096 000

[(h6M? + keMg) + 175(h4k2M‘g + 122/(WEN (7.1.2.2)

7.1.3 ERROR ESTIMATES
Stroud [49] uses the Radon-AIbrecht-Collatz fomiula which he designates C225-1 (and which

we call Radon) to approximate

i i
4
I J, V3+x+y dxdy = 7511 -18\/3 + 25V5) (7.1.3.1)
-1 -1
z 6.859942 64033465

and then gives several error bounds.

128

 

100m T

 

 

 

I‘fl “
.. \ 2.1 2.1
\ f f 921v. . ...2 2.1
r "V
\ 1 1
\ a
_ \ ‘
" \\ \ o MQUAD
I. \ \ A BOOLE
" kt “ \9 o GAUSS
' \ \
MINTOV
\ \ °
\ \ 0 CSA
\
\ .. v xsscs
1000 _ . \
I \ \
i- n \ \
I— § 0 \
'3 \
o " \
r: \ GK
3 - \
\
g \
m - \
5 \ \
p.
U 1‘
- ‘ \
8 \ \
9 “ \
1oo - \
_ . \ \9
.. t t
7 \
— \‘ 0 L‘
h \‘
\-
a
)—
\
v E
‘0 1 1 l l l L i i i i 1 i i
10'" 10"2 10‘10 to" 10" 10" 10"

LOG (ABSOLUTE ERROR)

Figure 7.1.2.1 Performance Comparison in Approximating 1.12.1 f12'1(xy)'ldxdv

 

 

 

 

 

 

10 000 _
,_ A
_ “ 2.1 2.1
f f ——Md = In2 2.1
_. xv
-. 1 1
Cl SCSCSS
o TANIMOTO
O LYNESS
a RADON
1 000 __ .\
*- \\
P'— \‘\ ‘
r- \\
g “\
9 h \\
i- \
< ’— \
3 .\
a‘
h— \‘
a \e
Z \\
9 r—
p.
o
z
3
IL
& I—-
O
o'
z
100 _
10 I 1 I l I I I I
10’“ 10"2 10"° 10" 10" 10" 10" 10°

LOG (ABSOLUTE ERROR)

Figure 7.1.2.2 Performance Comparison in Approximating [12" f12.1 (xy)'ldxd_v

130

Using (7.1.2), the MINT 0V truncation error estimate for this example is

(7.13.2)

[(h" +k6) + 3501“]:2 + h2k4)]

l

10 240

 

IE(f)|<

=k

orinthecaseh

(7.133)

9116
1280

IE(f)|<

The results of applying various cubature formulas and error estimates to the integral (7.1.3.1) are

tabulated in Table 7.1.3.1.

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Figure 7.1.3.1 Graph of z

3+x+y on [-l,l]2

131

Table 7.1.3.1 Error Estimates for Approximating [.11 1:11 \/3 + x + y dxdy =

 

 

 

 

 

 

 

 

 

 

 

6.859 942 640 334 65 when h = k = 2 (n = m = 1).

Rule nfe 2:2; Error affirm , lEst /Errl Order
Trapezoidal 4 .001 1 .60—1 6 .67-1 4 1
Midpoint 1 .000 -6.83-2 333—1 5 1
EM143 5 .002 775—4 250—1 323 3
Ewing 5 .001 775-3 250-1 32 3
DF543S 9 .003 —2.o3—3 417—2 20 3
Simpson 9 .002 142—3 417—2 29 ' 3
MINTOV i7 .007 —2.61-3 450-1 172 5
CSA 21 .009 727—4 1.25-2 i7 5
SC9CSS 21 .008 153-4 125-2 81 5
Lyness 9 .002 451-4 1.10-1 244 5
Radon 7 .002 4.394 1 .48-1IL 1 065 5
Radon 7 .002 -1.39—4 3.56-21 256 5
Radon 7 .002 4394 5.73.11” 14 696 5

 

 

 

 

*See Table 6.4.1.
iComputed by Stroud [49].

At first glance, Radon appears more accurate and efficient than MINTOV. Indeed, for this
example, one application of Radon does win over MINTOV. In fact, Simpson and Lyness also win
over MINTOV.

However, as the grid size is decreased, the efficiency of MINTOV becomes evident. This may be
seen in Table 7.1.3 .2 where n = m = 10. Here Radon uses 2.6 times as many function evaluations as
MINTOV and yet produces approximately the same error. MINTOV is now more efficient and
accurate than the third-order Simpson rule.

The accuracy and efficiency of the third-order cubature formula EM143 are somewhat surprising.
For n = m = 10, EM143 uses only 1/3 as many function evaluations as Simpson and yet produces a

smaller error.

132

Table 7.1.3.2 Error Estimates for Approximating L: L: V3 + x +y dxdy =
6.859 942 640 334 65 when h = k = 1/5 (n = m = 10).

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Rule rife (12:; I Error [5:332:16 lEst/Errl Order
Trapezoidal 121 .029 152—3 667—3 4 1
Midpoint 100 .025 -7.S9—4 333-3 4 1
EM I 43 140 .046 192-7 250-5 130 3
Ewing 221 .054 1.16-6 250—5 3
DF543S 225 .056 186-7 4.17—6 22 3
Simpson 441 .106 1.95-7 4.17-6 21 3
MINTOV 269 .079 —6.69-—9 450-7 67 5
CSA 309 .100 -l.28—10 1.25-8 98 S
SC9C5S 489 .131 -1.70-10 125—8 74 5
Lyness 621 .136 162-9 1.10-7 68 5
Radon 700 .148 -5.55-10 [ — — 5

 

 

 

The error estimates given by Stroud [49] for Radon have no provision for estimating the grid size
h X k which guarantees a prescribed error. The MINTOV error estimate (7.1.3.2) as well as those
given in Table 6.4.1 enjoy this advantage over those given by Stroud. Moreover, our error estimates are
much easier to obtain and apply than Stroud’s estimates.

Now we compute the grid size h X k which minimizes the number of function evaluations,
2(nm) + 3(n +m)+ 9, and which guarantees that the MINTOV truncation error (7.1.3.2) is smaller than
10"“. The results are presented in Table 7.1.3.3.

For a function such as g(x,y) = e’y sin (100x) which is changing more rapidly with respect to one
variable than the other, one would take h < k.

However, in the case of (7.1.3.1), it is more convenient to set h = k in (7.1.3.2). With this

simplification, the MINTOV truncation error is guaranteed to be smaller than a prescribed E > 0 by taking

11 = m > 0.876 541/6). (7.1.3.4)

133

Table 7.1.3.3 I}, 1.11 ¢3 +3: + y dxdy = 6.859 942 640 334 65.
Guaranteed MINTOV Error = 10‘“, a = 1(1)] 2. (Grid size h =# k.)

 

 

 

 

 

Absolute Time _ - Error
Error n m nfe Approxrmation Error _ lEst/Errl
Guaranteed (sec) Estimate
10'1 l 2 22 .009 6.860 4.73 009 882 55 -530-4 747-2 141.
10'2 2 2 29 .010 6.860 014 219 003 29 -7.l6-5 703-3 98.
10'3 3 3 45 .016 6.859 95g 207 911 96 -7.57-6 617-4 82.
10'4 4 5 76 .026 6.859 943 _428 653 38 -7.88-7 580-5 74.
10‘5 6 6 117 .037 6.859 942178 121 94 -1 .38-7 9.65-6 70.
10‘6 8 10 223 .066 6.859 942 653 717 00 -1 .34-8 906-? 68.
10‘7 12 14 423 .119 6.859 942 64 .1. 780 92 445-9 9.63-8 67.
10'8 18 20 843 .227 6.859 942 640 481 64 -1.47-10 9.71.9 66.
10'9 25 31 1 727 .450 6.859 942 640 339 67 -1.50-1l 994-10 66.
10"10 37 45 3 585 .910 6.859 942 640 336 07 -142-12 998-11 70.
10'11 56 64 7 537 1.882 6.859 942 640 334 21* 4.43-13 988-12 22.
10'12 82 94 15 953 3.936 6.859 942 640 33_3_ 46* 1.20-12 994-13 1.

 

 

 

 

 

 

 

 

 

 

 

 

 

*Contaminated by roundoff error.

The results are given in Table 7.1.3 .4 and are similar to those presented in Table 7.1.3.3, except
the cost of setting h = k requires a few more function evaluations. In Figure 7.1.3.2, the grid size
h = k versus the log of the absolute error is plotted for several cubature rules.

Thus for n = m = 41 or I: = k as .05, the actual MINTOV error is -1.35 X 10'”. The inequality

(7.1.3.3) provides a guaranteed error estimate of 9.47 X 10"“.

134

Table 7.1.3.4 1:115/3 +x+y dxdy = 6.859 94264033465

Guaranteed MINTOV Error = 10'“, a = 1(1)12. (Grid Size h = k.)

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Absolute Time Grid Error
Error n=m nfe Size Error . IEst/Errl
Guaranteed (sec) h=k Estimate
10'1 2 29 .011 l -7.16-5 7.03-3 98
10'2 2 29 .01 1 1 -7.l6-S 7.03 -3 98
10'3 3 45 .016 .667 -7.57—6 617-4 82
10'4 5 89 .030 .400 -4.01-7 288-5 7 2
10‘5 6 1 i7 .037 .333 -1.38-7 965-6 70
10‘6 9 225 .068 .222 -1.25-8 847—7 68
10'7 13 425 .120 .154 4.409 9.32-8 67
10‘8 19 845 .227 .105 -1.45-10 957—9 66
10’9 28 i 745 .456 .0714 -1 .41-11 9.34-10 66
10'10 41 3 617 .919 .0488 —i.35-1 2 9.47-11 70
10'11 60 7 569 1.890 .0333 466-13“ 9.65-12 21
10‘12 88 16 025 3.956 .0227 122-12* 969—13 1

 

*
Accuracy affected by machine limitation of 15 significant digits.

 

135

 

 

 

 

 

10‘
_ ‘l 1
T f f (‘3 +——_x +y dxdy = 6.859 942640 334 65
_ -1 -1
10° _.
x F
II
.C
11.1
E p
U)
9
C
0
10“ ..
- 3 617
o MINTOV
13 5119011 ESTIMATE
; FUNCTION EVALUATIONS
10-2 i i 1 i i i i l i 1 i i i
10'“ 10“? 10“° 10"3 10“ 10‘4 10‘2

LOG (ABSOLUTE ERROR)

. . . . 1
Figure 7.1.3.2 Error Curves in Approx1mating f: 1:1 V 3 + x +y dxdy

136

A RATIONAL FUNCTION WITH A SINGULARITY APPROACHING THE DOMAIN OF

INTEGRATION

7.1.4

Consider the double integral

=0.

1/4 [(w +4)ln(w+ 4) - 2(w+2)ln(w+2) + W In (w)] , w > 0
ln(2) (Cauchy Principal Value), w

l

The errors obtained when various cubature formulas are applied to this integral with h

dxdy
4(w + 2 + x + y)

I 1
III

1

(7.1.4.1)

k =1/50 or

and 0 are given in Table 7.1.4.1.

Q

100 as the parameter w takes on the values 1, .5. .1, .05, .01, .001

n=m=

Squire and the derivative corrected formula, EM143, perform quite well as w approaches 0. Simpson,

function evaluations for the error

Tanimoto, Lyness, Radon, Gauss, and Boole require too many

returned.

 

2

%(2 +x +_v)'l on (-l.1]

Graph of 2

Figure 7.1.4.1

137

Table 7.1.4.1 Application of Various Cubature Formulas to (111.: [4(w+ 2+x+y)] 'ldxdy

withh =k = 1/50 (r1 =m= 100).

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Avg w
Rule rife Time Order
(sec) 1.0 0.1 0.01 0.001 0.0
Trapezoidal 10 201 1.755 -8.89-6 -155—4 -1.81-3 -231-2 * l
Midpoint 10 000 1.771 4.44—6 7.71—5 670-4 187-3 249-3 1
EM143 10 400 1.976 -5.1 8-11 -5.32-8 ~9.48—6 125-4 444-4 3
Squire 20 200 3.491 -2.22-6 -3.84-5 ~2.87-4 -3.12-4 3.12—6 3
MM443 20 600 3.696 1.43-10 1.49-7 5.29-5 560—4 103-3 3
Ewing 20 201 3.527 -3.1 1-10 -3.27-7 -1.57-4 -6.46-3 * 3
DF543S 20 205 3.529 -5.17-11 -4.94-8 1.20-4 2.71—1 :1: 3
DM543A 20 601 3.731 -1.38-10 -1.45-7 -5.88-5 -2.07-3 * 3
Tyler 30 200 5.263 779-1] 814-8 321-5 4.15-4 8.33—4 3
XM543T 30 600 5.468 1.30-10 135—7 487-5 5.31-4 989-4 3
Miller 30 401 5.247 4.67-10 490-7 222-4 7.29-3 * 3
Simpson 40 401 7.019 -5.l8-1 1 -5.48-8 ~3.1 1-5 —1 .88-3 * 3
MINTOV 20 609 3.735 1.26-13 437-9 109—4 135—1 * 5
CSA 21 009 3.940 821-14 525-11 -4.83-5 —2.66-1 * 5
SC9C5S 40 809 7.227 8.28-14 103-10 -1 .17—5 -8.36-2 * 5
Tanimoto 41 209 7.390 8.41-14 .23-10 3.00—6 -1.08-2 * 5
Lyness 60 201 9.140 106-14 -1 .05-9 --1 .55-5 -1.78-3 * 5
Radon 70 000 10.034 506-14 359-10 4.43—6 184—4 509-4 5
Gauss 90 000 12.458 392-14 368-1] 844-7 871-5 338-4 5
Boole 160 801 24.881 4.95-14 -3.79-1 l -9.66-7 —2.73-4 * 5

 

*
Cubature rule not applicable due to singularity at (-1. -1) when w = 0

7.1.5 THE COMPARISON OF 45 NEW CUBATURE FORMULAS WITH 12 CONVENTIONAL RULES
In this example we present the results of applying each of the 52 cubature formulas listed in

Table 6.4.1 to the integral

1 1
I J. 112“} + 1)sin (7(1')(].\'(1_1' : e/n. (7.151)
O 0

138

 

Figure 7.1.5.] Graph ofz = 111(e" + 1) sin(1ry) on [0,1]2

For the sake of comparison, several additional cubature formulas are included. Error estimates are
also provided.

The approximations were computed by a single computer program which contained a subroutine
to calculate the 6 cubature elements, F0, FV, FM, FVI, FMl, F V1 1, using 4nm + 6(n + m) + 9 fe.
The cubature elements were then combined as indicated in Table 6.4.1 to produce the 52 approxima-
tions to [(f). The results are tabulated in Table 7.1.5.1.

For comparison we include the results of 5 additional fifth-order cubature formulas: Tanimoto,

Lyness, Radon, Gauss, and Boole.

139

 

 

 

 

 

 

 

Table 7.1.5.1 The Comparison of 45 New Cubature Formulas with 12 Conventional Rules for
the Double Integral fol fol 12(e" + l)sin (11y)dxdy = e/rr
h=k=1(n=m=l) h=k=1/10 (n=m=10) Or-
No“ Rule . der
rife Time Error Error Estl nfe Time Error Em” I3331'
(sec) Estimate"I Err (SOC) Estimate‘ E11

1 Midpoint 1 .000 -4.59-1 821-1 2 100 .047 -3.34-3 821-3 2 l
2 EM143 5 .002 -1.48-1 245-1 2 140 .063 -l.l3-5 245—5 2 3
3 ET183 9 .003 -8.48-2 3.38-1 4 144 .065 -5.64-6 3.38-5 6 3
4 EX183S 9 .003 -1.39-1 2.22-1 2'. 144 .065 -1.03-S 222—5 2 3
5 EC1C3S 13 .004 -l.32-l 222-1 2 148 .066 —l.03-5 2.22-5 2 3
6 ESIC3S 13 .005 —1.38-1 222-1 2 184 .081 -1.03-5 222-5 2 3
7 Trapezoidal 4 .001 865-1 1.64 2 121 .054 668—3 1.64-2 2 1
8 TM443 8 .002 243-1 4.40-1 2 161 .070 1.93-5 4.40-5 2 3
9 TT483 12 .004 1.17—1 3.47-1 3 165 .071 806-6 3.47-5 4 3
10 TX483S 12 .003 1.68-1 253-1 2 165 .071 1.18-5 2.53-5 2 3
11 TC4C38 16 .005 1.54-1 2.5 3-1 2 169 .073 1.18—5 253-5 2 3
12 TS4C3S 16 .005 1.59—1 2.5 3-1 2 205 .087 1.18-5 253-5 2 3
13 Squire 4 .001 1.50-1 4.11-1 3 220 .100 166-3 4.11-3 2 l
14 MM443 8 .003 -4.99—3 3.91-2 8 260 .1 l6 -2.03-7 391-6 19 3
15 MT483 12 .004 -3.67-2 109—1 3 264 .118 -3.02-6 109-5 4 3
16 MX483S 12 .004 4.38-3 1.5 8-2 4 264 .118 7.34-7 1.58-6 2 3
l7 MC4C3S 16 .005 8.28-4 1.5 8-2 19 268 .119 732-7 1.5 8-6 2 3
18 MS4C38 16 .006 5.5 7-3 1.58-2 3 304 .134 735-7 158-6 2 3
l9 Ewing 5 .001 -l.77-2 1.10-1 6 221 .101 -l.08—6 1.10-5 10 3
20 DF543S 9 .002 -3.64-2 6.34-2 2 225 .103 -2.95-6 634-6 2 3
21 DM54 3A 9 .003 —1.05-1 169—1 2 261 .117 -7.87-6 169-5 2 3
22 DM543B 9 .003 3.46-2 7.45-2 2 261 .117 3.00-6 7.45-6 2 3
23 DT583A 13 .004 271-2 4.22-2 2 265 .119 1.97-6 4.22-6 2 3
24 DT583B 13 .004 923-3 186-2 2 265 .119 7.51-7 186-6 2 3
25 DX585 13 .004 4.5 7-3 1.77-2 4 265 .119 3.489 1.77 -8 5 5
26 MINTOV 17 .005 1.73-3 1.14-2 7 269 .120 1.40-9 1.14-8 8 5
27 DSSCS 17 .006 7.81-4 929-3 12 305 .135 704-10 929—9 13 5
28 CSA 21 .007 -2.06-3 296-3 1 309 .136 -l.38-9 296-9 2 5
29 Tyler 5 .002 -5.27-2 8.66-2 2 320 .148 -3.89—6 866-6 2 3
30 XF543S 9 .003 -4.33-2 6.34-2 2 324 .149 -2.96-6 6.34-6 2 3
31 XM54 3T 9 .003 -1.4S-2 233-2 2 360 .164 -9.41-7 233-6 2 3
32 XT583A 13 .005 -5.81-2 897-2 2 364 .165 -4.l9-6 897-6 2 3
33 XTS 83B 13 .005 -3.99—2 732-2 2 364 . 165 -3 .19-6 7 .9 2-6 2 3
34 XX585 13 .004 —5.16-3 103-2 2 364 .165 -3.81-9 103-8 3 5
35 XCSCS 17 .006 -8.01-3 1.67-2 2 368 .167 -5.89-9 167-8 3 5
36 XSSCS 17 .006 -3.98-3 7.70—3 2 404 .181 -2.94-9 7.70-9 3 5
37 XHSGSS 21 .007 —1 .85-3 2.96-3 2 408 .183 -1.38-9 236-9 2 5

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

"' See Table 6.4.1

 

140

Table 7.1.5.] (cont’d.)

 

 

 

 

 

 

h=k=l (n=m=l) ll=k=l/10 1n=m=10)
No” Rule 3;;
Time . Error list - Time Error 1351
nfe (sec) I-.rror Estimate" It—rrI nie (sec) Error Estimate" E—rr
38 Miller 8 .002 -8.78-2 157-1 2 341 .154 -6.71-6 1.57-5 2 3
39 OF843S 12 .003 -5.03-2 634-2 1 345 .156 -2.96-6 6.34-6 2 3
40 OM843A 12 .004 226-2 4.22-2 2 381 .170 1.97-6 4.22-6 2 3
41 OM84BB 12 .004 -2.15-2 3.73-2 2 381 .170 -l.50-6 3.73-6 2 3
42 OT883T 16 .005 -4.69—2 932-2 2 385 .172 -3-76-6 932-6 2 3
43 OX885 16 .005 -6.55-3 135-2 2 385 .172 485-9 1.35-8 3 5
44 OC8C5 20 .006 - -9.40-3 1.98-2 2 389 .173 -6.93-9 138-8 3 5
45 OSSCS 20 .007 -4.66-3 9.29-3 2 425 .188 -3.46—9 929-9 3 5
46 OH9GSS 24 .008 -1.81-3 296-3 2 429 .189 -1.38-9 296-9 2 5
47 Simpson 9 .003 -4.10-2 6.34-2 2 441 .202 -2.95-6 6.34-6 2 3
48 SM945 13 .004 -2.85-3 507-3 2 481 .213 -2.07-9 507-9 2 5
49 ST985 17 .006 426-4 7.18-3 57 485 .219 991-12 7.18-9 724 5
50 SX985S I7 .005 -1.92-3 296-3 2 485 .219 -l.38-9 296-9 2 5
51 SC9CSS 21 .007 -l.98-3 2.96-3 1 489 .221 -1.38-9 2.96-9 2 5
52 SS9C‘SS 21 .007 -1.94—3 296-3 2 525 .235 -l.38-9 296-9 2 5
53 Tanimoto 25 .010 -1.95-3 — — 529 .236 -1.38-9 —- — 5
54 Lyness 9 .004 —8.604 226-3 3 621 .272 -6.26-10 226-9 4 5
55 Radon 7 .003 879-4 _ — 700 .303 6.38-10 — — 5
56 Gauss 9 .004 -6.01-4 — — 900 .385 4.14-10 - — 5
57 Boole 25 .010 6.18-4 — — 1681 .739 4-31-10 — — 5

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

*See Table 6.4.1

These results are included in the sense of a benchmark test to reveal possible errors in the cubature
formulas of Table 6.4.1. These and other results confirm the accuracy of the entries in Table 6.4.].

Finally, we observe that for the double integrals considered, for any reasonable grid size the
composite fifth-order cubature rules, CSA, Tanimoto, Lyness, Radon, Gauss, and Boole generally all
produced approximately the same error; however, CSA required far less function evaluations and thus is
the best fifth-order method.

Similar comparisons can be made for the powerful third-order derivative corrected rule EM 143.

141
7.2 TRIPLE INTEGRALS
Denoting the step size in each dimension by h]- = (bi - ai)/nl-, i = 1(1)3, the 3-dimensional

formulation of MINTOV is

b3 b2 61
f f j f(x, y, z)dxdydz
”3 02 01

=—h h2h3Z Z Ema, +h (11 My.) az+h (12-14), a3+h3 (13-14))

13:1 12=li1=l
"3 "2 "I
7 fit 7 "1' . . .
+ 1—2—6-h1h2h3 Z Z L f(l11 +11’111‘12 +12h2~ ”3 +’3h3)
i3=012=011=0
1 "3 ”2
_ _ 2 ' YT“ . .
240111112113 Z 2‘ Ifx(b1~02 +12h2~03+l3h31
I3=OIZ=O
— fx(al. (12 + izhzr a3 +1'3h3)]

"3 ”1
I 2 ' ' . .
'- jib—(1115,13 Z Z [f,(al+11hl. b2. 03 +l3h3)
13:0 11:0

' fy(al +1.1hl302~a3 +1'3’13)I (72.1)

"2 "I
I 2 ' ' . .
12:0 11:0

- fz(a1+ilh1. 02 'I'Izhz. 03)]

"3
l 2 2 ' . .
- fill—O- hlh2h3 ZO [fxy(al.a2.a3 +13h3)- fxy(al,b2.a3 +l3h3)
I3:
- fxy(b1,a2.a3 +i3h3) + fxy(b1,b2,a3 +i3h3)]
"2
.. I 2 2 ’ . .
144Oh1h2h3 Z [fxz(01,02 +12h2,a3) - fxz(al,02+12h2,b3)
12:0
_ fxz(b11a2 +i2h21a3) + fxz(b1.02 +i2h2,b3)]
"1
1,12 2
- mill}? 23h 20 [f),z(al +ilh1,02,a3)- fyz(al +iIhI 02,123)
I]:

The truncation error may be estimated by

(b1 ‘ 00092 ‘ ”211.193 ' ‘13)
604 800

 

lE(f)I < [(th? + 713/143

6 6 4 2 42 2 4 24
4 2 42 2 4 24 4 2 42
+ h1h3M13 + hlh3M13 + h2h3M23

2 4 24 2 2 2 222
+ h2h3 M23) + 280hlh2h3M123]

Table 7.2.1 lists the number of function evaluations required by various composite multiple

quadrature fomiulas. In this reSpect, MINTOV is superior to the conventional formulas.

Table 7.2.1 Nfe for Several Multiple Quadrature Formulas Repeated 711112173 Times

 

 

 

 

 

 

_ _ MINTOV Simpson Lyness Gauss Boole
"1 ”12-,” 2n3+9n2+27n+ i9 (2n+1)3 (n+1)3+7n3 (3n)3 (411311)3

1 57 27 15 27 64

2 125 125 83 216 729

4 399 729 573 l 728 4 913

8 1835 4913 4313 13824 35 937

16 10 947 35 937 33 585 110 592 274 625

32 75 635 274 625 265 313 884 736 2 146 689

64 562 899 2 146 689 2 109 633 7 077 888 16 974 593

100 2092 719 8120601 8 030 301 27 000 000 64 481 201

 

 

 

 

 

 

 

7.2.1 A FUNCTION WITH VANISHING MIXED HIGHER-ORDER PARTIAL DERIVATIVES

We wish to approximate the triple integral

2 2 2
J‘ J J ln(xyz)dxdydz = In 64 - 3.
I l 1

(7.2.1.1)

 

143

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Figure 7.2.1.1 Graph ofz =ln (xy) on [1,2]2

 

 

 

 

 

 

 

 

 

 

 

 

 

 

144

 

 

 
 
 
 
 
 
  

 

 

10
"' 2 2 2
_ If] ln(xyz)dxdydz = 1.158 883 08335967
1 1 1
10° _ o MINTOV
: c1 EWING
L—
o MIDPOINT
" o TRAPEZOIDAL
M
a: r-
I
N
8
I h
‘P
I“
'3
In
9 I-
E
0
10‘1 1.
P
P
F75 635
* FUNCTION
EVALUATIONS
10'2 1 1 1 1 1 i i 1 1 1 i i
10'“ 10"? 10"° 10'8 10'6 10“ 10'2 10°

LOG (ABSOLUTE ERROR)

Figure 7.2.1.2 Error Curves in Approximating (2 If f: ln(xyz)dxdydz

145

The first-order partial derivatives of the function f(x,y) = ln (xyz) are much less costly to
evaluate on a computer than is f. Moreover. the mixed second-order partials vanish everywhere. Thus it

is of interest to apply MINTOV to this function. The results are presented in Table 7.2.1.].

Table 7.2.1.1 [12112112 In (xyz)dxd_vdz = 1.158 883 083 359 67

 

 

 

 

 

 

 

 

 

hl=h2=h3=1 hl=h2=h3=I/IO
Or-
Rule der
rife Time (sec) Error nfe Time (sec) Error

Trapezoidal 8 .002 1.19-1 1331 .457 125-3 1
Midpoint l .000 —5.75-2 1000 .361 .6,24.4 1
Ewing 9 .003 138-3 2331 .818 132.7 3
MINTOV 57 .009 -6.41-5 3189 .933 _1 .1440 5

 

 

 

 

 

 

 

 

 

For hi = h2 = 113 = 1/10, the MINTOV approximation is better than Ewing’s result by 3 orders
of magnitude. This confirms the theoretical result that the addition of partial derivative correction

terms can greatly enhance a numerical integration formula.

7.2.2 MINTOV vs. JPL’s MQUAD

Bunton, Diethelm, and Haigler [8] give the following example:

11/2 11/2 1r/2
J. I f cos (x) cos (13) cos (z)dxdydz = 8.0. (7.2.2.1)
-1r/2 -1r/2 -11/2

The results of the DC-MQUAD algorithm using n1 = n2 = 713 = 2,3,5,8,° ° 3 , (cf. the Fibonnaci

sequence) and JPL’s routine MQUAD are presented in Table 7.2.2.1. Again, our method based on the

3-dimensional MINTOV formula (7.2.1) is superior to MQUAD.

146

Table 7.2.2.1 MINTOV vs. MQUAD for the Integral 1:; : [Zr/22 L1];

Relative Errors Requested = 10’“, a = 1(1)10.

cos (x) cos (y) cos (z)dxdydz = 8.

 

 

 

 

 

 

 

 

Relative MINTOV (DC'MQUAD Algorithm) MQUAD (JPL Routine)
Error .

Requested rife Time (sec) R2128 rife Time (sec) R2132?
10"1 360 .198 -1.1 1-3 2 197 1.634 8.74-5
10'2 989 .546 -5.07-5 2 197 1.635 8.74-5
10“3 2 824 1.529 -3.00—6 2 197 1.627 874-5
10'4 2 824 1.529 -3.00-6 24 389 17.641 4.13-6
10‘5 9 109 4.971 -1.63-7 35 937 25.693 9.70—8
10'6 32 186 17.558 -9.14-9 117 649 82.090 1.63-9
10'7 122 135 66.611 -5.06-10 614 125 421.889 259-11
10'8 122 135 66.611 -5.06-10 4 173 281 2 844.152 407-13
10’9 483 614 264.343 -2.69—1l 4 173 281 2 844.114 407—13
10'10 1 967 263 1 076.794 155-11 31 855 013 21 647.323 8.40-15

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

.//LL“11L

Figure 7.2.2.1 Graph of z = cos (x) cos (v) on [-1r/2, 11/2] 2

 

 

 

 

 

 

 

 

“a

 

 

 

 

 

 

 

 

LLrul

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

147

7.2.3 AN EXAMPLE ILLUSTRATING THE PERSISTENCE OF FORM

Consider the integral

 

f'lzf'flf'” 1 “Kim Wadydz (7 23 1)

sin (x) sin (y) sin (z)e -

= 1.531 670 226 93.

If we let w2 = x2 + y2 + 22 then the first-order partial with respect to x is

 

 

sin sin sin “9
fx(x.y.z) = [(1 +w) cos(x) - (l +w+x2) x(x)] 0: (z)e (723.2)
yz
and the mixed second-order partial with respect to x and y is
+ 2+ 2+ 2 + 2 2
from.» = [(1 + W)c08(x)costy) + w “’ ”(3” y) " ” sin(x)sin(y)

(7.23.3)

 

 

.. l + w: x2 sin (x) cos (y) - 1 +wy+y2 cos (x) sin(y)]——L-$in (z 640'

xyz

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Figure 7.2.3.1 Graph ofz = sin(x) sin (y) e' V ‘2 ”2 on [-11/2,1r/2] 2

l +\/.1:I -1-yI
xy

148

The results of applying MINTOV with various grid sizes are presented 111 Table 7.2.3.1.

_ 1+ , .
Table 7.2.3.1 10”” 10”” 10"” “ sin(x)sin(1')sin(:)e’“dxd_l'dz =

 

1.531670 226 93.

xyz

 

 

 

 

 

Ewing MINTOV
nl=n3=n3
nfe Time (sec) Error rife Time (sec) Error
1 9 .007 -1 .54-2 57 .048 -5 .30-3
2 35 .031 -1.04-3 125 .111 -8.75-5
4 189 .174 —6.59-5 399 .364 -1.36-6
1241 1.159 -4.13-6 1835 1.712 -2.l3-8
16 9 009 8.486 -2.58-7 10 947 10.314 -3.66- 10
32 68 705 64.978 -1 .62-8 75 635 71.579 -3.76-1 I
64 536 769 509.407 -1.03-9 562 899 534.424 -2 .00-1 1

 

 

 

 

 

 

 

 

 

 

 

This example illustrates the persistence of form concept, namely, the cost of evaluating a

derivative is approximately the same as the cost of evaluating the function.

In Table 7.2.3.1 for a given mesh size. the difference between the number of function evaluations
required by MINTOV and Ewing is the number of partial derivative evaluations required by MINTOV.
For example, for n] = n2 = r13 = 16, MINTOV uses 1938 first- and second-order partial derivative
evaluations (pde) and takes 1.841 seconds, or an average of 950 microseconds per partial derivative
evaluation (us/pde). This is to be compared with the conventional Ewing’s formula which uses 9009 fe
in 8.438 seconds, or an average of 937 us/fe. MINTOV averages 939 microseconds per evaluation,
whether function or partial derivative.

Thus throughout this study we‘have weighted a partial derivative evaluation the same as a function

evaluation, and refer to either simply as a function evaluation.

149

 

         

 

 

 

10‘
h .3. 1 1
’- 2 2 21+ (x2+y§+22
._ f f f ———————sin(x)sin(y)sin(z)e' I” T V I z dxdydz = 1.531 67022693
xyz
0 0 0
14
10° _ 0 MINTOV
: 13 EWING
_ o MIDPOINT
:9 ” FUNCTION
I EVALUATIONS
.cN _
II
zp
... L
'2‘
(I)
9
t
o F-
10"_
~ 75635
)—
536769
L 4243841 2097152
10‘2 1 l i 1 1 1 i i 1 1 i 1 i
10“ to"2 10'” 10" to" to" 10‘2 10°

LOG (ABSOLUTE E RROR)

Figure 7.2.3.2 Error Curves

8. CONCLUSIONS AND RECOMMENDATIONS

We have studied the problem of enhancing the accuracy of conventional formulas for evaluating
multiple integrals numerically over d-dimensional rectangles by the addition of partial derivative
correction terrr'is evaluated on the boundary of the domain of integration.

The formulas in Chapter 5 were based on the double Euler-Maclaurin Summation formula (5.2.1)
and were found somewhat cumbersome to apply to practical situations because of the different weights
for different nodes.

The formulas in Chapter 6, based on the finite Taylor Series expansion of the integrand. are
much easier to apply in practice. We have constructed MINTOV (6.2.17), a fifth-order multi-
dimensional integration formula which may be considered as the d-dimensional generalization of the
l-dimensional Lanczos quadrature formula (6.2.19) or of Ewing’s cubature formula D0503 given in
Table 6.4.1.

For a single integral, the derivative correction terms are evaluated only at the end points Of the
interval Of integration. The situation is somewhat more complicated in higher dimensions since, as the
dimension increases, the boundary becomes increasingly more complex.

Of greater significance is the fact that in higher dimensions. most of the volume of a d-rcctangle
lies near the boundary. We have accounted for this by constructing multidimensional integration
formulas with boundary partial derivative correction terms. the number of which increases as the
dimension increases.

Indeed, as can be seen in (6.3.2), the d-dimensional MINTOV with n subdivisions requires 71" fe
at the centroids of the subregions, (n + I)“ re at lattice points. 2d(n + l)"‘1 first-order partial
derivative evaluations at the lattice points of the 2d “faces" or (d—1)-dimensional hyperplanes. B,
which bound the domain of integration, and 2d(d-1)(n + l)‘1"2 mixed second-order partial derivative .
evaluations at the lattice points of the 2‘1 “edges” or (d - 2)-dimensional hyperplanes which bound B.

In Chapter 6. 47 hitherto unpublished cubature formulas with boundary partial derivative

correction temis were given. Computable bounds for the truncation errors were also given.

150

151

Three integration formulas, DF543S, XF543S, and OF843S (see Table 6.4.1) were discovered
which require mixed second-order partial derivative correction terms evaluated only at the corners Of
the rectangular domain Of integration. These formulas have the same error bound as Simpson’s rule;
however, they require far fewer function evaluations than Simpson’s rule. Of the three, DF543S is the
most efficient, requiring only half as many function evaluations as Simpson’s rule.

OF543S has no interior nodes and yet it is not as efficient as DF543S which has one node interior
to each cell. Thus, it is not always advisable to take as many points as possible on the boundary.

In cases where a third-order rule is considered adequate, DF543S should be preferred to Simpson’s
rule.

The fifth-order MINTOV was compared with JPL’s routine MQUAD, and on the two test integrals
considered, MINTOV required fewer function evaluations to achieve a prescribed relative error than did
MQUAD.

The numerical results presented indicate that formula CSA or DH5G5S is an efficient and accurate
composite fifth-order integration rule. For a sufficiently small grid size (n = m = 4 for a double integral,
n = m = 2 for a triple integral. etc.) CSA requires fewer function evaluations than Simpson's rule.

Therefore, in situations where first- and mixed second-order partial derivatives are easily calculated.
we have shown that the use of nodes satisfying the inclusion property coupled with the use of derivative
correction terms exhibiting the equal weight-alternate sign property can improve the accuracy and
efficiency Of integration formulas. That is, whenever the first- and mixed second-order partial derivatives
are easily computed. the formulas of Chapter 6 are to be preferred to conventional composite integration
rules of comparable degrees.

Finally, future investigations should include the addition of weight. functions. more general domains

ofintegration. and the use of “difference correction terms."

LIST OF REFERENCES

152

GENERAL REFERENCES

Abramowitz, Milton
On the practical evaluation of integrals. J. Soc. lndust. Appl. Math 2 (1954) 20-35.

Abramowitz, Milton and Stegun, Irene A.
Handbook of Mathematical Functions. National Bureau of Standards. Applied Math. Series. V55 (1964)

Aitken. A.C. and Frewin. G.L.
The Numerical Evaluation of Double Integrals. Proc. Edinburgh Math. Soc. 42 (1923) 2-13.

Anders, Edward B.
An Extension of Romberg Integration Procedures to N-Variables. J. ACM 13 (1966) 505-510.

Bailey, Carl B., Jones, Randall E.
Usage and Argument Monitoring of Mathematical Library Routines. ACM Trans. 011 Math. Software. 1
(1975) 196-209.

Baker, Christopher TH. and Hodgson, Graham S.
Asymptotic Expansions for Integration Formulas in One or More Dimensions. SIAM J. Numer. Anal.
8(197l)473-480.

Barnes, E.W.
The Generalisation of the Maclaurin Sum Formula, and the Range Ofits Applicability. Quart. J. of
Pure and Appl. Math. 35 (1904) 175-188.

Barrett. W.
On the Convergence of Cote’s Quadrature Formulae. J. London Math. Soc. 39 (1964) 296-302.

Baten, William D.
A Remainder for the Euler-Maclaurin Summation Formula in Two Independent Variables. Amer. J.
Math. 54 (1932) 265-275.

Bierens de Haan, David
Supplément aux Tables D’lnte’grales Définies. Amsterdam: G.G. Van Der Post (1864).

Camp. C.C.
Note on Numerical Evaluation of Double Series. Ann. Math. Stat. 8 (1937) 72-75.

Chakravarti, P.C.
Integrals and Sums. Some New Formulae for their Numerical Evaluation. London: The Anthlone Press
(1970).

Cody, WJ.
The Construction of Numerical Subroutine Libraries. SIAM Rev. 16 (1974) 3646.

Cranley, R. and Patterson, T. N. L.
The Evaluation Of Multidimensional Integrals. Comp. J. 1 1 (1968-69) 102-110.

Davis, Philip J. and Rabinowitz, Philip
Methods of Numerical Integration. New York: Academic Press (1975).

Frame, J.S.
Numerical Integration. Amer. Math. Monthly 50 (1943) 244-250.

153

Frank, Irving
An Application of the Euler-Maclaurin Sum Formula to Operational Mathematics. Quart. Appl. Math.
20 (1962) 89-91.

Gauss. C.F.
Metliodus Nova Integralium Valores per Approximationem lnveniendi. Carl Friedrich Gauss Werke.
GOItingen: K611igliclien Gesellschaft der Wissenschaften 3 ( 1866) 163-196.

Ghizzetti. A., and Ossicini. A.
Quadrature Formulae. New York: Academic Press (1970).

Gibb. David
Interpolation and Numerical Integration. London: G. Bell & Sons. Ltd. (1915).

Gould, H.W. and Squire, William
Maclaurin's Second Formula and Its Generalization. Amer. Math. Monthly 70 (1963) 44-52.

Gradsliteyn, LS. and Ryzhik, I.M.
Tables of Integrals. Series. and Products. (Trans. from 4th Russian ed. by Alan Jeffrey). New York:
Academic Press (1965 ).

Grant. J.A.
Derivation of Correction Terms for General Quadrature Formulae. Proc. Cambridge Phil. Soc. 66 (1969)
571-586.

Guenther. RB. and Roetman. E.L.
Newton-Cotes Formulae in n-Dimensions. Numer. Math. 14 (1970) 330-345.

Haber. Seymour
Numerical Evaluation of Multiple Integrals. SIAM Rev. 12 (l970)481-526.

Hammer, Preston C. and Wicke, Howard H.
Quadrature Formulas Involving Derivatives ofthe Integrand. Math. Comp. I4 (1960) 3-7.

Hamming, R.W. and Pinkham, R.S.
A Class of Integration Formulas. J. ACM 13 (1966) 430-438.

Hgvie, T.
Derivation of Explicit Expressions for the Error Terms in the Ordinary and the Modified Romberg
Algorithms. BIT 9 (1969) 18-29.

Hfivie, T.
Some Algorithms for Numerical Quadrature Using the Derivatives of the Integrand in the Integration
Interval. BIT 10 (1970) 277-294.

Hildebrand, F.B.
Introduction to Numerical Analysis. New York: McGraw-Hill Book Company. 211d Ed.. (1974).

Hillstrom, K.E.
Comparison of Several Adaptive Newton-Cotes Quadrature Routines in Evaluating Definite Integrals
with Peaked Integrands. ANL-7511. Argonne Natl. Lab. Argonne, II. (1968).

Hummel, PM. and Seebeck, Jr., CL.
A Generalization Of Taylor’s Expansion. Amer. Math. Monthly 56 (1949) 243-246.

154

Irwin, .10.
On Quadrature and Cubature or 011 Methods of Determining Approximately Single and Double Integrals.
London: Tracts for Computers. No. 10. Cambridge: Cambridge Univ. Press. (1923).

Isaacson, E. and Keller. [1.8.
Analysis of Numerical Methods. New York: John Wiley & Sons, Inc. (1966).

Johnson. W. Woolsey
011 Cotesian Numbers: Their Ilistory. Computation and Values to ”=20. Quart. J. Pure Appl. Math.
46 (1915) 52-65.

Kahaner. David K.
Los Alamos Workshop on Quadrature Algorithms. SIGNUM Newsletter 11 (1976) 4-26.

Knuth. Donald E. and Bucklioltz. Thomas J.
Computation of Tangent. Euler. and Bernoulli Numbers. Math. Comp. 21 (1967) 663-688.

Krylov, VI and Pal‘tsev. A.A.
Tables for Numerical Integration of Functions with logarithmic and Power Singularities. (Trans. frorr
lst Russian ed., 1967. by the Israel Program for Scientific Translations Staff). Jerusalem: Keter Press
(1971).

Krylov. VI and Si'llgina. LT.
Handbook ofNumericaJ Integration. (In Russian). Moscow: Izdat. Nauka. (1966).

Lambert, JD. and Mitchell. A.R.
The Use of Higher Derivatives in Quadrature Formulae. Comp. J. 5 (1962-63) 322-327.

Lambert, John D. and Mitchell, Andrew R.
Repeated Quadratures Using Derivatives of the Integrand. Z. Aiigew. Math. Phys. 15(1964)84-90.

Lehmer. D.H.
An Extension of the Table of Bernoulli Numbers. Duke Math. J. 2 (1936) 460-464.

Lipow. Peter R. and Stenger. Frank
Ilow Slowly Can Quadrature Formulas Converge? Math. Comp. 26 (1972) 917-922.

Lyness. IN. and McHugh, B.J.J.
Integration Over Multidimensional Hypercubes. I. A Progressive Procedure. Comp. J. 6 (1963-64)
264-270.

Lyness. J.N. and Ninham. B.W.
Numerical Quadrature and Asymptotic Expansions. Math. Comp. 21 (1967) 162-178.

McNamee. John
Error-Bounds for the Evaluation of Integrals by the Euler-Maclaurin Formula and by Gauss-Type
Formulae. Math. Comp. 18 (1964) 368-381.

Maxwell, .1. Clerk
On Approximate Multiple Integration between Limits of Summation. Proc. Cambridge Phil. Soc. 3
(1877) 3947.

Miller, J.C.P.
Quadrature in Terms of Equally-Spaced Function Values. U.S. Army Math. Research Center. Madison.

Wisc. Rep. I67,July. 1960. 1-91.

Munro, W.D.
Note on the Euler-Maclaurin Formula. Amer. Math. Monthly 65 (1958) 201-203.

Mysovskih. I.P.
Lectures on Numerical Methods. Groningen, The Netherlands; Wolters-Noordhopf Publishing (1969).

Nikol‘skii. S.M.
Quadrature Formulae. (Trans. from the lst Russian ed., 1958) Delhi: Hindustan Publishing Corp. (1964).

Nortliam. Jack 1.
Certain Summation and Cubature Formulas. Iiast Lansing: Michigan State University. Master‘s thesis
(1939).

Ohm. M.
Etwas Uber die Bernoullischen Zalilen. J. Reine Angew. Math. 20 (1840) I l-l2.

Patterson. T.N.L.
Integration Formulae Involving Derivatives. Math. Comp. 23 (1969)411412.

Philips. G.M.
Numerical Integration in Two and Three Dimensions. Comp. J. 10(1967) 202-204.

Ralston. A.
A Family of Quadrature Formulas Which Achieve High Accuracy in Composite Rules. J. ACM 6(1959)

384-394.

Rosenstock. Herbert B.
Euler-Maclaurin Formula in Three Dimensions. J. Math. and Phys. 43 (1964) 342-346.

Sack. R.A.
Newton-Cotes Type Quatrature Formulas with Terminal Corrections. Comp. J. 5 {1962-63) 230-237.

Sadowsky. Michael
A Formula For Approximate Computation of a Triple Integral. Amer. Math. Monthly 47 (1940) 539-

543.

Saidel. Frank
Some Interpolation Formulas in Two Variables. East Lansing: Michigan State University. Masters thesis
(I941).

Scarborough. J.B.
On the Relative Accuracy of Simpson’s Rules and Weddle’s Rule. Amer. Math. Monthly 34 (1927)

135-139.

Scarborough. James B.
Numerical Mathematical Analysis. Baltimore: The Johns Hopkins Press. 5111 ed. (1962).

Schoenberg. H. and Shanna. A.
The Interpolatory Background of the Euler-Maclaurin Quadrature Formula. Bull. AMS 77 (1971 ) 1034-

1038.

156

Shampine, Lawrence F.
Quadrature Formulas Using Derivatives. Math. Comp. I9 (1965) 481-482.

Smith, Francis J.
Quadrature Methods Based on the Euler-Maclaurin Formula and on the Clenshaw-Curtis Method of
Integration. Numer. Math. 7 (1965) 406-41 1.

Smith, J.M.
Recent Developments in Numerical Integration. J. Dynamic Systems. Measurement. and Control. March
(1974) 61-70.

Squire. William
Some Applications of Quadrature by Differentiation. J. Soc. Indust. Appl. Math. 9(1961)94-108.

Steffensen, J.F.
Interpolation. New York: Chelsea Publ. Co., 2nd ed. (1950).

Strb’m, Torsten
Strict Error Bounds in Romberg Quadrature. BIT 7 (1967) 314-321.

Stroud. A.H.
Quadrature Methods for Functions of More Than One Variable. Ann. NY Acad. Sci. 86 (1960) 776-791.

Stroud, A.H.
A Bibliography on Approximate Integration. Math. Comp. 15 (1961) 52-80.

Stroud, A.H.
Numerical Quadrature and Solution of Ordinary Differential Equations. New York: Springer-Verlag
(1974).

Stroud. A.H. and Secrest. D.
Gaussian Quadrature Formulas. Englewood Cliffs: Prentice-Hall, Inc. (1966).

Strubble. George
Tables for Use in Quadrature Formulas Involving Derivatives of the Integrand. Math. Comp. I4 (1960)
8-12.

Takeyama. Hisao
Expressions for Interpolation and Numerical Integration of High Accuracy. TOhOku Univ. Technology
Reports. 23 (1958) 47-70.

Thacher, Jr., Henry C.
An Efficient Composite Formula for Multidimensional Quadrature. Comm. ACM 7 (1964) 23-25.

Tyler, G.W.
Numerical Integration of Functions of Several Variables. Canadian J. Math. 5 (1953) 393-412.

USpensky, J .V.
On an Expansion ofthe Remainder in the Gaussian Quadrature Formula. Bull. AMS 40 (1934) 871-876.

Weddle, Thomas
On a New and Simple Rule for Approximating to the Area Ofa Figure by Means of Seven Equidistant
Ordinates. Cambridge Math. J. 9 (1854) 79-80.

f0

BIBLIOGRAPHY

Adams. J.C.
Table of the values of the first sixty-two numbers of Bernoulli. J. Reine Angew. Math. 85 (1878)
269-272.

Alilin, A.C.
On Error Bounds for Gaussian Cubature. SIAM Rev. 4 (1962) 25-39.

Albrecht. J., and Collatz, L.
Zur Numerischen Auswertung h’lellrdimensionaler Integrale, Z. Angew. Math. Mech. 38 (1958)
1-15.

Barnes. E.W.
The Maclaurin Sum-Formula. Proc. London Math. Soc. 3 (1905) 253-272.

Bauer, F.L., Rutishauser, H.. and Stiefel. E.
New Aspects in Numerical Quadrature. Proc. Symp. Appl. Math. Providence: Amer. Math. Soc.
15 (1963) 199-218.

Becker. George F.
Some New Mechanical Quadratures. Philos. Mac. 22 (191 ”342-353.

Bickley. W.G.
Finite Difference Formulae for the Square Lattice. Quart. J. Mech. Appl. Math. 1 (1948) 3542.

Bunton. Wiley R.. Diethelm, Michael. and Haigler. Karen
Romberg Quadrature Subroutines for Single and Multiple Integrals. Jet Propulsion Lab.. Pasadena.
Calif. TM-314-22l. July I, 1969,1-50.

Bunton, Wiley R., and Diethelm, Michael
Modifications to the JPL Romberg Subroutines. Jet Propulsion Lab., Pasadena. Calif. TM-314-247.
September 1. 1970, 1-9.

Bunton, Wiley R.. Diethelm, Michael. and Winje. Gilbert L.
Modified Romberg Quatrature: A Subroutine to Support General Scientific Computing. Jet
Propulsion Lab. Pasadena. Calif. TM-3l4-258. April I. 1970. 1-35.

Burnside, W.
An Approximate Quadrature Formula. Messenger ofMath. 37 (1908) 166-167.

Caselleto, J.. Pickett. M., and Rice. J.
A Comparison of Some Numerical Integration Programs. SIGNUM Newsletter. 4 (1969) 3040.

Davis. Philip J.. and Rabinowitz. Philip
Methods of Numerical Integration. New York: Academic Press (1975).

26

27

158

DeDoncker. Flise. and Piessens. Robert
A Bibliography on Automatic Integration. J. Comput. and Appl. Math. (to appear).

Euler. Leonhard
Methodus Generalis Summandi Progressiones. Commentarii Acad. Sci. Imp. Petropolitanae. Vol. 6
(1732-33, published 1738) St. Petersburg.

Ewing. G.M.
Oii Approximate Cubature. Amer. Math. Monthly. 48 (1941) 134-136.

Forsythe. George E.. and Moler. Cleve B.
Computer Solution of Linear Algebraic Systems. Englewood Cliffs: Prentice-Hall (1967).

Frame. J. Sutherland
Numerical Integration and the Euler-Maclaurin Summation Formula. East Lansing: Michigan
State University (to appear).

Fritsch. F.N.
A Bibliography on Approximate Multidimensional Integration 1960-1968. Lawrence Radiation
Lab., Livermore. Calif. UCRL-50610. March 5, 1969. 1-20.

Good. I.J.. and Gaskins, RA.
The Centroid Method of Numerical Integration. Numer. Math. 16 (1971) 343-359.

Hammer. Preston C.
Numerical Evaluation of Multiple Integrals. On Numerical Approximation. Langer. R.E. (Ed.)
(Proceedings of a symposium conducted by the U.S. Army Math. Research Center). Madison:
The University of Wisconsin Press (1959). 99-1 15.

Hammer. Preston C.. and Wymore. A. Wayne
Numerical Evaluation of Multiple Integrals I. MTAC II (1957) 59-67.

Ionescue, D.V.
Generalization of the Quadrature Formula of N. Ob reschkoff for Double Integrals (In Romanian).
Stud. Cerc. Mat. 17 (1965) 831-841.

Joyce, D.C.
Survey of Extrapolation Process in Numerical Analysis. SIAM Rev. 13 (I97l)435490.

Kahaner, D.K.
Comparison of Numerical Quadrature Formulas. Mathematical Software. I. R. Rice (ed).
New York: Academic Press (1970) 229-259.

Knopp, K.
Theory and Application of Infinite Series. London: Blackie and Son (1951).

Krylov, V.I.
Approximate Calculation of Integrals. (Trans. from lst Russian ed., 1959, by A. H. Stroud)
New York: Macmillan Company (1962 ).

Lanczos, C.
Applied Analysis. Englewood Cliffs: Prentice-Hall (1956).

30

31

33

34

35

36

37

38

39

40

41

43

44

159

Lyness. J.N.
Symmetric Integration Rules for Ilypercubes I. Iirror Coefficients. Math. Comp. I9 (1965)
260-276.

Lyness. J.N... and Kaganove. J.J.
A Technique for Comparing Automatic Quadrature Routines. Private communication (1975).

Lyness. J.N., and Kaganove, J.J.
Comments on the Nature of Automatic Quadrature Routines. ACM Trans. on Math. Software.
2(1976)65-81.

Lyness, J.N., and McHugh, B.J.J.
On the Remainder Term in the N-I)imensional Euler Maclaurin Expansion. Numer. Math.
15 (1970) 333-344.

Maclaurin. Colin
A Treatise on Fluxions. Edinburgh ( 1742). p. 672.

Meister, Bernd
011 a Family of Cubature Formulae. Comp. J. 8 (1966) 368-371.

Miller, J.C.P.
Numerical Quadrature over a Rectangular Domain in Two or More Dimensions. I Quadrature
over 3 Square using up to Sixteen Equally Spaced Points. Math. Comp. 14 (1960) 13-20.

Milne, W.E.
Numerical Calculus. Princeton: Princeton University Press (1949).

Mustard. D.. Lyness. J.N.. and Blatt. J.M.
Numerical Quadrature in N Dimensions. Comp. J. 6 (1963-64) 7587.

Ob reschkoff. N.
Neue Quadraturaturformeln. Abhanrllungen der Preussisclien Akademie der Wissenschaften.
Berlin. 4 (1940) 1-20.

Oliver. F.W.J.
Asymptotics and Special Functions. New York: Academic Press (1974).

Oliver. J.
The Evaluation of Definite Integrals Using High-Order Formulae. Comp. J. 14 (1971) 301-306.

Price, J.F.
Examples and Notes on Multiple Integration. Boeing Scientific Research Lab.. Seattle. Wash.
Mathematical Note No. 285, Dl-82-023l (1963). 1-34.

Radon. Johann
Zur Meclianischen Kubatur. Monatsh. Math. 52 (1948) 286-300.

Richardson. L.F.. and Gaunt. J.A.
The Deferred Approach to the Limit. Trans. Roy. Soc. London 226A (1927) 299-361.

Romberg, Werner
Vereinfaclite Numerische Integration. Norske Vid. Selsk. Forli. Trondheim 28 (1955) 30-36.

46

47

48

49

50

51

I60

Sheppard. W.E. ,
Some Quadrature-Formulae. Proc. London Math. Soc. 32 (1900) 258-277.

Simpson. T.
Mathematical Dissertations. London (1743).

Squire. William
Integration for Engineers and Scientists. New York: Amer. Elsevier Pub. Company (1970').

Squire, William
Numerical Evaluation of a Class of Singular Double Integrals by Symmetric Pairing. Intl. J.
for Numer. Meth. in Engr. 10 (1976) 703-708.

Stroud, A.H.
Approximate Calculation of Multiple Integrals. Englewood Cliffs: Prentice-Hall (1971 ).

Tanimoto, B.
An Efficient Modification of Euler-Maclaurin‘s Formula. Trans. Japan Soc. Civil Engrs.
24 (1955) 1-5.

Tyler, George W.
The Experimental Evaluation of Definite Integrals. Blacksburg: Virginia Polytechnic Institute
and State University. Ph.D. dissertation (I949).

Uspensky, J.V.
On the Expansion of the Remainder in the Newton-Cotes Formula. Trans. AMS 37 (1935)
381-396.

.‘r‘

 

 

   

111111191111111111111111I