大学数学实验之蒙特卡洛方法
《数学实验》报告
班级:序号:姓名:
1. 问题描述
I、用蒙特卡罗方法计算以下函数在区间上的积分,并改变随机点数目观察对结果的影响。
(1)y=1/(1+x), 0=
(2)y= (exp(3*x))*sin(2*x), 0=
(3)y=(1+x^2)^0.5, 0=
(4)y=(1/(2*pi)^0.5)*exp(-x(i)^2/2),0=
(5)y=exp(x(i)/2)*(sin(x(i)))^2, 0=
(6)f(x,y)=exp(-x^2-y^2)0=
II、用蒙特卡罗法求解全局最优化及约束问题并通过图形做出评论,求下列函数的最大值。
(1) f(x)=(1-x.^2).*sin(3*x),-2*pi=
(2) max
f(x)=x1*x2*x3,s.t.:-x1+2x2+2x3>=0,x1+2x2+2x3
=x2
(3) f(x,y)=(X.^2+2*(Y.^2)+X.*Y).*exp(-X.^2-Y.^2),
abs(x)
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)
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)
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
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
结果分析:可以看到,蒙特卡罗法求出的最大值接近于3400,而从图中可以看出最大值是3437,求出的结果比较符合。
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,而从图中可以看出最大值是0.8025,求出的结果比较符合。
3.实验总结和实验感悟
这次蒙特卡洛法令我印象比较深刻,特别是可以利用多次模拟实验的方法来求圆周率,这是我以前没有接触过的。蒙特卡洛法可以理解成一种思想,就是多次随机的实验来求近似值。不过这种方法比较适合电脑模拟,模拟次数足够高才可以保证误差不过大,而且某些可以直接求解的问题并不需要用蒙特卡罗法来做。