实验六 蒙特卡罗方法及程序.ppt

*/10 ? ? ? ? 实验六:蒙特卡罗方法实验 二维填充图绘制方法 面积、体积计算问题 冰淇淋锥的体积计算 思考题与练习题 fill(x,y,c) 用c所指定的颜色对多边形填充.其中,多边形的顶点由x,y确定. c 指定颜色不允许省略;x,y确定的点必须形成封闭的多边形. 二维多边形填充图 fill() 使用格式 一.二维填充图绘制方法 clf;n=10; dt=2*pi/n;t=0:dt:2*pi; t=[t,t(1)]; x=sin(t);y=cos(t); fill(x,y,c);axis off rand 产生一个0到 1 之间均匀随机数 rand(m,n) 产生m×n个0到 1之间均匀随机数 X=rand(10000,1); hist(X) (1)蒙特卡罗方法——利用随机试验做近似计算 10000个随机数较均匀地分布在各个小区间上， 随机变量X落入小区间的概率仅与小区间长度有关,而与小区间位置无关 二. 面积、体积计算问题 例1. 计算下面两条曲线所围区域面积 function S=area2(N) if nargin==0,N=2000;end X=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 例2. 计算两个半径为1的直交圆柱面所围成体积 x2 + y2 = 1 , x2 + z2 = 1 function V=mlab4(N) if nargin==0,N=2000;end P=rand(N,3); x=P(:,1);y=P(:,2);z=P(:,3); II=find(x.^2+y.^2=1x.^2+z.^2=1); n=length(II);V=8*n/N h=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([0 0 1]) axis off view(119,34) V = 16/3 D: (2)定积分数值计算方法:quad()的使用格式 quad(F,a,b) 返回被积函数F(X)从a 到 b的定积分值,F是被积函 数名构成的字符串.相关命令:dblquad() —重积分计算 例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)-u^2,0,1) S = 0.3327 S0 = 1/3 三.实验:冰淇淋锥的体积计算 x=2*rand-1产生– 1到1之间的随机数 y=2*rand-1产生– 1到1之间的随机数 z=2*rand;产生0到2之间的随机数 冰淇淋锥含于体积 = 8 的六面体 2 2 由于rand 产生0 到1之间的随机数,所以 N个点均匀分布于六面体中,锥体中占有m个,则锥体与六面体体积之比近似为 m : N function data=icecream(L) if nargin==0,L=7;end N=10000; for k=1:L P=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=Rz=1+sqrt(1-R2)); m=length(II); q(k)=8*m/N; end data=[q; q-pi]; 实验任务:蒙特卡罗方法计算冰淇淋锥的体积 function V=icecream(N) if nargin==0,N=10000;end P=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=RZ=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