0% found this document useful (0 votes)
10 views15 pages

A Method For Numerical Integration: by C. B. Haselgrove

This document describes numerical integration methods based on Diophantine approximation theory. Specifically, it discusses methods to estimate multidimensional integrals using sums of the integrand evaluated at points defined by irrational numbers. It shows that certain Césaro mean-based methods can provide estimates with errors that decrease as O(N-r) for r=1,2,3,4 as the number of evaluation points N increases, assuming the integrand satisfies appropriate conditions. These methods are generally superior to Monte Carlo and direct product Gauss methods, requiring fewer function evaluations for a given error.

Uploaded by

Abid Mulla
Copyright
© © All Rights Reserved
We take content rights seriously. If you suspect this is your content, claim it here.
Available Formats
Download as PDF, TXT or read online on Scribd
0% found this document useful (0 votes)
10 views15 pages

A Method For Numerical Integration: by C. B. Haselgrove

This document describes numerical integration methods based on Diophantine approximation theory. Specifically, it discusses methods to estimate multidimensional integrals using sums of the integrand evaluated at points defined by irrational numbers. It shows that certain Césaro mean-based methods can provide estimates with errors that decrease as O(N-r) for r=1,2,3,4 as the number of evaluation points N increases, assuming the integrand satisfies appropriate conditions. These methods are generally superior to Monte Carlo and direct product Gauss methods, requiring fewer function evaluations for a given error.

Uploaded by

Abid Mulla
Copyright
© © All Rights Reserved
We take content rights seriously. If you suspect this is your content, claim it here.
Available Formats
Download as PDF, TXT or read online on Scribd
You are on page 1/ 15

A Method For Numerical Integration

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

\ g(y) dy = h^Ci g(a + «¿Ä) + R.


Ja 1-1

Received October 31, 1960.


323
324 C. B. HASELGROVE

The remainder R is 0(A2") as A varies for functions g(y) satisfying certain condi-
tions. We may estimate the integral

/ / ••• / /(x) dxi dxi ■■■dxk


Jax Ja2 Jah

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

OXl dx2 ■ ■ -oxk

but they are probably more than is necessary.


If we divide the hypercube 0 £ x¡ £ 1 into n smaller cubes with side A = 1/n,
then using the v-point Gauss formula we obtain a formula for / which requires
N = {vn)k values of /(x). This gives the value of / with an error which is
0(n-2r+1) = 0{N{-"+,)lk)

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.

(i) ¿s [k •• • j[ä/(x) asi • • • dxt = (1 - ifc)/(0, 0, • • •, 0)

+ i[f(h, 0, • • • ,0) + /( -A, 0, • • •, 0) + /(0, h, 0, • • •, 0)


+ /(0, -A,0, •••,0) + ...} +0(A4).

(H) ¿* L " ' L/(x) dXi" ' dXk= f/(0' °' ' ' ' ' 0)

+ ~ {/(A,h, ■■■, h) + /(-A, A, • • •, A) + /(A, -A, h, ■■■, A)

+ /( -A, -A, A, • • •, A) + • • • ! + 0(A4).


METHOD FOR NUMERICAL INTEGRATION 325

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. The Integration of Periodic Functions.

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

/ / ** * / /(x) dxi dx2 - • • dxk.


(2r)*
These estimates will be based on sums of the type
(1) s(N) = Yl cNmf{2irmai, 2irmai, • ■• , 2wmak),
m

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

SÁN) = (N + 1)2W + 3) {SÁ2N+l)- 2S*iN)]


and

S*(Ar)= {N + 1)4 i^(2JV) - 4S^ - DI

give substantially improved estimates of the integral. Under these conditions we


shall show that for large N,
8r(N) - I = 0{N~r)

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

hm,0s 1 sin2 UN + 1)0


kN {0) " (N + l)2 ~ sin2 he ~ '

w,as _ 1 sin2 \{N + 1)0 sin (N + f)0


(N + l)2(2iV + 3) sin2 ¿0
METHOD FOR NUMERICAL INTEGRATION 327

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 ,

where the constant c0 is chosen so that

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

(2) ¿(mMn*) • • • <t>(nk)\n« - n | á l/(* + 1)


for all sets of integers n\ , n2, ■■■ , nk, n not all zero. Further, the measure of the set
of ai, a2, ■■■ , ak with 0 á <«¡á 1 such that
(3) lbd¿(niMws) ■•■4>(nk)\n-a - n\ g S/(k + 1)
is less than 5.
We may suppose without loss of generality that 0 ^ a¿ ^ 1. We observe that
for fixed n, n the measure of the set of a such that
(4) «Kni)<¿>(n2)• • • 4>(nk)\ n-a - n \ g 8/(k + 1)
is less than or equal to
g _L_1_
(k + 1) max | nt \ <Kni)0(ii2) ' ' ' 0(n*) '
But n is a non-negative integer less than or equal to k max | n¿ | . Hence the measure
328 C. B. HASELGROVE

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

= zZ ' ' • ¿2 a„¡n2...nickN(2irn-a).


ni nk

Now ao,o,--,o = J, which is the integral which we wish to evaluate. We shall


suppose that the coefficients cNmare scaled so that kN(0) = 1. Then
(6) s(N) - I = 2Z'anin2---nkkN(2Trn-a),
where the prime ' means that the term (0, 0, • • • , 0) is omitted.
We suppose that the method is of order r so that there exists a constant K such
that
(7) \kN(6)\ ^ K\ Nsmhe\~r
for all N and 0. Now
IBin«Í | fc 2 || € ||
where || £ || denotes the distance of £ from the nearest integer. Thus if the numbers
ai satisfy the condition (2) we have

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

||(n' - n") ■« || ^ l/(t + 1)*« - V)4>(w2' - n2") ■■■<t>{nk'- nk")

^ !/(*+ lMNiMN,) ■■■4>(Nk)


so that
| || n'-a || - || n"-o || | ^ \/{k + IMN^N,) ■■■<t>(Nk).
We deduce that the number of sets of numbers n in the zone with

IIn-a || £ v/(k + lMNiMN*) ■■■4>{Nk)


330 C. B. HASELGROVE

is at most 2v + 2. Now for any set of numbers n,


|| n-a || ^ l/(k + l)0(ni)0(w2) • • • <t>(nk)
so that the contribution of that set to the sum (6) is less than or equal to
2~r(k + íyMJCMrhMn*) ■■■4>{nk)Y\ nxn2 ■■■nk\ ~sN~r.
If we assume that {<t>(n)}r\n\"s decreases monotonically for n > 0 we deduce
that the contribution of this term to the sum is less than
2~r(k + lYMJCMNiMN,) ■■■<t>(Nk)\r| N,N2 ■■■Nk\ -HT.
Thus, the sum of the contributions of all terms in the zone Ni < n,- ^ 2Ar,:to the
sum (6) is less than
2~r(k + iyMtKU(Ni)<l>(N2) ••• d>(Nk)V
(//¡^•■■fk \
2+2 g v-r\ AT'.

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.

3. Numerical Estimates for the Error.

3.1. A Method for Obtaining Numerical Estimates. In this section we shall be


especially concerned with the methods using the sums s2(N) and s4(Ar) to estimate
the integral of a periodic function
/(x) = f(xi ,x2,---,xk) = ]C • • • 2Z a„in2...„kemx.

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

(9) I an¡„2...„k I g Ms I nxn2 ■■■nk\~".

We shall- be particularly concerned with the cases s = 2 and s = 4. In order to


obtain a convenient form for the final results we shall replace zero factors in the
denominator of the right-hand side of (9) by (6/ir2) when s = 2 and by (360/77T4)
when s = 4.
METHOD FOR NUMERICAL INTEGRATION 331

We recall that the functions fcjvr(0)for r = 2 and r = 4 take the form

,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).

We now define the functions

/s(x) = M* Z ■■■E Im«* ■• • nk |-Vito*

= m*n ( e i»r«"*-')
J=l \n—oo /

where Ms* is chosen so that/„(0, 0, • • • , 0) = 1 and where, when n = 0, the fac-


tor n~a is replaced by the appropriate fraction for s = 2 or s = 4.
By summation of the trigonometric series we deduce that

/2(x) = ¿{(i -I*,!)2!

and

/4(x) = fi ¡2- (1 -|^|)2Î(1 -l^l)2


We denote the values of the sums sr(N) corresponding to the functions /»(x)
by sr¡s(N), and the mean value of /,(x) over the periodic cell by /„.
Then we deduce that

| Sr(iV) - I\â (Ma/M*)(Sr,,(N) - Is).

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

(N + l)2{s2A(N) - I<\ è 2B{1 + 0(NTus)\

and such that | a' - a | = 0(ATn/8).

4. Application to the Integration of Non-Periodic Functions. We shall now show


how the methods of Section 2 may be applied to obtain estimates of integrals of
non-periodic functions of several variables. Suppose that it is required to evaluate
the integral

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

I=¡ [ ••• [ F(xux2, ■ xk) dxi dx2 ■• ■ dxk.


Jo Jo Jo

This integral may be written in the form

, | zk | ) dzj dz2 ■■• dzk.

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

of bounded variation it may be proved that the Fourier coefficients of / satisfy,


for some M2,
anin = M21 nxn2 ■■• nk

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

S,(A0 = F(0, 0, ■■■, 0) + 2Ë F(2 HJn«,}! ,2 Ken«*}


|)
METHOD FOR NUMERICAL INTEGRATION 335

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

I = ¡ j ••• j F(zi » 22» ■■■,zk) dzi dz2 ■■■dzk.

We make the substitution a,- = Pr(xi) where r = 1 and

Pr(xi) = Ar Í' (1 - u2)r du


Jo
and Ar is chosen so that Pr(±l) = ±1. Then

I = [ • • j[ FiPrfa), ■■■,PÁXMPXXJ ■■■Pr'(xk) dxX• • ■dxk

- ôî j / - - - J /(^i >x2. ""i**) dxi dx2 •■• dxk,

say. Then / = 0 if any Xi = ±1, and F is bounded. Thus if F is continuous / may


be regarded as a continuous periodic function of the x<. If F possesses a partial
derivative
dkrF
(dXidx2 • ■• dxk)r

of bounded variation then so will /. It is then possible to obtain an estimate


Sr(N) - I = 0(ATr)

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

Si*(N) = E \f(ma) - In*


and
S2*(N) £&*(«).

Table 3
Estimates for So So • ■■ So e~z\x*"x<,dxi dx2 dxs
S2Í.N) nW Monte Carlo

1000 0.97062580 0.97062392 0.96763166


2000 0.97063927 0.97082902 0.96870265
3000 0.97066765 0.97054070 0.96885258
4000 0.97066383 0.97068153 0.96944396
5000 0.97065630 0.97065925 0.96950137
6000 0.97065761 0.97061983 0.96990269
7000 0.97065639 0.97068925 0.97018578
8000 0.97065632 0.97064881 0.97030504
9000 0.97065706 0.97063833 0.97038771
10000 0.97065854 0.97066307 0.97032729
11000 0.97065860 0.97065947 0.97029480
12000 0.97065744 0.97067426 0.97048290
Exact value = 0.97065719
METHOD FOR NUMERICAL INTEGRATION 337

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

if ••• f r"»"»•"*•dxidxi ■■■


dxk
Jo Jo Jo
was calculated on the Mercury computer. Some results for k = 5 are given in Table
3. We give the sum s2(N) for N = 1000( 1000) 12000. The method described above
was used to reduce rounding errors, with the estimate In* calculated at every 10th
value of N up to N = 100 and every 100th value thereafter. For comparison, the
sum Si(N) is also given, and an estimate using the ordinary Monte Carlo method
with N random sampling points.

6. Acknowledgments. The author is grateful for many helpful discussions with


Dr. J. C. P. Miller, Dr. D. J. Wheeler, Mr. H. P. F. Swinnerton-Dyer, and others,
and to the staff of the Manchester University Computing Machine Laboratory.
Department of Mathematics
The University
Manchester, 13, England
1. P. Davis & P. Rabinowitz, "Some Monte Carlo experiments in computing multiple
integrals," MTAC, v. 10, 1956,p. 1-8.
2. P. C. Hammer, "Numerical evaluation of multiple integrals," On Numerical Approxi-
mation, éd. R. E. Langer. Proceedings of symposium, University of Wisconsin Press, Madison,
1959, p. 99-115.
3. J. M. Hammersley & K. W. Morton, "A new Monte Carlo technique: antithetic
variâtes," Proc. Camb. Philos. Soc, v. 52, 1956, p. 449-475.
4. J. M. Hammersley & J. G. Mauldon, "General principles of antithetic variâtes,"
Proc. Camb. Philos. Soc, v. 52, 1956, p. 476-481.
5. J. C. P. Miller, "Numerical quadrature over a rectangular domain in two or more
dimensions, Part 1. Quadrature over a square, using up to sixteen equally spaced points,"
Math. Comp., v. 14, 1960, p. 13-20.
6. J. C. P. Miller, "Numerical quadrature over a rectangular domain in two or more
dimensions, Part 2. Quadrature in several dimensions, using special points," Math. Comp.,
v. 14, 1960, p. 130-138.
7. J. C. P. Miller, "Numerical quadrature over a rectangular domain in two or more
dimensions, Part 3. Quadrature of a harmonic integral," Math. Comp., v. 14, 1960, p. 240-248.

You might also like