蒙特卡羅算法舉例

1、    首先看看上面這個(gè)問題。這個(gè)問題是我在一個(gè)MATLAB交流群里碰到的提問，計(jì)算陰影部分面積。    什么是蒙特卡羅在這里我就不多做介紹了，感興趣的朋友可以自己去查閱相關(guān)資料，相信可以得到全面的解釋，在這里我只介紹如果用蒙特卡羅方法來計(jì)算上圖中陰影部分的面積，注意這只是蒙特卡羅方法的一個(gè)應(yīng)用而已。    記得第一次接觸到蒙特卡羅是在一次數(shù)學(xué)建模培訓(xùn)中。當(dāng)時(shí)我們老師給我們講了一個(gè)故事，故事的全部我已經(jīng)記不清了，大概內(nèi)容是：    一個(gè)古人要求一個(gè)圖形的面積，他把圖形畫在一塊

2、方形布上，然后找來一袋豆子，然后將所有豆子灑在布上，落在圖形內(nèi)豆子的重量比上那塊布上所有豆子的重量再乘以布的面積就是他所要求的圖形的面積。    這確實(shí)是一個(gè)求面積的好方法，這是我聽到這個(gè)故事后的第一反應(yīng)。從此我就記住了這個(gè)方法，記得很深刻。所以當(dāng)群里有人問如何求上面這個(gè)圖形的面積的時(shí)候我馬上就回想起用蒙特卡羅方法來計(jì)算。    仔細(xì)思考后，以我的知識(shí)面我能找到兩種編程思路來計(jì)算這個(gè)面積：方法一：將整個(gè)坐標(biāo)軸看成一個(gè)邊長(zhǎng)為12的正方形，然后均勻的這個(gè)正方形分成N（N的大小取決于劃分的步長(zhǎng)）個(gè)點(diǎn)，然后找出N個(gè)點(diǎn)中有多少個(gè)點(diǎn)是屬于陰影部分

3、中，假設(shè)這個(gè)值為k，則陰影部分的面積為：k/N*122方法二：將整個(gè)坐標(biāo)軸看成一個(gè)邊長(zhǎng)為12的正方形，然后在（-6，6）中隨機(jī)出N（N越大越好，至少超過1000）個(gè)點(diǎn)，然后找出這N個(gè)點(diǎn)中有多少個(gè)點(diǎn)在陰影區(qū)域內(nèi)，假設(shè)這個(gè)值為k，則陰影部分的面積為：k/N*122。然后重復(fù)這個(gè)過程100次，求出100次面積計(jì)算結(jié)果的均值，這個(gè)均值為陰影部分面積。對(duì)比分析：以上兩個(gè)方法都是利用蒙特卡羅方法計(jì)算陰影部分面積，只是在處理的細(xì)節(jié)有一點(diǎn)區(qū)別。前者是把豆子均勻分布在布上；后者則是隨機(jī)把豆子仍在布上。就計(jì)算結(jié)果的精度而言，前者取決點(diǎn)的分割是否夠密，即N是否夠大；后者不僅僅通過N來控制精度，因?yàn)殡S機(jī)的因素會(huì)造成單

4、次計(jì)算結(jié)果偏高和偏小，所以進(jìn)行反復(fù)多次計(jì)算最后以均值來衡量陰影部分面積。附上MATLAB程序：方法一：clearx=-6:0.01:6;y=x;s=size(x);zs=s(1,2)2;k=0;for i=1: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;    &#

5、160;   a3=(x(i)-2)2+(y(j)+1)2;        if a1<1            if a2<1                if a3<9     

6、;               k=k+1;                end            end        e

7、nd    endendmj=(122)*k/zs;運(yùn)行結(jié)果：mj =    7.2150方法二：clearN=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;    a

8、2=(x(i)2)/36+y(i)2;    a3=(x(i)-2)2+(y(i)+1)2;    if a1<1       if a2<1          if a3<9             k=k+1;           end