1、 计 算 方 法实验报告实验序号：实验三 实验名称：高维积分数值计算的MonteCarlo 方法 实 验 人：杨莉专业年级：06 计算机科学与技术教 学 班：C01学 号：0061436实验时间：2009.04.10江西财经大学信息管理学院1实验三 高维积分数值计算的Monte-Carlo 方法一、实验目的学会计算机上高维随机点的选取；学会计算区域测度的随机抽样方法；研究 Monte-Carlo 方法的计算过程和数值性质二、实验内容对于一般的区域 ，计算其测度（只要理解为平面上的面积和空间中的体积）的一般D方法是：先找一个规则的区域 包含 ，且 的测度 是已知的生成 中 个均A)(AmAN匀分
2、布的点 ，如果其中有 个落在区域 中，则区域 的测度近似为)()2(1,Nxx nD)(m函数 在区域 上的积分可以近似为区域 的测度与函数 在区域 中的 个随)(xfD)(xfn机值的平均值的乘积，即DxkDkfndxf)(1)()(xkfNAm)()(三、实验过程假设冰淇淋的下部为一锥体而上面为一半球，考虑冰淇淋锥的体积问题：计算锥面 上方和球面 的内部区域的体积如果使用球面坐标，22yxz1)(22zyx该区域可以表示为如下的积分：40cos202ind用 Monte-Carlo 方法计算该积分；图 1 冰淇淋图2表 1 体积对比 1精确体积 V v1=yl1(1000) v2=yl1(
3、2000) v3=yl1(3000) v4=yl1(4000) v5=yl1(5000)3.1416 3.7417 3.5124 4.2150 4.1665 3.8812表 2 体积对比 2精确体积 V v6=yl1(6000) v7=yl1(7000) v8=yl1(8000) v9=yl1(9000) v10=yl1(10000)3.1416 4.0116 4.0196 3.7674 4.1042 3.98490 2 4 6 8 1033.544.5值值值图 2 体积散点图显然这样的冰淇淋可以装在如下立方体的盒子里,1,zyx而该立方体的体积为生成这个盒子里均匀分布的随机点，落入冰淇淋锥中
4、的点数与总点数之比再乘以得到的就是冰淇淋锥的体积将该结果与中的结果进行比较；表 3 随机点体积对比精确体积 V v=yl2(10) v=yl2(10) v=yl2(10)个数 m 1000 531 396 94体积 V 8 4.2480 3.1680 0.7520比较：由以上表 1、2、3 可知 Monte-Carlo 方法所求体积要优于表 3 求出的体积，且体积比较固定，应为所取为平均值，而表 3 的随机性太大，导致所求体积不太精确。类似的方法可以计算复杂区域的测度（面积或体积） 试求由下列关系所界定区域的测度：（） 0sin1,2,0zyexx图 3.1 区域图 13利用 MATLAB 编
5、程求得区域的侧度为：表 4 区域 1 体积v=yl31(10) v=yl31(10) v=yl31(10) v=yl31(10) v=yl31(10)个数 m 136 150 464 265 258体积 V 1.0880 1.2000 3.7120 2.1200 2.0640（） 2941,3xey图 3.2 区域图 2 利用 MATLAB 编程求得区域的侧度为：表 5 区域 2 体积v=yl32(10) v=yl32(10) v=yl32(10) v=yl32(10)个数 m 4 6 22 12体积 V 0.4000 0.6000 2.200 1.2000（） 1sin0,02yezxz图
6、3.3 区域图 3利用 MATLAB 编程求得区域的侧度为：4表 6 区域 3 体积v=yl33(10) v=yl33(10) v=yl33(10) v=yl33(10)个数 m 844 788 660 549体积 V 0.8440 0.7880 0.6600 0.5490四、实验结论通过以上实验可知：Monte-Carlo 方法中通过样本均值求得的体积与利用概率所求体积相差较大，且概率法的随机性太大，所求的体积不确定性较大，变化范围太广，也许将用概率所求体积在进行平均可能效果会好一点。5附录：1、Monte-Carlo 方法 1-样本均值法 MATLAB 程序：function v=yl1(
7、n)x=rand(1,n);y=rand(1,n);f_g=fcnchk(1+sqrt(1-x.2-y.2)-sqrt(x.2+y.2);k=0;s=0;for i=1:nif x(i).2+y.2=0m=m+1;endendendendmv=m/n.3*8;5、第三题第一问 MATLAB 画图程序：function huatu2x=linspace(0,1,100);y=linspace(1,2,100);x,y=meshgrid(x,y);z=0*x-1;mesh(x,y,z)hold onz=0*x+3;mesh(x,y,z)z=0*x;mesh(x,y,z)z=linspace(-1,
8、3,100);z=exp(x)-y;mesh(x,y,z)hold off76、第三题第二问MATLAB程序：function v=yl32(n)x=1+2*rand(1,n);y=-1+5*rand(1,n);m=0;x y=meshgrid(x,y);s=x.3+y.3;p=exp(x)-2-y;for i=1:nfor j=1:nif s(i,j)=0endend8endendmv=m/n.3;9、第三题第三问MATLAB画图程序：function huatu4x=linspace(0,1,100);y=linspace(0,1,100);x,y=meshgrid(x,y);z=0*x;mesh(x,y,z)hold onz=x.2+sin(y);mesh(x,y,z)z=x+exp(y)-1;mesh(x,y,z)z=0*x+1;mesh(x,y,z)hold off
