蒙特卡羅方法課件

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

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

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

理論依據(jù)：

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

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

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

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

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

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

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

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

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

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

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

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

用來作為隨機(jī)數(shù)發(fā)生器的物理源主要有兩種：一種是根據(jù)放射性物質(zhì)的放射性，另一種是利用計(jì)算機(jī)的固有噪聲。用物理方法產(chǎn)生的隨機(jī)數(shù)序列無法重復(fù)實(shí)現(xiàn)，不能進(jìn)行程序復(fù)算，給驗(yàn)證結(jié)果帶來很大困難。而且，需要增加隨機(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ù)為例，平方取中法是把一個2S位的十進(jìn)制自平方后，去頭截尾只保留中間2S個數(shù)字，然后用102S來除，這樣就可以得到在[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)，很少使用。§2偽隨機(jī)數(shù)的產(chǎn)生和檢驗(yàn)---數(shù)學(xué)方法由Lehmer在1951年提出來的，它的一般形式是：對于任一初始值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)等長度子區(qū)間中隨機(jī)數(shù)的數(shù)量是一樣的。按先后順序出現(xiàn)的隨機(jī)數(shù)中，每個隨機(jī)數(shù)的取值與其相距一定間隔的隨機(jī)數(shù)取值之間無關(guān)?！?偽隨機(jī)數(shù)序列的統(tǒng)計(jì)檢驗(yàn)判斷偽隨機(jī)數(shù)序列是否滿足均勻和相互獨(dú)立的要求，要靠統(tǒng)計(jì)檢驗(yàn)的方法實(shí)現(xiàn)。對于偽隨機(jī)數(shù)的統(tǒng)計(jì)檢驗(yàn)，一般包括兩大類：均勻性檢驗(yàn)和獨(dú)立性檢驗(yàn)偽隨機(jī)數(shù)的均勻性將區(qū)間[0,1]分為K個子區(qū)間，統(tǒng)計(jì)隨機(jī)數(shù)落在第k

個子區(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

?！?隨機(jī)數(shù)的產(chǎn)生和檢驗(yàn)---統(tǒng)計(jì)檢驗(yàn)K的取值不能太大也不能太小，太大反映不出“小區(qū)間”的均勻性，太小反映不出“大區(qū)間”的均勻性。限制條件（1）概率論中的Pearson定理說明，上式的極限概率分布是χ2分布給出了χ2≤x的概率。整數(shù)ν是系統(tǒng)自由度表示獨(dú)立測量的次數(shù)，由于存在一個限制條件，故ν=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ì)算出來的χ小于χα，則認(rèn)為在此置信度下，偽隨機(jī)數(shù)序列在[0,1]中是均勻分布的。（1）順序相關(guān)法它用相鄰兩個隨機(jī)數(shù)的自相關(guān)函數(shù)（或相關(guān)系數(shù)）來標(biāo)識偽隨機(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個隨機(jī)數(shù)作為S維空間中的一個點(diǎn)的坐標(biāo)值，于是可以構(gòu)成一個點(diǎn)序列。（2）把S維空間中的單位方體分成為K個子方體，方體邊長（3）統(tǒng)計(jì)落在第k個子方體中的實(shí)際頻數(shù)nk，它應(yīng)當(dāng)趨近于理論頻數(shù)：例如將2N

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

,

ξ

3,…}和{ξ

2,ξ

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

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

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

，則§2隨機(jī)數(shù)產(chǎn)生及檢驗(yàn)---獨(dú)立性統(tǒng)計(jì)檢驗(yàn)用連續(xù)的兩個隨機(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個值x1,x2,x3，它們出現(xiàn)的幾率分別為2/8,5/8,1/8，則隨機(jī)數(shù)小于2/8時實(shí)現(xiàn)x1

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

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

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

σ2

σ3

反應(yīng)幾率f：ε1=σ1/σε2=σ2/σε3=σ3/σ歸一化：例如Poisson分布是離散型分布對此分布進(jìn)行抽樣得到第n個事件發(fā)生的條件為電子對產(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ù)。這對一些簡單的幾率密度函數(shù)解析表達(dá)式是很容易做到的。如粒子隨機(jī)運(yùn)動的自由程分布為指數(shù)分布：則求反函數(shù)后得§3隨機(jī)變量的抽樣---連續(xù)型累積函數(shù)注意（1-ξ）和ξ同樣服從[0，1]的均勻分布§3隨機(jī)變量的抽樣---變換抽樣法一維：變換抽樣法的基本思想是將一個比較復(fù)雜的分布p(x)的抽樣，變換為已知的簡單分布g(y)的抽樣我們希望找到x?y

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

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

。二維：有兩個變量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)|

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

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

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

或y都滿足正態(tài)分布?？梢?，為了得到滿足一個復(fù)雜分布的隨機(jī)抽樣，這里用了兩個滿足簡單分布的隨機(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足夠大時，In就是積分I的一個無偏估計(jì)值。相應(yīng)的方差為：可見，當(dāng)f(x)在其定義域內(nèi)變化較大時，方差較大?！?-1蒙特卡洛方法在積分計(jì)算中的應(yīng)用重要抽樣法：當(dāng)f(x)在其定義域內(nèi)有顯著的起伏變化時，可采用重要抽樣法?？蓪⒎e分稍作變換：適當(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í)際物理問題中，被積函數(shù)在超立方體區(qū)域內(nèi)可能強(qiáng)烈變化。若在積分區(qū)域內(nèi)均勻抽樣，積分貢獻(xiàn)可能主要來自少數(shù)僅僅只有幾個蒙特卡洛投點(diǎn)的小區(qū)域，從而導(dǎo)致很大的統(tǒng)計(jì)誤差。所以采用重要抽樣法，使得隨機(jī)點(diǎn)更多地投在取值大的區(qū)間§4-1蒙特卡洛方法在積分計(jì)算中的應(yīng)用選取偏倚分布密度函數(shù)，并定義（要求方差較小）有：按照偏倚分布密度函數(shù)在區(qū)域抽樣N個子樣則：積分的近似值§4-1蒙特卡洛方法在積分計(jì)算中的應(yīng)用一維積分的擲點(diǎn)法一維標(biāo)準(zhǔn)積分：定義：則有：在單位正方形內(nèi)投N個點(diǎn)，落在曲線f(x)下的有M個，則由于對y的積分可以解析計(jì)算，故此法的誤差較平均值法大?！?-2蒙特卡洛方法在量子力學(xué)中的應(yīng)用量子力學(xué)回顧薛定諤方程定義：費(fèi)曼傳播子量子力學(xué)的基本理論告訴我們：系統(tǒng)的所有信息，包括基態(tài)、激發(fā)態(tài)的能量、波函數(shù)等論可由費(fèi)曼傳播子給出。特別是和基態(tài)有關(guān)的信息，可以很方便的得到。§4-2蒙特卡洛方法在量子力學(xué)中的應(yīng)用例如：基態(tài)波函數(shù)的模方可表示為費(fèi)曼傳播子可以表示為路徑積分的形式，對于簡單系統(tǒng)，即有其中，作用量§4-2蒙特卡洛方法在量子力學(xué)中的應(yīng)用路徑積分量子蒙特卡洛方法實(shí)際數(shù)值計(jì)算中，費(fèi)曼傳播子表達(dá)式為

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

格林函數(shù)歸一，且與ｘ，ｔ無關(guān)。擴(kuò)散方程馬爾科夫過程（格林函數(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)解虛時薛定諤方程（?。?/p>

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

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

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

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

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

3.計(jì)算步驟

4.結(jié)果及分析

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

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

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

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

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

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

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

1.1蒙特卡洛處理的兩類問題1.2求圓周率的數(shù)學(xué)模型1.3計(jì)算結(jié)果1.4氫原子電子的分布密度1.5模擬結(jié)果目錄確定性問題1.1蒙特卡洛處理的兩類問題隨機(jī)性問題原子核物理問題、運(yùn)籌學(xué)中的庫存問題、隨機(jī)服務(wù)系統(tǒng)中的排隊(duì)問題，動物的生態(tài)競爭和傳染病的蔓延計(jì)算多重積分、求逆矩陣、解線性代數(shù)方程組、解積分方程、解某些偏微分方程邊值問題和計(jì)算微分算子的本征值基本思想：針對待求問題，根據(jù)物理現(xiàn)象本身的統(tǒng)計(jì)規(guī)律，或人為構(gòu)造一合適的依賴隨機(jī)變量的概率模型，使某些隨機(jī)變量的統(tǒng)計(jì)量為待求問題的解，進(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個隨機(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ù)，與和無關(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