計算材料學(xué)概述-之-蒙特卡洛方法.詳解

計算材料學(xué)概述第三章蒙特卡羅方法（MonteCarlo）主要內(nèi)容MonteCarlo模擬發(fā)展簡介MonteCarlo模擬基本原理MonteCarlo模擬典型算法MonteCarlo模擬典型應(yīng)用蒙特卡洛法是什么？蒙特卡洛(MonteCarlo)方法，是在簡潔的理論準則基礎(chǔ)上，接受反復(fù)隨即抽樣的方法，解決困難系統(tǒng)的問題。其實質(zhì)是一種概率和統(tǒng)計的問題。蒙特·卡羅方法（MonteCarlomethod），也稱統(tǒng)計模擬方法，是二十世紀四十年頭中期由于科學(xué)技術(shù)的發(fā)展和電子計算機的獨創(chuàng)，而被提出的一種以概率統(tǒng)計理論為指導(dǎo)的一類特別重要的數(shù)值計算方法。是指運用隨機數(shù)（或更常見的偽隨機數(shù)）來解決很多計算問題的方法。MC的基本思想MC基本思想很早以前就被人們所發(fā)覺和利用。17世紀，人們就知道用事務(wù)發(fā)生的“頻率”來確定事務(wù)的“概率”。但要真正實現(xiàn)隨機抽樣是很困難的，甚至幾乎是不行能的。高速計算機的出現(xiàn)，使得用數(shù)學(xué)方法在計算機上大量、快速地模擬這樣的試驗成為可能。確定性系統(tǒng)隨機性系統(tǒng)模擬自然界Monte-Carlo模擬，即隨機模擬（重復(fù)“試驗”）重復(fù)試驗計算機模擬MonteCarlo方法：亦稱統(tǒng)計模擬方法，statisticalsimulationmethod

利用隨機數(shù)進行數(shù)值模擬的方法MonteCarlo名字的由來：是由Metropolis在二次世界大戰(zhàn)期間提出的：Manhattan支配，探討與原子彈有關(guān)的中子輸運過程；MonteCarlo是摩納哥（monaco)的首都，該城以賭博著名NicholasMetropolis(1915-1999)Monte-Carlo,MonacoMonteCarlo方法簡史簡潔地介紹一下MonteCarlo方法的發(fā)展歷史1、Buffon投針試驗：18世紀，法國數(shù)學(xué)家ComtedeBuffon利用投針試驗估計的值dL1777年法國科學(xué)家布豐提出的一種計算圓周率的方法——隨機投針法，即著名的布豐投針問題。這一方法的步驟是：

1)取一張白紙，在上面畫上很多條間距為d的平行線。

2)取一根長度為l（l<d）的針，隨機地向畫有平行直線的紙上擲n次，視察針與直線相交的次數(shù)，記為m

3）計算針與直線相交的概率．

布豐本人證明白，這個概率是:

p=2l/(πd)，π為圓周率：

利用這個公式可以用概率的方法得到圓周率的近似值。下面是一些資料試驗者年頭投擲次數(shù)相交次數(shù)圓周率估計值沃爾夫1850500025313.1596史密斯1855320412193.1554德摩137福克斯188410304893.1595拉澤里尼1901340818083.1415929賴納192525208593.1795布豐投針試驗是第一個用幾何形式表達概率問題的例子，他首次運用隨機試驗處理確定性數(shù)學(xué)問題，為概率論和蒙特卡羅方法的發(fā)展起到確定的推動作用。

MonteCarlo方法之隨機數(shù)的產(chǎn)生很多計算機系統(tǒng)都有隨機數(shù)生成函數(shù)F90:callrandom_seedcallrandom_number(a)2、ISEED=RTC()X=RAN(ISEED)Y=RAN(ISEED）Matlab:x=rand(N)產(chǎn)生元素在(0,1)間隨機分布的N*N矩陣

s=rand(‘state’,0)重設(shè)該生成函數(shù)到初始狀態(tài)留意：上述隨機數(shù)序列均具周期性，如上頁random子程序的周期約230。實例一、計算π值計算過程：1、構(gòu)造或描述問題的概率過程2、從概率密度函數(shù)動身進行隨機抽樣，實現(xiàn)從已知概率分布的抽樣，得到特征量的一些模擬結(jié)果——計算均值MonteCarlo方法之典型算法與應(yīng)用考慮平面上的一個邊長為1的正方形及其內(nèi)部的一個形態(tài)不規(guī)則的“圖形”，如何求出這個“圖形”的面積呢？MonteCarlo方法是這樣一種“隨機化”的方法：向該正方形“隨機地”投擲N個點，若有M個點落于“圖形”內(nèi)，則該“圖形”的面積近似為M/N。用該方法計算π的基本思路是：

1

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

s=πR^2

,當(dāng)R=1時，S=π。

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

3、假如在1*1的正方形中勻整地落入隨機點，則落入1/4圓中的點的概率就是1/4圓的面積。其4倍，就是圓面積。由于半徑為1，該面積的值為π的值。

REALR,R1,R2,PI ISEED=RTC() N0=0 N=300000 DOI=1,N R1=RAN(ISEED) R2=RAN(ISEED) R=SQRT(R1*R1+R2*R2) IF(R<1.0)N0=N0+1 ENDDO PI=4.0*N0/N WRITE(*,*)PI END面積的計算f(x)x辛普遜方法I=ΣSn蒙特-卡洛方法f(x)x在長方形中均勻投N0組(x,y)如y<f(x),則N=N+1I=(N/N0)×S0SS011MC的優(yōu)點MC與傳統(tǒng)數(shù)學(xué)方法相比，具有直觀性強，簡便易行的優(yōu)點，該方法能處理一些其他方法無法解決的負責(zé)問題，并且簡潔在計算機上實現(xiàn)，在很大程度上可以代替很多大型、難以實現(xiàn)的困難試驗和社會行為。無污染、無危急、能擺脫試驗誤差。MonteCarlo方法利用隨機抽樣的方法來求解物理問題;數(shù)值解法:從一個物理系統(tǒng)的數(shù)學(xué)模型動身,通過求解一系列的微分方程來的導(dǎo)出系統(tǒng)的未知狀態(tài);MonteCarlo方法并非只能用來解決包含隨機過程的問題:例如:用MonteCarlo方法計算定積分.

對這樣的問題可將其轉(zhuǎn)換成相關(guān)的隨機過程,然后用MonteCarlo方法進行求解留意以下兩點：MC的應(yīng)用自然現(xiàn)象的模擬：宇宙射線在地球大氣中的傳輸過程；高能物理實驗中的核相互作用過程；數(shù)值分析：數(shù)學(xué)問題，求積分，求逆矩陣，解線性代數(shù)等經(jīng)濟學(xué)模擬：庫存問題，隨機服務(wù)系統(tǒng)中排隊問題人口問題：人口的誕生，傳染病的擴散；乃至動物的生態(tài)競爭金屬學(xué)：擴散、組織長大、相變過程蒙特-卡洛模擬的意義能探討不同邊界、不同材料的影響理論不行能、試驗耗費太大用于試驗設(shè)計無污染反應(yīng)堆防護核彈爆炸能擺脫試驗誤差作理論和試驗的橋梁MonteCarlo模擬的步驟：依據(jù)欲探討的物理系統(tǒng)的性質(zhì)，建立能夠描述該系統(tǒng)特性的理論模型，導(dǎo)出該模型的某些特征量的概率密度函數(shù)；（即構(gòu)造或描述問題的概率過程）從概率密度函數(shù)動身進行隨機抽樣，實現(xiàn)從已知概率分布的抽樣，得到特征量的一些模擬結(jié)果；有了明確的概率過程后，為了實現(xiàn)過程的數(shù)值模擬，必需實現(xiàn)從已知概率分布的隨機數(shù)的抽樣，進行大量的隨機模擬試驗，從中獲得隨機變量的大量試驗值。產(chǎn)生已知概率分布的隨機變量，是實現(xiàn)MC方法的關(guān)鍵步驟，其中最基本的是（0，1）勻整分布。對模擬結(jié)果進行分析總結(jié)，預(yù)言物理系統(tǒng)的某些特性。4.模擬結(jié)果的檢驗MonteCarlo算法的主要組成部分概率密度函數(shù)(pdf)—必需給出描述一個物理系統(tǒng)的一組概率密度函數(shù);隨機數(shù)產(chǎn)生器—能夠產(chǎn)生在區(qū)間[0,1]上（勻整）分布的隨機數(shù)抽樣規(guī)則—如何從在區(qū)間[0,1]上勻整分布的隨機數(shù)動身,隨機抽取聽從給定的pdf的隨機變量;模擬結(jié)果記錄—記錄一些感愛好的量的模擬結(jié)果誤差估計—必需確定統(tǒng)計誤差（或方差）隨模擬次數(shù)以及其它一些量的變更；削減方差的技術(shù)—利用該技術(shù)可削減模擬過程中計算的次數(shù)；并行和矢量化—可以在先進的并行計算機上運行的有效算法實例二定積分計算事實上，不少的統(tǒng)計問題，如計算概率、各階距等，最終都歸結(jié)為定積分的近似計算問題。下面考慮一個簡潔的定積分!計算x**2在(0,1)上積分計算過程：1、構(gòu)造或描述問題的概率過程:產(chǎn)生聽從分布f（x）的隨機變量Xi()（i=1,2,···,N）2、從概率密度函數(shù)動身進行隨機抽樣，實現(xiàn)從已知概率分布的抽樣，得到特征量的一些模擬結(jié)果——計算均值()

REALY Y=0 N=300000 ISEED=RTC() DOI=1,N X=RAN(ISEED) Y=Y+X**2/NENDDO WRITE(*,*)Y END

lim∑x2dx(dx0)MonteCarlo方法另一個重要問題：隨機數(shù)隨機數(shù)：由單位矩陣分布中所產(chǎn)生的簡潔子樣稱為隨機數(shù)序列，其中的每一個個體稱為隨機數(shù)。但真正的隨機數(shù)的不適合電子計算機上運用，因為它須要很大的存儲量。利用某些物理現(xiàn)象可以在電子計算機上產(chǎn)生隨機數(shù)，且其產(chǎn)生的序列無法重復(fù)實現(xiàn)，使程序無法復(fù)算，結(jié)果無法驗證，同時須要增加隨機數(shù)發(fā)生器和電路聯(lián)系等附加設(shè)備。偽隨機數(shù)：是有數(shù)學(xué)遞推公式所產(chǎn)生的隨機數(shù)。（近似的具備隨機數(shù)的性質(zhì)。）An+1=T(A);An+1=An+k+1偽隨機的優(yōu)點和缺點：推斷偽隨機數(shù)好壞的方法：1、它能夠有較好的勻整性和獨立性；2、它的費用大小，即指所消耗計算機的時間；3、容量要求盡可能大。隨機數(shù)產(chǎn)生的方法產(chǎn)生勻整分布隨機數(shù)的幾種方法；（1）物理方法；（2）數(shù)學(xué)方法。偽隨機數(shù)產(chǎn)生方法：加同余法乘同余法乘加同余法取中方法逆變換法合成法篩選法。。。。。。關(guān)于隨機數(shù)的幾點留意注1由于勻整分布的隨機數(shù)的產(chǎn)生總是接受某個確定的模型進行的，從理論上講，總會有周期現(xiàn)象出現(xiàn)的。初值確定后，全部隨機數(shù)也隨之確定，并不滿足真正隨機數(shù)的要求。因此通常把由數(shù)學(xué)方法產(chǎn)生的隨機數(shù)成為偽隨機數(shù)。但其周期又相當(dāng)長，在實際應(yīng)用中幾乎不行能出現(xiàn)。因此，這種由計算機產(chǎn)生的偽隨機數(shù)可以當(dāng)作真正的隨機數(shù)來處理。注2應(yīng)對所產(chǎn)生的偽隨機數(shù)作各種統(tǒng)計檢驗，如獨立性檢驗，分布檢驗，功率譜檢驗等等。FORTRAN語言產(chǎn)生隨機數(shù)的實例random_number(x)

產(chǎn)生一個0到1之間的隨機數(shù)（x可以是向量），但是每次總是那幾個數(shù)。

用了random_seed

()后，系統(tǒng)依據(jù)日期和時間隨機地供應(yīng)種子，使得隨機數(shù)更隨機了。

program

random

real

::

x

call

random_seed

()

!

系統(tǒng)依據(jù)日期和時間隨機地供應(yīng)種子

call

random_number

(x)

!

每次的隨機數(shù)就都不一樣了

write(*,*)

x

stop

end

program

randomMonteCarlo方法之隨機數(shù)的產(chǎn)生很多計算機系統(tǒng)都有隨機數(shù)生成函數(shù)F90:callrandom_seedcallrandom_number(a)2、ISEED=RTC()X=RAN(ISEED)Y=RAN(ISEED）Matlab:x=rand(N)產(chǎn)生元素在(0,1)間隨機分布的N*N矩陣

s=rand(‘state’,0)重設(shè)該生成函數(shù)到初始狀態(tài)留意：上述隨機數(shù)序列均具周期性，如上頁random子程序的周期約230。FinitedifferenceapproximationofdifferentialequationsAdifferentialequationcanbeapproximatedbyafinitedifferencescheme.ForexampleForwardtimeBackwardspaceCentralspaceForwardtaime-CentralspaceFICK其次定律Fick其次定律——穩(wěn)態(tài)擴散解

REALC0(0:1000+1),C(0:1000+1)!C(DISTANCE)

C0=0.1 C0(0)=0.8!BOUNDARYCONDITION C0(1001)=0.1!BOUNDARYCONDITION C=C0OPEN(1,FILE="F:\DIF.dat") DOJT=1,50000！Time DOIX=1,1000!distance C(IX)=C0(IX)+0.45*(C0(IX+1)+C0(IX-1)-2.0*C0(IX))! C(0)=0.8! C(1001)=0.1 IF(JT==50000)WRITE(1,*)IX,C(IX) ENDDO C0=CENDDO END應(yīng)用之二生日問題MC模擬假設(shè)有n個人在一起，各自的生日為365天之一，依據(jù)概率理論，與很多人的直覺相反，只需23個人便有大于50％的幾率人群中至少有2個人生日相同。n理論幾率模擬幾率0.1170.1100.4110.4120.5270.5200.7060.6920.9410.936500.9860.987INTEGERM(1:10000),

NUMBER1(0:364),

NUMBER2REALX,Y

ISEED=RTC() DOJ=1,10000 NUMBER1=0 X=RAN(ISEED) NUMBER1(0)=INT(365*X+1)

JJJ=1 DOI=1,365 Y=RAN(ISEED) NUMBER2=INT(365*Y+1) ETR=COUNT(NUMBER1.EQ.NUMBER2) IF(ETR==1)THEN EXIT ELSE JJJ=JJJ+1 M(J)=JJJ NUMBER1(I)=NUMBER2ENDIFENDDO ENDDO

DOI=1,10000IF(M(I).LE.23)SUM=SUM+1 ENDDOPRINT*,SUM/10000ENDMC在材料學(xué)領(lǐng)域的應(yīng)用

——隨機行走背景如，布朗運動－－－最簡潔、無限制隨機行走(Unrestrictedrandonwalk,RW)startend平均平方端-端位移:,自然科學(xué)和社會生活中很多現(xiàn)象都與隨機運動有關(guān)可以模擬的內(nèi)容？擴散；分子運動；。。。。。。如圖所示，第i個分子在經(jīng)過N步隨機行走后距原點距離為R，對n個分子每步的位移平方求和后取平均值就得到了全部分子距原點的方均距離<R2>：

！MonteCarloSimulationofOneDimensionalDiffusion

INTEGERX,XX(1:1000,1:1000) REALXXM(1:1000)! X:INSTANTANEOUSPOSITIONOFATOM! XX(J,I)：X*X,J:第幾天試驗,I:第幾步跳動! XXM(I):THEMEANOFXX WRITE(*,*)"試驗天數(shù)JMAX，試驗次數(shù)IMAX" READ(*,*)JMAX,IMAX ISEED=RTC() DOJ=1,JMAX!第幾天試驗 X=0!!! DOI=1,IMAX!第幾步跳動 RN=RAN(ISEED) IF(RN<0.5)THEN X=X+1 ELSE X=X-1 ENDIF XX(J,I)=X*X ENDDO ENDDO OPEN(1,FILE=“f:\DIF1.DAT") DOI=1,IMAX XXM=0.0 XXM(I)=1.0*SUM(XX(1:JMAX,I))/JMAX!! WRITE(1,*)I,XXM(I) ENDDO CLOSE(1) END! MonteCarloSimulationofTwoDimensionalDiffusion

INTEGERX,Y,XY(1:1000,1:1000) REALXYM(1:1000)! X:INSTANTANEOUSPOSITIONOFATOM! XY(J,I)：X*Y,J:第幾天試驗,I:第幾步跳動! XYM(I):THEMEANOFXY WRITE(*,*)"試驗天數(shù)JMAX，試驗次數(shù)IMAX" READ(*,*)JMAX,IMAX ISEED=RTC() DOJ=1,JMAX!第幾天試驗 X=0!!! Y=0!!! DOI=1,IMAX!第幾步跳動 RN=RAN(ISEED) IF(RN.LT.0.25)THEN x=x y=y-1 ENDIF IF(RN.LT.0.5.AND.RN.GE.0.25)THEN x=x y=y+1 ENDIF IF(RN.LT.0.75.AND.RN.GE.0.5)THEN x=x-1 y=y ENDIF IF(RN.GE.0.75)THEN x=x+1 y=y ENDIF XY(J,I)=X*X+Y*Y ENDDO ENDDO OPEN(1,FILE=“f:\DIF2.DAT") DOI=1,IMAX XYM=0.0 XYM(I)=1.0*SUM(XY(1:JMAX,I))/JMAX!! WRITE(1,*)I,XYM(I) ENDDO CLOSE(1) END! MonteCarloSimulationofTwoDimensionalDiffusion

INTEGERX,XY(1:1000,1:1000),y,XN(1:4),YN(1:4),RN REALXYM(1:1000)! X:INSTANTANEOUSPOSITIONOFATOM! XY(J,I)：X*Y,J:第幾天試驗,I:第幾步跳動! XYM(I):THEMEANOFXY WRITE(*,*)"試驗天數(shù)JMAX，試驗次數(shù)IMAX" READ(*,*)JMAX,IMAX XN=(/0,0,-1,1/) YN=(/-1,1,0,0/) ISEED=RTC() DOJ=1,JMAX!第幾天試驗 X=0!!! Y=0!!! DOI=1,IMAX!第幾步跳動 RN=4*RAN(ISEED)+1 X=X+XN(RN) Y=Y+YN(RN) XY(J,I)=X*X+Y*Y ENDDO ENDDO OPEN(1,FILE="C:\DIF2.DAT") DOI=1,IMAX XYM=0.0 XYM(I)=1.0*SUM(XY(1:JMAX,I))/JMAX!! WRITE(1,*)I,XYM(I) ENDDO CLOSE(1)！做三維空間隨機行走？ END實際環(huán)境是困難的:可變，缺陷，風(fēng)，季節(jié)，電場等。

空間中有隨機性缺陷不退行走自回避行走例如模擬高分子的位形用隨機行走方法模擬高分子位形是用隨機行走的軌跡代表高分子的位形，行走過的位置代表的是構(gòu)成分子的原子或官能團，因此，無限制隨機行走忽視了體斥效應(yīng)。不退行走就是禁止在每一步行走后馬上倒退，可以解決剛走的一步與上一步重疊的問題。但不退行走沒有完全解決高分子的體斥效應(yīng)。自回避行走就是全部已走過的位置不能再走，這樣就完全解決了體斥效應(yīng)問題。

氣體分子的隨機行走

假設(shè)在一個空曠封閉的房間中心滴上一滴香水，揮發(fā)的香水分子隨機的與空氣中的粒子發(fā)生碰撞，最終會擴散到整個房間。這里我們將通過計算機模擬400個分子在二維平面內(nèi)的隨機運動探討熵變現(xiàn)在來探討一下熵在我們這個氣體擴散模型中是什么意思。在剛起先的時候，我們的系統(tǒng)處在一個特別有序的狀態(tài)：全部的分子都處于原點。這時系統(tǒng)的熵為零，系統(tǒng)不存在任何的無序度。隨著時間的推移，分子起先不斷的向外擴散，這時系統(tǒng)出現(xiàn)了無序狀態(tài)，熵起先漸漸增加。系統(tǒng)的熵和我們預(yù)料的一樣在隨時間增大，但當(dāng)全部分子布滿整個邊界以內(nèi)區(qū)域時，系統(tǒng)的熵起先趨于穩(wěn)定。從這一結(jié)果中我們可以了解：當(dāng)全部的分子隨機的布滿整個區(qū)域時，雖然當(dāng)我們跟蹤某一個確定分子時，它還是在區(qū)域內(nèi)到處亂竄，但每個小區(qū)域內(nèi)分子的密度卻不會再變更了。所以，一旦氣體分子擴散到整個區(qū)域以后，不管我們再等上多少時間，系統(tǒng)的熵都不會再有太大的起伏。換句話說，讓系統(tǒng)自動回到起先的狀態(tài)，即全部分子都在原點的狀態(tài)，已經(jīng)不行能了。畫一個圓晶粒

USEMSFLIB! INTEGERXR,YR!在的區(qū)域中畫一個圓

PARAMETERXR=400,YR=400 INTEGERR,S(1:XR,1:YR) X0=XR/2!圓心位置X0，YO Y0=YR/2 R=MIN(X0-10,Y0-10)!圓半徑

S=0!像素的初始狀態(tài)（顏色）

DOI=1,XR DOJ=1,YR IF((I-X0)**2+(J-Y0)**2<=R**2)S(I,J)=10 IER=SETCOLOR(S(I,J)+3) IER=SETPIXEL(I,J) ENDDO ENDDO END畫晶界! 畫一個圓 USEMSFLIB INTEGERXR，YR!在的區(qū)域中畫一個圓 PARAMETERXR=400,YR=400 INTEGERR,S(0:XR+1,0:YR+1),XN(1:4),YN(1:4),SNS XN=(/0,0,-1,1/) YN=(/-1,1,0,0/) X0=XR/2!圓心位置X0，YO Y0=YR/2 R=MIN(X0-10,Y0-10)!圓半徑 S=0!像素的初始狀態(tài)（顏色） DOI=1,XR DOJ=1,YR IF((I-X0)**2+(J-Y0)**2<=R**2)S(I,J)=10 IER=SETCOLOR(S(I,J)) IER=SETPIXEL(I,J) ENDDO ENDDO DOI=1,XR!畫晶界 DOJ=1,YR NDS=0 DOK=1,4 IF(S(I,J).NE.S(I+XN(K),J+YN(K)))NDS=NDS+1 ENDDO IF(NDS>0)THEN IER=SETCOLOR(9) ELSE IER=SETCOLOR(8) ENDIF IER=SETPIXEL(I,J) ENDDO ENDDO END如何畫有確定寬度的晶界？例題例求前N個自然數(shù)倒數(shù)和。即求程序須要幾個變量：累加的項數(shù)I(整型)，存放通項值的T(實型)，存放總和的S(實型)，待輸入的總項數(shù)N(整型)。PROGRAMMAINIMPLICITNONEINTEGERN,IREALT,SREAD*,NS=0.0;I=1DOT=1.0/IS=S+TI=I+1IF(I>N)EXITENDDOPRINT*,"SUM=",SENDPROGRAMMAIN程序的運行結(jié)果如下15

SUM=3.182290用DO結(jié)構(gòu)處理循環(huán)問題，最關(guān)鍵的是要處理好初值、循環(huán)終止條件和循環(huán)體語句的關(guān)系。FinitedifferenceapproximationofdifferentialequationsAdifferentialequationcanbeapproximatedbyafinitedifferencescheme.ForexampleForwardtimeBackwardspaceCentralspaceForwardtime-CentralspaceFICK其次定律Fick其次定律——穩(wěn)態(tài)擴散解

REALC0(0:1000+1),C(0:1000+1)!C(DISTANCE)

C0=0.1 C0(0)=0.8!BOUNDARYCONDITION C0(1001)=0.1!BOUNDARYCONDITION C=C0OPEN(1,FILE="F:\DIF.dat")

DOJT=1,50000！Time DOIX=1,1000!distance C(IX)=C0(IX)+0.45*(C0(IX+1)+C0(IX-1)-2.0*C0(IX)) C(0)=0.8 C(1001)=0.1 IF(JT==50000)WRITE(1,*)IX,C(IX) ENDDO C0=CENDDO PRINT*,"PROGRAMMEISFINISHED"CLOSE(1)END蒙特卡羅(MonteCarlo)方法1、簡介2、相關(guān)幾個例子MonteCarlo方法：亦稱統(tǒng)計模擬方法，statisticalsimulationmethod

利用隨機數(shù)進行數(shù)值模擬的方法MonteCarlo名字的由來：是由Metropolis在二次世界大戰(zhàn)期間提出的：Manhattan支配，探討與原子彈有關(guān)的中子輸運過程；MonteCarlo是摩納哥（monaco)的首都，該城以賭博著名NicholasMetropolis(1915-1999)Monte-Carlo,MonacoMonteCarlo方法的基本思想很早以前就被人們所發(fā)覺和利用。早在17世紀，人們就知道用事務(wù)產(chǎn)生的“頻率”來近似事務(wù)的“概率”。18世紀人們用投針試驗的方法來確定圓周率π。本世紀40年頭電子計算機的出現(xiàn)，特殊是近年來高速電子計算機的出現(xiàn)，使得用數(shù)學(xué)方法在計算機上大量、快速地模擬這樣的試驗成為可能。1930年，EnricoFermi利用MonteCarlo方法探討中子的擴散，并設(shè)計了一個MonteCarlo機械裝置，F(xiàn)ermiac,用于計算核反應(yīng)堆的臨界狀態(tài)VonNeumann是MonteCarlo方法的正式奠基者,他與StanislawUlam合作建立了概率密度函數(shù)、反累積分布函數(shù)的數(shù)學(xué)基礎(chǔ)，以及偽隨機數(shù)產(chǎn)生器。在這些工作中，StanislawUlam意識到了數(shù)字計算機的重要性合作起源于Manhattan工程：利用ENIAC(ElectronicNumericalIntegratorandComputer)計算產(chǎn)額MonteCarlo模擬的應(yīng)用：自然現(xiàn)象的模擬：宇宙射線在地球大氣中的傳輸過程；高能物理實驗中的核相互作用過程；試驗探測器的模擬數(shù)值分析：利用MonteCarlo方法求積分材料學(xué)中的應(yīng)用：擴散、晶粒的長大、再結(jié)晶等MonteCarlo模擬在物理探討中的作用MonteCarlo模擬的步驟：依據(jù)欲探討的物理系統(tǒng)的性質(zhì)，建立能夠描述該系統(tǒng)特性的理論模型，導(dǎo)出該模型的某些特征量的概率密度函數(shù)；從概率密度函數(shù)動身進行隨機抽樣，得到特征量的一些模擬結(jié)果；對模擬結(jié)果進行分析總結(jié)，預(yù)言物理系統(tǒng)的某些特性。留意以下兩點：MonteCarlo方法與數(shù)值解法的不同:MonteCarlo方法利用隨機抽樣的方法來求解物理問題;數(shù)值解法:從一個物理系統(tǒng)的數(shù)學(xué)模型動身,通過求解一系列的微分方程來的導(dǎo)出系統(tǒng)的未知狀態(tài);MonteCarlo方法并非只能用來解決包含隨機的過程的問題:很多利用MonteCarlo方法進行求解的問題中并不包含隨機過程例如:用MonteCarlo方法計算定積分.對這樣的問題可將其轉(zhuǎn)換成相關(guān)的隨機過程,然后用MonteCarlo方法進行求解MonteCarlo算法的主要組成部分概率密度函數(shù)—必需給出描述一個物理系統(tǒng)的一組概率密度函數(shù);隨機數(shù)產(chǎn)生器—能夠產(chǎn)生在區(qū)間[0,1]上勻整分布的隨機數(shù)抽樣規(guī)則—如何從在區(qū)間[0,1]上勻整分布的隨機數(shù)動身,隨機抽取聽從給定的pdf的隨機變量;模擬結(jié)果記錄—記錄一些感愛好的量的模擬結(jié)果誤差估計—必需確定統(tǒng)計誤差（或方差）隨模擬次數(shù)以及其它一些量的變更；削減方差的技術(shù)—利用該技術(shù)可削減模擬過程中計算的次數(shù)；并行和矢量化—可以在先進的并行計算機上運行的有效算法確定性系統(tǒng)隨機性系統(tǒng)模擬自然界Monte-Carlo模擬，即隨機模擬（重復(fù)“試驗”）重復(fù)試驗計算機模擬蒙特卡羅方法的基本原理及思想當(dāng)所要求解的問題是某種事務(wù)出現(xiàn)的概率，或者是某個隨機變量的期望值時，它們可以通過某種“試驗”的方法，得到這種事務(wù)出現(xiàn)的頻率，或者這個隨機變數(shù)的平均值，并用它們作為問題的解。把蒙特卡羅解題歸結(jié)為三個主要步驟： 構(gòu)造或描述概率過程； 實現(xiàn)從已知概率分布抽樣； 建立各種估計量。對于本身就具有隨機性質(zhì)的問題，主要是正確描述和模擬這個概率過程；對于不是隨機性質(zhì)的確定性問題，比如計算定積分，就必需事先構(gòu)造一個人為的概率過程。最簡潔、最基本、最重要的一個概率分布是(0,1)上的勻整分布（或稱矩形分布）。實現(xiàn)從已知概率分布抽樣(模擬試驗)：構(gòu)造了概率模型以后，產(chǎn)生已知概率分布的隨機變量（或隨機向量），就成為實現(xiàn)蒙特卡羅方法模擬試驗的基本手段，這也是蒙特卡羅方法被稱為隨機抽樣的緣由。建立各種估計量：構(gòu)造了概率模型并能從中抽樣后，即實現(xiàn)模擬試驗后，就要確定一個隨機變量，作為所要求的問題的解，我們稱它為無偏估計。收斂性:

由大數(shù)定律,Monte-Carlo模擬的收斂是以概率而言的．誤差:

用頻率估計概率時誤差的估計，可由中心極限定理，給定置信水平的條件下，有：

模擬次數(shù):由誤差公式得隨機數(shù)的產(chǎn)生原理-計算機抽樣試驗隨機數(shù)的產(chǎn)生是實現(xiàn)MC計算的先決條件。而大多數(shù)概率分布的隨機數(shù)的產(chǎn)生都是基于勻整分布U(0,1]的隨機數(shù)。首先，介紹聽從勻整分布U(0,1]的隨機數(shù)的產(chǎn)生方法。其次，介紹聽從其他各種分布的隨機數(shù)的產(chǎn)生方法。以及聽從正態(tài)分布的隨機數(shù)的產(chǎn)生方法。最終，關(guān)于隨機數(shù)的幾點注。(0,1]勻整分布隨機數(shù)的產(chǎn)生最常用的是在(0,1]區(qū)間內(nèi)勻整分布的隨機數(shù)，其他分布的隨機數(shù)可利用勻整分布隨機數(shù)來產(chǎn)生。勻整分布隨機數(shù)的產(chǎn)生：乘同余法產(chǎn)生區(qū)間(0,1]內(nèi)勻整分布的隨機數(shù)的遞推公式是：

其中，是乘因子，M是模數(shù)。第一式稱做以M為模數(shù)（modulus）的同余式，即以M除后得到的余數(shù)記為。當(dāng)給定了一個初值（稱為種子），計算出的即在(0,1]上勻整分布的隨機數(shù)。

例1：R=RAND()ISEED=RTC()R=RAN(ISEED)其他各種分布的隨機數(shù)的產(chǎn)生基本方法有如下三種：逆變換法合成法篩選法逆變換法設(shè)隨機變量的分布函數(shù)為，定義定理設(shè)隨機變量聽從上的勻整分布，則的分布函數(shù)為。因此，要產(chǎn)生來自的隨機數(shù)，只要先產(chǎn)生來自的隨機數(shù)，然后計算即可。其步驟為

合成法

合成法的應(yīng)用最早見于Butlter的書中。構(gòu)思如下：假如的密度函數(shù)難于抽樣，而關(guān)于的條件密度函數(shù)以及的密度函數(shù)均易于抽樣，則的隨機數(shù)可如下產(chǎn)生：可以證明由此得到的聽從。篩選抽樣

假設(shè)我們要從抽樣，假如可以將表示成，其中是一個密度函數(shù)且易于抽樣，而，是常數(shù)，則的抽樣可如下進行：定理設(shè)的密度函數(shù)，且，其中，，是一個密度函數(shù)。令和分別聽從和，則在的條件下，的條件密度為

標(biāo)準正態(tài)分布的隨機數(shù)的隨機數(shù)產(chǎn)生方法很多。簡要介紹三種。法1、變換法（Box和Muller1958）設(shè)，是獨立同分布的變量，令則與獨立，均聽從標(biāo)準正態(tài)分布。法2、結(jié)合合成法與篩選法。（略）法3、近似方法（利用中心極限定理）即用個變量產(chǎn)生一個變量。其中是抽自的隨機數(shù)，可近似為一個變量。由于正態(tài)隨機數(shù)的獨立性，當(dāng)很小時，按(a)、(b)所構(gòu)成的過程的相關(guān)時間很短，從而具有很高的截止頻率。當(dāng)然，過小，樣本的計算量過大。因此，選取適當(dāng)小即可。關(guān)于隨機數(shù)的幾點注注1由于勻整分布的隨機數(shù)的產(chǎn)生總是接受某個確定的模型進行的，從理論上講，總會有周期現(xiàn)象出現(xiàn)的。初值確定后，全部隨機數(shù)也隨之確定，并不滿足真正隨機數(shù)的要求。因此通常把由數(shù)學(xué)方法產(chǎn)生的隨機數(shù)成為偽隨機數(shù)。但其周期又相當(dāng)長，在實際應(yīng)用中幾乎不行能出現(xiàn)。因此，這種由計算機產(chǎn)生的偽隨機數(shù)可以當(dāng)作真正的隨機數(shù)來處理。注2應(yīng)對所產(chǎn)生的偽隨機數(shù)作各種統(tǒng)計檢驗，如獨立性檢驗，分布檢驗，功率譜檢驗等等。MonteCarlo方法相關(guān)引例:1、Buffon投針試驗：18世紀，法國數(shù)學(xué)家ComtedeBuffon利用投針試驗估計的值dLProblemofBuffon’sneedle:Ifaneedleoflengthlisdroppedatrandomonthemiddleofahorizontalsurfaceruledwithparallellinesadistanced>lapart,whatistheprobabilitythattheneedlewillcrossoneofthelines?Solution:Thepositioningoftheneedlerelativetonearbylinescanbedescribedwitharandomvectorwhichhascomponents:Therandomvectorisuniformlydistributedontheregion[0,d)×[0,).Accordingly,ithasprobabilitydensityfunction1/d.Theprobabilitythattheneedlewillcrossoneofthelinesisgivenbytheintegral圓周率的值π=3.

14159265358979323846264338327950288419716939937510

58209749445923078164062862089986280348253421170679

82148086513282306647093844609550582231725359408128

48111745028410270193852110555964462294895493038196

44288109756659334461284756482337867831652712019091

45648566923460348610454326648213393607260249141273

72458700660631558817488152092096282925409171536436

78925903600113305305488204665213841469519415116094

33057270365759591953092186117381932611793105118548

07446237996274956735188575272489122793818301194912

98336733624406566430860213949463952247371907021798

60943702770539217176293176752384674818467669405132

00056812714526356082778577134275778960917363717872

14684409012249534301465495853710507922796892589235

420199561121290219608640344181598136297747713.....什么是蒙特-卡洛模擬？計算機試驗物理試驗——用探測器取數(shù)據(jù)5RHIC-STAR實驗探測器蒙特-卡洛模擬

——

用計算機取數(shù)據(jù)一個例子——道爾頓板理論——分布曲線試驗——丟鐵球模擬——計算機丟球6模擬的動身點——模型依據(jù)球和釘子碰撞的規(guī)律追蹤每個球的運動依據(jù)落入每個格子的幾率給出這一幾率的每次實現(xiàn)89令n為鐵球數(shù)，m為格子數(shù)。p(i)是鐵球落入第i個格子中的概率再產(chǎn)生n個勻整分布的隨機數(shù)： ra(j),j=1,2,...,n將p(i)累加s=0.r(1)=0.do2i=1,ms=s+p(i)r(i+1)=senddodo3k=1,mnp(k)=0do3j=1,nif(ra(j).ge.r(k).and.ra(j).lt.r(k+1))[r(k)≤ra(j)<r(k+1)]np(k)=np(k)+1Enddonp(k)是落入第k個盒子中的粒子數(shù)。r(m)Σp(i)=1r(3)p(1)+p(2)r(2)p(1)r(1)0···i=1,m···通過比較第j個鐵球的ra(j)和r(i)來確定這個球落入那一格子中：例一：道爾頓板的蒙特-卡洛模擬p3p4p5p2p1p6p7例二面積的計算f(x)x辛普遜方法I=ΣSn蒙特-卡洛方法f(x)x在長方形中均勻投N0組(x,y)如y<f(x),則N=N+1I=(N/N0)×S0SS011考慮平面上的一個邊長為1的正方形及其內(nèi)部的一個形態(tài)不規(guī)則的“圖形”，如何求出這個“圖形”的面積呢？MonteCarlo方法是這樣一種“隨機化”的方法：向該正方形“隨機地”投擲N個點，若有M個點落于