Modelling and Simulation: Theory and Examples

4. Monte Carlo Simulation

Sometimes it is difficult to obtain an analytical solution to a problem. Monte Carlo methods rely on repeated random sampling to obtain numerical results. Randomness may be used to solve some deterministic problems.

4.1. Calculate the approximate value of \(\pi\).

Example 1

Use the Monte Carlo method to evaluate

(a) the approximate value of \(\pi\).
(b) the approximate area of a triangle without \(\frac{1}{2} bh\).
(c) \(\int_{-3}^{3} e^{-x^2}\)
(d) \(\int_0^{15} \Big( \big(7x^3 + 20x^2 + 45x + 5 \big)^{\frac{1}{3}} \Big) e^{-\frac{1}{5}x} dx \)

Solution:

(a) the approximate value of \(\pi\)

\[\begin{align*} \frac{A_{square}} {A_{circle}} &= \frac{(2rad)^2}{\pi \times rad^2} \\ \pi &= \frac{(2rad)^2}{rad^2} \times \frac {A_{circle}} {A_{square}} \\ \pi &= \frac{(2rad)^2}{rad^2} \times \frac {N_{circle}} {N_{square}} \end{align*}\]

\(N\) is the number of dots in a square or a circle

You may use the equation of a circle \(r^2 = x^2 + y^2\) to know when a point is inside the circle.

 1import numpy as np
 2import matplotlib.pyplot as plt
 3
 4inCirc = 0               # inside circle counts 
 5trial = 1000             # number of trials 
 6rad = 2                  # circle radius 
 7x0, y0 = 1, 15           # circle center
 8x_circ = []; y_circ = [] # inside circle points
 9x_sq = []; y_sq = []     # square points 
10
11for i in range(1,trial):
12    x = np.random.uniform(x0-rad,x0+rad)
13    y = np.random.uniform(y0-rad,y0+rad)
14    if ((x-x0)**2 + (y-y0)**2 <= rad**2):
15        inCirc += 1
16        x_circ.append(x); y_circ.append(y)
17    else:
18        x_sq.append(x); y_sq.append(y)
19
20print("Estimated Prob = %.5f" %(inCirc/trial))
21print("Real Prob = %.5f" %(np.pi*rad**2/((2*rad)**2))) 
22
23print("Estimated Area = %.5f" %( ((2*rad)**2)*(inCirc/trial) ))
24print("Real Area = %.5f" %(np.pi*rad**2))
25
26print("Estimated PI = %.5f" %(((2*rad)**2)*(inCirc/trial)/rad**2))
27print("Real PI = %.5f" %(np.pi)) 
28
29fig, ax = plt.subplots()
30ax.scatter(x_sq,y_sq, marker='o', color='k')
31ax.scatter(x_circ,y_circ, marker='o', color='b')
32ax.set_xlabel('X')
33ax.set_ylabel('Y')
34ax.set_aspect('equal')
35#plt.savefig(r'../images/monte_pi.png')
36plt.show()

Estimated Prob = 0.79600
Real Prob = 0.78540
Estimated Area = 12.73600
Real Area = 12.56637
Estimated PI = 3.18400
Real PI = 3.14159

Figure 4. Monte Carlo \(\pi\) calculation

(b) the area of a triangle without \(\frac{1}{2} bh\)

\[\begin{align*} \frac{Triangle_A}{Square_A} &= \frac{Triangle_{dots}}{Square_{dots}} \\ Triangle_A &= \frac{Triangle_{dots}}{Square_{dots}} \times Square_A \end{align*}\]

You may use the equation of a straight line \(y = mx + c\) or just the gradient to know when a point is inside or outside the desired triangle.

 1import numpy as np
 2import matplotlib.pyplot as plt 
 3
 4xt0, xtl = 0, 2      # x width
 5yt0, ytl = 0, 2      # y height
 6xta = []; yta = []   # points inside triangle
 7x = []; y = []       # all points
 8xtat = []; ytat = [] # points outside triangle
 9Tdots = 0            # dots count inside triangle 
10Sdots = int(1e3)     # number of trials 
11
12def h():
13    return np.sqrt(xt**2 + yt**2)
14
15for i in range(Sdots):
16    xt = np.random.uniform(xt0,xtl) 
17    yt = np.random.uniform(yt0,ytl) 
18    x.append( np.random.uniform(xt0,xtl) )
19    y.append( np.random.uniform(yt0,ytl) )
20    if (yt/xt) <= (ytl/xtl) : # gradient 
21        xta.append(xt)
22        yta.append(yt)
23        Tdots += 1
24    else:
25        xtat.append(xt)
26        ytat.append(yt)
27
28fig, ax = plt.subplots()
29ax.scatter(xtat,ytat, marker='.', color='blue')
30ax.scatter(xta,yta, marker='.', color='black')
31ax.set_xlabel('X')
32ax.set_ylabel('Y')
33ax.set_aspect('equal')
34TA = (Tdots/Sdots) * (xtl*ytl)
35print("Triangle Area = %.4f" %TA)
36#plt.savefig(r'../images/monte_area.png')
37plt.show()

Triangle Area = 2.0000

Figure 5. Monte Carlo area calculation

\[\begin{align*} I &= \int_a^b g(x) dx \\ I_n &\approx \frac{1}{n} \sum_{i=1}^n I_i \\ &\approx \frac{b-a}{n} \sum_{i=1}^n g(a + (b-a) U_i) \end{align*}\]

 1import numpy as np
 2
 3def monte(fxn, a=-3, b=3, N=1000):
 4    p = []
 5    for i in range(N):
 6        p.append(np.random.uniform(a,b))
 7
 8    integral = 0
 9    for i in p:
10        integral += fxn(i)
11
12    return ((b-a)/N) * integral 
13
14def f(x):
15    return np.exp(-(x)**2) 
16
17print("Monte Integral %.2e" %monte(f, -3, 3))

Monte Integral 1.80e+00

evaluating \(\int_{-3}^3 e^{-x^2}\) using np.scipy

 1import numpy as np
 2import matplotlib.pyplot as plt
 3from scipy import integrate
 4
 5x = np.arange(-3,3+.1, step=0.1)
 6y = ( np.exp(-(x)**2)) 
 7
 8def f(x):
 9    return np.exp(-(x)**2)
10
11val, err = integrate.quad(f, -3, 3)
12print("Scipy Integral = %.2e" %val)
13
14fig, ax = plt.subplots()
15ax.set_xlabel('X')
16ax.set_ylabel('Y')
17ax.fill_between(x,y, hatch="x", color="none", edgecolor="blue")
18ax.plot(x,y, color="red", lw=3)
19
20#plt.savefig(r'../images/scipy_int.png')
21plt.show()

Scipy Integral = 1.77e+00

Figure 6. Monte Carlo area calculation

Add Comment

Comments

No comments yet. Be the first!