第五章蒙特卡羅方法在計算機上的實現(xiàn)

1、1.源分布抽樣過程源分布抽樣過程2.空間、能量和運動方向的隨機游動過程空間、能量和運動方向的隨機游動過程3.記錄貢獻和分析結(jié)果過程記錄貢獻和分析結(jié)果過程4.核截面數(shù)據(jù)的引用核截面數(shù)據(jù)的引用5.蒙特卡羅程序結(jié)構(gòu)蒙特卡羅程序結(jié)構(gòu)作作 業(yè)業(yè) 蒙特卡羅方法是隨著計算機的出現(xiàn)和發(fā)展而逐步發(fā)展起來的。在計算機上能夠產(chǎn)生符合要求的隨機數(shù)，實現(xiàn)對已知分布的抽樣，奠定了蒙特卡羅方法在計算機上得以實現(xiàn)的基礎(chǔ)。在計算機上使用蒙特卡羅方法解粒子輸運問題大致包括三個過程：源分布抽樣過程，空間、能量和運動方向的隨機游動過程以及記錄、分析結(jié)果過程 。 源分布抽樣的目的是產(chǎn)生粒子的初始狀態(tài) 。下面我們介紹一些常見的特定類型的

2、源分布抽樣方法。),(0000rSE(1) 圓內(nèi)均勻分布設(shè)圓半徑為R0，粒子在圓內(nèi)均勻分布時，從發(fā)射點到中心的距離 r 的分布密度函數(shù)為： r 的抽樣方法為：其它當(dāng)002)(020RrRrrf),max(210 Rr(2) 圓環(huán)內(nèi)均勻分布設(shè)圓環(huán)的內(nèi)半徑為R0，外半徑為R1，則粒子在該圓環(huán)內(nèi)均勻分布時，從發(fā)射點到中心的距離 r 的分布密度函數(shù)為： r 的抽樣方法為：其它當(dāng)02)(102021RrRRRrrf02010320101011)(),max()(RRRrRRRrRRRR(3) 球內(nèi)均勻分布設(shè)球的半徑為R，粒子在球內(nèi)均勻分布時，從發(fā)射點到中心的距離 r 的分布密度函數(shù)為： r 的抽樣方法為

3、：在直角坐標系下，抽樣方法為：其它當(dāng)003)(32RrRrrf),max(321 Rr302010232221,1RzRyRx(4) 球殼內(nèi)均勻分布設(shè)球殼的內(nèi)半徑為R0，外半徑為R1，在均勻分布時，從發(fā)射點到中心的距離 r 的分布密度函數(shù)為： r 的抽樣方法為：其它當(dāng)03)(1030312RrRRRrrf001432322211020101211020201)(),max(),max(33RxRRrxxxRRRRRRRRRRR在直角坐標系下，球殼內(nèi)點的坐標為： 其中，r 由前面的抽樣方法確定，、服從各向同性分布，其抽樣方法為：cossinsincossin000rzryrx2322222123

4、222221023222221310232222212101223222221cos2sinsin2cossin)(AAAArrzAAArryAAArrxAA(5) 圓柱內(nèi)均勻分布圓柱內(nèi)均勻分布是指粒子發(fā)射點均勻地分布在底半徑為 R，高為 2H 的圓柱內(nèi)。若固定圓柱的中心為原點，圓柱的軸向為 z 軸，則分布密度函數(shù)為： 抽樣方法為：其它當(dāng)0| ,21),(2222HzRyxHRzyxf3020102221,1HzRyRx(6) 點源分布 點源分布是指粒子由一固定點 發(fā)射，其分布密度函數(shù)為：其中， 為狄拉克函數(shù)，源粒子的抽樣方法為：在球坐標系中，粒子發(fā)射點到球心的距離 r 的分布密度函數(shù)為：其中

5、， 為點源到球心的距離。源粒子的位置抽樣為：)()()(),(*0*0*0zzyyxxzyxf),(*0*0*0zyx*0*0*0,zzyyxx*0rr *0r)()(*0rrrf)(7) 球外平行束源分布球外平行束源分布是指粒子平行入射到半徑為 R 的球面上，或球外點源距離球很遠，可以近似地看作平行束源。設(shè) r 為粒子發(fā)射點到球心的距離 , 其分布密度函數(shù)為： r 的抽樣方法為：在直角坐標系中，抽樣方法為：)()(RrrfRr 20202020102221,1zyRxRzRy(1) 單能源分布單能源分布是指粒子的發(fā)射能量為一固定值 E0 ，其分布密度函數(shù)為 ：源粒子的能量為：)()(0EEE

6、f0EE (2) 裂變中子譜分布裂變中子譜分布的一般形式為： 其中A，B，C，Emin，Emax 均為與元素有關(guān)的量。對于鈾-235，A=0.965，B=2.29，C=0.453，Emin=0，Emax=。maxmin,sh)(EEEBEeCEfAE采用近似修正抽樣，抽樣方法為：其中，m0.8746，M10.2678，0.5543。此外，裂變譜分布也有以數(shù)值曲線形式給出的，此時，用數(shù)值曲線抽樣方法抽取 E 。3321121ln1)ln()(AAEEAEMEHmEAmEBECAAEHexpsh1)(21(3) 麥克斯韋(Maxwell) 譜分布麥克斯韋譜分布的一般形式為：該分布的抽樣方法為0,2

7、)(23EeEEfE22221ln23lnEe(1) 各向同性分布各向同性分布密度函數(shù)為： 其中，cos，為運動方向與 z 軸的夾角，為方位角。21)(,21)(41)()()(2121fffff在直角坐標系下，各方向余弦 u，v，w 為：其抽樣方法為：cossinsincossinwvu2322222123222221232222213123222221211223222221cos2sinsin2cossin)(AAAAwAAAvAAAuAA(2) 半面各向同性分布不妨設(shè)在 x0 的半面方向上各向同性發(fā)射粒子，則在前述各向同性分布的抽樣方法中，用2代替2就能得到所需分布的抽樣。對于其它方向

8、的情況，可用類似的方法處理。(3) 球外平行束源分布令cos，為粒子運動方向的徑向夾角，則分布密度函數(shù)為：的抽樣方法為：01,2)(f),max(21u(4) 球外各向同性點源分布設(shè)球外點源 S 到球心的距離為D0。點源 S 到球的最大張角為*，則球外各向同性點源分布的抽樣方法是：先抽樣確定 ，再轉(zhuǎn)換成。0220*cosDRD RDR2202*sincos)cos1 (1cos在直角坐標系下，取 OS 為 z 軸，抽樣方法為：cos0sinwvu 在有關(guān)次級粒子（如裂變中子，中子生成光子，光子生成中子）的輸運過程中，次級粒子源分布的抽樣方法，主要可分為以下兩種：(1) 直接生成法 可將生成的次

9、級粒子的位置、能量、方向、權(quán)重等參數(shù)直接作為源分布的抽樣結(jié)果。也就是直接對生成的次級粒子進行跟蹤。這種方法比較簡單、直觀。(2) 離散分布法 將生成的次級粒子的權(quán)重，按空間位置、能量、方向分別記錄，得到次級粒子的空間、能量、運動方向的離散的近似分布。再根據(jù)該分布，利用各種抽樣技巧，得到源分布的抽樣，對抽樣的源粒子進行跟蹤、記錄。 當(dāng)一個問題需要用兩個以上的蒙卡程序處理時，可采用這種方法。 粒子由狀態(tài)Sm到狀態(tài)Sm+1時，需要確定粒子的空間位置 rm+1，能量 Em+1和運動方向m+1。 設(shè) rm 為粒子第 m 次碰撞點的位置，m 為碰撞后的運動方向，則粒子第 m+1 次碰撞點的位置 rm+1

10、為：即其中 為 的方向余弦，L 為兩次碰撞點間的距離。mmmmmmmmmwLzzvLyyuLxx111mmmL rr1mmmmwvu),(L 的分布密度函數(shù)為：由 f (L) 抽樣確定 L 的方法通常有三種：(1) 直接抽樣方法確定 L 的直接抽樣方法是：首先由自由程分布中抽取再由下列關(guān)系式解出 L 。0,),(exp),()(01LdlElELfLmmmtmmtrrLmmmtdlEl0),(r ef)(ln對于均勻介質(zhì)，有對于多層介質(zhì)，如果則其中， 為粒子由 rm 出發(fā)，沿m 方向在順序經(jīng)過的第 i 個介質(zhì)區(qū)域內(nèi)走過的距離， 為第 i 個介質(zhì)區(qū)域的宏觀總截面 ( i =1，2，Imax )。

11、 當(dāng)時，意味著粒子穿出系統(tǒng)。)(ln)(mtmtEEL)()()()(,10,100,10,mItIimitiIiiIimitiIimitiEELLLELELiL)(,mitEmax0,)(IimitiEL(2) 最大截面法對于多層介質(zhì)，或其他介質(zhì)密度與位置有關(guān)的問題，在求 ( i =1，2，Imax ) 時，如果系統(tǒng)形狀復(fù)雜，計算是非常煩雜的。在這種情況下，使用最大截面法更方便。最大截面抽樣方法為：其中iL),(max)(max,EEttrr1max,12max,1111)(),()(ln0LLEELELLLmtmmmtmtr(3) 限制抽樣法 當(dāng)介質(zhì)區(qū)域很小時，如使用直接抽樣法抽取輸運長度

12、，粒子很容易穿出介質(zhì)，此時使用限制抽樣法確定自由程個數(shù)較好，的分布密度函數(shù)為：其中 Dm 為粒子由rm 出發(fā)，沿m 方向到達區(qū)域邊界的自由程個數(shù)。的抽樣方法是：然后用直接抽樣法中根據(jù)計算 L 的方法計算輸運長度 L 。此時，粒子的權(quán)重需乘以糾偏因子 。其它當(dāng)001)(mDDeefm)1 (1lnmDe)1 (mDe 粒子在介質(zhì)中發(fā)生碰撞后，首先要確定與哪種原子核發(fā)生何種反應(yīng)。粒子發(fā)生碰撞后（吸收除外）的能量 Em+1 一般只與其碰撞前后運動方向的夾角（散射角）有關(guān)。 粒子碰撞后常見的能量分布有下面幾種情況。(1) 裂變中子譜 中子引起原子核裂變反應(yīng)時，裂變中子的能量服從裂變譜分布。其抽樣方法可

13、參考以前的介紹。(2) 中子彈性散射后能量的確定 中子彈性散射后，能量與質(zhì)心系散射角C的關(guān)系是：能量與實驗室系散射角L的關(guān)系是：其中，A 為碰撞核的質(zhì)量， ?；?確定后，即可求出 Em+1。222211) 1(LLmmAAEE12) 1(221CmmAAAEELLCCcos,cosLC(3) 中子非彈性散射后能量的確定 中子非彈性散射后，能量與質(zhì)心系散射角C的關(guān)系是：其中， 為第 K 個能級的閾能， 為第 K 個能級的激發(fā)態(tài)能量。 如果確定了實驗室系散射角L，則根據(jù)下式確定 后，再計算出 Em+1。KKCmKmKmmAAEAEAAEE11121) 1(221KKC11111222LLmKLmK

14、CEAEA(4) 光子康普頓（Compton）散射后能量的確定光子發(fā)生康普頓散射后，其能量分布密度函數(shù)為：其中， K() 為歸一因子。 ， 和 分別為光子散射前后的能量，以 m0c2 為單位，m0為電子靜止質(zhì)量，c 為光速。211,1111)(1)(322xxxxxxKxf22)21 (21421)21ln() 1(21)(Kx光子康普頓散射能量分布的抽樣方法為：x 的抽樣確定后，散射后的能量為：232223221121321) 1(427212942711212121xxxxxxxxxxEcmxcmEmm 20201粒子碰撞后運動方向m+1的確定，一般與散射角有關(guān)。由已知分布抽樣確定散射角后

15、，再確定m+1。常見的散射角分布有如下幾種：(1) 質(zhì)心系各向同性分布散射角在質(zhì)心系服從各向同性分布時，其抽樣方法為 。質(zhì)心系散射角C抽樣確定后，需轉(zhuǎn)換成實驗室系散射角L。12 C在中子彈性散射情況下，轉(zhuǎn)換公式為：其中 A 為碰撞核質(zhì)量， 。在中子非彈性散射情況下，轉(zhuǎn)換公式為：其中， 為第 K 個能級的閾能。CCLAAA2112CmKmKCmKLEAEAEA1211112LLCCcos,cosK(2) 中子彈性散射勒讓德 (Legendre) 多項式分布 中子彈性散射角分布常以勒讓德多項式的展開形式給定。散射角余弦 xcos的分布密 度函數(shù)為：其中 Pl(x) 為 l 階勒讓德多項式。該分布即

16、為 n 階勒讓德近似展開。勒讓德多項式由以下遞推公式確定：其它當(dāng)01|)(212)(0 xxPflxfnlllxxPxPxnPxxPnxPnnnn)(1)(0)()() 12()() 1(1011考慮新的分布：當(dāng)選取 x0，x1， xn 為 Pn+1(x)0 的根，且時，fa(x) 依照勒讓德多項式展開的前 n 項與 f (x) 的展開形式相同。因此，可以用 fa(x) 作為 f (x) 的近似分布。nlklnlkllkxPlxPfl020)(212)(212nkkkaxxxf0)()(在實際問題中，由于勒讓德多項式展開項數(shù)不夠，可能出現(xiàn)某個 為負值的現(xiàn)象。此時可以采用如下近似分布：其中：對于

17、該近似分布，可用加抽樣方法進行抽樣：此時，由于偏倚抽樣而引起的糾偏因子為 wK ，也就是說，粒子的權(quán)重要乘上wK。nkkkkknkkkknkkka|w|xxxf000*|,|)()(kKkkKkkKxx010,當(dāng)(3) 光子康普頓散射角分布光子的康普頓散射角與其散射前后的能量有關(guān) , 它的分布密度函數(shù)為：抽樣方法為：)111 ()(LLf111L實驗室系散射角L確定后，依據(jù)不同的坐標系的表現(xiàn)形式，有不同的確定方法。(1) 確定方向余弦 um+1，vm+1，wm+1mmmmmmmmmmmmmmmmmmawvubcwavvubduvbcwvauvubdvubcwu221221221其中，方位角 在

18、 0, 2 上均勻分布。 當(dāng) 時，不能使用上述公式，可用下面的簡單公式：sin,cos,1sin,cos2dcabaLL022mmvummmmawwbdvbcu111(2) 確定m+1的球坐標 (m+1，m+1)設(shè)m的球坐標分別為 (m,m)，其中，為粒子運動方向與 z 軸的夾角， 為粒子運動方向在 x y 平面上投影的方位角。則m+1的球坐標 (m+1,m+1) 分別由下式確定：111111sinsincoscoscos)cos(sinsinsin)sin(cossinsincoscoscosmmmmLmmmLmmLmLmm 一般幾何的隨機游動公式可以應(yīng)用到球形幾何，而對球?qū)ΨQ問題，使用特殊

19、形式更為方便。(1) 下次碰撞點的徑向位置 rm+1的確定 兩次碰撞點間的距離 L 確定之后，下次碰撞點的徑向位置 rm+1的計算公式為：設(shè)系統(tǒng)的外半徑為R，如 rm+1R，則粒子逃出系統(tǒng)。mmmmLLcos2221rrr(2) 粒子碰撞后瞬時運動方向的確定 在球?qū)ΨQ系統(tǒng)中，粒子運動方向用其與徑向夾角余弦來描述。使用球面三角公式，粒子碰撞后瞬時運動方向與徑向夾角余弦 cosm+1的計算公式為：其中， 為在 0, 2 上均勻分布的方位角， 為在 rm+1 點進入碰撞前瞬時運動方向與 rm+1 徑向之間的夾角。111coscossinsincossinsincoscoscosmmmmmmmmLmL

20、mmLrrrrm在抽樣確定輸運距離、判斷粒子是否穿透系統(tǒng)時， 常遇到求由 rm 出發(fā)，沿m 方向到達某個區(qū)域表面的距離問題。在記錄對結(jié)果的貢獻時，也常使用類似的量。區(qū)域表面通常是平面或二次曲面。 求到達區(qū)域表面的距離問題，實際上是求直線（或半射線）與平面或二次曲面的交點問題。這是 蒙特卡羅方法解粒子輸運的各種實際問題時 , 所遇到的基本幾何問題。(1) 點到平面的距離點 沿方向 的直線方程為：該直線到達方程為的平面的距離為：當(dāng) 與平面平行時，即直線與平面無交點。如果 l 為負值，直線與平面也無交點。這時，粒子的運動方向是背離平面的。),(0000wvu0000cwbvau000000)(cwb

21、vauczbyaxdldczbyax00rrl),(0000zyxr(2) 點到球面的距離在三維直角坐標系中，設(shè)球心為 rc(xc, yc, zc) ，球半徑為 R，則球面方程為：將直線方程代入該球面方程，得到點 r0沿0方向到達球面的距離 l ：其中2222)()()(Rzzyyxxcccl2020202020220200000000)()()(|)()()()(cccccccczzyyxxRRRwzzvyyuxxrrrr當(dāng) R0R 時，即 r0 點在球內(nèi) ，2，l 只有一個正根：當(dāng) R0R 時，即 r0 點在球外，分以下三種情況：a) 若0，l 無正實根，直線與球面無交點。b) 若0，0，

22、l 無實根，直線與球面無交點。c) 若0，0，l 有兩個正實根，直線與球面有兩個交點。ll在球坐標系中，不失一般性，設(shè)球心為 rc0，則球面方程為 rR。當(dāng) r0R 時，即 r0 點在球內(nèi) ，有一個交點：其中0為0與 r0 的徑向夾角。當(dāng) r0R 時，即 r0 點在球外 ，令當(dāng) cos00 時，直線與球面無交點。當(dāng) cos00 時，若 dR，則直線與球面無交點。若 dR，則有兩個交點：200200)sin(cosrRrl200200)sin(cosrRrl00sin rd(3) 點到圓柱面的距離設(shè)圓柱面的方程為：其中 (xc, yc, 0) 為圓柱的中心，R 為圓柱底半徑。點 r0沿0方向到達

23、圓柱面的距離 l 為：其中222)()(Ryyxxcc201wl2020202022020000)()()1)()()(ccccyyxxRwRRvyyuxx 當(dāng) R0R 時，r0 點在圓柱內(nèi) ，如果 ，則 l 有一個正根： 如果 ，即0平行于圓柱的對稱軸，直線與圓柱面無交點。 當(dāng) R0R 時，r0 點在圓柱外，分以下三種情況：a) 若0，l 無正實根，直線與圓柱面無交點。b) 若0，0，l 無實根，直線與圓柱面無交點。c) 若0，0 且 ，l 有兩個正實根，直線與圓柱面有兩個交點。d)在 的情況下，直線與圓柱面不相交。201wl201wl120w120w120w120w(4) 點到圓錐面的距離

24、 設(shè)圓錐頂點在原點，以 z 軸為對稱軸，則圓錐面的方程為：點 r0沿0方向到達圓錐面的距離 l 為：其中如果0與錐面某一母線平行，即 ，則2222zcyx202)1 (1wcl)1 (1)(202202202020020000wczcyxwzcvyux22022020zcyxl1)1 (202wc(5) 空腔處理 在粒子輸運問題中，所考慮的系統(tǒng)常有空腔存在，如中空的球殼 , 平板間的空隙等。粒子輸運時，可有兩種處理空腔的方法：a) 將空腔作為宏觀總截面t0 的區(qū)域 , 按通常的方法輸運。b) 設(shè) 分別為由 rm 出發(fā)，沿m 方向到空腔區(qū)域的近端和遠端的交點，當(dāng)粒子超過 時，以 為新的起點，重新

25、開始輸運。 顯然，這兩種方法在統(tǒng)計上是等價的。mmrr 、mrmr (1) 全反射邊界在反應(yīng)堆活性區(qū)中，元件盒常常按正方形或六角形排列。假定元件盒足夠多，每個盒結(jié)構(gòu)相同，那么活性區(qū)中每個盒所占的柵胞的物理情況，可以代表整個活性區(qū)中的狀況。 進一步假定，元件盒是圓對稱的，那么每個柵胞中情況，可以用更小的單位（柵元）來反映。比如對六角形柵胞可取其 1/12 的OAB 來做代表；正方形柵胞可用其 1/8 的OAB 來做代表。這樣一來問題就大大簡化了。 現(xiàn)在的問題是怎樣計算直角三角形柵元的物理量（如通量）。用蒙特卡羅方法如何模擬中子在柵元內(nèi)的運動，反映出整個活性區(qū)對它的影響。我們可把OA、OB、AB

26、作為全反射邊界來處理。所謂全反射邊界，就是當(dāng)中子打到該邊界上時，按鏡面反射的方式，從邊界 上全部反射回來，中子的能量與權(quán)重均不改變。 在這種邊界上的反射條件，稱之為全反射條件，就是通常的鏡面反射條件。 在全反射邊界條件下，一條通過活性區(qū)若干個區(qū)域的中子徑跡，可以用柵元OAB 中的一條折線軌 道來反映出來。 反過來，在直角三角形柵元OAB 中任一條反射成的折線軌道，都代表了中子在活性區(qū)內(nèi)一條直線軌道的作用。由于系統(tǒng)的對稱性，在活性區(qū)內(nèi)，凡是與柵元內(nèi)位置相當(dāng)?shù)牡胤?，都有相同的物理情況，因此柵元內(nèi)各處的情況，當(dāng)然代表了整個活性區(qū)的情況。(2) 一般曲面全反射條件 對于一般曲面的全反射，設(shè)入射方向為，

27、入射點的內(nèi)法線方向為 n ，則反射方向 為：其中設(shè)則nnncos22ncos),(zyxnnnn)(coscos2cos2cos2zyxzyxnwnvnunwwnvvnuu(3) 平面全反射條件 設(shè)三角形柵元的橫截面OAB 在 X-Y 平面上，OAB。則邊界 OA、OB、AB 上的反射都是平面全反射。在任一與 X-Y 平面垂直且與 X 軸成角的平面上，全反射條件為：由此就可得到OA、OB 和 AB 邊上的全反射條件，對于 OB 邊，=；對于 OA 邊，= 0；對于 AB 邊，=/ 2。wwvuvvuu2cos2sin2sin2cos(4) 反射層邊界條件 對于具有大反射層的系統(tǒng)，如存放，運輸和

28、生產(chǎn)裂變物質(zhì)的倉 庫、車廂和車間等，當(dāng)中子從里面打到四周墻上或反射層時，還要繼續(xù)對它進行跟蹤。這種跟蹤常常要花費很大的計算量，并且在結(jié)果中引起的方差也比較大。如果在計算這種系統(tǒng)的不同方案中，反 射層條件不變，那么這種大量重復(fù)的計算是很不經(jīng)濟的。 中子射入反射層后，一部分被介質(zhì)吸收，只有一部分返回，由于中子的散射慢化，損失一部分能量，因此反射回來的中子有一個能量方向分布。顯然，對這種反射層，不能應(yīng)用全反射條件。不過，我們?nèi)匀豢梢园阉?dāng)做邊界，在邊界上按反射層的物理作用來處理。 比如，如果反射層是一種平板幾何，我們可以用數(shù)值方法或蒙特卡羅方法，預(yù)先算好在各種不同入射能量 E 下的反照率(E)，反射

29、中子的能量分布 RE(EE )。于是代替在反射層中眼蹤中子，我們可在反射層邊界上作如下處理：一旦中子打入反射 層，立即返回，反射后權(quán)重為其中，E 為射入反射層中子的能量，W 為中子的權(quán)重。反射后的能量 E 由反射能譜 RE(EE) 中抽樣產(chǎn)生。反射后的方向 由半平面各向同性分布或余弦分布中抽樣。反射后的中子位置為入射時的位置。 計算表明，對于大尺寸的反射層來說，這樣的近似，引 起的結(jié)果上的誤差是可以忽略的，卻能帶來計算量的大量節(jié)省。)(EWW 在粒子輸運問題中，除了得到某些量的積分結(jié)果外，還需要得到這些量的方差、協(xié)方差、以及這些量的空間、能量、方向和時間的分布。這些量可以利用分類記錄手續(xù)同時得

30、到。 為了得到所求量的估計值，在粒子輸運過程中需進行記錄，即求每個粒子對所求量的貢 獻。 設(shè)模擬了 N 個粒子，所求量的估計值為：其中 gi 為第 i 個粒子的總貢獻。NiiNgNg11 記錄的貢獻由所求量決定。對于同一個所求量，又隨所用的蒙特卡羅技巧的不同而不同。 例如，所求量是粒子穿透屏蔽概率，使用直接模擬法時，如粒子穿透屏蔽，在疊加記錄單元加“1” ( 初始值為零 )，否則沒有貢獻。使用加權(quán)法時，如粒子穿透屏蔽，在疊加記錄單元加粒子的權(quán)重，否則沒有貢獻。使用統(tǒng)計估計法時，粒子每發(fā)生一次碰撞 (包括零次碰撞)，都要記錄貢獻，等等。估計量 g 和 g 的方差和協(xié)方差為：它們可以用下式估計：N

31、iiNiiNiiiggNiiNiiggNgNggNgNgN11122112211111 gEEgggEEggEggg2222因此，要得到 和 的估計，只要對每一個歷史記錄結(jié)果的 和 進行記錄，并加以累加即可。方差估計值 確定后，可得到誤差其中 為置信限，它隨置信水平 而定。在通常情況下，取 。22gggiiiggg22gNg96. 1,95. 011 在前面已經(jīng)提到，用蒙特卡羅方法求某種量的空間、能量、方向和時間分布，實質(zhì)上是得到這種分布的階梯函數(shù)近似的估計值。而求這種估計值是很方便的，只要將跟蹤過程中所得到的感興趣量，按其狀態(tài)的空間、能量、方向、時間特征，分別記錄其權(quán) 重，最后將這些記錄結(jié)果適當(dāng)處理即可。 事先，將問題的空間、能量、方向（常按相對于某個方向的夾角余弦）、時