清华数学与实验实验五 蒙特卡罗方法 .pptVIP

*/10 ? ? ? 实验五:蒙特卡罗方法实验 面积、体积计算问题 冰淇淋锥的体积计算 思考题与练习题 蒙特卡罗方法——随机投点试验求近似解 引例. 给定曲线y =2 – x2 和曲线y3 = x2，曲线的交点为：P1( – 1，1 )、P2( 1，1 )。曲线围成平面有限区域，用蒙特卡罗方法计算区域面积。 P=rand(10000,2); x=2*P(:,1)-1; y=2*P(:,2); II=find(y=2-x.^2y.^3=x.^2); M=length(II); S=4*M/10000 plot(x(II),y(II),g.) S = 2.1136 例5.14 计算 其中D为y= x – 2与y2 = x 所围 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.^2x=y+2z=x.*(y.^2)); M=length(II); V=192*M/10000 例5.15 用蒙特卡罗方法计算 其中,积分区域是由 和 z = 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(zsqrt(x.^2+y.^2)z=1u=x.^2+y.^2+z.^2); M=length(II); V=8*M/10000 实验:蒙特卡罗方法计算体积 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 [q,error]=MonteC(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 error=q-pi; 实验参考程序 蒙特卡罗方法计算体积 半球体积 圆锥体积 实验任务一：记录L次实验的实验数据及误差 实验任务二：修改实验程序MonteC计算L次实验数据均值及均值误差( mean 计算平均值 ) 序号 1 2 3 4 5 6 7 数据 误差 L 8 16 32 64 128 256 均值 误差 function icecream(m,n) if nargin==0,m=20;n=100;end t=linspace(0,2*pi,n); r=linspace(0,1,m); x=r*cos(t);y=r*sin(t); z1=sqrt(x.^2+y.^2); z2=1+sqrt(1+eps-x.^2-y.^2); X=[x;x];Y=[y;y]; Z=[z1;z2]; mesh(X,Y,Z) view(0,-18) colormap([0 0 1]),axis off 冰淇淋锥体积 冰淇淋锥图形绘制程序