搜索
《数学实验》报告班级：序号：姓名：1．问题描述I、用蒙特卡罗方法计算以下函数在区间上的积分，并改变随机点数目观察对结果的影响。（1）y=1/(1+x),0=x=1;（2）y=(exp(3*x))*sin(2*x),0=x=2;（3）y=(1+x^2)^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(-x^2-y^2)0=x=pi,0=y=sin(x);II、用蒙特卡罗法求解全局最优化及约束问题并通过图形做出评论，求下列函数的最大值。（1）f(x)=(1-x.^2).*sin(3*x),-2*pi=x=2*pi;（2）maxf(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)使用均值估计法程序：functionp=shell1(a,b,n)z=0;x=unifrnd(a,b,1,n);fori=1:nu=(x(i)+1)^(-1);z=z+u;endp=(b-a)*z/n;运行结果：p=shell1(0,1,1000)p=0.6975p=shell1(0,1,10000)p=0.6922p=shell1(0,1,100)p=0.7001p=shell1(0,1,500)p=0.6890结果分析：改变了四次随机点数，结果都趋近于0.69，说明积分值约等于0.69，但是点数越多，值越接近。I、(2)使用均值估计法程序：functionp=shell2(a,b,n)z=0;x=unifrnd(a,b,1,n);fori=1:nu=(exp(3*x(i)))*sin(2*x(i));z=z+u;endp=(b-a)*z/n;运行结果：p=shell2(0,2,1000)p=-24.4911p=shell2(0,2,100)p=-43.8720p=shell2(0,2,10000)p=-30.8699p=shell2(0,2,500)p=-23.2955p=shell2(0,2,100000)p=-30.0058结果分析：改变了5次随机点数，结果变化较大，但是点数越多，值越接近真实积分值。所以积分值近似于-30。I、(3)使用均值估计法程序：functionp=shell3(a,b,n)z=0;x=unifrnd(a,b,1,n);fori=1:nu=(1+x(i)^2)^0.5;z=z+u;endp=(b-a)*z/n;运行结果：p=shell3(0,2,100)p=2.9293p=shell3(0,2,1000)p=2.9516p=shell3(0,2,10000)p=2.9512p=shell3(0,2,100000)p=2.9600结果分析：改变了四次随机点数，结果都趋近于2.95，说明积分值约等于2.95，而且点数越多，值越接近真实积分值。I、(4)使用均值估计法程序：functionp=shell4(a,b,n)z=0;x=unifrnd(a,b,1,n);fori=1:nu=(1/(2*pi)^0.5)*exp(-x(i)^2/2);z=z+u;endp=(b-a)*z/n;运行结果：p=shell4(0,2,100000)p=0.4783p=shell4(0,2,10000)p=0.4777p=shell4(0,2,1000)p=0.4765p=shell4(0,2,100)p=0.4432结果分析：改变了四次随机点数，结果都趋近于0.47，说明积分值约等于0.47，而且点数越多，值越接近真实积分值。I、(5)使用均值估计法程序：functionp=shell5(a,b,n)z=0;x=unifrnd(a,b,1,n);fori=1:nu=exp(x(i)/2)*(sin(x(i)))^2;z=z+u;endp=(b-a)*z/n;运行结果：p=shell5(0,2*pi,100)p=22.0140p=shell5(0,2*pi,1000)p=20.2718p=shell5(0,2*pi,10000)p=20.9394p=shell5(0,2*pi,100000)p=20.7968结果分析：改变了四次随机点数，结果都趋近于20.8，说明积分值约等于20.8，而且点数越多，值越接近真实积分值。I、(6)使用均值估计法程序：functionp=shell6(a1,b1,a2,b2,n)z=0;x=unifrnd(a1,b1,1,n);y=unifrnd(a2,b2,1,n);fori=1:nify(i)=sin(x(i));u=exp(-x(i)^2-y(i)^2);z=z+u;endendp=(b1-a1)*(b2-a2)*z/n;运行结果：p=shell6(0,pi,0,1,100)p=0.4368p=shell6(0,pi,0,1,1000)p=0.3378p=shell6(0,pi,0,1,10000)p=0.3674p=shell6(0,pi,0,1,100000)p=0.3610结果分析：改变了四次随机点数，结果都趋近于0.36，说明积分值约等于0.36，而且点数越多，值越接近真实积分值。II、(1)使用蒙特卡罗法分析：将x在它被允许的范围内生成多个随机的数值，利用max函数可以近似地求出结果。然后做出图像，进行结果的比较。程序：functionf81(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.3293f81(10000)ans=32.4002f81(100000)ans=32.4006做出函数的图像，并且标出最高点的值结果分析：可以看到，蒙特卡罗法求出的最大值接近于32.4，而从图中可以看出最大值是32.33，求出的结果比较符合。II、(2)使用均值估计法分析：由于x1=x2+10,所以可以消元，使其变为两个自变量x2和x3。x2,x3在它们被允许的范围内生成多个随机的数值，利用max函数可以近似地求出结果。然后做出图像，进行结果的比较。程序：functionf82(n)x2=unifrnd(10,20,1,n);x1=10+x2;x3=unifrnd(-10,20,1,n);fori=1:nif-x1(i)+2*x2(i)+2*x3(i)=0ifx1(i)+2*x2(i)+2*x3(i)=72y(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*Y10;err2=3*X+2*Y62;X(err1)=nan;Y(err2)=nan;Z=X.*Y.*(X+10);surf(X,Y,Z)运行结果：f82(1000)ans=3.3889e+03f82(10000)ans=3.4357e+03f82(100)ans=3.3726e+03f82(100000)ans=3.4441e+03结果分析：可以看到，蒙特卡罗法求出的最大值接近于3400，而从图中可以看出最大值是3437，求出的结果比较符合。II、(3)使用蒙特卡罗法分析：x,y在它们被允许的范围内生成多个随机的数值，利用max函数可以近似地求出结果。然后做出图像，进行结果的比较。程序：functionf83(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.8105f83(10000)ans=0.8117作出函数图，并且标出最大值结果分析：可以看到，蒙特卡罗法求出的最大值接近于0.81，而从图中可以看出最大值是0.8025，求出的结果比较符合。3.实验总结和实验感悟这次蒙特卡洛法令我印象比较深刻，特别是可以利用多次模拟实验的方法来求圆周率，这是我以前没有接触过的。蒙特卡洛法可以理解成一种思想，就是多次随机的实验来求近似值。不过这种方法比较适合电脑模拟，模拟次数足够高才可以保证误差不过大，而且某些可以直接求解的问题并不需要用蒙特卡罗法来做。