A Method For Numerical Integration: by C. B. Haselgrove
A Method For Numerical Integration: by C. B. Haselgrove
By C. B. Haselgrove
1. Introduction. In this paper we shall give an account of some methods de-
veloped for the numerical evaluation of multidimensional integrals. These methods
are based on the theory of Diophantine approximation. They are suitable for some
problems for which the Monte Carlo method is commonly used and, like the Monte
Carlo method, are well fitted for use with an electronic digital computer. We shall
show, however, that they are superior to the Monte Carlo method provided that
the integrand satisfies certain conditions. We shall also show that they are superior,
for integrals in space of several dimensions, to formulas typified by those of Gauss
and Simpson; they may be superior even to certain new integration formulas
specially constructed for the evaluation of multiple integrals (see for example
Hammer [2], who gives a bibliography, and Miller [5], [6], [7]).
The method of antithetic variâtes which is described by Hammersley and others
[3], [4] may be used to obtain better estimates than the Monte Carlo method but
the author thinks that the method described in the present paper is simpler to apply
and gives better results.
Various authors have suggested methods which are particular cases of those
described in this paper but without the underlying theory. See for example Davis
and Rabinowitz [1].
In this section we shall give a short account of the behavior of the error in the
Monte Carlo method and the direct-product Gauss-type methods so that we can
compare these with the errors of the new methods. We shall not give an account of
the method of antithetic variâtes.
Suppose that we wish to estimate the integral
1 = / / " ' / f(Xl >Xl> " ' ' x*) dxi d%2 • • ■ dxk.
Jo Jo Jo
We shall denote the vector (xi, x2, • ■■ , xk) by x. Numerical methods for the
evaluation of I involve the calculation of /(x) at a number N of points x¿. The
most desirable of such methods for use on an electronic computer are those which
require the evaluation of/(x) at the smallest number of points x, in order to obtain
an estimate with a sufficiently small error.
The Monte Carlo method gives as an estimate for I the sum
where the points x¿ are chosen at random in the range of integration. The error
of such an estimate has standard deviation 0{N~m) provided that the function
/(x) satisfies certain conditions. It is sufficient that the function be bounded.
A Gauss-type formula for a one-dimensional integral takes the form
r>a+h v
The remainder R is 0(A2") as A varies for functions g(y) satisfying certain condi-
tions. We may estimate the integral
by applying the formula in each dimension, obtaining the direct product estimate
v v
m=i >*=i
This gives the value of the integral with an error which is 0(h +2"~1) ag h varies,
provided that the function /(x) possesses partial derivatives of order up to 2 v with
respect to each variable and that these derivatives are bounded. These conditions
are very strong in that /(x) must have the derivative
for fixed v and k. Thus if 2v — 1 < \k, v and k being fixed, the error will not de-
crease as N increases so rapidly as in the Monte Carlo method. This suggests that
v should be chosen large. However, it is not desirable to use a high-order Gauss
formula since large errors may then arise if high-order partial derivatives of the
function are large or if the function does not possess such derivatives. The high-
order partial derivatives are large if the function has singularities in the 2fc-dimen-
sional complex domain near the region of integration.
Similar arguments to those above may be applied to multidimensional integra-
tion formulas obtained as the direct product of rules such as Simpson's rule, which
gives an error that is 0(N~ilk). They may in fact be applied to any case where
the domain of integration is divided into smaller regions as above.
In this connection the following two integration formulas should be mentioned.
(H) ¿* L " ' L/(x) dXi" ' dXk= f/(0' °' ' ' ' ' 0)
Formula (i) requires (2k + 1) values of /(x) to estimate the integral over the
hypercube of side 2A, and formula (ii) requires (2 +1) values. However, if the
hypercube 0 2i Xj ¿ 1 is divided into n smaller ones of side 2A then to evaluate
/ by using (i) requires (k + \)n + 0{n'~l) values of/(x) and by using (ii) re-
quires 2n + 0(n ~ ) values. Simpson's rule used in this way requires (k + 2)n -\-
0(n _1) values. Thus, for fixed k the formulas (i) and (ii) and Simpson's rule give
errors which are 0(N~ilk), but (ii) is the most efficient.
In this paper we shall describe methods which give errors which are 0(A^_1),
0(N~2), 0(N~ ) and 0(N~ ) when applied in any number of dimensions to func-
tions which satisfy appropriate conditions. The points x¿ at which the values of
the function /(x) are required can be calculated by a computer more easily than
the random numbers required for the Monte Carlo method. No coefficients, such
as are required in the case of the Gauss-type methods, need be calculated, nor is it
necessary to decide in advance the number of points to be used.
In the following sections we first describe the methods as applied to the integra-
tion of functions periodic in each of the variables. In Section 3 we describe calcula-
tions leading to numerical estimates of some of the errors and to suitable sets of
points x, at which to calculate the function. In Section 4 we describe the practical
application of the methods to non-periodic functions, and in Section 5 we give an
example.
2.1. General Theory. Let f{x\, x2, ■■■xk) = /(x) be a periodic function with
period 2ir in each of the k variables Xi, x2, • • • , xk. We shall describe a number of
methods for estimating the integral
where the a» are certain linearly independent irrational numbers and the cNm are
certain coefficients chosen so that s(N) —>/ as N —» ». We shall be particularly
concerned with sets of coefficients cNmsuch that the number of non-zero coefficients
in the sum ( 1 ) is finite and of order N. We shall give estimates for the difference
s(N) —Im terms of N and of bounds for the derivatives of/(x).
Certain formulas which are particularly simple to use are based on sums s(JV)
which can be expressed in terms of the repeated sums
Si(N) = E /(2*ma)
and
Sr(N) = ¿ Sr-rim)
m=0
326 C. B. HASELGROVE
where r ^ 2. The Sr(N) are modified Césaro means of the sequence/(2xma), terms
with positive and negative m being taken together. We shall show that the means
•AN) = 2^ SÁN)
and
SÁN)= whySÁN)'
which are of the form ( 1 ), give good estimates of / for large N. However, although
for functions satisfying certain conditions the estimate for the error s2(N) — I is
better than the estimate for Si(iV) — /, we find that the higher-order Césaro means
do not give better estimates of /. We shall show that if /(x) satisfies certain condi-
tions then the means
for r = 1, 2, 3 and 4.
In the theory of the general method we are led to the consideration of the
function
ktf(d) = Y^cNmem\
m
in terms of which we shall give estimates of the error. If for a particular set of
coefficients cNm ,
m»i*wem
for some fixed constant K = K(r) and for all N and 0, we shall say that the method
has order r with constant K. For the sums Si(iV), s2(N), s3(Ar) and st(N) the
corresponding functions kN(6) are
u(i)/^ .-
kiP(O) 1 sin (iV + |)0
(2N + 1) sin %0
and
{«,,„x . 1 sin4 i(iV + l)0
(N + l)4 sin4 |0
Hence these sums give methods of order 1, 2, 3 and 4 respectively. We shall prove
that a method of order r gives an error which is 0{N~r) provided that the inte-
grand satisfies certain conditions.
A method may have order r for several different values of r. In particular, if
-m2/2A'l!
Cum — CoC ,
2—1C"m = 1,
m
it follows from the theory of elliptic modular functions that the method is of order
r for all positive r. However, this method has certain disadvantages in practice.
2.2. The Numbers «¿. We shall now prove a lemma which will enable us to
prove that there exist sets of irrational numbers a,- which may be used to obtain
good estimates of integrals. This lemma is a standard result in the theory of Dio-
phantine approximation but we include the proof here for completeness. The lemma
also shows that for most sets of numbers a,- the error is not substantially worse
than in the estimates given. We show in Section 3 how suitable sets of numbers
ai may be found.
We shall denote a set of integers (ni , n2, • • • , nk) by n and we shall write the
scalar product (niax + rha2 + • • • + nkak) in the form n-a.
Lemma. Let <j>{n)be any positive function such that
Z l/*(n) = 1.
—oo
Then there exist irrational numbers ai, a2, ■■■ , ak such that
of the set of a which satisfy (4) for given n and some n is less than or equal to
5/<K«i)4>(tt2)• • • <f>(nk).
Thus the measure of the set of a which satisfy (3) is less than
00 OO -i
8 V ... V_-_ = 3
-»o <t>(ni)<t>(n2) ■ ■■ <t>(nk)
(we have strict inequality since the regions overlap). Taking 5 = 1 we obtain the
first part of the lemma.
2.3. Estimates for the Error in the Integration of Periodic Functions. We now
consider the problem of estimating the integral of a periodic function of k variables
which can be represented by an absolutely convergent multiple Fourier series
/(X) = Z •••Za«,n2---n/n'X.
«I nk
We shall suppose that for some positive fixed number s there exists a constant
Ms such that if none of the n, = 0,
(5) | a„,„2...„t | ^ M, | n<n2 ■■■ nk\~".
We suppose further that if any of the nt = 0 the same inequality holds with the
zero factors omitted from the denominator on the right hand side.
Since the Fourier series is absolutely convergent we may treat the contributions
to the sum s(N) from the individual terms of the Fourier series separately. We
obtain
s(N) = £ <W(2irroa)
m
i sin
• irn ■a i - 2 1
k + 1 </>(n,)0(n2) • • • <f>(nk)'
Thus if we apply the bound for the Fourier coefficient and the condition (7) we
METHOD FOR NUMERICAL INTEGRATION 329
deduce that the contribution of the term % , n2, • • • , nk to the sum is less than
KMs2-r(k + l)r{<KrH)¿(rc2) ■■■<t>(nk)VI nxn2 ■■■nk\ ~,N~r.
We have supposed that we may drop any of the ?i¿ from the denominator if it
is zero. If we suppose that the sum
E 'í*(»)r
i«r*
converges we deduce that there exists a constant Cr,s depending only on r and s
such that
I s(N) - 11 á KMs2~\k + \yCkr,sN-r.
We have the bound
C, £ Í0(O)}r+ ¿'{*(n)}r|n|-
—00
for the coefficients Cr,„. If s > r + 1 we can choose an appropriate function <f>(n)
which leads to a numerical bound for the error of the method. However, this bound
for Cr,s is by no means the best possible and we shall show how it may be improved.
We shall now prove the theorem.
Theorem. // the numbers ax, a2, • ■• , ak are chosen appropriately and if the
function f(xi, x2, • • • , xk) is periodic with period 2x in each of the variables Xi and
its Fourier coefficients satisfy the condition (5) then if we apply any method of order
r < s, where r è 1, with constant K there exist numbers C,,e such that
(8) | s(N) -./ | S KM,2-'(k + l)TAr.s.kCkr,sN-r.
The numbers Cr,s depend only on the numbers a, and not on the particular method or
the function to be integrated.
If r = s it can be proved that | s(N) - 11 = 0(N~r+e) for « > 0, but we
shall not give the proof here.
In order to prove the theorem we suppose that the function <j>(n) is chosen to
be monotonie increasing for n ^ 0 and such that 4>(—n) = 4>{n). We divide the
range of summation of the variables nx, n2, ■• • , nk into zones defined by
Ni < Hi ^ 2Ni
(or if n¿ is negative — 2JV<^ fi; < — A", ; zones for which some of the variables are
zero are also allowed). Now if n' and n" are two sets of numbers in the same zone
we deduce that | n/ — n" \ £ AT,-. Hence
Similar results hold for those zones for which some of the n,- are negative or zero.
If s > r and r > 1 we may choose <t>(n) so that the sum over all zones con-
verges. We obtain the inequality
I s(N) - I |
S 2-\k + iyMsK{2r(r) + 2}({0(O)}r+ 2^ '{0(Ar1)}rAT,)tATr,
where Ni runs through the powers of 2 and f (r) is the Riemann zeta function.
This inequality is of the form (8) with Ar,s,k = {2f(r) + 2).
If r = 1 and s > 1 we obtain an expression with a different Cr,s and a coefficient
Ar.s.k depending on s and k.
These estimates for the errors are very crude but give some indication of the
power of the methods. In the next section we shall describe how sets of numbers
a, can be obtained and how numerical bounds for the errors in the estimates of the
integrals can be calculated.
Note that this function has period 2 in each of the variables Xi ■As before, we shall
assume that the Fourier coefficients satisfy an inequality
,mu,
Ú2)(e) 1 sm2h(N+l)d
(N + l)2 sin2 \&
and
sin4 è(A" + 1)0
^4,(0)
(N + l)4 sin410
Both these functions are positive.
Thus
| sr(N) - I | = Z •••E 'an¡ni...nkk¿r)(Trn-a)
á MsY. • • • Z I Wi«2• • • n» | ""s/cir)(7rn-o).
= m*n ( e i»r«"*-')
J=l \n—oo /
and
Thus, we may say that the functions /2(x) and Mx) are the worst functions for
numerical integration whose coefficients satisfy the inequality (9) with s = 2 and
s = 4.
If, therefore, we find sets of numbers <*¿which make (s,,,(N) — /„) for a prac-
tical range of N as small as possible, these a¿ will give good estimates for the in-
tegral of any periodic function whose Fourier coefficients satisfy (9).
3.2. The Determination of the a¿. The Ferranti Mercury computer at Man-
chester University was used to find good sets of a¿ for (r, s) = (2, 2) and (2, 4).
This was done by minimizing the upper bounds of (N + l)2[sr,a(N) — Is] for
0 < N ^ Ni. The minimization was performed by a random walk in the «-space,
a step being taken only if it decreased the upper bound. The size of the step was
332 C. B. HASELGROVE
reduced when the upper bound could not be decreased otherwise and the number
Ni was increased when the step was reduced.
For s = 2, Is = ($)* and for s = 4, 7S = (if)k. In both cases
/(0, 0, • • • , 0) = 1
and/(x) è 0 for all x, so that Si(m) ^ 1 and S2(N) äJV+1. Hence, for s = 2
we have
(N + 1)2{S2,2(N) - I2] è (N + 1) - (AT + l)2(i)fc.
For a certain value A7™of A^ the right-hand side of this inequality approaches its
maximum value of \ •3*. Therefore, for values of Ni > Nm we cannot expect the
upper bound lub (N + l)2{s2,2(N) — I2\ to be less than J-3*. The corresponding
value for s = 4 is \{xr-)k.
It was found that the points a converged to points for which the upper bound
was reasonably small. In Figure 1 the logarithms of the upper bounds for s = 2
and s — 4 are plotted against k. For s = 2, Ni was taken up to 1500; for s = 4,
up to 1000. The values j-3* for s = 2 and K^)* for s = 4 are shown on the same
logarithmic graph as straight lines. The point for s = 2, k = 8 lies below the ap-
propriate line; for this case Ni was less than the number Nm . The numbers a¿ ob-
tained for s = 2 and s — 4 are given in Tables 1 and 2 respectively. It should be
realized that although these are good sets of a¿ they are not claimed to be the best
possible.
The convergence of the points a suggests that there may be limit points a' such
that the bounds
lub (N + íyiSr.s(N) - Is]
0<JV<oo
Fig. 1.—Errors associated with the numbers on for (r, s) = (2, 2) and (2, 4). The points
are the logarithms to base 10 of the computed errors and the lines are those below which, for
large enough JV, the points cannot be expected to fall.
METHOD FOR NUMERICAL INTEGRATION 333
Table 1
Numbers ai for (r, s) = (2, 2)
k = 1 k = 7
0.73258893 0.95734608 0.80638723
k = 2 0.86730270 0.22584927
0.09724025 0.72510075
0.62055505 0.31301950 0.51310685
0.22610245 0.48476582 0.11080509
k= 3 k= 6 0.60161858
0.92715171
0.96498949 0.43657951
0.81091316 0.59185199 k= 8
0.46960090 0.05024400 0.73750248
k= 4 0.84373919 0.08314415
0.38104000 0.84753682
0.62366851 0.75808683 0.88989711
0.04150108 0.80254484
0.48574769 0.27951501
0.27210703 0.67340402
0.53040927
for these points are not substantially larger than the values given in the figure, but
the author has not been able to prove this.
In the case k = l,r = 2, s = 4, however, he has been able to prove the following
result.
//, for some a, (N + l)2{s2,i(N) — h\ for 0 g N ^ N, then there exists an
irrational a such that for all N
I = I I " ' I F(xi¡ xi, ■■■, xk) dxi dx2 ■■• dxk
where R is some region of the ft-dimensional space and F(x) is some given function.
We shall show how under certain special circumstances it is possible to transform
the problem to that of the integration of a periodic function.
We have seen that we can obtain good estimates for the error in the integration
of a periodic function if its Fourier coefficients tend to zero rapidly. Now, the
Fourier coefficients tend to zero rapidly if the function is continuous and possesses
partial derivatives of high orders. We shall, therefore, attempt to construct
"smooth" periodic functions. The methods given below for doing this are intended
only as examples. It will be apparent that each problem may require special treat-
ment to obtain the best results. Also, if the function F(x) possesses unknown dis-
334 C. B. HASELGROVE
Table 2
Numbers at for (r, s) (2,4)
k = 1 k = 5 k= 7
0.83969144 0.44810200 0.58505729
k = 2 0.53589831 0.50196855
0.56039410 0.77797734
0.59734470 0.83630131 0.60504620
0.92828094 0.22148205 0.62193588
k= 3 k= 6 0.84244165
0.74235492 0.64543976
0.10613747
0.57387033 0.40278232 k= 8
0.32279917 0.88772556 0.23975940
k = 4 0.43554826 0.01544979
0.17219381 0.57794809
0.17665781 0.63794472 0.81182909
0.71327190 0.78068912
0.98875216 0.62319488
0.60299793 0.70710061
0.60389317
continuities, it will not be possible to remove these from the periodic function
constructed.
We first consider the integral
The function f(zx, z2, • • • , zk) = F(\zx\ ,\z2\ , ■■■ ,\zk\) may be regarded as a
periodic function with period 2 in each of the variables z¿. Then / is continuous if
F is continuous, and if F possesses a derivative
dkF
dxi dx2 dxk
This implies that the means sx(N) and s2(N) will give estimates for I with errors
OiN-1) and 0(N~2+t) respectively. These methods have the particular advantage
that since/(z) = /(— z) the means Si(N) and s2(N) maybe evaluated from the
values of F at N + 1 points. The sum Si(N) may be expressed in the form
where the brackets { } denote the fractional part, lying in the range (— \, \). For
such a function the basic points of Table 1 are suitable.
We now consider the application of higher-order methods to the integral
so that the application of a method of order r will give an error of order A7"' and
the application of a method of order r + 1 will give an error of order N~r~1+t.
The numbers a¿ of Table 2 can be used to calculate the sums s2(N) when the
partial derivative d kf/{dx\dx2 ■• ■dxk) of the function has bounded variation.
It is also possible to apply these methods to the more difficult problem of the
evaluation of an integral taken over some more general regions. It is, of course,
possible that there may exist some transformation of the variables which reduces
the integral to periodic form. For example, in the evaluation of the integral of a
function taken over the interior of a circle we may take polar coordinates. The
function is then a periodic function of the angular variable. It may be made into
a periodic function of the radius with derivatives of appropriate orders by one of
the techniques above.
In the case of the integration of a function F(x) over the infinite range ( — °°, » )
we may at once consider the periodic function
*
E F(x + 2m)
—00
but this may mean that we have to take a very large number of points if F(x) does
not tend to zero rapidly asx-> ± °°.
A more general method is the method of extrusion. Suppose that a region R with
a smooth boundary is contained in a sphere S with center O in R. Then, if P is
336 c. b. haselgrove
any point in the sphere we may construct the line OP. If R is a star body about 0
this line will meet the boundary in one point Q. We may then evaluate the function
F at a point P' at a distance OP-OQ/r from O along the line OP, where r is the
radius of the sphere. If this value is multiplied by (OQ/r) we obtain a function
defined in the sphere which has the same integral over the sphere as F has over R.
This process of extrusion may be generalized to the extrusion of any star body
to any other star body containing it. It will generally be better to use a sphere
rather than a cube if R has a smooth boundary, otherwise discontinuities are in-
troduced into the derivatives.
5. Example^of Use. The method with the most practical application is probably
that using the sum s2(N). The sum st(N) may be more suitable if higher accuracy
is required but for the results to be better than those for s2(N), the function / must
have higher-order partial derivatives.
The values of s2(N) can be calculated at chosen values of N while the repeated
sum S2(N) is being accumulated. The accuracy of the results can be estimated by
examining the convergence of the successive values of s2(N).
A refinement may be introduced in order to decrease the rounding errors. The
sum S2(N) becomes very large and for a floating point computer with a fixed
number of significant figures the rounding error introduced at each addition is
proportional to the larger of the two quantities added. Thus, the sum s2(N) will
not be accurate to the full capacity of the computer. A method of overcoming this
difficulty is to calculate an estimate /„* = S2(N)/(N + l)2 at intervals in the
computation, and to subtract this value from each subsequent calculation of the
integrand, so that the accumulated sums are
Table 3
Estimates for So So • ■■ So e~z\x*"x<,dxi dx2 dxs
S2Í.N) nW Monte Carlo
The sum S2*(N) does not then become very large. At each successive re-estimate
it+i = h* + S2*(N)/(N + l)2 we replace the current S2*(N) by zero and the
current &*(#) by <Si*(A0+ (2N + 1)(J„* - ft+1). The same method used in a
fixed-point computer will reduce the amount of scaling necessary.
As an example the integral