计算物理讲义2

  • 格式:pdf
  • 大小:911.11 KB
  • 文档页数:10

第一章 Monte Carlo方法基础 §1.1 随机数产生器 1-1第一章 Monte Carlo方法基础 Monte Carlo方法起源于二次大战中的原子弹设计研究。中子在非均匀介质(如反应堆芯或原子弹)中的散射和输运问题由一组复杂的积分-微分方程描述,当时Metropolis和Ulam等人采用了一种数值的实验方法来计算而非直接求解输运方程。为了保密起见,他们将该方法以著名的赌城-摩纳哥王国的Monte Carlo城(该城现在主要以F-1车赛的主办地之一而闻名)而命名(另一种说法是因为Ulam的叔叔每年要去Monte Carlo赌博)。实际上,它是采用随机数于各种物理计算和模拟实验,以类似于赌博中掷骰子的方法来随机决定其中某个单独事件的结果,因此模拟本身颇有游戏的色彩。对于模拟产生的大量随机事件发生后的最终结局应是服从基本物理规律和统计规律的,从而得到物理问题的正确答案。通常,这类物理问题本身相当复杂,很难用解析的方法进行分析和计算,而Monte Carlo方法特别适用于如基于统计力学和量子力学等领域中的复杂问题,是一很好的数值计算和数值模拟的方法,在凝聚态物理(表面物理、临界现象、非晶态)、应用物理(金属学、扩散和偏析)、理论物理(量子场论、统计理论、基本粒子、核物理)、化学(化学反应、高分子物理)以及非线性现象(分形、逾渗)等研究领域中起作非常重要的作用,成为科学研究的一种标准手段。除了用于研究数学物理问题以外,它的应用还已渗透到了经济和技术等其它各个领域,充分显示出其独到之处。 学习Monte Carlo模拟方法的意义,不仅在于它是科学研究的重要和普遍的手段,而且由于它的简单易操作性可以使学生更容易理解和掌握理论课程中的知识,如果再借助于计算机绘图方法,更使得一些深奥的理论概念成为直观上可以想象和理解的。

§1.1 随机数产生器 许多数值计算都需要用随机数来构造初始构型或直接用于数值模拟。尽管人们会认为计算机可以产生随机数,但实际上在计算机中没有什么东西是随机的,计算机是按照给定的程序严格执行的,如果程序起始是一样的话,其结果也应该是一样的,即可预测的。随机数产生器其实是一伪随机数产生器。而随机数的随机性是指单个随机数之间的关联很小;当用两个不同的随机数产生器产生的随机数序列代入某一实际的应用程序后,其造成的统计输出结果值应该是一致的。如果不是的话,至少其中一个随机数产生器是不好的。有一些统计测试标准可以用来判断随机数的质量。

1.1.1 均匀随机数产生器 1.1.1.1 伪随机数的要求 最有用的随机数是[0,1]区间内均匀分布的随机数。对于一个随机数产生器的第一章 Monte Carlo方法基础 §1.1 随机数产生器 1-2好坏,有三个判断标准,好的随机数产生器要求: 1、 一个随机数序列并非是永不重复的,在经过一个周期后,其序列将再次

形成一个循环重复。因此需要有长的循环周期。如在32位机上,要求周期接近最大可取的正整数31212147483647−=(整数的取值区间是31312,21⎡⎤−−⎣⎦),也可设计出周期更长的随机数序列。

2、 有良好的随机性,也就是说,随机数序列{}nx中的所有随机数之间没有

大的关联,或自相关函数()()1,2,3,Cll=…的值很小。说明这个关联的一个办法是在xy平面上绘出点(),nnlxx+。好的随机数产生器是,对于任何l值,产生的平面上数据点的分布都是非常均匀的。而差的情况是有条带、格点等非均匀分布出现。 3、 由于Monte Carlo计算中需要大量使用随机数,因此要求算法的速度非

常快。为了得到好的统计结果,随机数产生器的速度也是一个非常重要的考虑因素。

1.1.1.2 Lehmer线性同余法 最简单的均匀随机数的产生法是线性同余法,随机数序列{}nx是按线性关系 ()1modnnIaIbm+=+, (1.1.1.2-1)

nnxIm=, (1.1.1.2-2) 得到。整数12,,II…的取值是0至1m−中的一个值,故实数nx是严格小于1的,但偶尔会为0。也因为nx的取值个数为m,因此希望m足够大以使得[0,1]区间中的随机数分布足够稠密。这里的常数a是乘子,b是增量,m是模数,如何选择它们就决定了产生器的质量。0b=时称为乘同余法,0b≠时为混合同余法。 线性同余法的优点是速度快,每次调用只要几个操作步骤即可,因此是最为常用的。不好之处是序列相关性方面还差一些,如果用k个连续随机数在k维空间画点(一个点的k个坐标值)的话,那么这些点不会均匀充满k维空间,而是布满在1k−维平面上,形成一种类晶体结构(2k=时,用连续的两个随机数在xy平面上绘出点,则点是沿着直线段的),最多会有()1!kkm个平面。如果常数,,abm选择得不合适的话,则平面数远远小于此值。m的值常取为32位计算机所可能取的最大值

312∼。因此连续使用随机数于3维空间中,则

空间的平面数的上界为2344,所以当你关注的体积中的一小块时,平面的规则离散性对计算结果有很大影响。 例如,以前一个有名的子程序RANDU在IBM主流机型上使用了多年,并且移植到了多

种其他系统上,它采用的常数值是:3122147483648m==,162365539a

=+=,

图1.1.1.2-1 RANDU形成的3维晶体

结构。 第一章 Monte Carlo方法基础 §1.1 随机数产生器 1-30b=。如果以连续的3个随机数作为(),,xyz坐标画图的话,可以发现它产生的平面数仅为15,平面()965,4,,8,9xyzkk−+==−−󰀢以等间隔排列的(图1.1.1.2-1),这是因为随机数满足方程 {}2

2169655396655399nnnnxxxx++−+=−×+

3131220mod2nx=×= (1.1.1.2-3)

有趣的是,某个计算中心的顾问对此的解释是:“我们保证每一个随机数是单独随机的,但不能同时保证其它也是随机的”。 (1.1.1.2-1)式将自身循环,显见周期不会超过m。如果适当选择常数,,abm,可使周期有最大的长度m,此时,0至1m−之间的所有整数可出现一次,因此随意选择的任一初始值0I与其它初始值有相同的效果。从理论和实际效果看,0b=

和0b≠是差不多的。 对混合同余法,可取m为最大整数312,而a和b的选择要经过仔细研究,如当a满足()mod85a=,m满足100mamm<<−且二进制表示上没有特别的形态,b是奇数满足12360.21132bm−≈∼的情形下可以给出好的随机数。

1.1.1.3 16807产生器 有一个被认为是“最低标准”的产生器是: 531716807,0,212147483647abm====−=。 (1.1.1.3-1)

它已经通过了许多理论测试和大量实用的考验,被认为是32位机上最有效的产生器,可以作为其他产生器与之比较的标准,但它仍然不是十全十美的。 要把(1.1.1.3-1)式直接变成计算机程序并不是一件易事,因为32位机的整数区间是3131[2,21]−−,naI在取模之前就可能超过最大整数值。如果计算机可以自动取模的话,则编程就很容易了,只需给一初始值然后连续产生116807nnII+=即可,但这样的话,物理问题可能其值已超界但计算机仍能持续运行,故实际上计算机系统不会自动取模,而必须设计取模的方法。 Schrage方法:设m可表示为

,[/],modmaqrqmarma

=+==,

(1.1.1.3-2)

如:()2147483647168071277732836,127773,2836qr=×+==。则对:rq01zm<<−

,()modazq和[]rzq取值区间是[]0,1m−,有

()

[]()()[]()[]{}()[]()[]

modmodmodmodmod

modmodmod,0,mod,

zrzazmaqrm

qq

zqzqzqaqraqrrzqrm

qq

azqrzqmazqrzqifazqrzqmotherwise

⎧⎫=+−⎨⎬

⎩⎭

⎧⎫=+++−−⎨⎬

⎩⎭

=−

⎧−≥

=⎨

−+⎪⎩

。 (1.1.1.3-3) 第一章 Monte Carlo方法基础 §1.1 随机数产生器 1-4另外一种简便但不严格的编程方法是直接利用整数的区间。当4字节的整数naI超过最大值3121−时,计算机将保留后32位有效位数,而自动扔掉前面多余

的位数(因此,随机数序列不等同于Schrage方法得到的结果),剩下的位数中最左边的第一位决定了该整数的符号(0为正,1为负)。该位有一半的可能性为1,即整数值成为负的。这时可取绝对值,或加312使它变为正的,但312已经超

过最大值,故要加上31212147483647−=再加1,最后再乘以311024.65661210−−=×

以得到[]0,1区间中均匀分布的随机数。最后,应取奇数种子值使整数值恒为奇数。

[作业]:用Schrage方法编写随机数子程序,用连续两个随机数作为点的xy坐标

值绘出若干点的平面分布图;再用下面的2χ统计测试其2维独立性(总点数N大于7210×)。再取106,1288,6075abm===,重复以上计算,总结和评价产生器的随机性质量。

1.1.1.4 种子值 在用同余法时,有时要求每次使用随机数序列时初始的种子值也是不一样的,因此需要一个系统产生种子值的方法,几乎每个计算机系统都有一个报告当前时间的整数值,可以用它来构造一个初始值。例如:大多数计算机上都有 年: 099yi≤≤, 月: 112mi≤≤, 日: 131di≤≤ 时: 023hi≤≤, 分: 059ni≤≤, 秒: 059si≤≤ 则可设种子值为:(){}()07012312359ymdhnsIiiiiii=+++++⎡⎤⎣⎦,它的值约在区间310,21⎡⎤−⎣⎦内,第二部分的括号在100年内不会重复。

1.1.1.5 Tausworthe位移计数器法 Tausworthe方法或称一般反馈位移计数器法可以用以消除同余法中的关联

问题。人们曾经认为它要强于16807法,但后来发现情况并非那么简单,但至少该方法提供了随机数产生的另一种途径。该方法中即对整数进行位操作:首先用其它方法产生一个随机的整数序列,然后对两个整数进行XOR(与或)操作以产生一个新的随机整数(需再除以m得到[]0,1区间的随机数):

nnpnqIII−−=⊕

(1.1.1.5-1)

其中的[],pq是一对整数,最佳的选择是满足条件 221pqprime++=, (1.1.1.5-2)

的Mersine素数,如:[31,3]、[98,27]、[250,103]、[1279,(216,418)]、[9689,(84,471,1836,2444,4187)]等,其中的R250 (250p=,103q=) 是最为常用的产生