CALCULATE INTEGRALS BY MONTE-CARLO METHODDuong Ngoc HaoUniversity of Information TechnologyVietnam 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) <=dhaodn[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:restart;
read `D:/haodn_e.txt`;
with(haodn);Suppose we want to determine LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYoLUkjbWlHRiQ2JVEkaW50RicvJSdpdGFsaWNHUSV0cnVlRicvJSxtYXRodmFyaWFudEdRJ2l0YWxpY0YnLUkjbW9HRiQ2LVEiKEYnL0YzUSdub3JtYWxGJy8lJmZlbmNlR0YxLyUqc2VwYXJhdG9yR1EmZmFsc2VGJy8lKXN0cmV0Y2h5R0YxLyUqc3ltbWV0cmljR0Y/LyUobGFyZ2VvcEdGPy8lLm1vdmFibGVsaW1pdHNHRj8vJSdhY2NlbnRHRj8vJSdsc3BhY2VHUSwwLjE2NjY2NjdlbUYnLyUncnNwYWNlR0ZMLUYsNiVRInhGJ0YvRjItRjY2LVEnJnNkb3Q7RidGOS9GPEY/Rj0vRkFGP0ZCRkRGRkZIL0ZLUSYwLjBlbUYnL0ZORlgtRiw2JVEmeWR4ZHlGJ0YvRjJGOQ==LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYtLUkobXN1YnN1cEdGJDYnLUkjbW9HRiQ2LVEmJmludDtGJy8lLG1hdGh2YXJpYW50R1Enbm9ybWFsRicvJSZmZW5jZUdRJmZhbHNlRicvJSpzZXBhcmF0b3JHRjcvJSlzdHJldGNoeUdRJXRydWVGJy8lKnN5bW1ldHJpY0dGNy8lKGxhcmdlb3BHRjwvJS5tb3ZhYmxlbGltaXRzR0Y3LyUnYWNjZW50R0Y3LyUnbHNwYWNlR1EmMC4wZW1GJy8lJ3JzcGFjZUdGRy1GIzYnLUkjbW5HRiQ2JFEiMEYnRjIvJSdpdGFsaWNHRjwvJStmb3JlZ3JvdW5kR1EsWzEyOCwwLDEyOF1GJy8lLHBsYWNlaG9sZGVyR0Y8L0YzUSdpdGFsaWNGJy1GIzYnLUZNNiRRIjFGJ0YyRlBGUkZVRlcvJTFzdXBlcnNjcmlwdHNoaWZ0R1EiMEYnLyUvc3Vic2NyaXB0c2hpZnRHRmpuLUYvNi1RIn5GJ0YyRjVGOC9GO0Y3Rj0vRkBGN0ZBRkNGRUZILUYvNi1RMCZEaWZmZXJlbnRpYWxEO0YnRjJGNUY4RmBvRj1GYW9GQUZDRkVGSC1JI21pR0YkNidRInhGJ0ZQRlJGVUZXRl1vRistRmZvNiVRI3h5RidGUEZXRl1vRmJvLUZmbzYlUSJ5RidGUEZXRjI=, we we will use procedure intm2 as followresult:=intm2([x*y,x,y,0,1,0,1,0,1,1000]);JSFHIt's is easy to determine the exact value of the integral above to compare with the approximated value, exact_value:=int(int(x*y,y=0..1),x=0..1);JSFHBy this way, we can easily compute more complicated integral like LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYtLUkobXN1YnN1cEdGJDYnLUkjbW9HRiQ2LVEmJmludDtGJy8lLG1hdGh2YXJpYW50R1Enbm9ybWFsRicvJSZmZW5jZUdRJmZhbHNlRicvJSpzZXBhcmF0b3JHRjcvJSlzdHJldGNoeUdRJXRydWVGJy8lKnN5bW1ldHJpY0dGNy8lKGxhcmdlb3BHRjwvJS5tb3ZhYmxlbGltaXRzR0Y3LyUnYWNjZW50R0Y3LyUnbHNwYWNlR1EmMC4wZW1GJy8lJ3JzcGFjZUdGRy1JI21pR0YkNidRImFGJy8lJ2l0YWxpY0dGPC8lK2ZvcmVncm91bmRHUSxbMTI4LDAsMTI4XUYnLyUscGxhY2Vob2xkZXJHRjwvRjNRJ2l0YWxpY0YnLUZLNidRImJGJ0ZORlBGU0ZVLyUxc3VwZXJzY3JpcHRzaGlmdEdRIjBGJy8lL3N1YnNjcmlwdHNoaWZ0R0Zmbi1GLDYnRi4tRiM2Jy1GSzYlUSJjRidGTkZVRk5GUEZTRlUtRiM2JS1GSzYlUSJkRidGTkZVRk5GVUZaRmduLUZLNiVRInhGJ0ZORlUtRi82LVEnJnNkb3Q7RidGMkY1RjgvRjtGN0Y9L0ZARjdGQUZDRkVGSC1GSzYlUSRleHBGJy9GT0Y3RjItSShtZmVuY2VkR0YkNiQtRiM2LC1GLzYtUSomdW1pbnVzMDtGJ0YyRjVGOEZbcEY9RlxwRkFGQy9GRlEsMC4yMjIyMjIyZW1GJy9GSUZqcC1JJW1zdXBHRiQ2JS1GSzYlRk1GTkZVLUYjNiUtSSNtbkdGJDYkUSIyRidGMkZORlVGWi1GXXE2JUZlb0ZhcUZaLUYvNi1RKCZtaW51cztGJ0YyRjVGOEZbcEY9RlxwRkFGQ0ZpcEZbcS1GXXE2JS1GSzYlRllGTkZVRmFxRlpGaG8tRks2JVEkY29zRidGYHBGMi1GXXE2JS1GYnA2JC1GIzYkLUZLNiVRInlGJ0ZORlVGMkYyRmFxRlotRl1xNiVGZW8tRiM2JS1GZHE2JFEiNEYnRjJGTkZVRlpGMkYyLUYvNi1RIn5GJ0YyRjVGOEZbcEY9RlxwRkFGQ0ZFRkgtRi82LVEwJkRpZmZlcmVudGlhbEQ7RidGMkY1RjhGW3BGPUZccEZBRkNGRUZILUZLNidGZ29GTkZQRlNGVS1GSzYlUSNkeUYnRk5GVUYy.JSFH