計算材料學概述之蒙特卡洛方法詳解課件

1、計算材料學概述之蒙特卡洛方法詳解計算材料學概述之蒙特卡洛方法詳解主要內(nèi)容Monte Carlo模擬發(fā)展簡介Monte Carlo模擬基本原理Monte Carlo模擬典型算法Monte Carlo模擬典型應用 主要內(nèi)容Monte Carlo模擬發(fā)展簡介蒙特卡洛法是什么？蒙特卡洛(Monte Carlo)方法，是在簡單的理論準則基礎上，采用反復隨即抽樣的方法，解決復雜系統(tǒng)的問題。其實質(zhì)是一種概率和統(tǒng)計的問題。 蒙特卡羅方法（Monte Carlo method），也稱統(tǒng)計模擬方法，是二十世紀四十年代中期由于科學技術的發(fā)展和電子計算機的發(fā)明，而被提出的一種以概率統(tǒng)計理論為指導的一類非常重要的數(shù)值計

2、算方法。是指使用隨機數(shù)（或更常見的偽隨機數(shù)）來解決很多計算問題的方法。 蒙特卡洛法是什么？蒙特卡洛(Monte Carlo)方法，是MC的基本思想 MC基本思想很早以前就被人們所發(fā)現(xiàn)和利用。17世紀，人們就知道用事件發(fā)生的“頻率”來決定事件的“概率”。但要真正實現(xiàn)隨機抽樣是很困難的，甚至幾乎是不可能的。 高速計算機的出現(xiàn)，使得用數(shù)學方法在計算機上大量、快速地模擬這樣的試驗成為可能。MC的基本思想 MC基本思想很早以前就被人們所發(fā)現(xiàn)和確定性系統(tǒng)隨機性系統(tǒng)模擬自然界Monte-Carlo模擬，即隨機模擬（重復“試驗”）重復試驗計算機模擬確定性系統(tǒng)隨機性系統(tǒng)模擬自然界Monte-Carlo模擬，即M

3、onte Carlo方法：亦稱統(tǒng)計模擬方法，statistical simulation method 利用隨機數(shù)進行數(shù)值模擬的方法Monte Carlo名字的由來：是由Metropolis在二次世界大戰(zhàn)期間提出的：Manhattan計劃，研究與原子彈有關的中子輸運過程；Monte Carlo是摩納哥（monaco)的首都，該城以賭博聞名Nicholas Metropolis (1915-1999)Monte-Carlo, MonacoMonte Carlo方法：亦稱統(tǒng)計模擬方法，statistMonte Carlo方法簡史簡單地介紹一下Monte Carlo方法的發(fā)展歷史1、Buffon投針

4、實驗：18世紀，法國數(shù)學家Comte de Buffon利用投針實驗估計的值dLMonte Carlo方法簡史簡單地介紹一下Monte Ca1777年法國科學家布豐提出的一種計算圓周率的方法隨機投針法，即著名的布豐投針問題。這一方法的步驟是： 1) 取一張白紙，在上面畫上許多條間距為d的平行線。 2) 取一根長度為l（ld） 的針，隨機地向畫有平行直線的紙上擲n次，觀察針與直線相交的次數(shù)，記為m 3）計算針與直線相交的概率 布豐本人證明了，這個概率是: p=2l/(d) ，為圓周率 ：利用這個公式可以用概率的方法得到圓周率的近似值。下面是一些資料 1777年法國科學家布豐提出的一種計算圓周率的

5、方法隨機 實驗者 年代 投擲次數(shù) 相交次數(shù) 圓周率估計值 沃爾夫 1850 5000 2531 3.1596 史密斯 1855 3204 1219 3.1554 德摩根 1880 600 383 3.137 福克斯 1884 1030 489 3.1595 拉澤里尼 1901 3408 1808 3.1415929 賴納 1925 2520 859 3.1795 布豐投針實驗是第一個用幾何形式表達概率問題的例子，他首次使用隨機實驗處理確定性數(shù)學問題，為概率論和蒙特卡羅方法的發(fā)展起到一定的推動作用。 實驗者 年代 投擲次數(shù) 相交次數(shù) 圓周率估計值Monte Carlo方法之隨機數(shù)的產(chǎn)生 許多計算

6、機系統(tǒng)都有隨機數(shù)生成函數(shù)F90: call random_seed call random_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ù)到初始狀態(tài)注意：上述隨機數(shù)序列均具周期性，如上頁random子程序的周期約230。Monte Carlo方法之隨機數(shù)的產(chǎn)生F90: 實例一、計算值計算過程：1、構(gòu)造或描述問題的概率過程2、從概率密度函數(shù)出發(fā)進行隨機抽樣，實現(xiàn)從已知概率分布的抽樣，得到特征量的一些模擬結(jié)果計算均值Mon

7、te Carlo方法之典型算法與應用實例一、計算值Monte Carlo方法之典型算法與應用考慮平面上的一個邊長為1的正方形及其內(nèi)部的一個形狀不規(guī)則的“圖形”，如何求出這個“圖形”的面積呢？Monte Carlo方法是這樣一種“隨機化”的方法：向該正方形“隨機地”投擲N個點，若有M個點落于“圖形”內(nèi)，則該“圖形”的面積近似為M/N?？紤]平面上的一個邊長為1的正方形及其內(nèi)部的一個形狀不規(guī)則的“用該方法計算的基本思路是：1、根據(jù)圓面積的公式：s=R2,當R=1時，S=。2、由于圓的方程是：x2+y2=1(x2為x的平方的意思),因此1/4圓面積為x軸、y軸和上述方程所包圍的部分。3、如果在1*1的

8、正方形中均勻地落入隨機點，則落入1/4圓中的點的概率就是1/4圓的面積。其4倍，就是圓面積。由于半徑為1，該面積的值為的值。用該方法計算的基本思路是：REAL R,R1,R2,PIISEED=RTC()N0=0N=300000DO I=1,NR1=RAN(ISEED)R2=RAN(ISEED)R=SQRT(R1*R1+R2*R2)IF(R1.0)N0=N0+1END DOPI=4.0*N0/NWRITE(*,*)PIENDREAL R,R1,R2,PI面積的計算f (x)x辛普遜方法I = Sn蒙特-卡洛方法f (x)x在長方形中均勻投N0組(x,y)如 yf(x), 則 N=N+1I =(N

9、/N0)S0SS011面積的計算f (x)x辛普遜方法I = Sn蒙特-卡洛方法MC的優(yōu)點MC與傳統(tǒng)數(shù)學方法相比，具有直觀性強，簡便易行的優(yōu)點，該方法能處理一些其他方法無法解決的負責問題，并且容易在計算機上實現(xiàn)，在很大程度上可以代替許多大型、難以實現(xiàn)的復雜實驗和社會行為。無污染、無危險、能擺脫實驗誤差。Monte Carlo方法利用隨機抽樣的方法來求解物理問題;數(shù)值解法:從一個物理系統(tǒng)的數(shù)學模型出發(fā),通過求解一系列的微分方程來的導出系統(tǒng)的未知狀態(tài);Monte Carlo方法并非只能用來解決包含隨機過程的問題: 例如:用Monte Carlo方法計算定積分. 對這樣的問題可將其轉(zhuǎn)換成相關的隨機過

10、程, 然后用Monte Carlo方法進行求解注意以下兩點：MC的優(yōu)點MC與傳統(tǒng)數(shù)學方法相比，具有直觀性強，簡便易行的優(yōu)MC的應用自然現(xiàn)象的模擬：宇宙射線在地球大氣中的傳輸過程；高能物理實驗中的核相互作用過程；數(shù)值分析：數(shù)學問題，求積分，求逆矩陣，解線性代數(shù)等經(jīng)濟學模擬：庫存問題，隨機服務系統(tǒng)中排隊問題 人口問題：人口的出生，傳染病的蔓延；乃至動物的生態(tài)競爭金屬學：擴散、組織長大、相變過程MC的應用自然現(xiàn)象的模擬：宇宙射線在地球大氣中的傳輸過程；高蒙特-卡洛模擬的意義能研究不同邊界、不同材料的影響 理論不可能、實驗耗費太大 用于實驗設計無污染 反應堆防護 核彈爆炸能擺脫實驗誤差 作理論和實驗的

11、橋梁蒙特-卡洛模擬的意義能研究不同邊界、不同材料的影響Monte Carlo模擬的步驟：根據(jù)欲研究的物理系統(tǒng)的性質(zhì)，建立能夠描述該系統(tǒng)特性的理論模型，導出該模型的某些特征量的概率密度函數(shù)； （即構(gòu)造或描述問題的概率過程）從概率密度函數(shù)出發(fā)進行隨機抽樣，實現(xiàn)從已知概率分布的抽樣，得到特征量的一些模擬結(jié)果； 有了明確的概率過程后，為了實現(xiàn)過程的數(shù)值模擬，必須實現(xiàn)從已知概率分布的隨機數(shù)的抽樣，進行大量的隨機模擬實驗，從中獲得隨機變量的大量試驗值。產(chǎn)生已知概率分布的隨機變量，是實現(xiàn)MC方法的關鍵步驟，其中最基本的是（0，1）均勻分布。對模擬結(jié)果進行分析總結(jié)，預言物理系統(tǒng)的某些特性。4. 模擬結(jié)果的檢驗

12、Monte Carlo模擬的步驟：根據(jù)欲研究的物理系統(tǒng)的性質(zhì)Monte Carlo算法的主要組成部分概率密度函數(shù)(pdf) 必須給出描述一個物理系統(tǒng)的一組概率密度函數(shù);隨機數(shù)產(chǎn)生器能夠產(chǎn)生在區(qū)間0,1上（均勻）分布的隨機數(shù)抽樣規(guī)則如何從在區(qū)間0,1上均勻分布的隨機數(shù)出發(fā),隨機抽取服從給定的pdf的隨機變量;模擬結(jié)果記錄記錄一些感興趣的量的模擬結(jié)果誤差估計必須確定統(tǒng)計誤差（或方差）隨模擬次數(shù)以及其它一些量的變化；減少方差的技術利用該技術可減少模擬過程中計算的次數(shù)；并行和矢量化可以在先進的并行計算機上運行的有效算法Monte Carlo算法的主要組成部分概率密度函數(shù)(pdf實例二 定積分計算事實上

13、，不少的統(tǒng)計問題，如計算概率、各階距等，最后都歸結(jié)為定積分的近似計算問題。下面考慮一個簡單的定積分實例二 定積分計算事實上，不少的統(tǒng)計問題，如計算概率、各階距! 計算x*2在(0,1)上積分計算過程：1、構(gòu)造或描述問題的概率過程:產(chǎn)生服從分布f（x）的隨機變量Xi( )（i=1,2, ,N）2、從概率密度函數(shù)出發(fā)進行隨機抽樣，實現(xiàn)從已知概率分布的抽樣，得到特征量的一些模擬結(jié)果計算均值( )! 計算x*2在(0,1)上積分REAL YY=0N=300000ISEED=RTC()DO I=1,N X=RAN(ISEED) Y=Y+X*2/N END DOWRITE(*,*)YEND limx2dx

14、 (dx 0)REAL Ylimx2dx (dx 0)Monte Carlo方法另一個重要問題：隨機數(shù) 隨機數(shù)：由單位矩陣分布中所產(chǎn)生的簡單子樣稱為隨機數(shù)序列，其中的每一個個體稱為隨機數(shù)。 但真正的隨機數(shù)的不適合電子計算機上使用，因為它需要很大的存儲量。利用某些物理現(xiàn)象可以在電子計算機上產(chǎn)生隨機數(shù)，且其產(chǎn)生的序列無法重復實現(xiàn)，使程序無法復算，結(jié)果無法驗證，同時需要增添隨機數(shù)發(fā)生器和電路聯(lián)系等附加設備。Monte Carlo方法另一個重要問題：隨機數(shù) 隨機數(shù)：偽隨機數(shù)：是有數(shù)學遞推公式所產(chǎn)生的隨機數(shù)。（近似的具備隨機數(shù)的性質(zhì)。） An+1=T(A);An+1= An+k+1偽隨機的優(yōu)點和缺點：判

15、斷偽隨機數(shù)好壞的方法： 1、它能夠有較好的均勻性和獨立性； 2、它的費用大小，即指所消耗計算機的時間； 3、容量要求盡可能大。偽隨機數(shù)：是有數(shù)學遞推公式所產(chǎn)生的隨機數(shù)。（近似的具備隨機數(shù)隨機數(shù)產(chǎn)生的辦法產(chǎn)生均勻分布隨機數(shù)的幾種方法； （1）物理方法；（2）數(shù)學方法。偽隨機數(shù)產(chǎn)生方法：加同余法乘同余法乘加同余法取中方法逆變換法合成法篩選法 。隨機數(shù)產(chǎn)生的辦法產(chǎn)生均勻分布隨機數(shù)的幾種方法； 關于隨機數(shù)的幾點注意注1 由于均勻分布的隨機數(shù)的產(chǎn)生總是采用某個確定的模型進行的，從理論上講，總會有周期現(xiàn)象出現(xiàn)的。初值確定后，所有隨機數(shù)也隨之確定，并不滿足真正隨機數(shù)的要求。因此通常把由數(shù)學方法產(chǎn)生的隨機數(shù)成

16、為偽隨機數(shù)。 但其周期又相當長，在實際應用中幾乎不可能出現(xiàn)。因此，這種由計算機產(chǎn)生的偽隨機數(shù)可以當作真正的隨機數(shù)來處理。注2 應對所產(chǎn)生的偽隨機數(shù)作各種統(tǒng)計檢驗，如獨立性檢驗，分布檢驗，功率譜檢驗等等。 關于隨機數(shù)的幾點注意注1 由于均勻分布的隨機數(shù)的產(chǎn)生總是采用FORTRAN 語言產(chǎn)生隨機數(shù)的實例random_number(x)產(chǎn)生一個0到1之間的隨機數(shù)（x可以是向量），但是每次總是那幾個數(shù)。用了random_seed()后，系統(tǒng)根據(jù)日期和時間隨機地提供種子，使得隨機數(shù)更隨機了。programrandomreal:xcallrandom_seed()!系統(tǒng)根據(jù)日期和時間隨機地提供種子call

17、random_number(x)!每次的隨機數(shù)就都不一樣了write(*,*)xstopendprogramrandom FORTRAN 語言產(chǎn)生隨機數(shù)的實例random_numbeMonte Carlo方法之隨機數(shù)的產(chǎn)生 許多計算機系統(tǒng)都有隨機數(shù)生成函數(shù)F90: call random_seed call random_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ù)到初始狀態(tài)注意：上述隨機數(shù)序列均具周期性，如上頁ra

18、ndom子程序的周期約230。Monte Carlo方法之隨機數(shù)的產(chǎn)生F90: Finite difference approximation of differential equationsA differential equation can be approximated by a finite difference scheme. For exampleForward timeBackward spaceCentral spaceForward t aime-Central spaceFinite difference approximatioFICK 第二定律FICK 第二定律Fi

19、ck 第二定律穩(wěn)態(tài)擴散解 REAL C0(0:1000+1),C(0:1000+1) !C(DISTANCE)C0=0.1C0(0)=0.8 !BOUNDARY CONDITIONC0(1001)=0.1 !BOUNDARY CONDITIONC=C0 OPEN(1,FILE=F:DIF.dat) DO JT=1,50000 ！Time DO IX=1,1000 ! distanceC(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) END D

20、O C0=C END DOENDFick 第二定律穩(wěn)態(tài)擴散解 REAL C0(0:應用之二 生日問題MC模擬假設有n個人在一起，各自的生日為365天之一，根據(jù)概率理論，與很多人的直覺相反，只需23個人便有大于50的幾率人群中至少有2個人生日相同。n 理論幾率 模擬幾率0.117 0.1100.411 0.4120.527 0.5200.706 0.6920.941 0.93650 0.986 0.987應用之二 生日問題MC模擬假設有n個人在一起，各自的生日為INTEGER M(1:10000), NUMBER1(0:364), NUMBER2 REAL X,Y ISEED=RTC()DO J

21、=1,10000NUMBER1=0X=RAN(ISEED)NUMBER1(0)=INT(365*X+1)JJJ=1 DO I=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)=NUMBER2 END IF END DOEND DO DO I=1,10000 IF(M(I).LE.23) SUM=SUM+1END DO PRINT *,SUM/10000 END INTEGER M(1:10000

22、), NUMBER1(0MC在材料學領域的應用 隨機行走MC在材料學領域的應用 隨機行走背景如，布朗運動最簡單、無限制隨機行走(Unrestricted randon walk,RW )startend平均平方端-端位移:,自然科學和社會生活中很多現(xiàn)象都與隨機運動有關可以模擬的內(nèi)容？擴散；分子運動；。背景如，布朗運動最簡單、無限制隨機行走(Unrestr如圖所示，第i個分子在經(jīng)過N步隨機行走后距原點距離為R，對n個分子每步的位移平方求和后取平均值就得到了所有分子距原點的方均距離：如圖所示，第i個分子在經(jīng)過N步隨機行走后距原點距離為R，對計算材料學概述之蒙特卡洛方法詳解計算材料學概述之蒙特卡洛方

23、法詳解！Monte Carlo Simulation of One Dimensional Diffusion INTEGER X,XX(1:1000,1:1000)REAL XXM(1:1000)!X:INSTANTANEOUS POSITION OF ATOM!XX(J,I)：X*X ,J:第幾天實驗,I:第幾步跳躍!XXM(I): THE MEAN OF XX WRITE(*,*) 實驗天數(shù)JMAX，實驗次數(shù) IMAXREAD(*,*) JMAX,IMAXISEED=RTC()DO J=1,JMAX !第幾天實驗X=0 !DO I=1,IMAX !第幾步跳躍RN=RAN(ISEED)IF

24、(RN0.5)THEN X=X+1ELSE X=X-1 END IF XX(J,I)=X*XEND DOEND DOOPEN(1,FILE=“f:DIF1.DAT)DO I=1,IMAXXXM=0.0XXM(I)=1.0*SUM(XX(1:JMAX,I)/JMAX !WRITE(1,*) I, XXM(I)END DOCLOSE(1)END！Monte Carlo Simulation of On!Monte Carlo Simulation of Two Dimensional Diffusion INTEGER X,Y,XY(1:1000,1:1000)REAL XYM(1:1000)!X

25、:INSTANTANEOUS POSITION OF ATOM!XY(J,I)：X*Y ,J:第幾天實驗,I:第幾步跳躍!XYM(I): THE MEAN OF XY WRITE(*,*) 實驗天數(shù)JMAX，實驗次數(shù) IMAXREAD(*,*) JMAX,IMAXISEED=RTC()DO J=1,JMAX !第幾天實驗X=0 !Y=0 !DO I=1,IMAX !第幾步跳躍RN=RAN(ISEED)IF(RN.LT.0.25)THEN x=xy=y-1 END IFIF(RN.LT.0.5.AND.RN.GE.0.25)THEN x=xy=y+1 END IFIF(RN.LT.0.75.AN

26、D.RN.GE.0.5)THEN x=x-1y=y END IFIF(RN.GE.0.75)THEN x=x+1y=y END IFXY(J,I)=X*X+Y*YEND DOEND DOOPEN(1,FILE=“f:DIF2.DAT)DO I=1,IMAXXYM=0.0XYM(I)=1.0*SUM(XY(1:JMAX,I)/JMAX !WRITE(1,*) I, XYM(I)END DOCLOSE(1)END!Monte Carlo Simulation of Tw!Monte Carlo Simulation of Two Dimensional Diffusion INTEGER X,XY

27、(1:1000,1:1000),y,XN(1:4),YN(1:4),RNREAL XYM(1:1000)!X:INSTANTANEOUS POSITION OF ATOM!XY(J,I)：X*Y ,J:第幾天實驗,I:第幾步跳躍!XYM(I): THE MEAN OF XY WRITE(*,*) 實驗天數(shù)JMAX，實驗次數(shù) IMAXREAD(*,*) JMAX,IMAXXN=(/0,0,-1,1/)YN=(/-1,1,0,0/)ISEED=RTC()DO J=1,JMAX !第幾天實驗X=0 !Y=0 !DO I=1,IMAX !第幾步跳躍RN=4*RAN(ISEED)+1X=X+XN(RN)

28、 Y=Y+YN(RN)XY(J,I)=X*X+Y*YEND DOEND DOOPEN(1,FILE=C:DIF2.DAT)DO I=1,IMAXXYM=0.0XYM(I)=1.0*SUM(XY(1:JMAX,I)/JMAX !WRITE(1,*) I, XYM(I)END DOCLOSE(1) ！做三維空間隨機行走？END!Monte Carlo Simulation of Tw實際環(huán)境是復雜的:可變，缺陷，風， 季節(jié)，電場等。 空間中有隨機性缺陷不退行走自回避行走實際環(huán)境是復雜的:可變，缺陷，風，空間中有隨機性缺陷不退行走例如模擬高分子的位形 用隨機行走方法模擬高分子位形是用隨機行走的軌跡代

29、表高分子的位形，行走過的位置代表的是構(gòu)成分子的原子或官能團，因此，無限制隨機行走忽略了體斥效應。 不退行走就是禁止在每一步行走后立即倒退，可以解決剛走的一步與上一步重疊的問題。但不退行走沒有完全解決高分子的體斥效應。 自回避行走就是所有已走過的位置不能再走，這樣就完全解決了體斥效應問題。例如模擬高分子的位形氣體分子的隨機行走 假設在一個空曠封閉的房間中心滴上一滴香水，揮發(fā)的香水分子隨機的與空氣中的粒子發(fā)生碰撞，最終會擴散到整個房間。這里我們將通過計算機模擬400個分子在二維平面內(nèi)的隨機運動 氣體分子的隨機行走 假設在一個空曠封閉的房間中心滴上一滴計算材料學概述之蒙特卡洛方法詳解討論熵變現(xiàn)在來討

30、論一下熵在我們這個氣體擴散模型中是什么意思。在剛開始的時候，我們的系統(tǒng)處在一個非常有序的狀態(tài)：所有的分子都處于原點。這時系統(tǒng)的熵為零，系統(tǒng)不存在任何的無序度。隨著時間的推移，分子開始不斷的向外擴散，這時系統(tǒng)出現(xiàn)了無序狀態(tài)，熵開始逐漸增加。 討論熵變現(xiàn)在來討論一下熵在我們這個氣體擴散模型中是什么意思。系統(tǒng)的熵和我們預測的一樣在隨時間增大，但當所有分子布滿整個邊界以內(nèi)區(qū)域時，系統(tǒng)的熵開始趨于穩(wěn)定。從這一結(jié)果中我們可以了解：當所有的分子隨機的布滿整個區(qū)域時，雖然當我們跟蹤某一個確定分子時，它還是在區(qū)域內(nèi)到處亂竄，但每個小區(qū)域內(nèi)分子的密度卻不會再變化了。所以，一旦氣體分子擴散到整個區(qū)域以后，不管我們再

31、等上多少時間，系統(tǒng)的熵都不會再有太大的起伏。換句話說，讓系統(tǒng)自動回到開始的狀態(tài)，即所有分子都在原點的狀態(tài)，已經(jīng)不可能了。 系統(tǒng)的熵和我們預測的一樣在隨時間增大，但當所有分子布滿整個計算材料學概述之蒙特卡洛方法詳解計算材料學概述之蒙特卡洛方法詳解畫一個圓晶粒USE MSFLIB !INTEGER XR,YR !在的區(qū)域中畫一個圓PARAMETER XR=400,YR=400 INTEGER R,S(1:XR,1:YR) X0=XR/2 ! 圓心位置X0，YOY0=YR/2R=MIN(X0-10,Y0-10) !圓半徑S=0 !像素的初始狀態(tài)（顏色）DO I=1,XRDO J=1,YRIF(I-X

32、0)*2+(J-Y0)*2=R*2)S(I,J)=10IER=SETCOLOR(S(I,J)+3) IER=SETPIXEL(I,J)END DOEND DOEND畫一個圓晶粒畫晶界!畫一個圓USE MSFLIB INTEGER XR，YR !在的區(qū)域中畫一個圓PARAMETER XR=400,YR=400 INTEGER R,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，YOY0=YR/2R=MIN(X0-10,Y0-10) !圓半徑S=0 !像素的初始狀態(tài)（顏色）DO

33、I=1,XRDO J=1,YRIF(I-X0)*2+(J-Y0)*20)THENIER=SETCOLOR(9)ELSEIER=SETCOLOR(8)END IF IER=SETPIXEL(I,J)END DOEND DOEND如何畫有一定寬度的晶界？畫晶界!畫一個圓例題例 求前N個自然數(shù)倒數(shù)和。即求 程序需要幾個變量：累加的項數(shù)I(整型)，存放通項值的T(實型)，存放總和的S(實型)，待輸入的總項數(shù)N(整型)。PROGRAM MAIN IMPLICIT NONE INTEGER N,I REAL T,S READ *,N S=0.0; I=1 DO T=1.0/I S=S+T I=I+1 IF

34、(IN)EXIT END DO PRINT *,SUM=,SEND PROGRAM MAIN 程序的運行結(jié)果如下15 SUM= 3.182290 用DO結(jié)構(gòu)處理循環(huán)問題，最關鍵的是要處理好初值、循環(huán)終止條件和循環(huán)體語句的關系。例題例 求前N個自然數(shù)倒數(shù)和。即求 程序需要幾個變量Finite difference approximation of differential equationsA differential equation can be approximated by a finite difference scheme. For exampleForward timeBackwa

35、rd spaceCentral spaceForward time-Central spaceFinite difference approximatioFICK 第二定律FICK 第二定律Fick 第二定律穩(wěn)態(tài)擴散解 REAL C0(0:1000+1),C(0:1000+1) !C(DISTANCE)C0=0.1C0(0)=0.8 !BOUNDARY CONDITIONC0(1001)=0.1 !BOUNDARY CONDITIONC=C0 OPEN(1,FILE=F:DIF.dat) DO JT=1,50000 ！Time DO IX=1,1000 ! distanceC(IX)=C0(I

36、X)+0.45*(C0(IX+1)+C0(IX-1)-2.0*C0(IX)C(0)=0.8C(1001)=0.1 IF(JT=50000)WRITE(1,*) IX,C(IX) END DO C0=C END DOPRINT *, PROGRAMME IS FINISHED CLOSE(1) ENDFick 第二定律穩(wěn)態(tài)擴散解 REAL C0(0:蒙特卡羅(Monte Carlo)方法1、簡介2、相關幾個例子蒙特卡羅(Monte Carlo)方法1、簡介Monte Carlo方法：亦稱統(tǒng)計模擬方法，statistical simulation method 利用隨機數(shù)進行數(shù)值模擬的方法Mont

37、e Carlo名字的由來：是由Metropolis在二次世界大戰(zhàn)期間提出的：Manhattan計劃，研究與原子彈有關的中子輸運過程；Monte Carlo是摩納哥（monaco)的首都，該城以賭博聞名Nicholas Metropolis (1915-1999)Monte-Carlo, MonacoMonte Carlo方法：亦稱統(tǒng)計模擬方法，statistMonte Carlo方法的基本思想很早以前就被人們所發(fā)現(xiàn)和利用。早在17世紀，人們就知道用事件產(chǎn)生的“頻率”來近似事件的“概 率”。 18世紀人們用投針試驗的方法來決定圓周率。本世紀40年代電子計算機的出現(xiàn)，特別是近年來高速電子計算機的出

38、現(xiàn)，使得用數(shù)學方法在計算機上大量、快速地模擬這樣的試驗成為可能。Monte Carlo方法的基本思想很早以前就被人們所發(fā)現(xiàn)和1930年，Enrico Fermi利用Monte Carlo方法研究中子的擴散，并設計了一個Monte Carlo機械裝置，F(xiàn)ermiac,用于計算核反應堆的臨界狀態(tài)Von Neumann是Monte Carlo方法的正式奠基者,他與Stanislaw Ulam合作建立了概率密度函數(shù)、反累積分布函數(shù)的數(shù)學基礎，以及偽隨機數(shù)產(chǎn)生器。在這些工作中， Stanislaw Ulam意識到了數(shù)字計算機的重要性合作起源于Manhattan工程：利用ENIAC(Electronic

39、Numerical Integrator and Computer)計算產(chǎn)額1930年，Enrico Fermi利用Monte CarlMonte Carlo模擬的應用：自然現(xiàn)象的模擬：宇宙射線在地球大氣中的傳輸過程；高能物理實驗中的核相互作用過程；實驗探測器的模擬數(shù)值分析：利用Monte Carlo方法求積分材料學中的應用：擴散、晶粒的長大、再結(jié)晶等Monte Carlo模擬的應用：自然現(xiàn)象的模擬：宇宙射線在Monte Carlo模擬在物理研究中的作用Monte Carlo模擬在物理研究中的作用Monte Carlo模擬的步驟：根據(jù)欲研究的物理系統(tǒng)的性質(zhì)，建立能夠描述該系統(tǒng)特性的理論模型，導

40、出該模型的某些特征量的概率密度函數(shù)；從概率密度函數(shù)出發(fā)進行隨機抽樣，得到特征量的一些模擬結(jié)果；對模擬結(jié)果進行分析總結(jié)，預言物理系統(tǒng)的某些特性。Monte Carlo模擬的步驟：根據(jù)欲研究的物理系統(tǒng)的性質(zhì)注意以下兩點：Monte Carlo方法與數(shù)值解法的不同:Monte Carlo方法利用隨機抽樣的方法來求解物理問題;數(shù)值解法:從一個物理系統(tǒng)的數(shù)學模型出發(fā),通過求解一系列的微分方程來的導出系統(tǒng)的未知狀態(tài);Monte Carlo方法并非只能用來解決包含隨機的過程的問題:許多利用Monte Carlo方法進行求解的問題中并不包含隨機過程 例如:用Monte Carlo方法計算定積分. 對這樣的問題

41、可將其轉(zhuǎn)換成相關的隨機過程, 然后用Monte Carlo方法進行求解注意以下兩點：Monte Carlo方法與數(shù)值解法的不同:MMonte Carlo算法的主要組成部分概率密度函數(shù) 必須給出描述一個物理系統(tǒng)的一組概率密度函數(shù);隨機數(shù)產(chǎn)生器能夠產(chǎn)生在區(qū)間0,1上均勻分布的隨機數(shù)抽樣規(guī)則如何從在區(qū)間0,1上均勻分布的隨機數(shù)出發(fā),隨機抽取服從給定的pdf的隨機變量;模擬結(jié)果記錄記錄一些感興趣的量的模擬結(jié)果誤差估計必須確定統(tǒng)計誤差（或方差）隨模擬次數(shù)以及其它一些量的變化；減少方差的技術利用該技術可減少模擬過程中計算的次數(shù)；并行和矢量化可以在先進的并行計算機上運行的有效算法Monte Carlo算法的

42、主要組成部分概率密度函數(shù) 必須確定性系統(tǒng)隨機性系統(tǒng)模擬自然界Monte-Carlo模擬，即隨機模擬（重復“試驗”）重復試驗計算機模擬確定性系統(tǒng)隨機性系統(tǒng)模擬自然界Monte-Carlo模擬，即蒙特卡羅方法的基本原理及思想當所要求解的問題是某種事件出現(xiàn)的概率，或者是某個隨機變量的期望值時，它們可以通過某種“試驗”的方法，得到這種事件出現(xiàn)的頻率，或者這個隨機變數(shù)的平均值，并用它們作為問題的解。把蒙特卡羅解題歸結(jié)為三個主要步驟： 構(gòu)造或描述概率過程； 實現(xiàn)從已知概率分布抽樣； 建立各種估計量。 蒙特卡羅方法的基本原理及思想當所要求解的問題是某種事件出現(xiàn)的對于本身就具有隨機性質(zhì)的問題，主要是正確描述和

43、模擬這個概率過程；對于不是隨機性質(zhì)的確定性問題，比如計算定積分，就必須事先構(gòu)造一個人為的概率過程。最簡單、最基本、最重要的一個概率分布是(0,1)上的均勻分布（或稱矩形分布）。 實現(xiàn)從已知概率分布抽樣(模擬試驗)：構(gòu)造了概率模型以后，產(chǎn)生已知概率分布的隨機變量（或隨機向量），就成為實現(xiàn)蒙特卡羅方法模擬實驗的基本手段，這也是蒙特卡羅方法被稱為隨機抽樣的原因。建立各種估計量： 構(gòu)造了概率模型并能從中抽樣后，即實現(xiàn)模擬實驗后，就要確定一個隨 機變量，作為所要求的問題的解，我們稱它為無偏估計。對于本身就具有隨機性質(zhì)的問題，主要是正確描述和模擬這個概率過收斂性: 由大數(shù)定律, Monte-Carlo模擬

44、的收斂是以概率而言的誤差: 用頻率估計概率時誤差的估計，可由中心極限定理，給定置信水平 的條件下，有： 模擬次數(shù):由誤差公式得收斂性: 由大數(shù)定律, Monte-Carlo模擬的收斂是以隨機數(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)生方法。最后，關于隨機數(shù)的幾點注。隨機數(shù)的產(chǎn)生原理-計算機抽樣實驗 隨機數(shù)的產(chǎn)生是實現(xiàn)MC計(0,1均勻分布隨機數(shù)的產(chǎn)生最常用的是在(0,1區(qū)間內(nèi)均

45、勻分布的隨機數(shù)，其他分布的隨機數(shù)可利用均勻分布隨機數(shù)來產(chǎn)生。 均勻分布隨機數(shù)的產(chǎn)生：乘同余法 產(chǎn)生區(qū)間(0,1內(nèi)均勻分布的隨機數(shù)的遞推公式是： 其中，是乘因子，M是模數(shù)。第一式稱做以M為模數(shù)（modulus）的同余式，即以M 除后得到的余數(shù)記為。當給定了一個初值（稱為種子）， 計算出的 即在(0,1上均勻分布的隨機數(shù)。(0,1均勻分布隨機數(shù)的產(chǎn)生最常用的是在(0,1區(qū)間內(nèi)均例1：例1：R=RAND()ISEED=RTC() R=RAN(ISEED)R=RAND()其他各種分布的隨機數(shù)的產(chǎn)生基本方法有如下三種： 逆變換法 合成法 篩選法 其他各種分布的隨機數(shù)的產(chǎn)生基本方法有如下三種：逆變換法設隨

46、機變量 的分布函數(shù)為 ，定義 定理 設隨機變量 服從 上的均勻分布，則 的分布函數(shù)為 。因此，要產(chǎn)生來自 的隨機數(shù)，只要先產(chǎn)生來自 的隨機數(shù)，然后計算 即可。其步驟為 逆變換法設隨機變量 的分布函數(shù)為 ，定義 合成法 合成法的應用最早見于Butlter 的書中。構(gòu)思如下： 如果 的密度函數(shù) 難于抽樣，而 關于 的條件密度函數(shù) 以及 的密度函數(shù) 均易于抽樣，則 的隨機數(shù)可如下產(chǎn)生：可以證明由此得到 的服從 。合成法 合成法的應用最早見于Butlter 的書中。篩選抽樣 假設我們要從 抽樣，如果可以將 表示成 ，其中 是一個密度函數(shù)且易于抽樣，而 ， 是常數(shù)，則 的抽樣可如下進行： 定理 設 的密

47、度函數(shù) ，且 ，其中 ， ， 是一個密度函數(shù)。令 和 分別服從 和 ，則在 的條件下， 的條件密度為 篩選抽樣 假設我們要從 抽樣，如果可以將 表示成 標準正態(tài)分布的隨機數(shù) 的隨機數(shù)產(chǎn)生方法很多。簡要介紹三種。 法1、 變換法（Box 和Muller 1958）設 ， 是獨立同分布的 變量，令 則 與 獨立，均服從標準正態(tài)分布。法2、 結(jié)合合成法與篩選法。（略）法3、 近似方法（利用中心極限定理）即用 個 變量產(chǎn)生一個 變量。其中 是抽自 的隨機數(shù)， 可近似為一個 變量。標準正態(tài)分布的隨機數(shù) 的隨機數(shù)產(chǎn)生方法很多。簡要介紹三由于正態(tài)隨機數(shù)的獨立性，當 很小時，按(a)、(b)所構(gòu)成的過程的相關

48、時間很短，從而具有很高的截止頻率。當然， 過小，樣本的計算量過大。因此，選取 適當小即可。由于正態(tài)隨機數(shù)的獨立性，當 很小時，按(a)、(b)關于隨機數(shù)的幾點注注1 由于均勻分布的隨機數(shù)的產(chǎn)生總是采用某個確定的模型進行的，從理論上講，總會有周期現(xiàn)象出現(xiàn)的。初值確定后，所有隨機數(shù)也隨之確定，并不滿足真正隨機數(shù)的要求。因此通常把由數(shù)學方法產(chǎn)生的隨機數(shù)成為偽隨機數(shù)。 但其周期又相當長，在實際應用中幾乎不可能出現(xiàn)。因此，這種由計算機產(chǎn)生的偽隨機數(shù)可以當作真正的隨機數(shù)來處理。注2 應對所產(chǎn)生的偽隨機數(shù)作各種統(tǒng)計檢驗，如獨立性檢驗，分布檢驗，功率譜檢驗等等。 關于隨機數(shù)的幾點注注1 由于均勻分布的隨機數(shù)的

49、產(chǎn)生總是采用某Monte Carlo方法相關引例:1、Buffon投針實驗：18世紀，法國數(shù)學家Comte de Buffon利用投針實驗估計的值dLMonte Carlo方法相關引例:1、Buffon投針實驗Problem of Buffons needle:If a needle of length l is dropped at random on the middle of a horizontal surface ruled with parallel lines a distance dl apart, what is the probability that the needle

50、 will cross one of the lines?Problem of Buffons needle:If Solution:The positioning of the needle relative to nearby lines can be described with a random vector which has components:The random vector is uniformly distributed on the region 0,d)0,). Accordingly, it has probability density function 1/d.

51、The probability that the needle will cross one of the lines is given by the integralSolution:The positioning of th計算材料學概述之蒙特卡洛方法詳解計算材料學概述之蒙特卡洛方法詳解計算材料學概述之蒙特卡洛方法詳解圓周率的值 = 3.14159 26535 89793 23846 26433 83279 50288 41971 69399 3751058209 74944 59230 78164 06286 20899 86280 34825 34211 7067982148 0865

52、1 32823 06647 09384 46095 50582 23172 53594 0812848111 74502 84102 70193 85211 05559 64462 29489 54930 3819644288 10975 66593 34461 28475 64823 37867 83165 27120 1909145648 56692 34603 48610 45432 66482 13393 60726 02491 4127372458 70066 06315 58817 48815 20920 96282 92540 91715 3643678925 90360 011

53、33 05305 48820 46652 13841 46951 94151 1609433057 27036 57595 91953 09218 61173 81932 61179 31051 1854807446 23799 62749 56735 18857 52724 89122 79381 83011 9491298336 73362 44065 66430 86021 39494 63952 24737 19070 2179860943 70277 05392 17176 29317 67523 84674 81846 76694 0513200056 81271 45263 56

54、082 77857 71342 75778 96091 73637 1787214684 40901 22495 34301 46549 58537 10507 92279 68925 8923542019 95611 21290 21960 86403 44181 59813 62977 47713 . 圓周率的值 = 3.14159 26535 89793什么是蒙特-卡洛模擬？計算機實驗物理實驗 用探測器取數(shù)據(jù)5RHIC-STAR實驗探測器什么是蒙特-卡洛模擬？計算機實驗物理實驗 用探測器取數(shù) 蒙特-卡洛模擬 用計算機取數(shù)據(jù)計算材料學概述之蒙特卡洛方法詳解一個例子 道爾頓板理論 分布曲線實

55、驗 丟鐵球模擬 計算機丟球6一個例子 道爾頓板理論 分布曲線6模擬的出發(fā)點 模型根據(jù)球和釘子碰撞的規(guī)律 追蹤每個球的運動根據(jù)落入每個格子的幾率 給出這一幾率的每次實現(xiàn)8模擬的出發(fā)點 模型根據(jù)球和釘子碰撞的規(guī)律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.do 2 i=1, ms=s+p (i)r(i+1)=senddodo 3 k=1, mnp(k)=0do 3 j=1,nif(ra( j ).ge.r(k).and.ra( j ).lt.r(k+

56、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) 0i=1,m通過比較第 j 個鐵球的 ra(j) 和 r(i) 來決定這個球落入那一格子中：例一：道爾頓板的蒙特-卡洛模擬p3p4p5p2p1p6p79令 n 為鐵球數(shù)，m 為格子數(shù)。再產(chǎn)生 n 個均勻分布的隨例二 面積的計算f (x)x辛普遜方法I = Sn蒙特-卡洛方法f (x)x在長方形中均勻投N0組(x,y)如 yf(x), 則 N=N+1I =(N/N0)S0SS011例二 面積的計算f (x)x辛普遜方法I = Sn蒙特-考慮平面上的一個邊長為1的正方形及其內(nèi)部的一個形狀不規(guī)則的“圖形”，如何求出這個“圖形”的面積呢？Monte Carlo方法是這樣一種“隨機化”的方法：向該正方形“隨機地”投擲N個點，若有M個點落于“圖形”內(nèi)，則該“圖形”的面積近似為M/N。考慮平面上的一個邊長為1的正方形及其內(nèi)部的一個形狀不規(guī)則的“用該方法計算的基本思路是：1根據(jù)圓面積的公式：s=R2,當R=1時，S=。由于圓的方程是：x2+y2=1(x2為x的平方的意思),