蒙特卡羅(Monte Carlo)模擬講義課件

17.

蒙特卡羅(MonteCarlo)模擬蒙特卡羅方法，或稱計算機隨機模擬方法，是一種基于“隨機數(shù)”的計算方法。該方法源于美國在WWI后研制原子彈的“曼哈頓計劃”。S.M.Ulam(1908-1984)和該計劃的主持人之一、數(shù)學(xué)家馮·諾伊曼（J.vonNeumann）用馳名世界的賭城—摩納哥的MonteCarlo—來命名這種方法。蒙特卡羅模擬其基本思想很早以前就被人們所發(fā)現(xiàn)和利用：17世紀，用事件發(fā)生的“頻率”來決定事件的“概率”；19世紀人們用投針試驗的方法來決定π。高速計算機：用數(shù)學(xué)方法在計算機上大量模擬這樣的試驗必要性科技計算及社會中的問題比計算π要復(fù)雜得多。比如金融衍生產(chǎn)品（期權(quán)、期貨、掉期等）的定價及交易風(fēng)險估算，問題的維數(shù)（即變量的個數(shù)）可能高達數(shù)百甚至數(shù)千。對這類問題，難度隨維數(shù)的增加呈指數(shù)增長，這就是所謂的“維數(shù)的災(zāi)難”(CurseofDimensionality)，傳統(tǒng)的數(shù)值方法難以對付（即使使用速度最快的計算機）。MonteCarlo方法能很好地用來對付維數(shù)的災(zāi)難，因為該方法的計算復(fù)雜性不再依賴于維數(shù)。以前那些本來是無法計算的問題現(xiàn)在也能夠計算量。為提高方法的效率，科學(xué)家們提出了許多所謂的“方差縮減”技巧。另一類形式與MonteCarlo方法相似，但理論基礎(chǔ)不同的方法—“擬蒙特卡羅方法”(Quasi-MonteCarlo方法)—近年來也獲得迅速發(fā)展。我國數(shù)學(xué)家華羅庚、王元提出的“華—王”方法即是其中的一例。這種方法的基本思想是“用確定性的超均勻分布序列(數(shù)學(xué)上稱為LowDiscrepancySequences)代替MonteCarlo方法中的隨機數(shù)序列。對某些問題該方法的實際速度一般可比MonteCarlo方法提出高數(shù)百倍，并可計算精確度。/gkjqy/rkx/rd2.htm擬蒙特卡羅模擬許多計算機系統(tǒng)都有隨機數(shù)生成函數(shù)F90:callrandom_seedcallrandom_number(a)計算程序產(chǎn)生的隨機數(shù)不是真正的隨機數(shù)，它們是確定的，但看上去是隨機的，且能通過一些隨機性的檢驗，故常稱為偽隨機數(shù)若采用當(dāng)前時刻作隨機數(shù)種子，則也是不可重復(fù)的隨機數(shù)如對32位字長的計算機realprocedurerandom((xi))integerarray(li)nrealarray

(xi)nl0anyintegerthat1<l0<231-1fori=1tondoli=(231-1)除以16807li-1的余數(shù)xi

＝li/(231-1)endfor其它算法1.取初值x0(0,1)，letxi

bethefractionalpartof(π+xi-1)52.Letu0=1.Fori=1,2,3,…,letui

betheremainderof(8t-3)ui-1dividedby28,andxi

=ui/2i.Heretcanbeanylargeinteger注意：上述隨機數(shù)序列均具周期性，如上頁random子程序的周期約230。Maplehasacollectionofrandomnumbergenerators.如：With(stats)x

:=uniform(0..1):seq(x(),i=1..10);Matlab:x=rand(10,1)產(chǎn)生10個隨機數(shù)

x=rand(N)產(chǎn)生元素在(0,1)間隨機分布的N*N矩陣

s=rand(‘state’,0)重設(shè)該生成函數(shù)到初始狀態(tài)如(b-a)r+a

x

(a,b)integer((n+1)r)x

{0,1,2,……n}試產(chǎn)生1000個在橢圓x2+4y2=4內(nèi)均勻分布的隨機點：方法：在中均勻產(chǎn)生足夠多隨機點，當(dāng)位于橢圓內(nèi)的點數(shù)為1000時停止?？捎呻S機數(shù)序列r

(0,1)構(gòu)造一系列隨機數(shù)序列有時隨機數(shù)產(chǎn)生函數(shù)通不過嚴格的隨機性測試。而在一些計算（如MC積分）中隨機性非常重要。故使用更長字節(jié)的計算機更好。應(yīng)用之一估計面積和體積Numericalintegration選取(0,1)中隨機數(shù)序列x1，x2，x3，……

xn。則誤差約，它并不能和一些高級的數(shù)值積分算法比擬，但對多維情況，MC方法卻很有吸引力。我們可產(chǎn)生一系列隨機數(shù)可簡單取3個隨機數(shù)構(gòu)成一個隨機點，即相應(yīng)地，一般地，例如：求其中在該區(qū)域中取5000個點，可得該積分約0.57.計算體積Pseudo-codeProgramvolume_regionIntegerparametern5000,iprt1000Realarrayintegeri,mRealvol,x,y,zCallrandom((rij))ForI=1tondoxri1yri2zri3Ifthenmm+1Ifmod(I,iprt)=0thenVolreal(m)/real(I)OutputI,volEndifEndforEndprogram計算體積之作業(yè)：冰淇淋的體積應(yīng)用之二MC模擬生日問題假設(shè)有n個人在一起，各自的生日為365天之一，根據(jù)概率理論，與很多人的直覺相反，只需23個人便有大于50％的幾率人群中至少有2個人生日相同。n理論幾率模擬幾率0.1170.1100.4110.4120.5070.5200.7060.6920.9410.936500.9860.987Realfunctionprob(n,npts)Integeri,nptsLogicalbirthRealsumSum=0ForI=1tonptsdoIfbirth(n)thensumsum+1EndforProbsum/real(npts)EndfunctionLogicalfunctionbirth(n)……蒲豐的投針問題Buffon’sNeedleproblem蒲豐證明在相距d的一系列平行線上投長為l針。針與線的交點數(shù)除以投的次數(shù)等于2l/(πd)。Realfunctionprob(n,npts)Integeri,nptsLogicalbirthRealsumSum=0ForI=1tonptsdoIfbirth(n)thensumsum+1EndforProbsum/real(npts)Endfunction不妨取d＝l。產(chǎn)生隨機數(shù)u和，其中u為針中心距最近的線的距離，故u小于等于0.5，為夾角，根據(jù)對稱性，可取0到π/2。算時需要用上是否不自洽（姜冰）作業(yè)中子屏蔽問題NeutronShieldingproblem假設(shè)鉛墻長為5，中子在鉛中的平均自由程為1，中子與鉛原子碰撞后各向同性散射。令碰撞8次后中子能量耗盡，試求穿透鉛墻的中子的比例。暫不考慮垂直紙面的運動，則中子的水平位移是。入口鉛墻（長為5）21點問題Blackjack21problem牛刀小試參數(shù)擬合問題Parameterfittingproblem輻射轉(zhuǎn)移問題Radiativetransferproblem[1]Ahn,S.-H.,Lee,H.-W.&Lee,H.-M.2000,JKAS,33,p29.[2]

ZhengZheng&JordiMiralda-Escude.2002,ApJ,578,p33.[3]Watson,A.M.&Henney,W.J.2001,RMxAA,37,p221.[4]LorneW.Avery&LewisL.H.,1968,ApJ,152,p493.[5]LawrenceJ.C.,PeterD.N.&JefferyD.S.1972,ApJ,176,439.[6]DavidA.Neufeld,1990,ApJ,350,p216.輻射轉(zhuǎn)移問題對三維輻射轉(zhuǎn)移問題,蒙特卡羅模擬幾乎是必須的計算方法將空間劃分為Nx*Ny*Nz個網(wǎng)格選取空間中的隨機點產(chǎn)生光子包(photonpackage),并此向隨機方向傳播考慮光子與粒子的隨機碰撞產(chǎn)生具有某種分布的隨機數(shù)試驗粒子數(shù)值模擬研究粒子在電磁場中的運動研究粒子加速等全粒子試驗粒子研究微觀