計算材料學概述-之-蒙特卡洛方法.詳解_第1頁
計算材料學概述-之-蒙特卡洛方法.詳解_第2頁
計算材料學概述-之-蒙特卡洛方法.詳解_第3頁
計算材料學概述-之-蒙特卡洛方法.詳解_第4頁
計算材料學概述-之-蒙特卡洛方法.詳解_第5頁
已閱讀5頁,還剩107頁未讀 繼續免費閱讀

下載本文檔

版權說明:本文檔由用戶提供并上傳,收益歸屬內容提供方,若內容存在侵權,請進行舉報或認領

文檔簡介

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

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

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

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

3)計算針與直線相交的概率.

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

p=2l/(πd),π為圓周率:

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

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

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

1

、依據圓面積的公式:

s=πR^2

,當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的優點MC與傳統數學方法相比,具有直觀性強,簡便易行的優點,該方法能處理一些其他方法無法解決的負責問題,并且簡潔在計算機上實現,在很大程度上可以代替很多大型、難以實現的困難試驗和社會行為。無污染、無危急、能擺脫試驗誤差。MonteCarlo方法利用隨機抽樣的方法來求解物理問題;數值解法:從一個物理系統的數學模型動身,通過求解一系列的微分方程來的導出系統的未知狀態;MonteCarlo方法并非只能用來解決包含隨機過程的問題:例如:用MonteCarlo方法計算定積分.

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

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

產生一個0到1之間的隨機數(x可以是向量),但是每次總是那幾個數。

用了random_seed

()后,系統依據日期和時間隨機地供應種子,使得隨機數更隨機了。

program

random

real

::

x

call

random_seed

()

!

系統依據日期和時間隨機地供應種子

call

random_number

(x)

!

每次的隨機數就都不一樣了

write(*,*)

x

stop

end

program

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

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

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應用之二生日問題MC模擬假設有n個人在一起,各自的生日為365天之一,依據概率理論,與很多人的直覺相反,只需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在材料學領域的應用

——隨機行走背景如,布朗運動---最簡潔、無限制隨機行走(Unrestrictedrandonwalk,RW)startend平均平方端-端位移:,自然科學和社會生活中很多現象都與隨機運動有關可以模擬的內容?擴散;分子運動;。。。。。。如圖所示,第i個分子在經過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(*,*)"試驗天數JMAX,試驗次數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(*,*)"試驗天數JMAX,試驗次數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(*,*)"試驗天數JMAX,試驗次數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實際環境是困難的:可變,缺陷,風,季節,電場等。

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

氣體分子的隨機行走

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

USEMSFLIB! INTEGERXR,YR!在的區域中畫一個圓

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!像素的初始狀態(顏色)

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!在的區域中畫一個圓 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!像素的初始狀態(顏色) 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個自然數倒數和。即求程序須要幾個變量:累加的項數I(整型),存放通項值的T(實型),存放總和的S(實型),待輸入的總項數N(整型)。PROGRAMMAINIMPLICITNONEINTEGERN,IREALT,SREAD*,NS=0.0;I=1DOT=1.0/IS=S+TI=I+1IF(I>N)EXITENDDOPRINT*,"SUM=",SENDPROGRAMMAIN程序的運行結果如下15

SUM=3.182290用DO結構處理循環問題,最關鍵的是要處理好初值、循環終止條件和循環體語句的關系。FinitedifferenceapproximationofdifferentialequationsAdifferentialequationcanbeapproximatedbyafinitedifferencescheme.ForexampleForwardtimeBackwardspaceCentralspaceForwardtime-CentralspaceFICK其次定律Fick其次定律——穩態擴散解

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、相關幾個例子MonteCarlo方法:亦稱統計模擬方法,statisticalsimulationmethod

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

由大數定律,Monte-Carlo模擬的收斂是以概率而言的.誤差:

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

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

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

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

合成法

合成法的應用最早見于Butlter的書中。構思如下:假如的密度函數難于抽樣,而關于的條件密度函數以及的密度函數均易于抽樣,則的隨機數可如下產生:可以證明由此得到的聽從。篩選抽樣

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

標準正態分布的隨機數的隨機數產生方法很多。簡要介紹三種。法1、變換法(Box和Muller1958)設,是獨立同分布的變量,令則與獨立,均聽從標準正態分布。法2、結合合成法與篩選法。(略)法3、近似方法(利用中心極限定理)即用個變量產生一個變量。其中是抽自的隨機數,可近似為一個變量。由于正態隨機數的獨立性,當很小時,按(a)、(b)所構成的過程的相關時間很短,從而具有很高的截止頻率。當然,過小,樣本的計算量過大。因此,選取適當小即可。關于隨機數的幾點注注1由于勻整分布的隨機數的產生總是接受某個確定的模型進行的,從理論上講,總會有周期現象出現的。初值確定后,全部隨機數也隨之確定,并不滿足真正隨機數的要求。因此通常把由數學方法產生的隨機數成為偽隨機數。但其周期又相當長,在實際應用中幾乎不行能出現。因此,這種由計算機產生的偽隨機數可以當作真正的隨機數來處理。注2應對所產生的偽隨機數作各種統計檢驗,如獨立性檢驗,分布檢驗,功率譜檢驗等等。MonteCarlo方法相關引例:1、Buffon投針試驗:18世紀,法國數學家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.....什么是蒙特-卡洛模擬?計算機試驗物理試驗——用探測器取數據5RHIC-STAR實驗探測器蒙特-卡洛模擬

——

用計算機取數據一個例子——道爾頓板理論——分布曲線試驗——丟鐵球模擬——計算機丟球6模擬的動身點——模型依據球和釘子碰撞的規律追蹤每個球的運動依據落入每個格子的幾率給出這一幾率的每次實現89令n為鐵球數,m為格子數。p(i)是鐵球落入第i個格子中的概率再產生n個勻整分布的隨機數: 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個盒子中的粒子數。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的正方形及其內部的一個形態不規則的“圖形”,如何求出這個“圖形”的面積呢?MonteCarlo方法是這樣一種“隨機化”的方法:向該正方形“隨機地”投擲N個點,若有M個點落于

溫馨提示

  • 1. 本站所有資源如無特殊說明,都需要本地電腦安裝OFFICE2007和PDF閱讀器。圖紙軟件為CAD,CAXA,PROE,UG,SolidWorks等.壓縮文件請下載最新的WinRAR軟件解壓。
  • 2. 本站的文檔不包含任何第三方提供的附件圖紙等,如果需要附件,請聯系上傳者。文件的所有權益歸上傳用戶所有。
  • 3. 本站RAR壓縮包中若帶圖紙,網頁內容里面會有圖紙預覽,若沒有圖紙預覽就沒有圖紙。
  • 4. 未經權益所有人同意不得將文件中的內容挪作商業或盈利用途。
  • 5. 人人文庫網僅提供信息存儲空間,僅對用戶上傳內容的表現方式做保護處理,對用戶上傳分享的文檔內容本身不做任何修改或編輯,并不能對任何下載內容負責。
  • 6. 下載文件中如有侵權或不適當內容,請與我們聯系,我們立即糾正。
  • 7. 本站不保證下載資源的準確性、安全性和完整性, 同時也不承擔用戶因使用這些下載資源對自己和他人造成任何形式的傷害或損失。

評論

0/150

提交評論