Monte Carlo simulation offers a simple numerical method for calculating the area under a curve where one has the equation of the curve, and the limits of the range for which we wish to calculate the area.
For example, imagine we wish to perform the integral:
| \int\limits_0^{10} x^3 ~ sin ^ 2 (x) ~ exp ~ \bigg(\frac{-1}{x}\bigg)dx |
Figure 1 shows the curve x^3 ~ sin^2 (x) ~ exp ~ \bigg(\frac {-1}{x}\bigg)between 0 and 10. Note that y <500 over this range.
We can estimate the area under the curve as follows:
Generate a random value for x on [0,10]: =Uniform(0,10)
Generate a random value for y on [0,500]: =Uniform(0,500)
Determine the value of the curve (ycurve) at x: =x^3*SIN(x)^2*EXP(-1/x)
Test whether the random y is below the curve (test): =IF(y<ycurve,1,0)
Determine the fraction of random samples falling below the curve (fraction): =CB.GetForeStatFN(test,2)
Multiply this fraction by the area of the rectangle (0,0; 10;500) = fraction*(500-0)*(10-0)
This is the required area. You can check the model here: Function to Integrate - 1.
The links to the Function to Integrate - 1 software specific models are provided here:
Functions crossing the x-axis
Where the function y can be negative (e.g. Figure 2) we often want the area between the curve and the x-axis, which is easily achieved with a simple ABS( ) modification to the IF statement.
So, for example, determining
| \int\limits_0^{10} \frac{\sin(x)}{log(x)} dx |
proceeds as follows (noting |y| does not exceed 25 over x = [0,10]):
Generate a random value for x on [0,10]: =Uniform(0,10)
Generate a random value for y on [0,25]: =Uniform(0,25)
Determine the value of the curve ycurve at x: =SIN(x)/LOG(x)
Test whether the random y is below the curve (test): =IF(y<ABS(ycurve),1,0)
Determine the fraction of random samples falling below the curve (fraction): =CB.GetForeStatFN(test,2)
Multiply this fraction by the area of the rectangle (0,0; 10;25) = fraction*(25-0)*(10-0)
This is the required area. You can check the model here: Function to Integrate - 2.
The links to the Function to Integrate - 2 software specific models are provided here:
A refinement to this method is to split the x-range up into increments each with a different ymax so that the simulation is more efficient. For example, Figure 1 could be split into three separate areas as shown in Figure 3. Then the model simulates each area and adds the results as shown below:
You can check the model here: Function to Integrate - 3.
The links to the Function to Integrate - 3 software specific models are provided here:
This method can produce any required level of accuracy (in principle, at least) by simply performing more iterations. It has the advantage of being able to integrate any area where the formula for y can be coded easily, but it is of course quite slow and needs rerunning if one changes the values of the equation's parameters. It also has the advantage of being unaffected by singularities at the limits of the definite integral because any Uniform(a,b) distribution has zero probability of generating exactly a or b.
Integrating over three dimensions and more
The method can also be extended to definite integrals over three dimensions very easily. For example, the equation for the surface of a sphere with center (0,0,0) and radius r is:
x2+y2+z2=r2
The volume of a sphere V is \frac {4}{3} \pi r^3 .
For r = 4, that gives V = 268.08257.. If we didn't have the formula, we could get to this value by randomly generating a point in a cube of length 8 centered at the origin, i.e.:
Within =
IF (x^2 + y^2 + z^2 < 4^2, 1, 0)
The mean of the variable 'Within' equals the fraction of iterations that a generated point falls within the cube. Finally we determine the volume of the sphere, equal to the product of the fraction of generated points falling within the sphere, and the volume of the cube:
Where 8^3 is the volume of the cube within which we are sampling.
