[分享]关于蒙特卡罗方法的一些简单介绍，我们广大的新手...

DoubleX

最近买了一本书，书上有关于蒙特卡罗（mc）方法的简单介绍，下面我就分几天时间把这书上的介绍给发上来吧。让我们广大的新手能从最基本的学起什么是蒙特卡罗方法。

1、概述

mc方法又叫做随机取样、统计模拟或统计实验方法。mc方法是一种利用随机数的统计规律来进行计算和模拟的方法。它可用于数值计算，也可用于数字仿真。在数值计算中常用于多重积分、线性方程组求解、矩阵求逆以及方程求解，包括常微分方程，偏微分方程，本征方程、非齐次线性积分方程和非线性方程等。在数字仿真方面，常用于核系统临界条件模拟、反应堆模拟以及高能物理，实验核物理、统计物理、真空、地震、生物物理和信息物理等领域。mc方法的缺点是收敛速度慢，计程长。为对mc方法有个初步的认识先看个用mc方法求圆周率的例子。

设园的半径是r，圆心位于xoy平面的（r,r）处，且内切于边长为2r的正方形中，那么显然正方形的面积是S1＝4r×r，用mc方法计算圆的面积的基本思想是随机的正方形范围内部画点，若共画了N个点，而落在圆内部的有M个点，当N足够大的时候，圆的面积可为s=M/N*S1,基本方法是首先产生 两个随机数x和y其值域为[0，2r]，然后进行判断（x,y）是否落在圆内

（x－r）2＋（y-r）2<r2

记录下总的点数N和落在圆中的点数M，则圆的面积为S＝4r2M/N，从而得pi＝4×M/N

这里要注意两点：第一，两个随机数要是在矩形范围内均匀分布得。大部分计算机语言提供了随机数发生器，可产生0－1之间得均匀随机数，经过适当得变换不难变成任意值域得均匀随机数。第二，随机数得个头N必须要足够大。已确保一定得精度。负责由于统计涨落过大而使得误差很大，这是统计方法得自身得规律。

mc方法不仅能用于物理方程得数值计算，也可以用于物理过程得数字模拟。氢原子的电子云模拟就是一个简单的例证，由量子力学可以得知氢原子的s态波函数ψ＝ψ（r)只是半径的函数，于θ，φ无关，而氢原子中的电子沿半径r的分布密度即电子在半径r处单位厚度球壳内出现的 几率

D=4×pi×r2×ψ2

习惯上把这种分布称作为电子云。

氢的基态1s态（n=1,l=0,m=0），有

D=1/(a1*a1*a1)4×r2×e-2r/a1

Dmax=1.1

r0=0.25nm

其中a1＝5。29×10－2 nm，是D的最大值Dmax处的r值，其值是与波尔半径相同，r0是Dmax处的r值，其值与波尔半径相同。r0是D收敛处的r值，即D的收敛点。

氢的2s态（n=2,l=0,m=0）有

D=1/(8×a1*a1*a1)×（2－r/a1)2×e-r/a1

Dmax=0.14

r0=1.0nm

氢原子的3s态（n=3,l=0,m=0）有

D=1/(81×81×3×a1*a1*a1)[27-18*r/a1+2*(r/a)2]×e-2r/3a1

Dmax=0.2

r0=2.0nm

氢的电子云的模拟是依据上述的分布函数用绘图点的密度来描绘电子的概率分布函数，在氢原子例子中设随机数发生器可产生0－1的随机数，记做rand(k)（k=1,2,3……），首先利用一个随机数rand(1)来产生一个随机的电子半径轨道

r=r0*rand(1)

显然0<r<r0，由r计算出D（r)。再产生一个随机的概率判据D0＝Dmax* rand(2)进行判断，如果D(r）<D0，则从新开始，否则再产生一个 随机的角度值

θ＝2*pi*rand(3)。最后计算要描绘的点的坐标x=srcosθ y=srsinθ。其中s是控制图形的尺寸因子，与每纳米点数相对应。绘出点（x,y)后再从新上述的过程，经过足够多的次数电子云的图貌就由点的疏密反应出来了。

DoubleX

2、假随机数的产生
真的随机数如同掷骰子一样，产生1－6范围内的随机整数，抽奖用的摇号码机产生的0－9范围内的随机整数，这些真的随机数除了统计规律外无其他规律可循。假随机数就是用某种算法给出的似乎随机的数，既然数的给出是按照某种算法的那么就有某种规律，那么这些数必然有一定的周期，设周期是n那么n+1就一定是第一个数了，此后均要一次重复出现。当然如果周期足够大，可使在整个使用过程中不至于呈现周期性，那么假随机数也是使用的。例如，让计算机的假随机数发生器其周期大于机器的记忆单元数。此外假随机数的统计性质也是表征随机数品质的又一个重要的指标。均匀分布的随机数，即要求数的出现是随机均匀的，也要求数的分布是均匀的。至于评价数的均匀，这里不作介绍。
2、1 均匀随机数的产生
产生均匀随机数常用幂剩余法，例如可按如下的公式产生随机数
xn=cx n-1(modN)
或者写成xn=cx n-1－Nint（cx n-1/N)
其中c、N为给定的常数，给出x0后，就可以用上式依次给出x1，x2……等等，如何确定常数c、N、X0是十分关键的问题，人们仍然在不断的研究中，现面给出c、N和X0的一般原则。N的取值一般取N＝2 m-1，其中m为计算机的二进制数的字长，N-1一般取计算机能表现出来的最大整数，如字长16位的N=2^15=32768;字长32位的N=2^31=2147483648,c的取值一般取c=8*M+-3，其中M位任一正整数，如有取c=16897，c=65539，c=397204099等，建议取c~N 1/2，这样统计性比较好，关于x0的取值，一般取x0为奇数，如x0=13。可以检验其他条件不变，x0是奇数时，周期是T,偶数时周期时T/2。如N=64,c＝5，x0＝2则得
x1=10
x2=50
x3=58
x4=34
x5=42
x6=18
x7=26
x8=2
如此可见x8=x0，可见这列数字得周期是8，若令x0=1，其余不变，则数列为
5，25，51，49，53，9，45，33，37，57，29，17，21，41，13，1
则周期为16
若按照上式产生得随机序列其值域为0～N-1，若要产生0～1之间的得随机数只要把每个随机数在除以N即可。
2、2随机性统计检验
一个好得随机数发生器或一个好得随机数生成程序要有两个条件
1、生成的随机数有足够长的周期。2、生成的随机数列应当有真正随机序列的统计性质。其周期的长短比较容易测试和判断，尔统计性质的优劣则不那么简单了。下面着重讨论统计性质的两种常用的检验方法，频数分布检验和行程频数检验。
先讨论频数分布检验。对于一个均匀分布的随机数发生器，设所产生的随机数序列的值域为[0,1]，则所产生的随机数字应与从0～1之间均匀的频数分布一致，为检验频数分布情况，可按画统计直方图的方法，将整个至于分成M个宽度相等的子区间，设xi是第i子区间内出现的随机数的个 数，即第i个子区间的频数，则所有子区间中随机数的个数的平均频数
郁闷，没有公式编辑器，这后面的公式就不介绍了，大概了解下算了
就是每一个子区间中出现x个数字的几率要服从高斯分布规律
行程频数检验中的行程是指数字序列中单调上升或单调下降的连续数字串中数字的个数。由统计规律可知，对一个真正随机数序列其中行程长度k出现次数的期望值为E(K)
k=1,E(K)＝1/12*(5N+1)
k=2,E(K)=1/60*(11N-14)
……
k<N-1,E(K)=2/(k+3)!*[N(k2+3k+1)-(k3+3k2-k-4)]
行程长度k为任意值E(K)＝1/3(2N-1)
检验的方法就是输出给定随机数序列中出现的各种行程长度的次数，并与上述的期望值想比较，符合程度好的随机性好。
上述的两种检验方法都获得满意结果的随机数，一般 来说是可以信任的均匀分布的随机数序列。但仍然不是真正的随机数，因为并不能保证他们满足真正随机数的其他统计性质。


[em01][em01]手敲的酸死了。今天先到这里吧，剩下的以后在发

DoubleX

3、用mc方法计算积分
已知函数f(x)在区间[a,b]上连续，A点x=a，B点x=b，其函数曲线f(x)－x，如图。
下面介绍用mc方法来计算对f(x)的在[a,b]上的定积分
如图在图上做一个矩形，其宽度是b-a，高度是f(x)在[a,b]上的最大值f(m)，那么矩形的面积就是s1=(b-a)*f(m)，给出两个随机数xi和yi，且满足
a<=xi<=b,0<=yi<=f(m)，则所有的随机点都落在图中的矩形内。如果不等式yi<f(xi)成立，那么随即点（xi,yi）还落在矩形中的阴影区域内。
总共设产生了N个随机点，落在阴影中的点数为M，那么当N足够大的时候对f(x)求区间[a,b]上的积分就换算成了求阴影的面积s2=M/N*(b-a)*f(m)。其实这和前面计算圆的面积没有本质的差别。对于比较复杂的问题，其基本处理方法还是相同的。如用mc方法计算一个具有空洞的球体质量和质心坐标等。

DoubleX

下面利用mc方法来计算一个具有空洞球体得质量和质心坐标
设已知物体得轮廓满足曲面方程z=f(x,y)
定义域为x[-a/2,a/2]；y[-b/2,b/2]值域为z[-c/2;c/2]。物体密度是ρ＝ρ（x,y,z)，体积元dV＝dxdydz，则质量元为dm=ρ（x,y,z)dxdydz
那么物体质量为m=∫∫∫ρ（x,y,z)dxdydz
对于均匀分布得物体，即密度为常量，那么物体的密度为m=ρ∫∫∫dxdydz，用mc方法计算物体的体积可以做一个中心位于原点的长方体，将上述物体包含其中。长方体的边长分别是a、b、c。则长方体的体积Vs＝abc，假设整个长方体均充满与物体相同的物质，则其质量ms＝ρabc
然后产生3个随机数xi,yi,zi其值域分别是[-a/2,a/2]、[-b/2,b/2]和[-c/2;c/2]，其中i的取值1～N。这一组随机数代表长方体空间内的一个点。如果一******生了N组随机数，则对应于长方体中的N个点。随机数是均匀分布的，则对应的点再长方体内的分布也是均匀的。显然可以理解为将长方体均等分N个小长方体，每个小长方体的体积是ΔV＝1/N×Vs＝1/N×abc，小长方体质量或称为质量元Δm＝1/N*ρabc
也就是说，一组坐标（xi,yi,zi）就代表一个小体积，也代表一个质量元，统计所有位于物体内的点数M，当N足够大的时候，可以得到物体的体积
V＝MΔV＝M/N*VS=M/N*abc
物体的质量m＝MΔm＝M/N*ms＝M/N*ρabc
对于非均匀物体，质量元Δmj＝ρjΔV＝1/N*ρjVs
其中ρj＝ρ（xi，yi，zi），1<=j<=N，但j只取位于物体内的随机点的i值。所以物体的质量

【例】一个球体半径为R＝0.5m球上有一个半径r=0.3m的圆柱形空洞，其轴线与球的直径相重合设球体的密度ρ＝1kg/m3，试用mc方法来求实体的体积
【解】图略。令球的中心位于坐标系的原点O处，做边长为2R＝1m的正方体，其中心和球心相重合，则正方体的体积Vs＝1m3，产生一组随机数（xi,yi,zi），他们的值域均为        [-0.5m,0.5m]。然后判断改随机点是否再实体内，判断依据是

若******生了N组随机数，而满足上述判据的有M组，那么球体的质量m＝M/N*VSkg，为了提高测量精度，可重复上述计算过程，例如用同样的方法做3次，取平均值作为结果，并求出其方均根偏差
4、链式反应的模拟
放射性物质的链式反应式一个随机过程，可借助计算机用mc方法来模拟和研究。由原子核物理学可以得知U235的原子核本质式不稳定的，会自发的发生裂变。裂变的激烈程度可用放射性物质的半衰期来描述，半衰期是指大量核中有1/2的核发生裂变所需要的时间。U235的半衰期为7亿多年。因此任何时刻发生裂变的核只是相对很小的一部分，其释放的能量只能使其本身微微温热。但是，再一定条件下，自发裂变放出的两个种子轰击其他U235核而被吸收，引起新的裂变，依次类推，可迅速释放出大量的能量，甚至引起爆炸，这就是链式反应。
设开始有N个U235核发生裂变，每个核放出两个中子，称为第一代中子，共有2N个。2N个中子有感生新一轮的裂变，产生第二代中子，为4N个。……如此下去，直至第n次裂变，产生n代中子为2nN个，如此计算，30代可产生裂变的核数为2^30N＝10N亿，即为第一次裂变核数的10亿倍，现在的问题是，在什么条件下才能产生链式反应呢？其基本要求是裂变产生的两个中子至少有一个能使第二个铀块发生裂变。为此要求核材料中杂质的含量，包括U238的含量要足够小。以避免中子被其他杂质或U238吸收。另外，由于热中子使U235里裂变的机会很大，所以在U堆里面还要加入减速剂，如重水或石墨等，以使快中子减速到热中子，最后十分重要的条件是U堆的体积要足够大，以避免裂变所放出的中子过多的未与U核相遇而飞出U块体外。这就涉及到临界体积和临界质量的概念。所谓临界质量是指可裂变物质能发生自续链式反应的最小质量。由于U核体积很小，一U核裂变放出的中子在和另外一块U核作用并使之发生裂变之前，平均来说要金国一定相对很长的距离，约为厘米数量级，因此，假定有N0个核自发裂变而放出2N0个中子，其中有N个中子在U块中引起另外的核发生裂变，其余的中子未与其他核碰撞而飞出U块。为描述一次裂变能引起下次裂变的众寡程度，定义裂变过程的倍增系数
k=N/N0
不难理解，维持自续联是反应的条件是N>=N0
即k>= 1
倍增系数k=1是临界质量Mc的条件。k的值与前面论技的诸多因素有关，本节将只讨论k核铀块的质量形状有关的问题。用计算机程序来模拟具有一定大小和形状的铀块众大量随机的裂变过程，统计算出相应的倍增系数k。
设铀块为长方体a*a*b，发生裂变的游客位于铀块内随机点（x0,y0,z0）处如图所示。随机点坐标的值域x0 [－0.5a,0.5a]；y0 [-0.5a,0.5a]；z0 [-0.5b,0.5b]；该核子裂变反应产生两个中子，其运动方向可用两个角坐标θ和φ来表示。如图。释放的每个中子按飞行方向的几率均等来考虑，或者说中子飞行方向的几率是按以（x0,y0,z0）为定点的立体角均匀分布的。立体角元可表示为
dΩ＝sinθdθdφ=-dφdcosθ

可见立体角均匀分布是按照φ角均匀分布和按cosθ均匀分布，而并非按θ角均匀分布。一次对应的两个随机数的值域为φ(0,2pi)，cosθ(-1,1)，平均的说能否击中另一个核只取决于中子在铀块体内飞行的距离，假设在0－1cm距离之间经过任何一段相同距离击中另一铀块的几率均等。或者说，中子在击中另一铀块之前飞行的距离为0～1之间均匀分布的随机数。一次与飞行距离相应的随机数为0<d<1，由此可计算出倍击中的铀核的位置
x1=x0+dsinθcosφ
y1=y0+dsinθcosφ
z1=z0+dcosθ
最后家查计算处的碰撞点（x1,y1,z1）是否位于铀块体内。若在铀块体内，则累计入引起新裂变中子数N0按照上述原则，归纳计算k的具体步骤为：
1.给定铀块质量M、铀块边长比s＝a/b和用于计算k的随机自发裂变核的个数，即旧裂变核个数N0，并设所选约化单位可使铀块的密度为1，体积为V，则得
M＝V=a2b=a3/s=b3s2
或a=(Ms）1/3
b=(Ms-2)1/3
2.产生九个0～1之间得随机数
r1,r2,r3……r9
3.旧裂变核得位置
x0=a(r1-0.5)
y0=a(r2-0.5)
z0=b(r3-0.5)
4.旧裂变放出得两个中子的方向
φ＝2pi×r4
cosθ1＝2r3-1
φ2＝2pi×r6
cosθ2＝2r7－1
5.中子飞行距离
d1=r8
d2=r9
6.可能发生新裂变得位置
x1=x0+d1sinθ1cosφ1
y1=y0+d1sinθ1sinφ1
z1=z0+d1cosθ1
和
x2=x0+d2sinθ2cosφ2
y2=y0+d2sinθ2sinφ2
z2=z0+d2cosθ2
7.检查上述位置（x1,y1,z1)和（x2,y2,z2)是否在铀块体内。如果
-0.5a<=x1<=0.5a
-0.5a<=y1<=0.5a
-0.5a<=z1<=0.5a
均满足，则N的值增加1。同样如果
-0.5a<=x2<=0.5a
-0.5a<=y2<=0.5a
-0.5a<=z2<=0.5a
也满足，N的值也增加1
8.重复步骤2～7，共执行N0次，然后计算
k=N/N0
若计算结果k<>1，则调整M和s的值后在进行上述步骤，直至k=1此时M即为临界质量。显然临界质量和s有关，与核材料的形状有关。

DoubleX

5、趋向平衡态
1。宏观过程的方向性
气体的自由膨胀是表征宏观过程方向性的典型例子。如图所示，一个密闭的盒子ei 分割成A和B两个体积相等的空间。起初A空间充满某种气体，B空间为真空，打开隔板后气体分子将充满整个盒子空间。这是一个自动进行的宏观过程，而相反的过程确从来未见过自动进行过。因此这一过程称为不可逆过程。不可逆扩成是各式各样的，但它们是相顾关联的。热力学第二定律的各种叙述都解释孤立系统内宏观过程的方向性。
在热力学众，任一宏观状态所包含的微观状态数，定义为该宏观状态的热力学概率，用Ω来表示。热力学概率是分子热运动混乱程度的量度。热力学系统的另一个状态函数熵S的微观意思也是分子热运动的无序性的量度。波尔兹曼关系
S=klnΩ
给出了熵和热力学概率之间的数量关系，其中k是波尔兹曼常数。
熵增加原理指出，孤立系统的熵永远不会减少。换言之，孤立系统内所发生的宏观过程总是朝着熵增加的方向进行，即对于孤立系统内的宏观过程
ΔS>0
由波尔兹曼关系看出，孤立系统内所发生的过程，总是由热力学嘎率小的宏观状态相热力学概率大的宏观张太进行，或者说朝着无序性增大的 方向进行；也就是说，总是由包含微观状态数少的宏观状态向包含微观状态数多的宏观状态进行，总之是朝着平衡态进行。熵增加原理是热力学第二定律的数学表示。他们均建立在统计力学的基础上。因此，孤立系统内必须包含大量分子。
2。趋向平衡态模拟概述
现仍利用上图的盒子来进行趋向平衡态的模拟。设有N个分子，其编号为1，2……，N。开始时所有分子均在A空间内，然后随机的 抽取某一编号的分子，令其穿过隔板进入另一个空间。当然，第一次抽取的分子肯定是由A进入B空间。以后每次抽取的分子则可能是由A进入B空间，也可能是由B进入A空间。令n为第x次抽取分子后B众的分子数。在
下次抽取时，选中B中分子而由B空间进入A空间的几率

PA=n/N
而选中A空间的分子由A进入B空间的几率
PB＝1－n/N
则B空间中分子数的增量
Δn＝PB-PA=1-2n/N
考察Δx个分子时，可近似利用这一关系，有
Δn＝（1－2n/N）Δx
用 连续量替代离散量可得
dn/(1-2n/N)=dx
作积分
得B空间中分子数占总分子数得比率
n/N＝1/2（1－e－2x/N)
可见当x→时,n/N＝1/2，即盒内气体 处于平衡态
3。趋向平衡态模拟程序流程
趋向平衡模拟程序得流程图可如下，其中N是分子总数，MX是抽取分子得总次数，B是空间B中的分子数。数组NU（I）＝1，（I＝1，2……,N)，表示出世状态错有N个分子均在空间A中，而如果某NU(I)=－1，则表示编号为I的分子在空间B中。因此若抽取到编号K的分子，只要将对应的NU(K)值取反即可表示该分子穿过隔板。RAND(I)代表值域为[0，1]的随机数发生器。INT(X)表示将X按去尾法取整。
本帖转自52mc论坛