《計算物理學(xué)》課件 第五章 蒙特卡羅模擬方法

第五章蒙特卡羅模擬方法2確定性方法和隨機性方法確定性方法：分子動力學(xué)模擬--明確的相互作用方程模擬整個系統(tǒng)的行為--具有確定性，雖然與粒子運動的隨機過程有相互關(guān)聯(lián)微觀粒子的運動：是否具有確定性呢？本身是個隨機過程，具有概率的特征--例如，一個O2分子從教室一角擴散到教室的另外一角--不具有確定的描述需要一種隨機性方法來描述上述過程蒙特卡羅方法（MonteCarlo，MC）3蒙特卡羅方法1、蒙特卡羅方法概述2、隨機數(shù)3、對概率分布函數(shù)抽樣4、蒙特卡羅方法的應(yīng)用5、動力學(xué)蒙特卡羅方法1、蒙特卡羅方法概述蒙特卡洛方法的基本概念蒙特卡洛方法的數(shù)學(xué)基礎(chǔ)蒙特卡洛方法的主要構(gòu)成蒙特卡洛方法的優(yōu)缺點45蒙特卡羅方法MonteCarlomethods，或稱蒙特卡羅實驗MonteCarloexperiments，是一大類計算算法的集合，依靠重復(fù)的隨機抽樣來獲得數(shù)值結(jié)果?；靖拍钍抢秒S機性來解決理論上可能是確定性的問題。這類方法通常用于解決物理和數(shù)學(xué)問題，當面對棘手問題而束手無策時，往往它們可以大顯身手。蒙特卡羅方法主要用于解決3類問題：最優(yōu)化，數(shù)值積分，依據(jù)概率分布生成圖像。蒙特卡洛方法的基本概念6在物理學(xué)相關(guān)問題中，蒙特卡羅方法可用于模擬具有多個耦合自由度的系統(tǒng)，如流體、無序材料、強耦合固體和細胞結(jié)構(gòu)（參見細胞波茨模型CellularPottsModel、相互作用粒子系統(tǒng)InteractingParticleSystems、麥肯-弗拉索夫過程McKean-VlasovProcesses、氣體動力學(xué)模型KineticModelsofGases）。對輸入中具有重大不確定性的現(xiàn)象進行建模，如商業(yè)中的風險計算，以及在數(shù)學(xué)中對具有復(fù)雜邊界條件的多維定積分進行評估。在系統(tǒng)工程問題（空間、石油勘探、飛機設(shè)計等）的應(yīng)用中，基于蒙特卡羅的故障預(yù)測、成本超支和進度超支通常比人類的直覺或其他的“軟性”方法更有效。7理論上，蒙特卡羅方法可以用來解決任何具有概率解釋的問題。根據(jù)大數(shù)定律LawofLargeNumbers，用某個隨機變量的期望值描述的積分可以用該變量獨立樣本的經(jīng)驗均值（即樣本均值）來近似。當變量的概率分布被參數(shù)化時，數(shù)學(xué)家們經(jīng)常使用馬爾可夫鏈蒙特卡羅MarkovchainMonteCarlo(MCMC)采樣器。其中心思想是設(shè)計一個具有給定穩(wěn)態(tài)概率分布StationaryProbabilityDistribution的有效馬爾可夫鏈模型。也就是說，在極限情況下，馬爾科夫鏈蒙特卡洛方法生成的樣本將成為來自期望（目標）分布的樣本。通過遍歷定理ErgodicTheorem，穩(wěn)態(tài)分布可以用馬爾科夫鏈蒙特卡洛采樣器隨機狀態(tài)的經(jīng)驗測度來近似。

8蒙特卡羅方法各不相同，但趨于遵循一個特定的模式：1、定義可能輸入的域2、從域上的概率分布隨機生成輸入3、對輸入進行確定性計算4、匯總結(jié)果9例如，考慮一個單位正方形內(nèi)嵌的四分之一圓?？紤]到它們的面積比是π/4，π的值可以用蒙特卡羅方法來近似：1、畫一個正方形，然后在其中劃出一個四分之一圓2、在正方形上均勻散布給定數(shù)量的點3、計算四分之一圓內(nèi)的點數(shù)，即滿足距離原點小于1的4、四分之一圓內(nèi)部計數(shù)與總樣本計數(shù)之比是兩個區(qū)域之比的估計值，π/4。把結(jié)果乘以4就可以估算出π的值。在這個過程中，輸入域是限定四分之一圓的正方形。我們通過將顆粒散射到正方形上來產(chǎn)生隨機輸入，然后對每個輸入執(zhí)行計算（測試它是否在四分之一圓內(nèi)）。匯總這些結(jié)果會產(chǎn)生最終的結(jié)果——π的近似值。10兩個重要的考慮因素：1、如果這些點不是均勻分布的，那么近似效果就會很差。2、這一過程需要很多點。如果整個正方形中只有幾個點是隨機放置的，那么這個近似值通常是很差的。平均而言，隨著放置更多的點，近似值精度會提升。應(yīng)用蒙特卡羅方法需要大量的隨機數(shù)，這也就刺激了偽隨機數(shù)生成器的發(fā)展11蒙特卡羅方法的早期變種被設(shè)計來解決

布豐投針Buffon'sNeedleProblem問題，在布豐投針問題中，π可以通過將針落在由平行等距條組成的地板上來估計。20世紀30年代，恩里科·費米EnricoFermi在研究中子擴散時首次嘗試了蒙特卡羅方法，但他沒有發(fā)表這項工作。20世紀40年代末，斯坦尼斯拉夫·烏拉姆StanislawUlam在洛斯阿拉莫斯國家實驗室研究核武器項目時，發(fā)明了現(xiàn)代版的馬爾可夫鏈蒙特卡羅方法在烏拉姆的突破之后，約翰·馮·諾伊曼JohnvonNeumann立即意識到了它的重要性。馮·諾伊曼為ENIAC（人類第一臺電子數(shù)字積分計算機）編寫了程序來進行蒙特卡羅計算。蒙特卡羅方法簡史12他后來回憶當初靈感產(chǎn)生過程：

我最初構(gòu)想和嘗試蒙特卡洛法是在1946年，當時我正從疾病中康復(fù)，時常玩單人紙牌游戲。那時我會思考這樣一個問題：一盤52張的加菲爾德紙牌成功出牌的幾率有多大？在花了大量時間嘗試通過純粹的組合計算來估計它們之后，我想知道是否有一種比“抽象思維”更實際的方法，可能不是將它展開100次，然后簡單地觀察和計算成功的游戲數(shù)量。在快速計算機新時代開始時，這已經(jīng)是可以想象的了，我立刻想到了中子擴散和其他數(shù)學(xué)物理的問題，以及更一般的情形—如何將由某些微分方程描述的過程轉(zhuǎn)換成可解釋為一系列隨機操作的等價形式。后來（1946年），我向約翰·馮·諾伊曼描述了這個想法，然后我們開始計劃實際的計算。1946年，洛斯阿拉莫斯的核武器物理學(xué)家正在研究中子在可裂變材料中的擴散。盡管擁有大部分必要的數(shù)據(jù)，例如中子在與原子核碰撞之前在物質(zhì)中的平均運行距離，以及碰撞后中子可能釋放出多少能量，但洛斯阿拉莫斯的物理學(xué)家們無法用傳統(tǒng)的、確定性的數(shù)學(xué)方法解決這個問題。此時烏拉姆建議使用隨機實驗。13馮·諾依曼和烏拉姆的工作是秘密進行的，需要一個代號。馮·諾依曼和烏拉姆的一位同事，尼古拉斯·梅特羅波利斯NicholasMetropolis建議使用蒙特卡羅這個名字，這個名字指的是摩納哥的蒙特卡羅賭場，烏拉姆的叔叔會和親戚借錢然后去那里賭博。使用“真正隨機”的隨機數(shù)列表是非常慢的，然而馮·諾依曼使用平方取中法Middle-SquareMethod開發(fā)了一種計算偽隨機數(shù)生成器的方法。雖然許多人一直批評這種方法較為粗糙原始，但是馮·諾依曼也意識到這一點：他證明這種方法比任何其他方法都快，并指出當它出錯時，人們可以輕易發(fā)現(xiàn)，不像其他方法產(chǎn)生的錯誤可能會不易察覺。盡管受到當時的計算工具嚴重限制，蒙特卡羅方法依然是曼哈頓計劃ManhattanProject所需模擬的核心關(guān)鍵。20世紀50年代，它們在洛斯阿拉莫斯用于與氫彈開發(fā)有關(guān)的早期工作，并在物理學(xué)、物理化學(xué)和運籌學(xué)領(lǐng)域得到普及。14更復(fù)雜的平均場型粒子蒙特卡羅方法的理論產(chǎn)生于20世紀60年代中期，最初來自于小亨利·麥基恩HenryP.McKeanJr.研究流體力學(xué)中出現(xiàn)的一類非線性拋物型偏微分方程的馬爾可夫解釋。西奧多·愛德華·哈里斯TheodoreE.Harris和赫曼·卡恩HermanKahn在1951年發(fā)表了一篇開創(chuàng)性文章，使用平均場遺傳型蒙特卡羅方法來估計粒子傳輸能量。這一方法在演化計算中也被用作啟發(fā)式自然搜索算法（又稱元啟發(fā)式）。這些平均場計算技術(shù)的起源可以追溯到1950年和1954年，當時阿蘭·圖靈AlanTuring在基因類型突變-選擇學(xué)習機器上的工作，以及來自新澤西州普林斯頓高等研究院的尼爾斯·阿爾·巴里里利NilsAallBarricelli的文章。15量子蒙特卡羅方法，更具體地說，擴散蒙特卡羅方法也可以解釋為費曼-卡茨路徑積分Feynman–KacPathIntegrals的平均場粒子蒙特卡羅近似。量子蒙特卡羅方法的起源通常歸功于恩里科·費米EnricoFermi和羅伯特·里希特邁耶RobertRichtmyer于1948年開發(fā)了中子鏈式反應(yīng)的平均場粒子解釋，但是用于估計量子系統(tǒng)的基態(tài)能量（在簡化矩陣模型中）的第一個類啟發(fā)式和遺傳型粒子算法（也稱為重取樣或重構(gòu)蒙特卡洛方法）則是由杰克·H·海瑟林頓在1984年提出。16Buffon投針實驗1768年，法國數(shù)學(xué)家ComtedeBuffon利用投針實驗估計的值dL蒙特卡羅方法簡史17ProblemofBuffon’sneedle:

蒙特卡羅方法簡史18Solution:Thepositioningoftheneedlerelativetonearbylinescanbedescribedwitharandomvectorwhichhascomponents:Therandomvectorisuniformlydistributedontheregion[0,d)×[0,).Accordingly,ithasprobabilitydensityfunction1/d.Theprobabilitythattheneedlewillcrossoneofthelinesisgivenbytheintegral(1)蒙特卡羅方法簡史19

20EnricoFermi1930年，EnricoFermi利用MonteCarlo方法研究中子的擴散，并設(shè)計了一個MonteCarlo機械裝置，用于計算核反應(yīng)堆的臨界狀態(tài)。蒙特卡羅方法簡史21StanislawUlam,aPolishbornmathematicianwhoworkedforJohnvonNeumannontheUnitedStates’ManhattanProjectduringWorldWarII.UlamisprimarilyknownfordesigningthehydrogenbombwithEdwardTellerin1951.StanislawUlam

(1909-1984)HeinventedtheMonteCarlomethodin1946whileponderingtheprobabilitiesofwinningacardgameofsolitaire.

蒙特卡羅模擬的先驅(qū)22VonNeumann是MonteCarlo方法的正式奠基者,他與StanislawUlam合作建立了概率密度函數(shù)、反累積分布函數(shù)的數(shù)學(xué)基礎(chǔ)，以及偽隨機數(shù)產(chǎn)生器。在這些工作中，StanislawUlam意識到了數(shù)字計算機的重要性.

合作起源于Manhattan工程：利用ENIAC(ElectronicNumericalIntegratorandComputer)計算產(chǎn)額VonNeumann蒙特卡羅方法簡史23ThealgorithmbyMetropolis(andARosenbluth,MRosenbluth,ATellerandETeller,1953)hasbeencitedasamongthetop10algorithmshavingthe"greatestinfluenceonthedevelopmentandpracticeofscienceandengineeringinthe20thcentury."NicholasMetropolis(1915-1999)24方法名字的由來Monte-Carlo,MonacoMetropoliscoinedthename“MonteCarlo”,fromitsgamblingCasino.蒙特卡羅方法的基本思想基本思想：針對待求問題，根據(jù)物理現(xiàn)象本身的統(tǒng)計規(guī)律，或人為構(gòu)造一合適的依賴隨機變量的概率模型，使某些隨機變量的統(tǒng)計量為待求問題的解，進行大統(tǒng)計量N→∞的統(tǒng)計實驗方法或計算機隨機模擬方法。理論依據(jù)：大數(shù)定理：均勻分布的算術(shù)平均收斂于真值中心極限定理：置信水平下的統(tǒng)計誤差一個例子：

Buffen投針實驗求25實驗結(jié)果根據(jù)上式，一些人進行了實驗，其結(jié)果列于下表：實驗者年份投針次數(shù)π的實驗值Wolf185050003.1596Smith185532043.1553Fox189411203.1419Lazzarini190134083.1415929距今越近的實驗結(jié)果越準確幾千次的抽樣已經(jīng)可以得到很好的結(jié)果了測量方法、實驗設(shè)計等對最終實驗結(jié)果有較大的影響。2627薩維羅斯基Sawilowsky區(qū)分了模擬、蒙特卡羅方法和蒙特卡羅模擬：模擬是對現(xiàn)實的一種虛構(gòu)的表現(xiàn)，蒙特卡羅方法是一種可以用來解決數(shù)學(xué)或統(tǒng)計問題的技術(shù)，蒙特卡羅模擬使用重復(fù)抽樣來獲得某些現(xiàn)象（或行為）的統(tǒng)計特性。例如：模擬：從區(qū)間[0,1]中繪制一個偽隨機均勻變量可以用來模擬拋硬幣：如果值小于或等于0.50，則結(jié)果為正面，但如果值大于0.50，則結(jié)果為反面。這是一個模擬，但不是蒙特卡洛模擬。蒙特卡洛方法：將一盒硬幣倒在桌子上，然后計算正面與反面落地的硬幣比例。這是一種確定重復(fù)擲硬幣行為的蒙特卡洛方法，但它不是模擬。蒙特卡羅模擬法：一次或多次從區(qū)間[0,1]中繪制“大量”偽隨機均勻變量，賦值小于或等于0.50作為正面，大于0.50為反面。這是一個多次擲硬幣的“蒙特卡羅模擬”行為。28蒙特卡羅和隨機數(shù)MonteCarloandrandomnumbers這種方法的主要思想是基于重復(fù)隨機抽樣和統(tǒng)計分析來計算結(jié)果。蒙特卡洛模擬實際上是一種隨機實驗，在這種情況下，這些實驗的結(jié)果并不為人所知。蒙特卡羅模擬的典型特征是有許多未知參數(shù)，其中許多參數(shù)很難通過實驗獲得。蒙特卡羅模擬方法并不總是要求真正的隨機數(shù)是有用的（盡管對于一些應(yīng)用程序，如質(zhì)數(shù)測試，不可預(yù)測性是至關(guān)重要的）。許多最有用的技術(shù)是使用確定性的偽隨機序列，使測試和重新運行模擬變得很容易。偽隨機序列在某種意義上表現(xiàn)地“足夠隨機”，這是進行良好模擬唯一必需的性質(zhì)。29薩維羅斯基列出了高質(zhì)量蒙特卡羅模擬的特點：（偽隨機數(shù)）生成器具有某些特征（例如，序列重復(fù)之前有一個很長的“周期”）（偽隨機數(shù)）生成器可以產(chǎn)生能通過隨機性測試的值有足夠的樣本來確保準確的結(jié)果使用適當?shù)娜臃椒ㄊ褂玫乃惴▽?nèi)容是有效的它可以對問題中的現(xiàn)象進行模擬30蒙特卡羅方法的數(shù)學(xué)基礎(chǔ)大數(shù)定理：均勻分布的算術(shù)平均值收斂于真值。

31

32

減小方差的各種技巧

3334蒙特卡羅模擬的基本步驟根據(jù)問題本身，構(gòu)造或者確定一個概率模型對于本身不是隨機性質(zhì)的確定性問題，必須根據(jù)問題本身的特點和規(guī)律，建立一個概率模型，把確定性問題轉(zhuǎn)化成具有概率特征的問題，是所求的解就是這個模型的概率分布或數(shù)學(xué)期望。定義一個合適的隨機變量確定概率模型后，要定義一個對應(yīng)于這一物理問題的隨機變量，使這個概率模型下的隨機變量的概率分布或數(shù)學(xué)期望就是所模擬問題的解選擇隨機抽樣方法，通過模擬獲得子樣根據(jù)模型，選取隨機抽樣方法，實現(xiàn)由已知概率分布的隨機抽樣，得到具有給定分布的隨機數(shù)。用給定分布的隨機數(shù)進行模擬，相當于制造符合所定概率模型的一系列實驗數(shù)據(jù)，稱為“子樣”，具有原概率母體的性質(zhì)。統(tǒng)計計算用上述抽樣的算術(shù)平均值替代其數(shù)學(xué)期望（統(tǒng)計平均值），給出近似解35蒙特卡羅算法的主要組成部分概率密度函數(shù)(pdf)隨機數(shù)產(chǎn)生器抽樣規(guī)則模擬結(jié)果記錄

記錄一些感興趣的量的模擬結(jié)果

如何從在區(qū)間[0,1]上均勻分布的隨機數(shù)出發(fā),隨機抽取服從給定的pdf的隨機變量;

能夠產(chǎn)生在區(qū)間[0,1]上均勻分布的隨機數(shù)

必須給出描述一個物理系統(tǒng)的一組概率密度函數(shù);誤差估計減少方差的技術(shù)

利用該技術(shù)可減少模擬過程中計算的次數(shù)；

必須確定統(tǒng)計誤差（或方差）隨模擬次數(shù)以及其它一些量的變化；蒙特卡羅方法的優(yōu)缺點36蒙特卡羅方法能夠比較逼真地描述具有隨機性質(zhì)的事物的特點及物理實驗過程。受幾何條件限制小，收斂速度與問題的維數(shù)無關(guān)。具有同時計算多個方案與多個未知量的能力。誤差容易確定。程序結(jié)構(gòu)簡單，易于實現(xiàn)。

描述具有隨機性質(zhì)的事物能夠比較逼真地描述具有隨機性質(zhì)的事物的特點及物理實驗過程從這個意義上講，蒙特卡羅方法可以部分代替物理實驗，甚至可以得到物理實驗難以得到的結(jié)果。用蒙特卡羅方法解決實際問題，可以直接從實際問題本身出發(fā)，而不從方程或數(shù)學(xué)表達式出發(fā)。它有直觀、形象的特點。37受幾何條件限制小在計算s維空間中的任一區(qū)域Ds上的積分時，無論區(qū)域Ds的形狀多么特殊，只要能給出描述Ds的幾何特征的條件，就可以從Ds中均勻產(chǎn)生N個點，得到積分的近似值。其中Ds為區(qū)域Ds的體積。這是數(shù)值方法難以作到的。另外，在具有隨機性質(zhì)的問題中，如考慮的系統(tǒng)形狀很復(fù)雜，難以用一般數(shù)值方法求解，而使用蒙特卡羅方法，不會有原則上的困難。38收斂速度和問題的維數(shù)無關(guān)由誤差定義可知，在給定置信水平情況下，蒙特卡羅方法的收斂速度為，與問題本身的維數(shù)無關(guān)。維數(shù)的變化，只引起抽樣時間及估計量計算時間的變化，不影響誤差。也就是說，使用蒙特卡羅方法時，抽取的子樣總數(shù)N與維數(shù)s無關(guān)。維數(shù)的增加，除了增加相應(yīng)的計算量外，不影響問題的誤差。這一特點，決定了蒙特卡羅方法對多維問題的適應(yīng)性。而一般數(shù)值方法，比如計算定積分時，計算時間隨維數(shù)的冪次方而增加，而且，由于分點數(shù)與維數(shù)的冪次方成正比，需占用相當數(shù)量的計算機內(nèi)存，這些都是一般數(shù)值方法計算高維積分時難以克服的問題。39同時計算多個方案與多個未知量對于那些需要計算多個方案的問題，使用蒙特卡羅方法有時不需要像常規(guī)方法那樣逐個計算，而可以同時計算所有的方案，其全部計算量幾乎與計算一個方案的計算量相當。例如，對于屏蔽層為均勻介質(zhì)的平板幾何，要計算若干種厚度的穿透概率時，只需計算最厚的一種情況，其他厚度的穿透概率在計算最厚一種情況時稍加處理便可同時得到。另外，使用蒙特卡羅方法還可以同時得到若干個所求量。例如，在模擬粒子過程中，可以同時得到不同區(qū)域的通量、能譜、角分布等，而不像常規(guī)方法那樣，需要逐一計算所求量。40誤差容易確定對于一般計算方法，要給出計算結(jié)果與真值的誤差并不是一件容易的事情，而蒙特卡羅方法則不然。根據(jù)蒙特卡羅方法的誤差公式，可以在計算所求量的同時計算出誤差。對干很復(fù)雜的蒙特卡羅方法計算問題，也是容易確定的。一般計算方法常存在著有效位數(shù)損失問題，而要解決這一問題有時相當困難，蒙特卡羅方法則不存在這一問題。41不足之處：收斂速度慢如前所述，蒙特卡羅方法的收斂速度為，一般不容易得到精確度較高的近似結(jié)果。對于維數(shù)少（三維以下）的問題，不如其他方法好。誤差具有概率性由于蒙特卡羅方法的誤差是在一定置信水平下估計的，所以它的誤差具有概率性，而不是一般意義下的誤差。42不足：計算結(jié)果與系統(tǒng)大小有關(guān)例如在粒子輸運問題中:經(jīng)驗表明，只有當系統(tǒng)的大小與粒子的平均自由程可以相比較時（一般在十個平均自由程左右），蒙特卡羅方法計算的結(jié)果較為滿意。但對于大系統(tǒng)或小概率事件的計算問題，計算結(jié)果往往比真值偏低。因此，在使用蒙特卡羅方法時，可以考慮把蒙特卡羅方法與解析（或數(shù)值）方法相結(jié)合，取長補短。既能解決解析（或數(shù)值）方法難以解決的問題，也可以解決單純使用蒙特卡羅方法難以解決的問題。這樣，可以發(fā)揮蒙特卡羅方法的特長，使其應(yīng)用范圍更加廣泛。43主要應(yīng)用范圍根據(jù)前面分析的特點，蒙特卡羅方法特別適用于處理維數(shù)高、邊界復(fù)雜而精度要求不是很高的問題。大量的科學(xué)、工程問題均具有上述特點，因此其有著廣泛的應(yīng)用。它的主要應(yīng)用范圍包括：典型數(shù)學(xué)問題，粒子輸運問題，統(tǒng)計物理，量子力學(xué)，真空技術(shù)，激光技術(shù)以及醫(yī)學(xué)，生物，探礦等方面。此外，其在金融、經(jīng)濟、社會學(xué)等領(lǐng)域也有著廣泛的應(yīng)用。隨著科學(xué)技術(shù)的發(fā)展，其應(yīng)用范圍將更加廣泛。442、隨機數(shù)與偽隨機數(shù)隨機數(shù)的定義及產(chǎn)生方法偽隨機數(shù)產(chǎn)生偽隨機數(shù)的各種方法4546用蒙特卡羅方法在計算機上模擬一個隨機過程，其首要任務(wù)是產(chǎn)生滿足相應(yīng)概率分布的隨機變量。在連續(xù)型隨機變量的分布中，最簡單、最基本的分布是單位均勻分布。由該分布抽取的簡單子樣稱為隨機數(shù)序列，其中每一個體稱為隨機數(shù)，且在[0,1]區(qū)間上均勻分布。因此，隨機數(shù)是指一個數(shù)列，其中的每一個體稱為隨機數(shù)，其值與數(shù)列中的其它數(shù)無關(guān)。在一個均勻分布的隨機數(shù)中，每一個體出現(xiàn)的概率是均等的。將[0,1]區(qū)間上均勻分布的隨機數(shù)作為已知量，用適當?shù)臄?shù)學(xué)方法可以產(chǎn)生具有任意已知分布的簡單子樣。

蒙特卡羅方法與隨機數(shù)的關(guān)系47隨機數(shù)的定義和特性什么是隨機數(shù)？單個的數(shù)字不是隨機數(shù)是指一個數(shù)列，其中的每一個體稱為隨機數(shù)，其值與數(shù)列中的其它數(shù)無關(guān)；在一個均勻分布的隨機數(shù)中，每一個體出現(xiàn)的概率是均等的；例如：在[0,1]區(qū)間上均勻分布的隨機數(shù)序列中，0.00001與0.5出現(xiàn)的機會均等隨機數(shù)的定義及性質(zhì)在連續(xù)型隨機變量的分布中，最簡單而且最基本的分布是單位均勻分布。由該分布抽取的簡單子樣稱為隨機數(shù)序列，其中每一個體稱為隨機數(shù)。單位均勻分布也稱為［0，1］上的均勻分布，其分布密度函數(shù)為：分布函數(shù)為：4849隨機數(shù)應(yīng)具有的基本特性考慮一個對高能粒子反應(yīng)過程的模擬：需用隨機數(shù)確定：出射粒子的屬性：能量、方向、…粒子與介質(zhì)的相互作用對這一過程的模擬應(yīng)滿足以下要求（相空間產(chǎn)生過程）：出射粒子的屬性應(yīng)是互不相關(guān)的，即每一粒子的屬性的確定獨立于其它的粒子的屬性的確定；粒子的屬性的分布應(yīng)滿足物理所要求的理論分布；所模擬的物理過程要求隨機數(shù)應(yīng)具有下列特性：隨機數(shù)序列應(yīng)是獨立的、互不相關(guān)的(uncorrelated)：即序列中的任一子序列應(yīng)與其它的子序列無關(guān)；50均勻分布的隨機數(shù)應(yīng)滿足均勻性(Uniformity)：隨機數(shù)序列應(yīng)是均勻的、無偏的，即：如果兩個子區(qū)間的“面積”相等，則落于這兩個子區(qū)間內(nèi)的隨機數(shù)的個數(shù)應(yīng)相等。例如：對[0,1)區(qū)間均勻分布的隨機數(shù)，如果產(chǎn)生了足夠多的隨機數(shù)，而有一半的隨機數(shù)落于區(qū)間[0,0.1]

不滿足均勻性如果均勻性不滿足，則會出現(xiàn)序列中的多組隨機數(shù)相關(guān)的情況

均勻性與互不相關(guān)的特性是有聯(lián)系的所模擬的物理過程要求隨機數(shù)應(yīng)具有下列特性：隨機數(shù)序列應(yīng)是獨立的、互不相關(guān)的(uncorrelated)：即序列中的任一子序列應(yīng)與其它的子序列無關(guān)；51隨機數(shù)的產(chǎn)生[0,1]區(qū)間上均勻分布的隨機數(shù)是蒙特卡洛模擬的基礎(chǔ)服從任意分布的隨機數(shù)序列可以用[0,1]區(qū)間均勻分布的隨機數(shù)序列作適當?shù)淖儞Q或舍選后求得。因此，下面著重討論[0,1]均勻分布的隨機數(shù)的產(chǎn)生方法。

由于隨機數(shù)在蒙特卡羅方法中所處的特殊地位，它們雖然也屬于由具有已知分布的總體中產(chǎn)生簡單子樣的問題，但就產(chǎn)生方法而言，卻有著本質(zhì)上的差別。為了產(chǎn)生隨機數(shù)，可以使用隨機數(shù)表。隨機數(shù)表是由0，1，…，9十個數(shù)字組成，每個數(shù)字以0.1的等概率出現(xiàn)，數(shù)字之間相互獨立。這些數(shù)字序列叫作隨機數(shù)字序列。如果要得到n位有效數(shù)字的隨機數(shù)，只需將表中每n個相鄰的隨機數(shù)字合并在一起，且在最高位的前邊加上小數(shù)點即可。例如，某隨機數(shù)表的第一行數(shù)字為86397584102…，要想得到三位有效數(shù)字的隨機數(shù)依次為0.863，0.975，0.841。因為隨機數(shù)表需在計算機中占有很大內(nèi)存，而且也難以滿足蒙特卡羅方法對隨機數(shù)需要量非常大的要求，因此，該方法不適于在計算機上使用。隨機數(shù)的產(chǎn)生方法——隨機數(shù)表52隨機數(shù)的產(chǎn)生方法——物理方法（1）用物理方法產(chǎn)生隨機數(shù)的基本原理是：利用某些物理現(xiàn)象，在計算機上增加些特殊設(shè)備，可以在計算機上直接產(chǎn)生隨機數(shù)。這些特殊設(shè)備稱為隨機數(shù)發(fā)生器。用來作為隨機數(shù)發(fā)生器的物理源主要有兩種：一種是根據(jù)放射性物質(zhì)的放射性，另一種是利用計算機的固有噪聲。一般情況下，任意一個隨機數(shù)在計算機內(nèi)總是用二進制的數(shù)表示的：其中εi（i=1，2，…，m）或者為0，或者為1。53隨機數(shù)的產(chǎn)生方法——物理方法（2）因此，利用物理方法在計算機上產(chǎn)生隨機數(shù)，就是要產(chǎn)生只取0或1的隨機數(shù)字序列，數(shù)字之間相互獨立，每個數(shù)字取0或1的概率均為0.5。用物理方法產(chǎn)生的隨機數(shù)序列無法重復(fù)實現(xiàn)，不能進行程序復(fù)算，給驗證結(jié)果帶來很大困難。而且，需要增加隨機數(shù)發(fā)生器和電路聯(lián)系等附加設(shè)備，費用昂貴。因此，該方法也不適合在計算機上使用。5455

隨機數(shù)的產(chǎn)生方法——數(shù)學(xué)方法存在的問題

56第一個問題，不能從本質(zhì)上加以改變，但只要遞推公式選得比較好，隨機數(shù)間的相互獨立性是可以近似滿足的。第二個問題，則不是本質(zhì)的。因為用蒙特卡羅方法解任何具體問題時，所使用的隨機數(shù)的個數(shù)總是有限的，只要所用隨機數(shù)的個數(shù)不超過偽隨機數(shù)序列出現(xiàn)循環(huán)現(xiàn)象時的長度就可以了。用數(shù)學(xué)方法產(chǎn)生的偽隨機數(shù)容易在計算機上得到，可以進行復(fù)算，而且不受計算機型號的限制。因此，這種方法雖然存在著一些問題，但仍然被廣泛地在計算機上使用，是在計算機上產(chǎn)生偽隨機數(shù)的主要方法?？梢越鉀Q否？57計算程序產(chǎn)生的隨機數(shù)并不是真正的隨機數(shù)，它們是確定的，但看上去是隨機的，且能通過一些隨機性的檢驗，故常稱為偽隨機數(shù)。發(fā)生周期性循環(huán)現(xiàn)象的偽隨機數(shù)的個數(shù)稱為偽隨機數(shù)的周期。從偽隨機數(shù)序列的初始值開始，到出現(xiàn)循環(huán)現(xiàn)象為止，所產(chǎn)生的偽隨機數(shù)的個數(shù)稱為偽隨機數(shù)的最大容量。隨機數(shù)發(fā)生器5859線性同余法（LinearCongruentialMethod)1948年由Lehmer提出的一種產(chǎn)生偽隨機數(shù)的方法,是最常用的方法。遞推公式：

該方法產(chǎn)生整型的隨機數(shù)序列，隨機性來源于取模運算。如果產(chǎn)生的是[0,1]之間的隨機數(shù)，則只需要每個數(shù)除以m即可

60

通常將m取為計算機所能表示的最大的整型量，在32位計算機上，m=231=2x109a,c,m的選擇：用線性同余方法產(chǎn)生的隨機數(shù)序列具有周期m的條件是：c和m為互質(zhì)數(shù)；a-1是m的任一奇數(shù)因子的倍數(shù)；如果m是4的倍數(shù)，a-1也是4的倍數(shù)。例：a=5,c=1,m=16,I0=1周期=m=161,6,15,12,13,2,11,8,9,14,7,4,5,10,3,0,1,6,15,12,13,2,..線性同余方法在計算機上的使用為了便于在計算機上使用，通常取：

m=2s其中s為計算機中二進制數(shù)的最大可能有效位數(shù)

x1=奇數(shù)

a=52k+1其中k為使52k+1在計算機上所能容納的最大整數(shù)，即a為計算機上所能容納的5的最大奇次冪。一般地，s=32時，a=513；s=48，a=515等。偽隨機數(shù)序列的最大容量λ(m)=2s-2。乘同余方法是使用的最多、最廣的方法，在計算機上被廣泛地使用。6162乘同余法（multiplicativeLinearCongruentialMethod)

63從概率分布函數(shù)的抽樣

(SamplingfromProbabilityDistributionFunctions)64對概率分布函數(shù)的抽樣MonteCarlo算法的一個重要組成部分：描述所要模擬的物理系統(tǒng)的一些概率密度函數(shù)（PDF)MonteCarlo模擬的主要任務(wù):通過對這些概率密度函數(shù)的隨機抽樣來模擬物理系統(tǒng)的狀態(tài);為描述系統(tǒng)的演化所必需的一些附加運算.物理過程的描述從描述物理系統(tǒng)的pdf出發(fā)，隨機抽取系統(tǒng)的可能狀態(tài)。

描述整個系統(tǒng)在空間、能量、時間或多維相空間中的發(fā)展和演化；65直接抽樣法（反函數(shù)法）2.舍選抽樣法3.變換抽樣方法4.重要抽樣法66注意：pdf

f(x)必須是歸一化的設(shè)y=F(x)為隨機變量x的累積分布函數(shù)

x和y是一一對應(yīng)的先隨機抽取y，然后通過求F(x)的反函數(shù)F-1(y)得到隨機變量x的值隨機變量y在區(qū)間[0,1]上均勻分布利用[0,1]區(qū)間上均勻分布隨機數(shù)產(chǎn)生器抽取直接抽樣方法是根據(jù)概率分布進行采樣。對一個已知概率密度函數(shù)與累積概率密度函數(shù)的概率分布，我們可以直接從累積分布函數(shù)（cdf）進行采樣使用累積分布函數(shù)進行采樣看似簡單，但是由于很多分布我們并不能寫出概率密度函數(shù)與累積分布函數(shù)，所以這種方法的適用范圍較窄。直接抽樣法（反函數(shù)法）67p3=0.2b3+c3p2=0.3b2+c2p1=0.5b1+c1a

例、粒子衰變末態(tài)的隨機抽樣設(shè)粒子a有三種衰變方式，其分支比如下隨機選取每次衰變的衰變方式（衰變道）直接抽樣法對指數(shù)分布的直接抽樣分布密度函數(shù)為：積分得到分布函數(shù)：令則指數(shù)分布隨機變量的抽樣為：注意（1-ξ）和

ξ同樣服從[0，1]的均勻分布68舍選抽樣法（acceptance-rejectionsampling)直接抽樣法的困難：許多隨機變量的累積分布函數(shù)無法用解析函數(shù)給出;有些隨機變量的累積分布函數(shù)的反函數(shù)不存在或難以求出；即使反函數(shù)存在，但計算困難

簡單舍選抽樣法舍選法抽樣步驟：

設(shè)隨機變量x的取值區(qū)間為x

[a,b],其概率密度函數(shù)f(x)有界，即

y

f(x)X=x

>

abxf(x)c幾何解釋:在二維圖上，隨機選取位于矩形abef內(nèi)的點[x,y]；選取位于曲線f(x)下的那些點，則這些點將服從概率密度為f(x)的分布ef抽樣效率:對舍選抽樣法：欲產(chǎn)生m個隨機變量x的值需產(chǎn)生n對(x,y)，顯然，mn

abxf(x)cefd改進的舍選抽樣法

改進的舍選抽樣法簡單舍選抽樣法的問題：如果f(x)曲線下的面積占矩形面積的比例很小，則抽樣效率很低，這是因為隨機數(shù)x和y是在區(qū)間[a,b]和[0,c]內(nèi)均勻分布，所產(chǎn)生的大部分投點不會落在f(x)曲線下改進方法：構(gòu)造一個新的概率密度函數(shù)g(x)，使它的形狀接近f(x),且有式中C為常數(shù)，而g(x)的抽樣相對比較容易。這樣抽樣的效率會大大提高。xcf(x)Cgg(x)

改進的舍選抽樣法抽樣方法：1.產(chǎn)生兩個隨機數(shù)產(chǎn)生分布為g(x)

的隨機數(shù)x

，x[a,b];產(chǎn)生[0,Cg(x)]

區(qū)間上均勻分布的隨機數(shù)y，y=Cg

(x)

,[0,1].2.接收或舍棄取樣值

x.如果

y>f(x)，舍棄，返回到1,重復(fù)上述過程;否則，接受;改進的舍選抽樣法幾何解釋:在二維圖上，隨機選取位于曲線Cg(x)下的點[x,y]；選取位于曲線f(x)下的那些點，則這些點將服從概率密度為f(x)的分布xcf(x)Cg(x)例1:標準正態(tài)分布的抽樣,x[-a,a]無法用直接抽樣法,累積分布函數(shù)無解析表達式

77

78變換抽樣法變換抽樣法的基本思想是將一個比較復(fù)雜的分布函數(shù)的抽樣變換為一個已知的簡單分布函數(shù)的抽樣。

79

80

81

82

83

這兩個分布的抽樣可以采用直接抽樣獲得。

84

85舍選抽樣法完美的解決了累積分布函數(shù)不可求時的抽樣問題。但是舍選抽樣非常依賴于提議分布(proposaldistribution)的選擇。如果提議分布選擇的不好，可能抽樣時間很長卻獲得很少滿足分布的粒子。而重要性采樣就解決了這一問題Importance

sampling86重要抽樣法重要抽樣法的基本思想是通過改變隨機變量的權(quán)重，改變樣本空間的概率分布，使用加權(quán)平均的方法得到期望值。

87

88

89

MC模擬生日問題假設(shè)有n個人在一起，各自的生日為365天之一，根據(jù)概率理論，與很多人的直覺相反，只需23個人便有大于50％的幾率人群中至少有2個人生日相同。n理論幾率模擬幾率

0.1170.1100.4110.4120.5070.5200.7060.6920.9410.936500.9860.987中子屏蔽問題NeutronShieldingproblem假設(shè)鉛墻長為5，中子在鉛中的平均自由程為1，中子與鉛原子碰撞后各向同性散射。令碰撞8次后中子能量耗盡，試求穿透鉛墻的中子的比例。暫不考慮垂直紙面的運動，則中子的水平位移是。入口鉛墻（長為5）排隊系統(tǒng)：排隊系統(tǒng)模擬：M/G/1排隊系統(tǒng)服務(wù)規(guī)則：先到先服務(wù)（Firstcome,firstservice:FCFS）假設(shè)：（１）顧客到達遵循Poisson分布；（２）服務(wù)時間服從一般分布；（３）到達間隔與服務(wù)時間相互獨立．關(guān)心的指標：

（１）時刻t時，系統(tǒng)中的顧客數(shù)；即隊長的分布；（２）顧客的等待時間；（３）服務(wù)的忙碌程度；（４）．．．．．．．用最樸素的Monte-Carlo方法可以得到這些指標的估計．排隊論起源于20世紀初的電話通話。1909-1920丹麥數(shù)學(xué)家，電氣工程師A.K.Erlang用概率論方法研究電話通話問題。排隊論的應(yīng)用非常廣泛。適用于一切服務(wù)系統(tǒng)。尤其在通信系統(tǒng)，交通系統(tǒng)，計算機，存儲系統(tǒng)，生產(chǎn)管理系統(tǒng)等方面應(yīng)用得最多。連續(xù)系統(tǒng)模擬實例:追逐問題

狀態(tài)隨時間連續(xù)變化的系統(tǒng)稱為連續(xù)系統(tǒng)。對連續(xù)系統(tǒng)的計算機模擬只能是近似的，只要這種近似達到一定的精度，也就可以滿足要求。例追逐問題:

如圖,正方形ABCD的四個頂點各有一人.在某一時刻,四人同時出發(fā)以勻速v=1米/秒按順時針方向追逐下一人,如果他們始終保持對準目標,則最終按螺旋狀曲線于中心點O.試求出這種情況下每個人的行進軌跡.OBCDA1.建立平面直角坐標系:A(x1,y1),B(x2,y2),C(x3,y3),D(x4,y4).2.取時間間隔為Δt,計算每一點在各個時刻的坐標.4.對每一個點，連接它在各時刻的位置,即得所求運動軌跡.求解過程:設(shè)某點在t時刻的坐標為:則在時刻的坐標為:其中

v=1;dt=0.05;x=[001010];y=[010100];fori=1:4plot(x(i),y(i),'.'),holdonendd=20;while(d>0.1)x(5)=x(1);y(5)=y(1);

fori=1:4d=sqrt((x(i+1)-x(i))^2+(y(i+1)-y(i))^2);x(i)=x(i)+v*dt*(x(i+1)-x(i))/d;y(i)=y(i)+v*dt*(y(i+1)-y(i))/d;plot(x(i),y(i),'.'),holdon

endend計算程序:編一個福利彩票電腦選號的程序。

某設(shè)備上安裝有四只型號規(guī)格完全相同的電子管，已知電子管壽命為1000--2000小時之間的均勻分布。當電子管損壞時有兩種維修方案，一是每次更換損壞的那一只；二是當其中一只損壞時四只同時更換。已知更換時間為換一只時需1小時，4只同時換為2小時。更換時機器因停止運轉(zhuǎn)每小時的損失為20元，又每只電子管價格10元，試用模擬方法決定哪一個方案經(jīng)濟合理？

導(dǎo)彈追蹤問題：設(shè)位于坐標原點的甲艦向位于x軸上點A(1,0)處的乙艦發(fā)射導(dǎo)彈，導(dǎo)彈頭始終對準乙艦.如果乙艦以最大的速度(是常數(shù))沿平行于y軸的直線行駛，導(dǎo)彈的速度是5，模擬導(dǎo)彈運行的軌跡.又乙艦行駛多遠時，導(dǎo)彈將它擊中？98MC解決隨機性問題物理研究中經(jīng)常遇到一些隨機性問題，例如粒子輸運/擴散，晶體生長等等。MC方法可以對這些過程進行直接的模擬，因此這些問題本身就提供了一個概率模型（省去了抽象的困難）通過MC將一個隨機過程和多次過程之后的效果聯(lián)系起來了，具體應(yīng)用分為兩種：一是隨機過程的概率分別已知，可以直接使用MC對復(fù)雜過程的結(jié)果做出預(yù)言；二是概率未知，則需要假定一個模型，通過與實驗結(jié)果比較，判斷模型的正確性。兩種情況在解決實際問題中都有應(yīng)用99庫存管理問題某種商品進貨價格為a元，出售價格為b元，假設(shè)對該商品每天早晨進貨配齊n個(必須當天賣掉)，每日顧客相互獨立地到來，平均每日m人，且服從Poisson分布，

其中k為每日到來的人數(shù)。問當a=2、b＝3、m=10時，每天早晨該商品備齊多少個可得到最大利潤。假設(shè)購買此商品的顧客每人只購一個，而且備齊的商品如果當天賣不出去則不可在第二天出售，若用a(n,k)表示顧客為k人時的日利潤，則100程序?qū)崿F(xiàn)N=365;%mc次數(shù)n=15;%進貨，考慮最多15個m=10;%平均人流量a=2;b=3;forj=1:ntotal=0;%總利潤randnumber=poissrnd(m,N,1);%Poisson分布隨機數(shù)fori=1:N ifrandnumber(i)<j total=total+randnumber(i)*b-j*a; else total=total+j*b-j*a; endendy(j)=total/N%平均利潤end101計算結(jié)果12345678Try11.0001.9922.9673.9754.9265.5236.3676.907Try21.0001.9842.9753.9674.8855.6716.3346.7019101112131415Try16.8146.3925.8054.8162.8820.858-0.0417.0116.0305.5514.5872.1841.786-0.625每日進貨8或者9個是利潤量最大利潤分辨率還不高，可以考慮前期的波動則有選擇8或者9102自然界中的隨機過程103隨機游走問題隨機游走問題最早是Pearson在1905年提出的。假設(shè)有個醉漢從一根電線桿的位置出發(fā)（其坐標為x=0，x坐標向右為正，向左為負），假定醉漢的步長為l，他走的每一步的取向是隨機的，與前一步的方向無關(guān)。如果醉漢在每個時間間隔內(nèi)向右行走一步的幾率為p，則向左走一步的幾率為q=1-p。我們記錄醉漢向右走了nr

步，向左走了nl

步，即總共走了N=nr+nl

步。那末醉漢在行走了N步以后，離電線桿的距離為x=(nr-nl)l，其中-Nl=<x=<Nl。然而我們更感興趣的是醉漢在行走N步以后，離電線桿的距離為x的概率PN(x)104方差估計計算醉漢在走了N步后的位移和方差的平均值<xN>,<ΔxN2>其中公式中的求平均是指對步中所有可能的行走過程的平均。注意到在向左、向右對稱的情況下，即p=q=1/2，得到<xN>=0105求解上述問題：查點法在查點法中，對給定的行走總步數(shù)N及總位移x，要求把游走時可能的每一步的坐標和幾率都確定下來。這是可以用概率理論精確計算的。例:對于N=3,l=1的醉漢一維行走問題，由概率理論可得由此可以算出則:查點法只有在總步數(shù)N較小時才可以使用。N比較大時用起來就比較困難了。106蒙特卡羅方法求解蒙特卡洛方法就可以克服在游走中的這個困難，具有更廣泛的可操作性。蒙特卡洛方法可以對許多步的游走過程進行抽樣。例如。我們可以按照正確的概率，對確定的產(chǎn)生出各種可能的行走樣本。對于一維情況：假設(shè)p、q分別為向左、右行走的概率，則產(chǎn)生隨機數(shù)ξ若0=<ξ=<p，則x-1若p<ξ=<1，則x+1原則上只要我們增加抽樣的個數(shù)，要達到較高的精度總是可能的。107兩次獨立隨機行走的路徑步長為+-1500次的平均結(jié)果D:

diffusion

constant與自由粒子做對比：距離隨時間線性增長隨機行走模擬108步長不定的隨機行走兩次獨立隨機行走的路徑步長是隨機的，[-1,+1]500次的平均結(jié)果109三維隨機行走（500次的平均結(jié)果，步長為+-1）110Self-avoiding

walkA

random

walk

that

is

subject

to

the

constraint

of

avoiding

intersect

itself

is

called

a

self-avoiding

walk.

Each

distinct

SAW

of

the

same

number

of

steps

must

occur

with

equal

probability.

The

collection

of

all

such

SAWs

is

often

called

the

SAW

ensemble.Simulation:

We

must

keep

track

of

all

prior

steps

and

make

sure

that

configurations

that

would

revisit

a

previously

trampled

site

are

not

included.我們可以從一點開始，隨機選擇可走的近鄰位置，直到?jīng)]有可選的近鄰位置為止。這種算法有問題嗎？有什么問題？111想象一個在平衡溶液中的高分子鏈鏈的所有可能構(gòu)型恰恰與SAW的所有可能構(gòu)型一致，這些構(gòu)型應(yīng)該是等概率出現(xiàn)的。上面的算法不滿足等概率條件。采用上述算法產(chǎn)生的系綜是不滿足所有可能的SAW等概率出現(xiàn)怎么辦？112最簡單的方式：從所有可能的行走方向上隨機選擇下一步，如果遇到了self-intersection，就放棄整個行走過程，重新開始。這樣可以保證產(chǎn)生的SAW構(gòu)型出現(xiàn)的概率是相等的。

1132-dimensional

self-avoiding

walksn=3/4(2D),

3/5(3D)114Cream-in-coffee

problem400個粒子（20X20）限制在x=+-100,

y=+-100范圍之內(nèi)115Pi

is

the

probability

of

finding

a

particle

in

cell

i.1161953年，Metropolis等研究了由具有相互作用分子組成的物質(zhì)系統(tǒng)的性質(zhì)。在利用蒙特卡洛方法計算高維積分時，提出了一種新的相空間的采樣算法，后來被普遍稱為Metropolis算法。簡單來說，就是通過采用轉(zhuǎn)移概率從前一個態(tài)產(chǎn)生新的態(tài)，利用細致平衡條件得到目標概率分布。因此，Metropolis算法是一種根據(jù)玻爾茲曼分布生成系統(tǒng)狀態(tài)的馬爾科夫鏈蒙特卡洛方法。之后，Hastings改進了該算法，提高了采樣率，能夠模擬隨機變量序列，更精確地模擬了期望分布為平穩(wěn)分布的馬爾科夫鏈，特別是在許多隨機變量的分布無法直接模擬的情況下，稱為Metropolis-Hastings算法。該算法已經(jīng)被廣泛應(yīng)用于物理系統(tǒng)的蒙特卡洛模擬。Metropolis方法117

118

119

這樣，Metropolis算法通過構(gòu)造一個滿足細致平衡條件的馬爾科夫鏈的演化過程，產(chǎn)生一個狀態(tài)系列。

120隨機生長過程的模擬-DLA1981年,美國密捷安大學(xué)Witten和Sander開創(chuàng)性地提出了一個擴散限制凝聚(DiffusionLimitedAggregation)的分形生長模型,簡稱DLA模型。最初該模型提出時主要是為了研究懸浮在大氣中的煤灰、金屬粉末或煙塵擴散的凝聚問題。后來受到不同領(lǐng)域的學(xué)者重視,被引入不同的學(xué)科,用來解釋各種與分形形態(tài)有關(guān)的生長和凝聚現(xiàn)象.121DLA模型DLA模型采用了在晶格原點上置一個初始粒子作為種子,以該點為圓心(原點),r為半徑作一個大圓。在該圓附近隨機產(chǎn)生一個微粒在圓上作近似于Brown運動的隨機行走,即以1/4的概率向上下左右方向行走。

若微粒與種子顆粒相碰,則令其附著于種子微粒之上,與之結(jié)合形成凝聚集團;若微粒行走到圓的邊界或離開此圓,令其被邊界吸收而消失接著再隨機地產(chǎn)生第二個微粒并重復(fù)以上步驟,直至附著于凝聚集團或被邊界吸收。如此不斷地進行,當種子集團長大到一定程度后,可以發(fā)現(xiàn),它就形成了以種子微粒為中心的無規(guī)分叉圖形。122DLA模型的性質(zhì)和特點DLA模型的主要性質(zhì)和特點是:這類分形結(jié)構(gòu)在形成過程中都遵從可動邊界的拉普拉斯方程以凝聚生長為例,凝聚粒子在拉普拉斯?jié)舛葓鲋凶鳠o規(guī)擴散運動,生長集團的界面具有復(fù)雜的性質(zhì)和不穩(wěn)定的性質(zhì),生長過程是一個遠離平衡的動力學(xué)過程,是一個非線性的、非平衡態(tài)的進化過程可以說,DLA模型的出現(xiàn),深化了人們對非平衡生長現(xiàn)象的認識,吸引了大批數(shù)學(xué)、物理、化學(xué)、生物科學(xué)家投入其中,研究有關(guān)物理的、化學(xué)的、醫(yī)學(xué)的生長現(xiàn)象,并已經(jīng)取得了令人鼓舞的進展123計算機實現(xiàn)構(gòu)造二維格子初始化種子外圍隨機產(chǎn)生粒子粒子進行隨機游走如果碰到種子則（以一定幾率）吸附吸附后重復(fù)產(chǎn)生粒子吸附幾率的控制，隨機游走是否幾率均等。。。124DLA模型深化了人們對非平衡生長現(xiàn)象的認識，被廣泛應(yīng)用于研究時間依賴的生長、擴散、輸運、吸附等物理和化學(xué)過程。同時，DLA模型也不斷地被推廣和發(fā)展。例如，在真實的生長過程中，粒子不斷被沉積到基底表面，在表面上同時有多個粒子在隨機行走。在行走的過程中，可能會產(chǎn)生新的聚集中心。這就導(dǎo)致在基底表面會形成多個集團構(gòu)型。此外，粒子與其他粒子接觸之后可能也不會立即穩(wěn)定下來，而是有一定概率離開或者沿所形成的的集團的邊界行走。當襯底上形成了集團構(gòu)型，新的粒子可能沉積在集團之上，這樣逐漸在襯底上形成三維的集團構(gòu)型。這些推廣和擴展為理解真實系統(tǒng)的生長、擴散等提供了堅實的理論基礎(chǔ)。動力學(xué)蒙特卡洛方法是模擬這類問題的強有力工具。IsingModelSpinmodelforaferromegnetCollectionofMagneticmomentsSpinangularmomentumThesimplestIsingmodelassumesaninteractiononlybetweennearestneighborsJ(>0)isexchangeconstantAssumethateachspinisabletopointalongeither+zor–z;(upordown).Forconveniencewetaketobesi=1or-1.Theenergyofthespinsystemislowestifallofthespinsareparalleltooneanother.125126早在1920年，德國物理學(xué)家WilhelmLenz為解釋鐵磁相變，提出一個包含小箭頭的網(wǎng)格的簡單模型。1924年Lenz的學(xué)生ErnstIsing證明，當空間維數(shù)為1時，模型沒有相變。之后，眾多學(xué)者對二維Ising模型進行了研究，并發(fā)展了平均場理論。1944年，美國物理學(xué)家LarsOnsager發(fā)表了二維Ising模型的嚴格解。然而，精確求解三維Ising模型仍是困擾物理學(xué)家的一個未解難題。IsingModelDisorderingeffectoftemperatureAssumethatthespinsystemisinequilibriumwithaheatbathatT,sothatovertimedifferentspinsflipbackandforth,andthesystemwillthereby“move”intodifferentspinconfigurations.Forasysteminequilibriumwithaheatbath,theprobabilityoffindingthesysteminanyparticularstateisproportionaltotheBoltzmannfactorTheprobabilityoffindingthesysteminstate(microstate)Thereare2NdifferentpossiblemicrostatesofasystemwithparticleN.Themeasuredmagnetization127Mean-FieldTheoryofIsingModelAllspinsmusthavethesameaverageproperties.ThetotalmagneticationatTforsystemofNspinsisIfwecancalculate<si>,wehaveMimmediately.Anexactcomputationof<si>requirestheprobabilitiesofallpossiblemicrostates.IfweaddamagneticfieldH128Mean-FieldTheoryofIsingModelForasystemwithasinglespin129Mean-FieldTheoryofIsingModelThemeanfieldapproximationisbasedontheassumptionthattheinteractionofaspinsiwithitsneighboringspinsisequivalenttoaneffectivemagneticfieldactingonsi.H=0130Mean-FieldTheoryofIsingModel<S>Freeenergyvs<s>131Mean-FieldTheoryofIsingModelForsmallx132MonteCarloMethodforIsingModelSimulationofSpinflipsAspinischosenandtheenergyrequiredtomakeitflipEflipiscalculated.IfEflip<0,thespinisflippedandthesystemmovesintoadifferentmicrostate.IfEflip>0,whattodo?ArandomnumberRthatisdistributeduniformlybetween0and1isgenerated.IfP>R,thespinisflipped;otherwise,thespinisleftundisturbed.Metropolisalgorithm133MonteCarloalgorthmfortheIsingmodelonanLxLsquarelattice134

135L=10,H=0,Jasenergyunit,sotheunitofTisJ/kB136MonteCarloSimulationforIsingModel蒙特卡洛模擬Ising模型能量（左）和磁化強度（右）隨溫度的演化。以上結(jié)果是1000次的平均。137在Metropolis算法中，當系統(tǒng)接近或者處于平衡態(tài)時，產(chǎn)生新的構(gòu)型或者狀態(tài)的概率變得非常小。為了提高模擬效率，Bortz，Karlos和Lebowitz于1975年提出N-foldway算法來模擬二維Ising模型。簡單來說，在二維正方格子上，根據(jù)每個自旋及其近鄰自旋的取向，所有自旋可以分為10類構(gòu)型，每一類都可以計算出確定的反轉(zhuǎn)率N-foldway算法138構(gòu)型中心自旋4個近鄰自旋狀態(tài)中心自旋的反轉(zhuǎn)率1

2

，

，

，

3

，

，

，

，

4

，

，

，

5

6

7

，

，

，

8

，

，

，

，

9

，

，

，

10

139

N-foldway算法的特點140對于傳統(tǒng)的蒙特卡洛方法，一旦系統(tǒng)達到平衡態(tài)，時間尺度并沒有明確的物理意義，因此也沒有必要確定每一步所對應(yīng)的物理時間尺度。但是，對于一個處于非平衡態(tài)的系統(tǒng)，系統(tǒng)的演化以及動力學(xué)過程就需要考慮相應(yīng)的時間尺度，而基于N-foldway算法的蒙特卡洛模擬方法能夠?qū)γ商乜迥M的時間尺度與真實的時間尺度之間建立聯(lián)系。因此，基于N-foldway算法發(fā)展起來的蒙特卡洛模擬方法不僅能夠模擬處于平衡態(tài)的系統(tǒng)，還能夠模擬處于遠離平衡態(tài)下系統(tǒng)的動力學(xué)過程。20世紀90年代，該算法被廣泛應(yīng)用于研究晶體核薄膜生長、表面擴散和吸附、催化等很多領(lǐng)域。基于該算法的蒙特卡洛模擬方法也被統(tǒng)稱為動力學(xué)蒙特卡洛方法（kineticMonteCarlo）。動力學(xué)蒙特卡洛方法（kineticMonteCarlo）141動力學(xué)蒙特卡羅(KineticMC，KMC)在原子模擬領(lǐng)域內(nèi)，分子動力學(xué)(moleculardynamics,MD)具有突出的優(yōu)勢。它可以非常精確的描述體系演化的軌跡。一般情況下MD的時間步長在飛秒(fs)量級，因此足以追蹤原子振動的具體變化。但是這一優(yōu)勢同時限制了MD在大時間尺度模擬上的應(yīng)用。現(xiàn)有的計算條件足以支持MD到10ns，運用特殊的算法可以達到10μs的尺度。即便如此，很多動態(tài)過程，如表面生長或材料老化等，時間跨度均在s以上，大大超出了MD的應(yīng)用范圍。有什么方法可以克服這種局限呢？當體系處于穩(wěn)定狀態(tài)時，我們可以將其描述為處于維勢能函數(shù)面的一個局域極小值（阱底）處。有限溫度下，雖然體系內(nèi)的原子不停的進行熱運動，但是絕大部分時間內(nèi)原子都是在勢能阱底附近振動。偶然情況下體系會越過不同勢阱間的勢壘從而完成一次“演化”，這類小概率事件才是決定體系演化的重點。142KMC-續(xù)因此，如果我們將關(guān)注點從“原子”升格到“體系”，同時將“原子運動軌跡”粗化為“體系組態(tài)躍遷”，那么模擬的時間跨度就將從原子振動的尺度提高到組態(tài)躍遷的尺度。這是因為這種處理方法擯棄了與體系穿越勢壘無關(guān)的微小振動，而只著眼于體系的組態(tài)變化。因此，雖然不能描繪原子的運動軌跡，但是作為體系演化，其“組態(tài)軌跡”仍然是正確的。此外，因為組態(tài)變化的時間間隔很長，體系完成的連續(xù)兩次演化是獨立的，無記憶的，所以這個過程是一種典型的馬爾可夫過程(Markovprocess)，即體系從組態(tài)到組態(tài)，這一過程只與其躍遷速率有關(guān)。143KMC-續(xù)如果精確地知道，我們便可以構(gòu)造一個隨機過程，使得體系按照正確的軌跡演化。這里``正確''的意思是某條給定演化軌跡出現(xiàn)的幾率與MD模擬結(jié)果完全一致(假設(shè)我們進行了大量的MD模擬，每次模擬中每個原子的初始動量隨機給定)。這種通過構(gòu)造隨機過程研究體系演化的方法即為動力學(xué)蒙特卡洛方法(kineticMonteCarlo,KMC)KineticMonteCarlo:ACoarse-Grained,Atomistic,Lattice-BasedTechniqueforCondensed-MatterDynamics∞ContinuumEquationsKMCMDTime(s)10-1510-1210-910-6

∞Length(m)10-610-810-11KineticMonteCarlo:Coarse-GrainingMDRareEvents:MDofCoonCu(001):TheWholeTrajectoryKMC:Coarse-GrainedHopsRotateRareEventsInfrequenttransitionsfromalocalminimumtoanotherlocalminimumTransitionstatetheory

tounderstandchemicalreactionrates:(MichaelPolanyi&HenryEyringin1920’s&1930’s)MasterEquation:ProbabilitytobeatStateatTimet:TransitionProbabilityperUnitTimefromtoKineticMonteCarlosimulationKMCTransitionProbabilitiesandDetailedBalance

Energy

U*

UiUfUifbTSTSatisfiesDetailedBalanceandKineticsp(if)=

0

exp(Ubif/kBT)MetropolisMCSatisfiesDetailedBalance,ButNotKineticsAGenericKMCAlgorithmInitializeLatticeFinished?IdentifyAllProcessesandRatesRiDoProcessa,IncrementTimeYesNoChooseaProcessa01?=a1)(iiP?-=11)(aiiP1.Surfacediffusion2.Surfacegrowth3.Vacancydiffusioninalloys4.Coarseningofdomainevolution5.Defectmobilityandclusteringinionorneutronirradiatedsolids6.ViscoelasticityofphysicallycrosslinkednetworksApplicationtoKMCSimulationsofThin-FilmEpitaxyDeposition

AggregationNucleationTerraceDiffusion

EdgeDiffusionGerdBinnigHeinrichRohrerTheysharedthehalfoftheNobelPrizeinPhysicsin1986forthedesignofthescanningtunnelingmicroscope(STM).(theotherhalfofthePrizewasawardedtoErnstRuskaforhisachievementsinelectronopticsincludingthedesignofthefirstelectronmicroscope).Vsample=1.00VItunne