实验六 蒙特卡罗方法.pptVIP

*/10 ? ? ? 实验六:蒙特卡罗方法实验 面积、体积计算问题 冰淇淋锥的体积计算 思考题与练习题 rand 产生一个0到 1 之间均匀随机数 rand(m,n) 产生m×n个0到 1之间均匀随机数 X=rand(10000,1); hist(X) 蒙特卡罗方法——通过随机变量的统计试验求近似解 10000个随机数较均匀地分布在各个小区间上， 随机变量X落入小区间的概率仅与小区间长度有关,而与小区间位置无关 例1 计算两条抛物线y = x2, x = y2所围面积. 在正方形[0,1]×[0,1]区域投入2000个均匀随机点 则随机点落入抛物线所围区域的概率为所求面 积与正方形面积之比 function S=area1(N) if nargin==0,N=2000;end X=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 定积分数值计算方法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) syms u, S0=int(sqrt(u)-u^2,0,1) S = 0.3327 S0 = 1/3 例3. 计算两个半径为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); m=length(II); V=8*m/N V = 16/3 D: 实验:冰淇淋锥的体积计算 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(x,y,z1),hold on mesh(x,y,z2) colormap([0 0 1]) axis off view(0,-18)