电子科技大学数学实验实验五 Monte Carlo方法计算Viviani体体积
合集下载
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
10/12
思考题与练习题
1.蒙特卡罗方法计算面积和计算体积方法有何差异？ 2.说明蒙特卡罗方法计算Viviani的体积误差与哪些因 素有关。 3.概率论中的贝努里大数定律在本次实验中体现如何
4.叙述概率论中著名独立同分布大数定律，并以这一 大数定律解释实验的数据变换规律。 5. 说明L次蒙特卡罗实验平均值计算Viviani的体积误 差与实验次数之间关系
theta=(-180:10:1800)*pi/180; r=0:0.1:1; X=r'*(1+cos(theta));Y=r'*sin(theta); Z=sqrt(4-X.^2-Y.^2); meshz(X,Y,Z) axis off colormap([0 0 0]) view(130,30)
8/12
实验任务一：记录L次实验的实验数据及误差 序号 数据 误差 1 2 3 4 5 6 7
实验任务二：修改实验程序MonteC计算L次实验数 据均值及均值误差( mean 计算平均值 ) L 8 16 32 64 128 256
均值
误差
9/12
Viviani体绘制程序
圆柱面在 x y 平面的投影用极坐 标表示为 x 0.5R(1 cos ), y 0.5R sin , [0, 2 ]
N个点均匀分布于六面体中, Viviani体中 占有m个,则锥体与六面体体积之比近似为 m:N
V 8
m N
5/12
求球体
2 2 x 2 y 2 z 2 R 2 被圆柱面 x y Rx
所割下部分的体积 (称为维维安尼 (Viviani)体 ). 解 由所求立体的对称性，只要求出在第 一卦限内的部分体积,再乘以4,即得所求立体的体
实验五:蒙特卡罗方法实验
面积、体积计算问题
Viviani体的体积计算
思考题与练习题
1/12
蒙特卡罗方法——随机投点试验求近似解 引例. 给定曲线y =2 – x2 和曲线y3 = x2，曲线的交点为： P1( – 1，1 )、P2( 1，1 )。曲线围成平面有限区域，用 蒙特卡罗方法计算区域面积。 P=rand(10000,2); x=2*P(:,1)-1;
11/12
5.下面程序绘出二维图形填 充图(右图)。分析每条语句功 能,给程序中语句写注记 x1=-1:0.1:1;
y1=x1.^2.^(1/3); %函数1 x2=1:-0.1:-1; y2=2-x2.^2; %函数2 fill([x1,x2],[y1,y2],‘g’)%填 充函数1曲线和函数2曲线之 间的区域
–1≤ x ≤1，–1≤ y ≤1， 0≤ z ≤1，0≤ u ≤2 P=rand(10000,4); x=-1+2*P(:,1); y=-1+2*P(:,2); z=P(:,3);u=2*P(:,4); II=find(z>sqrt(x.^2+y.^2)&z<=1&u<=x.^2+y.^2+z.^2); M=length(II); V=8*M/10000
D的边界曲线交点为：(1，-1)，(4，2)，被积函数在求 积区域内的最大值为16。积分值是三维体积，该三维 图形位于立方体区域 0≤ x ≤4，–1≤ y ≤2，0 ≤ z ≤16 内，立方体区域的体积为192。 data=rand(10000,3); x=4*data(:,1); y=-1+3*data(:,2); z=16*data(:,3); II=find(x>=y.^2&x<=y+2&z<=x.*(y.^2)); M=length(II); V=192*M/10000
4/12
实验:蒙特卡罗方法计算Viviani
体体积
x 2 y 2 z 2 R 2 & x 2 y 2 Rx， R 2
含于体积 = 8 的六面体 {( x, y , z ) | 0 x 2, 1 y 1, 0 z 2} 由于rand 产生0 到1之间的随机数,所以 x=2*rand产生0到2之间的随机数 y=2*rand-1产生-1到1之间的随机数 z=2*rand产生0到2之间的随机数
axis off%去掉坐标系
12/12
xy 2
xy dxdy
2
1dzdxdy
D
D
0
3/12
( x
例5.15 用蒙特卡罗方法计算
其中,积分区域是由
2
y z ) dxdydz
x2 y2 z2
z x y和 z = 1 所围成。 0
1 dudxdydz
2
被积函数在求积区域上的最大 值为2。所以有四维超立方体
y=2*P(:,2);
II=find(y<=2-x.^2&y.^3>=x.^2);
M=length(II);
S=4*M/10000 plot(x(II),y(II),'g.')
1 1
S=
2
2.1136
3 2
S ( 2 x x ) dx
32 15
2/12
例5.14 计算
D
xy 2 dxdy 其中D为y= x – 2与y2 = x 所围
2 2
(x
R 2
) y
2 2
R2 4
V
R
0
Rx x 2 Rx x 2
2 R x y dy dx R . 3 2 3
2 2 2
4
3
前页 后页 返回
实验参考程序 蒙特卡罗方法计算体积 syms x y; f=sqrt(4-x^2-y^2); y1=-sqrt(2*x-x^2);y2=sqrt(2*x-x^2); S1=int(f,y,y1,y2);S2=int(S1,x,0,2) V=double(S2) P=rand(10000,3); X=2*P(:,1);Y=2*P(:,2)-1;Z=2*P(:,3); II=find((X1).^2+Y.^2<=1&Z<=sqrt(4-X.^2Y.^2)); V1=8*length(II)/10000 Error=V-V1
z
y
R
x
前页 后页 返回
积. 在第一卦限内的立体是一个曲顶柱体, 其底为 xy 平面内由 y 0 和 x 2 y 2 Rx 所确定的区域 D (图21-32),而曲顶的方程为 z
R 2 x 2 y 2 . 所以
V 4 R 2 x 2 y 2 d ,
D
其中 D ( x , y ) y 0 , x y Rx .
思考题与练习题
1.蒙特卡罗方法计算面积和计算体积方法有何差异？ 2.说明蒙特卡罗方法计算Viviani的体积误差与哪些因 素有关。 3.概率论中的贝努里大数定律在本次实验中体现如何
4.叙述概率论中著名独立同分布大数定律，并以这一 大数定律解释实验的数据变换规律。 5. 说明L次蒙特卡罗实验平均值计算Viviani的体积误 差与实验次数之间关系
theta=(-180:10:1800)*pi/180; r=0:0.1:1; X=r'*(1+cos(theta));Y=r'*sin(theta); Z=sqrt(4-X.^2-Y.^2); meshz(X,Y,Z) axis off colormap([0 0 0]) view(130,30)
8/12
实验任务一：记录L次实验的实验数据及误差 序号 数据 误差 1 2 3 4 5 6 7
实验任务二：修改实验程序MonteC计算L次实验数 据均值及均值误差( mean 计算平均值 ) L 8 16 32 64 128 256
均值
误差
9/12
Viviani体绘制程序
圆柱面在 x y 平面的投影用极坐 标表示为 x 0.5R(1 cos ), y 0.5R sin , [0, 2 ]
N个点均匀分布于六面体中, Viviani体中 占有m个,则锥体与六面体体积之比近似为 m:N
V 8
m N
5/12
求球体
2 2 x 2 y 2 z 2 R 2 被圆柱面 x y Rx
所割下部分的体积 (称为维维安尼 (Viviani)体 ). 解 由所求立体的对称性，只要求出在第 一卦限内的部分体积,再乘以4,即得所求立体的体
实验五:蒙特卡罗方法实验
面积、体积计算问题
Viviani体的体积计算
思考题与练习题
1/12
蒙特卡罗方法——随机投点试验求近似解 引例. 给定曲线y =2 – x2 和曲线y3 = x2，曲线的交点为： P1( – 1，1 )、P2( 1，1 )。曲线围成平面有限区域，用 蒙特卡罗方法计算区域面积。 P=rand(10000,2); x=2*P(:,1)-1;
11/12
5.下面程序绘出二维图形填 充图(右图)。分析每条语句功 能,给程序中语句写注记 x1=-1:0.1:1;
y1=x1.^2.^(1/3); %函数1 x2=1:-0.1:-1; y2=2-x2.^2; %函数2 fill([x1,x2],[y1,y2],‘g’)%填 充函数1曲线和函数2曲线之 间的区域
–1≤ x ≤1，–1≤ y ≤1， 0≤ z ≤1，0≤ u ≤2 P=rand(10000,4); x=-1+2*P(:,1); y=-1+2*P(:,2); z=P(:,3);u=2*P(:,4); II=find(z>sqrt(x.^2+y.^2)&z<=1&u<=x.^2+y.^2+z.^2); M=length(II); V=8*M/10000
D的边界曲线交点为：(1，-1)，(4，2)，被积函数在求 积区域内的最大值为16。积分值是三维体积，该三维 图形位于立方体区域 0≤ x ≤4，–1≤ y ≤2，0 ≤ z ≤16 内，立方体区域的体积为192。 data=rand(10000,3); x=4*data(:,1); y=-1+3*data(:,2); z=16*data(:,3); II=find(x>=y.^2&x<=y+2&z<=x.*(y.^2)); M=length(II); V=192*M/10000
4/12
实验:蒙特卡罗方法计算Viviani
体体积
x 2 y 2 z 2 R 2 & x 2 y 2 Rx， R 2
含于体积 = 8 的六面体 {( x, y , z ) | 0 x 2, 1 y 1, 0 z 2} 由于rand 产生0 到1之间的随机数,所以 x=2*rand产生0到2之间的随机数 y=2*rand-1产生-1到1之间的随机数 z=2*rand产生0到2之间的随机数
axis off%去掉坐标系
12/12
xy 2
xy dxdy
2
1dzdxdy
D
D
0
3/12
( x
例5.15 用蒙特卡罗方法计算
其中,积分区域是由
2
y z ) dxdydz
x2 y2 z2
z x y和 z = 1 所围成。 0
1 dudxdydz
2
被积函数在求积区域上的最大 值为2。所以有四维超立方体
y=2*P(:,2);
II=find(y<=2-x.^2&y.^3>=x.^2);
M=length(II);
S=4*M/10000 plot(x(II),y(II),'g.')
1 1
S=
2
2.1136
3 2
S ( 2 x x ) dx
32 15
2/12
例5.14 计算
D
xy 2 dxdy 其中D为y= x – 2与y2 = x 所围
2 2
(x
R 2
) y
2 2
R2 4
V
R
0
Rx x 2 Rx x 2
2 R x y dy dx R . 3 2 3
2 2 2
4
3
前页 后页 返回
实验参考程序 蒙特卡罗方法计算体积 syms x y; f=sqrt(4-x^2-y^2); y1=-sqrt(2*x-x^2);y2=sqrt(2*x-x^2); S1=int(f,y,y1,y2);S2=int(S1,x,0,2) V=double(S2) P=rand(10000,3); X=2*P(:,1);Y=2*P(:,2)-1;Z=2*P(:,3); II=find((X1).^2+Y.^2<=1&Z<=sqrt(4-X.^2Y.^2)); V1=8*length(II)/10000 Error=V-V1
z
y
R
x
前页 后页 返回
积. 在第一卦限内的立体是一个曲顶柱体, 其底为 xy 平面内由 y 0 和 x 2 y 2 Rx 所确定的区域 D (图21-32),而曲顶的方程为 z
R 2 x 2 y 2 . 所以
V 4 R 2 x 2 y 2 d ,
D
其中 D ( x , y ) y 0 , x y Rx .