蒙特卡羅方法1課件

計算物理3/lesson/ComputationalPhysics蒙特卡羅方法計算物理3/less1蒙特卡羅方法蒲豐投針收斂性、誤差和優(yōu)缺點任意分布的隨機數(shù)粒子輸運問題隨機過程模擬梅氏抽樣√蒙特卡羅方法蒲豐投針√2蒲豐投針(1/5)蒙特卡羅方法又稱隨機抽樣技巧或統(tǒng)計試驗方法以概率統(tǒng)計理論為基礎(chǔ)的能夠比較逼真地描述事物的特點及物理實驗過程解決一般數(shù)值方法難以解決的問題隨著電子計算機的發(fā)展而發(fā)展首先在核武器的試驗與研制中得到了應(yīng)用蒲豐投針法國數(shù)學家蒲豐的1777年出版的著作：“在平面上畫有一組間距為d的平行線，將一根長度為

l

(l<d)

的針任意擲在這個平面上，求此針與平行線中任一條相交的概率?！薄唐沿S投針(1/5)蒙特卡羅方法蒲豐投針√3蒲豐投針(2/5)步驟在桌面上畫出間距為

2d的平行線準備長度為

2l

(l<d)

的針向桌面隨機投針如果針與平行線相交，則計數(shù)器n加

1計算：計數(shù)器n與總投針數(shù)N的比例(視作相交概率P)概率分析

P

=

?各條平行線地位等同，僅考慮某條平行線附近的情況平行線方向的x坐標對概率沒影響針的中點的

y坐標在線之間等概率落入(均勻分布在

[0,d])，僅當

yl才可能針-線相交針-線的夾角

q

均勻分布在

[0,

p]，q

與

y獨立xy√蒲豐投針(2/5)步驟概率分析P=?xy√4蒲豐投針(3/5)xy概率

P

=

2l/

(pd)，可計算圓周率實驗者時間針長總投數(shù)相交數(shù)p

值Wolf18500.85,0002,5323.1596Smith18550.63,2041,2183.1554DeMorgan,C18601.06003823.137Fox18840.751,0304893.1595Lezzerini19010.833,4081,8083.1415929Reina19250.54192,5208593.1795√蒲豐投針(3/5)xy概率P=2l/(pd)，5蒲豐投針(4/5)關(guān)于蒙特卡羅方法的分析和總結(jié)基本思想確定所求問題的解是某事件的概率(或某隨機變量的數(shù)學期望、或與概率/數(shù)學期望有關(guān)的量，如

p

)通過試驗方法，得出事件發(fā)生的頻率(或該隨機變量若干個具體觀察值的算術(shù)平均值，如

P)，求解數(shù)學期望與概率：當隨機變量的取值僅為1或0時，它的數(shù)學期望就是事件的概率；反之亦然數(shù)學期望與算術(shù)平均值用隨機試驗的方法計算積分，將積分看作服從分布密度函數(shù)f(r)的隨機變量

g(r)的數(shù)學期望通過試驗，得到N個觀察值

r1,

r2,

…,rN

(從

f(r)

中抽取N個子樣

r1,

r2,

…,rN

)，求

g

(r)

的算術(shù)平均值√蒲豐投針(4/5)關(guān)于蒙特卡羅方法的分析和總結(jié)通過試驗，得到6蒲豐投針(5/5)試驗方法和次數(shù)試驗方法不一定可行精確的近似解需要巨量的試驗次數(shù)，但人工方法作大量的試驗相當困難，甚至是不可能的用計算機模擬隨機試驗過程，完成巨量的次數(shù)抽象x的分布密度函數(shù)：q

的分布密度函數(shù)：產(chǎn)生任意的(x,q)

=

由

f1(x)抽樣

x

+

由

f2(q)抽樣

q對應(yīng)的隨機變量：數(shù)學期望與算術(shù)平均值新的問題：誤差？收斂？√蒲豐投針(5/5)試驗方法和次數(shù)q的分布密度函數(shù)：產(chǎn)生任意7收斂性、誤差和優(yōu)缺點(1/4)收斂性求解：以隨機變量X

的簡單子樣X1,

X2,…,

XN

的算術(shù)平均值，作為求解的近似值近似值的收斂性大數(shù)定理：當試驗次數(shù)足夠多時，事件發(fā)生的頻率無窮接近于該事件發(fā)生的概率如果

X1,…,

XN

獨立分布，且期望值有限(

E(X)<

)，那么

當隨機變量X

的簡單子樣數(shù)N充分大時，其算術(shù)平均值以概率

1

收斂于期望值

E(X)√收斂性、誤差和優(yōu)缺點(1/4)收斂性近似值的收斂性當隨機8收斂性、誤差和優(yōu)缺點(2/4)誤差概率論的中心極限定理：如果隨機變量序列X1,

X2,…,

XN獨立分布，且具有有限非零的方差s

2，即其中的

f(x)是X的分布密度函數(shù)，則當N足夠大時(蒙特卡羅方法)不等式的概率約

1-a誤差定義為，收斂速度為a0.5000.0500.003la0.6741.9602.968√收斂性、誤差和優(yōu)缺點(2/4)誤差其中的f(x)是9收斂性、誤差和優(yōu)缺點(3/4)蒙特卡羅方法的誤差為概率誤差均方差

s

是未知的，估計值為減少誤差的技巧(在確定的置信度

a前提下)誤差

e

與試驗次數(shù)的開根號

N1/2

成反比：精度一個數(shù)量級，次數(shù)N兩個數(shù)量級——巨大的代價誤差

e

與均方差

s

成正比：精度一個數(shù)量級，均方差

s

一個數(shù)量級——可接受的代價效率降低方差增加觀察子樣的時間固定時間內(nèi)樣本數(shù)減少代價增加蒙特卡羅方法中的效率：由均方差

s

和觀察一個子樣的費用(計算機時)c衡量

=

s2

c√收斂性、誤差和優(yōu)缺點(3/4)蒙特卡羅方法的誤差為概率誤差減10收斂性、誤差和優(yōu)缺點(4/4)優(yōu)點較逼真地描述具有隨機性質(zhì)的事物的特點及物理實驗過程受幾何條件限制小收斂速度與問題的維數(shù)無關(guān)具有同時計算多個方案與多個未知量的能力誤差容易確定程序結(jié)構(gòu)簡單，易于實現(xiàn)缺點收斂速度慢誤差具有概率性在某些問題中，計算結(jié)果與系統(tǒng)大小有關(guān)主要應(yīng)用范圍粒子輸運，統(tǒng)計物理，典型數(shù)學，真空技術(shù)，激光技術(shù)，醫(yī)學，生物，探礦等√收斂性、誤差和優(yōu)缺點(4/4)優(yōu)點√11任意分布的隨機數(shù)(1/11)隨機數(shù)(蒙特卡羅方法的關(guān)鍵)基本的：均勻地在

(0,1)

內(nèi)分布任意的：非均勻地在

(a,b)

內(nèi)分布方法：直接的(離散、連續(xù))、舍選的(簡單、乘、乘加)直接抽樣方法離散型：產(chǎn)生隨機量

x

的抽樣值

xi

，概率為Pi

(i=

1,

2,…)方法：先計算x的分布函數(shù)

對于產(chǎn)生的R，如果滿足Fi-1<

R

Fi，則令抽樣x=

xi

證明：√任意分布的隨機數(shù)(1/11)隨機數(shù)(蒙特卡羅方法的關(guān)鍵)方法12任意分布的隨機數(shù)(2/11)例：二項分布

直接抽樣方法：

產(chǎn)生R，如果滿足，則令x

=

n例：泊松分布

直接抽樣方法：

產(chǎn)生R，如果滿足，則令x=n√任意分布的隨機數(shù)(2/11)例：二項分布直接抽樣13任意分布的隨機數(shù)(3/11)連續(xù)型：產(chǎn)生隨機量

x

的抽樣值，概率密度函數(shù)為f(x)方法：先計算x的分布函數(shù)

對于隨機數(shù)R'，解方程

R'

=

F(x')，得到

x'

=

F-1(R')，則令x=

x'

證明：√任意分布的隨機數(shù)(3/11)連續(xù)型：產(chǎn)生隨機量x的抽樣值14任意分布的隨機數(shù)(4/11)例：均勻介質(zhì)中，粒子運動的自由程S是隨機量，其概率密度函數(shù)為

f(S)

=

Se

-SS；其中的S

=

ns是宏觀總截面，s是原子截面，n為介質(zhì)中的原子數(shù)密度分布函數(shù)：介質(zhì)中粒子的碰撞過程

隨機量

R中的抽樣過程例：散射方位角余弦分布√任意分布的隨機數(shù)(4/11)例：均勻介質(zhì)中，粒子運動的自由程15舍任意分布的隨機數(shù)(5/11)舍選抽樣方法直接法的特點：簡單，但分布函數(shù)的反函數(shù)不一定有解析解簡單分布定理：如果Z是

[a,b]

上均勻分布的隨機數(shù)，那么利用條件f(Z)/MR

(M是f(Z)

的上界)選出的Z將是

[a,b]

上概率密度函數(shù)為f(Z)

的隨機數(shù)

證明：如右圖所示dZ(Z',R)RZf(Z)

/

Mab10Z為橫坐標，R為縱坐標，實曲線為函數(shù)f(Z)/MR在(0,1)內(nèi)、Z在[a,b]內(nèi)均勻分布隨機點(Z',R)在虛框內(nèi)均勻分布隨機點落入窄條的概率=兩面積之比選已知概率密度函數(shù)抽樣方法√舍任意分布的隨機數(shù)(5/11)舍選抽樣方法證明：如右圖所16任意分布的隨機數(shù)(6/11)簡單分布抽樣方法的流程產(chǎn)生隨機數(shù)x和R計算Z

=

a+(b-a)xf(Z)

/

MR?Z'

=

ZYN效率(有一部分隨機數(shù)被舍棄)舍dZ(Z',R)RZf(Z)

/

Mab10選例：受限的散射方位角余弦分布

效率√任意分布的隨機數(shù)(6/11)簡單分布抽樣方法的流程產(chǎn)生隨機數(shù)17任意分布的隨機數(shù)(7/11)乘分布：f(x)

有銳鋒時效率很低，需要改進的簡單分布方法：將函數(shù)寫成

f(Z)

=

h(Z)

f1(Z)，由容易抽樣的

f1(Z)抽樣出Z，代入h(Z)，如果滿足條件

h(Z)/MR

(M是

h(Z)

的上界)，那么得到概率密度函數(shù)為f(Z)

的抽樣值舍

證明：如右圖所示dF1(F1',R)RF1(Z)h(Z)

/

M110f1(Z)的分布函數(shù)F1(Z)

為橫坐標，R為縱坐標，實曲線為函數(shù)

h(Z)/MR和F1

在(0,1)內(nèi)均勻分布隨機點(F1',R)在虛框內(nèi)均勻分布隨機點落入窄條的概率=兩面積之比選效率√任意分布的隨機數(shù)(7/11)乘分布：f(x)有銳鋒時效率很18任意分布的隨機數(shù)(8/11)乘分布抽樣方法的流程產(chǎn)生隨機數(shù)F1和R利用F1，由f1(Z)抽樣Zh(Z)

/

MR?Z'

=

ZYN例：半正態(tài)分布，概率密度函數(shù)√任意分布的隨機數(shù)(8/11)乘分布抽樣方法的流程產(chǎn)生隨機數(shù)利19任意分布的隨機數(shù)(9/11)乘加分布方法：如果隨機量x具有以下乘加形式

其中f1(x)和f2(x)滿足密度函數(shù)的要求，即

并且h1(x)和h2(x)滿足

那么從f(x)的以下變形，可以得到乘加分布

其中

和

為加權(quán)因子

證明：參照直接法和舍選法(略)√任意分布的隨機數(shù)(9/11)乘加分布其中f1(x)和20任意分布的隨機數(shù)(10/11)乘加分布抽樣方法的流程產(chǎn)生隨機數(shù)g,

F,R利用F由f1(x)抽樣x1gg1

?YNh1(x1)/M1

R?利用F由f2(x)抽樣x2h2(x2)/M2

R?YYNNx

=

x1x

=

x2效率√任意分布的隨機數(shù)(10/11)乘加分布抽樣方法的流程產(chǎn)生隨機21任意分布的隨機數(shù)(11/11)例：散射光子能量抽樣。能量為

E0

的入射光子，經(jīng)原子散射后的能量E按某種概率分布。如果令

x

=

E0/E，那么概率密度函數(shù)

f(x)

為(其中K=K(E0)為歸一因子)√任意分布的隨機數(shù)(11/11)例：散射光子能量抽樣。能量為22粒子輸運問題(1/7)蒙特卡羅模擬粒子輸運是隨機過程，運動規(guī)律是大量粒子運動的統(tǒng)計蒙特卡羅模擬：模擬一定數(shù)量粒子在介質(zhì)中的運動，再現(xiàn)粒子運動的統(tǒng)計規(guī)律例(平板介質(zhì)模型)：由單一物質(zhì)組成的均勻介質(zhì)，厚度為H，能量為

E00

的平行光子束垂直射入板內(nèi)，求光子對板的投射率H抽樣：自由程、作用類型、散射能量、散射方向自由程抽樣(PPT15)：粒子運動的自由程S是隨機量，其概率密度函數(shù)為

f(S)

=

Se

-SS；其中的S

=

ns是宏觀總截面，s是原子截面，n為介質(zhì)中的原子數(shù)密度分布函數(shù)：F(S')

=

1

-

e-SS'粒子的碰撞過程

隨機量

R中的抽樣過程：S'

=

-lnR'

/

S√粒子輸運問題(1/7)蒙特卡羅模擬例(平板介質(zhì)模型)：由單一23粒子輸運問題(2/7)作用類型抽樣光子-介質(zhì)：散射(康普頓散射)和吸收(光電效應(yīng))

散射截面

吸收截面

總截面：St

=

Ss

+

Sa光子隨機地被散射或吸收的過程隨機量

R中的抽樣過程：如果

Ss

/

St

R，則散射，否則為吸收√粒子輸運問題(2/7)作用類型抽樣吸收截面總截面：S24粒子輸運問題(3/7)散射能量抽樣(PPT22)：能量為

E0

的入射光子，經(jīng)原子散射后的能量E按某種概率分布。如果令

x

=

E0/E，那么概率密度函數(shù)

f(x)

為(其中K=K(E0)為歸一因子)√粒子輸運問題(3/7)散射能量抽樣(PPT22)：能量為E25粒子輸運問題(4/7)散射方向抽樣坐標系：以入射方向(由

q

和

f

確定)為參考系，散射方向

由散射角

q'

和散射方位角

f'

確定xyz流程圖輸入g0,E0,E計算m

=1+1/E0-1/E計算g產(chǎn)生隨機數(shù)R計算f'

=

2pR√粒子輸運問題(4/7)散射方向抽樣坐標系：以入射方向(26粒子輸運問題(5/7)直接模擬方法模擬粒子在介質(zhì)中的真實物理過程粒子在介質(zhì)中的狀態(tài)：空間位置，能量和運動方向碰撞點的狀態(tài)參數(shù)

(Zm,

Em,

gm)

表示從源出發(fā)的粒子在介質(zhì)中經(jīng)過m次碰撞后的狀態(tài)ZHOZ0E0q0Z1E1q1Z2E2q2粒子的運動過程

=

碰撞點的狀態(tài)序列

模擬粒子的運動過程

=確定狀態(tài)序列的問題√粒子輸運問題(5/7)直接模擬方法碰撞點的狀態(tài)參數(shù)(Zm,27粒子輸運問題(6/7)記錄和計算的內(nèi)容穿透率和誤差估計(設(shè)N為入射粒子數(shù)，N1

為透射數(shù))透射粒子的能量和方向分布離散化能量：Emin

=

E0

<E1

<

<Ei

<

<

EI離散化角度：0

=

q0

<

q1

<

<

qj

<

<

qJ

=

p

/

2能量和方向分布粒子的軌跡粒子運動關(guān)于Z軸對稱，只需要記錄(Zim,

qim)流程圖√粒子輸運問題(6/7)記錄和計算的內(nèi)容透射粒子的能量和方向分28粒子輸運問題(7/7)√輸入N,

H,E00,

Emin,g00光子i的初態(tài)m=0,Z=0,

E=E00,g=g00利用E計算截面Ss,St利用Ss抽樣自由程S計算碰撞點位置ZZ+Sg利用Ss,St抽樣作用類型計數(shù)器N1N1+1利用E0

E抽樣散射能量E利用g0g,

E0

,

E抽樣散射方向gmm+1ii+1輸出P=N1/N,

s,

DPZ

0

?類型=吸收

?E

Emin

?YNNYNYm>M?NYi>N?NYZ

H

?NY粒子輸運問題(7/7)√輸入N,H,光子i的初態(tài)m=029隨機過程模擬(1/3)蒙特卡羅模擬隨機過程的兩種情況已知隨機過程的概率分布函數(shù)隨機過程的統(tǒng)計特征建立與隨機過程的模型，形成隨機量，并使其數(shù)字特征(概率、平均值、方差等)是問題的解由已知的概率分布函數(shù)進行大量的抽樣統(tǒng)計處理抽樣結(jié)果，給出問題的解及其誤差已知隨機現(xiàn)象的觀測數(shù)據(jù)概率分布函數(shù)分析現(xiàn)象的特征，假設(shè)隨機量服從某種分布，建立概率模型由觀測數(shù)據(jù)推斷分布中的參數(shù)按推斷的分布進行大量的抽樣比較抽樣值和觀測值，根據(jù)誤差，判斷假設(shè)的準確性√隨機過程模擬(1/3)蒙特卡羅模擬隨機過程的兩種情況√30隨機過程模擬(2/3)例：a

粒子衰變的蒙特卡羅模擬隨機現(xiàn)象的觀測數(shù)據(jù)概率分布函數(shù)觀測數(shù)據(jù)：每隔

Dt<<半衰期

測量一次放射的

a

粒子數(shù)，共測量

N=2608

次，測得

k

個粒子的次數(shù)為

nkknkknk05762731203