1、例題,例 求前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(IN)EXIT END DO PRINT *,SUM=,S END PROGRAM MAIN,程序的運(yùn)行結(jié)果如下 15 SUM= 3.182290,用DO結(jié)構(gòu)處理循環(huán)問題，最關(guān)鍵的是要處理好初值、循環(huán)終止條件和循環(huán)體語句的關(guān)系。,Finite differen

2、ce approximation of differential equations,A differential equation can be approximated by a finite difference scheme. For example,Forward time,Backward space,Central space,Forward time-,Central space,FICK 第二定律,Fick 第二定律穩(wěn)態(tài)擴(kuò)散解,REAL C0(0:1000+1),C(0:1000+1) !C(DISTANCE) C0=0.1 C0(0)=0.8 !BOUNDARY CONDI

3、TION C0(1001)=0.1 !BOUNDARY CONDITION C=C0 OPEN(1,FILE=F:DIF.dat) DO JT=1,50000 ！Time DO IX=1,1000 ! distance C(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 DO C0=C END DO PRINT *, PROGRAMME IS FINISHED CLOSE(1) END,蒙特卡羅(Monte Carlo)方法,1、簡介 2

4、、相關(guān)幾個(gè)例子,Monte Carlo方法：,亦稱統(tǒng)計(jì)模擬方法，statistical simulation method 利用隨機(jī)數(shù)進(jìn)行數(shù)值模擬的方法,Monte Carlo名字的由來：,是由Metropolis在二次世界大戰(zhàn)期間提出的：Manhattan計(jì)劃，研究與原子彈有關(guān)的中子輸運(yùn)過程；,Monte Carlo是摩納哥（monaco)的首都，該城以賭博聞名,Nicholas Metropolis (1915-1999),Monte-Carlo, Monaco,Monte Carlo方法的基本思想很早以前就被人們所發(fā)現(xiàn)和利用。早在17世紀(jì)，人們就知道用事件產(chǎn)生的“頻率”來近似事件的“概

5、率”。 18世紀(jì)人們用投針試驗(yàn)的方法來決定圓周率。本世紀(jì)40年代電子計(jì)算機(jī)的出現(xiàn)，特別是近年來高速電子計(jì)算機(jī)的出現(xiàn)，使得用數(shù)學(xué)方法在計(jì)算機(jī)上大量、快速地模擬這樣的試驗(yà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ī)的重要性,合作起源于Ma

6、nhattan工程：利用ENIAC(Electronic Numerical Integrator and Computer)計(jì)算產(chǎn)額,Monte Carlo模擬的應(yīng)用：,自然現(xiàn)象的模擬：,實(shí)驗(yàn)探測器的模擬,數(shù)值分析：,利用Monte Carlo方法求積分,材料學(xué)中的應(yīng)用：,擴(kuò)散、晶粒的長大、再結(jié)晶等,Monte Carlo模擬在物理研究中的作用,Monte Carlo模擬的步驟：,根據(jù)欲研究的物理系統(tǒng)的性質(zhì)，建立能夠描述該系統(tǒng)特性的理論模型，導(dǎo)出該模型的某些特征量的概率密度函數(shù)；,從概率密度函數(shù)出發(fā)進(jìn)行隨機(jī)抽樣，得到特征量的一些模擬結(jié)果；,對模擬結(jié)果進(jìn)行分析總結(jié)，預(yù)言物理系統(tǒng)的某些特性。,注

7、意以下兩點(diǎn)：,Monte Carlo方法與數(shù)值解法的不同:,Monte Carlo方法利用隨機(jī)抽樣的方法來求解物理問題;,數(shù)值解法:從一個(gè)物理系統(tǒng)的數(shù)學(xué)模型出發(fā),通過求解一系列的微分方程來的導(dǎo)出系統(tǒng)的未知狀態(tài);,Monte Carlo方法并非只能用來解決包含隨機(jī)的過程的問題:,許多利用Monte Carlo方法進(jìn)行求解的問題中并不包含隨機(jī)過程 例如:用Monte Carlo方法計(jì)算定積分. 對這樣的問題可將其轉(zhuǎn)換成相關(guān)的隨機(jī)過程, 然后用Monte Carlo方法進(jìn)行求解,Monte Carlo算法的主要組成部分,概率密度函數(shù) 必須給出描述一個(gè)物理系統(tǒng)的一組概率密度函數(shù);,隨機(jī)數(shù)產(chǎn)生器能夠產(chǎn)生

8、在區(qū)間0,1上均勻分布的隨機(jī)數(shù),抽樣規(guī)則如何從在區(qū)間0,1上均勻分布的隨機(jī)數(shù)出發(fā),隨機(jī)抽取服從給定的pdf的隨機(jī)變量;,模擬結(jié)果記錄記錄一些感興趣的量的模擬結(jié)果,誤差估計(jì)必須確定統(tǒng)計(jì)誤差（或方差）隨模擬次數(shù)以及其它一些量的變化；,減少方差的技術(shù)利用該技術(shù)可減少模擬過程中計(jì)算的次數(shù)；,并行和矢量化可以在先進(jìn)的并行計(jì)算機(jī)上運(yùn)行的有效算法,確定性系統(tǒng),隨機(jī)性系統(tǒng),模擬,自然界,Monte-Carlo模擬，即隨機(jī)模擬（重復(fù)“試驗(yàn)”）,重復(fù)試驗(yàn),計(jì)算機(jī)模擬,蒙特卡羅方法的基本原理及思想,當(dāng)所要求解的問題是某種事件出現(xiàn)的概率，或者是某個(gè)隨機(jī)變量的期望值時(shí)，它們可以通過某種“試驗(yàn)”的方法，得到這種事件出現(xiàn)的

9、頻率，或者這個(gè)隨機(jī)變數(shù)的平均值，并用它們作為問題的解。 把蒙特卡羅解題歸結(jié)為三個(gè)主要步驟： 構(gòu)造或描述概率過程； 實(shí)現(xiàn)從已知概率分布抽樣； 建立各種估計(jì)量。,對于本身就具有隨機(jī)性質(zhì)的問題，主要是正確描述和模擬這個(gè)概率過程；對于不是隨機(jī)性質(zhì)的確定性問題，比如計(jì)算定積分，就必須事先構(gòu)造一個(gè)人為的概率過程。最簡單、最基本、最重要的一個(gè)概率分布是(0,1)上的均勻分布（或稱矩形分布）。 實(shí)現(xiàn)從已知概率分布抽樣(模擬試驗(yàn))：構(gòu)造了概率模型以后，產(chǎn)生已知概率分布的隨機(jī)變量（或隨機(jī)向量），就成為實(shí)現(xiàn)蒙特卡羅方法模擬實(shí)驗(yàn)的基本手段，這也是蒙特卡羅方法被稱為隨機(jī)抽樣的原因。 建立各種估計(jì)量： 構(gòu)造了概率模型并能

10、從中抽樣后，即實(shí)現(xiàn)模擬實(shí)驗(yàn)后，就要確定一個(gè)隨 機(jī)變量，作為所要求的問題的解，我們稱它為無偏估計(jì)。,收斂性: 由大數(shù)定律, Monte-Carlo模擬的收斂是以概率而言的 誤差: 用頻率估計(jì)概率時(shí)誤差的估計(jì)，可由中心極限定理，給定置信水平 的條件下，有： 模擬次數(shù):由誤差公式得,隨機(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)注。,(0,1均勻分布

11、隨機(jī)數(shù)的產(chǎn)生,最常用的是在(0,1區(qū)間內(nèi)均勻分布的隨機(jī)數(shù)，其他分布的隨機(jī)數(shù)可利用均勻分布隨機(jī)數(shù)來產(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ù)。,例1：,R=RAND() ISEED=RTC() R=RAN(ISEED),其他各種分布的隨機(jī)數(shù)的產(chǎn)生,基本方法有如下三種： 逆變換法 合成法 篩選法,逆變換法,設(shè)隨機(jī)變量 的分布函數(shù)為 ，定義 定理 設(shè)隨機(jī)變量 服從 上的均勻分布，則

12、 的分布函數(shù)為 。 因此，要產(chǎn)生來自 的隨機(jī)數(shù)，只要先產(chǎn)生來自 的隨機(jī)數(shù)，然后計(jì)算 即可。 其步驟為,合成法,合成法的應(yīng)用最早見于Butlter 的書中。 構(gòu)思如下： 如果 的密度函數(shù) 難于抽樣，而 關(guān)于 的條件密度函數(shù) 以及 的密度函數(shù) 均易于抽樣，則 的隨機(jī)數(shù)可如下產(chǎn)生： 可以證明由此得到 的服從 。,篩選抽樣,假設(shè)我們要從 抽樣，如果可以將 表示成 ，其中 是一個(gè)密度函數(shù)且易于抽樣，而 ， 是常數(shù)，則 的抽樣可如下進(jìn)行： 定理 設(shè) 的密度函數(shù) ，且 ，其中 ， ， 是一個(gè)密度函數(shù)。令 和 分別服從 和 ，則在 的條件 下， 的條件密度為,標(biāo)準(zhǔn)正態(tài)分布的隨機(jī)數(shù),的隨機(jī)數(shù)產(chǎn)生方法很多。簡要介

13、紹三種。 法1、 變換法（Box 和Muller 1958） 設(shè) ， 是獨(dú)立同分布的 變量，令 則 與 獨(dú)立，均服從標(biāo)準(zhǔn)正態(tài)分布。 法2、 結(jié)合合成法與篩選法。（略） 法3、 近似方法（利用中心極限定理） 即用 個(gè) 變量產(chǎn)生一個(gè) 變量。 其中 是抽自 的隨機(jī)數(shù)， 可近似為一個(gè) 變量。,由于正態(tài)隨機(jī)數(shù)的獨(dú)立性，當(dāng) 很小時(shí)，按(a)、(b)所構(gòu)成的過程的相關(guān)時(shí)間很短，從而具有很高的截止頻率。 當(dāng)然， 過小，樣本的計(jì)算量過大。因此，選取 適當(dāng)小即可。,關(guān)于隨機(jī)數(shù)的幾點(diǎn)注,注1 由于均勻分布的隨機(jī)數(shù)的產(chǎn)生總是采用某個(gè)確定的模型進(jìn)行的，從理論上講，總會(huì)有周期現(xiàn)象出現(xiàn)的。初值確定后，所有隨機(jī)數(shù)也隨之確定，

14、并不滿足真正隨機(jī)數(shù)的要求。因此通常把由數(shù)學(xué)方法產(chǎn)生的隨機(jī)數(shù)成為偽隨機(jī)數(shù)。 但其周期又相當(dāng)長，在實(shí)際應(yīng)用中幾乎不可能出現(xiàn)。因此，這種由計(jì)算機(jī)產(chǎn)生的偽隨機(jī)數(shù)可以當(dāng)作真正的隨機(jī)數(shù)來處理。 注2 應(yīng)對所產(chǎn)生的偽隨機(jī)數(shù)作各種統(tǒng)計(jì)檢驗(yàn)，如獨(dú)立性檢驗(yàn)，分布檢驗(yàn)，功率譜檢驗(yàn)等等。,Monte Carlo方法相關(guān)引例:,1、Buffon投針實(shí)驗(yàn)：,18世紀(jì)，法國數(shù)學(xué)家Comte de Buffon利用投針實(shí)驗(yàn)估計(jì)的值,Problem of Buffons needle:,If a needle of length l is dropped at random on the middle of a horizon

15、tal surface ruled with parallel lines a distance dl apart, what is the probability that the needle will cross one of the lines?,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 t

16、he region 0,d)0,). Accordingly, it has probability density function 1/d.,The probability that the needle will cross one of the lines is given by the integral,圓周率的值 = 3.14159 26535 89793 23846 26433 83279 50288 41971 69399 3751058209 74944 59230 78164 06286 20899 86280 34825 34211 7067982148 08651 32

17、823 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 01133 0

18、5305 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 56082

19、77857 71342 75778 96091 73637 1787214684 40901 22495 34301 46549 58537 10507 92279 68925 8923542019 95611 21290 21960 86403 44181 59813 62977 47713 .,什么是蒙特-卡洛模擬？,計(jì)算機(jī)實(shí)驗(yàn),物理實(shí)驗(yàn) 用探測器取數(shù)據(jù),5,RHIC-STAR 實(shí)驗(yàn)探測器,蒙特-卡洛模擬 用計(jì)算機(jī)取數(shù)據(jù),一個(gè)例子 道爾頓板,理論 分布曲線 實(shí)驗(yàn) 丟鐵球 模擬 計(jì)算機(jī)丟球,6,模擬的出發(fā)點(diǎn) 模型,根據(jù)球和釘子碰撞的規(guī)律 追蹤每個(gè)球的運(yùn)動(dòng) 根據(jù)落入每個(gè)格子的幾率 給出這一幾率

20、的每次實(shí)現(xiàn),8,9,令 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, m s=s+p (i) r(i+1)=s enddo,do 3 k=1, m np(k)=0 do 3 j=1,n if(ra( j ).ge.r(k).and.ra( j ).lt.r(k+1) r(k) ra (j ) r(k+1) np(k)=np(k)+1 Enddo np(k) 是落入第 k 個(gè)盒子中的粒子數(shù)。,r(m) p(i)=1 r(3)

21、 p(1)+p(2) r(2) p(1) r(1) 0, ,i=1,m, ,通過比較第 j 個(gè)鐵球的 ra(j) 和 r(i) 來決定這個(gè)球落入那一格子中：,例一：道爾頓板的蒙特-卡洛模擬,p3,p4,p5,p2,p1,p6,p7,例二 面積的計(jì)算,f (x),x,辛普遜方法,I = Sn,蒙特-卡洛方法,11,考慮平面上的一個(gè)邊長為1的正方形及其內(nèi)部的一個(gè)形狀不規(guī)則的“圖形”，如何求出這個(gè)“圖形”的面積呢？Monte Carlo方法是這樣一種“隨機(jī)化”的方法：向該正方形“隨機(jī)地”投擲N個(gè)點(diǎn)，若有M個(gè)點(diǎn)落于“圖形”內(nèi)，則該“圖形”的面積近似為M/N。,用該方法計(jì)算的基本思路是： 1根據(jù)圓面積的

22、公式：s=R2,當(dāng)R=1時(shí)，S=。由于圓的方程是：x2+y2=1(x2為x的平方的意思),因此1/4圓面積為x軸、y軸和上述方程所包圍的部分。如果在1*1的正方形中均勻地落入隨機(jī)點(diǎn)，則落入1/4圓中的點(diǎn)的概率就是1/4圓的面積。其4倍，就是圓面積。由于半徑為1，該面積的值為的值。,例1 算圓周率,REAL X,Y,R,PI ISEED=RTC() N=3000000 DO 10 I=1,N X=RAN(ISEED) Y=RAN(ISEED) R=SQRT(X*X+Y*Y) IF(R.LT.1.0)N0=N0+1 10CONTINUE PI=4.0*N0/N WRITE(*,*)YUAN ZHO

23、U LV IS,PI END,例、 定積分（面積）的MC計(jì)算 隨機(jī)投點(diǎn)法 樣本平均值法,定積分的MC計(jì)算,事實(shí)上，不少的統(tǒng)計(jì)問題，如計(jì)算概率、各階距等，最后都?xì)w結(jié)為定積分的近似計(jì)算問題。 下面考慮一個(gè)簡單的定積分 為了說明問題，我們首先介紹兩種求 的簡單的MC方法，然后給出幾種較為復(fù)雜而更有效的MC方法。,在計(jì)算積分上，MC的實(shí)用場合是計(jì)算重積分 其中 是 維空間的點(diǎn)，當(dāng) 較大時(shí)，用MC方法比一般的數(shù)值方法有優(yōu)點(diǎn)，主要是它的誤差與維數(shù) 無關(guān)。,面積正比于點(diǎn)數(shù)?？傸c(diǎn)數(shù)n。點(diǎn)（x,y）,那么我們可以得到 的一個(gè)估計(jì) 具體試驗(yàn)步驟為,注1 隨機(jī)投點(diǎn)法的思想簡單明了，且每次投點(diǎn)結(jié)果服從二項(xiàng)分布，故 ，其中 注2 可證 是 的無偏估計(jì)。若用估計(jì)的標(biāo)準(zhǔn)差來衡量其精度，則估計(jì) 的精度的階為 。 注3 這里，定積分的解，就對應(yīng)我們選定的隨機(jī)變量的概率值。,樣本平均值法,基本原理：對積分 ，設(shè) 是 上的一個(gè)密度函數(shù)，改寫 可見，任一積分均可以表示為某個(gè)隨機(jī)變量 （函數(shù)）的期望。由矩法，若有n個(gè)來自a,b的觀測值，則可給出 的一個(gè)矩估計(jì)。 最簡單的，若n有限，可取 。,設(shè) 是來自 的隨機(jī)數(shù)，則 的一個(gè)估計(jì)為 具體步驟為 注 可證 是 的無偏估計(jì)