計(jì)算材料學(xué)概述之蒙特卡洛方法詳解課件

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

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

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

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

5、方法隨機(jī) 實(shí)驗(yàn)者 年代 投擲次數(shù) 相交次數(shù) 圓周率估計(jì)值 沃爾夫 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í)驗(yàn)是第一個(gè)用幾何形式表達(dá)概率問(wèn)題的例子，他首次使用隨機(jī)實(shí)驗(yàn)處理確定性數(shù)學(xué)問(wèn)題，為概率論和蒙特卡羅方法的發(fā)展起到一定的推動(dòng)作用。 實(shí)驗(yàn)者 年代 投擲次數(shù) 相交次數(shù) 圓周率估計(jì)值Monte Carlo方法之隨機(jī)數(shù)的產(chǎn)生 許多計(jì)算

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

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

8、正方形中均勻地落入隨機(jī)點(diǎn)，則落入1/4圓中的點(diǎn)的概率就是1/4圓的面積。其4倍，就是圓面積。由于半徑為1，該面積的值為的值。用該方法計(jì)算的基本思路是：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面積的計(jì)算f (x)x辛普遜方法I = Sn蒙特-卡洛方法f (x)x在長(zhǎng)方形中均勻投N0組(x,y)如 yf(x), 則 N=N+1I =(N

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

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

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

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

13、，不少的統(tǒng)計(jì)問(wèn)題，如計(jì)算概率、各階距等，最后都?xì)w結(jié)為定積分的近似計(jì)算問(wèn)題。下面考慮一個(gè)簡(jiǎn)單的定積分實(shí)例二 定積分計(jì)算事實(shí)上，不少的統(tǒng)計(jì)問(wèn)題，如計(jì)算概率、各階距! 計(jì)算x*2在(0,1)上積分計(jì)算過(guò)程：1、構(gòu)造或描述問(wèn)題的概率過(guò)程:產(chǎn)生服從分布f（x）的隨機(jī)變量Xi( )（i=1,2, ,N）2、從概率密度函數(shù)出發(fā)進(jìn)行隨機(jī)抽樣，實(shí)現(xiàn)從已知概率分布的抽樣，得到特征量的一些模擬結(jié)果計(jì)算均值( )! 計(jì)算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方法另一個(gè)重要問(wèn)題：隨機(jī)數(shù) 隨機(jī)數(shù)：由單位矩陣分布中所產(chǎn)生的簡(jiǎn)單子樣稱為隨機(jī)數(shù)序列，其中的每一個(gè)個(gè)體稱為隨機(jī)數(shù)。 但真正的隨機(jī)數(shù)的不適合電子計(jì)算機(jī)上使用，因?yàn)樗枰艽蟮拇鎯?chǔ)量。利用某些物理現(xiàn)象可以在電子計(jì)算機(jī)上產(chǎn)生隨機(jī)數(shù)，且其產(chǎn)生的序列無(wú)法重復(fù)實(shí)現(xiàn)，使程序無(wú)法復(fù)算，結(jié)果無(wú)法驗(yàn)證，同時(shí)需要增添隨機(jī)數(shù)發(fā)生器和電路聯(lián)系等附加設(shè)備。Monte Carlo方法另一個(gè)重要問(wèn)題：隨機(jī)數(shù) 隨機(jī)數(shù)：偽隨機(jī)數(shù)：是有數(shù)學(xué)遞推公式所產(chǎn)生的隨機(jī)數(shù)。（近似的具備隨機(jī)數(shù)的性質(zhì)。） An+1=T(A);An+1= An+k+1偽隨機(jī)的優(yōu)點(diǎn)和缺點(diǎn)：判

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

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

17、random_number(x)!每次的隨機(jī)數(shù)就都不一樣了write(*,*)xstopendprogramrandom FORTRAN 語(yǔ)言產(chǎn)生隨機(jī)數(shù)的實(shí)例random_numbeMonte Carlo方法之隨機(jī)數(shù)的產(chǎn)生 許多計(jì)算機(jī)系統(tǒng)都有隨機(jī)數(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)間隨機(jī)分布的N*N矩陣 s= rand(state,0) 重設(shè)該生成函數(shù)到初始狀態(tài)注意：上述隨機(jī)數(shù)序列均具周期性，如上頁(yè)ra

18、ndom子程序的周期約230。Monte Carlo方法之隨機(jī)數(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)擴(kuò)散解 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)擴(kuò)散解 REAL C0(0:應(yīng)用之二 生日問(wèn)題MC模擬假設(shè)有n個(gè)人在一起，各自的生日為365天之一，根據(jù)概率理論，與很多人的直覺(jué)相反，只需23個(gè)人便有大于50的幾率人群中至少有2個(gè)人生日相同。n 理論幾率 模擬幾率0.117 0.1100.411 0.4120.527 0.5200.706 0.6920.941 0.93650 0.986 0.987應(yīng)用之二 生日問(wèn)題MC模擬假設(shè)有n個(gè)人在一起，各自的生日為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在材料學(xué)領(lǐng)域的應(yīng)用 隨機(jī)行走M(jìn)C在材料學(xué)領(lǐng)域的應(yīng)用 隨機(jī)行走背景如，布朗運(yùn)動(dòng)最簡(jiǎn)單、無(wú)限制隨機(jī)行走(Unrestricted randon walk,RW )startend平均平方端-端位移:,自然科學(xué)和社會(huì)生活中很多現(xiàn)象都與隨機(jī)運(yùn)動(dòng)有關(guān)可以模擬的內(nèi)容？擴(kuò)散；分子運(yùn)動(dòng)；。背景如，布朗運(yùn)動(dòng)最簡(jiǎn)單、無(wú)限制隨機(jī)行走(Unrestr如圖所示，第i個(gè)分子在經(jīng)過(guò)N步隨機(jī)行走后距原點(diǎn)距離為R，對(duì)n個(gè)分子每步的位移平方求和后取平均值就得到了所有分子距原點(diǎn)的方均距離：如圖所示，第i個(gè)分子在經(jīng)過(guò)N步隨機(jī)行走后距原點(diǎn)距離為R，對(duì)計(jì)算材料學(xué)概述之蒙特卡洛方法詳解計(jì)算材料學(xué)概述之蒙特卡洛方

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:第幾天實(shí)驗(yàn),I:第幾步跳躍!XXM(I): THE MEAN OF XX WRITE(*,*) 實(shí)驗(yàn)天數(shù)JMAX，實(shí)驗(yàn)次數(shù) IMAXREAD(*,*) JMAX,IMAXISEED=RTC()DO J=1,JMAX !第幾天實(shí)驗(yàn)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:第幾天實(shí)驗(yàn),I:第幾步跳躍!XYM(I): THE MEAN OF XY WRITE(*,*) 實(shí)驗(yàn)天數(shù)JMAX，實(shí)驗(yàn)次數(shù) IMAXREAD(*,*) JMAX,IMAXISEED=RTC()DO J=1,JMAX !第幾天實(shí)驗(yàn)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:第幾天實(shí)驗(yàn),I:第幾步跳躍!XYM(I): THE MEAN OF XY WRITE(*,*) 實(shí)驗(yàn)天數(shù)JMAX，實(shí)驗(yàn)次數(shù) IMAXREAD(*,*) JMAX,IMAXXN=(/0,0,-1,1/)YN=(/-1,1,0,0/)ISEED=RTC()DO J=1,JMAX !第幾天實(shí)驗(yàn)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) ！做三維空間隨機(jī)行走？END!Monte Carlo Simulation of Tw實(shí)際環(huán)境是復(fù)雜的:可變，缺陷，風(fēng)， 季節(jié)，電場(chǎng)等。 空間中有隨機(jī)性缺陷不退行走自回避行走實(shí)際環(huán)境是復(fù)雜的:可變，缺陷，風(fēng)，空間中有隨機(jī)性缺陷不退行走例如模擬高分子的位形 用隨機(jī)行走方法模擬高分子位形是用隨機(jī)行走的軌跡代

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

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

31、等上多少時(shí)間，系統(tǒng)的熵都不會(huì)再有太大的起伏。換句話說(shuō)，讓系統(tǒng)自動(dòng)回到開(kāi)始的狀態(tài)，即所有分子都在原點(diǎn)的狀態(tài)，已經(jīng)不可能了。 系統(tǒng)的熵和我們預(yù)測(cè)的一樣在隨時(shí)間增大，但當(dāng)所有分子布滿整個(gè)計(jì)算材料學(xué)概述之蒙特卡洛方法詳解計(jì)算材料學(xué)概述之蒙特卡洛方法詳解畫一個(gè)圓晶粒USE MSFLIB !INTEGER XR,YR !在的區(qū)域中畫一個(gè)圓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畫一個(gè)圓晶粒畫晶界!畫一個(gè)圓USE MSFLIB INTEGER XR，YR !在的區(qū)域中畫一個(gè)圓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如何畫有一定寬度的晶界？畫晶界!畫一個(gè)圓例題例 求前N個(gè)自然數(shù)倒數(shù)和。即求 程序需要幾個(gè)變量：累加的項(xiàng)數(shù)I(整型)，存放通項(xiàng)值的T(實(shí)型)，存放總和的S(實(shí)型)，待輸入的總項(xiàng)數(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 程序的運(yùn)行結(jié)果如下15 SUM= 3.182290 用DO結(jié)構(gòu)處理循環(huán)問(wèn)題，最關(guān)鍵的是要處理好初值、循環(huán)終止條件和循環(huán)體語(yǔ)句的關(guān)系。例題例 求前N個(gè)自然數(shù)倒數(shù)和。即求 程序需要幾個(gè)變量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)擴(kuò)散解 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)擴(kuò)散解 REAL C0(0:蒙特卡羅(Monte Carlo)方法1、簡(jiǎn)介2、相關(guān)幾個(gè)例子蒙特卡羅(Monte Carlo)方法1、簡(jiǎn)介Monte Carlo方法：亦稱統(tǒng)計(jì)模擬方法，statistical simulation method 利用隨機(jī)數(shù)進(jìn)行數(shù)值模擬的方法Mont

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

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

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

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

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

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

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

44、的收斂是以概率而言的誤差: 用頻率估計(jì)概率時(shí)誤差的估計(jì)，可由中心極限定理，給定置信水平 的條件下，有： 模擬次數(shù):由誤差公式得收斂性: 由大數(shù)定律, Monte-Carlo模擬的收斂是以隨機(jī)數(shù)的產(chǎn)生原理-計(jì)算機(jī)抽樣實(shí)驗(yàn) 隨機(jī)數(shù)的產(chǎn)生是實(shí)現(xiàn)MC計(jì)算的先決條件。而大多數(shù)概率分布的隨機(jī)數(shù)的產(chǎn)生都是基于均勻分布U(0,1的隨機(jī)數(shù)。 首先，介紹服從均勻分布U(0,1的隨機(jī)數(shù)的產(chǎn)生方法。其次，介紹服從其他各種分布的隨機(jī)數(shù)的產(chǎn)生方法。以及服從正態(tài)分布的隨機(jī)數(shù)的產(chǎn)生方法。最后，關(guān)于隨機(jī)數(shù)的幾點(diǎn)注。隨機(jī)數(shù)的產(chǎn)生原理-計(jì)算機(jī)抽樣實(shí)驗(yàn) 隨機(jī)數(shù)的產(chǎn)生是實(shí)現(xiàn)MC計(jì)(0,1均勻分布隨機(jī)數(shù)的產(chǎn)生最常用的是在(0,1區(qū)間內(nèi)均

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

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

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

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

49、產(chǎn)生總是采用某Monte Carlo方法相關(guān)引例:1、Buffon投針實(shí)驗(yàn)：18世紀(jì)，法國(guó)數(shù)學(xué)家Comte de Buffon利用投針實(shí)驗(yàn)估計(jì)的值dLMonte Carlo方法相關(guān)引例:1、Buffon投針實(shí)驗(yàn)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計(jì)算材料學(xué)概述之蒙特卡洛方法詳解計(jì)算材料學(xué)概述之蒙特卡洛方法詳解計(jì)算材料學(xué)概述之蒙特卡洛方法詳解圓周率的值 = 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什么是蒙特-卡洛模擬？計(jì)算機(jī)實(shí)驗(yàn)物理實(shí)驗(yàn) 用探測(cè)器取數(shù)據(jù)5RHIC-STAR實(shí)驗(yàn)探測(cè)器什么是蒙特-卡洛模擬？計(jì)算機(jī)實(shí)驗(yàn)物理實(shí)驗(yàn) 用探測(cè)器取數(shù) 蒙特-卡洛模擬 用計(jì)算機(jī)取數(shù)據(jù)計(jì)算材料學(xué)概述之蒙特卡洛方法詳解一個(gè)例子 道爾頓板理論 分布曲線實(shí)

55、驗(yàn) 丟鐵球模擬 計(jì)算機(jī)丟球6一個(gè)例子 道爾頓板理論 分布曲線6模擬的出發(fā)點(diǎn) 模型根據(jù)球和釘子碰撞的規(guī)律 追蹤每個(gè)球的運(yùn)動(dòng)根據(jù)落入每個(gè)格子的幾率 給出這一幾率的每次實(shí)現(xiàn)8模擬的出發(fā)點(diǎn) 模型根據(jù)球和釘子碰撞的規(guī)律89令 n 為鐵球數(shù)，m 為格子數(shù)。p(i) 是鐵球落入第 i 個(gè)格子中的概率再產(chǎn)生 n 個(gè)均勻分布的隨機(jī)數(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 個(gè)盒子中的粒子數(shù)。r(m) p(i)=1r(3) p(1)+p(2)r(2) p(1)r(1) 0i=1,m通過(guò)比較第 j 個(gè)鐵球的 ra(j) 和 r(i) 來(lái)決定這個(gè)球落入那一格子中：例一：道爾頓板的蒙特-卡洛模擬p3p4p5p2p1p6p79令 n 為鐵球數(shù)，m 為格子數(shù)。再產(chǎn)生 n 個(gè)均勻分布的隨例二 面積的計(jì)算f (x)x辛普遜方法I = Sn蒙特-卡洛方法f (x)x在長(zhǎng)方形中均勻投N0組(x,y)如 yf(x), 則 N=N+1I =(N/N0)S0SS011例二 面積的計(jì)算f (x)x辛普遜方法I = Sn蒙特-考慮平面上的一個(gè)邊長(zhǎng)為1的正方形及其內(nèi)部的一個(gè)形狀不規(guī)則的“圖形”，如何求出這個(gè)“圖形”的面積呢？Monte Carlo方法是這樣一種“隨機(jī)化”的方法：向該正方形“隨機(jī)地”投擲N個(gè)點(diǎn)，若有M個(gè)點(diǎn)落于“圖形”內(nèi)，則該“圖形”的面積近似為M/N?？紤]平面上的一個(gè)邊長(zhǎng)為1的正方形及其內(nèi)部的一個(gè)形狀不規(guī)則的“用該方法計(jì)算的基本思路是：1根據(jù)圓面積的公式：s=R2,當(dāng)R=1時(shí)，S=。由于圓的方程是：x2+y2=1(x2為x的平方的意思),