3.4: Numerical Approximation of Multiple Integrals
As you have seen, calculating multiple integrals is tricky even for simple functions and regions. For complicated functions, it may not be possible to evaluate one of the iterated integrals in a simple closed form. Luckily there are numerical methods for approximating the value of a multiple integral. The method we will discuss is called the Monte Carlo method . The idea behind it is based on the concept of the average value of a function, which you learned in single-variable calculus. Recall that for a continuous function \(f (x)\), the average value \(\bar f \text{ of }f\) over an interval \([a,b]\) is defined as
\[\bar f = \dfrac{1}{b-a}\int_a^b f (x)\,dx \label{Eq3.11}\]
The quantity \(b − a\) is the length of the interval \([a,b]\), which can be thought of as the “volume” of the interval. Applying the same reasoning to functions of two or three variables, we define the average value of \(f (x, y)\) over a region \(R\) to be
\[\bar f = \dfrac{1}{A(R)} \iint\limits_R f (x, y)\,d A \label{Eq3.12}\]
where \(A(R)\) is the area of the region \(R\), and we define the average value of \(f (x, y, z)\) over a solid \(S\) to be
\[\bar f = \dfrac{1}{V(S)} \iiint\limits_S f (x, y, z)\,dV \label{Eq3.13}\]
where \(V(S)\) is the volume of the solid \(S\). Thus, for example, we have
\[\iint\limits_R f (x, y)\,d A = A(R) \bar f \label{Eq3.14}\]
The average value of \(f (x, y)\) over \(R\) can be thought of as representing the sum of all the values of \(f\) divided by the number of points in \(R\). Unfortunately there are an infinite number (in fact, uncountably many) points in any region, i.e. they can not be listed in a discrete sequence. But what if we took a very large number \(N\) of random points in the region \(R\) (which can be generated by a computer) and then took the average of the values of \(f\) for those points, and used that average as the value of \(\bar f\) ? This is exactly what the Monte Carlo method does. So in Formula \ref{Eq3.14} the approximation we get is
\[\iint\limits_R f (x, y)\,d A \approx A(R) \bar f \pm A(R) \sqrt{\dfrac{\bar {f^2}-(\bar f)^2}{N}}\label{Eq3.15}\]
where
\[\bar f = \dfrac{\sum_{i=1}^{N}f(x_i,y_i)}{N}\text{ and } \bar {f^2} = \dfrac{\sum_{i=1}^{N}(f(x_i,y_i))^2}{N}\label{Eq3.16}\]
with the sums taken over the \(N\) random points \((x_1 , y_1), ..., (x_N , y_N )\). The \(\pm\) “error term” in Formula \ref{Eq3.15} does not really provide hard bounds on the approximation. It represents a single standard deviation from the expected value of the integral. That is, it provides a likely bound on the error. Due to its use of random points, the Monte Carlo method is an example of a probabilistic method (as opposed to deterministic methods such as Newton’s method, which use a specific formula for generating points).
For example, we can use Formula \ref{Eq3.15} to approximate the volume \(V\) under the plane \(z = 8x + 6y\) over the rectangle \(R = [0,1] \times [0,2]\). In Example 3.1 in Section 3.1, we showed that the actual volume is 20. Below is a code listing (montecarlo.java) for a Java program that calculates the volume, using a number of points \(N\) that is passed on the command line as a parameter.
Listing 3.1 Program listing for montecarlo.java
The results of running this program with various numbers of random points (e.g. java montecarlo 100) are shown below:
As you can see, the approximation is fairly good. As \(N \to \infty\), it can be shown that the Monte Carlo approximation converges to the actual volume (on the order of \(O( \sqrt{N})\), in computational complexity terminology).
In the above example the region \(R\) was a rectangle. To use the Monte Carlo method for a nonrectangular (bounded) region \(R\), only a slight modification is needed. Pick a rectangle \(\tilde R \text{ that encloses }R\), and generate random points in that rectangle as before. Then use those points in the calculation of \(\bar f\) only if they are inside \(R\). There is no need to calculate the area of \(R\) for Equation \ref{Eq3.15} in this case, since the exclusion of points not inside \(R\) allows you to use the area of the rectangle \(\tilde R\) instead, similar to before.
For instance, in Example 3.4 we showed that the volume under the surface \(z = 8x + 6y\) over the nonrectangular region \(R = {(x, y) : 0 ≤ x ≤ 1, 0 ≤ y ≤ 2x^2 }\) is 6.4. Since the rectangle \(\tilde R = [0,1] \times [0,2] \text{ contains }R\), we can use the same program as before, with the only change being a check to see if \(y < 2x^2\) for a random point \((x, y) \text{ in }[0,1] \times [0,2]\). Listing 3.2 below contains the code (montecarlo2.java):
Listing 3.2 Program listing for montecarlo2.java
The results of running the program with various numbers of random points (e.g. java montecarlo2 1000) are shown below:
To use the Monte Carlo method to evaluate triple integrals, you will need to generate random triples \((x, y, z)\) in a parallelepiped, instead of random pairs \((x, y)\) in a rectangle, and use the volume of the parallelepiped instead of the area of a rectangle in Equation \ref{Eq3.15} (see Exercise 2). For a more detailed discussion of numerical integration methods, see PRESS et al.