<< . .

. 58
( : 95)



. . >>

where f is a continuous function in „¦. For this purpose, in Sections 9.9.1
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.

<< . .

. 58
( : 95)



. . >>