大学数学实验之蒙特卡洛方法

1、数学实验报告班级： 序号： 姓名： 1 问题描述I、用蒙特卡罗方法计算以下函数在区间上的积分，并改变随机点数目观察对结果的影响。（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, 0=y=sin(x);II、用蒙特卡罗法求解全局最优化及约束问题并通过图形做出评论，求下列函数

2、的最大值。（1） f(x)=(1-x.2).*sin(3*x), -2*pi=x=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) p=shell1(0,1,10000)p = 0.6922 p=shell1(0,1,100)p = 0.7001 p=shell1(0,1,500)p =0.6890 结果分析：改变了四次随机点数，结果都趋近于0.69，说明积分值约等于0.69，但是点数越多，值越接近。 I、(2)使用均值估计法程序：function p

3、=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;运行结果： 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,)p = -30.0058结果分析：改变了5次随机点数，结果变化较大，但是点数越多，值越接近真实积分值。所以积分值近似于-30。

4、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;运行结果： 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=shell3(0,2,)p = 2.9600结果分析：改变了四次随机点数，结果都趋近于2.95，说明积分值约等于2.95，而且点数越多，值越接近真实积分值。I、(4)使用均值估计法程序：fun

5、ction 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;运行结果： p=shell4(0,2,)p = 0.4783 p=shell4(0,2,10000)p = 0.4777 p=shell4(0,2,1000)p = 0.4765 p=shell4(0,2,100)p = 0.4432结果分析：改变了四次随机点数，结果都趋近于0.47，说明积分值约等于0.47，而且点数越多，值越接近真实积分值。I、(5)使用均值估计法程序：functi

6、on 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;运行结果： p=shell5(0,2*pi,100)p = 22.0140 p=shell5(0,2*pi,1000)p = 20.2718 p=shell5(0,2*pi,10000)p = 20.9394 p=shell5(0,2*pi,)p = 20.7968结果分析：改变了四次随机点数，结果都趋近于20.8，说明积分值约等于20.8，而且点数越多，值越接近真实积分值。I、(6)使用均值估计法程

7、序：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) 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,)p = 0.3610结果分析：改变了四次随机点数，结果都趋近于0.36，说明积分值约等于0.36，而且点数越多，值越接近真实积分值。II、(1)使用蒙特卡罗法分析

8、：将x在它被允许的范围内生成多个随机的数值，利用max函数可以近似地求出结果。然后做出图像，进行结果的比较。程序：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);运行结果： f81(1000)ans = 32.3293 f81(10000)ans = 32.4002 f81()ans = 32.4006做出函数的图像，并且标出最高点的值结果分析：可以看到，蒙特卡罗法求出的最大

9、值接近于32.4，而从图中可以看出最大值是32.33，求出的结果比较符合。II、(2)使用均值估计法分析：由于x1=x2+10,所以可以消元，使其变为两个自变量x2和x3。x2,x3在它们被允许的范围内生成多个随机的数值，利用max函数可以近似地求出结果。然后做出图像，进行结果的比较。程序：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)+2*x2(i)+2*x3(i)=72 y(i)=(x1(i)*(x2(i)*(

10、x3(i); endendendmax(y)x2=10:0.1:20;x3=-5:21/100:16;X,Y=meshgrid(x2,x3);err1 = X+2*Y62;X(err1) = nan;Y(err2) = nan;Z=X.*Y.*(X+10);surf(X,Y,Z)运行结果： f82(1000)ans = 3.3889e+03 f82(10000)ans = 3.4357e+03 f82(100)ans = 3.3726e+03 f82()ans = 3.4441e+03 结果分析：可以看到，蒙特卡罗法求出的最大值接近于3400，而从图中可以看出最大值是3437，求出的结果比较符

11、合。II、(3)使用蒙特卡罗法分析：x,y在它们被允许的范围内生成多个随机的数值，利用max函数可以近似地求出结果。然后做出图像，进行结果的比较。程序：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)运行结果： f83(1000)ans = 0.8105 f83(10000)ans =0.8117 作出函数图，并且标出最大值 结果分析：可以看到，蒙特卡罗法求出的最大值接近于0.81，而从图