大學(xué)數(shù)學(xué)實驗之蒙特卡洛方法

1、數(shù)學(xué)實驗報告班級： 序號： 姓名： 1 問題描述I、用蒙特卡羅方法計算以下函數(shù)在區(qū)間上的積分，并改變隨機點數(shù)目觀察對結(jié)果的影響。（1） y=1/(1+x), 0=<x=<1;（2） y= (exp(3*x)*sin(2*x), 0=<x=<2 ;（3） y=(1+x2)0.5, 0=<x=<2;（4） y=(1/(2*pi)0.5)*exp(-x(i)2/2), 0=<x=<2;（5） y=exp(x(i)/2)*(sin(x(i)2, 0=<x=<2*pi;（6） f(x,y)=exp(-x2-y2) 0=<x=<pi,

2、 0=<y=<sin(x);II、用蒙特卡羅法求解全局最優(yōu)化及約束問題并通過圖形做出評論，求下列函數(shù)的最大值。（1） f(x)=(1-x.2).*sin(3*x), -2*pi=<x=<2*pi;（2） max f(x)=x1*x2*x3,s.t.:-x1+2x2+2x3>=0,x1+2x2+2x3<=72,10<=x2<=20,x1-x2=10;（3） f(x,y)=(X.2+2*(Y.2)+X.*Y).*exp(-X.2-Y.2), abs(x)<1.5,abs(y)<1.5;2 問題分析與實驗過程I、(1)使用均值估計法程序：f

3、unction p=shell1(a,b,n)z=0;x=unifrnd(a,b,1,n);for i=1:n u=(x(i)+1)(-1); z=z+u;endp=(b-a)*z/n;運行結(jié)果：p=shell1(0,1,1000)p = 0.6975>> p=shell1(0,1,10000)p = 0.6922>> p=shell1(0,1,100)p = 0.7001>> p=shell1(0,1,500)p =0.6890 結(jié)果分析：改變了四次隨機點數(shù)，結(jié)果都趨近于0.69，說明積分值約等于0.69，但是點數(shù)越多，值越接近。 I、(2)使用均值估計法

4、程序：function p=shell2(a,b,n)z=0;x=unifrnd(a,b,1,n);for i=1:n u=(exp(3*x(i)*sin(2*x(i); z=z+u;endp=(b-a)*z/n;運行結(jié)果：>> p=shell2(0,2,1000)p = -24.4911>> p=shell2(0,2,100)p = -43.8720>> p=shell2(0,2,10000)p = -30.8699>> p=shell2(0,2,500)p = -23.2955>> p=shell2(0,2,100000)p =

5、-30.0058結(jié)果分析：改變了5次隨機點數(shù)，結(jié)果變化較大，但是點數(shù)越多，值越接近真實積分值。所以積分值近似于-30。I、(3)使用均值估計法程序：function p=shell3(a,b,n)z=0;x=unifrnd(a,b,1,n);for i=1:n u=(1+x(i)2)0.5; z=z+u;endp=(b-a)*z/n;運行結(jié)果：>> p=shell3(0,2,100)p = 2.9293>> p=shell3(0,2,1000)p = 2.9516>> p=shell3(0,2,10000)p = 2.9512>> p=shell

6、3(0,2,100000)p = 2.9600結(jié)果分析：改變了四次隨機點數(shù)，結(jié)果都趨近于2.95，說明積分值約等于2.95，而且點數(shù)越多，值越接近真實積分值。I、(4)使用均值估計法程序：function p=shell4(a,b,n)z=0;x=unifrnd(a,b,1,n);for i=1:n u=(1/(2*pi)0.5)*exp(-x(i)2/2); z=z+u;endp=(b-a)*z/n;運行結(jié)果：>> p=shell4(0,2,100000)p = 0.4783>> p=shell4(0,2,10000)p = 0.4777>> p=shel

7、l4(0,2,1000)p = 0.4765>> p=shell4(0,2,100)p = 0.4432結(jié)果分析：改變了四次隨機點數(shù)，結(jié)果都趨近于0.47，說明積分值約等于0.47，而且點數(shù)越多，值越接近真實積分值。I、(5)使用均值估計法程序：function p=shell5(a,b,n)z=0;x=unifrnd(a,b,1,n);for i=1:n u=exp(x(i)/2)*(sin(x(i)2; z=z+u;endp=(b-a)*z/n;運行結(jié)果：>> p=shell5(0,2*pi,100)p = 22.0140>> p=shell5(0,2*

8、pi,1000)p = 20.2718>> p=shell5(0,2*pi,10000)p = 20.9394>> p=shell5(0,2*pi,100000)p = 20.7968結(jié)果分析：改變了四次隨機點數(shù)，結(jié)果都趨近于20.8，說明積分值約等于20.8，而且點數(shù)越多，值越接近真實積分值。I、(6)使用均值估計法程序：function p=shell6(a1,b1,a2,b2,n)z=0;x=unifrnd(a1,b1,1,n);y=unifrnd(a2,b2,1,n);for i=1:n if y(i)<=sin(x(i); u=exp(-x(i)2-y(

9、i)2); z=z+u; endendp=(b1-a1)*(b2-a2)*z/n;運行結(jié)果：>> p=shell6(0,pi,0,1,100)p = 0.4368>> p=shell6(0,pi,0,1,1000)p = 0.3378>> p=shell6(0,pi,0,1,10000)p = 0.3674>> p=shell6(0,pi,0,1,100000)p = 0.3610結(jié)果分析：改變了四次隨機點數(shù)，結(jié)果都趨近于0.36，說明積分值約等于0.36，而且點數(shù)越多，值越接近真實積分值。II、(1)使用蒙特卡羅法分析：將x在它被允許的范圍內(nèi)生

10、成多個隨機的數(shù)值，利用max函數(shù)可以近似地求出結(jié)果。然后做出圖像，進行結(jié)果的比較。程序：function f81(n)x=unifrnd(-2*pi,2*pi,1,n);y=(1-x.2).*sin(3*x);max(y)x=-2*pi:0.001:2*pi;y=(1-x.2).*sin(3*x);plot(x,y)xlabel('x');ylabel('y');運行結(jié)果：>> f81(1000)ans = 32.3293>> f81(10000)ans = 32.4002>> f81(100000)ans = 32.4006

11、做出函數(shù)的圖像，并且標(biāo)出最高點的值結(jié)果分析：可以看到，蒙特卡羅法求出的最大值接近于32.4，而從圖中可以看出最大值是32.33，求出的結(jié)果比較符合。II、(2)使用均值估計法分析：由于x1=x2+10,所以可以消元，使其變?yōu)閮蓚€自變量x2和x3。x2,x3在它們被允許的范圍內(nèi)生成多個隨機的數(shù)值，利用max函數(shù)可以近似地求出結(jié)果。然后做出圖像，進行結(jié)果的比較。程序：function f82(n)x2=unifrnd(10,20,1,n);x1=10+x2;x3=unifrnd(-10,20,1,n);for i=1:nif -x1(i)+2*x2(i)+2*x3(i)>=0 if x1(i

12、)+2*x2(i)+2*x3(i)<=72 y(i)=(x1(i)*(x2(i)*(x3(i); endendendmax(y)x2=10:0.1:20;x3=-5:21/100:16;X,Y=meshgrid(x2,x3);err1 = X+2*Y<10;err2 = 3*X+2*Y>62;X(err1) = nan;Y(err2) = nan;Z=X.*Y.*(X+10);surf(X,Y,Z)運行結(jié)果：>> f82(1000)ans = 3.3889e+03>> f82(10000)ans = 3.4357e+03>> f82(100

13、)ans = 3.3726e+03>> f82(100000)ans = 3.4441e+03 結(jié)果分析：可以看到，蒙特卡羅法求出的最大值接近于3400，而從圖中可以看出最大值是3437，求出的結(jié)果比較符合。II、(3)使用蒙特卡羅法分析：x,y在它們被允許的范圍內(nèi)生成多個隨機的數(shù)值，利用max函數(shù)可以近似地求出結(jié)果。然后做出圖像，進行結(jié)果的比較。程序：function f83(n)x=unifrnd(-1.5,1.5,1,n);y=unifrnd(-1.5,1.5,1,n);z=(x.2+2*(y.2)+x.*y).*exp(-x.2-y.2);max(z)x=-1.5:0.1:1.5;y=-1.5:0.1:1.5;X,Y=meshgrid(x,y);Z=(X.2+2*(Y.2)+X.*Y).*exp(-X.2-Y.2);surf(X,Y,Z)運行結(jié)果：>> f83(1000)ans = 0.8105>> f83(10000)ans =0.8117 作出函數(shù)圖，并且標(biāo)出