《蒙特卡罗算法举例》由会员分享，可在线阅读，更多相关《蒙特卡罗算法举例（5页珍藏版）》请在金锄头文库上搜索。
1、首先看看上面这个问题。这个问题是我在一个MATLAB交流群里碰到的提问，计算阴影部分面积。什么是蒙特卡罗在这里我就不多做介绍了，感兴趣的朋友可以自己去查阅相关资料，相信可以得到全面的解释，在这里我只介绍如果用蒙特卡罗方法来计算上图中阴影部分的面积，注意这只是蒙特卡罗方法的一个应用而已。记得第一次接触到蒙特卡罗是在一次数学建模培训中。当时我们老师给我们讲了一个故事，故事的全部我已经记不清了，大概内容是：一个古人要求一个图形的面积，他把图形画在一块方形布上，然后找来一袋豆子，然后将所有豆子洒在布上， 落在图形内豆子的重量比上那块布上所有豆子的重量再乘以布的面积就是他所要求的图形的面积。这确实是一个
2、求面积的好方法，这是我听到这个故事后的第一反应。从此我就记住了这个方法，记得很深刻。 所以当群里有人问如何求上面这个图形的面积的时候我马上就回想起用蒙特卡罗方法来计算。求圆周率 , 方法一：随机投点clearl=100;n=10000;x = 2*rand(n,1)-1; y = 2*rand(n,1)-1; p=zeros(l,1);for j=1:lk=0;for i=1:nd = x(i)2+y(i)2;if d=1k=k+1; endendp(j)=k*4/n;endpi=mean(p) 或者m = 10000; pivalue = zeros(m,1); x = 2*rand(m,1
3、)-1; y = 2*rand(m,1)-1; d = x.2+y.2;pivalue= 4*sum(d = 1)/m 方法二：均匀分布点仔细思考后，以我的知识面我能找到两种编程思路来计算这个面积：方法一：将整个坐标轴看成一个边长为12 的正方形，然后均匀的这个正方形分成N（N 的大小取决于划分的步长）个点，然后找出N 个点中有多少个点是属于阴影部分中，假设这个值为k ，则阴影部分的面积为：k/N*122 方法二：将整个坐标轴看成一个边长为12 的正方形，然后在（-6 ，6 ）中随机出N（N 越大越好，至少超过1000 ）个点，然后找出这N 个点中有多少个点在阴影区域内，假设这个值为k ，则阴
4、影部分的面积为：k/N*122。然后重复这个过程100次，求出 100次面积计算结果的均值，这个均值为阴影部分面积。对比分析： 以上两个方法都是利用蒙特卡罗方法计算阴影部分面积，只是在处理的细节有一点区别。前者是把豆子均匀分布在布上；后者则是随机把豆子仍在布上。就计算结果的精度而言，前者取决点的分割是否够密，即N 是否够大；后者不仅仅通过N 来控制精度，因为随机的因素会造成单次计算结果偏高和偏小，所以进行反复多次计算最后以均值来衡量阴影部分面积。附上 MATLAB程序：方法一：clear x=-6:0.01:6; y=x; s=size(x); zs=s(1,2)2; k=0; for i=1
5、:s(1,2) for j=1:s(1,2) a1=(x(i)2)/9+(y(j)2)/36; a2=(x(i)2)/36+y(j)2; a3=(x(i)-2)2+(y(j)+1)2; if a11 if a21 if a39 k=k+1; end end end end end mj=(122)*k/zs; 运行结果：mj = 7.2150 方法二：clear N=10000; n=100; for j=1:n k=0; for i=1:N a=12*rand(1,2)-6; x(i)=a(1,1); y(i)=a(1,2); a1=(x(i)2)/9+(y(i)2)/36; a2=(x(i
6、)2)/36+y(i)2; a3=(x(i)-2)2+(y(i)+1)2; if a11 if a21 if a39 k=k+1; end end end end m(j)=(122)*k/N; end mj=mean(m); 运行结果：mj = 7.2500 改进方法二：m=length(x);n=1000;for i=1:nx=12*rand(10000,1)-6;y=12*rand(10000,1)-6; a1=(x.2)/9+(y.2)/36;a2=(x.2)/36+y.2;a3=(x-2).2+(y+1).2;k=sum(a11s(i)=122*k/m;endmj=mean(s)绘出
7、简单阴影图：clearfigure h1=ezplot(x2/9+y2/36-1);set(h1,color, r);hold onh2=ezplot(x2/36+y2-1);set(h2,color, r);hold onh3=ezplot(x-2)2+(y+1)2-9);set(h3,color, r);hold on%x=-6:0.1:6;%均匀分布的点% y=-6:0.1:6; x=12*rand(1000,1)-6;% 随机投点y=12*rand(1000,1)-6;% n=length(x);for j=1:nfor i=1:na1=(x(j)2)/9+(y(i)2)/36;a2=(x(j)2)/36+y(i)2;a3=(x(j)-2)2+(y(i)+1)2;if a11if a21if a39plot(x(j),y(i),b);hold onendendendendend
