實驗六蒙特卡羅方法第一頁，共十四頁，2022年，8月28日rand

產生一個0到

1之間均勻隨機數rand(m,n)

產生m×n個0到

1之間均勻隨機數X=rand(10000,1);hist(X)蒙特卡羅方法——利用隨機試驗做近似計算10000個隨機數較均勻地分布在各個小區(qū)間上，隨機變量X落入小區(qū)間的概率僅與小區(qū)間長度有關,而與小區(qū)間位置無關unifrnd(a,b)

產生一個a到

b之間均勻隨機數X=unifrnd(2,10)Y=2+(10-2)*rand第二頁，共十四頁，2022年，8月28日例1

計算兩條拋物線y=x2,x=y2所圍面積.在正方形[0,1]×[0,1]區(qū)域投入2000個均勻隨機點則隨機點落入拋物線所圍區(qū)域的概率為所求面積與正方形面積之比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第三頁，共十四頁，2022年，8月28日定積分數值計算方法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第四頁，共十四頁，2022年，8月28日例3.計算下面兩條曲線所圍區(qū)域面積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第五頁，共十四頁，2022年，8月28日例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:第六頁，共十四頁，2022年，8月28日實驗:冰淇淋錐的體積計算&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第七頁，共十四頁，2022年，8月28日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)第八頁，共十四頁，2022年，8月28日思考題與練習題1.如何用極坐標變換處理重積分二維多邊形填充圖

fill()使用格式fill(x,y,c)用c所指定的顏色對多邊形填充.其中,多邊形的頂點由x,y確定.c指定顏色不允許省略;x,y確定的點必須形成封閉的多邊形.第九頁，共十四頁，2022年，8月28日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大圓小圓第十頁，共十四頁，2022年，8月28日/

第十一頁，共十四頁，2022年，8月28日

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)第十二頁，共十四頁，2022年，8月28日

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第十三頁，共十四頁，202