Module for Monte Carlo Integration

Numerical Methods Using Mathematica 
Complementary Supplement 
for 
Numerical Analysis - Numerical Methods 
(c) 
John H. Mathews and Kurtis Fink, 2005

Background

	We will be discussing how to approximate the value of an integral based on the average of function values.  The following concept is useful.

Theorem  (Mean Value Theorem for Integrals).  If  f (x)  is continuous over  [a,b],  then there exists a number  c,  with  a < c < b,  such that  

		(1/(b  -  a))(∫_a^b) f (x) dx  =  f (c). 

This can be written in the equivalent form:

		(∫_a^b) f (x) dx  = (b - a)* f (c).  

Remark.  This computation shows that the area under the curve is the base width  (b - a)  times the "average height"  f (c).

Example 1.  Let  f (x)=sin (x)+(1/3) sin(3 x).  Find  c  so that  

		(∫_0^π)(sin (x) + 
                (1/3) sin (3 x)) dx  =  (π - 0)*f (c).

Solution 1.

f[x_] = Sin[x] + (1/3) Sin[3 x];  
g[
          x_] = ∫f[
              x] \[DifferentialD]x;  
a = 0;  
b = π;  
val = (∫_a^b) f[x] \[DifferentialD]x;  
ave = (((∫_a^b) f[x] \[DifferentialD]x)/(b - a));  
Print["     f[x]  =  ", f[x] ]; 
Print["  ∫f[x]\[DifferentialD]x  =  ", 
        g[x]];  
\(\(Print["  ", "∫"\+a\%b, "f[x]\[DifferentialD]x  =  ", "∫"\+a\%b, "(", f[x], ")\[DifferentialD]x = ", 
        val\ ];\)\)
\(\(Print[
        1\/\(b - a\), "∫"\+a\%b, "f[x]\[DifferentialD]x  =  ", "∫"\+a\%b, "(", f[x], ")\[DifferentialD]x = ", 
        ave\ ];\)\)
\(\(Print["Solve  f[c] = ", 
        1\/\(b - a\), "∫"\+a\%b, "f[x]\[DifferentialD]x", "  for c"];\)\ \ \)
Print["f[c] = ", ave];  
Clear[c];
eqn  =  f[c] == ave; 
Print[ eqn ];
solset = Solve[eqn, c];
c = 
        (Subscript[solset, [[1, 1, 2]]]);
Print["c = ", c];  
c = N[Re[c]];
Print["c = ", c];

Needs["Graphics`Colors`"];  
Needs["Graphics`FilledPlot`"];  
gr1  =  
        FilledPlot[{f[c]}, {x, a, b}, Fills -> {{{1, Axis}, Honeydew}}, 
          DisplayFunction -> Identity];  
grf = 
        Plot[f[x], {x, a, b}, 
          DisplayFunction -> Identity];  
dots = 
        Graphics[{Green, PointSize[0.025], Point[{0, f[c]}], Red, 
            PointSize[0.025], Point[{c, 0}], 
            Point[{c, f[c]}]} ]; 
line1 = 
        Graphics[{Red, Thickness[0.007], 
            Line[{{c, 0}, {c, f[c]}}]}];  
line2 = 
        Graphics[{Green, Thickness[0.008], Line[{{a, 0}, {a, f[c]}}], 
            Line[{{a, f[c]}, {b, f[c]}}], 
            Line[{{b, f[c]}, {b, 0}}]}];  
Show[gr1, grf, line1, line2, dots, 
        PlotRange -> {{0, π}, {0, 1}}, 
        Ticks -> {{0, {c, "c"}, (π/2), π}, {0.5, {f[
                  c], "f(c)"}, 1.0}}, AxesLabel -> {x, y}, 
        DisplayFunction -> $DisplayFunction];  
Print["y = f[x] = ", f[x] ];  
Print["The value of the integral is over [", a, ",", 
        b, "] is" ];  
\(\(Print["∫"\+a\%b, "f[x]\[DifferentialD]x  =  ", "∫"\+a\%b, "(", 
        f[x], ")\[DifferentialD]x"\ ];\)\)
\(\(Print["∫"\+a\%b, "f[x]\[DifferentialD]x  =  ", 
        val, "  =  ", 
        NumberForm[N[val], 17]];\)\)
Print["The average height of the function is"]; 
\(\(Print[\*""(Overscript[f, ∧])"", " = f[c] = ", 
        1\/\(b - a\), "∫"\+a\%b, "f[x]\[DifferentialD]x = ", 1\/\(b - a\), "∫"\+a\%b, "(", 
        f[x], ")\[DifferentialD]x", " = ", 
        N[ave]\ ];\)\)
Print["      c  =  ", c ];
\(\(Print[\*""\!\(f\&∧\)"", " = f[c] = f[", 
        c, "]  =  ", f[c]\ ];\)\)
Print["The area under the curve is the same as"];  
Print["the area of the rectangle of height  f[c] = ", 
        f[c] ];  
\(\(Print["∫"\+a\%b, "f[x]\[DifferentialD]x  =  ", "(b - a)*", \*""(Overscript[f, ∧])""];\)\)
\(\(Print["∫"\+a\%b, "f[x]\[DifferentialD]x  =  ", "(b - a)*f[c]"];\)\)
\(\(Print["∫"\+a\%b, "f[x]\[DifferentialD]x  =  (", 
        b, " - ", a, ")*f[", 
        c, "]"];\)\)
\(\(Print["∫"\+a\%b, "f[x]\[DifferentialD]x  =  (", 
        b - a, ")*", NumberForm[f[c], 17]];\)\)
\(\(Print["∫"\+a\%b, "f[x]\[DifferentialD]x  =  ", 
        NumberForm[\((b - a)\) f[c], 17]\ ];\)\)

Composite Midpoint Rule

An intuitive method of finding the area under a curve  y = f (x)  is to approximate that area with a series of rectangles that lie above the intervals  (Subscript[{[(Subscript[x, k - 1]),(Subscript[x, k])]}, k = 1]^n).

Theorem (Composite Midpoint Rule).  Consider  y = f (x)  over  [a,b].  Let the interval  [a,b]  be subdivided into  n  subintervals  (Subscript[{[(Subscript[x, k - 1]),(Subscript[x, k])]}, k = 1]^n)  of equal width  h = ((b  -  a)/n).  Form the equally spaced nodes  (Subscript[c, k]) = a + (k - (1/2)) h  for  k = 1,2,...,n.  The composite midpoint rule for  n  subintervals is  

		(∫_a^b) f (x) dx  ≈  h (∑_(k = 1)^n) f ((Subscript[c, k])).  

This can be written in the equivalent form  

		(∫_a^b) f (x) dx  ≈  (b - a)* (Overscript[f, ∧]),     where     (Overscript[f, ∧])    =  (1/n) (∑_(k = 1)^n) f ((Subscript[c, k])).

Module  Midpoint Rule

Corollary (Remainder term for the Midpoint Rule)  The composite midpoint rule  M(f,h)=(b - a) (1/n) (∑_(k = 1)^n) f ((Subscript[c, k]))  is an numerical approximation to the integral, and  

		(∫_a^b) f (x) dx  =  M(f,h)  +  (Subscript[E, M])(f,h).   

Furthermore, if  f (x)∈(C^2)[a,b],  then there exists a value  c  with  a < c < b  so that the error term  (Subscript[E, M])(f,h)  has the form  

		(Subscript[E, M])(f,h)  =  (((b  -  a) (f^2) (c))/24) (h^2).  

This is expressed using the "big O" notation  (Subscript[E, M])(f,h) = O((h^2)).

Module  Midpoint Rule

Animations  Midpoint Rule

Algorithm Composite Midpoint Rule.  To approximate the integral  

		(∫_a^b)f (x) dx  ≈ (b - a) (1/n) (∑_(k = 1)^n) f ((Subscript[c, k])),  

by sampling  f (x)  at the  n  equally spaced points  (Subscript[c, k]) = a + (k - (1/2)) h  for  k = 1,2,...,n,  where  h = ((b  -  a)/n).

Computer Programs  Midpoint Rule

Mathematica Subroutine (Midpoint Rule).

\(\(\(MidpointRule[
          n0_] := 
Module[{h, i}, 
n = 
            n0; 
h\  = \ \(b - a\)\/n; 
X = Table[a + \((i - 0.5)\)\ h, {i, 1, 
                n}]; \ 
f\&^\  = \ \(1\/n\) \(∑\+\(i = 1\)\%n f[X\_\(\([[\)\(i\)\(]]\)\)]\); 
approx = \((b - a)\)*f\&^; \ 
Return[{n, 
              approx}];\ ];\)\(\ \)\)\)

Example 2.  Let  f (x) = (Sqrt[x]).  Use the midpoint rule to calculate approximations to the integral  (∫_0^4) (Sqrt[x]) dx.

Solution 2.

f[x_] = (Sqrt[x]);
a = 0; 
b = 4; 
Needs["Graphics`Colors`"]; 
gticks = {Range[0, 4, 1], 
          Range[0, 2, 1]}; 
\(\(gtick[
          h_] := {gticks\_\(\([[\)\(1\)\(]]\)\), 
          Append[gticks\_\(\([[\)\(2\)\(]]\)\), {h, \*""(Overscript[f, ∧])""}]};\)\ \)
graph = 
        Plot[f[x], {x, a, b}, PlotStyle -> Magenta, 
          PlotRange -> All, 
          Ticks -> gticks];
Print["y = f[x] = ", f[x] ];

Needs["Graphics`Colors`"]; 
\(\(MidpointRuleResults[n0_] := 
        Module[{n = n0}, 
MidpointRule[
            n]; \ 
Y = f[X]; \ 
P = 
            Map[Point, Transpose[{X, Y}]]; \ 
dots = 
            Graphics[{Red, PointSize[0.01], P}\ ]; \ 
L = 
            Map[Line, 
              Transpose[{Transpose[{X, 0  Y}], 
                  Transpose[{X, Y}]}]]; \ 
lines = 
            Graphics[{Red, Thickness[0.005], 
                L}\ ]; \ 
dots1 = 
            Graphics[{Green, PointSize[0.025], 
                Point[{0, f\&^}]}\ ]; \ 
line1 = 
            Graphics[{Green, Thickness[0.007], Line[{{a, 0}, {a, f\&^}}], 
                Line[{{a, f\&^}, {b, f\&^}}], 
                Line[{{b, f\&^}, {b, 0}}]}]; \ \ 
Show[
            graph, dots, lines, dots1, line1, 
            Ticks -> 
              gtick[f\&^]]; \ 
Print["Points Generated  n  =  ", \ 
            n]; \ 
Print[\*""The average of {f[(Subscript[x, i])]"", "}"\+\("i=1"\)\%n\ , "  is"]; \ 
Print[\*""(Overscript[f, ∧])"", " = ", \ 
            1\/n, "∑"\+\("i=1"\)\%n, \*""f[(Subscript[x, i])]"", " =  ", 
            f\&^\ \ ]; \ 
Print["Approximation for the integral"]; \ 
Print["∫"\+a\%b, "(", 
            f[x], ")\[DifferentialD]x", \*""  ≈  (b-a)*((1/n)(∑_(i = 1)^n) f[(Subscript[x, i])])""]; \ 
Print["∫"\+a\%b, "(", 
            f[x], ")\[DifferentialD]x", "  ≈  (b-a)*", \ 
            1\/n, "∑"\+\("i=1"\)\%n, \*""f[(Subscript[x, i])]""]; \ 
Print["∫"\+a\%b, "(", 
            f[x], ")\[DifferentialD]x", "  ≈  (b-a)*", \*""(Overscript[f, ∧])""]; \ 
Print["∫"\+a\%b, "(", f[x], ")\[DifferentialD]x", "  ≈  (", 
            b, " - ", a, ")*(", 
            f\&^, ")"]; \ 
Print["∫"\+a\%b, "(", 
            f[x], ")\[DifferentialD]x", "  ≈  (", 
            b - a, ")*(", 
            f\&^, ")"]; \ 
Print["∫"\+a\%b, "(", 
            f[x], ")\[DifferentialD]x", "  ≈  ", 
            N[approx]\ ]; \ 
Return[{n, 
              approx}];\ \ ];\)\)

val1 = MidpointRuleResults[10];

val2 = MidpointRuleResults[100];

val3 = MidpointRuleResults[1000];

val4 = MidpointRuleResults[10000];

We are done.  

Aside.  The analytic value of the integral can be found.

f[x_] = (Sqrt[x]); 
F[
          x_]  =  ∫f[
              x] \[DifferentialD]x;  
Fa = 
        Limit[F[x], x -> a, 
          Direction -> 1]; 
Needs["Graphics`FilledPlot`"]; 
FilledPlot[f[x], {x, a, b}, PlotStyle -> Blue, 
        Fills ->  Cyan, PlotRange -> {{a, b}, {0, 2}}, 
        Ticks -> {Range[0, 4, 1], 
            Range[0, 2, 1]}]; 
Print["f[x] = ", f[x] ]; 
Print["F[x] = ∫(", 
        f[x], ")\[DifferentialD]x"];  
Print["F[x] = ", F[x] ];  
\(\(Print["∫"\+a\%b, "(", 
        f[x], ")\[DifferentialD]x = F[", b, "] - F[", 
        a, "]"\ ];\)\ \ \)
\(\(Print["∫"\+a\%b, "(", 
        f[x], ")\[DifferentialD]x = (", F[b]\ , ") - (", 
        F[a]\ , ")"\ ];\)\ \ \)
Print["Final answer"]; 
\(\(Print["∫"\+a\%b, "(", 
        f[x], ")\[DifferentialD]x = ", 
        F[b]\  - \ F[a]\ ];\)\ \ \)
\(\(Print["∫"\+a\%b, "(", 
        f[x], ")\[DifferentialD]x = ", 
        NumberForm[N[F[b]\  - \ F[a]], 
          10]\ ];\)\ \ \)
\(\(Print["The approximations obtained with the midpoint rule:"];\)\ \)
Print["Using n =", 
        PaddedForm[
          (Subscript[val1, [[1]]]), 
          5], "    area ≈ ", 
        NumberForm[
          (Subscript[val1, [[2]]]), 
          10] ]; 
Print["Using n =", 
        PaddedForm[
          (Subscript[val2, [[1]]]), 
          5], "    area ≈ ", 
        NumberForm[
          (Subscript[val2, [[2]]]), 
          10] ]; 
Print["Using n =", 
        PaddedForm[
          (Subscript[val3, [[1]]]), 
          5], "    area ≈ ", 
        NumberForm[
          (Subscript[val3, [[2]]]), 
          10] ]; 
Print["Using n =", 
        PaddedForm[
          (Subscript[val4, [[1]]]), 
          5], "    area ≈ ", 
        NumberForm[
          (Subscript[val4, [[2]]]), 
          10] ];

Monte Carlo Method

	Monte Carlo methods can be thought of as statistical simulation methods that utilize a sequences of random numbers to perform the simulation. The name "Monte Carlo'' was coined by Nicholas Constantine Metropolis (1915-1999) and inspired by Stanslaw Ulam (1909-1986), because of the similarity of statistical simulation to games of chance, and because Monte Carlo is a center for gambling and games of chance.

Approximation for an Integral

	The Monte Carlo method can be used to numerically approximate the value of an integral.  For a function of one variable the steps are:  

(i)	Pick n randomly distributed points  (Subscript[x, 1]),(Subscript[x, 2]),(Subscript[x, 3]),...,(Subscript[x, n]) in the interval  [a, b].    

(ii)	Determine the average value of the function  

		 \(f\&^\ \  = \ \ \(1\/n\) \(∑\+\(i = 1\)\%n\ f \((x\_i)\)\)\).  

(iii)	Compute the approximation to the integral  

		\(∫\_a\%b\ 
          f \((x)\)\ dx\ \  ≈ \ \ \((b - a)\)*
          f\&^\). 

(iv)	An estimate for the error is 

		Error  ≈ \(\((b - a)\) \@\(\(\(f\^2\)\&^ - \((f\&^)\)\^2\)\/n\)\),     where     \(\(f\^2\)\&^\  = \ 
        1\/n\ \(∑\+\(i = 1\)\%n\ \(f\^2\) \((x\_i)\)\)\).

Every time a Monte Carlo simulation is made using the same sample size it will come up with a slightly different value.  Larger values of  n  will produce more accurate approximations.  The values converge very slowly of the order  O((n^(-1/2))).  This property is a consequence of the Central Limit Theorem.

Proof  Monte Carlo Integration

Computer Programs  Monte Carlo Integration

Mathematica Subroutine (Monte Carlo for 1 Dimensional Integrals).

\(\(\(MonteCarloIntegral1D[
          n0_] := 
Module[{i}, 
n = 
            n0; 
X = 
            Table[Random[Real, {a, b}], {i, 1, 
                n}]; \ 
f\&^\  = \ \(1\/n\) \(∑\+\(i = 1\)\%n f[X\_\(\([[\)\(i\)\(]]\)\)]\); 
v = \((b - a)\); 
approx = 
            v*f\&^; \ 
Return[{n, 
              approx}];\ ];\)\(\ \)\)\)

The above subroutine is all we need to "do the math."  The following subroutine presents the results in a nice format.

Needs["Graphics`Colors`"]; 
\(\(MonteCarloResults1D[n0_] := 
        Module[{n = n0}, 
MonteCarloIntegral1D[
            n]; \ 
Y = f[X]; \ 
P = 
            Map[Point, Transpose[{X, Y}]]; \ 
dots = 
            Graphics[{Red, PointSize[0.01], P}\ ]; \ 
L = 
            Map[Line, 
              Transpose[{Transpose[{X, 0  Y}], 
                  Transpose[{X, Y}]}]]; \ 
lines = 
            Graphics[{Red, Thickness[0.005], 
                L}\ ]; \ 
dots1 = 
            Graphics[{Green, PointSize[0.025], 
                Point[{0, f\&^}]}\ ]; \ 
line1 = 
            Graphics[{Green, Thickness[0.007], Line[{{a, 0}, {a, f\&^}}], 
                Line[{{a, f\&^}, {b, f\&^}}], 
                Line[{{b, f\&^}, {b, 0}}]}]; \ \ 
Show[
            graph, dots, lines, dots1, line1, 
            Ticks -> 
              gtick[f\&^]]; \ 
Print["Points Generated  n  =  ", \ 
            n]; \ 
Print[\*""The average of {f[(Subscript[x, i])]"", "}"\+\("i=1"\)\%n\ , "  is"]; \ 
Print[\*""(Overscript[f, ∧])"", " = ", \ 
            1\/n, "∑"\+\("i=1"\)\%n, \*""f[(Subscript[x, i])]"", " =  ", 
            f\&^\ ]; \ 
Print["Approximation for the integral"]; \ 
Print["∫"\+a\%b, "(", 
            f[x], ")\[DifferentialD]x", \*""  ≈  (b-a)*((1/n)(∑_(i = 1)^n) f[(Subscript[x, i])])""]; \ 
Print["∫"\+a\%b, "(", 
            f[x], ")\[DifferentialD]x", "  ≈  (b-a)*", \ 
            1\/n, "∑"\+\("i=1"\)\%n, \*""f[(Subscript[x, i])]""]; \ 
Print["∫"\+a\%b, "(", 
            f[x], ")\[DifferentialD]x", "  ≈  (b-a)*", \*""(Overscript[f, ∧])""]; \ 
Print["∫"\+a\%b, "(", f[x], ")\[DifferentialD]x", "  ≈  (", 
            b, " - ", a, ")*(", 
            f\&^, ")"]; \ 
Print["∫"\+a\%b, "(", 
            f[x], ")\[DifferentialD]x", "  ≈  (", 
            b - a, ")*(", 
            f\&^, ")"]; \ 
Print["∫"\+a\%b, "(", 
            f[x], ")\[DifferentialD]x", "  ≈  ", 
            N[approx]\ ]; 
Print[""]; \ 
F2 = 
            Y\^2; \ 
\(f\^2\)\&^\  = \ \(1\/n\) 
              ReplaceAll[F2, 
                List -> Plus]; \ 
error\  = \ \((b - 
                  a)\) \@\(\(\(f\^2\)\&^ - \((f\&^)\)\^2\)\/n\); \ 
Print["The 'error estimate'  ≈  ", \ 
            error]; \ 
Print["Actual |Area-approx|  ≈  ", 
            Abs[N[∫\_a\%b f[x] \[DifferentialD]x] - 
                approx]\ ]; \ 
Return[{n, 
              approx}];\ ];\)\)

Example 3.  Let  f (x) = (Sqrt[x]).  Use the Monte Carlo method to calculate approximations to the integral  (∫_0^4) (Sqrt[x]) dx.

Solution 3.

f[x_] = (Sqrt[x]);
a = 0; 
b = 4; 
Needs["Graphics`Colors`"]; 
gticks = {Range[0, 4, 1], 
          Range[0, 2, 1]}; 
\(\(gtick[
          h_] := {gticks\_\(\([[\)\(1\)\(]]\)\), 
          Append[gticks\_\(\([[\)\(2\)\(]]\)\), {h, \*""(Overscript[f, ∧])""}]};\)\ \)
graph = 
        Plot[f[x], {x, a, b}, PlotStyle -> Magenta, 
          PlotRange -> All, 
          Ticks -> gticks];
Print["y = f[x] = ", f[x] ];

val1 = MonteCarloResults1D[10];

val2 = MonteCarloResults1D[100];

val3 = MonteCarloResults1D[1000];

val4 = MonteCarloResults1D[10000];

We are done!

Aside.  We can use Mathematica's subroutine NIntegrate to perform Monte Carlo integration, and we can verify that it is computing the same results as our subroutine MonteCarloIntegral1D.

We will need to seed the random number generator with the same integer so that we get the same sequence of random numbers.

SeedRandom[5];
approx1 = MonteCarloIntegral1D[100]; 
\(\(Print["Area  ≈  ", \ 
        approx1\_\(\([[\)\(2\)\(]]\)\)];\)\ \)

approx2 = 
        NIntegrate[f[x], {x, a, b}, Method -> MonteCarlo[5], 
          MaxPoints -> 100 ];
Print["Area  ≈  ",  approx2];

Observe that the results are the same.

We are really done!  

Aside.  The analytic value of the integral can be found.

f[x_] = (Sqrt[x]); 
F[
          x_]  =  ∫f[
              x] \[DifferentialD]x;  
Fa = 
        Limit[F[x], x -> a, 
          Direction -> 1]; 
Needs["Graphics`FilledPlot`"]; 
FilledPlot[f[x], {x, a, b}, PlotStyle -> Blue, 
        Fills ->  Cyan, PlotRange -> {{a, b}, {0, 2}}, 
        Ticks -> {Range[0, 4, 1], 
            Range[0, 2, 1]}]; 
Print["f[x] = ", f[x] ];  
Print["F[x] = ∫(", 
        f[x], ")\[DifferentialD]x"];  
Print["F[x] = ", F[x] ];  
\(\(Print["∫"\+a\%b, "(", 
        f[x], ")\[DifferentialD]x = F[", b, "] - F[", 
        a, "]"\ ];\)\ \ \)
\(\(Print["∫"\+a\%b, "(", 
        f[x], ")\[DifferentialD]x = (", F[b]\ , ") - (", 
        F[a]\ , ")"\ ];\)\ \ \)
Print["Final answer"]; 
\(\(Print["∫"\+a\%b, "(", 
        f[x], ")\[DifferentialD]x = ", 
        F[b]\  - \ F[a]\ ];\)\ \ \)
\(\(Print["∫"\+a\%b, "(", 
        f[x], ")\[DifferentialD]x = ", 
        NumberForm[N[F[b]\  - \ F[a]], 
          10]\ ];\)\ \ \)
Print["The approximations obtained with Monte Carlo \
simulation:"];
Print["Using n =", 
        PaddedForm[
          (Subscript[val1, [[1]]]), 
          5], "    area ≈ ", 
        NumberForm[
          (Subscript[val1, [[2]]]), 
          10] ]; 
Print["Using n =", 
        PaddedForm[
          (Subscript[val2, [[1]]]), 
          5], "    area ≈ ", 
        NumberForm[
          (Subscript[val2, [[2]]]), 
          10] ]; 
Print["Using n =", 
        PaddedForm[
          (Subscript[val3, [[1]]]), 
          5], "    area ≈ ", 
        NumberForm[
          (Subscript[val3, [[2]]]), 
          10] ]; 
Print["Using n =", 
        PaddedForm[
          (Subscript[val4, [[1]]]), 
          5], "    area ≈ ", 
        NumberForm[
          (Subscript[val4, [[2]]]), 
          10] ];

Example 4.  Let  f (x) = (4/(1 + (x^2))).  Use the Monte Carlo method to calculate approximations to the integral  (∫_0^1)  (4/(1 + (x^2))) dx.

Solution 4.

f[x_] = (4/(1 + (x^2)));
a = 0; 
b = 1; 
Needs["Graphics`Colors`"]; 
gticks = {Range[0, 1, 0.25], 
          Range[0, 4, 1]}; 
\(\(gtick[
          h_] := {gticks\_\(\([[\)\(1\)\(]]\)\), {1, 2, {h, \*""(Overscript[f, ∧])""}, 
            4}};\)\ \)
graph = 
        Plot[f[x], {x, a, b}, PlotRange -> {{0, 1}, {0, 4}}, 
          PlotStyle -> Magenta, 
          Ticks -> gticks];
Print["y = f[x] = ", f[x] ];

val1 = MonteCarloResults1D[10];

val2 = MonteCarloResults1D[100];

val3 = MonteCarloResults1D[1000];

val4 = MonteCarloResults1D[10000];

We are done!

Aside.  We can use Mathematica's subroutine NIntegrate to perform Monte Carlo integration, and we can verify that it is computing the same results as our subroutine MonteCarloIntegral1D.

We will need to seed the random number generator with the same integer so that we get the same sequence of random numbers.

SeedRandom[5];
approx1 = MonteCarloIntegral1D[100]; 
\(\(Print["Area  ≈  ", \ 
        approx1\_\(\([[\)\(2\)\(]]\)\)];\)\ \)

approx2 = 
        NIntegrate[f[x], {x, a, b}, Method -> MonteCarlo[5], 
          MaxPoints -> 100 ];
Print["Area  ≈  ",  approx2];

Observe that the results are the same.

We are really done!  

Aside.  The analytic value of the integral can be found.

f[x_] = (4/(1 + (x^2))); 
F[
          x_]  =  ∫f[
              x] \[DifferentialD]x; 
Needs["Graphics`FilledPlot`"];
FilledPlot[f[x], {x, a, b}, PlotStyle -> Blue, 
        Fills ->  Cyan, PlotRange -> {{a, b}, {0, 4}}, 
        Ticks -> {Range[0, 1, 0.25], 
            Range[0, 4, 1]}];
Print["f[x] = ", f[x] ]; 
Print["F[x] = ∫(", 
        f[x], ")\[DifferentialD]x"]; 
Print["F[x] = ", F[x] ]; 
\(\(Print["∫"\+a\%b, "(", 
        f[x], ")\[DifferentialD]x = F[", b, "] - F[", 
        a, "]"\ ];\)\ \)
\(\(Print["∫"\+a\%b, "(", 
        f[x], ")\[DifferentialD]x = (", F[b]\ , ") - (", 
        F[a]\ , ")"\ ];\)\ \)
Print["Final answer"];
\(\(Print["∫"\+a\%b, "(", 
        f[x], ")\[DifferentialD]x = ", 
        F[b]\  - \ F[a]\ ];\)\ \)
\(\(Print["∫"\+a\%b, "(", 
        f[x], ")\[DifferentialD]x = ", 
        NumberForm[N[F[b]\  - \ F[a]], 10]\ ];\)\ \)
Print["The approximations obtained with Monte Carlo \
simulation:"];
Print["Using n =", 
        PaddedForm[
          (Subscript[val1, [[1]]]), 
          5], "    area ≈ ", 
        NumberForm[
          (Subscript[val1, [[2]]]), 
          10] ];
Print["Using n =", 
        PaddedForm[
          (Subscript[val2, [[1]]]), 
          5], "    area ≈ ", 
        NumberForm[
          (Subscript[val2, [[2]]]), 
          10] ];
Print["Using n =", 
        PaddedForm[
          (Subscript[val3, [[1]]]), 
          5], "    area ≈ ", 
        NumberForm[
          (Subscript[val3, [[2]]]), 
          10] ];
Print["Using n =", 
        PaddedForm[
          (Subscript[val4, [[1]]]), 
          5], "    area ≈ ", 
        NumberForm[
          (Subscript[val4, [[2]]]), 
          10] ];

Example 5.  Let  f (x) = (Sqrt[x + (Sqrt[x])]).  Use the Monte Carlo method to calculate approximations to the integral  (∫_0^1) (Sqrt[x + (Sqrt[x])]) dx.

Solution 5.

f[x_] = (Sqrt[x + (Sqrt[x])]);
a = 0; 
b = 1; 
Needs["Graphics`Colors`"]; 
gticks = {Range[0, 1, 0.25], 
          Range[0, 1.25, 0.25]}; 
\(\(gtick[
          h_] := {gticks\_\(\([[\)\(1\)\(]]\)\), {0.25, 0.5, 
            0.75, {h, \*""(Overscript[f, ∧])""}, 
            1.25}};\)\ \)
graph = 
        Plot[f[x], {x, a, b}, PlotStyle -> Magenta, 
          PlotRange -> All, 
          Ticks -> gticks];
Print["y = f[x] = ", f[x] ];

val1 = MonteCarloResults1D[10];

val2 = MonteCarloResults1D[100];

val3 = MonteCarloResults1D[1000];

val4 = MonteCarloResults1D[10000];

We are done!

Aside.  We can use Mathematica's subroutine NIntegrate to perform Monte Carlo integration, and we can verify that it is computing the same results as our subroutine MonteCarloIntegral1D.

We will need to seed the random number generator with the same integer so that we get the same sequence of random numbers.

SeedRandom[5];
approx1 = MonteCarloIntegral1D[100]; 
\(\(Print["Area  ≈  ", \ 
        approx1\_\(\([[\)\(2\)\(]]\)\)];\)\ \)

approx2 = 
        NIntegrate[f[x], {x, a, b}, Method -> MonteCarlo[5], 
          MaxPoints -> 100 ];
Print["Area  ≈  ",  approx2];

Observe that the results are the same.

We are really done!  

Aside.  The analytic value of the integral can be found.

f[x_] = (Sqrt[x + (Sqrt[x])]);
F[
          x_]  =  ∫f[
              x] \[DifferentialD]x; 
Fa = Limit[F[x], x -> a];
Fa = 
        Limit[F[x], x -> a, 
          Direction -> 1]; 
Needs["Graphics`FilledPlot`"];
FilledPlot[f[x], {x, a, b}, PlotStyle -> Blue, 
        Fills ->  Cyan, PlotRange -> {{a, b}, {0, f[b]}}, 
        Ticks -> {Range[0, 1, 0.25], 
            Range[0, 2, 0.25]}];
Print["f[x] = ", f[x] ]; 
Print["F[x] = ∫(", 
        f[x], ")\[DifferentialD]x"]; 
Print["F[x] = ", F[x] ]; 
\(\(Print["∫"\+a\%b, "(", 
        f[x], ")\[DifferentialD]x = F[", b, "] - F[", 
        a, "]"\ ];\)\ \)
\(\(Print["∫"\+a\%b, "(", 
        f[x], ")\[DifferentialD]x = (", F[b]\ , ") - (", 
        Fa\ , ")"\ ];\)\ \)
Print["Final answer"];
\(\(Print["∫"\+a\%b, "(", 
        f[x], ")\[DifferentialD]x = ", 
        F[b]\  - \ Fa\ ];\)\ \)
\(\(Print["∫"\+a\%b, "(", 
        f[x], ")\[DifferentialD]x = ", 
        NumberForm[N[F[b]\  - \ Fa], 10]\ ];\)\ \)
Print["The approximations obtained with Monte Carlo \
simulation:"];
Print["Using n =", 
        PaddedForm[
          (Subscript[val1, [[1]]]), 
          5], "    area ≈ ", 
        NumberForm[
          (Subscript[val1, [[2]]]), 
          10] ];
Print["Using n =", 
        PaddedForm[
          (Subscript[val2, [[1]]]), 
          5], "    area ≈ ", 
        NumberForm[
          (Subscript[val2, [[2]]]), 
          10] ];
Print["Using n =", 
        PaddedForm[
          (Subscript[val3, [[1]]]), 
          5], "    area ≈ ", 
        NumberForm[
          (Subscript[val3, [[2]]]), 
          10] ];
Print["Using n =", 
        PaddedForm[
          (Subscript[val4, [[1]]]), 
          5], "    area ≈ ", 
        NumberForm[
          (Subscript[val4, [[2]]]), 
          10] ];

Approximation for a Double Integral

	The Monte Carlo method can be used to numerically approximate the value of a double integral.  For a function of two variables the steps are:  

(i)	Pick n randomly distributed points  ((Subscript[x, 1]),(Subscript[y, 1])),((Subscript[x, 2]),(Subscript[y, 2])),((Subscript[x, 3]),(Subscript[y, 3])),...,((Subscript[x, n]),(Subscript[y, n])) in the rectangle  [a, b]×[c, d].    

(ii)	Determine the average value of the function  

		 \(f\&^\)  =  (1/n)(∑_(i = 1)^n) f((Subscript[x, i]),(Subscript[y, i])).  

(iii)	Compute the approximation to the integral  

		\(∫\_a\%b\((\ ∫\_c\%d\ \ f \((x, 
                    y)\) dy\ )\) dx\ \  ≈ \ \ \((b - a)\)*\((d - c)\)*f\&^\). 

(iv)	An estimate for the error is

		Error  ≈ \(\((b - a)\)*\((d - 
              c)\) \@\(\(\(f\^2\)\&^ - \((f\&^)\)\^2\)\/n\)\),     where     \(\(f\^2\)\&^\) = (1/n) (∑_(i = 1)^n) (f^2)((Subscript[x, i]),(Subscript[y, i])).

Mathematica Subroutine (Monte Carlo for 2 Dimensional Integrals).

\(\(\(MonteCarloIntegral2D[
          n0_] := 
Module[{i}, 
n = 
            n0; 
X = 
            Table[{Random[Real, {a, b}], Random[Real, {c, d}]}, {i, 1, 
                n}]; \ 
f\&^\  = \ \(1\/n\) \(∑\+\(i = 1\)\%n f[X\_\(\([[\)\(i\)\(]]\)\)]\); 
v = \((b - a)\) \((d - c)\); 
approx = 
            v*f\&^; \ 
Return[{n, 
              approx}];\ ];\)\(\ \)\)\)

The above subroutine is all we need to "do the math."  The following subroutine presents the results in a nice format.

Needs["Graphics`Colors`"]; 
\(\(MonteCarloResults2D[n0_] := 
        Module[{n = n0}, 
MonteCarloIntegral2D[
            n]; \ 
Z = Map[f, X]; \ 
P = 
            Map[Point, 
              Transpose[{X\_\(\([[\)\(All, 1\)\(]]\)\), 
                  X\_\(\([[\)\(All, 2\)\(]]\)\), Z}]]; \ 
Q = 
            Map[Point, 
              Transpose[{X\_\(\([[\)\(All, 1\)\(]]\)\), 
                  X\_\(\([[\)\(All, 2\)\(]]\)\), 0*Z}]]; \ 
dots = 
            Graphics3D[{Red, PointSize[0.01], P}\ ]; \ 
L = 
            Map[Line, 
              Transpose[{Transpose[{X\_\(\([[\)\(All, 1\)\(]]\)\), 
                      X\_\(\([[\)\(All, 2\)\(]]\)\), Z}], 
                  Transpose[{X\_\(\([[\)\(All, 1\)\(]]\)\), 
                      X\_\(\([[\)\(All, 2\)\(]]\)\), 0*Z}]}]]; \ 
lines = 
            Graphics3D[{Red, Thickness[0.005], 
                L}\ ]; \ 
Show[dots, 
            PlotRange -> grange, BoxRatios -> gratio, 
            Ticks -> gticks, 
            AxesLabel -> {"x", "y", "z"}, 
            AxesEdge -> {{1, \(-1\)}, {1, \(-1\)}, {1, \(-1\)}}, 
            ViewPoint -> 
              gview]; 
Print["Points Generated  n  =  ", \ n]; \ 
Print[\*""The average of {f[(Subscript[x, i])]"", "}"\+\("i=1"\)\%n\ , "  is"]; \ 
Print[\*""(Overscript[f, ∧])"", " = ", \ 
            1\/n, "∑"\+\("i=1"\)\%n, \*""f[(Subscript[x, i])]"", " =  ", 
            f\&^\ ]; \ 
Print["Approximation for the integral"]; \ 
Print["∫"\+a\%b, "(", "∫"\+c\%d, "(", 
            f[{x, y}], ")\[DifferentialD]y", ")\[DifferentialD]x", \*""  ≈  (b-a)*(d-c)*((1/n)(∑_(i = 1)^n) f[(Subscript[x, i]),(Subscript[y, i])])""]; \ 
Print["∫"\+a\%b, "(", "∫"\+c\%d, "(", 
            f[{x, y}], ")\[DifferentialD]y", ")\[DifferentialD]x", "  ≈  (b-a)*(d-c)*", \ 
            1\/n, "∑"\+\("i=1"\)\%n, \*""f[(Subscript[x, i])]""]; \ 
Print["∫"\+a\%b, "(", "∫"\+c\%d, "(", 
            f[{x, y}], ")\[DifferentialD]y", ")\[DifferentialD]x", "  ≈  (b-a)*(d-c)*", \*""(Overscript[f, ∧])""]; \ 
Print["∫"\+a\%b, "(", "∫"\+c\%d, "(", 
            f[{x, y}], ")\[DifferentialD]y", ")\[DifferentialD]x", "  ≈  (", \((b - a)\)*\((d - c)\), ")*(", 
            f\&^, ")"]; \ 
Print["∫"\+a\%b, "(", "∫"\+c\%d, "(", 
            f[{x, y}], ")\[DifferentialD]y", ")\[DifferentialD]x", "  ≈  ", 
            N[approx]\ ]; 
Print[""]; \ 
F2 = 
            Z\^2; \ 
\(f\^2\)\&^\  = \ \(1\/n\) 
              ReplaceAll[F2, 
                List -> Plus]; \ 
error\  = \ \((b - 
                  a)\) \((d - 
                  c)\) \@\(\(\(f\^2\)\&^ - \((f\&^)\)\^2\)\/n\); \ 
Print["The   'error estimate'  ≈  ", \ 
            error]; \ 
Print["Actual |Volume-approx|  ≈  ", 
            Abs[N[∫\_a\%b\((∫\_c\%d 
                              f[{x, y}] \[DifferentialD]y)\) \[DifferentialD]x] - approx]\ ]; \ 
Return[{n, 
              approx}];\ ];\)\)

Example 6.  Let  f (x, y) = 4 - (x^2) - (y^2).  Use the Monte Carlo method to calculate approximations to the double integral  

		(∫_0^(5/4))((∫_0^(5/4)) (4 - (x^2) - 
                    (y^2)) dy) dx.

Solution 6.

f[{x_, y_}] = 4 - (x^2) - (y^2);
a = 0; 
b = (5/4);  
c = 0;  
d = (5/4);  
e = 4;
Needs["Graphics`Colors`"]; 
grange = {{a, b}, {c, d}, {0, e}};
gratio = {1, 1, 1};
gticks = {Range[a, b, 0.5], Range[c, d, 0.5], 
          Range[0, 10, 2]};
gview  = {2, 1, 1};
graph  =  
        Plot3D[f[{x, y}], {x, a, b}, {y, c, d}, PlotRange -> grange, 
          BoxRatios -> gratio, Ticks -> gticks, 
          AxesLabel -> {"x", "y", "z"}, 
          AxesEdge -> {{1, -1}, {1, -1}, {1, -1}}, 
          ViewPoint -> gview];
Print["z = f[{x,y}] = ", f[{x, y}] ];

val1 = MonteCarloResults2D[100];

val2 = MonteCarloResults2D[1000];

val3 = MonteCarloResults2D[10000];

We are done!

Aside.  We can use Mathematica's subroutine NIntegrate to perform Monte Carlo integration, and we can verify that it is computing the same results as our subroutine MonteCarloIntegral2D.

We will need to seed the random number generator with the same integer so that we get the same sequence of random numbers.

SeedRandom[5];
approx1 = MonteCarloIntegral2D[100]; 
\(\(Print["Volume  ≈  ", \ 
        approx1\_\(\([[\)\(2\)\(]]\)\)];\)\ \)

approx2 = 
        NIntegrate[Evaluate[f[{x, y}]], {x, a, b}, {y, c, d}, 
          Method -> MonteCarlo[5], 
          MaxPoints -> 100 ];
Print["Volume  ≈  ",  approx2];

Observe that the results are the same.

We are really done!  

Aside.  The analytic value of the double integral can be found.

f[{x_, y_}] = 4 - (x^2) - (y^2);
a = 0; 
b = (5/4);  
c = 0;  
d = (5/4);  
e = 4;
Needs["Graphics`Colors`"]; 
grange = {{a, b}, {c, d}, {0, e}};
gratio = {1, 1, 1};
gticks = {Range[a, b, 0.5], Range[c, d, 0.5], 
          Range[0, 10, 2]};
gview  = {2, 1, 1};
graph  =  
        Plot3D[f[{x, y}], {x, a, b}, {y, c, d}, PlotRange -> grange, 
          BoxRatios -> gratio, Ticks -> gticks, 
          AxesLabel -> {"x", "y", "z"}, 
          AxesEdge -> {{1, -1}, {1, -1}, {1, -1}}, 
          ViewPoint -> gview];
Print["     z = f[{x,y}] = ", 
        f[{x, y}] ];
\(\(Print["     ", "∫"\+a\%b, "(", "∫"\+c\%d, "(", 
        f[{x, y}], ")\[DifferentialD]y)\[DifferentialD]x = ", "∫"\+a\%b, "(", " (", ∫f[{x, 
                y}] \[DifferentialD]y, ")", "|"\+\(y == c\)\%\(y == d\), " )\[DifferentialD]x"\ ];\)\ \)
\(\(Print["     ", "∫"\+a\%b, "(", "∫"\+c\%d, "(", 
        f[{x, y}], ")\[DifferentialD]y)\[DifferentialD]x = ", "∫"\+a\%b, "(", ∫\_c\%d 
              f[{x, y}] \[DifferentialD]y, ")\[DifferentialD]x"\ ];\)\ \)
\(\(Print["     ", "∫"\+a\%b, "(", "∫"\+c\%d, "(", 
        f[{x, y}], ")\[DifferentialD]y)\[DifferentialD]x = ", "(", ∫\((∫\_c\%d 
                    f[{x, y}] \[DifferentialD]y)\) \[DifferentialD]x, ")", "|"\+\(x == a\)\%\(x == b\)\ ];\)\ \)
\(\(Print["     ", "∫"\+a\%b, "(", "∫"\+c\%d, "(", 
        f[{x, y}], ")\[DifferentialD]y)\[DifferentialD]x = ", \ ∫\_a\%b\((∫\_c\%d 
                    f[{x, y}] \[DifferentialD]y)\) \[DifferentialD]x\ ];\)\ \)
\(\(Print["     ", "∫"\+a\%b, "(", "∫"\+c\%d, "(", 
        f[{x, y}], ")\[DifferentialD]y)\[DifferentialD]x = ", \ 
        NumberForm[
          N[∫\_a\%b\((∫\_c\%d 
                        f[{x, y}] \[DifferentialD]y)\) \[DifferentialD]x\ ], 
          10]\ ];\)\ \)
Print["The approximations obtained with Monte Carlo \
simulation:"];
Print["Using n =", 
        PaddedForm[
          (Subscript[val1, [[1]]]), 
          5], "    volume ≈ ", 
        NumberForm[
          (Subscript[val1, [[2]]]), 
          10] ];
Print["Using n =", 
        PaddedForm[
          (Subscript[val2, [[1]]]), 
          5], "    volume ≈ ", 
        NumberForm[
          (Subscript[val2, [[2]]]), 
          10] ];
Print["Using n =", 
        PaddedForm[
          (Subscript[val3, [[1]]]), 
          5], "    volume ≈ ", 
        NumberForm[
          (Subscript[val3, [[2]]]), 
          10] ];

Example 7.  Let  f (x, y) = (Sqrt[4 - (x^2) - (y^2)]).  Use the Monte Carlo method to calculate approximations to the double integral  

		(∫_0^(5/4))((∫_0^(5/4)) (Sqrt[4 - (x^2) - 
                      (y^2)]) dy) dx.

Solution 7.

f[{x_, y_}] = (Sqrt[4 - (x^2) - (y^2)]);
a = 0; 
b = (5/4);  
c = 0;  
d = (5/4);  
e = 2;
Needs["Graphics`Colors`"]; 
grange = {{a, b}, {c, d}, {0, e}};
gratio = {1, 1, 1};
gticks = {Range[a, b, 0.5], Range[c, d, 0.5], 
          Range[0, 3, 1]};
gview  = {2, 1, 1};
graph  =  
        Plot3D[f[{x, y}], {x, a, b}, {y, c, d}, PlotRange -> grange, 
          BoxRatios -> gratio, Ticks -> gticks, 
          AxesLabel -> {"x", "y", "z"}, 
          AxesEdge -> {{1, -1}, {1, -1}, {1, -1}}, 
          ViewPoint -> gview];
Print["z = f[{x,y}] = ", f[{x, y}] ];

val1 = MonteCarloResults2D[100];

val2 = MonteCarloResults2D[1000];

val3 = MonteCarloResults2D[10000];

We are done!

Aside.  We can use Mathematica's subroutine NIntegrate to perform Monte Carlo integration, and we can verify that it is computing the same results as our subroutine MonteCarloIntegral2D.

We will need to seed the random number generator with the same integer so that we get the same sequence of random numbers.

SeedRandom[5];
approx1 = MonteCarloIntegral2D[100]; 
\(\(Print["Volume  ≈  ", \ 
        approx1\_\(\([[\)\(2\)\(]]\)\)];\)\ \)

approx2 = 
        NIntegrate[Evaluate[f[{x, y}]], {x, a, b}, {y, c, d}, 
          Method -> MonteCarlo[5], 
          MaxPoints -> 100 ];
Print["Volume  ≈  ",  approx2];

Observe that the results are the same.

We are really done!  

Aside.  The analytic value of the double integral can be found.

f[{x_, y_}] = (Sqrt[4 - (x^2) - (y^2)]);
a = 0; 
b = (5/4);  
c = 0;  
d = (5/4);  
e = 2;
Needs["Graphics`Colors`"]; 
grange = {{a, b}, {c, d}, {0, e}};
gratio = {1, 1, 1};
gticks = {Range[a, b, 0.5], Range[c, d, 0.5], 
          Range[0, 10, 2]};
gview  = {2, 1, 1};
graph  =  
        Plot3D[f[{x, y}], {x, a, b}, {y, c, d}, PlotRange -> grange, 
          BoxRatios -> gratio, Ticks -> gticks, 
          AxesLabel -> {"x", "y", "z"}, 
          AxesEdge -> {{1, -1}, {1, -1}, {1, -1}}, 
          ViewPoint -> gview];
Print["z = f[{x,y}] = ", 
        f[{x, y}] ];
\(\(Print["∫"\+a\%b, "(", "∫"\+c\%d, "(", f[{x, 
            y}], ")\[DifferentialD]y)\[DifferentialD]x  = ", "∫"\+a\%b, "(", " (", ∫f[{x, 
                y}] \[DifferentialD]y, ")", "|"\+\(y == c\)\%\(y == d\), " )\[DifferentialD]x"\ ];\)\ \)
\(\(Print["∫"\+a\%b, "(", "∫"\+c\%d, "(", f[{x, 
            y}], ")\[DifferentialD]y)\[DifferentialD]x  = ", "∫"\+a\%b, "(", ∫\_c\%d 
              f[{x, y}] \[DifferentialD]y, ")\[DifferentialD]x"\ ];\)\ \)
\(\(Print["∫"\+a\%b, "(", "∫"\+c\%d, "(", f[{x, 
            y}], ")\[DifferentialD]y)\[DifferentialD]x  = ", "(", ∫\((∫\_c\%d 
                    f[{x, 
                      y}] \[DifferentialD]y)\) \[DifferentialD]x, ")", "|"\+\(x == a\)\%\(x == b\)\ ];\)\ \)
\(\(Print["∫"\+a\%b, "(", "∫"\+c\%d, "(", f[{x, 
            y}], ")\[DifferentialD]y)\[DifferentialD]x  = ", \ ∫\_a\%b\((∫\_c\%d 
                    f[{x, y}] \[DifferentialD]y)\) \[DifferentialD]x\ ];\)\ \)
\(\(Print["∫"\+a\%b, "(", "∫"\+c\%d, "(", f[{x, y}], ")\[DifferentialD]y)\[DifferentialD]x  =  ", \ 
        NumberForm[
          N[∫\_a\%b\((∫\_c\%d 
                        f[{x, y}] \[DifferentialD]y)\) \[DifferentialD]x\ ], 
          10]\ ];\)\ \)
Print["The approximations obtained with Monte Carlo \
simulation:"];
Print["Using n =", 
        PaddedForm[
          (Subscript[val1, [[1]]]), 
          5], "   volume ≈ ", 
        NumberForm[
          (Subscript[val1, [[2]]]), 
          10] ];
Print["Using n =", 
        PaddedForm[
          (Subscript[val2, [[1]]]), 
          5], "   volume ≈ ", 
        NumberForm[
          (Subscript[val2, [[2]]]), 
          10] ];
Print["Using n =", 
        PaddedForm[
          (Subscript[val3, [[1]]]), 
          5], "   volume ≈ ", 
        NumberForm[
          (Subscript[val3, [[2]]]), 
          10] ];

Iterated Integrals in Higher Dimensions

	Sometimes we are given integrals which cannot be done analytically, especially in higher dimensions where the standard  methods of discretization can become computationally expensive.  For example, the error in the composite midpoint rule (and the composite trapezoidal rule) of an d-dimensional integral has the order of convergence   O((n^(-2/d))).  We can apply the inequality  

		O((n^(-1/2)))< O((n^(-2/d)))  when  d > 4  

to see that  Monte-Carlo integration will usually converge faster for quintuple multiple integrals and higher, i.e. \(\(∫\(∫\(∫\(∫∫\)\)\)\)\+V\ \ f\ \ dw\ du\ dz\ dy\ dx\), etc.

Approximation for a Multiple Integral

	The Monte Carlo method can be used to numerically approximate the value of a multiple integrals.  For a function of d variables the steps are:  

(i)	Pick n randomly distributed points  (Subscript[((Overscript[X,  ⇀]) ), 1]),(Subscript[((Overscript[X,  ⇀]) ), 2]),...,(Subscript[((Overscript[X,  ⇀]) ), n])  in the "volume"   V = [(Subscript[a, 1]), (Subscript[b, 1])]×[(Subscript[a, 2]), 
                (Subscript[b, 2])]× ... ×[(Subscript[a, d]), (Subscript[b, d])].    

(ii)	Determine the average value of the function  

		 \(f\&^\ \  = \ \ \(1\/n\) \(∑\+\(i = 1\)\%n\ 
            f \((\(X\& ⇀ \)\_i)\)\)\).  

(iii)	Compute the approximation to the integral  

		\(\(\(\(∫∫ ... \) ∫\)\+V\ \ f \((\(X\& ⇀ \))\)\ \ dV\ \  ≈ \ \ \((b\_1 - 
                  a\_1)\) \((b\_2 - a\_2)\) ... \) \((b\_d - a\_d)\)*
        f\&^\). 

(iv)	An estimate for the error is 

		Error  ≈ ((Subscript[b, 1]) - (Subscript[a, 1])) ((Subscript[b, 2]) - (Subscript[a, 2]))...((Subscript[b, d]) - 
            (Subscript[a, d]))\(\@\(\(\(f\^2\)\&^ - \((f\&^)\)\^2\)\/n\)\),     where     \(\(f\^2\)\&^\  = \ 
        1\/n\ \(∑\+\(i = 1\)\%n\ \(f\^2\) \((\(X\& ⇀ \)\_i)\)\)\).

Mathematica Subroutine (Monte Carlo for 3 Dimensional Integrals).

\(\(\(MonteCarloIntegral3D[
          n0_] := 
Module[{i}, 
n = 
            n0; 
X = 
            Table[{Random[Real, {a1, b1}], Random[Real, {a2, b2}], 
                Random[Real, {a3, b3}]}, {i, 1, 
                n}]; \ 
f\&^\  = \ \(1\/n\) \(∑\+\(i = 1\)\%n f[X\_\(\([[\)\(i\)\(]]\)\)]\); 
v = \((b1 - a1)\) \((b2 - a2)\) \((b3 - 
                  a3)\); 
approx = 
            v*f\&^; \ 
Return[{n, 
              approx}];\ ];\)\(\ \)\)\)

Example 8.  Let  f (x, y, z)  =  4 - (x^2) - (y^2) - (z^2).  Use the Monte Carlo method to calculate approximations to the triple integral  

		\(∫\_0\%\(9/10\)\((∫\_0\%1\((∫\_0\%\(11/10\)\ \ \((4 - x\^2 - y\^2 - 
                          z\^2)\) dz)\) dy)\) dx\).

Solution 8.

f[{x_, y_, z_}]  =  
        4 - (x^2) - (y^2) - (z^2); 
\(\(int = \*""(∫_0^(9/10))((∫_0^1)((∫_0^(11/10)) (4-(x^2)-(y^2)-(z^2))\[DifferentialD]z)\[DifferentialD]y)\[DifferentialD]x"";\)\ \)
val  =  (∫_0^(9/10))((∫_0^1)((∫_0^(11/10)) 
                      f[{x, y, 
                            z}] \[DifferentialD]z) \[DifferentialD]y) \[DifferentialD]x; 
app  =  
        NIntegrate[
          f[{x, y, z}], {x, 0, (9/10)}, {y, 0, 1}, {z, 0, 
            1  (11/10)}];  
Print["f[{x,y,z}] = ", 
        f[{x, y, z}] ]; 
Print[int, " = ", val, " = ",  
        NumberForm[N[val], 10] ]; 
Print["Mathematica's approximation using 'NIntegrate' is"];
Print[int, " ≈ ", 
        NumberForm[app, 10] ];

a1 = a2 = a3 = 0; 
b1 = (9/10);
b2 = 1;
b3 = (11/10); 
val1 = MonteCarloIntegral3D[100]; 
val2 = MonteCarloIntegral3D[1000]; 
val3 = MonteCarloIntegral3D[10000]; 
Print["f[{x,y,z}] = ", 
        f[{x, y, z}] ];
Print[int, " = ", val ];
Print["Mathematica's approximation using 'NIntegrate' is"];
Print[int, " ≈ ", 
        NumberForm[app, 10] ];
Print["The approximations obtained with Monte Carlo \
simulation:"];
Print["Using n =", 
        PaddedForm[
          (Subscript[val1, [[1]]]), 
          5], "    integral ≈ ", 
        NumberForm[
          (Subscript[val1, [[2]]]), 
          10] ];
Print["Using n =", 
        PaddedForm[
          (Subscript[val2, [[1]]]), 
          5], "    integral ≈ ", 
        NumberForm[
          (Subscript[val2, [[2]]]), 
          10] ];
Print["Using n =", 
        PaddedForm[
          (Subscript[val3, [[1]]]), 
          5], "    integral ≈ ", 
        NumberForm[
          (Subscript[val3, [[2]]]), 
          10] ];

We are done!

Aside.  We can use Mathematica's subroutine NIntegrate to perform Monte Carlo integration, and we can verify that it is computing the same results as our subroutine MonteCarloIntegral3D.

SeedRandom[5];
approx1 = MonteCarloIntegral3D[100]; 
\(\(Print["Integral  ≈  ", \ 
        approx1\_\(\([[\)\(2\)\(]]\)\)];\)\ \)

approx2  =  
        NIntegrate[
          Evaluate[f[{x, y, z}]], {x, 0, (9/10)}, {y, 0, 1}, {z, 0, 
            1  (11/10)}, Method -> MonteCarlo[5], 
          MaxPoints -> 100 ];
Print["Integral  ≈  ",  approx2];

Observe that the results are the same.

We are really done!  

Aside.  The analytic value of the integral can be found and involves the following steps.

∫ (4 - (x^2) - (y^2) - (z^2)) \[DifferentialD]z

(∫_0^(11/10)) (4 - (x^2) - (y^2) - 
            (z^2)) \[DifferentialD]z

∫((∫_0^(11/10)) (4 - (x^2) - (y^2) - 
                  (z^2)) \[DifferentialD]z) \[DifferentialD]y

(∫_0^1)((∫_0^(11/10)) (4 - (x^2) - (y^2) - 
                  (z^2)) \[DifferentialD]z) \[DifferentialD]y

∫((∫_0^1)((∫_0^(11/10)) (4 - 
                        (x^2) - (y^2) - 
                        (z^2)) \[DifferentialD]z) \[DifferentialD]y) \[DifferentialD]x

(∫_0^(9/10))((∫_0^1)((∫_0^(11/10)) (4 - (x^2) - (y^2) - 
                        (z^2)) \[DifferentialD]z) \[DifferentialD]y) \[DifferentialD]x

Mathematica Subroutine (Monte Carlo for 4 Dimensional Integrals).

\(\(\(MonteCarloIntegral4D[
          n0_] := 
Module[{i}, 
n = 
            n0; 
X = 
            Table[{Random[Real, {a1, b1}], Random[Real, {a2, b2}], 
                Random[Real, {a3, b3}], Random[Real, {a4, b4}]}, {i, 1, 
                n}]; \ 
f\&^\  = \ \(1\/n\) \(∑\+\(i = 1\)\%n f[X\_\(\([[\)\(i\)\(]]\)\)]\); 
v = \((b1 - a1)\) \((b2 - a2)\) \((b3 - a3)\) \((b4 - 
                  a4)\); 
approx = 
            v*f\&^; \ 
Return[{n, 
              approx}];\ ];\)\(\ \)\)\)

Example 9.  Let  f (x, y, z, u)  =  5 - (x^2) - (y^2) - (z^2) - (u^2).  Use the Monte Carlo method to calculate approximations to the integral  

		\(∫\_0\%\(4/5\)\((∫\_0\%\(9/10\)\((∫\_0\%1\((∫\_0\%\(11/10\)\ \ \((5 - x\^2 - y\^2 - z\^2 - 
                                u\^2\ )\) du)\) dz)\) dy)\) dx\).

Solution 9.

f[{x_, y_, z_, u_}] = 
        5 - (x^2) - (y^2) - (z^2) - (u^2); 
\(\(int = \*""(∫_0^(4/5))((∫_0^(9/10))((∫_0^1)((∫_0^(11/10)) (5-(x^2)-(y^2)-(z^2)-(u^2))\[DifferentialD]u)\[DifferentialD]z)\[DifferentialD]y)\[DifferentialD]x"";\)\ \)
val  =  (∫_0^(4/5))((∫_0^(9/10))((∫_0^1)((∫_0^(11/10)) 
                            f[{x, y, z, 
                                  u}] \[DifferentialD]u) \[DifferentialD]z) \[DifferentialD]y) \[DifferentialD]x; 
app  =  
        NIntegrate[
          f[{x, y, z, u}], {x, 0, (4/5)}, {y, 0, (9/10)}, {z, 0, 1}, {u, 0, 
            (11/10)}];  
Print["f[{x,y,z,u}] = ", 
        f[{x, y, z, u}] ]; 
Print[int, " = ", val, " = ",  
        NumberForm[N[val], 10] ]; 
Print["Mathematica's approximation using 'NIntegrate' is"];
Print[int, " ≈ ", 
        NumberForm[app, 10] ];

a1 = a2 = a3 = a4 = 0; 
b1 = (4/5);
b2 = (9/10);
b3 = 1;
b4 = (11/10); 
val1 = MonteCarloIntegral4D[100];  
val2 = MonteCarloIntegral4D[1000];  
val3 = MonteCarloIntegral4D[10000];  
Print["f[{x,y,z,u}] = ", 
        f[{x, y, z, u}] ];  
Print[int, " = ", val ]; 
Print["Mathematica's approximation using 'NIntegrate' is"];
Print[int, " ≈ ", 
        NumberForm[app, 10] ];
Print["The approximations obtained with Monte Carlo \
simulation:"]; 
Print["Using n =", 
        PaddedForm[
          (Subscript[val1, [[1]]]), 
          5], "    integral ≈ ", 
        NumberForm[
          (Subscript[val1, [[2]]]), 
          10] ];
Print["Using n =", 
        PaddedForm[
          (Subscript[val2, [[1]]]), 
          5], "    integral ≈ ", 
        NumberForm[
          (Subscript[val2, [[2]]]), 
          10] ];
Print["Using n =", 
        PaddedForm[
          (Subscript[val3, [[1]]]), 
          5], "    integral ≈ ", 
        NumberForm[
          (Subscript[val3, [[2]]]), 
          10] ];

We are done!

Aside.  We can use Mathematica's subroutine NIntegrate to perform Monte Carlo integration, and we can verify that it is computing the same results as our subroutine MonteCarloIntegral4D.

SeedRandom[5];
approx1 = MonteCarloIntegral4D[100]; 
\(\(Print["Integral  ≈  ", \ 
        approx1\_\(\([[\)\(2\)\(]]\)\)];\)\ \)

approx2  =  
        NIntegrate[
          Evaluate[f[{x, y, z, u}]], {x, 0, (4/5)}, {y, 0, (9/10)}, {z, 0, 
            1}, {u, 0, (11/10)}, Method -> MonteCarlo[5], 
          MaxPoints -> 100 ];
Print["Integral  ≈  ",  approx2];

Mathematica Subroutine (Monte Carlo for 5 Dimensional Integrals).

\(\(\(MonteCarloIntegral5D[
          n0_] := 
Module[{i}, 
n = 
            n0; 
X = 
            Table[{Random[Real, {a1, b1}], Random[Real, {a2, b2}], 
                Random[Real, {a3, b3}], Random[Real, {a4, b4}], 
                Random[Real, {a5, b5}]}, {i, 1, 
                n}]; \ 
f\&^\  = \ \(1\/n\) \(∑\+\(i = 1\)\%n f[X\_\(\([[\)\(i\)\(]]\)\)]\); 
v = \((b1 - a1)\) \((b2 - a2)\) \((b3 - a3)\) \((b4 - 
                  a4)\) \((b5 - a5)\); 
approx = 
            v*f\&^; \ 
Return[{n, 
              approx}];\ ];\)\(\ \)\)\)

Example 10.  Let  f (x, y, z, u, w)  =  6 - (x^2) - (y^2) - (z^2) - (u^2) - (w^2).  Use the Monte Carlo method to calculate approximations to the integral  

		\(∫\_0\%\(7/10\)\((∫\_0\%\(4/5\)\((∫\_0\%\(9/10\)\((∫\_0\%1\((∫\_0\%\(11/10\)\ \ \((6 - x\^2 - 
                                      y\^2 - z\^2 - u\^2 - 
                                      w\^2)\) dw)\) du)\) dz)\) dy)\) dx\).

Solution 10.

f[{x_, y_, z_, u_, w_}] = 
        6 - (x^2) - (y^2) - (z^2) - (u^2) - (w^2); 
\(\(int = \*""(∫_0^1)((∫_0^1)((∫_0^1)((∫_0^1)((∫_0^1) (6-(x^2)-(y^2)-(z^2)-(u^2)-(w^2))\[DifferentialD]w)\[DifferentialD]u)\[DifferentialD]z)\[DifferentialD]y)\[DifferentialD]x"";\)\ \)
val  =  (∫_0^(7/10))((∫_0^(4/5))((∫_0^(9/10))((∫_0^1)((∫_0^(11/10)) 
                                  f[{x, y, z, u, 
                                        w}] \[DifferentialD]w) \[DifferentialD]u) \[DifferentialD]z) \[DifferentialD]y) \[DifferentialD]x; 
app  =  
        NIntegrate[
          f[{x, y, z, u, w}], {x, 0, (7/10)}, {y, 0, (4/5)}, {z, 0, (9/10)}, {u, 
            0, 1}, {w, 0, (11/10)}];  
Print["f[{x,y,z,u,w}] = ", 
        f[{x, y, z, u, w}] ]; 
Print[int, " = ", val, " = ",  
        NumberForm[N[val], 10] ]; 
Print["Mathematica's approximation using 'NIntegrate' is"];
Print[int, " ≈ ", 
        NumberForm[app, 10] ];

a1 = a2 = a3 = a4 = a5 = 
                0;
b1 = (7/10);
b2 = (4/5);
b3 = (9/10);
b4 = 1;
b5 = (11/10); 
val1 = MonteCarloIntegral5D[100]; 
val2 = MonteCarloIntegral5D[1000];  
val3 = MonteCarloIntegral5D[10000]; 
Print["f[{x,y,z,u,w}] = ", 
        f[{x, y, z, u, w}] ]; 
Print[int, " = ", val ];
Print["Mathematica's approximation using 'NIntegrate' is"];
Print[int, " ≈ ", 
        NumberForm[app, 10] ];
Print["The approximations obtained with Monte Carlo \
simulation:"];
Print["Using n =", 
        PaddedForm[
          (Subscript[val1, [[1]]]), 
          5], "    integral ≈ ", 
        NumberForm[
          (Subscript[val1, [[2]]]), 
          10] ];
Print["Using n =", 
        PaddedForm[
          (Subscript[val2, [[1]]]), 
          5], "    integral ≈ ", 
        NumberForm[
          (Subscript[val2, [[2]]]), 
          10] ];
Print["Using n =", 
        PaddedForm[
          (Subscript[val3, [[1]]]), 
          5], "    integral ≈ ", 
        NumberForm[
          (Subscript[val3, [[2]]]), 
          10] ];

We are done!

Aside.  We can use Mathematica's subroutine NIntegrate to perform Monte Carlo integration, and we can verify that it is computing the same results as our subroutine MonteCarloIntegral5D.

SeedRandom[5];
approx1 = MonteCarloIntegral5D[100]; 
\(\(Print["Integral  ≈  ", \ 
        approx1\_\(\([[\)\(2\)\(]]\)\)];\)\ \)

approx2  =  
        NIntegrate[
          Evaluate[f[{x, y, z, u, w}]], {x, 0, (7/10)}, {y, 0, (4/5)}, {z, 0, 
            (9/10)}, {u, 0, 1}, {w, 0, (11/10)}, Method -> MonteCarlo[5], 
          MaxPoints -> 100 ];
Print["Integral  ≈  ",  approx2];

Example 11.  Let  f (x, y, z)  =  (Sqrt[4 - (x^2) - (y^2) - (z^2)]).  Use the Monte Carlo method to calculate approximations to the integral  

		\(∫\_0\%\(9/10\)\((∫\_0\%1\((∫\_0\%\(11/10\)\ \ \(\@\(4 - x\^2 - y\^2 - 
                            z\^2\)\) dz)\) dy)\) dx\).

Solution 11.

f[{x_, y_, 
            z_}]  =  (Sqrt[4 - (x^2) - (y^2) - (z^2)]); 
\(\(int = \*""(∫_0^(9/10))((∫_0^1)((∫_0^(11/10)) (Sqrt[4 - (x^2) - (y^2) - (z^2)]) \[DifferentialD]z)\[DifferentialD]y)\[DifferentialD]x"";\)\ \)
app  =  
        NIntegrate[
          f[{x, y, z}], {x, 0, (9/10)}, {y, 0, 1}, {z, 0, 
            1  (11/10)}];  
Print["f[{x,y,z}] = ", 
        f[{x, y, z}] ]; 
Print["Mathematica's approximation using 'NIntegrate' is"];
Print[int, " ≈ ", 
        NumberForm[app, 10] ];

a1 = a2 = a3 = 0; 
b1 = (9/10);
b2 = 1;
b3 = (11/10); 
val1 = MonteCarloIntegral3D[100]; 
val2 = MonteCarloIntegral3D[1000]; 
val3 = MonteCarloIntegral3D[10000]; 
Print["f[{x,y,z}] = ", 
        f[{x, y, z}] ];
Print["Mathematica's approximation using 'NIntegrate' is"];
Print[int, " ≈ ", 
        NumberForm[app, 10] ];
Print["The approximations obtained with Monte Carlo \
simulation:"];
Print["Using n =", 
        PaddedForm[
          (Subscript[val1, [[1]]]), 
          5], "    integral ≈ ", 
        NumberForm[
          (Subscript[val1, [[2]]]), 
          10] ];
Print["Using n =", 
        PaddedForm[
          (Subscript[val2, [[1]]]), 
          5], "    integral ≈ ", 
        NumberForm[
          (Subscript[val2, [[2]]]), 
          10] ];
Print["Using n =", 
        PaddedForm[
          (Subscript[val3, [[1]]]), 
          5], "    integral ≈ ", 
        NumberForm[
          (Subscript[val3, [[2]]]), 
          10] ];

We are done!

Aside.  We can use Mathematica's subroutine NIntegrate to perform Monte Carlo integration, and we can verify that it is computing the same results as our subroutine MonteCarloIntegral3D.

SeedRandom[5];
approx1 = MonteCarloIntegral3D[100]; 
\(\(Print["Integral  ≈  ", \ 
        approx1\_\(\([[\)\(2\)\(]]\)\)];\)\ \)

approx2  =  
        NIntegrate[
          Evaluate[f[{x, y, z}]], {x, 0, (9/10)}, {y, 0, 1}, {z, 0, 
            1  (11/10)}, Method -> MonteCarlo[5], 
          MaxPoints -> 100 ];
Print["Integral  ≈  ",  approx2];

We are really done!  

Aside.  The analytic value of the integral can be found and involves the following steps.

∫ (Sqrt[4 - (x^2) - (y^2) - 
              (z^2)]) \[DifferentialD]z

(∫_0^(11/10)) (Sqrt[4 - (x^2) - (y^2) - 
              (z^2)]) \[DifferentialD]z

∫((∫_0^(11/10)) (Sqrt[4 - (x^2) - (y^2) - 
                    (z^2)]) \[DifferentialD]z) \[DifferentialD]y

(∫_0^1)((∫_0^(11/10)) (Sqrt[4 - (x^2) - (y^2) - 
                    (z^2)]) \[DifferentialD]z) \[DifferentialD]y

∫((∫_0^1)((∫_0^(11/10)) (Sqrt[4 - 
                          (x^2) - (y^2) - 
                          (z^2)]) \[DifferentialD]z) \[DifferentialD]y) \[DifferentialD]x

todo  =  (∫_0^(9/10))((∫_0^1)((∫_0^(11/10)) (Sqrt[4 - (x^2) - (y^2) - 
                            (z^2)]) \[DifferentialD]z) \[DifferentialD]y) \[DifferentialD]x

todo=NumberForm[N[todo], 10]

todo = Chop[todo]

Example 12.  Let  f (x, y, z, u)  =  (Sqrt[5 - (x^2) - (y^2) - (z^2) - (u^2)]).  Use the Monte Carlo method to calculate approximations to the integral  

		\(∫\_0\%\(4/5\)\((∫\_0\%\(9/10\)\((∫\_0\%1\((∫\_0\%\(11/10\)\ \ \(\@\(5 - x\^2 - y\^2 - z\^2 - 
                                  u\^2\)\) du)\) dz)\) dy)\) dx\).

Solution 12.

f[{x_, y_, z_, 
            u_}] = (Sqrt[5 - (x^2) - (y^2) - (z^2) - (u^2)]); 
\(\(int = \*""(∫_0^(4/5))((∫_0^(9/10))((∫_0^1)((∫_0^(11/10)) (Sqrt[5 - (x^2) - (y^2) - (z^2) - (u^2)]) \[DifferentialD]u)\[DifferentialD]z)\[DifferentialD]y)\[DifferentialD]x"";\)\ \)
app  =  
        NIntegrate[
          f[{x, y, z, u}], {x, 0, (4/5)}, {y, 0, (9/10)}, {z, 0, 1}, {u, 0, 
            (11/10)}];  
Print["f[{x,y,z,u}] = ", 
        f[{x, y, z, u}] ]; 
Print["Mathematica's approximation using 'NIntegrate' is"];
Print[int, " ≈ ", 
        NumberForm[app, 10] ];

a1 = a2 = a3 = a4 = 0; 
b1 = (4/5);
b2 = (9/10);
b3 = 1;
b4 = (11/10); 
val1 = MonteCarloIntegral4D[100];  
val2 = MonteCarloIntegral4D[1000];  
val3 = MonteCarloIntegral4D[10000];  
Print["f[{x,y,z,u}] = ", 
        f[{x, y, z, u}] ];  
Print["Mathematica's approximation using 'NIntegrate' is"];
Print[int, " ≈ ", 
        NumberForm[app, 10] ]; 
Print["The approximations obtained with Monte Carlo \
simulation:"]; 
Print["Using n =", 
        PaddedForm[
          (Subscript[val1, [[1]]]), 
          5], "    integral ≈ ", 
        NumberForm[
          (Subscript[val1, [[2]]]), 
          10] ];
Print["Using n =", 
        PaddedForm[
          (Subscript[val2, [[1]]]), 
          5], "    integral ≈ ", 
        NumberForm[
          (Subscript[val2, [[2]]]), 
          10] ];
Print["Using n =", 
        PaddedForm[
          (Subscript[val3, [[1]]]), 
          5], "    integral ≈ ", 
        NumberForm[
          (Subscript[val3, [[2]]]), 
          10] ];

We are done!

Aside.  We can use Mathematica's subroutine NIntegrate to perform Monte Carlo integration, and we can verify that it is computing the same results as our subroutine MonteCarloIntegral4D.

SeedRandom[5];
approx1  =  
        MonteCarloIntegral4D[100]; 
\(\(Print["Integral  ≈  ", \ 
        approx1\_\(\([[\)\(2\)\(]]\)\)];\)\ \)

approx2  =  
        NIntegrate[
          Evaluate[f[{x, y, z, u}]], {x, 0, (4/5)}, {y, 0, (9/10)}, {z, 0, 
            1}, {u, 0, (11/10)}, Method -> MonteCarlo[5], 
          MaxPoints -> 100 ];
Print["Integral  ≈  ",  approx2];

Example 13.  Let  f (x, y, z, u, 
            w)  =  (Sqrt[6 - (x^2) - (y^2) - (z^2) - (u^2) - (w^2)]).  Use the Monte Carlo method to calculate approximations to the integral  

		\(∫\_0\%\(7/10\)\((∫\_0\%\(4/5\)\((∫\_0\%\(9/10\)\((∫\_0\%1\((∫\_0\%\(11/10\)\ \ \(\@\(6 - x\^2 - 
                                        y\^2 - z\^2 - u\^2 - 
                                        w\^2\)\) dw)\) du)\) dz)\) dy)\) dx\).

Solution 13.

f[{x_, y_, z_, u_, 
            w_}] = (Sqrt[6 - (x^2) - (y^2) - (z^2) - (u^2) - (w^2)]); 
\(\(int = \*""(∫_0^(7/10))((∫_0^(4/5))((∫_0^(9/10))((∫_0^1)((∫_0^(11/10)) (Sqrt[6 - (x^2) - (y^2) - (z^2) - (u^2) - (w^2)]) \[DifferentialD]w)\[DifferentialD]u)\[DifferentialD]z)\[DifferentialD]y)\[DifferentialD]x"";\)\ \)
app  =  
        NIntegrate[
          f[{x, y, z, u, w}], {x, 0, (7/10)}, {y, 0, (4/5)}, {z, 0, (9/10)}, {u, 
            0, 1}, {w, 0, (11/10)}];  
Print["f[{x,y,z,u,w}] = ", 
        f[{x, y, z, u, w}] ]; 
Print["Mathematica's approximation using 'NIntegrate' is"];
Print[int, " ≈ ", 
        NumberForm[app, 10] ];

a1 = a2 = a3 = a4 = a5 = 
                0;
b1 = (7/10);
b2 = (4/5);
b3 = (9/10);
b4 = 1;
b5 = (11/10); 
val1 = MonteCarloIntegral5D[100]; 
val2 = MonteCarloIntegral5D[1000]; 
val3 = MonteCarloIntegral5D[10000]; 
Print["f[{x,y,z,u,w}] = ", 
        f[{x, y, z, u, w}] ]; 
Print["Mathematica's approximation using 'NIntegrate' is"];
Print[int, " ≈ ", 
        NumberForm[app, 10] ];
Print["The approximations obtained with Monte Carlo \
simulation:"];
Print["Using n =", 
        PaddedForm[
          (Subscript[val1, [[1]]]), 
          5], "    integral ≈ ", 
        NumberForm[
          (Subscript[val1, [[2]]]), 
          10] ];
Print["Using n =", 
        PaddedForm[
          (Subscript[val2, [[1]]]), 
          5], "    integral ≈ ", 
        NumberForm[
          (Subscript[val2, [[2]]]), 
          10] ];
Print["Using n =", 
        PaddedForm[
          (Subscript[val3, [[1]]]), 
          5], "    integral ≈ ", 
        NumberForm[
          (Subscript[val3, [[2]]]), 
          10] ];

We are done!

Aside.  We can use Mathematica's subroutine NIntegrate to perform Monte Carlo integration, and we can verify that it is computing the same results as our subroutine MonteCarloIntegral5D.

SeedRandom[5];
approx1  =  
        MonteCarloIntegral5D[100]; 
\(\(Print["Integral  ≈  ", \ 
        approx1\_\(\([[\)\(2\)\(]]\)\)];\)\ \)

approx2  =  
        NIntegrate[
          Evaluate[f[{x, y, z, u, w}]], {x, 0, (7/10)}, {y, 0, (4/5)}, {z, 0, 
            (9/10)}, {u, 0, 1}, {w, 0, (11/10)}, Method -> MonteCarlo[5], 
          MaxPoints -> 100 ];
Print["Integral  ≈  ",  approx2];

Example 14.  Let  f (x, y, z, u, 
            w)  =  (Sqrt[6 - (x^2) - (y^2) - (z^2) - (u^2) - (w^2)]).  Use the Monte Carlo method to calculate approximations to the integral  

		\(∫\_0\%\(7/10\)\((∫\_0\%\(4/5\)\((∫\_0\%\(9/10\)\((∫\_0\%1\((∫\_0\%\(11/10\)\ \ \(\@\(6 - x\^2 - 
                                        y\^2 - z\^2 - u\^2 - 
                                        w\^2\)\) dw)\) du)\) dz)\) dy)\) dx\).  
		
Compare the Monte Carlo method with the midpoint rule.

Solution 14.

f[{x_, y_, z_, u_, 
            w_}] = (Sqrt[6 - (x^2) - (y^2) - (z^2) - (u^2) - (w^2)]); 
\(\(int = \*""(∫_0^(7/10))((∫_0^(4/5))((∫_0^(9/10))((∫_0^1)((∫_0^(11/10)) (Sqrt[6 - (x^2) - (y^2) - (z^2) - (u^2) - (w^2)]) \[DifferentialD]w)\[DifferentialD]u)\[DifferentialD]z)\[DifferentialD]y)\[DifferentialD]x"";\)\ \)
app  =  
        NIntegrate[
          f[{x, y, z, u, w}], {x, 0, (7/10)}, {y, 0, (4/5)}, {z, 0, (9/10)}, {u, 
            0, 1}, {w, 0, (11/10)}];  
Print["f[{x,y,z,u,w}] = ", 
        f[{x, y, z, u, w}] ]; 
Print["Mathematica's approximation using 'NIntegrate' is"];
Print[int, " ≈ ", 
        NumberForm[app, 10] ];

m = 6;
a1 = a2 = a3 = a4 = a5 = 
                0;
b1 = (7. /10);
b2 = (4. /5);
b3 = (9. /10);
b4 = 1.0;
b5 = (11. /10); 
xx = 
        Table[b1 (((i - 0.5))/m), {i, 1, m}];
yy = 
        Table[b2 (((i - 0.5))/m), {i, 1, m}];
zz = 
        Table[b3 (((i - 0.5))/m), {i, 1, m}];
uu = 
        Table[b4 (((i - 0.5))/m), {i, 1, m}];
ww = 
        Table[b5 (((i - 0.5))/m), {i, 1, m}];
\(\(f\&^\  = \ \(1\/m\^5\) \(∑\+\(i = 1\)\%m\(∑\+\(ii = 1\)\%m\(∑\+\(iii = 1\)\%m\(∑\+\(iiii = 1\)\%m\(∑\+\(iiiii = 1\)\%m f[{xx\_\(\([[\)\(i\)\(]]\)\), 
                        yy\_\(\([[\)\(ii\)\(]]\)\), 
                        zz\_\(\([[\)\(iii\)\(]]\)\), 
                        uu\_\(\([[\)\(iiii\)\(]]\)\), 
                        ww\_\(\([[\)\(iiiii\)\(]]\)\)}]\)\)\)\)\);\)\ \)
v = b1 b2 b3 b4 b5; 
\(\(midpoint = v*f\&^;\)\ \)
val4 = MonteCarloIntegral5D[7776]; 
Print["f[{x,y,z,u,w}] = ", 
        f[{x, y, z, u, w}] ]; 
Print["Mathematica's approximation using 'NIntegrate' is"];
Print[int, " ≈ ", 
        NumberForm[app, 10] ];
\(\(Print[\*""Using  \!\(m\^5\) = "", 
        m\^5, "  function evaluations for the midpoint rule"];\)\)
Print[int, " ≈ ", 
        NumberForm[midpoint, 10] ];
\(\(Print["Using  n =", 
        PaddedForm[
          val4\_\(\([[\)\(1\)\(]]\)\), 
          5], "  function evaluations for the Monte Carlo method"\ ];\)\)
Print[int, " ≈ ", 
        NumberForm[
          (Subscript[val4, [[2]]]), 
          10] ];
Print[""]; 
Print["The midpoint error = ", 
        PaddedForm[app - midpoint] ];
\(\(Print["Monte Carlo  error = ", 
        PaddedForm[
          app - val4\_\(\([[\)\(2\)\(]]\)\), 8]\ ];\)\)

m = 7;
b1 = (7. /10);
b2 = (4. /5);
b3 = (9. /10);
b4 = 1.0;
b5 = (11. /10); 
xx = 
        Table[b1 (((i - 0.5))/m), {i, 1, m}];
yy = 
        Table[b2 (((i - 0.5))/m), {i, 1, m}];
zz = 
        Table[b3 (((i - 0.5))/m), {i, 1, m}];
uu = 
        Table[b4 (((i - 0.5))/m), {i, 1, m}];
ww = 
        Table[b5 (((i - 0.5))/m), {i, 1, m}];
\(\(f\&^\  = \ \(1\/m\^5\) \(∑\+\(i = 1\)\%m\(∑\+\(ii = 1\)\%m\(∑\+\(iii = 1\)\%m\(∑\+\(iiii = 1\)\%m\(∑\+\(iiiii = 1\)\%m f[{xx\_\(\([[\)\(i\)\(]]\)\), 
                        yy\_\(\([[\)\(ii\)\(]]\)\), 
                        zz\_\(\([[\)\(iii\)\(]]\)\), 
                        uu\_\(\([[\)\(iiii\)\(]]\)\), 
                        ww\_\(\([[\)\(iiiii\)\(]]\)\)}]\)\)\)\)\);\)\ \)
v = b1 b2 b3 b4 b5; 
\(\(midpoint = v*f\&^;\)\ \)
val5 = MonteCarloIntegral5D[16807]; 
Print["f[{x,y,z,u,w}] = ", 
        f[{x, y, z, u, w}] ]; 
Print["Mathematica's approximation using 'NIntegrate' is"];
Print[int, " ≈ ", 
        NumberForm[app, 10] ];
\(\(Print[\*""Using \!\(m\^5\) = "", 
        m\^5, "  function evaluations for the midpoint rule"];\)\)
Print[int, " ≈ ", 
        NumberForm[midpoint, 10] ];
\(\(Print["Using  n =", 
        PaddedForm[
          val5\_\(\([[\)\(1\)\(]]\)\), 
          5], "  function evaluations for the Monte Carlo method"\ ];\)\)
Print[int, " ≈ ", 
        NumberForm[
          (Subscript[val5, [[2]]]), 
          10] ];
Print[""]; 
Print["The midpoint error = ", 
        PaddedForm[app - midpoint] ];
\(\(Print["Monte Carlo  error = ", 
        PaddedForm[
          app - val5\_\(\([[\)\(2\)\(]]\)\), 8]\ ];\)\)

Remark.  For this quintuple multiple integral,  most of the time the Monte Carlo simulation has less error than the midpoint rule.  This is what we expect to happen.

Research Experience for Undergraduates

Monte Carlo Integration  Internet hyperlinks to web sites and a bibliography of articles.