第二章隨機(jī)數(shù)第四 方法解粒子輸運(yùn)問

第一 方法概第一 方法概，并首先在核的試驗(yàn)與研制中得到了應(yīng)用。蒙特很大區(qū)別。它是以概率統(tǒng)計理論為基礎(chǔ)的法。由于方法能夠比較真地描述事物的特點(diǎn) 電子計算機(jī)的發(fā)明，方法作為一種獨(dú)立的方法被提出來，并首先在核的試驗(yàn)與研制中得到了例1. 例2.射擊問題（打靶游戲例 P

2l2l 思例2.射擊問題（打靶游戲 g0g(r)f g現(xiàn)假設(shè)該運(yùn)動員進(jìn)行了Ｎ次射擊，每次射擊的彈著點(diǎn)依次為r1r2rNg(r1，N 1

g(ri由以上兩個例子可以看出，當(dāng)所求問題的解是某個的概率，或者是某個隨量的數(shù)學(xué)期望，或者是與概率、數(shù)學(xué)期望有關(guān)的量時，通過某種試驗(yàn)的方法，得出該發(fā)生的頻率，或者該隨量若干個具體觀察值的算術(shù)平均值，通過它得到問題的解。這就是方法的基本思想。當(dāng)隨量的取值僅為1或0時，它的數(shù)學(xué)期望就是某個的概率。或者說，某種的概率也是隨因此，可以通俗地說，方法是用隨機(jī)試種分布密度函數(shù)f(r)的隨量ｇ(r)的數(shù)學(xué)期望g0g(r)f率語言來說，從分布密度函數(shù)f(r)中抽?。蝹€子樣r1，r2，…，rN，），將相應(yīng)的Ｎ個隨量的值g(r1)，N 1

g(ri 計算機(jī)模擬試驗(yàn)過程，就是將試驗(yàn)過程（如投針，射擊）化為數(shù)學(xué)問題，在計算機(jī)上實(shí)現(xiàn)。以上述例 例2.射擊問題（打靶游戲由上面兩個例題看出，方法常以一個“概布抽樣的方法，得到部分試驗(yàn)結(jié)果的觀察值，求得問xl

勻 勻

1a, 0xaf1x)0, 1/ 0f2() xs(x,)

當(dāng)xls

N1s(x,N NPs(x,)f1(x)f2(

dx

2l 每次投針試驗(yàn)，實(shí)際上變成在計算機(jī)上從兩個均勻分布的隨量中抽樣得到（x,θ），然后定義描述針與平行線相交狀況的隨量每次投針試驗(yàn)，實(shí)際上變成在計算機(jī)上從兩個均勻分布的隨量中抽樣得到（x,θ），然后定義描述針與平行線相交狀況的隨量s(x,θ)，為789

N動動

1

g(ri N NNi 1NiN（E(X)<∞），PlimX

E(X)N 樣的算術(shù)平均值XN 方法的近似值與真值的誤差問題，概率論的中心極限定理給出了答案。該定理 ，如果隨機(jī)零的方差σ2，即02(xE(X))2f(x)dxN N

t2/

XNE(X

x N

t2/XNE(X)

dt1 N 這表明，不等式 E(X)近似地以概 1－α成立，且誤差收斂速度的階為O(N1/2)。 N上式中與置信度α是一一對應(yīng)的，根據(jù)問題的確定出。α31NN 1NN 21NN2iXi顯然，當(dāng)給定置信度α后，誤差ε由σ和N決定。要減小ε，或者是增大N，或者是減小方差σ2。在σ固定的情況下，要把精度提高一個數(shù)量級，試驗(yàn)次數(shù)N需增加兩個數(shù)量級。因此，單純增大N不是一個一般來說，降低方差的技巧，往往會使觀察一個子樣的時間增加。在固定時間內(nèi)，使觀察的樣本數(shù)減少。所以， 法的優(yōu)劣，需要由方差和觀察一個子樣的費(fèi)用（使用計算機(jī)的時間）是 方法中效率的概念。它定義為2c是觀察一個子樣的平均費(fèi)用。顯然2c 越小，方法有效。能夠比較真地描述具有隨

能夠比較真地描述具有隨機(jī)性質(zhì) 在計算 g g(x, sx d(x(i),x(i) ,(i

g(x(i),x(i) ,x )

i

卡羅方法，不會有原則上的。卡羅方法的收斂速度為O(N12)，與問題本身的維數(shù)間的變化，不影響誤差。也就是說，使用方特點(diǎn)，決定了方法對問題的適應(yīng)性。而般數(shù)值方法計算積分時難以克服的問題。方案的計算量相當(dāng)。例如，對于層為均勻介質(zhì)的平板幾何，要計算若干種厚度的概率時，只需計算最厚的一種情況，其他厚度的概率在計算最厚另外，使用方法還可以同時得到若干個差并不是一件容易的事情，而方法則不然。根據(jù)方法的誤差，可以在計算所求量的同時計算出誤差。對干很復(fù)雜的方法計算問解決這一問題有時相當(dāng)，方法則不存在

O(N1/2 由于方法的誤差是在一定置信水平下估概率的計算問題，計算結(jié)果往往比真值偏低。而 方法所特有的優(yōu)點(diǎn)，使得它的應(yīng)用范圍越來越廣。它的主要應(yīng)用范圍包括：粒子輸運(yùn)問題，統(tǒng)計物理，典型數(shù)學(xué)問題，真空技術(shù)，激光技術(shù)以及醫(yī)學(xué)，生物，探礦等方面。隨著科學(xué)技術(shù)的發(fā)展，其方法在粒子輸運(yùn)問題中的應(yīng)用范圍主要 第二章 f(x)

0x

F(x)

x0xx 上的點(diǎn)在

0ai1，i1,2, ssP(niai i1 ,s) 也難以滿足方法對隨機(jī)數(shù)需要量非常大的要12122 m能進(jìn)行程序復(fù)算，給驗(yàn)證結(jié)果帶來很大。而且， nk

T(n,n1 ,nk

n確定ξ，ｎ=1，2，…。經(jīng)常使用的是k=1的情況，其遞推為：nk

T(n遞推和初始值ξ1,ξ2…，ξk確定后，整個隨機(jī)數(shù)由于隨機(jī)數(shù)序列是由遞推確定的，而在計算機(jī)上出現(xiàn)這樣的n＇，n(n＇<n)，使得下面等式成

i1,2, k隨機(jī)數(shù)序列便出現(xiàn)了周期性的循環(huán)現(xiàn)象。對于k=1的要遞推選得比較好，隨機(jī)數(shù)間的相互獨(dú)立性是可為用方法解任何具體問題時，所使用的隨機(jī)發(fā)生周期性循環(huán)現(xiàn)象的偽隨機(jī)數(shù)的個數(shù)稱為偽隨機(jī)數(shù)的周期。對于前面介紹的情況，偽隨機(jī)數(shù)的周期從偽隨機(jī)數(shù)序列的初始值開始，到出現(xiàn)循環(huán)現(xiàn)象為止，所產(chǎn)生的偽隨機(jī)數(shù)的個數(shù)稱為偽隨機(jī)數(shù)的最大容量。前面的例子中，偽隨機(jī)數(shù)的最大容量為n″。乘同余方法是由Lehmer在1951年提出來的，它的遞推確定：xi1axi (modM

xi1M

iMP0 P0=2P1，…(M){(P

),(P1 ( 0(P0

2)2

當(dāng)00或當(dāng)0當(dāng)0(Pi1)Pi

i1,2, ra

(mod 當(dāng)0(mod 當(dāng)0 當(dāng) i對于0n(Piiian (modPiix1a=其中k為使52k+1在計算機(jī)上所能容納的最大整數(shù)，即a為計算機(jī)上所能容納的5的最大奇次冪。一般地，容量λ(M)=2s-2。乘同余方法是使用的最多、最廣的方法，在計算機(jī)產(chǎn)生偽隨機(jī)數(shù)的乘加同余方法是由Rotenberg于1960年提出來的，由于這個方法有很多優(yōu)點(diǎn)，已成為隨機(jī)數(shù)序列由下面遞推確定：xi1axi

xi1M

i1,2,關(guān)于乘加同余方法的最大容量問題，有如下結(jié)論：如果對于正整數(shù)M的所有素數(shù)因子P，下式均成a (moda (modM=a=2b+ c=這樣在計算中可以使用移位和指令加法，提高計算樣時的均勻性問題。其中n為偽隨機(jī)數(shù)序列的最大容對于任意的0≤x≤1Nn(x表示偽隨機(jī)數(shù)序列

(n)sup|Nn(x)x (n)

|i|,| imax n (n) 1當(dāng)(n)2n

2i

i1,2, n(n)(n)n對于任意0x,y ，令Nn(x,y)表示 x,i1 (n)sup|Nn(x,y)Nn(x)Nn(y)

則ε(n)標(biāo)志偽隨機(jī)數(shù)序列ξ1,ξ2…，ξn的獨(dú)立程度，(n)

nnn(n)

nnnmax(1,2)

第三章作業(yè)中抽取簡樣。隨機(jī)數(shù)序列是由單位均勻分布的總體中抽取的簡樣，屬于一種特殊的由已知分布的簡樣，是在假設(shè)隨機(jī)數(shù)為已知量的前提下，使用簡樣的。對于連續(xù)型分布，常用分布密度函函數(shù)f(x)產(chǎn)生的簡樣的。另外，在抽樣過程中對于任意給定的分布函數(shù)F(x)，直接抽樣方法如Xn n1,2, ,F(t)XFinf FX(x)P(Xnx)P( t P(nF(x))F對于任意的n成立，因此 量序列X1，，ξ2，…，ξ是相互獨(dú)立的，而直接抽樣 Measuretheory,N.Y.VonNosrtand,1950]§45定理2）。F(x)xi其中x1x2，為離散型分布函數(shù)的跳躍點(diǎn)，P1，P2，…為相應(yīng)的概率，根據(jù)前述直接抽樣法，有 XFxI 當(dāng)PiPi 例1P(xn)

CnPn(1P)NNN XF 當(dāng)PiPi 例2.泊松(Possion)P(xn)

eXF

當(dāng)

ei＝0 i＝0例3.擲 P(Xn)6n1 則XFXF[6]例4iiPit

i1,2, n Pi

Pi 例5依次為：

P PP

對于連續(xù)型分布，如果分布函數(shù) 的反函XFF1例6.在［a，b］bF(x)b

當(dāng)x當(dāng)ax當(dāng)x則XFa(ba)例7βf(x) 0xxF(x)x

f(t)dt

x2tdtx2 0x0則XF例8f(x)aeax xxxF(x)

f(t)dt

x0 1ln(1 1 反函數(shù)存在且容易實(shí)現(xiàn)的情況，使用起來是很方便的。但是對于以下幾種情況，直接抽樣法是不合適分布函數(shù)可以給出其解析形式，但是反函數(shù)給不出 為了實(shí)現(xiàn)從己知分布密度函數(shù)f(x)抽樣，選取與M

f(x)xff(Xh Mh(XhXfXh

f(xhMh(xh

P(x

xdx)

x

f(Xh P

Mh(Xh) f(Xh PxXhxdx,Mh(X

f(Xh Mh(Xh)

f(Xh Mh(Xh)h(Xh)dXh f(XhMh(Xh

h(Xh)dX

x

f(XhMh(X

h(Xh)dX

f(XhMh(Xh

h(Xh)dXx

f(Xh

ff(Xh)dXh時要使得h(x)容易抽樣且M的值要盡量小。因?yàn)镸提高抽樣效率。抽樣效率是指在挑選抽樣方法中進(jìn)行

f(Xh Mh(Xh)

f(XhMh(Xh

h(Xh)dXh1M1f(x在［0，1］h(x)=1，Xh=ξ，Msupf(x) MXf例9f(r)

當(dāng)0r

F(r)

RrR20rfR0由于開方運(yùn)算在計算機(jī)上很費(fèi)時間，該方法不是 h(r)1，f(r)2r，M2，rR ≤rfR0取rfR01就可以了，亦即rfR0max(1,2布F(r)r。

與max(1,2在實(shí)際問題中，經(jīng)常有這樣的隨量，它服從分布的隨量，稱這樣的隨量服從復(fù)合分布。f(x)Pnfn F(y) f(x)f2(xy)dF1(其中f2(x/y表示與參數(shù)y有關(guān)的條件分布密度函數(shù)，數(shù)f2(x/YF1)中抽樣確定Xf2(x/YF)f Xf(x f

p(xX xdx)p(xXf(x/ )x f2(xY)dxdF1(Y)f(x 例10 eEn(x)

當(dāng)xynnyf1(y)yef2(xy)

當(dāng)y1當(dāng)x則En(x) f2(xy)f1(nf nf1

)

,

max(

f(x)H(x,y)f2(xy)dF1( H(Xf(x),YF> H(Xf(x),YF>Xf

f2(x/YF1

P(x

xdx)

x

H(Xf,YF)2 P H(Xf,YF)MPxX xdx, 1M2 H(Xf,YF)P 1

H(x,y

f2(xy)dxdF1( H(x,y

f2(xy)dxdF1(

H(x,y)f(xy)dxdF( H(x,y)f(xy)dxdF(

H(x,y)f2(xy)dF1(y)dxf 量y的抽樣，將其表 量x1，x2，…，xn的函數(shù)yg(x1,x2 ,xnYfg(X1,X2 ,Xn例111f(x) 1f(x) sinsincoscosxyx2x2x2x2(x,y)(x,y)在上半個單coscos2cos2sin2

x2x2sinsin22sincos

x2(x,y)的問題。x y

12

222

2 cos 2,sin 122112211

2

2為實(shí)現(xiàn)散射方位角余弦分布抽樣，最重要的是在上半個單位圓內(nèi)產(chǎn)生均勻分布點(diǎn)。下面這種方法，首先在單位圓的半個外切正六邊形內(nèi)產(chǎn)生均勻分布點(diǎn)，1 311 11,21132

222

1212cos 2,sin 1212

32 323E 3例12f(x)

ex2 f(x,y)1e(x2y2)

x ysinf(,)e2eef1()2f() 22XfYf

22~N(μ,σ2)的抽樣，~~XfX~YfYff(x)

(k1)!(n

xk1(1

0x一組簡單的相互獨(dú)立隨量的函數(shù)，通過這些簡單隨量的抽樣，實(shí)現(xiàn)βx1x2，kRk(x1,x2 ,xnF(x)=x

(x)k

CCnnnn

(x)k

Ci

xi

不難驗(yàn)證，ζk的分布密度函數(shù)為βXfRk(1,2 ,nf(x)Pnfn 首先抽樣確定n’，然后由fn’(x)n－1 XfXf

當(dāng)PnPn 例14if(x)axii f(x)

i

(i1)xi Pf Xf ,n

當(dāng)PiPi 例15 f(r)

當(dāng)R0r

r3R3F(r) R3 r(R1R0)xf(x)

3R0

2x

0R2RR ≤1

00≤ 1≤

>Xf

Xfmax(2,3

Xfmax(2,3,4rf(R1R0)Xff(x)A1f1(x)A2f2f(x)

f2(x)f1(x)

2f(x) Af2(Xf Af2(Xf 1 f1(Xf1E

1XfX11

f f(x)f2(x)A1 A2 f2

f1(Xf

2 2 f2(Xf E

2XfX2 A1mA2例16βf(x)2(1 0x或11 ≤Xf

21 ≤Xf≤XfXfmin(1,2f(x)H(x) H H(Xf M E1XfX 1例17f(x) 1 1x f aelnf x1fi(x)i(a1i x 1xiXf[(a1i1)iif(x)H(x)fi

i(a1i i(a1i

H(x)

M

x1i1[(a1i1)21]1Xf[(a1

≤2≤

E

lnai(a1i1) 23f(x) x

x,xf1(x)

e2x x33 3 H(x)31H(x)31x3,xM則12

≤ 3

E

f(x)Hn(x)fnf(x)H1(x)f1(x)H2(x)f2f(x)

H1

f1(x)

f2 Pf*(x)Pf 1H1(x)f1(x)dxP2H2(x)f2(x)dx1f*(x)1

H1

f112f*(x)2

f2≤>1≤> 2

H( f11f1M1

2 H2(H2(X2) 212XfX XfX12P

E11 1P2 2 2 M M 的。下面的方法可以不用計算P1：對于任意小于1的正數(shù)P1，令P2=1－P1H1(x) H(x,y)

1y

f

y1y

Mmax(M1,M2 2f(xy) 2f2

y

E 1, 2M1 M P

P M M1M

M12E2

M1M≤≤1 H1(Xf1)

M1

M1M1M>H2(Xf2) M212XfX XfX12 E1E2 2 M M1MM(MM)P2M

M)P2M

M2

M2P2MM(P2P2)M M1M2(M1M2M2

M2P22MMPP (M21M1P2 例19令光子散射前后的能量分別為和（m0c2為單位，m0為電子靜止質(zhì)量，c為光速），x，x的分布密度函數(shù)為：f(x)

1x

1

1

1x1K

x3 K()2(1)ln(12)14

f(x)

(x1)2

f(x)12 f(x) 1x2H1(x)K()(12)

(xH2(x)K f(x)H1(x)f1(x)H2(x)f2(x)X

11X 12H(x) M H(x) M 27K ≤≤X112>11f1 >X 1227(X >3

2

X4fX42 XfX1

XfX2E 27(12)KM1M f(x)H1(x)f1(x)H2(x)f2其中H1(x、H2(x)為非負(fù)函數(shù)，f1(x)、f2(x為任意分布f(xf(x)f1(x)H1(x)

H2(x)f2(x) H(x)f(x 令H(x)的上界為

H2(x)f2

H1(x)f1(x)H(H()H2(Xf)f2(Xf)11 M(1

1

f( 1XfX1

f(xf(x)

(x)H

(x)

H1(x)

H(x)f H1(XH1(Xf)f1(Xf H(X) M(1m)

f(

2

2XfX2例20f(E)CeEA

EminE其中A，B，C，Emin，Emaxf1(E)

(E)

eE E

H(E)

H2(E)

exp f(E)H1(E)f1(E)H2(E)f2 2exp41A1 11

AB

8

m＝011

Emin1 ln 1 1(ln21

≤

22

EfEf1E

ef(x)f1(x)H件，H(x)為任意奇函數(shù)，即對任意x滿足：f1(x)f1(x)H(x)H≤H(Xf11f1(Xf1XfX1

XfX122 1 H(Xf)

1 H(Xf) x)P x, 1 1P

x,

2 f(

2 f(

1 H(Xf)

1 H(Xf)P x, 1 1P

x,

1 1

2 f(

2 f(

x1

H(Xf)

1

H(Xf)

1 1f(

)dXf

1 1f(

)dX 2

f( )

x2

f( ) x 1f( )H( x

1f( )H(

x

x 1f( )H( x1f( )H( x

xx f( )H( xx

f( )

.f()1 1 CcosC2 LcosL

11A22 1 A2122 依照質(zhì)心系散射各向同性的假定，可得到系散射角余弦μL的分布如下：1A212 f(L)

L2L

1LL2AL

A21 該密度函數(shù)中的第一項為偶函數(shù)，第二項為奇函數(shù)，1f1(L)1

A21 L L

1L2

1A21 A22 1A21 1≤L

A212 LA21L P(a)P(2a2) A21 A2a 1A21 1A2

2

(A212211A 111A 另一方面，在2a2即1 a1的條件下，η/a 2,A21 2于是，得到質(zhì)心系各向同性散射角余弦分布的抽12A21 (A2122 12A 2A

1A211≤ A21 L

2

Lf(x)

f0(x,

f0(x,YH(X) XfXfP(x

xdx)Px

xdx H( ff

x H( f

PYfH(Xf

f0(x,

f

f0(x,例22.各向同性散射方向的抽樣為了確定各向同性散射方向uivjwk usinvsinsinw1cos211cos21 x

A22 y A22 23可以證明，當(dāng)221時， 23

xy

f(x,y)1x

1 2y1 2y1

2

1

1x(1x)12

f(x, f(x) (1x)12

f(x, A2條件2

33333311

2)2 2≤2usincosvsinsin

2 323222

2A223wcos 321221

E12A2 Ig(x)fI

g(x)f(x)f*(x)dxf*(x)

g*(x)f*(x)dxN是可以這樣計算積分I： 1g*(X)

f(Xi

g(X)

)g(X Nf*(X

N

f(x)抽樣，現(xiàn)改為由另一f*(x)抽樣，并附帶一個權(quán)重糾偏因子W(x)f f*(x)f(xXf，滿足P(xXfddx)f(x)dxff

*)P(x

*ddx)W(x)

*(x)dxff一般情況下，Xf是具有分布f(x)總體的簡 ，只代表一個。Xf*是具有分布f*(x)總體的簡 ，但不代表一個，而是代表W(Xf*)個，ffa(xf(x)fa(xf(x的一個近似分布密fa(x)

f(x)dxfi*, 當(dāng)xi1xX xi1

Fa

(x

當(dāng)F

)

(x F(x)

或

i1

f

當(dāng)f

ff fi

ii f(x)

xx

f f

i1x 其中，cfi＝f(xi，x0，x1，，xn為任Fa(xi11Fa(xi)F

≤

>X xi1(xixi1)a

X xi1(xixi1)(12a理，以及用具有布的隨量。例23.正態(tài)分布的近似抽樣

ni1in n n16n＝12時，有6X122i2i1m fa(

ffaf(x)mfa(x)H1(x)f1(x)H1(x為非負(fù)函數(shù)，f1(x為一分布密度函數(shù)。 ≤1 H1(Xf ≤XfX XfX 抽樣效率Em(1 1由上述近似-修正抽樣方法可以看出，如果近似分fa(x)選得好，m1，這時有很大可能fa(x)Xfa，而只有很少的情況需要計算與f(x)有關(guān)的函數(shù)H1(Xf1)。在乘抽樣方法中，每一次都H(Xfa)＝f(Xfa)／fa(Xfa)f(x)比較復(fù)雜例24裂變中子譜分布的近似-f(E)CeEA

EminE其中A，B，C，Emin，Emax f(E)EexpEA A f(E)1Aexp1A

H(E)

1A

A2 f(E)mfa(E)H1(E)f1(E)2m fa(E

f(E)fa

fa(E

ACE

Ef Aln(23

Ef 1

a1a

1m

＞ 2

H1(Xf1>1>Xf

X

A1

對于鈾-235，m≈0.8746，M≈0.2678，λ≈0.5543，抽樣效率E≈0.9333。而且近似修正抽樣對數(shù)。因此，較之乘減抽樣方法了計算時 f(x,y)f1(x)f2(yf1(x)f(x,2f(yx)f(x,y)2

f(x,

f(x,對 分布密度函數(shù)，也可直接采用類似于一f(x,y)H(x,y)f1(x,H(x,y)為非負(fù)函數(shù)，f1(x,y)為任意二維分布密度函MH(x,y)的上界，則有二維分布的乘抽樣方H(H(Xf,Yf11>M≤XfXf1 YfYf1例25ef(x,y) x1,yxf(x,yf(x,y)f1(x)f2(y f(x) f2(yx)xXf11 Xf1 f(x)ex xXfln

ni>i>≤NY

ii

nnXfn

ni0，YYY>0>≤NY

ii

nnXfn f(x)

1

1

2

2K

f(x)12 f(x) 1x2H1(x)K()(12)

(xH2(x)K >≤第四 方法解粒子輸運(yùn)問 作業(yè)第四 方法解輻 問輻射（光子和中子）問題是理直觀出發(fā)，說明方法解決這類粒子巧對于諸如輻射、多次散射和通量計算等常要用一些吸收材料做成物擋住光子或中子。我們所關(guān)心的是經(jīng)過后射線的強(qiáng)度及其能量分布，這就是問題。當(dāng)物的形狀復(fù)雜，散射各向異性，材,核反應(yīng)截面與能量、位置有關(guān)時，難以用數(shù)值方法求解，用方法能 模型，在厚度為能量、方向分布的輻射源S。求粒子 rE和運(yùn)動方向ΩS＝rE,Ω表示。的權(quán)重W，這時S'＝rE,Ωt,W。S＝zEZS＝rE或

(rm,Em,?m,tm,Wm它表示一個由源發(fā)出的粒子，在介質(zhì)中經(jīng)過m次碰撞一個由源發(fā)出的粒子在介質(zhì)中運(yùn)動，經(jīng)過若干次碰撞后，直到其運(yùn)動歷史結(jié)束（如逃出系統(tǒng)或被吸收等）。假定粒子在兩次碰撞之間按直線運(yùn)動，其運(yùn)動方向與能量均不改變，則粒子在介質(zhì)中的運(yùn)動過程可S0，S1，…，SM-1

rM M M

1

MS0SMM這樣的序列稱為粒子隨機(jī)運(yùn)動的歷史，模擬一個 說明，這時狀態(tài)參數(shù)取S＝(z,E,cosα)。 f(z0,E0,cos0)f1(z0)f2(E0)f3(cos0

f1(z0)(z0已知狀態(tài)Sm-1，要確定狀態(tài)Sm，首先要確定下一個碰撞點(diǎn)的位置zm。在相鄰兩次碰撞之間，中子的輸運(yùn)長度l服從如下分布：f(l) l , )l ,

對于平板模型，lf(l) l , )

l l

,

ll ,

flll0t(rm1l?m1Em1)dl l

l zmzm1l

如果zm≥a，則中子 ，若zm≤0,則中子被反 通常介質(zhì)由幾種原子核組成，中子與核碰撞時，要確定與哪一種核碰撞。設(shè)介質(zhì)由A、B、C三種原子核組成，其核密度分別為NA、NB、NC，則介質(zhì)的宏觀 )A 其中ABC分別為核A、B、C () )N() )、N )、N、

)分別表示(·)核的宏觀總截 AAB AB

CP C

t(Em1 >PA >與C

撞類型。中子與核的反應(yīng)類型有彈性散射、非彈性散射、(n,2n)反應(yīng)，裂變和俘獲等，它們的微觀截面分別為B

(n,2n)

)B

(n,2n)

)B(E

m )B

)B(E(n,2n)

(n,2n)

m )B(E

m )B(E

m 問題中，中子與核反應(yīng)常只有彈性散射和 )B

)B

B Pel，則為彈性散射；否則為吸收，發(fā)生吸收反如果中子被碰撞核吸收，則其輸運(yùn)歷史結(jié)束。如果發(fā)生彈性散射，需要確定散射后中子的能量和運(yùn)動方向。中子能量Em為：Em

C2C

r(AA C fC(μC)中抽樣產(chǎn)生 系散射角θL的余弦

1C 1A22C如果給出 系散射角余弦分布fL(μL)，可直接從fL(μL)中抽取μL，此時能量Em與μL的關(guān)系式

A212 (A cosmcosm1cosLsinm1sinLcos至此，由Sm-1完全可以確定Sm。因此，當(dāng)中子由源出發(fā)后，即S0確定后，重復(fù)步驟(2)～(5)游動歷史終止。于是得到了一個中子的隨機(jī)游動歷史S0，S1，…，SM-1，SM z0

zMEM

cosM,,以上模擬過程可分為兩大步：第一步確定粒子的步又分為兩個過程：第一個過程是確定碰撞點(diǎn)位置 能量及運(yùn)動方向，稱為碰撞過程。對于中子而言，碰撞過程是先確定散射角，進(jìn)而確定能量和運(yùn)動方向；而對于光子，碰撞過程是先確定能量，再確定散射角以及運(yùn)動方向。重復(fù)這兩個過程，直至粒子的歷史終這種模擬過程，是解任何類型的粒子輸運(yùn)問題所 算的物理量進(jìn)行估計。對于問題，我們要計算中子的率。每個中子的隨機(jī)游動歷史，它可能

了NN個中子， N1 ?(1)

1N N

N?(1)N?(1)

2P(1

?(1)?(1) 為得到中子 的能量、角分布，將能量、Emin

E1E0001 J 蔽的中子按其能量、方向分間隔記錄。設(shè)一的中子能量為EM，其運(yùn)動方向與Z軸夾角為αM量EM屬于第i個能量間隔ΔEi，角度αM屬于第j個角度間隔Δαjij個角度計數(shù)器中加"1"。N?(1)

i1,2, ?N2,

IN N2,

j1,2, J分別為中子的能量分布和角分布。其中N和N2,i分別為第i個能量和第j個角度間隔的中?(1)*

?

i1,2,

N

?(1)*

?N2,

N2,

j1,2, N2,

?N

N1 從模擬物理過程來說，直接模擬法是最簡單、也 都對沒有貢獻(xiàn)。而在許多問題中，率的數(shù)量級在10-6到10-8。進(jìn)一步，如果我們要求率達(dá)|?P| P(1P)

N N N那么，N要大到驚人的數(shù)量級 的

B )B )散射，一部分吸收。彈性散射這部分繼續(xù)；吸收部分則停止。也就是說，我們利用中子權(quán)重的變化來反應(yīng)繼續(xù)彈性散射的部分。這就是簡單法的顯然， 子的權(quán)重W已成為中子狀

zM

EM

cos

cos

M,

cosM Wm

B因子 B

)B ) 當(dāng)zM 如果我們 了N個中子， 率P的無偏估計 ?(2) 1 N N

n類似地，可以得到 中子的能量分布和角分布。只不過在對各計數(shù)器進(jìn)行的加"1"操作改為加WM。

1

(?(2) NNNNNN

1

(2 nnnn

N 法的思想在方法中用途很廣泛。例得用直接模擬法模擬每一個中子是非常的。這種情況可以利用法來處理。 1el1

2(n,2n)3t

f

：3

f

t1el1

2(n,2n)3

f法的思想，還可以應(yīng)用到連續(xù)分布情況和偏法雖然改進(jìn)了直接模擬法，但它同樣只關(guān)心中子是否這一信息，因此對每一個中子歷史一個中子，可能在介質(zhì)內(nèi)不發(fā)生碰撞而直接，也可能在介質(zhì)內(nèi)發(fā)生一次碰撞后 率P可表示為：P

S0

其中Pmm次碰撞而的概率。這表明，可以用求Pm(m=0,1,…)的方法得到P。這樣，中子對概率的貢獻(xiàn)就不只限于末次碰撞了。

zM

EM

cos

cos

M,

cosM m過m次碰撞后 的部分?m

m0,1, ,M顯然，具有初態(tài)S0＝(0,E0,cosα0,W0)的中子， ?Wexp(E 0 0

cosm次碰撞后的中子具有狀態(tài)Sm＝(zm,Em,cosαm,Wm)，其可能的部分，正好m次碰撞的部分： (E

a

，

m m0,1 ,M這里的這種估計技巧，由于是對每次碰撞后的狀態(tài)，求其后碰撞直接的貢獻(xiàn)，因此該方法也 PP m如果我們 了N個中子， ?(3)N

N

NN

?(n)(?) NN

?(3) N P(1P) 其中

N N N P*K*

N NN*K

N*K這意味著，達(dá)到同樣的相對誤差，粒子的數(shù)KK倍的計算量。指數(shù)變換法就 *(E)(E)C C EminE

散射截面仍為Σel(E)Emin、Emax分別為能量的下限和上限，α為粒子的運(yùn)動方向與Z軸的夾角?？梢宰C明這個偽過程的概率P*與原過程的概率P之間有如下關(guān)系:P*eCa顯然，KeCa1e指,可以明顯看出P*增大的原cosα=1Z軸方向一致，其截面最小，粒子沿Z軸方向輸運(yùn)的距離較遠(yuǎn)；而當(dāng)cosα=－1Z軸方向，這時其截面最大，粒子向后輸運(yùn)的距離較短。因此，截面變換的結(jié)果是加強(qiáng)了粒子向前運(yùn)動的能力，因而使概偽過程的構(gòu)造與幾何形狀及所考慮的問題有關(guān)。,所構(gòu)造的宏觀總截面與平板的情況不同，粒子的模衡量一種技巧的好壞，除了看其方差大Ef（方差與費(fèi)用乘積的倒數(shù)）全面考慮： 2其中σ2為方差，T為所需費(fèi)用。Ef大時，所用方在一般情況下，有些方法雖然減小了方差，卻增 f(x)

1

1

2

2

xf(x)12 f(x) 1x2H1(x)K()(12)

(xH2(x)

K ≤第五 作業(yè)第五 了方法在計算機(jī)上得以實(shí)現(xiàn)的基礎(chǔ)。在計算機(jī)上使用方法解粒子輸運(yùn)問題析結(jié)果過程。r的分布密度函數(shù)為：f(r)

當(dāng)0rr

rR0max(1,2r的分布密f(r)

rr

≤

R1R0 R1R1r(R1R0)max(2,3)

r(R1R0)2r的分布密度函數(shù)為：33f(r)

當(dāng)0rrrRmax(1,2,321221

3≤3

1 x0R

y0R2

z0Rr的分布密度函數(shù)為：f(r)

當(dāng)R0r0r0

≤

> R2RR 1 ≤1

0R0R1R12>R12x

xmax(2,3r(R1R0)x

xmax(2,3,4x0rsincosy0rsinsinz0r(2A22A22)2 ≤x0rsincosr

22A221rsinsinr

2 2A22 rcosr

2A22 原點(diǎn)，圓柱的軸向?yàn)閦軸，則分布密度函數(shù)為： 當(dāng)x2y2R2|z|f(x,y,z)

2 21221

1 x0R

≤y0R2

z0H

(x*,y*,z*)

f(x,y,z)(xx*)(yy*)(zz*) x y z 0f(r)(rr*)00r*00 r0球外平行束源分布是指粒子平行入射到半徑為Rr其分布f(r)(rrr21221

1 y0R1x0

R2y2R2y2z00其分布密度函數(shù)為：f(E)(EE0Ef(E)CeEA

EminE其中A，B，C，Emin，Emax 1m

＞ 2

(E) EAln(

EE

1 H(E) 1A

A2此時，用數(shù)值曲線抽樣方法抽取E。f(E)

23

EeE

,E 2>1e2ln22>≤E3 fΩ)f1()f2()f()1 f() usincosvsinsinwcosA2) 21(2A2) 21≤usincos

2A1A2vsinsin

2A1232A223

A22A22A2wcos

2312A2231

A2 f() 1umax(1,2 DDcos* 先抽樣確定，再轉(zhuǎn)換成θ。cosR

OSz R2D R2D2sin2vwcos成的次級粒子進(jìn)行。這種方法比較簡單、直觀。粒子由狀態(tài)Sm到狀態(tài)Sm+1時，需要確定粒子的空間位置rm+1，能量Em+1和rm1即

Lxm1xmLym1ymzm1zm

其中(umvmwm?m的方向余弦，L為兩次碰撞點(diǎn)間Lf(L)t(rm1,Em) l?,E LL f()e中抽取

L l?,EL

L ln Lit,i(Em)Lit,i(Em L

Lit,i(Em t,I(Emi個介質(zhì)區(qū)域內(nèi)走過的距離，t,i(Em)i個介質(zhì)(i=1，2，…，Imax) Lit,i(Em在求Li(i=1，2，…，Imax)時，如果系統(tǒng)形狀復(fù)LL

>>t(rmL1?m,Em2≤

L其中t,maxE)maxt(rr e

當(dāng)0f()1 ln[1(1eDm)L。此時，粒子的權(quán)重需乘以糾偏因子(1eDmEm1

A212LC2A212LC2Em1

其中，A為碰撞核的質(zhì)量，CcosC

L 1 1 Em1

Em2

C A1 1K 1K KE1mL

2

f(x)

1x(4)光(4)光

1

1

1x1K

x3K()2(1)ln(12)14

xm0c2為單位，m0為電子靜止質(zhì)量，c x1

12 1211x2 1

32 3

x2127 223>42x

x mc2

mc2

法為C21。質(zhì)心系散射角θC抽樣確定后， 1C 1A22CA為碰撞核質(zhì)量，CcosC,LcosL 1A1K Em 1A2 E2A1

2l1fP

當(dāng)|xf(x) P0(x)1(x)nfa(x)Φk(xxknnn n

2l1 2l

Pl(xk

[P(x 。

nf*(x)Ψ(xx

|Φ nkΨ , w k|Φnk|Φk

|Φk

xxK

當(dāng)Ψk f(L)(111L L11

bcwmummmmmmmm

wm1

mmacosLccos

bsinLdsin

1a2方位角在[0,2π 當(dāng)u2v20時，不能使用上述 um1bcvm1bdcosm1cosmcosLsinmsinLcos )

sinLsin

)

cosLcosm

徑向位置rm+1的計算 rm1

2L22Lrr r余弦來描述。使用球面三角，粒子碰撞后瞬時運(yùn)cosθm+1的計算為：cosm1cosmcosLsinmsinLcos

rm

sincos

r

[0,2π]上均勻分布的方位角，m為在rm+1rm+1徑向之間的夾 子輸運(yùn)的各種實(shí)際問題時,所遇到的基本幾何問題。點(diǎn)r0x0y0z0沿方向?0u0v0w0 rr0l dl

by0cz0au0bv0 與平面平行時abvcw rc＝(xcyczc，R，則球面方程為：c(xx)2c

(yyc

(zzc

R2l：l(r0rc)?0(x0xc)u0(y0yc)v0(z0zc0cc02R20cc0R2|rr|2 x)2(

y)2

zll球面方程為r＝R。R2(rsin200R2(rsin200dr0sin0當(dāng)cosθ0≥0R2(rsin200R2(rsin200c(xx)2c

(yyc

R2(xcyc0)為圓柱的中心，R為圓柱底半徑。r0沿Ω0l為：l 010(x0xc)u0(y0yc2(R2R2)(1w2 R2 x)2(yy 0R0≤R時，r0w0l 010

000R0＞R時，r00若δ＜0，Δ≥0w21，l0l 010w201w2x2y2c2 l 01(1c20

yvc2z 2(x2y2c2z2)[1(1c2)w20 0 x2y2 l

(1c2

1/12的ΔOAB來做代表；正方形1/8的ΔO'A'B'來做代表。這樣一來問題 域的中子徑跡，可以用柵元ΔO'A'B'中的一條折線軌反過來，在直角三角形柵元ΔO'A'B'中任一條反n，則反射方向Ω為：?2cos cos? n(n,n,n 則uu2cosnxvv2cosnyww2coscos(unxv

設(shè)三角形柵元的橫截面ΔOABX-YX-YX軸成α角的平vusin2vcos2ww由此就可得到OA、OBAB邊上的全反射條件，對α=π/2。 能量E下的反照率β(E)，反射中子的能量分布RE(E→E')。于是代替在反射層中眼蹤中子，我們可WW其中，E為射入反射層中子的能量，W為中子的權(quán)EβRE(E→Eβ中抽樣產(chǎn)生。反射后的方向Ωβ由半平面各向同性分布或余分結(jié)果外，還需要得到這些量的方差、協(xié)方分布。這些量可以利用分類記錄手續(xù)同時得NNi 1NiN又隨所用的量是粒子概率，使用直接模擬法時，如粒子則沒有貢獻(xiàn)。使用法時，如粒子，在疊g2Eg2Egg2EggEg 2

g

g NN

gi

N

NgiNgi

因此，要得到2和22 2g錄結(jié)果的 和gi gg方差估計值g

NN其中為置信限，它隨置信水平1而定。在通常情況下，取10.95, 1.96。 向、時間特征，分別記錄其權(quán)重，最后將這 Emin

E1E0101 K1;0t0t1 tL 在粒子時，對粒子的每一種能量，先 特卡羅，只需要以下群截面：ti(r)tai(r)asij(r)sfs(|r,ij)

在電子計算機(jī)上，方法解粒子輸運(yùn)問題集集 作業(yè)第六 問題。相對來說，點(diǎn)通量的計算要一些。設(shè)(rE,?分別表示粒子的位置、能量和運(yùn)(rE,?)的定義為：(rE,?)dVdEdrdV(r0)dV在r0點(diǎn)的體積元dV內(nèi)，粒子的平均00

(r)dA

(r,E,?)dEd?00A0)ds沿曲面A0的法線方向增加厚度ds00

(r)dV

(r,E,?)dEd?00(V0) 在體V0內(nèi)，粒子的平均徑跡長通量(rE,?可用粒子各次散射對通量的貢獻(xiàn)和(r,E,?)n(r,n(rE,?)dVdEdnn＋1rdV內(nèi)，能量E和運(yùn)動方向Ω屬于dEdΩ的用方法計算通量的能譜與角分布，所采用的與計算其它物理量一樣，即把能量和方向分成若干個區(qū)間，分別按粒子狀態(tài)所處的區(qū)間累積記i(r) (r,E,?)d?i j(r) j E?(r)

(r,E,?)d?

(r,E(r,E,?)dE

當(dāng)E(r,E,?)

E?

E V0*(V)Ws2 s(rl?,E1 1區(qū)域V0的近端和遠(yuǎn)端的交點(diǎn)的距離。如果點(diǎn)rn在V0 *1

* ets1et *(V)

et(V0)(s2s1

(V*(V

10

l?,

sn次散射的通量貢獻(xiàn)為： s *(V)

s1s ss或與V *(V)f(s)ds

*(V) ,E)

s l?,E

*(V)d s l?,E W l?,Es2(ss)d s2(ss)dW l?,E2s1 s1

Wn(s2s1)expW s)

s2 l?,E

W 2 (rl?,E s2 s21sW (rl?,E 1 1sn次散射的通量貢獻(xiàn)為： *(V

ss ,E

*(V)f(s)ds

*(V) ,E)

s l?,E

,E)

s l?,E1 t(rn1,En1

W 2exp(rl?,E 1 1fn(s)s*n次散射通量貢獻(xiàn)的估計為： *(V

exp0t(rnl?n,EnW

fn(s*fn(s*)

s2

s1s

0fV(r)為在V0上定義的任一概率密度函數(shù)，則0

V(r)dV

(r) V V0

(r)dV*(V)

(r*V0V *(A)

|n?

s1(rl?,E

n有多個交點(diǎn)，則*An

交點(diǎn)，則*A0 sn次散射的通量貢獻(xiàn)為：

s*(A)

|n?

A0設(shè) (r)為在A上定義的任一概率密度函數(shù)，A00(A0)A(r)dA

(r) A A

(r)dA0

0*(A)

(r*A0AAA0

0)0) VA*(A)lim1 s0 110*(A)0

s(VA00計算最。這是因?yàn)?，在大量的模擬粒子n次散射后粒子的狀態(tài)為(rnEn,?n,Wn，n次碰撞的粒子的狀態(tài)為(rnEn1,?n1,Wn1，C(En1En,?n1C(EE,??r)

*(r*)

C(E

E ?*r)

|r*rn

|r*rn|

l?*,E)dl

r*n |r*rn當(dāng)n＝0時，用源分布密度函數(shù)S(r0,E0,?0) 核C(En1En,?n1?nrn)。C(EE,??r)K(EEr)111t(r,

其中光子能量E以電子靜止能量mc2≈0.511MeV為單位；K(E'→E/r)為Klein－Nishina K(EEr)πN(r)z(r)r

1

1

E2

當(dāng) EE2EKKn1*nrnnn

Wn1

(r, 2π|r*r exp

|r*rn|

l?*,E*)dl

nE*n

1 C(EE,??r)A,i(E)NA(r)A,i(E)(r,

(EE,? A,i(EA,i(E)分別表示能量為E'的中子與第A種i種反應(yīng)后產(chǎn)生的平均次級中子數(shù)和微觀截面；NA(r)r處第A種原子核的核密度；fA,i(EE,??表示能量為E'和方向?yàn)棣?的中子*(r*)

A,i(En1)NA(rn)A,i(En1

t(rn,En1

exp

|r*rn|

n n

|r*

r 概率方法的估計量中含有因子 |r*rn 狀態(tài)為(rE,?P(r,E,?)(r,E,?)|nPAP(r,E,?)dEd?A(r,E,?)|n?|dEd?

(r,E,?)t(r,E)(r,f(r,E,?)f(r,E)f(r,E)(r,

a(r,E,?)a(r,E)(r,r(r,E,?)r(r,E)(r,第七 第七 方法求積分的一般規(guī)則如下：任何一個積分，都可看作某個