CALCULATE INTEGRALS BY MONTE-CARLO METHOD

Duong Ngoc Hao

University of Information Technology

Vietnam national university, Ho Chi Minh city- Vietnam.

We build two procedures to determine integrals in one- dimenssion and two-dimenssion by Monte-Carlo method. These procedures are contained in

packadge "haodn_e.txt" and suppose to be stored in drive D:

# Calculate integrals by Monte-Carlo method.

# Written by Duong Ngoc Hao, University of Information Technology, Vietnam national university, Hochiminh city.

# Truong dai hoc cong nghe thong tin, dai hoc quoc gia Tp HCM.

# Updated 1 ngay 7.3.2009

# Updated 2 ngay 22.3.2010

# Determine integral of the function f(x) by Monte-carlo method.

# Syntax: intm1([f,x,a,b,c,d,n])

# where f is an express of the function, x is a variable

# n= number of random points

# a,b means a<= x <=b

# c,d means c<= y=f(x) <=d

haodn[intm1]:=proc(L::list(algebraic))

        local m,n,i,j,B,H,a,b,c,d,rf,MM,LN,hits,estimate,Data;

        n:=nops(L);

        m:=5;

        c:=L[n-3]-L[n-4];

        d:=L[n-1]-L[n-2];

        H:=unapply(L[1],L[2]);

        MM:=1000*L[n];

        for j from 1 to m do

                rf:=rand(0..MM);

                hits:=0;

                for i to L[n] do

                        a:=evalf(L[n-4]+c*rf()/MM);

                        b:=evalf(L[n-2]+d*rf()/MM);

                        B:=evalf(H(a));

                        if b<B then hits:=hits+1;

                        fi;

                od;

                estimate[j]:=evalf(c*d*hits/L[n]);

        od;

        Data:=[seq(estimate[j],j=1..m)];

        estimate:=(1/m)*sum(Data[k],k=1..m);

        end:

# Determine double integral of the function f(x,y) by Monte-Carlo method

# Syntax:  intm2([f,x,y,a1,a2,b1,b2,c1,c2,n])

# where f is an express of the function; x,y are two independent variables.

# n= so diem lay ngau nhien

#[a1,a2] is closed interval of x

#[b1,b2] is closed interval of y

#[c1,c2] is closed interval of z=f(x,y)

haodn[intm2]:=proc(L::list(algebraic))

        local m,n,k,i,j,B,H,a,b,c,aa,bb,cc,g,rf,MM,LN,hits,estimate,Data;

        n:=nops(L);

        aa:=L[5]-L[6];            #Do dai khoang tren Ox

        bb:=L[n-3]-L[n-4];        #Do dai khoang tren Oy

        cc:=L[n-1]-L[n-2];        #Do dai khoang tren Oz

        H:=unapply(L[1],L[2],L[3]);

        MM:=1000*L[n];

        m:=5;

        k:=L[n]^2;

        rf:=rand(0..MM);

        for j from 1 to m do

                hits:=0;

                for i to k do

                        a:=evalf(L[n-4]+aa*rf()/MM);

                        b:=evalf(L[n-2]+bb*rf()/MM);

                        g:=evalf(L[n-6]+cc*rf()/MM);

                        B:=evalf(H(a,b));

                        if g<B then hits:=hits+1;

                        fi;

                od;

                estimate[j]:=evalf(aa*bb*cc*hits/k);

        od;

        Data:=[seq(estimate[j],j=1..m)];

        estimate:=(1/m)*sum(Data[d],d=1..m);

        end:

Example:

Suppose we want to determine  , we we will use procedure intm2 as follow

It's is easy to determine the exact value of the integral above to compare with the approximated value,

By this way, we can easily compute more complicated integral like .