and 9.9.2 we address two methods.

The ¬rst method applies when „¦ is a normal domain with respect to a

coordinate axis. It is based on the reduction formula for double integrals

and consists of using one-dimensional quadratures along both coordinate

direction. The second method, which applies when „¦ is a polygon, consists

of employing composite quadratures of low degree on a triangular decom-

position of the domain „¦. Section 9.9.3 brie¬‚y addresses the Monte Carlo

9.9 Multidimensional Numerical Integration 403

method, which is particularly well-suited to integration in several dimen-

sions.

9.9.1 The Method of Reduction Formula

Let „¦ be a normal domain with respect to the x axis, as drawn in Figure

9.5, and assume for the sake of simplicity that φ2 (x) > φ1 (x), ∀x ∈ [a, b].

FIGURE 9.5. Normal domain with respect to x axis

The reduction formula for double integrals gives (with obvious choice of

notation)

b φ2 (x) b

I(f ) = f (x, y)dydx = Ff (x)dx. (9.56)

a φ1 (x) a

b

The integral I(Ff ) = a Ff (x)dx can be approximated by a composite

quadrature rule using Mx subintervals {Jk , k = 1, . . . , Mx }, of width H =

(k) (k)

(b ’ a)/Mx , and in each subinterval nx + 1 nodes {xk , i = 0, . . . , nx }.

i

Thus, in the x direction we can write

(k)

Mx nx

c

±i Ff (xk ),

k

Inx (f ) = i

k=1 i=0

k

where the coe¬cients ±i are the quadrature weights on each subinterval Jk .

For each node xk , the approximate evaluation of the integral Ff (xk ) is then

i i

carried out by a composite quadrature using My subintervals {Jm , m =

1, . . . , My }, of width hk = (φ2 (xk ) ’ φ1 (xk ))/My and in each subinterval

i i i

(m) (m)

i,k

ny + 1 nodes {yj,m , j = 0, . . . , ny }.

(k) (m)

In the particular case Mx = My = M , nx = ny = 0, for k, m =

1, . . . , M , the resulting quadrature formula is the midpoint reduction for-

mula

M M

0,k

c

hk f (xk , y0,m ),

I0,0 (f ) =H 0 0

m=1

k=1

404 9. Numerical Integration

where H = (b ’ a)/M , xk = a + (k ’ 1/2)H for k = 1, . . . , M and

0

0,k

y0,m = φ1 (x0 ) + (m ’ 1/2)hk for m = 1, . . . , M . With a similar procedure

k

0

the trapezoidal reduction formula can be constructed along the coordinate

(k) (m)

directions (in that case, nx = ny = 1, for k, m = 1, . . . , M ).

The e¬ciency of the approach can obviously be increased by employing

the adaptive method described in Section 9.7.2 to suitably allocate the

i,k

quadrature nodes xk and yj,m according to the variations of f over the

i

domain „¦. The use of the reduction formulae above becomes less and less

convenient as the dimension d of the domain „¦ ‚ Rd gets larger, due to

the large increase in the computational e¬ort. Indeed, if any simple integral

requires N functional evaluations, the overall cost would be equal to N d .

The midpoint and trapezoidal reduction formulae for approximating the

integral (9.56) are implemented in Programs 78 and 79. For the sake of

simplicity, we set Mx = My = M . The variables phi1 and phi2 contain

the expressions of the functions φ1 and φ2 which delimitate the integration

domain.

Program 78 - redmidpt : Midpoint reduction formula

function inte=redmidpt(a,b,phi1,phi2,m,fun)

H=(b-a)/m; xx=[a+H/2:H:b]; dim=max(size(xx));

for i=1:dim, x=xx(i); d=eval(phi2); c=eval(phi1); h=(d-c)/m;

y=[c+h/2:h:d]; w=eval(fun); psi(i)=h*sum(w(1:m)); end;

inte=H*sum(psi(1:m));

Program 79 - redtrap : Trapezoidal reduction formula

function inte=redtrap(a,b,phi1,phi2,m,fun)

H=(b-a)/m; xx=[a:H:b]; dim=max(size(xx));

for i=1:dim, x=xx(i); d=eval(phi2); c=eval(phi1); h=(d-c)/m;

y=[c:h:d]; w=eval(fun); psi(i)=h*(0.5*w(1)+sum(w(2:m))+0.5*w(m+1));

end; inte=H*(0.5*psi(1)+sum(psi(2:m))+0.5*psi(m+1));

9.9.2 Two-Dimensional Composite Quadratures

In this section we extend to the two-dimensional case the composite inter-

polatory quadratures that have been considered in Section 9.4. We assume

that „¦ is a convex polygon on which we introduce a triangulation Th of NT

triangles or elements, such that „¦ = T , where the parameter h > 0 is

T ∈Th

the maximum edge-length in Th (see Section 8.5.2).

Exactly as happens in the one-dimensional case, interpolatory composite

quadrature rules on triangles can be devised by replacing „¦ f (x, y)dxdy

with „¦ Πk f (x, y)dxdy, where, for k ≥ 0, Πk f is the composite interpolat-

h h

ing polynomial of f on the triangulation Th introduced in Section 8.5.2.

9.9 Multidimensional Numerical Integration 405

For an e¬cient evaluation of this last integral, we employ the property of

additivity which, combined with (8.38), leads to the following interpolatory

composite rule

Πk f (x, y)dxdy = Πk f (x, y)dxdy = T

c

Ik (f ) = Ik (f )

h T

T ∈Th T T ∈Th

„¦

(9.57)

dk ’1 dk ’1

f (˜T ) lj (x, y)dxdy =

T

±j f (˜T ).

T

= zj zj

T ∈Th j=0 T ∈Th j=0

T

(j) (j)

˜

The coe¬cients ±T and the points zT are called the local weights and

nodes of the quadrature formula (9.57), respectively.

(j) ˆ

The weights ±T can be computed on the reference triangle T of vertices

(0, 0), (1, 0) and (0, 1), as follows

(j)

lj,T (x, y)dxdy = 2|T | ˆj (ˆ, y )dˆdˆ, j = 0, . . . , dk ’ 1,

±T = l xˆ xy

T ˆ

T

(0)

where |T | is the area of T . If k = 0, we get ±T = |T |, while if k = 1 we

(j)

have ±T = |T |/3, for j = 0, 1, 2.

(j) (j)

3

Denoting respectively by aT and aT = j=1 (aT )/3, for j = 1, 2, 3, the

vertices and the center of gravity of the triangle T ∈ Th , the following

formulae are obtained.

Composite midpoint formula

|T |f (aT ).

c

I0 (f ) = (9.58)

T ∈Th

Composite trapezoidal formula

3

1 (j)

|T |

c

I1 (f ) = f (aT ). (9.59)

3

T ∈Th j=1

In view of the analysis of the quadrature error Ek (f ) = I(f ) ’ Ik (f ), we

c c

introduce the following de¬nition.

De¬nition 9.1 The quadrature formula (9.57) has degree of exactness

equal to n, with n ≥ 0, if Ik (p) = T pdxdy for any p ∈ Pn (T ), where

T

Pn (T ) is de¬ned in (8.35).

The following result can be proved (see [IK66], pp. 361“362).

406 9. Numerical Integration

Property 9.4 Assume that the quadrature rule (9.57) has degree of exact-

ness on „¦ equal to n, with n ≥ 0, and that its weights are all nonnegative.

Then, there exists a positive constant Kn , independent of h, such that

|Ek (f )| ¤ Kn hn+1 |„¦|Mn+1 ,

c

for any function f ∈ C n+1 („¦), where Mn+1 is the maximum value of the

modules of the derivatives of order n + 1 of f and |„¦| is the area of „¦.

The composite formulae (9.58) and (9.59) both have degrees of exactness

equal to 1; then, due to Property 9.4, their order of in¬nitesimal with

respect to h is equal to 2.

An alternative family of quadrature rules on triangles is provided by the so-

called symmetric formulae. These are Gaussian formulae with n nodes and

high degree of exactness, and exhibit the feature that the quadrature nodes

occupy symmetric positions with respect to all corners of the reference

triangle T or, as happens for Gauss-Radau formulae, with respect to the

straight line y = x.

Considering the generic triangle T ∈ Th and denoting by aT , j = 1, 2, 3,

(j)

the midpoints of the edges of T , two examples of symmetric formulae,

having degree of exactness equal to 2 and 3, respectively, are the following

3

|T |

f (aT ),

I3 (f ) = n = 3,

(j)

3 j=1

«

3 3

|T |

3 f (aT ) + 8 f (aT ) + 27f (aT ) ,

(i)

I7 (f ) = n = 7.

(j)

60 i=1 j=1

For a description and analysis of symmetric formulae for triangles, see

[Dun85], while we refer to [Kea86] and [Dun86] for their extension to tetra-

hedra and cubes.

The composite quadrature rules (9.58) and (9.59) are implemented in

Programs 80 and 81 for the approximate evaluation of the integral of f (x, y)

over a single triangle T ∈ Th . To compute the integral over „¦ it su¬ces

to sum the result provided by the program over each triangle of Th . The

coordinates of the vertices of the triangle T are stored in the arrays xv and

yv.

Program 80 - midptr2d : Midpoint rule on a triangle

function inte=midptr2d(xv,yv,fun)

y12=yv(1)-yv(2); y23=yv(2)-yv(3); y31=yv(3)-yv(1);

areat=0.5*abs(xv(1)*y23+xv(2)*y31+xv(3)*y12);

x=sum(xv)/3; y=sum(yv)/3; inte=areat*eval(fun);

Program 81 - traptr2d : Trapezoidal rule on a triangle

9.9 Multidimensional Numerical Integration 407

function inte=traptr2d(xv,yv,fun)

y12=yv(1)-yv(2); y23=yv(2)-yv(3); y31=yv(3)-yv(1);

areat=0.5*abs(xv(1)*y23+xv(2)*y31+xv(3)*y12); inte=0;

for i=1:3, x=xv(i); y=yv(i); inte=inte+eval(fun); end;

inte=inte*areat/3;

9.9.3 Monte Carlo Methods for Numerical Integration

Numerical integration methods based on Monte Carlo techniques are a valid

tool for approximating multidimensional integrals when the space dimen-

sion of Rn gets large. These methods di¬er from the approaches considered

thus far, since the choice of quadrature nodes is done statistically accord-

ing to the values attained by random variables having a known probability

distribution.

The basic idea of the method is to interpret the integral as a statistic

mean value

f (x)dx = |„¦| |„¦|’1 χ„¦ (x)f (x)dx = |„¦|µ(f ),

Rn

„¦

where x = (x1 , x2 , . . . , xn )T and |„¦| denotes the n-dimensional volume of

„¦, χ„¦ (x) is the characteristic function of the set „¦, equal to 1 for x ∈ „¦ and

to 0 elsewhere, while µ(f ) is the mean value of the function f (X), where

X is a random variable with uniform probability density |„¦|’1 χ„¦ over Rn .

We recall that the random variable X ∈ Rn (or, more properly, random

vector) is an n-tuple of real numbers X1 (ζ), . . . , Xn (ζ) assigned to every

outcome ζ of a random experiment (see [Pap87], Chapter 4).

Having ¬xed a vector x ∈ Rn , the probability P{X ¤ x} of the random

event {X1 ¤ x1 , . . . , Xn ¤ xn } is given by

x1 xn

P{X ¤ x} = ... f (X1 , . . . , Xn )dX1 . . . dXn

’∞ ’∞

where f (X) = f (X1 , . . . , Xn ) is the probability density of the random vari-

able X ∈ Rn , such that

f (X1 , . . . , Xn ) ≥ 0, f (X1 , . . . , Xn )dX = 1.

Rn

The numerical computation of the mean value µ(f ) is carried out by taking

N independent samples x1 , . . . , xN ∈ Rn with probability density |„¦|’1 χ„¦

and evaluating their average

N

1

fN = f (xi ) = IN (f ). (9.60)

N i=1

408 9. Numerical Integration

From a statistical standpoint, the samples x1 , . . . , xN can be regarded as

the realizations of a sequence of N random variables {X1 , . . . , XN }, mu-

tually independent and each with probability density |„¦|’1 χ„¦ .

For such a sequence the strong law of large numbers ensures with prob-

N

ability 1 the convergence of the average IN (f ) = i=1 f (Xi ) /N to

the mean value µ(f ) as N ’ ∞. In computational practice the sequence

of samples x1 , . . . , xN is deterministically produced by a random-number

generator, giving rise to the so-called pseudo-random integration formulae.

The quadrature error EN (f ) = µ(f ) ’ IN (f ) as a function of N can be

characterized through the variance

2

µ (IN (f ) ’ µ(f )) .

σ(IN (f )) =

Interpreting again f as a function of the random variable X, distributed

with uniform probability density |„¦|’1 in „¦ ⊆ Rn and variance σ(f ), we

have

σ(f )

σ(IN (f )) = √ , (9.61)

N

from which, as N ’ ∞, a convergence rate of O(N ’1/2 ) follows for the

statistical estimate of the error σ(IN (f )). Such convergence rate does not

depend on the dimension n of the integration domain, and this is a most

relevant feature of the Monte Carlo method. However, it is worth noting

that the convergence rate is independent of the regularity of f ; thus, un-

like interpolatory quadratures, Monte Carlo methods do not yield more

accurate results when dealing with smooth integrands.

The estimate (9.61) is extremely weak and in practice one does often

obtain poorly accurate results. A more e¬cient implementation of Monte

Carlo methods is based on composite approach or semi-analytical methods;

an example of these techniques is provided in [ NAG95], where a composite

Monte Carlo method is employed for the computation of integrals over

hypercubes in Rn .

9.10 Applications

We consider in the next sections the computation of two integrals suggested

by applications in geometry and the mechanics of rigid bodies.

9.10.1 Computation of an Ellipsoid Surface

Let E be the ellipsoid obtained by rotating the ellipse in Figure 9.6 around

the x axis, where the radius ρ is described as a function of the axial coor-

9.10 Applications 409

E

ρ(x )

x

- 1/ β 1/ β

FIGURE 9.6. Section of the ellipsoid

dinate by the equation

1 1

ρ2 (x) = ±2 (1 ’ β 2 x2 ), ’ ¤x¤ ,

β β

± and β being given constants, assigned in such a way that ±2 β 2 < 1.

√

We set the following values for the parameters: ±2 = (3 ’ 2 2)/100 and

√

β 2 = 100. Letting K 2 = β 2 1 ’ ±2 β 2 , f (x) = 1 ’ K 2 x2 and θ =

cos’1 (K/β), the computation of the surface of E requires evaluating the

integral

1/β

2π±

[(π/2 ’ θ) + sin(2θ)/2] .

I(f ) = 4π± f (x)dx = (9.62)

K

0

Notice that f (1/β) = ’100; this prompts us to use a numerical adaptive

formula able to provide a nonuniform distribution of quadrature nodes,

with a possible re¬nement of these nodes around x = 1/β.

Table 9.12 summarizes the results obtained using the composite midpoint,

trapezoidal and Cavalieri-Simpson rules (respectively denoted by (MP),

(TR) and (CS)), which are compared with Romberg integration (RO) and

with the adaptive Cavalieri-Simpson quadrature introduced in Section 9.7.2

and denoted by (AD).

In the table, m is the number of subintervals, while Err and flops denote

the absolute quadrature error and the number of ¬‚oating-point operations

required by each algorithm, respectively. In the case of the AD method,

we have run Program 77 taking hmin=10’5 and tol=10’8 , while for the

Romberg method we have used Program 76 with n=9.

The results demonstrate the advantage of using the composite adaptive

Cavalieri-Simpson formula, both in terms of computational e¬ciency and

accuracy, as can be seen in the plots in Figure 9.7 which allow to check

the successful working of the adaptivity procedure. In Figure 9.7 (left),

we show, together with the graph of f , the nonuniform distribution of

the quadrature nodes on the x axis, while in Figure 9.7 (right) we plot

the logarithmic graph of the integration step density (piecewise constant)

∆h (x), de¬ned as the inverse of the value of the stepsize h on each active

interval A (see Section 9.7.2).

410 9. Numerical Integration

Notice the high value of ∆h at x = 1/β, where the derivative of the

integrand function is maximum.