Slide1:
实验六:蒙特卡罗方法实验 面积、体积计算问题
冰淇淋锥的体积计算
思考题与练习题
二维填充图绘制方法
Slide2: 准备知识： (a,b)均匀分布的随机数函数：unifrnd
(0,1)均匀分布的随机数函数：rand，类似实现：unifrnd(0,1)
rand用法：
r=rand(n) 产生n维矩阵
r=rand(m,n) 产生m×n的矩阵
r=rand 产生1个在(0,1)均匀分别的随机数；
unifrnd用法：
r=unifrnd(a,b)产生(a,b)上均匀分布的随机数
r = unifrnd(a,b,m,n) returns an m by n matrix.
例子：
unifrnd(0,1) %产生(0,1)均匀分布
unifrnd(-1,1) %产生(-1,1)均匀分布
unifrnd(2,10,5,5) %what ’s this?
[a,b]均匀分布模拟：a+(b-a)*rand或unifrnd(a,b)
Slide3: rand 产生一个0到 1
之间均匀随机数
rand(m,n) 产生m×n个0到 1之间均匀随机数 X=rand(10000,1);
hist(X) 蒙特卡罗方法——利用随机试验做近似计算 10000个随机数较均匀地分布在各个小区间上，
随机变量X落入小区间的概率仅与小区间长度有关,而与小区间位置无关 rand功能的一个直观说明图
Slide4: 例1 计算两条抛物线y = x2, x = y2所围面积. 在正方形[0,1]×[0,1]区域投入2000个均匀随机点则随机点落入抛物线所围区域的概率p为所求面积与正方形面积之比 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 S0 S p=S/S0 ；
投N个点，m个落入S，则频率f=m/N；
频率近似概率：f≈p
→S =pS0≈S0m/N
Slide5: 定积分数值计算方法quad()的使用格式 quad('F',a,b)
返回被积函数F(X)从a 到 b的定积分值,F是被积函数名
构成的字符串. 例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
Slide6: 例3. 计算下面两条曲线所围区域面积 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;%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
Slide7: 实验:冰淇淋锥的体积计算 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
Slide8: 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>=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),hold on
mesh(x,y,z2)
colormap([0 0 1])
axis off
view(0,-18) ←参考代码；
要求：改用for语句实现体积估算
Slide9: 思考题与练习题 1. 如何用极坐标变换处理重积分 r in [0,1],seta in [0,2pi]
Slide10: 二维多边形填充图 fill() 使用格式 fill(x,y,c)
用c所指定的颜色对多边形填充.其中,多边形的顶点由x,y确定. c 指定颜色不允许省略;x,y确定的点必须形成封闭的多边形.
Slide11: 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')
axis off
辅助内容：: 辅助内容：
Slide13: 例4. 计算两个半径为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<=1&x.^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
与quad相关命令: 与quad相关命令 相关命令:dblquad() —— 重积分计算