蒙特卡羅(MonteCarlo)方法

——簡介及算法實現(xiàn)(C++版&Java版)物理學(xué)院0510267王洪彬蒙特卡洛法是什么？蒙特卡洛(MonteCarlo)方法，或稱計算機隨機模擬方法，是一種基于“隨機數(shù)”的計算方法。這一方法源于美國在第二次世界大戰(zhàn)中研制原子彈的“曼哈頓計劃”。該計劃的主持人之一、數(shù)學(xué)家馮·諾伊曼用馳名世界的賭城—摩納哥的MonteCarlo—來命名這種方法，為它蒙上了一層神秘色彩。

MonteCarlo方法的基本思想很早以前就被人們所發(fā)現(xiàn)和利用。早在17世紀，人們就知道用事件產(chǎn)生的“頻率”來近似事件的“概率”。19世紀人們用投針試驗的方法來決定圓周率π。本世紀40年代電子計算機的出現(xiàn)，特別是近年來高速電子計算機的出現(xiàn)，使得用數(shù)學(xué)方法在計算機上大量、快速地模擬這樣的試驗成為可能。

考慮平面上的一個邊長為1的正方形及其內(nèi)部的一個形狀不規(guī)則的“圖形”，如何求出這個“圖形”的面積呢？MonteCarlo方法是這樣一種“隨機化”的方法：向該正方形“隨機地”投擲N個點，若有M個點落于“圖形”內(nèi)，則該“圖形”的面積近似為M/N。圓周率的值π=3.

14159265358979323846264338327950288419716939937510

58209749445923078164062862089986280348253421170679

82148086513282306647093844609550582231725359408128

48111745028410270193852110555964462294895493038196

44288109756659334461284756482337867831652712019091

45648566923460348610454326648213393607260249141273

72458700660631558817488152092096282925409171536436

78925903600113305305488204665213841469519415116094

33057270365759591953092186117381932611793105118548

07446237996274956735188575272489122793818301194912

98336733624406566430860213949463952247371907021798

60943702770539217176293176752384674818467669405132

00056812714526356082778577134275778960917363717872

14684409012249534301465495853710507922796892589235

420199561121290219608640344181598136297747713.....

用該方法計算π的基本思路是：

1

根據(jù)圓面積的公式：

s=πR^2

,當R=1時，S=π。

由于圓的方程是：x^2+y^2=1(x^2為x的平方的意思),因此1/4圓面積為x軸、y軸和上述方程所包圍的部分。

如果在1*1的正方形中均勻地落入隨機點，則落入1/4圓中的點的概率就是1/4圓的面積。其4倍，就是圓面積。

由于半徑為1，該面積的值為π的值。

#include<iostream.h>#include<stdlib.h>#include<time.h>constlongN=2000000000;/*定義隨機點數(shù)*/voidmain(){

intn=0;doublex,y;/*坐標*/

srand(time(00));for(inti=1;i<=N;i++){ x=(double)rand()/RAND_MAX;

y=(double)rand()/RAND_MAX;/*在0~1之間產(chǎn)生一個隨機x,y坐標*/

if(x*x+y*y<=1.0)n++;/*統(tǒng)計落入單位圓中的點數(shù)*/}

cout<<"ThePIis"<<4*(double)n/N<<endl;/*計算出π的值*/}計算機模擬實驗數(shù)據(jù)處理1234513.1415333.1415193.1415273.1415363.14152923.1415283.1415283.1415093.1415533.14150633.1415273.1415213.1415373.1415273.1415380.0000120.00000323.14152790.000003314實驗簡評由實驗結(jié)果知此結(jié)果與我們前面給的值相差太大。方法不是很好，但可以通過擴大隨機測試值來獲得更精確的結(jié)果。注：實驗中我應(yīng)用了“iomanip.h”中的setprecision(int);來設(shè)置輸出精度。一道積分題我做過一道證明積分不等式的題：我想中間的積分值可以用蒙特卡洛法求得因為它是一個二重積分，其幾何直觀為一個立體的體積，很巧的是它可以完全包含于一個棱長為1的正方體中，我們在其中產(chǎn)生隨機點，其中落于所求體積的點數(shù)與正方體中產(chǎn)生的點數(shù)之比即為所求的積分值。下面是它的C++代碼：#include<iostream.h>#include<stdlib.h>#include<math.h>#include<time.h>constlongN=2000000000;/*定義隨機點數(shù)*/voidmain(){longn=0; doublex,y,z;/*坐標*/

srand(time(00));for(longi=1;i<=N;i++) {doublem; x=(double)rand()/RAND_MAX; y=(double)rand()/RAND_MAX; z=(double)rand()/RAND_MAX;/*在0~1之間產(chǎn)生一個隨機x,y,z坐標*/

m=0-(x*x+y*y); if(z<=exp(m))n++;/*統(tǒng)計落入所求體積中的點數(shù)*/ }

cout<<"Theintegralis"<<(double)n/N<<endl;/*計算出積分值*/}兩端的界值下面是一次運行結(jié)果：恰好界于0.496466325949717和0.766666666666667之間Theintegralis0.5577358下面是我編寫的Java程序，其實現(xiàn)方法與我編寫的C++程序是相同的：publicclassCalculatePI{publicstaticfinalintN=2000000000;publicstaticvoidmain(String[]args){doublex,y;intn=0;for(inti=0;i<N;i++){x=Math.random();y=Math.random();

if(x*x+y*y<=1)n++;}System.out.println("ThePIis"+4*(double)n/N);}}publicclassCalculateIntegral{publicstaticfinalintN=2000000000;

publicstaticvoidmain(String[]args){doublex,y,z,m;intn=0;for(inti=0;i<N;i++){x=Math.random();y=Math.random();z=Math.random();m=0-x*x-y*y;

if(z<=Math.exp(m))n++;}System.out.println("TheIntegralis"+(double)n/N);}}蒙特卡洛方法是一種比較粗糙的計算方法，它與其說是一種數(shù)學(xué)的方法，不如看成一種實驗方法；有的時候人們不得不依賴于此種方法。一來人們的才智不夠，二來它的確在某些方面給我們知識；但不可否認對于可控參數(shù)比較少的情況，蒙特卡洛法不失為一個對研究對象獲得認識的一個有效方法；從計算來說，用蒙特卡洛法計算值確實粗糙，我們一般用更好的方法來計算。比如說級數(shù)的方法，下面是我用級數(shù)的方法編的C++和Java程序：//CalculatePI_Ad.cpp#include<iostream.h>#include<iomanip.h>voidmain(){ doublesum=0.0,f=1; for(inti=1;;i+=2) { sum+=f/i; if(1.0/i<=1e-14) break; f=-f; } cout<<"ThePIis"<<setprecision(15)<<4*sum<<endl;}ThePIis3.14159265079415publicclassCalculatePI_Ad{publicstaticvoidmain(String[]args){doublesum=0.0,f=1;for(inti=1;;i+=2){sum+=f/i;if(1.0/i<=1e-15)