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

《数学实验》报告

班级：序号：姓名：

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） 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)使用均值估计法

程序：

function p=shell1(a,b,n)

z=0;

x=unifrnd(a,b,1,n);

fori=1:n

u=(x(i)+1)^(-1);

z=z+u;

end

p=(b-a)*z/n;

运行结果：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

结果分析：改变了四次随机点数，结果都趋近于0.69，说明积分值约等于0.69，但是点数越多，值越接近。

I、(2)使用均值估计法

程序：

function p=shell2(a,b,n)

z=0;

x=unifrnd(a,b,1,n);

fori=1:n

u=(exp(3*x(i)))*sin(2*x(i));

z=z+u;

end

p=(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,100000)

p =

-30.0058

结果分析：

改变了5次随机点数，结果变化较大，但是点数越多，值越接近真实积分值。所以积分值近似于-30。

I、(3)使用均值估计法

程序：

function p=shell3(a,b,n)

z=0;

x=unifrnd(a,b,1,n);

fori=1:n

u=(1+x(i)^2)^0.5;

z=z+u;

end

p=(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,100000)

p =

2.9600

结果分析：改变了四次随机点数，结果都趋近于2.95，说明积分值约等于2.95，而且点数越多，值越接近真实积分值。

I、(4)使用均值估计法

程序：

function p=shell4(a,b,n)

z=0;

x=unifrnd(a,b,1,n);

fori=1:n

u=(1/(2*pi)^0.5)*exp(-x(i)^2/2);

z=z+u;

end

p=(b-a)*z/n;

运行结果：

>> p=shell4(0,2,100000)

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)使用均值估计法

程序：

function p=shell5(a,b,n)

z=0;

x=unifrnd(a,b,1,n);

fori=1:n

u=exp(x(i)/2)*(sin(x(i)))^2;

z=z+u;

end

p=(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,100000)

p =

20.7968

结果分析：改变了四次随机点数，结果都趋近于20.8，说明积分值约等于20.8，而且点数越多，值越接近真实积分值。

I、(6)使用均值估计法

程序：

function p=shell6(a1,b1,a2,b2,n)

z=0;

x=unifrnd(a1,b1,1,n);

y=unifrnd(a2,b2,1,n);

fori=1:n

if y(i)<=sin(x(i));

u=exp(-x(i)^2-y(i)^2);

z=z+u;

end

end

p=(b1-a1)*(b2-a2)*z/n;

运行结果：

>> 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

结果分析：改变了四次随机点数，结果都趋近于0.36，说明积分值约等于0.36，而且点数越多，值越接近真实积分值。

II、(1)使用蒙特卡罗法

分析：将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(100000)

ans =

32.4006

做出函数的图像，并且标出最高点的值

结果分析：可以看到，蒙特卡罗法求出的最大值接近于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);

fori=1:n

if -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))*(x3(i));

end

end

end

max(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)

运行结果：

>>f82(1000)

ans =

3.3889e+03

>> f82(10000)

ans =

3.4357e+03

>> f82(100)

ans =

3.3726e+03

>> f82(100000)

ans =

3.4441e+03

结果分析：可以看 …… 此处隐藏：2274字，全部文档内容请下载后查看。喜欢就下载吧 ……