《实验六蒙特卡罗方法.ppt》由会员分享，可在线阅读，更多相关《实验六蒙特卡罗方法.ppt（14页珍藏版）》请在三一办公上搜索。
1、,实验六:蒙特卡罗方法实验,面积、体积计算问题冰淇淋锥的体积计算思考题与练习题二维填充图绘制方法,rand 产生一个0到 1之间均匀随机数rand(m,n)产生mn个0到 1之间均匀随机数,X=rand(10000,1);hist(X),蒙特卡罗方法利用随机试验做近似计算,10000个随机数较均匀地分布在各个小区间上，随机变量X落入小区间的概率仅与小区间长度有关,而与小区间位置无关,unifrnd(a,b)产生一个a到 b之间均匀随机数,X=unifrnd(2,10)Y=2+(10-2)*rand,例1 计算两条抛物线y=x2,x=y2所围面积.,在正方形0,10,1区域投入2000个均匀随机
2、点则随机点落入抛物线所围区域的概率为所求面积与正方形面积之比,function S=area1(N)if nargin=0,N=2000;endX=rand(N,1);Y=rand(N,1);II=find(Y=X.2);m=length(II);S=m/N;x1=0:0.01:1;x2=1:-0.01:0;y1=sqrt(x1);y2=x2.2;fill(x1,x2,y1,y2,c),S=0.3333,定积分数值计算方法quad()的使用格式,quad(F,a,b)返回被积函数F(X)从a 到 b的定积分值,F是被积函数名构成的字符串.相关命令:dblquad()重积分计算,例2.计算定积分
3、,fun=inline(sqrt(x)-x.2);S=quad(fun,0.01,1)t=0:0.01:1;y=fun(t);fill(0,t,0,y,c)syms u,S0=int(sqrt(u)-u2,0,1),S=0.3327S0=1/3,例3.计算下面两条曲线所围区域面积,function S=area2(N)if nargin=0,N=2000;endX=2*rand(N,1)-1;Y=2*rand(N,1);II=find(Y=abs(X);m=length(II);S=4*m/N;x1=0:0.01:1;y1=x1;x2=1:-0.01:-1;y2=1+sqrt(1-x2.2);
4、x3=-1:-0.01:0;y3=-x3;fill(x1,x2,x3,y1,y2,y3,c),S=2.5460,例4.计算两个半径为1的直交圆柱面所围成体积,x2+y2=1,x2+z2=1,function V=mlab4(N)if nargin=0,N=2000;endP=rand(N,3);x=P(:,1);y=P(:,2);z=P(:,3);II=find(x.2+y.2=1meshz(x,y,zz)colormap(0 0 1)axis offview(119,34),V=16/3,D:,实验:冰淇淋锥的体积计算,&,x=2*rand-1产生 1到1之间的随机数y=2*rand-1产生
5、 1到1之间的随机数z=2*rand;产生0到2之间的随机数,冰淇淋锥含于体积=8 的六面体,2,2,由于rand 产生0 到1之间的随机数,所以,N个点均匀分布于六面体中,锥体中占有m个,则锥体与六面体体积之比近似为 m:N,function V=icecream(N)if nargin=0,N=10000;endP=rand(N,3);X=2*P(:,1)-1;Y=2*P(:,2)-1;Z=2*P(:,3);R2=X.2+Y.2;R=sqrt(R2);II=find(Z=R,h=2*pi/100;t=0:h:2*pi;r=0:.05:1;x=r*cos(t);y=r*sin(t);z1=s
6、qrt(x.2+y.2);z2=1+sqrt(1+eps-x.2-y.2);mesh(x,y,z1),hold onmesh(x,y,z2)colormap(0 0 1)axis offview(0,-18),思考题与练习题,1.如何用极坐标变换处理重积分,二维多边形填充图 fill()使用格式,fill(x,y,c)用c所指定的颜色对多边形填充.其中,多边形的顶点由x,y确定.c 指定颜色不允许省略;x,y确定的点必须形成封闭的多边形.,2.圆 x2+(y 2)2=4内挖去小圆x2+(y 1)2=1后图形,h=2*pi/100;t=-pi/2:h:3*pi/2;x=cos(t);y=1+si
7、n(t);X=2*cos(t);Y=2+2*sin(t);n=length(x);x=x(n:-1:1);y=y(n:-1:1);fill(x,X,y,Y,c)axis off,大圆,小圆,/,function V=icecream(N)if nargin=0,N=10000;endP=rand(N,3);X=2*P(:,1)-1;Y=2*P(:,2)-1;Z=2*P(:,3);R2=X.2+Y.2;R=sqrt(R2);M=sum(Z=R,syms r tx=r*cos(t);y=r*sin(t);z1=sqrt(x.2+y.2);z2=1+sqrt(1+eps-x.2-y.2);f=r*(z2-z1)f=int(f,t,0,2*pi);f=int(f,r,0,1),function V=icecream(N)if nargin=0,N=10000;endm=0;for i=1:NX=-1+2*rand;Y=-1+2*rand;Z=2*rand;R2=X2+Y2;R=sqrt(R2);if Z=RendendV=8*m/N,function f=ff(t,r)x=r*cos(t);y=r*sin(t);z1=sqrt(x.2+y.2);z2=1+sqrt(1+eps-x.2-y.2);f=r*(z2-z1);,f=dblquad(ff,0,2*pi,0,1),
