Mathematics made simple.

Followers

Monte Carlo Integration using Python


 Monte Carlo Integration

Monte Carlo integration is a  powerful method of computing the value of complex integrals using probabilistic techniques. This technique uses random numbers to compute the definite integral of a function. Here we are going to use the python programs written in the previous post to generate pseudorandom numbers and approximate the value of the definite integral of a function.
Consider a function to be integrated, as shown below:

That is we need to evaluate the definite integral $\int_{a}^{b} f(x)\mathrm{d}x$. The integral is just the area under the curve. The width of the interval $(b- a)$ times the average value of the function is also the value of the integral, that is,

$\int_{a}^{b}f(x)\mathrm{d}x = (b - a)f_{average} = (b-a)\langle f \rangle$

So if we had some independent way of calculating the average value of the integrand, then we could evaluate the integral. This is where random numbers come in.
Imagine that we had a list of random numbers, $x_j$ uniformly distributed between $a$ and $b$. To calculate the function average, we simply evaluate $f(x_j)$ at each of the randomly selected points, add all of them together and divide the sum by the number of points.

$$\langle f \rangle _N = \frac{1}{N}\sum_{i=1}^{N}f(x_i)$$

As the number of points used in calculating the average increases, $\langle f \rangle_N$ approaches the true value, $\langle f \rangle$. Therefore, as a numerical approximation, we can write

$$\int_{a}^{b}f(x)\mathrm{d}x = \frac{(b-a)}{N}\sum_{i=1}^{N}f(x_i)$$

Alternatively, we can look at the Monte Carlo integration method in the following way:
  • Find some value $M$ such that $f(x) < M$ over the integral $[a,b]$.
  • Select random numbers $x$ and $y$ from the intervals $[a,b]$ and $[0,M]$ respectively.
  • Determine if $y > f(x)$ or $y \leq f(x)$.
  •  Repeat the process $N$ times, keeping track of the number of times $y \leq f(x)$ or under the curve (successes); call the total number of successes $S$.
The estimated probability of success is

$$\frac{S}{N} = \frac{Area \ under \ curve}{Total \ area \ inside \ rectangle} = \frac{\int_{a}^{b}f{x}\mathrm{d}x}{M(b-1)}$$

Thus the value of the integral can be calculated from the formula

$$\int_{a}^{b}f(x)\mathrm{d}x = M(b-a)\frac{S}{N}$$

The python program that implements the Monte Carlo method of integration is given below. Here the random numbers used in the integration is generated using Lehmer's method mention in the previous post.

Program

from sympy import *
def pseudo(a,m,x0):
    list1 = []
    for i in range(100):
        x0 = (a*x0)%m
        list1.append(x0)
    return(list1)

def integrate(f,var):
    f = sympify(f)
    var = Symbol(var)
    rand_values1 = pseudo(61,101,18)
    rand_values2 = pseudo(63,107,18)
    f_values1 = []
    f_values2 = []
    x = []
    y = []
    for i in rand_values1:
        x.append(a + (((b-1)/101)*i))
    i = a
    #to find maximum M
    while i <= b:
        f_values1.append(f.subs({var:i}))
        i += 0.01
    #to check whether y lies below or above the curve
    for i in x:
        f_values2.append(f.subs({var:i}))
    M = int(max(f_values1)) + 1
    for i in rand_values2:
        y.append((M/107)*(i))
    count = 0
    for i in range(len(y)):
        if y[i] <= f_values2[i]:
            count += 1
    return (b-a)*M*(count/100)

f = input("Enter the function: ")
var = input("Enter the variable: ")
a = int(input("Enter the lower limit: "))
b = int(input("Enter the upper limit: "))
print("The value of the definite integral is", integrate(f,var))


The value of some definite integrals obtained using the above program and their actual values are given below in the table:


0 comments:

Post a Comment

Ultimate Theorem. Powered by Blogger.

About Me

My photo
Mathematics has always fascinated me. I love the subject since childhood. This love towards the field helped me in completing a Masters Degree in Mathematics. As much as I love the subject, I love teaching it too.

Local Linear Approximation

You might have seen the famous approximation $\sin(x) \approx x$ for $x$ near 0. We can use derivatives to approximate non-linear functions ...

Search This Blog

ultimate theorem