蒙特卡洛方法在高分子材料中的应用

第六章高分子科学中的MonteCarlo方法MonteCarlo方法——一个十分独特的名字Monte-Carlo,MonacoMonteCarlo原为地中海沿岸Monaco的一个城市的地名，气候温和，景色怡人，人口不到一万，是世界闻名的大赌场。将MonteCarlo作为一种计算方法的命名固然已经赋予了新的内容。然而，顾名思义，MonteCarlo方法的随机抽样特征在它的命名上得到了反映。MC方法的发展归功于核武器早期工作期间LosAlamos（美国国家实验室中子散射研究中心）的一批科学家。vonNeumann,Metropolis,Ulam和Kahn等人在电子计算机上对中子行为进行随机抽样模拟，通过对大量中子行为的观察推断出所要求算的参数。LosAlamos小组的基础工作刺激了一次巨大的学科文化的迸发，并鼓励了MC在各种问题中的应用。学术界一般将Metropolis和Ulam在1949年发表的论文作为MonteCarlo方法诞生的标志。6.1MonteCarlo方法的基本思想MonteCarlo方法在数学上称其为随机模拟(randomsimulation)方法、随机抽样(randomsampling)技术或统计试验(statisticaltesting)方法．它的最基本思想是：为了求解数学、物理及化学等问题，建立一个概率模型或随机过程，使它的参数等于问题的解；当所解的问题本身属随机性问题时，则可采用直接模拟法，即根据实际物理情况的概率法则来构造MonteCarlo模型；然后通过对模型或过程的观察抽样试验来计算所求参数的统计特征，最后给出所求解的近似值。在高分子科学中的MonteCarlo模拟主要采用直接模拟方法。MonteCarlo方法的突出特点是，它的解是由试验得到的，而不是计算出来的。其程序结构简单，解题时受问题条件限制的影响较小，具有广泛的适应性。但不能解决精确度要求很高的问题。蒙特卡洛方法需要大量的随机数，计算量很大，人工计算需耗费大量的时间，利用计算机可大大减少计算时间，增加试验次数以提高计算精度，因此，蒙特卡洛方法的广泛应用与计算机技术的发展是不可分割的。设所要求的量x是随机变量ξ的数学期望E(ξ)，那么用MonteCarlo方法来近似确定x的方法是对ξ进行N次重复抽样，产生相互独立的ξ值的序列ξl，ξ2，…，ξN，并计算其算术平均值：根据Kolmogorov的大数定理则有：11NNiiN(lim)1NNPx即当N充分大时，成立的概率等于1，亦即可以用作为所求量x的估算值。()NExN例6-1用统计试验方法求圆周率π考虑边长为1的正方形，以其一角为圆心和边长为半径，在正方形内画一条1／4圆弧，如图所示。在正方形内等概率地产生n个随机点(xi，yi)，i=l，2，3…，n，设n个随机点中有k个点落在四分之一圆弧内，显然，当n→∞时有以下关系成立：因而，圆周率π的估值为：22/44krnr四分之一圆面积正方形面积4ˆkn判断随机点(xi，yi)是否位于圆内的判别式为：221iixy用一对(0，1)随机数Ul，U2分别模拟随机变量的取值xi和yi，当时，则计数器k值增1。这个判别式就是蒙特卡洛方法的概率模型。当试验次数n足够大时，所得的估值的精度也随之提高。22121UU例6-2.蒲丰氏问题ComtedeBuffon(1707-1788)FrenchNeedleexperiment,1777Buffon投针问题：平面上画很多平行线，间距为a。向此平面投掷长为l（l＜a）的针,求此针与任一平行线相交的概率p。alP2可以证明求出π值)(22nNalaPl其中Ｎ为投计次数，n为针与平行线相交次数。这就是古典概率论中著名的蒲丰氏问题。一些人进行了实验，其结果列于下表：试验者时间(年)针长投针次数相交次数π的估计值Wolf18500.80500025323.15956Smith18550.60320412183.15665Fox18840.7510304893.15951Lazzarini19250.83340818083.141592926.2MonteCarlo方法与高分子科学MonteCarlo模拟与高分子科学结下了不解之缘是由于高分子科学本身的特点所决定的，因为在高分子科学中存在着大量可供进行MonteCarlo直接模拟的随机性问题。如：由于聚合反应本身的随机性特点，高分子系综内各个成员之间存在着与其生成机理密切相关的特定分布，即体系中所生成的高分子链并非具有相同的分子量，而是存在着所谓的分子量分布问题；在多元聚合中，多元共聚物不仅具有分子量分布，而且导致了不同种单元在高分子链上的排列问题，即所谓的序列分布；在多官能团的聚合反应中的支化和凝胶化问题；高分子链的热降解和辐射降解等等，无一不是随机性问题。MonteCarlo方法在现代高分子科学中的应用主要具有以下特征:•由于高分子凝聚态物理的发展，高分子体系的MonteCarlo研究从对单链的研究转向对高浓度多链体系的研究。•由静态平衡态问题向动态和非平衡态问题发展也是当前高分子MonteCarlo模拟的重要特征。高分子链的分子运动学，尤其是高浓度多链体系的分子运动问题是当前研究的重要方面。•人们对共混和嵌段共聚物的界面、高分子和液晶的界面、高分子链的吸附、晶态和非晶态的界面性质和相互扩散问题开展了MonteCarlo模拟研究。•高分子MonteCarlo方法的新算法也是值得研究的。6.3随机数与伪随机数产生均匀分布随机数的方法可以采用物理方法和数学方法。最简单的产生随机数的物理方法是掷骰子游戏；采用电学噪声的变化也可产生随机数。但物理方法产生随机数的“费用”很高，且速度慢。因此，实际应用的随机数一般均在计算机上采用数学方法来产生。用数学方法产生的随机数一般均采用某种确定性的表达式来实现，因此其并非真正的随机，故通常称其为“伪随机数”。用数学方法产生伪随机数的优点是因为它借助于迭代公式，所以特别适合于计算机。而且其产生的速度快、费用低。目前，多数的计算机均附带有“随机数发生器”。用数学迭代方法产生的随机数存在两个问题：1、递推公式和初始值a1、a2、…、ak确定后，整个随机数序列便被唯一确定下来。即任意一个随机数被前面的随机数唯一确定了，不满足随机数相互独立的要求。2、既然随机数序列是用递推公式确定的，而在计算机上所能表示的[0，1]上的数又是有限多的，因此这样的随机数序列就不可能不出现重复地无限继续下去。这种随机数序列出现周期性的循环现象是与随机数的要求相矛盾的。对第一个问题不能从本质上改变，但只要递推公式选得好随机数的相互独立性是可近似满足；第二个问题，则不是本质的，因为用MonteCarlo方法解任何问题时，所用随机数个数总是有限的，只要保证不超过伪随机数序列出现循环现象的长度即可。用数学迭代方法产生随机数均存在周期现象，随着迭代过程的不同，其效果也各不相同。一般满足下列要求的产生方法才可被认为是好的：(1)随机性和统计独立性要好；(2)容易在计算机上实现；(3)省时，存贮量小；(4)伪随机数的周期长。乘同余法乘同余法由Lehmer首先提出。由于采用乘同余法具有在计算机上容易实现、快速等优点，因此乘同余法已被广泛采用。乘同余法的迭代公式为，1(mod)nnxxM作为[0，1]区间上均匀分布的伪随机数序列。（给出初始值x0及参数λ、M）当周期很大时，可用/nnrxM一个简单的例子1n6(mod11),/116,11nnnxxrxM（）01,x1,6,3,7,9,10,5,8,4当时得到序列:,1,6,3.,2.....003,1,,1,3,9........3,2,2,1,3,9,5,42,6,7,16....,,.08..xx如果令得到序列:如果令得到序列：上面的例子中，第一个随机数生成器的周期长度是10，而后两个的周期长度只有它的一半。我们自然希望随机数的周期越长越好，这样得到的分布就更接近于真实的均匀分布。0(Mx在给定的情况下，随机数的周期与和初值种子)选择有关。M(2s或其他)λ(52k+1或其他)M(2s)λ(32k+1)x0周期2305112303171或任意奇数2s-2231~234513231~2333191或任意奇数2s-2235~239515234~2363211或任意奇数2s-2241~244517237~2393231或任意奇数2s-2245~248519240~2423251或任意奇数2s-2243~2453271或任意奇数2s-2246~2483291或任意奇数2s-2108+1234759411858823521011751108101174k+1(k≠1)15×108231-116807任意整数231-1表：乘同余法的参数及周期MonteCarlo方法的核心就是随机数的使用，因此计算机模拟结果的优劣将强烈地依赖于伪随机数的质量。•伪随机数的均匀性•伪随机数的独立性对于已经产生的随机数质量的检验主要是：伪随机数的均匀性检验可用xn的矩来判别，均匀性好的随机数序列在N→∞时应满足下列要求：一阶矩二阶矩三阶矩四阶矩伪随机数独立性检验一般采用χ2检验。10111lim2NiNixxdxN1220111lim3NiNixxdxN1330111lim4NiNixxdxN1440111lim5NiNixxdxN随机变量的抽样：前面讨论了[0，1]均匀分布的伪随机数的产生，然而在实际应用中概率分布的形式是多种多样的。liilnlppp1)()0(),,2,1(,,,0并满足：niip11产生[0，1]随机数r，如果条件满足，则认为事件Ai发生。)()1(llprp一、从随机事件中抽样：假设随机事件的出现概率分别为Pi(i＝1，2，…n)。为了对随机事件Ai进行抽样，首先需构造累积概率：例6-3.掷骰子点数的抽样掷骰子点数X=n的概率为：选取随机数ξ，如则在等概率的情况下，可使用如下更简单的方法：其中［］表示取整数。61)(nXP661nnnXF1]6[FX二、连续型分布的抽样：连续型分布的一般形式如下：这里f(t)为分布的概率密度函数。如果分布函数的反函数存在，则连续型分布的一般抽样方法是通过其反函数直接抽样：这里r是[0，1]均匀分布的随机数，F-1为F(x)的反函数。()()xFxftdt1()Fr在［a,b］上均匀分布的分布函数为：例6-4.在［a,b］上均匀分布的抽样0()1xaxaFxaxbbaxb其抽样方法为：()abar这里r是[0,1]区间均匀分布的随机数。MonteCarlo方法的估值精度ε与试验次数N的平方根成反比，若精度提高10倍，则试验次数N要增加100倍。收敛速度慢是蒙特卡洛方法的主要缺点。蒙特卡洛方法的精度估算有概率性质，它并不断言精度一定好于ε，而只是表明，所算精度以接近于1的概率不超过某一界限，这是蒙特卡洛方法与其它确定性误差计算的根本区别之处。例6-5：中子扩散问题原子核反应堆的壁是铅制的，对中子起屏蔽作用。中子从反应堆内侧进入壁内与铅原子发生碰撞。求出穿透铅壁中子数的百分比，被吸收入铅壁中子数的百分比，以及重新返回反应堆中子数的百分比。1入口铅墙（长为3d）2d解：设壁厚为常量3d，中子是垂直进入壁内的，并设每个中子在壁内每次走过d（平均自由程）才与铅原子碰撞，碰撞后以随机的方向弹射，再走过d的距离，和第二个铅原子碰撞，如此继续下去。最后，有三种情况（1）中子穿透铅壁；（2）被铅壁吸收（假定经过8次碰撞后，没有穿透或返回，则认为被吸收；（3）重新返回反应堆。现在研究对中子运动的模拟：假设一个中子在壁内处于与壁内侧距离为x的位置上与铅原子碰撞，然后以θ角的方向弹射，那么θ是[0，2π]之间的均匀分布的随机数。中子经过弹射后，与壁内侧的距离x