蒙特卡羅方法課件

蒙特卡羅方法湘潭大學(xué)物理學(xué)專業(yè)基礎(chǔ)課計(jì)算物理及其應(yīng)用

材料與光電物理學(xué)院孟利軍2010年4月計(jì)算機(jī)模擬隨機(jī)性模擬方法確定性模擬方法直接蒙特卡羅模擬蒙特卡羅積分Metropolis蒙特卡羅模擬蒙特卡羅通過(guò)不斷產(chǎn)生隨機(jī)數(shù)序列來(lái)模擬過(guò)程通過(guò)數(shù)值求解一個(gè)個(gè)粒子的運(yùn)動(dòng)方程來(lái)模擬整個(gè)系統(tǒng)的行為計(jì)算機(jī)模擬和蒙特卡羅方法間接蒙特卡羅模擬目錄第一節(jié)蒙特卡羅方法概述第二節(jié)隨機(jī)數(shù)與偽隨機(jī)數(shù)第三節(jié)隨機(jī)變量的抽樣第四節(jié)蒙特卡羅方法的應(yīng)用實(shí)例§1蒙特卡羅方法概述---基本思想基本思想：

針對(duì)待求問(wèn)題，根據(jù)物理現(xiàn)象本身的統(tǒng)計(jì)規(guī)律，或人為構(gòu)造一合適的依賴隨機(jī)變量的概率模型，使某些隨機(jī)變量的統(tǒng)計(jì)量為待求問(wèn)題的解，進(jìn)行大統(tǒng)計(jì)量（N→∞）的統(tǒng)計(jì)實(shí)驗(yàn)方法或計(jì)算機(jī)隨機(jī)模擬方法。

理論依據(jù)：

大數(shù)定理：均勻分布的算術(shù)平均收斂于真值中心極限定理：置信水平下的統(tǒng)計(jì)誤差兩個(gè)例子：

Buffen投針實(shí)驗(yàn)求π射擊問(wèn)題（打靶游戲）Buffon投針實(shí)驗(yàn)(1777年）求π:３．各向同性隨機(jī)投針，則夾角α在[０,π]均勻分布，所以：４．設(shè)投針N次，相交次數(shù)為M，則相交概率的期望值：N→∞大數(shù)定理§1蒙特卡羅方法概述---基本思想２．針與平行線垂直方向夾角為α，則相交概率為：１．平行線間距＝針長(zhǎng)＝s一些人進(jìn)行了實(shí)驗(yàn)，其結(jié)果列于下表：實(shí)驗(yàn)者年份投計(jì)次數(shù)π的實(shí)驗(yàn)值沃爾弗(Wolf)185050003.1596史密思(Smith)185532043.1553?？怂?Fox)189411203.1419拉查里尼(Lazzarini)190134083.1415929§1蒙特卡洛方法的基本思想§1蒙特卡羅方法概述---基本思想1.設(shè)r表示射擊運(yùn)動(dòng)員彈著點(diǎn)到靶心的距離，ｇ(r)表示擊中r處相應(yīng)的得分?jǐn)?shù)（環(huán)數(shù)），f(r)為該運(yùn)動(dòng)員彈著點(diǎn)的分布密度函數(shù)，它們反映運(yùn)動(dòng)員的射擊水平。該運(yùn)動(dòng)員的射擊成績(jī)?yōu)椋海?用概率語(yǔ)言來(lái)說(shuō)，<g>是隨機(jī)變量ｇ(r)的數(shù)學(xué)期望，即３.假設(shè)該運(yùn)動(dòng)員進(jìn)行了Ｎ次射擊，每次射擊的彈著點(diǎn)依次為r1，r2，…，rN，則Ｎ次得分g(r1)，g(r2)，…，g(rN)的算術(shù)平均值代表了該運(yùn)動(dòng)員的成績(jī)射擊問(wèn)題（打靶游戲）§1蒙特卡羅方法概述---基本思想４.用Ｎ次試驗(yàn)所得成績(jī)的算術(shù)平均值作為數(shù)學(xué)期望<g>的估計(jì)值。例如，設(shè)射擊運(yùn)動(dòng)員的彈著點(diǎn)分布為環(huán)數(shù)78910概率0.10.10.30.5用計(jì)算機(jī)作隨機(jī)試驗(yàn)（射擊）的方法為，選取一個(gè)隨機(jī)數(shù)ξ，按右邊所列方法判斷得到成績(jī)。這樣，就進(jìn)行了一次隨機(jī)試驗(yàn)（射擊），得到了一次成績(jī)ｇ(r)，作Ｎ次試驗(yàn)后，得到該運(yùn)動(dòng)員射擊成績(jī)的近似值

§1蒙特卡羅方法概述---基本思想

收斂性:大數(shù)定理作為所求解的近似值。由大數(shù)定律可知，如果X1，X2，…，XN獨(dú)立同分布，且具有有限期望值（E(X)<∞），則

即隨機(jī)變量X的簡(jiǎn)單子樣的算術(shù)平均值，當(dāng)子樣數(shù)Ｎ充分大時(shí)，以概率1收斂于它的期望值E(X)。由前面介紹可知，蒙特卡羅方法是由隨機(jī)變量X的簡(jiǎn)單子樣X(jué)1，X2，…，XN的算術(shù)平均值：§1蒙特卡洛方法概述---大數(shù)定理f(x)是X的分布密度函數(shù)。則當(dāng)N充分大時(shí)，有如下的近似式蒙特卡羅方法的近似值與真值的誤差問(wèn)題，概率論的中心極限定理給出了答案。該定理指出，如果隨機(jī)變量序列X1，X2，…，XN獨(dú)立同分布，且具有有限非零的方差σ2，即§1蒙特卡洛方法概述---中心極限定理其中α稱為置信度，1－α稱為置信水平。這表明，不等式近似地以概率1－α成立，且誤差收斂速度的階為：O(N-1/2)上式中λα與置信度α是一一對(duì)應(yīng)的，根據(jù)問(wèn)題的要求確定出置信水平后，查標(biāo)準(zhǔn)正態(tài)分布表，就可以確定出λα。通常，蒙特卡羅方法的誤差ε定義為給出幾個(gè)常用的α與λα

的數(shù)值：α0.50.050.003λα

0.67451.963兩點(diǎn)說(shuō)明：(1)MC方法的誤差為概率誤差，這與其他數(shù)值計(jì)算方法是有區(qū)別的。(2)誤差中的均方差σ是未知的，必須使用其估計(jì)值來(lái)代替，在計(jì)算所求量的同時(shí)，可計(jì)算出估計(jì)值。§1蒙特卡洛方法概述---中心極限定理（2）減小估計(jì)的均方差σ，比如降低一半，那誤差就減小一半，這相當(dāng)于N增大四倍的效果。

減小方差的各種技巧:（1）增大試驗(yàn)次數(shù)N。在σ固定的情況下，要把精度提高一個(gè)數(shù)量級(jí)，試驗(yàn)次數(shù)N需增加兩個(gè)數(shù)量級(jí)。因此，單純?cè)龃驨不是一個(gè)有效的辦法。顯然當(dāng)給定置信度α（λα）后，誤差ε由σ和N決定。要減小ε：§1蒙特卡洛方法該概述---減小誤差

效率:一般來(lái)說(shuō)，降低方差的技巧，往往會(huì)使觀察一個(gè)子樣的時(shí)間增加。在固定時(shí)間內(nèi)，使觀察的樣本數(shù)減少。所以，一種方法的優(yōu)劣，需要由方差和觀察一個(gè)子樣的費(fèi)用（使用計(jì)算機(jī)的時(shí)間）兩者來(lái)衡量。這就是MC方法中效率的概念。它定義為σ2c，其中c是觀察一個(gè)子樣的平均費(fèi)用。顯然σ2c越小，方法越有效。（1）能夠比較逼真地描述具有隨機(jī)性質(zhì)的事物的特點(diǎn)及物理實(shí)驗(yàn)過(guò)程§1蒙特卡羅方法概述---MC優(yōu)點(diǎn)從這個(gè)意義上講，蒙特卡羅方法可以部分代替物理實(shí)驗(yàn)，甚至可以得到物理實(shí)驗(yàn)難以得到的結(jié)果。用蒙特卡羅方法解決實(shí)際問(wèn)題，可以直接從實(shí)際問(wèn)題本身出發(fā)，而不從方程或數(shù)學(xué)表達(dá)式出發(fā)。它具有直觀、形象的特點(diǎn)。（2）受幾何條件限制小計(jì)算s維空間中的任一區(qū)域Ds上的積分：無(wú)論區(qū)域Ds的形狀多么特殊，只要能給出描述Ds的幾何特征的條件，就可以從Ds中均勻產(chǎn)生N個(gè)點(diǎn)：因此，在具有隨機(jī)性質(zhì)的問(wèn)題中，如考慮的系統(tǒng)形狀很復(fù)雜，難以用一般數(shù)值方法求解，而使用蒙特卡羅方法，不會(huì)有原則上的困難。其中Ds為區(qū)域Ds的體積。這是數(shù)值方法難以作到的。得到積分的近似值：§1蒙特卡羅方法概述---MC優(yōu)點(diǎn)（3）收斂速度與問(wèn)題的維數(shù)無(wú)關(guān)由誤差定義可知，在給定置信水平情況下，MC方法的誤差為O(N-1/2)

，與問(wèn)題本身的維數(shù)無(wú)關(guān)。維數(shù)的變化，只引起抽樣時(shí)間及估計(jì)量計(jì)算時(shí)間的變化，不影響誤差。這一特點(diǎn)，決定了蒙特卡羅方法對(duì)多維問(wèn)題的適應(yīng)性。而一般數(shù)值方法，比如計(jì)算定積分時(shí)，計(jì)算時(shí)間隨維數(shù)的冪次方而增加，而且由于分點(diǎn)數(shù)與維數(shù)的冪次方成正比，需占用相當(dāng)數(shù)量的計(jì)算機(jī)內(nèi)存，這些都是一般數(shù)值方法計(jì)算高維積分時(shí)難以克服的問(wèn)題。（4）具有同時(shí)計(jì)算多個(gè)方案與多個(gè)未知量的能力（2）使用蒙特卡羅方法還可以同時(shí)得到若干個(gè)所求量。（1）對(duì)于那些需要計(jì)算多個(gè)方案的問(wèn)題，使用蒙特卡羅方法有時(shí)不需要像常規(guī)方法那樣逐個(gè)計(jì)算，而可以同時(shí)計(jì)算所有的方案，其全部計(jì)算量幾乎與計(jì)算一個(gè)方案的計(jì)算量相當(dāng)。例如對(duì)于屏蔽層為均勻介質(zhì)的幾何平板，要計(jì)算若干種厚度的穿透概率時(shí)，只需計(jì)算最厚的一種情況，其他厚度的穿透概率在計(jì)算最厚一種情況時(shí)稍加處理便可同時(shí)得到。例如在模擬粒子過(guò)程中，可以同時(shí)得到不同區(qū)域的通量、能譜、角分布等，而不像常規(guī)方法那樣，需要逐一計(jì)算所求量。§1蒙特卡羅方法概述---MC優(yōu)點(diǎn)（5）誤差容易確定根據(jù)蒙特卡羅方法的誤差公式，可以在計(jì)算所求量的同時(shí)計(jì)算出誤差（6）程序結(jié)構(gòu)簡(jiǎn)單，易于實(shí)現(xiàn)在計(jì)算機(jī)上進(jìn)行蒙特卡羅方法計(jì)算時(shí)，程序結(jié)構(gòu)簡(jiǎn)單，分塊性強(qiáng)，易于實(shí)現(xiàn)。（1）收斂速度慢（2）誤差具有概率性蒙特卡羅方法的收斂速度為O(N-1/2)

，一般不容易得到精確度較高的近似結(jié)果。對(duì)于維數(shù)少（三維以下）的問(wèn)題，不如其他方法好。由于蒙特卡羅方法的誤差是在一定置信水平下估計(jì)的，所以它的誤差具有概率性，而不是一般意義下的誤差?！?蒙特卡羅方法概述---MC缺點(diǎn)（3）計(jì)算結(jié)果與系統(tǒng)大小有關(guān)例如在粒子輸運(yùn)問(wèn)題中:經(jīng)驗(yàn)表明，只有當(dāng)系統(tǒng)的大小與粒子的平均自由程可以相比較時(shí)（一般在十個(gè)平均自由程左右），蒙特卡羅方法計(jì)算的結(jié)果較為滿意。但對(duì)于大系統(tǒng)或小概率事件的計(jì)算問(wèn)題，計(jì)算結(jié)果往往比真值偏低。在使用蒙特卡羅方法時(shí)，可以考慮把蒙特卡羅方法與解析（或數(shù)值）方法相結(jié)合，取長(zhǎng)補(bǔ)短，既能解決解析（或數(shù)值）方法難以解決的問(wèn)題，也可以解決單純使用蒙特卡羅方法難以解決的問(wèn)題。真隨機(jī)數(shù)：不可重復(fù)，物理方法產(chǎn)生。隨機(jī)數(shù)是實(shí)現(xiàn)由已知分布抽樣的基本量，在由已知分布的抽樣過(guò)程中，將隨機(jī)數(shù)作為已知量，用適當(dāng)?shù)臄?shù)學(xué)方法可以由它產(chǎn)生具有任意已知分布的簡(jiǎn)單子樣。由具有已知分布的總體中抽取簡(jiǎn)單子樣，在蒙特卡羅方法中占有非常重要的地位?？傮w和子樣的關(guān)系，屬于一般和個(gè)別的關(guān)系，或者說(shuō)屬于共性和個(gè)性的關(guān)系?！?隨機(jī)數(shù)的產(chǎn)生和檢驗(yàn)蒙特卡羅模擬就是邊產(chǎn)生隨機(jī)數(shù)邊進(jìn)行隨機(jī)模擬的方法準(zhǔn)隨機(jī)數(shù)：不具隨機(jī)性質(zhì)，只要處理問(wèn)題能得到正確結(jié)果。如放射性衰變，電子設(shè)備的熱噪音，宇宙射線的觸發(fā)時(shí)間等。偽隨機(jī)數(shù)：可重復(fù)，數(shù)學(xué)方法產(chǎn)生，必須通過(guò)統(tǒng)計(jì)檢驗(yàn)區(qū)分：數(shù)列的隨機(jī)性和隨機(jī)數(shù)的分布一個(gè)完美的隨機(jī)數(shù)序列可能具有某種分布（如均勻分布，高斯分布等），但具有某種分布的數(shù)列卻可能完全不是隨機(jī)的。良好統(tǒng)計(jì)分布，容易實(shí)現(xiàn)，效率高，周期長(zhǎng)，可移植性好等§2隨機(jī)數(shù)的產(chǎn)生和檢驗(yàn)分布函數(shù)為：最簡(jiǎn)單、最基本，也是最重要的隨機(jī)數(shù)是在單一區(qū)間[0,1]上的均勻分布的隨機(jī)數(shù)，其分布密度函數(shù)為由于隨機(jī)數(shù)在蒙特卡羅方法中占有極其重要的位置，我們用專門的符號(hào)ξ表示。用ξ1，ξ2

，…代表相互獨(dú)立且具有相同單位均勻分布的隨機(jī)數(shù)序列。獨(dú)立性、均勻性是隨機(jī)數(shù)必備的兩個(gè)特點(diǎn)。如：擲篩子游戲，投擲硬幣

用來(lái)作為隨機(jī)數(shù)發(fā)生器的物理源主要有兩種：一種是根據(jù)放射性物質(zhì)的放射性，另一種是利用計(jì)算機(jī)的固有噪聲。用物理方法產(chǎn)生的隨機(jī)數(shù)序列無(wú)法重復(fù)實(shí)現(xiàn)，不能進(jìn)行程序復(fù)算，給驗(yàn)證結(jié)果帶來(lái)很大困難。而且，需要增加隨機(jī)數(shù)發(fā)生器和電路聯(lián)系等附加設(shè)備，費(fèi)用昂貴。因此，該方法也不適合在計(jì)算機(jī)上使用。§2偽隨機(jī)數(shù)的產(chǎn)生和檢驗(yàn)---物理方法在計(jì)算機(jī)上用物理方法產(chǎn)生隨機(jī)數(shù)的基本原理是：利用某些物理現(xiàn)象，在計(jì)算機(jī)上增加某些特殊設(shè)備，可以在計(jì)算機(jī)上直接產(chǎn)生隨機(jī)數(shù)。這些特殊設(shè)備稱為隨機(jī)數(shù)發(fā)生器。（1）馮·諾伊曼平方取中法遞推公式：以十進(jìn)制數(shù)為例，平方取中法是把一個(gè)2S位的十進(jìn)制自平方后，去頭截尾只保留中間2S個(gè)數(shù)字，然后用102S來(lái)除，這樣就可以得到在[0,1]上均勻分布的偽隨機(jī)數(shù)序列。例如，設(shè)十進(jìn)制數(shù)的2S=4，并取x1=6406，則有：相應(yīng)的偽隨機(jī)數(shù)序列是0.6406，0.3680，0.1354，0.8333，0.4388等§2偽隨機(jī)數(shù)的產(chǎn)生和檢驗(yàn)---數(shù)學(xué)方法具有周期性，有些數(shù)甚至?xí)o接著重復(fù)出現(xiàn)，很少使用?！?偽隨機(jī)數(shù)的產(chǎn)生和檢驗(yàn)---數(shù)學(xué)方法由Lehmer在1951年提出來(lái)的，它的一般形式是：對(duì)于任一初始值x0，偽隨機(jī)數(shù)序列由下面遞推公式確定：（2）Lehmer線性同余法例如乘同余法x0稱為種子，改變它的值就得到基本序列的不同區(qū)段隨機(jī)數(shù)。a---乘子，c---增量，m---模乘同余法具有在計(jì)算機(jī)上容易實(shí)現(xiàn)、快速等特點(diǎn)，所以乘同余法已被廣泛采用。偽隨機(jī)數(shù)的均勻性偽隨機(jī)數(shù)的獨(dú)立性均勻性是指在[0,1]區(qū)間內(nèi)等長(zhǎng)度子區(qū)間中隨機(jī)數(shù)的數(shù)量是一樣的。按先后順序出現(xiàn)的隨機(jī)數(shù)中，每個(gè)隨機(jī)數(shù)的取值與其相距一定間隔的隨機(jī)數(shù)取值之間無(wú)關(guān)?！?偽隨機(jī)數(shù)序列的統(tǒng)計(jì)檢驗(yàn)判斷偽隨機(jī)數(shù)序列是否滿足均勻和相互獨(dú)立的要求，要靠統(tǒng)計(jì)檢驗(yàn)的方法實(shí)現(xiàn)。對(duì)于偽隨機(jī)數(shù)的統(tǒng)計(jì)檢驗(yàn)，一般包括兩大類：均勻性檢驗(yàn)和獨(dú)立性檢驗(yàn)偽隨機(jī)數(shù)的均勻性將區(qū)間[0,1]分為K個(gè)子區(qū)間，統(tǒng)計(jì)隨機(jī)數(shù)落在第k

個(gè)子區(qū)間的實(shí)際頻數(shù)nk，它應(yīng)當(dāng)趨近于理論頻數(shù)mk計(jì)算統(tǒng)計(jì)量如果χ2

值很大，表示遠(yuǎn)遠(yuǎn)偏離理想值，因此要求χ2值盡可能小，但如果它趨于0則有可能N已進(jìn)入循環(huán)。通常求和中的每一項(xiàng)的大小約為1，因此χ2的值約為K

。§2隨機(jī)數(shù)的產(chǎn)生和檢驗(yàn)---統(tǒng)計(jì)檢驗(yàn)K的取值不能太大也不能太小，太大反映不出“小區(qū)間”的均勻性，太小反映不出“大區(qū)間”的均勻性。限制條件（1）概率論中的Pearson定理說(shuō)明，上式的極限概率分布是χ2分布給出了χ2≤x的概率。整數(shù)ν是系統(tǒng)自由度表示獨(dú)立測(cè)量的次數(shù)，由于存在一個(gè)限制條件，故ν=K-1給出了χ2>x的概率余函數(shù)§2隨機(jī)數(shù)的產(chǎn)生和檢驗(yàn)---統(tǒng)計(jì)檢驗(yàn)因此，當(dāng)給定顯著水平α(或置信度1?α)后，由方程Q(χ2|υ)=α或P(χ2|υ）=1?α解出χα值，或從已有的χ2表中查得χα值，如果由（1）式計(jì)算出來(lái)的χ小于χα，則認(rèn)為在此置信度下，偽隨機(jī)數(shù)序列在[0,1]中是均勻分布的。（1）順序相關(guān)法它用相鄰兩個(gè)隨機(jī)數(shù)的自相關(guān)函數(shù)（或相關(guān)系數(shù)）來(lái)標(biāo)識(shí)偽隨機(jī)數(shù)序列的獨(dú)立性情況，間距為l的自相關(guān)函數(shù)是相關(guān)系數(shù)越小，獨(dú)立性越好?！?隨機(jī)數(shù)產(chǎn)生及檢驗(yàn)---獨(dú)立性統(tǒng)計(jì)檢驗(yàn)（2）多維頻率檢驗(yàn)（1）將偽隨機(jī)數(shù)序列用任意一種辦法進(jìn)行組合，每S個(gè)隨機(jī)數(shù)作為S維空間中的一個(gè)點(diǎn)的坐標(biāo)值，于是可以構(gòu)成一個(gè)點(diǎn)序列。（2）把S維空間中的單位方體分成為K個(gè)子方體，方體邊長(zhǎng)（3）統(tǒng)計(jì)落在第k個(gè)子方體中的實(shí)際頻數(shù)nk，它應(yīng)當(dāng)趨近于理論頻數(shù)：例如將2N

個(gè)隨機(jī)數(shù)序列分為兩組：{ξ1

,

ξ

3,…}和{ξ

2,ξ

4,…}，分別作為平面中N

個(gè)點(diǎn)的x和y坐標(biāo)值。在xy

平面中作K0×K0個(gè)小正方形網(wǎng)格區(qū)域，落在第(i,j)個(gè)網(wǎng)格區(qū)域中的實(shí)際頻數(shù)為nij

，則§2隨機(jī)數(shù)產(chǎn)生及檢驗(yàn)---獨(dú)立性統(tǒng)計(jì)檢驗(yàn)用連續(xù)的兩個(gè)隨機(jī)數(shù)作為點(diǎn)(x,y)的坐標(biāo)作圖可以直觀看出隨機(jī)數(shù)之間的關(guān)聯(lián)性。顯然左邊是不好的隨機(jī)數(shù)。左上顯示出條帶結(jié)構(gòu)，左下則是規(guī)則網(wǎng)格結(jié)構(gòu)。我們也可以把偽隨機(jī)數(shù)分為三列、四列等，用相似的方法進(jìn)行多維獨(dú)立性檢驗(yàn)§3隨機(jī)變量的抽樣---離散型隨機(jī)變量X：｛x1,x2,x3,···,xN｝例如：x可取3個(gè)值x1,x2,x3，它們出現(xiàn)的幾率分別為2/8,5/8,1/8，則隨機(jī)數(shù)小于2/8時(shí)實(shí)現(xiàn)x1

，在區(qū)間2/8,7/8中時(shí)實(shí)現(xiàn)x2

，大于7/8時(shí)實(shí)現(xiàn)x3.概率密度f(wàn)：｛p1,p2,p3,···,pN｝如果從[0,1]區(qū)間中均勻抽樣得到的隨機(jī)數(shù)ξ滿足下式時(shí)則物理量x取值為xn

。實(shí)際需要的大多數(shù)隨機(jī)變量并不是[0,1]區(qū)間均勻分布的，而是有各種不同形式分布密度函數(shù)的隨機(jī)變量。因此，對(duì)不均勻的隨機(jī)變量抽樣的關(guān)鍵問(wèn)題是如何從均勻分布的偽隨機(jī)變量樣本中，抽取符合所要求的分布密度函數(shù)的簡(jiǎn)單子樣。§3隨機(jī)變量的抽樣---離散型3MeV光子入射屏蔽鉛板的全吸收反射過(guò)程反應(yīng)類型X：光電效應(yīng)康普頓散射電子對(duì)產(chǎn)生反應(yīng)截面σ：σ1

σ2

σ3

反應(yīng)幾率f：ε1=σ1/σε2=σ2/σε3=σ3/σ歸一化：例如Poisson分布是離散型分布對(duì)此分布進(jìn)行抽樣得到第n個(gè)事件發(fā)生的條件為電子對(duì)產(chǎn)生光電效應(yīng)康普頓效應(yīng)NNYY設(shè)連續(xù)型變量x在區(qū)間[a,b]中取值，可視為將上述的離散情形取連續(xù)極限:顯然ξ(a)=0,ξ(b)=1且是單調(diào)增要求變量x，可由上式解析反解出x(ξ)的函數(shù)表達(dá)式，即求反函數(shù)。這對(duì)一些簡(jiǎn)單的幾率密度函數(shù)解析表達(dá)式是很容易做到的。如粒子隨機(jī)運(yùn)動(dòng)的自由程分布為指數(shù)分布：則求反函數(shù)后得§3隨機(jī)變量的抽樣---連續(xù)型累積函數(shù)注意（1-ξ）和ξ同樣服從[0，1]的均勻分布§3隨機(jī)變量的抽樣---變換抽樣法一維：變換抽樣法的基本思想是將一個(gè)比較復(fù)雜的分布p(x)的抽樣，變換為已知的簡(jiǎn)單分布g(y)的抽樣我們希望找到x?y

之間的對(duì)應(yīng)關(guān)系，使得幾率密度守恒：例如：黑體輻射的譜密度按頻率ω表示時(shí)為當(dāng)希望將譜密度用波長(zhǎng)λ=2πcω表示時(shí)，有§3隨機(jī)變量的抽樣---變換抽樣法顯然，當(dāng)g(y)取[0,1]均勻分布時(shí)，問(wèn)題即化為：尋找y(x)，使其導(dǎo)數(shù)為p(x)，然后在[0,1]區(qū)間中對(duì)變量y

抽樣得到均勻分布的隨機(jī)數(shù)，再由x(y)關(guān)系得到對(duì)應(yīng)幾率密度函數(shù)p(x)的隨機(jī)抽樣x

。二維：有兩個(gè)變量x和y的聯(lián)合分布密度函數(shù)為p(x,y)，欲變換至變量u和v，它們的聯(lián)合分布密度函數(shù)為g(u,v)取聯(lián)合分布密度函數(shù)g(u,v)為均勻分布：§3隨機(jī)變量的抽樣---變換抽樣法則我們的任務(wù)就是尋找變換式x=x(u,v),y=y(u,v)，以使p(x,y)=|?(u,v)/?(x,y)|

對(duì)均勻隨機(jī)變量(u,v)進(jìn)行抽樣，代入變換式得x

和y的抽樣。對(duì)于Gauss正態(tài)幾率分布的抽樣通過(guò)代換可以只考慮簡(jiǎn)單形式的分布令極角坐標(biāo)系下的角度為2πv

，半徑為現(xiàn)在我們?cè)噲D通過(guò)一個(gè)兩維聯(lián)合分布的抽樣獲得該一維分布的抽樣。u和v都是[0,1]區(qū)間中的均勻分布的隨機(jī)抽樣，則變換關(guān)系式為§3隨機(jī)變量的抽樣---變換抽樣法Jacobi行列式即兩維分布正為兩個(gè)獨(dú)立分布之積。顯然抽樣x

或y都滿足正態(tài)分布?？梢?jiàn)，為了得到滿足一個(gè)復(fù)雜分布的隨機(jī)抽樣，這里用了兩個(gè)滿足簡(jiǎn)單分布的隨機(jī)數(shù)?？傻梅醋儞Q§3隨機(jī)變量的抽樣---變換抽樣法§4蒙特卡羅方法在積分計(jì)算中的應(yīng)用一維積分平均值法：作變換：得標(biāo)準(zhǔn)積分：§4-1蒙特卡洛方法在積分計(jì)算中的應(yīng)用直接抽樣法：在x的定義域[0,1]上均勻隨機(jī)取點(diǎn)，該均勻分布的隨機(jī)變量記為ξ，定義隨機(jī)變量η為：則有因此，只要抽取足夠多的隨機(jī)點(diǎn)，即當(dāng)n足夠大時(shí)，In就是積分I的一個(gè)無(wú)偏估計(jì)值。相應(yīng)的方差為：可見(jiàn)，當(dāng)f(x)在其定義域內(nèi)變化較大時(shí)，方差較大?！?-1蒙特卡洛方法在積分計(jì)算中的應(yīng)用重要抽樣法：當(dāng)f(x)在其定義域內(nèi)有顯著的起伏變化時(shí)，可采用重要抽樣法。可將積分稍作變換：適當(dāng)選取偏倚分布密度函數(shù)，使得f*(x)在定義域內(nèi)變化比較平坦。然后產(chǎn)生[0,1]區(qū)間分布密度函數(shù)為g(x)的隨機(jī)變量ξ’，定義：則有：偏倚分布密度函數(shù)§4-1蒙特卡洛方法在積分計(jì)算中的應(yīng)用相應(yīng)的方差為：蒙特卡洛計(jì)算結(jié)果的方差為：直接抽樣重要抽樣§4-1蒙特卡洛方法在積分計(jì)算中的應(yīng)用高維積分平均值法標(biāo)準(zhǔn)形式：在實(shí)際物理問(wèn)題中，被積函數(shù)在超立方體區(qū)域內(nèi)可能強(qiáng)烈變化。若在積分區(qū)域內(nèi)均勻抽樣，積分貢獻(xiàn)可能主要來(lái)自少數(shù)僅僅只有幾個(gè)蒙特卡洛投點(diǎn)的小區(qū)域，從而導(dǎo)致很大的統(tǒng)計(jì)誤差。所以采用重要抽樣法，使得隨機(jī)點(diǎn)更多地投在取值大的區(qū)間§4-1蒙特卡洛方法在積分計(jì)算中的應(yīng)用選取偏倚分布密度函數(shù)，并定義（要求方差較?。┯校喊凑掌蟹植济芏群瘮?shù)在區(qū)域抽樣N個(gè)子樣則：積分的近似值§4-1蒙特卡洛方法在積分計(jì)算中的應(yīng)用一維積分的擲點(diǎn)法一維標(biāo)準(zhǔn)積分：定義：則有：在單位正方形內(nèi)投N個(gè)點(diǎn)，落在曲線f(x)下的有M個(gè)，則由于對(duì)y的積分可以解析計(jì)算，故此法的誤差較平均值法大?！?-2蒙特卡洛方法在量子力學(xué)中的應(yīng)用量子力學(xué)回顧薛定諤方程定義：費(fèi)曼傳播子量子力學(xué)的基本理論告訴我們：系統(tǒng)的所有信息，包括基態(tài)、激發(fā)態(tài)的能量、波函數(shù)等論可由費(fèi)曼傳播子給出。特別是和基態(tài)有關(guān)的信息，可以很方便的得到?！?-2蒙特卡洛方法在量子力學(xué)中的應(yīng)用例如：基態(tài)波函數(shù)的模方可表示為費(fèi)曼傳播子可以表示為路徑積分的形式，對(duì)于簡(jiǎn)單系統(tǒng)，即有其中，作用量§4-2蒙特卡洛方法在量子力學(xué)中的應(yīng)用路徑積分量子蒙特卡洛方法實(shí)際數(shù)值計(jì)算中，費(fèi)曼傳播子表達(dá)式為

實(shí)際計(jì)算中，Ｎ足夠大即可。積分常數(shù)時(shí)§4-2蒙特卡洛方法在量子力學(xué)中的應(yīng)用?。海袝r(shí)§4-2蒙特卡洛方法在量子力學(xué)中的應(yīng)用由此得到基態(tài)波函數(shù)：其中插入δ函數(shù)，得到§4-2蒙特卡洛方法在量子力學(xué)中的應(yīng)用采用Ｍetropolis方法計(jì)算基態(tài)波函數(shù)：首先，選擇任意的、連續(xù)Ｎ個(gè)時(shí)間間隔、且的一條路徑，計(jì)算相應(yīng)的能量。然后，再接著選一系列路徑，每條路徑與前一條路徑最多只有一個(gè)時(shí)刻（如）有不同的空間點(diǎn)。采用Ｍetropolis方法來(lái)確定滿足上述要求的新路徑。其中，將隨機(jī)定下的坐標(biāo)改變到的過(guò)渡概率為。其中，為兩條路徑的能量差。對(duì)于每條路徑，利用前述公式計(jì)算被積函數(shù)的估計(jì)值，并累加到求和之中。最終該求和所得值與抽樣路徑的總數(shù)相除得平均值，就得的數(shù)值結(jié)果。按上述方法，游走足夠多的步數(shù)后，我們就得到的值?！?-2蒙特卡洛方法在量子力學(xué)中的應(yīng)用變分量子蒙特卡洛方法對(duì)于任意的試探函數(shù)Ψ，其能量期望值滿足對(duì)于的計(jì)算，采用重要抽樣法。當(dāng)給定試探函數(shù)后，由Ｍetropolis方法產(chǎn)生分布的Ｎ個(gè)位形。對(duì)于每個(gè)位形，計(jì)算出相應(yīng)的局域能量，“局域能量”ε基態(tài)能量§4-2蒙特卡洛方法在量子力學(xué)中的應(yīng)用則：變分法步驟：（１）選擇一個(gè)物理上相對(duì)合理的基態(tài)試探波函數(shù)；（２）利用前述方法計(jì)算與之相應(yīng)的能量期望值；（３）改變?cè)囂胶瘮?shù)中的變分參數(shù)值，使試探函數(shù)改變一小量，記改變后的試探函數(shù)為，并計(jì)算相應(yīng)的能量期望值；（４）計(jì)算能量改變值，若改變量小于０，則接受試探函數(shù)，否則拒絕，并回到第三步；（５）反復(fù)循環(huán)，直至能量期望值不再有明顯變化為止。§4-2蒙特卡洛方法在量子力學(xué)中的應(yīng)用若上述循環(huán)至第Ｍ次終止，則格林函數(shù)量子蒙特卡洛方法一維擴(kuò)散方程：相應(yīng)的格林函數(shù)：§4-2蒙特卡洛方法在量子力學(xué)中的應(yīng)用

格林函數(shù)歸一，且與ｘ，ｔ無(wú)關(guān)。擴(kuò)散方程馬爾科夫過(guò)程（格林函數(shù)）（單步游走的概率分布）

分布§4-2蒙特卡洛方法在量子力學(xué)中的應(yīng)用Fokker-Planck方程：相應(yīng)的格林函數(shù)：

（Δｔ的一階近似）

構(gòu)造馬爾科夫鏈

分布力§4-2蒙特卡洛方法在量子力學(xué)中的應(yīng)用３Ｎ維多體定態(tài)薛定諤方程的基態(tài)解擴(kuò)散方程的定態(tài)解虛時(shí)薛定諤方程（?。?/p>

（具有勢(shì)函數(shù)的擴(kuò)散方程）擴(kuò)散項(xiàng)分支項(xiàng)格林函數(shù)：

（借用Ｄｉｒａｃ記號(hào)）§4-2蒙特卡洛方法在量子力學(xué)中的應(yīng)用上述格林函數(shù)正是量子力學(xué)中的時(shí)間演化算符在坐標(biāo)表象下的矩陣元，即費(fèi)曼傳播子。易知道當(dāng)，或者說(shuō)足夠大時(shí)，上述右邊的求和中只有基態(tài)才有貢獻(xiàn)，這個(gè)算符（時(shí)間演化算符）行為就如同作用在基態(tài)波函數(shù)上。解析計(jì)算格林函數(shù)，得到§4-2蒙特卡洛方法在量子力學(xué)中的應(yīng)用上面給出的是格林函數(shù)短時(shí)間近似結(jié)果。根據(jù)此結(jié)果，我們?cè)诟窳趾瘮?shù)蒙特卡洛模擬中，就必須進(jìn)行大量的短時(shí)間間隔的游走，最終使其分布近似滿足基態(tài)波函數(shù)。擴(kuò)散步游走分支步由游走到的權(quán)重需乘因子模擬效率不高。Ｒeynoids方法§4-3蒙特卡洛方法在統(tǒng)計(jì)物理中的應(yīng)用物理量的觀測(cè)值微觀粒子的某物理量在相空間的分布平均值高維積分相空間態(tài)矢熱力學(xué)平衡狀態(tài)下（恒溫Ｔ）：觀測(cè)量Ｈamilton量分布密度函數(shù)配分函數(shù)§4-3蒙特卡洛方法在統(tǒng)計(jì)物理中的應(yīng)用上述公式涉及的是高維積分問(wèn)題。只有理想氣體、諧振子系統(tǒng)、Ising模型等極少數(shù)類型的問(wèn)題可以解決求解。大多數(shù)情況下，只能借助近似方法求出。蒙特卡洛方法正則系綜，

除掉動(dòng)量以外的其它的相空間坐標(biāo)當(dāng)粒子間的相互作用與動(dòng)量無(wú)關(guān)時(shí)，動(dòng)量項(xiàng)的貢獻(xiàn)可以被積分，這相當(dāng)于將Hamilton量中的動(dòng)量項(xiàng)去掉。則，平衡態(tài)下的概率分布為Boltzmann分布?！?-3蒙特卡洛方法在統(tǒng)計(jì)物理中的應(yīng)用Boltzmann分布密度函數(shù)為：可見(jiàn)，所有對(duì)應(yīng)于大能量的狀態(tài)對(duì)應(yīng)觀測(cè)量積分的貢獻(xiàn)都很小。隨機(jī)選擇（均勻抽樣）ｎ各狀態(tài)，則有§4-3蒙特卡洛方法在統(tǒng)計(jì)物理中的應(yīng)用由于上式中大部分狀態(tài)對(duì)求和的貢獻(xiàn)很小，故抽樣的效率比較低。為有效的進(jìn)行計(jì)算，應(yīng)采用重要抽樣法。用Metropolis方法產(chǎn)生Boltzmann分布的ｎ個(gè)狀態(tài)，則湘潭大學(xué)物理學(xué)專業(yè)基礎(chǔ)課計(jì)算物理及其應(yīng)用

材料與光電物理學(xué)院變分量子蒙特卡洛方法求解簡(jiǎn)單諧振子的基態(tài)能量

1.引言2.變分量子蒙特卡洛方法

3.計(jì)算步驟

4.結(jié)果及分析

目錄1.引言經(jīng)典力學(xué)中，質(zhì)點(diǎn)在平衡位置附近的最基本最簡(jiǎn)單的運(yùn)動(dòng)是簡(jiǎn)諧振動(dòng)。在量子物理中與此對(duì)應(yīng)的微觀粒子的運(yùn)動(dòng)就是諧振子。簡(jiǎn)單諧振子的理論在應(yīng)用上有很大價(jià)值，因?yàn)榻?jīng)典力學(xué)告訴我們只要選擇適當(dāng)?shù)淖鴺?biāo)，任意粒子體系的微小振動(dòng)都可以認(rèn)為是一些相互獨(dú)立無(wú)關(guān)的振子的集合的運(yùn)動(dòng)。普朗克在他的輻射理論中將輻射物質(zhì)的中心當(dāng)作一些諧振子。從而得到和實(shí)驗(yàn)相符合的結(jié)果。在分子光譜中，可以把分子的振動(dòng)近似地當(dāng)作諧振子的波函數(shù)。另外在量子場(chǎng)論中電磁場(chǎng)的問(wèn)題也能歸結(jié)成諧振子的形式。因此，諧振子問(wèn)題的地位在物理學(xué)中非常重要。2.變分量子蒙特卡洛方法設(shè)量子體系的波函數(shù)為，則量子體系的Schr?dinger方程為簡(jiǎn)單諧振子體系的哈密頓量為V(rD)為諧振子勢(shì)，取其中，D表示的是維數(shù)，可取1、2和3維。取簡(jiǎn)單諧振子含調(diào)節(jié)參數(shù)的試探波函數(shù)形式2.變分量子蒙特卡洛方法使用量子單位，可得在計(jì)算過(guò)程中定義局部能量：波函數(shù)下的平均值為2.變分量子蒙特卡洛方法取N個(gè)含M個(gè)隨機(jī)狀態(tài)的系綜，則：標(biāo)準(zhǔn)方差：其中3.計(jì)算步驟

計(jì)算步驟：取試探波函數(shù)形式（在實(shí)際的模擬過(guò)程還取到了試探波函數(shù)的形式，用以來(lái)對(duì)比優(yōu)劣）；定義調(diào)節(jié)參數(shù)α的變化范圍，用循環(huán)的語(yǔ)句來(lái)控制變化；取N個(gè)系綜，第個(gè)系綜隨機(jī)取M個(gè)隨機(jī)數(shù)；N個(gè)系綜同步隨機(jī)演化，每個(gè)系綜變換按照Metropolis方法來(lái)產(chǎn)生馬爾柯夫鏈，判定接受的依據(jù)為3.計(jì)算步驟計(jì)算不同參數(shù)α?xí)r，EL平均值：和標(biāo)準(zhǔn)方差：4.結(jié)果及分析

計(jì)算結(jié)果：當(dāng)取一維諧振子基態(tài)試探波函數(shù)為時(shí)的結(jié)果如下：

圖2橫坐標(biāo)為試探波函數(shù)的調(diào)節(jié)參數(shù)，縱坐標(biāo)為變化調(diào)節(jié)參數(shù)后對(duì)應(yīng)能量的標(biāo)準(zhǔn)方差圖1橫坐標(biāo)為試探波函數(shù)的調(diào)節(jié)參數(shù)，縱坐標(biāo)為變化調(diào)節(jié)參數(shù)后對(duì)應(yīng)的能量值4.結(jié)果及分析當(dāng)取一維諧振子基態(tài)試探波函數(shù)為時(shí)的結(jié)果如下：

圖4橫坐標(biāo)為試探波函數(shù)的調(diào)節(jié)參數(shù)，縱坐標(biāo)為變化調(diào)節(jié)參數(shù)后對(duì)應(yīng)能量的標(biāo)準(zhǔn)方差圖3橫坐標(biāo)為試探波函數(shù)的調(diào)節(jié)參數(shù)，縱坐標(biāo)為變化調(diào)節(jié)參數(shù)后對(duì)應(yīng)的能量值4.結(jié)果及分析結(jié)果分析：諧振子系統(tǒng)具有精確的解析波函數(shù)，用變分量子蒙特卡洛方法計(jì)算出的結(jié)果與解析結(jié)果很好的吻合。用變分量子蒙特卡洛方法求解量子體系的基態(tài)能和基態(tài)波函數(shù)時(shí)，可以加入對(duì)局部能量的標(biāo)準(zhǔn)方差的統(tǒng)計(jì)，在結(jié)果中取標(biāo)準(zhǔn)方差最小的試探波函數(shù)形式，再由能量取最小來(lái)確定的α值，這樣，能夠準(zhǔn)確的從一系列試探波函數(shù)形式中選出最優(yōu)基態(tài)波函數(shù)，并同時(shí)確定比較精確的基態(tài)能。證實(shí)了變分量子蒙特卡洛方法方法在解決量子體系問(wèn)題方面的可行性和可靠性。統(tǒng)計(jì)力學(xué)：系綜時(shí)間平均統(tǒng)計(jì)力學(xué)：各態(tài)歷經(jīng)假設(shè)系綜平均統(tǒng)計(jì)力學(xué)：平均值隨機(jī)過(guò)程動(dòng)力學(xué)變量轉(zhuǎn)移概率分布函數(shù)歸一化條件Markov過(guò)程主方程Fokker-PlanckEquation布朗運(yùn)動(dòng)愛(ài)因斯坦模型概率演化的主方程差分形式微分形式平衡態(tài)簡(jiǎn)單抽樣方法BoltzmannStatisticalAverage重要抽樣方法BoltzmannStatisticalAverage經(jīng)典粒子系統(tǒng)經(jīng)典自旋系伊辛模型XY模型海森堡模型經(jīng)典自旋系量子蒙特卡羅方法湘潭大學(xué)物理學(xué)專業(yè)基礎(chǔ)課計(jì)算物理及其應(yīng)用

材料與光電物理學(xué)院蒙特卡洛方法的簡(jiǎn)單應(yīng)用

——求圓周率和模擬氫原子電子云

1.1蒙特卡洛處理的兩類問(wèn)題1.2求圓周率的數(shù)學(xué)模型1.3計(jì)算結(jié)果1.4氫原子電子的分布密度1.5模擬結(jié)果目錄確定性問(wèn)題1.1蒙特卡洛處理的兩類問(wèn)題隨機(jī)性問(wèn)題原子核物理問(wèn)題、運(yùn)籌學(xué)中的庫(kù)存問(wèn)題、隨機(jī)服務(wù)系統(tǒng)中的排隊(duì)問(wèn)題，動(dòng)物的生態(tài)競(jìng)爭(zhēng)和傳染病的蔓延計(jì)算多重積分、求逆矩陣、解線性代數(shù)方程組、解積分方程、解某些偏微分方程邊值問(wèn)題和計(jì)算微分算子的本征值基本思想：針對(duì)待求問(wèn)題，根據(jù)物理現(xiàn)象本身的統(tǒng)計(jì)規(guī)律，或人為構(gòu)造一合適的依賴隨機(jī)變量的概率模型，使某些隨機(jī)變量的統(tǒng)計(jì)量為待求問(wèn)題的解，進(jìn)行大統(tǒng)計(jì)量N→∞的統(tǒng)計(jì)實(shí)驗(yàn)方法或計(jì)算機(jī)隨機(jī)模擬方法。1.2求圓周率的數(shù)學(xué)模型y2rr2rOrx在正方形區(qū)域內(nèi)產(chǎn)生N個(gè)隨機(jī)點(diǎn){(x1,y1),(x2,y2),…,(xN,yN)}記錄落在圓內(nèi)的點(diǎn)數(shù)為M，如果取半徑r=1，那么圓周率為：1.3計(jì)算結(jié)果π值隨模擬次數(shù)的變化1.4氫原子電子的分布密度由原子物理理論和量子理論可知，氫原子態(tài)的波函數(shù)只是半徑的函數(shù)，與和無(wú)關(guān)。而氫原子中電子半徑的分布密度，即電子在半徑處單位厚度球殼內(nèi)出現(xiàn)的幾率習(xí)慣上把這種分布形象地稱為電子云。氫原子的基態(tài)即１s態(tài)（n=1,l=0,m=0），有，是D的最大值處的r值，其值與其中玻爾半徑相同。r0是D收斂處的r值，即D的收斂點(diǎn)。1.4氫原子電子的分布密度氫原子的2s態(tài)(n=2,l=0,m=0),有

氫原子的3s態(tài)(n=3,l=0,m=0),有

1.5模擬結(jié)果氫原子態(tài)電子云模擬圖氫原子態(tài)電子云模擬圖氫原子態(tài)電子云模擬圖氫原子組合電子云模擬圖附錄：程序具體FORTRAN程序如下：programmain

implicitnone

real(kind=8),parameter::m=1.0e8

integer::n,counter=0

real(kind=8)::x(1:m),y(1:m),Pi,r

open(unit=10,file='random_num.txt')

callrandom_seed()

CALLrandom_number(x)

callrandom_number(y)

don=1,m

r=sqrt(x(n)**2+y(n)**2)

if(r<=1)then

counter=counter+1

endif

enddo

pi=(counter*1.0)*4.0/(m*1.0)

writ