Introduction

For a while, while learning about Numerical Methods for integral, solutions of differential equation, and some other mathematical object, I've come across the Monte Carlo methods. It seemed like some distant method by which we could magically solve a problem. Of course, I knew better, but I avoided it, until I read the book "Partial Differential Equations for Scientists and Engineers" by Stanley J. Farlow.

He explained that Monte Carlo Methods, basically, are methods by which we numerically solve a deterministic problem using probabilistic means. For example, if we try to calculate the integral $\int_0^1 x^2 dx$. We can easily calculate it. $$\int_0^1 x^2 dx=\frac{x^3}{3}\bigg|_0^1=\frac{1}{3}$$

It can also be calculated by Monte Carlo Methods. If we plot the function $x^2$

Created using Octave 4.4.1

We know that the integral of the function is the area that is under the curve of the function. So, if we guess in a uniform manner the x and y coordinates on the plot between $0\le x\le 1$ and $0\le y \le 1$, where the y interval coordinates are based on the max and the min of the function in that interval, then, with sufficient guesses, we know that the $\frac{1}{3}$ of the guesses are going to be under the curve of the function.

Why?

The significance of this is the ability to numerically calculate integrals which are more difficult to calculate, such as the integral $\int_0^1 e^{x^2}$ which doesn't have an obvious analytical solution which doesn't use a special function such as the Error function, which is defined by $erf(x)=\frac{1}{\sqrt{\pi}} \int_{-x}^{x} e^{-t^2} dt$. We can however calculate the integral easily numerically using the Monte Carlo Method.
Monte Carlo Methods are used in various areas of engineering and physicical sciences such as microelectronic engineering, geostatistics molecular dynamics, fluid dynamics, and even marketing analytics. (link1, link2)


Code

I will be using Octave 4.4.1 which is a free alternative to matlab lacking some of the functionality, but for our uses it will work fine. Download link.

Let's first declare the $f(x)$ function as an anonymous function.

fun = @(x) exp(x.^2);

Set the number of guesses to 1000

% Number of guesses
n = 1000;

We need to set the range for the function on both axes. The x-axis is easy, it's the the interval for the integral; The y-axis in the other hand depends on the function. For this function, we know that the max is at the large end of the interval at $f(1)$.

% x axis
xlow = 0;
xhigh = 1;
% y axis
ylow = 0;
yhigh = fun(xhigh);

The variable to be used to count the number of guesses that are under the is

% Count of guesses under curve
j = 0;

Finally, the loop which will generate the $x$ and $y$ guess, as well as pass the $x$ guess throught the function $f(x)$. Lastly, guessed and actual values are compared so, if the guess is under the curve then $j$ is incremented by 1.

for i = 1:n
  % generate random number x
  x = xlow + (xhigh - xlow) * rand();
  % generate random guess y
  gy = ylow + (yhigh - ylow) * rand();
  % find real y
  fx = fun(x);
  
  % Compare
  if (gy < fx)
    j = j + 1;
  endif
endfor

Calculate the fraction under the curve and multiply it by the area where we were guessing.

% Calculate fraction
frac = j / n;
% Calculate Area
area = (xhigh - xlow) * (yhigh - ylow);
value = frac * area

Running the code I obtained $1.452$ as an answer compared to the actual $1.462$.

$\pi$

Another cool application is in the calculation of $pi$. We know that the equation for a circle is $x^2+y^2=r^2$, and set $r$, the radius, to $10$. By simplifying and solving for $y$, $y = \sqrt{100-x^2}$. Another given is the area of a circle $A=\pi r^2=100\pi$. If we use this monte carlo method to find the area under the function where $0\le x \le 10$ and $0\le x \le 10$ we can find the area of a quarter of a circle; I will call it $I$ for integral and $I=\frac{A}{4}$. So, $$\pi=\frac{I}{25}$$ where, in this case, $$f(x)=\sqrt{100-x^2}$$

Plot of $f(x)$, the quarter circle

Code

Let's modify out last code. We will first place the new $f(x)$

fun = @(x) sqrt(100 - x.^2);

The new ranges

xlow = 0;
xhigh 10;
ylow = 0;
yhigh = 10;

After calculating the value for the integral let's divide that by 25

% /\ The rest of the code /\
value = frac * area
value / 25

I obtained the value of $3.1440$ compared to the real value of $3.14159...$

Source Code