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

17.

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

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

x=rand(N)產(chǎn)生元素在(0,1)間隨機(jī)分布的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)均勻分布的隨機(jī)點(diǎn)：方法：在中均勻產(chǎn)生足夠多隨機(jī)點(diǎn)，當(dāng)位于橢圓內(nèi)的點(diǎn)數(shù)為1000時停止?？捎呻S機(jī)數(shù)序列r

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

xn。則誤差約，它并不能和一些高級的數(shù)值積分算法比擬，但對多維情況，MC方法卻很有吸引力。我們可產(chǎn)生一系列隨機(jī)數(shù)可簡單取3個隨機(jī)數(shù)構(gòu)成一個隨機(jī)點(diǎn)，即相應(yīng)地，一般地，例如：求其中在該區(qū)域中取5000個點(diǎn)，可得該積分約0.57.計(jì)算體積Pseudo-codeProgramvolume_regionIntegerparametern5000,iprt1000Realarrayintegeri,mRealvol,x,y,zCallrandom((rij))ForI=1tondoxri1yri2zri3Ifthenmm+1Ifmod(I,iprt)=0thenVolreal(m)/real(I)OutputI,volEndifEndforEndprogram計(jì)算體積之作業(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針。針與線的交點(diǎn)數(shù)除以投的次數(shù)等于2l/(πd)。Realfunctionprob(n,npts)Integeri,nptsLogicalbirthRealsumSum=0ForI=1tonptsdoIfbirth(n)thensumsum+1EndforProbsum/real(npts)Endfunction不妨取d＝l。產(chǎn)生隨機(jī)數(shù)u和，其中u為針中心距最近的線的距離，故u小于等于0.5，為夾角，根據(jù)對稱性，可取0到π/2。算時需要用上是否不自洽（姜冰）作業(yè)中子屏蔽問題NeutronShieldingproblem假設(shè)鉛墻長為5，中子在鉛中的平均自由程為1，中子與鉛原子碰撞后各向同性散射。令碰撞8次后中子能量耗盡，試求穿透鉛墻的中子的比例。暫不考慮垂直紙面的運(yùn)動，則中子的水平位移是。入口鉛墻（長為5）21點(diǎn)問題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)移問題,蒙特卡羅模擬幾乎是必須的計(jì)算方法將空間劃分為Nx*Ny*Nz個網(wǎng)格選取空間中的隨機(jī)點(diǎn)產(chǎn)生光子包(photonpackage),并此向隨機(jī)方向傳播考慮光子與粒子的隨機(jī)碰撞產(chǎn)生具有某種分布的隨機(jī)數(shù)試驗(yàn)粒子數(shù)值模擬研究粒子在電磁場中的運(yùn)動研究粒子加速等全粒子試驗(yàn)粒子研究微觀