計(jì)算物理MC_first

1、 Monte CarloStatistic of random variables Metropolis algorithm Implementation of the Metropolis MC method2對統(tǒng)計(jì)力學(xué)體系進(jìn)行計(jì)算機(jī)模擬時，需要確定對統(tǒng)計(jì)力學(xué)體系進(jìn)行計(jì)算機(jī)模擬時，需要確定體系的位形（組態(tài)）。按照產(chǎn)生位形變化的方體系的位形（組態(tài)）。按照產(chǎn)生位形變化的方法，可以將計(jì)算機(jī)模擬分成兩大類：法，可以將計(jì)算機(jī)模擬分成兩大類：一類是一類是( ( stochastic) )方法：方法：Monte Carlo: MC; Simulated annealing 一類是一類是( (determi

2、nistic方法：方法： 按照體系的動力學(xué)規(guī)律產(chǎn)生位形變化。按照體系的動力學(xué)規(guī)律產(chǎn)生位形變化。Molecular Dynamics: MDMonte Carlo方法（方法（MC）亦稱為隨機(jī)模擬（）亦稱為隨機(jī)模擬（Random Simulation）、隨機(jī)抽樣）、隨機(jī)抽樣(Random Sampling) 或統(tǒng)計(jì)試或統(tǒng)計(jì)試驗(yàn)（驗(yàn)（Statistical Testing）方法。）方法。Monte Carlo 方法名稱的固定、得到系統(tǒng)發(fā)展大方法名稱的固定、得到系統(tǒng)發(fā)展大約是二十世紀(jì)四十年代；但如就其方法特征可追約是二十世紀(jì)四十年代；但如就其方法特征可追溯到十八世紀(jì)著名的溯到十八世紀(jì)著名的Buffon

3、問題問題-用隨機(jī)投擲用隨機(jī)投擲縫針試驗(yàn)計(jì)算圓周率數(shù)值?？p針試驗(yàn)計(jì)算圓周率數(shù)值。The Comte de Buffon needle experiment, AD 1777dLp2dLNMdLp2MdNL2試驗(yàn)方案是：在平面內(nèi)劃一組相距試驗(yàn)方案是：在平面內(nèi)劃一組相距 d 的平行線，向此的平行線，向此平面隨意的投擲長度平面隨意的投擲長度 L 的細(xì)針，那末從針與平行線相的細(xì)針，那末從針與平行線相交的概率交的概率 p 可以得到可以得到 值。值。the Comte de Buffon needle experiment實(shí)驗(yàn)者實(shí)驗(yàn)者時間時間年份年份針長針長投針投針次數(shù)次數(shù)相交相交次數(shù)次數(shù)的的估計(jì)估計(jì)Wol

4、fWolf185018500.80.850005000253225323.15963.1596SmithSmith185518550.60.6320432041218.51218.53.15543.1554De Morgan CDe Morgan C186018601.01.0600600382.5382.53.1373.137FoxFox188418840.750.75103010304894893.15953.1595LezzeriniLezzerini190119010.830.8334083408180818083.14159293.1415929Reina19250.54192520

5、8593.1795歷史上一些有名用投針試驗(yàn)計(jì)算歷史上一些有名用投針試驗(yàn)計(jì)算值的結(jié)果值的結(jié)果Monte CarloMonte Carlo方法以概率統(tǒng)計(jì)為理論基礎(chǔ)，方法以概率統(tǒng)計(jì)為理論基礎(chǔ)，以隨機(jī)抽樣為主要手段。其基本思想是首先以隨機(jī)抽樣為主要手段。其基本思想是首先建立一個概率（或隨機(jī)過程）模型，使它的建立一個概率（或隨機(jī)過程）模型，使它的參數(shù)等于問題的解；然后通過對模型（或過參數(shù)等于問題的解；然后通過對模型（或過程）的抽樣試驗(yàn)來獲取有關(guān)參數(shù)的統(tǒng)計(jì)特征程）的抽樣試驗(yàn)來獲取有關(guān)參數(shù)的統(tǒng)計(jì)特征、解的近似值及精度估計(jì)。、解的近似值及精度估計(jì)。Monte Carlo的基本思想的基本思想NnnNN11設(shè)所要

6、求的量設(shè)所要求的量x是隨機(jī)變量是隨機(jī)變量的數(shù)學(xué)期望的數(shù)學(xué)期望E()，那么，那么用用Monte Carlo方法來近似確定方法來近似確定x的方法是對的方法是對進(jìn)行進(jìn)行N此重復(fù)抽樣，產(chǎn)生相互獨(dú)立的此重復(fù)抽樣，產(chǎn)生相互獨(dú)立的值的序列值的序列1，2，,N并計(jì)算其算術(shù)平均值：并計(jì)算其算術(shù)平均值：根據(jù)根據(jù)Kolmogorov的大數(shù)定理則有的大數(shù)定理則有即當(dāng)即當(dāng)N充分大時，充分大時， 成立的概率等于成立的概率等于1.亦即可以用亦即可以用 作為所求量作為所求量x的估計(jì)值。的估計(jì)值。1)lim(xPNNxEN)(NMeasuring ErrorN0Variance 方差方差: 2 = Var(X) = - 2St

7、andard deviation (standard error) of a random number is the square root of its variance2 = Var(X) = - 2Variance 實(shí)驗(yàn)者實(shí)驗(yàn)者時間時間年份年份針長針長投針投針次數(shù)次數(shù)相交相交次數(shù)次數(shù)的的估計(jì)估計(jì)WolfWolf185018500.80.850005000253225323.15963.1596SmithSmith185518550.60.6320432041218.51218.53.15543.1554De Morgan CDe Morgan C186018601.01.0600600

8、382.5382.53.1373.137FoxFox188418840.750.75103010304894893.15953.1595LezzeriniLezzerini190119010.830.8334083408180818083.14159293.1415929Reina19250.541925208593.1795歷史上一些有名用投針試驗(yàn)計(jì)算歷史上一些有名用投針試驗(yàn)計(jì)算值的結(jié)果值的結(jié)果N/1顯然這一早期的顯然這一早期的“古典古典Monte Carlo方法方法”及相應(yīng)的及相應(yīng)的抽樣實(shí)踐已寓示與這種數(shù)值計(jì)算方法相伴隨的巨大抽樣實(shí)踐已寓示與這種數(shù)值計(jì)算方法相伴隨的巨大工作量。所以可以理解

9、直到近半個世紀(jì)以來，工作量。所以可以理解直到近半個世紀(jì)以來，Monte Carlo方法的應(yīng)用范圍才不斷擴(kuò)展，形成計(jì)算方法的應(yīng)用范圍才不斷擴(kuò)展，形成計(jì)算數(shù)學(xué)的一個重要分支，是與電子計(jì)算機(jī)的相應(yīng)發(fā)展數(shù)學(xué)的一個重要分支，是與電子計(jì)算機(jī)的相應(yīng)發(fā)展不可分割的。不可分割的。 如果使如果使 值的精度達(dá)三位有效數(shù)字，則需數(shù)十萬次。值的精度達(dá)三位有效數(shù)字，則需數(shù)十萬次。實(shí)驗(yàn)者實(shí)驗(yàn)者時間時間年份年份針長針長投針投針次數(shù)次數(shù)相交相交次數(shù)次數(shù)的的估計(jì)估計(jì)WolfWolf185018500.80.850005000253225323.15963.1596SmithSmith185518550.60.6320432041

10、218.51218.53.15543.1554De Morgan CDe Morgan C186018601.01.0600600382.5382.53.1373.137FoxFox188418840.750.75103010304894893.15953.1595LezzeriniLezzerini190119010.830.8334083408180818083.14159293.1415929Reina19250.541925208593.1795歷史上一些有名用投針試驗(yàn)計(jì)算歷史上一些有名用投針試驗(yàn)計(jì)算值的結(jié)果值的結(jié)果N/1Monte Carlo模擬就是邊產(chǎn)生隨機(jī)數(shù)，邊在計(jì)算機(jī)上進(jìn)行隨模

11、擬就是邊產(chǎn)生隨機(jī)數(shù)，邊在計(jì)算機(jī)上進(jìn)行隨機(jī)過程模擬。機(jī)過程模擬。步驟如下步驟如下運(yùn)行運(yùn)行（1）在計(jì)算機(jī)上產(chǎn)生隨機(jī)數(shù)，并使隨機(jī)數(shù)）在計(jì)算機(jī)上產(chǎn)生隨機(jī)數(shù)，并使隨機(jī)數(shù)游動。游動。（2）進(jìn)行條件判斷，在滿足條件的情況下）進(jìn)行條件判斷，在滿足條件的情況下，按規(guī)則進(jìn)行作業(yè)。，按規(guī)則進(jìn)行作業(yè)。反復(fù)操作（反復(fù)操作（1）、（）、（2）。）。The Name of the GameMetropolis coined the name “Monte Carlo”, from its gambling Casino.Monte-Carlo, MonacoMonte Carlo 是位于歐洲地中海沿岸摩納哥是位于歐洲地中海

12、沿岸摩納哥（Monaco）國的一個重要城市。）國的一個重要城市。Stanislaw Ulam (1909-1984)S. Ulam is credited as the inventor of Monte Carlo method in 1940s, which solves mathematical problems using statistical sampling.大數(shù)定理是大數(shù)定理是Monte CaloMonte Calo模擬的理論基礎(chǔ)，模擬的理論基礎(chǔ)，我先做一下簡單的介紹我先做一下簡單的介紹大數(shù)定理大數(shù)定理 大數(shù)定律大數(shù)定律大數(shù)法則大數(shù)法則概率論歷史上第一個極限定理屬于伯努利，概率

13、論歷史上第一個極限定理屬于伯努利，后人稱之為后人稱之為“大數(shù)定律大數(shù)定律”。是概率論中討論。是概率論中討論隨機(jī)變量序列的隨機(jī)變量序列的算術(shù)平均值向常數(shù)收斂算術(shù)平均值向常數(shù)收斂的定的定律。是概率論與數(shù)理統(tǒng)計(jì)學(xué)的基本定律之一。律。是概率論與數(shù)理統(tǒng)計(jì)學(xué)的基本定律之一。又稱弱大數(shù)理論。又稱弱大數(shù)理論。數(shù)學(xué)家伯努利數(shù)學(xué)家伯努利 主要含義主要含義 有些隨機(jī)事件無規(guī)律可循，但不少卻是有規(guī)律有些隨機(jī)事件無規(guī)律可循，但不少卻是有規(guī)律的，這些的，這些“有規(guī)律的隨機(jī)事件有規(guī)律的隨機(jī)事件”中在大量重復(fù)中在大量重復(fù)出現(xiàn)的條件下，往往呈現(xiàn)幾乎必然的統(tǒng)計(jì)特性，出現(xiàn)的條件下，往往呈現(xiàn)幾乎必然的統(tǒng)計(jì)特性，這個規(guī)律就是大數(shù)定律。

14、這個規(guī)律就是大數(shù)定律。通俗地說，這個定理就是，在試驗(yàn)不變的條件通俗地說，這個定理就是，在試驗(yàn)不變的條件下，重復(fù)試驗(yàn)多次，隨機(jī)事件的下，重復(fù)試驗(yàn)多次，隨機(jī)事件的頻率頻率近似于它近似于它的的概率概率。比如，我們向上拋一枚硬幣，硬幣落下后哪一面朝上本比如，我們向上拋一枚硬幣，硬幣落下后哪一面朝上本來是偶然的，但當(dāng)我們上拋硬幣的次數(shù)足夠多后，達(dá)到來是偶然的，但當(dāng)我們上拋硬幣的次數(shù)足夠多后，達(dá)到上萬次甚至幾十萬幾百萬次以后，我們就會發(fā)現(xiàn)，硬幣上萬次甚至幾十萬幾百萬次以后，我們就會發(fā)現(xiàn)，硬幣每一面向上的次數(shù)約占總次數(shù)的二分之一。這種情況下，每一面向上的次數(shù)約占總次數(shù)的二分之一。這種情況下，偶然中包含著必然

15、。必然的規(guī)律與特性在大量的樣本中偶然中包含著必然。必然的規(guī)律與特性在大量的樣本中得以體現(xiàn)。得以體現(xiàn)。 又如稱量某一物體的重量，假如衡器不存在系統(tǒng)偏差，又如稱量某一物體的重量，假如衡器不存在系統(tǒng)偏差，由于衡器的精度等各種因素的影響，對同一由于衡器的精度等各種因素的影響，對同一物體重復(fù)稱物體重復(fù)稱量多次量多次，可能得到多個不同的重量數(shù)值，但它們的算術(shù)，可能得到多個不同的重量數(shù)值，但它們的算術(shù)平均值一般來說將隨稱量次數(shù)的增加而逐漸接近于物體平均值一般來說將隨稱量次數(shù)的增加而逐漸接近于物體的真實(shí)重量。的真實(shí)重量。 伯努利大數(shù)定律：伯努利大數(shù)定律：設(shè)設(shè) n 為為 n 重伯努利實(shí)驗(yàn)中事件重伯努利實(shí)驗(yàn)中事件

16、 A 發(fā)生的次數(shù)，發(fā)生的次數(shù)，p 為每次實(shí)驗(yàn)中為每次實(shí)驗(yàn)中 A 出現(xiàn)的概率，則對任意的出現(xiàn)的概率，則對任意的 0 有（有（2）成立。）成立。 limn1pnPn(2)這一定理指出，不論隨機(jī)變量的分布如何，只這一定理指出，不論隨機(jī)變量的分布如何，只要要n n足夠大，事件足夠大，事件A A出現(xiàn)的頻率趨向一個穩(wěn)定的出現(xiàn)的頻率趨向一個穩(wěn)定的值，這個值就是值，這個值就是A A的概率。的概率。What is a Random Number? Follow a definite distribution, usually uniform distribution Uncorrelated Unpredict

17、ableThe Roulette and DiceMechanical random number generators運(yùn)用運(yùn)用Monte Carlo方法首先要產(chǎn)生一組隨機(jī)數(shù)列。隨機(jī)數(shù)列方法首先要產(chǎn)生一組隨機(jī)數(shù)列。隨機(jī)數(shù)列的產(chǎn)生方法，一種是找一個實(shí)際的隨機(jī)物理過程，譬如說的產(chǎn)生方法，一種是找一個實(shí)際的隨機(jī)物理過程，譬如說擲骰子，紀(jì)錄放射輻射源中兩個衰變之間的時間值等等。擲骰子，紀(jì)錄放射輻射源中兩個衰變之間的時間值等等。另一種方法是由計(jì)算機(jī)自己產(chǎn)生隨機(jī)數(shù)列，由于這些隨機(jī)另一種方法是由計(jì)算機(jī)自己產(chǎn)生隨機(jī)數(shù)列，由于這些隨機(jī)數(shù)不是從實(shí)際過程中得來的，故稱之為數(shù)不是從實(shí)際過程中得來的，故稱之為“偽隨機(jī)偽

18、隨機(jī)數(shù)數(shù)” (Pseudo-random numbers) 。產(chǎn)生隨機(jī)數(shù)列的程序稱。產(chǎn)生隨機(jī)數(shù)列的程序稱為隨機(jī)數(shù)產(chǎn)生器為隨機(jī)數(shù)產(chǎn)生器(random number generator)Pseudo-Random Numbers & random number generator偽隨機(jī)偽隨機(jī)數(shù)數(shù) Truly random numbers can not be generated on a computer Pseudo-random numbers follows a well-defined algorithm, thus predictable and repeatable Pseudo-ra

19、ndom numbers have nearly all the properties of true random numbersPseudo-Random Numbers隨機(jī)數(shù)可以有各種分布形式，它們都可以從隨機(jī)數(shù)可以有各種分布形式，它們都可以從均勻分布均勻分布隨隨機(jī)數(shù)列中得出。而產(chǎn)生均勻分布隨機(jī)數(shù)列的數(shù)學(xué)方法又有機(jī)數(shù)列中得出。而產(chǎn)生均勻分布隨機(jī)數(shù)列的數(shù)學(xué)方法又有很多種，下面介紹幾種。很多種，下面介紹幾種。隨機(jī)性和統(tǒng)計(jì)獨(dú)立性要好隨機(jī)性和統(tǒng)計(jì)獨(dú)立性要好容易在計(jì)算機(jī)上實(shí)現(xiàn)容易在計(jì)算機(jī)上實(shí)現(xiàn)效率高效率高省時、占內(nèi)存小省時、占內(nèi)存小周期足夠大周期足夠大（1）平方取中）平方取中給定第一個數(shù)以后，則序

20、列中每個數(shù)可由前一個數(shù)按上給定第一個數(shù)以后，則序列中每個數(shù)可由前一個數(shù)按上式求得。式求得。0.6406, 0.0368, 0.135464061x4103683621xSnnx21110/把把2S位的數(shù)自乘，去頭截尾只保留中間的位的數(shù)自乘，去頭截尾只保留中間的2S位。位。13543x03682x0013542422x0183331623x)10),10/mod(int(221SSnnxx135403686406321xxx018333160013542441036836232221xxxSnnx21110/mod (a, p)a - INT(a / p) * pmod (7,3)3681000

21、0*)0368.41int(410368)10000,410368mod()10),10/mod(int(221SSnnxx)2),2/mod(int(221SSnnxxSnnx2112/平方取中平方取中)10),10/mod(int(221SSnnxxSnnx21110/)2),2/mod(int(212SSnnnxxxSnnx2222/乘積取中乘積取中二進(jìn)位二進(jìn)位)10),10/mod(int(221SSnnxxSnnx21110/)10),10/mod(int(212SSnnnxxxSnnx22210/乘積取中乘積取中平方取中平方取中十進(jìn)位十進(jìn)位（2）乘同余法）乘同余法(linear c

22、ongruential method)的疊代的疊代公式為公式為),mod(1cbaxxnn給定第一個數(shù)以后，則序列中每個數(shù)可由前一個數(shù)按上給定第一個數(shù)以后，則序列中每個數(shù)可由前一個數(shù)按上式求得。式求得。其中其中a, b and c are magic numbers: the choice of these numbers determines the quality of the generator.214748364712, 0,168077315cbacxranfn/()1MOD (a, p)a - INT(a / p) * pChoice of ParametersNameca (mu

23、ltiplier)bperiodANSIC rand()231110351524512345231Park-MillerNR ran0()231-1168070231-2drand48() 2482521490391711248Hayes 64-bit26463641362238467930051264(a x + b) mod cShort-Coming of LCGxnWhen (xn, xn+1) pairs are plotted for all n, a lattice structure is shown.xn+1),mod(1cbaxxnn1,1282, 0, 507xcbacx

24、ranfn/()1iiryrx212nnxx168071cxranfn/()1)2,2(313121474836471231c111cxxnn1*20iseedxIA/ICIQ )IA,IC(MODIR 12, 1 Iseed31Although popular, by virtue of the ease with which it can be programmed the linear congruential method does not satisfy all of the requirements that are now regarded as important in a r

25、andom number generator.隨機(jī)性和統(tǒng)計(jì)獨(dú)立性要好隨機(jī)性和統(tǒng)計(jì)獨(dú)立性要好容易在計(jì)算機(jī)上實(shí)現(xiàn)容易在計(jì)算機(jī)上實(shí)現(xiàn)效率高效率高省時、占內(nèi)存小省時、占內(nèi)存小周期足夠大周期足夠大mnnryrxxnxn+1（4） The Marsaglia random number generatorCombination generator：A period 2144.I) Lagged Fibonacci generator1yxyxyxyx snrnnxxx ii) Arithmetic sequence method16777216/16777213dcdcdcdc 16777216/7654

26、321cc1nn nnncxU Numerical integrationMonte Carlo integration1 定積分問題定積分問題 。dxxI10)exp(這個積分的值這個積分的值 I = - (e-1-1)0.63212為了用隨機(jī)抽樣方法來求解該積分為了用隨機(jī)抽樣方法來求解該積分, 我們先要構(gòu)造一我們先要構(gòu)造一個概率模型個概率模型. 設(shè)所要求的量設(shè)所要求的量x是隨機(jī)變量是隨機(jī)變量的數(shù)學(xué)期望的數(shù)學(xué)期望E()，那么用，那么用Monte Carlo方法來近似確定方法來近似確定x的方法是的方法是對對進(jìn)行進(jìn)行N次重復(fù)抽樣，產(chǎn)生相互獨(dú)立的次重復(fù)抽樣，產(chǎn)生相互獨(dú)立的值的序列值的序列1，2，,

27、N并計(jì)算其算術(shù)平均值：并計(jì)算其算術(shù)平均值：NnnNN11dxxI10)exp(圖中方框的總面積為圖中方框的總面積為1，而我們所求的積分值，而我們所求的積分值I即為圖中即為圖中陰影部分的面積。在正方形平面中均勻隨機(jī)的投點(diǎn)，則陰影部分的面積。在正方形平面中均勻隨機(jī)的投點(diǎn)，則落入陰影區(qū)的概率則為積分值落入陰影區(qū)的概率則為積分值I。NvI 擲點(diǎn)法擲點(diǎn)法Hit or miss method)exp( x00IENvI)(因此對于這里考慮的問題，我們擬構(gòu)造這樣的概率模型，即因此對于這里考慮的問題，我們擬構(gòu)造這樣的概率模型，即在在11的正方形平面中均勻隨機(jī)的投點(diǎn)，則落入陰影區(qū)的概率的正方形平面中均勻隨機(jī)的投

28、點(diǎn)，則落入陰影區(qū)的概率則為積分值則為積分值I。設(shè)在。設(shè)在N次投點(diǎn)試驗(yàn)中，落入陰影區(qū)的點(diǎn)為次投點(diǎn)試驗(yàn)中，落入陰影區(qū)的點(diǎn)為v次次,那那么觀察頻數(shù)么觀察頻數(shù)v也是一個隨機(jī)變量，其數(shù)學(xué)期望也是一個隨機(jī)變量，其數(shù)學(xué)期望 E(v)=NI。令。令表示觀察頻率，那么按照大數(shù)定理，當(dāng)表示觀察頻率，那么按照大數(shù)定理，當(dāng)N充分大時，頻率收斂充分大時，頻率收斂于概率，即以概率于概率，即以概率1成立。因此可由上述概率模型在成立。因此可由上述概率模型在N很大時很大時所得到的所得到的v/N等于所求的積分值等于所求的積分值I。圖中方框的總面積為圖中方框的總面積為1，而我們所求的積分值即為而我們所求的積分值即為圖中陰影部分的面

29、積。圖中陰影部分的面積。（1） 產(chǎn)生在產(chǎn)生在0，1區(qū)間上均勻分布獨(dú)立的隨機(jī)數(shù)區(qū)間上均勻分布獨(dú)立的隨機(jī)數(shù)r, r;NvI 計(jì)算機(jī)上的具體計(jì)算步驟為：計(jì)算機(jī)上的具體計(jì)算步驟為：（2） 令令r 和和r分別為所投點(diǎn)的分別為所投點(diǎn)的x, y坐標(biāo)值，若坐標(biāo)值，若re-r,則則表示所投點(diǎn)落在陰影區(qū)內(nèi)，因此表示所投點(diǎn)落在陰影區(qū)內(nèi)，因此v加上加上1, N也加上也加上1;若若re-r,則表示所投點(diǎn)落在陰影區(qū)外，因此則表示所投點(diǎn)落在陰影區(qū)外，因此v加上加上0,N加上加上1（3） 重復(fù)（重復(fù)（1），（），（2）直至）直至N足夠大足夠大.（4） 計(jì)算計(jì)算 以上稱為一維定積分計(jì)算的擲點(diǎn)法。以上稱為一維定積分計(jì)算的擲點(diǎn)法。

30、Nf1只要選取足夠多的隨機(jī)點(diǎn)，即當(dāng)只要選取足夠多的隨機(jī)點(diǎn)，即當(dāng)N充分大時充分大時)(iif)(NiifNI11（2）一維定積分計(jì)算的平均值法）一維定積分計(jì)算的平均值法若在若在x的定義域的定義域0,1上均勻隨機(jī)的取點(diǎn)，即上均勻隨機(jī)的取點(diǎn)，即選取隨機(jī)數(shù)選取隨機(jī)數(shù)，定義一個隨機(jī)變量，定義一個隨機(jī)變量1為為以概率以概率1成立。因此可由上述概率模型在成立。因此可由上述概率模型在N很大時得到所很大時得到所求的積分值求的積分值I。平均值法平均值法: Sample mean method)(NiifNI11)(11NiifNIiefi)(1) 產(chǎn)產(chǎn) 生在生在0，1區(qū)間上均勻分布獨(dú)立的隨機(jī)數(shù)區(qū)間上均勻分布獨(dú)立的

31、隨機(jī)數(shù)i。2) 令，令，, N加上加上1。3) 重復(fù)（重復(fù)（1），（），（2）直至）直至N足夠大足夠大.4) 計(jì)算計(jì)算一維定積分計(jì)算的平均值法具體計(jì)算步驟為一維定積分計(jì)算的平均值法具體計(jì)算步驟為dxxA1021Nf141412rMNA如果一維積分中的長度不是如果一維積分中的長度不是1，而是一般的情況（，而是一般的情況（a, b)， 那么，那么，f 在在(a, b)中中N個隨機(jī)點(diǎn)上的平均值是對個隨機(jī)點(diǎn)上的平均值是對dxxfba)(dyyfab10)()(作變換：作變換：x=a+(b-a)ydx=(b-a)dydxxf50)(NiiyfN1)(5dxxfba)(dyyfab10)()(作變換：作變

32、換：x=a+(b-a)yx=5yNiiyfNdxxf150)(5)(隨機(jī)點(diǎn)在區(qū)域（隨機(jī)點(diǎn)在區(qū)域（0, 1)中是均勻分布的。中是均勻分布的。dxxf52)(隨機(jī)點(diǎn)在區(qū)域（隨機(jī)點(diǎn)在區(qū)域（0,1)中是均勻分布的。中是均勻分布的。NiiyfN1)(3. 23 yxNf1x=5yNumerical Integration of Multi-dimensions 6152),(dxdyyxf NiiibafNdxdyyxf16152),(15),(AfdA用用Monte Carlo 方法計(jì)算積分方法計(jì)算積分, 最簡單的辦法最簡單的辦法就是在積分區(qū)域內(nèi)隨機(jī)地取一系列點(diǎn)就是在積分區(qū)域內(nèi)隨機(jī)地取一系列點(diǎn), 計(jì)算

33、被計(jì)算被積函數(shù)在這些點(diǎn)上的數(shù)值并取平均積函數(shù)在這些點(diǎn)上的數(shù)值并取平均, 然后乘以然后乘以積分區(qū)域的積分區(qū)域的. 這種方法稱為簡單抽樣法。這種方法稱為簡單抽樣法。2/11 NNfSuppose we need to integrate from x0 to x1. We shall subdivide this interval into n steps of size x=(x1x0)/n as shown in figurex)f(xf(x)dxI(f)xxx000 xxxf(x)dxI(f)00 xxxdxxx)(xfxx)(xf)f(x00)( 21)(200000302006121x)

34、(xfx)(xfx)f(xError Nx/ 12N2/11 NNf1N Trapezium rulexxf(x)f(xf(x)dxI(f)xxx)210000 xxxxxxdxxx)(xf)f(xf(x)dxI(f)0000)(000)()21) 6121212130020000 xOxxf(x)f(xxx)(xfx)(xf)f(x)f(x（Error 復(fù)化梯型公式復(fù)化梯型公式 Compound Trapezium rulexxf(x)f(xf(x)dxI(f)xxx)210000The Compound Trapezium Rule approximation to the integra

35、l is thereforeWhile the error for each step is O(x3), the cumulative error is N times this or O(x2) O(N-2)2N Simpsons rulex0 x0+ xx0+2 x)()(4)(6bfcfafabS4NNumerical Integration of Multi-dimensionsIn general if the error goes as O(Na) in one dimension, then the error in d dimensions goes as O(Na/d).多

36、重積分的計(jì)算多重積分的計(jì)算物理量的平均值的計(jì)算歸結(jié)為一個多重積分的計(jì)算物理量的平均值的計(jì)算歸結(jié)為一個多重積分的計(jì)算, 設(shè)系統(tǒng)由設(shè)系統(tǒng)由100個粒子構(gòu)成個粒子構(gòu)成, 每個粒子有每個粒子有6個自由度個自由度, 所所以需要計(jì)算以需要計(jì)算300重積分重積分. 現(xiàn)在考慮在每一維取現(xiàn)在考慮在每一維取10個點(diǎn)個點(diǎn), 總共有總共有10300個點(diǎn)個點(diǎn). 假設(shè)計(jì)算機(jī)每秒可以計(jì)算假設(shè)計(jì)算機(jī)每秒可以計(jì)算1012個點(diǎn)個點(diǎn), 計(jì)算計(jì)算 這個積分需要這個積分需要10288秒秒!問題比這個更嚴(yán)重問題比這個更嚴(yán)重! 如果如上述方式取點(diǎn)如果如上述方式取點(diǎn), 則積分區(qū)則積分區(qū)域的內(nèi)點(diǎn)數(shù)為域的內(nèi)點(diǎn)數(shù)為 8300, 在總的點(diǎn)中所占比例

37、為在總的點(diǎn)中所占比例為 (8/10)300300, 也就是說也就是說, 取的點(diǎn)基本上都在表面上取的點(diǎn)基本上都在表面上!定積分問題定積分問題Monte Carlo一維定積分計(jì)算一維定積分計(jì)算 X 多重定積分計(jì)算多重定積分計(jì)算 goodRandom variables 隨機(jī)變量隨機(jī)變量Probability distribution functions (PDF).Expectation ValueStatistic of a Random VariableStatistic of random variables0,1,1,0,0,0,1,0,0,1,1,0,1Random variables

38、隨機(jī)變量隨機(jī)變量Random variables 隨機(jī)變量隨機(jī)變量532463451）在相同的確定實(shí)驗(yàn)條件下，對的觀測無法給）在相同的確定實(shí)驗(yàn)條件下，對的觀測無法給出單一固定值出單一固定值; ;）必須依據(jù)遍舉測量原則，對所有可能取值給出）必須依據(jù)遍舉測量原則，對所有可能取值給出發(fā)生概率發(fā)生概率Random variables 隨機(jī)變量隨機(jī)變量 X Random variables are hence characterized by a domain which contains all possible valuesthat the random value may take. This d

39、omain has a corresponding probability distribution functions (PDF).Random variables 隨機(jī)變量隨機(jī)變量Discrete variable 離散變量舉例離散變量舉例: :DomainTo this domain we have the corresponding probabilitiesiif there is a continuous range of values, such as an angle between 0 and 2 Continuous variable 連續(xù)型連續(xù)型隨機(jī)變量隨機(jī)變量: :在連

40、續(xù)區(qū)間取值，其取某確定值的概率由分在連續(xù)區(qū)間取值，其取某確定值的概率由分布密度函數(shù)給出布密度函數(shù)給出f(x)dx gives probability that X falls between x and x +dx.Random Variable A variable X that takes “random” values. We assume that it follows a probability distribution, p(x). Discrete variable: p1, p2, Continuous variable: p(x)dx gives probability th

41、at X falls between x and x +dx.11iip( )1p t dtProbability distribution functions (PDF).Expectation Value If the probability distribution is known, the expectation value (average value) can be computed as (for continuous variable)f()(f()f( )p( )XEXxx dxExpectation Value 期望值期望值3.5Expectation Value7Sta

42、tistic of a Random VariableVariance 方差方差: 2 = - 2Correlation 相關(guān)相關(guān): C= - Standard deviation (standard error) of a random number is the square root of its variance隨機(jī)變量隨機(jī)變量x分布對期望值的離散程度分布對期望值的離散程度Variance 2 = - 2Variance 2 = - 22 =(1+4+9+16+25+36)/6-3.52=91/6-12.25=15.168-12.25=2.91VarianceVariance 2 = - 25.83Continuous Random VariablesFirst of all, if there is a continuous range of values, such as an angle between 0 and 2 , then the probability of getting exactly a speci