164
6 Computing Integrals and Testing Code
be implemented in half a page of code, but requires orders of magnitude more
function evaluations in double integrals compared to the midpoint rule.
Monte Carlo integration, however, is much more computationally efficient than
the midpoint rule when computing higher-dimensional integrals in more than three
variables over hypercube domains. Our ideas for double and triple integrals can
easily be generalized to handle an integral in m variables. A midpoint formula then
involves m sums. With n cells in each coordinate direction, the formula requires n m
function evaluations. That is, the computational work explodes as an exponential
function of the number of space dimensions. Monte Carlo integration, on the
other hand, does not suffer from this explosion of computational work and is the
preferred method for computing higher-dimensional integrals. So, it makes sense
in a chapter on numerical integration to address Monte Carlo methods, both for
handling complex domains and for handling integrals with many variables.
The Monte Carlo Integration Algorithm The idea of Monte Carlo integration
of
b
a f (x)dx is to use the mean-value theorem from calculus, which states that
the integral
b
a f (x)dx equals the length of the integration domain, here b − a,
times the average value of f , ¯
f , in [a, b]. The average value can be computed by
sampling f at a set of random points inside the domain and take the mean of the
function values. In higher dimensions, an integral is estimated as the area/volume
of the domain times the average value, and again one can evaluate the integrand
at a set of random points in the domain and compute the mean value of those
evaluations.
Let us introduce some quantities to help us make the specification of the integration algorithm more precise. Suppose we have some two-dimensional integral
Ω
f (x, y)dxdx,
where Ω is a two-dimensional domain defined via a help function g(x, y):
Ω = {(x, y) | g(x, y) ≥ 0}
That is, points (x, y) for which g(x, y) ≥ 0 lie inside Ω, and points for
which g(x, y) < Ω are outside Ω. The boundary of the domain ∂Ω is given
by the implicit curve g(x, y) = 0. Such formulations of geometries have been
very common during the last couple of decades, and one refers to g as a levelset function and the boundary g = 0 as the zero-level contour of the level-set
function. For simple geometries one can easily construct g by hand, while in more
complicated industrial applications one must resort to mathematical models for
constructing g.
Let A(Ω) be the area of a domain Ω. We can estimate the integral by this Monte
Carlo integration method:
1. embed the geometry Ω in a rectangular area R
2. draw a large number of random points (x, y) in R
6 Computing Integrals and Testing Code
be implemented in half a page of code, but requires orders of magnitude more
function evaluations in double integrals compared to the midpoint rule.
Monte Carlo integration, however, is much more computationally efficient than
the midpoint rule when computing higher-dimensional integrals in more than three
variables over hypercube domains. Our ideas for double and triple integrals can
easily be generalized to handle an integral in m variables. A midpoint formula then
involves m sums. With n cells in each coordinate direction, the formula requires n m
function evaluations. That is, the computational work explodes as an exponential
function of the number of space dimensions. Monte Carlo integration, on the
other hand, does not suffer from this explosion of computational work and is the
preferred method for computing higher-dimensional integrals. So, it makes sense
in a chapter on numerical integration to address Monte Carlo methods, both for
handling complex domains and for handling integrals with many variables.
The Monte Carlo Integration Algorithm The idea of Monte Carlo integration
of
b
a f (x)dx is to use the mean-value theorem from calculus, which states that
the integral
b
a f (x)dx equals the length of the integration domain, here b − a,
times the average value of f , ¯
f , in [a, b]. The average value can be computed by
sampling f at a set of random points inside the domain and take the mean of the
function values. In higher dimensions, an integral is estimated as the area/volume
of the domain times the average value, and again one can evaluate the integrand
at a set of random points in the domain and compute the mean value of those
evaluations.
Let us introduce some quantities to help us make the specification of the integration algorithm more precise. Suppose we have some two-dimensional integral
Ω
f (x, y)dxdx,
where Ω is a two-dimensional domain defined via a help function g(x, y):
Ω = {(x, y) | g(x, y) ≥ 0}
That is, points (x, y) for which g(x, y) ≥ 0 lie inside Ω, and points for
which g(x, y) < Ω are outside Ω. The boundary of the domain ∂Ω is given
by the implicit curve g(x, y) = 0. Such formulations of geometries have been
very common during the last couple of decades, and one refers to g as a levelset function and the boundary g = 0 as the zero-level contour of the level-set
function. For simple geometries one can easily construct g by hand, while in more
complicated industrial applications one must resort to mathematical models for
constructing g.
Let A(Ω) be the area of a domain Ω. We can estimate the integral by this Monte
Carlo integration method:
1. embed the geometry Ω in a rectangular area R
2. draw a large number of random points (x, y) in R