实验六蒙特卡罗方法

实验六蒙特卡罗方法第一页，共十四页，编辑于2023年，星期二rand

产生一个0到

1之间均匀随机数rand(m,n)

产生m×n个0到

1之间均匀随机数X=rand(10000,1);hist(X)蒙特卡罗方法——利用随机试验做近似计算10000个随机数较均匀地分布在各个小区间上，随机变量X落入小区间的概率仅与小区间长度有关,而与小区间位置无关unifrnd(a,b)

产生一个a到

b之间均匀随机数X=unifrnd(2,10)Y=2+(10-2)*rand第二页，共十四页，编辑于2023年，星期二例1

计算两条抛物线y=x2,x=y2所围面积.在正方形[0,1]×[0,1]区域投入2000个均匀随机点则随机点落入抛物线所围区域的概率为所求面积与正方形面积之比functionS=area1(N)ifnargin==0,N=2000;endX=rand(N,1);Y=rand(N,1);II=find(Y<=sqrt(X)&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第三页，共十四页，编辑于2023年，星期二定积分数值计算方法quad()的使用格式quad('F',a,b)返回被积函数F(X)从a

到

b的定积分值,F是被积函数名构成的字符串.相关命令:dblquad()——重积分计算例2.计算定积分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')symsu,S0=int(sqrt(u)-u^2,0,1)S=0.3327S0=1/3第四页，共十四页，编辑于2023年，星期二例3.计算下面两条曲线所围区域面积functionS=area2(N)ifnargin==0,N=2000;endX=2*rand(N,1)-1;Y=2*rand(N,1);II=find(Y<=1+sqrt(1-X.^2)&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);x3=-1:-0.01:0;y3=-x3;fill([x1,x2,x3],[y1,y2,y3],'c')S=2.5460第五页，共十四页，编辑于2023年，星期二例4.计算两个半径为1的直交圆柱面所围成体积x2+y2=1,x2+z2=1functionV=mlab4(N)ifnargin==0,N=2000;endP=rand(N,3);x=P(:,1);y=P(:,2);z=P(:,3);II=find(x.^2+y.^2<=1&x.^2+z.^2<=1);n=length(II);V=8*n/Nh=2*pi/100;t=0:h:pi;r=0:0.05:1;x=r'*cos(t);y=r'*sin(t);zz=sqrt(1-x.^2);meshz(x,y,zz)colormap([001])axisoffview(119,34)V=16/3D:第六页，共十四页，编辑于2023年，星期二实验:冰淇淋锥的体积计算&x=2*rand-1产生–1到1之间的随机数y=2*rand-1产生–1到1之间的随机数z=2*rand;产生0到2之间的随机数冰淇淋锥含于体积=8的六面体22由于rand产生0到1之间的随机数,所以N个点均匀分布于六面体中,锥体中占有m个,则锥体与六面体体积之比近似为m:N第七页，共十四页，编辑于2023年，星期二functionV=icecream(N)ifnargin==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&Z<=1+sqrt(1-R2));V=8*length(II)/N;h=2*pi/100;t=0:h:2*pi;r=0:.05:1;x=r'*cos(t);y=r'*sin(t);z1=sqrt(x.^2+y.^2);z2=1+sqrt(1+eps-x.^2-y.^2);mesh(x,y,z1),holdonmesh(x,y,z2)colormap([001])axisoffview(0,-18)第八页，共十四页，编辑于2023年，星期二思考题与练习题1.如何用极坐标变换处理重积分二维多边形填充图

fill()使用格式fill(x,y,c)用c所指定的颜色对多边形填充.其中,多边形的顶点由x,y确定.c指定颜色不允许省略;x,y确定的点必须形成封闭的多边形.第九页，共十四页，编辑于2023年，星期二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+sin(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')axisoff大圆小圆第十页，共十四页，编辑于2023年，星期二38/wlxt/ncourse/mathlab/web/shiyan//

第十一页，共十四页，编辑于2023年，星期二

functionV=icecream(N)ifnargin==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&Z<=1+sqrt(1-R2));V=8*M/N;symsrtx=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)第十二页，共十四页，编辑于2023年，星期二

functionV=icecream(N)ifnargin==0,N=10000;endm=0;fori=1:NX=-1+2*rand;Y=-1+2*rand;Z=2*rand;R2=X^2+Y^2;R=sqrt(R2);ifZ>=R&Z<=(1+sqrt(1-R2))m=m+1;endendV=8*m/N第十三页，共十四页，编辑于