计算物理MC_first

1、 Monte CarloStatistic of random variables Metropolis algorithm Implementation of the Metropolis MC method2对统计力学体系进行计算机模拟时，需要确定对统计力学体系进行计算机模拟时，需要确定体系的位形（组态）。按照产生位形变化的方体系的位形（组态）。按照产生位形变化的方法，可以将计算机模拟分成两大类：法，可以将计算机模拟分成两大类：一类是一类是( ( stochastic) )方法：方法：Monte Carlo: MC; Simulated annealing 一类是一类是( (determi

2、nistic方法：方法： 按照体系的动力学规律产生位形变化。按照体系的动力学规律产生位形变化。Molecular Dynamics: MDMonte Carlo方法（方法（MC）亦称为随机模拟（）亦称为随机模拟（Random Simulation）、随机抽样）、随机抽样(Random Sampling) 或统计试或统计试验（验（Statistical Testing）方法。）方法。Monte Carlo 方法名称的固定、得到系统发展大方法名称的固定、得到系统发展大约是二十世纪四十年代；但如就其方法特征可追约是二十世纪四十年代；但如就其方法特征可追溯到十八世纪著名的溯到十八世纪著名的Buffon

3、问题问题-用随机投掷用随机投掷缝针试验计算圆周率数值。缝针试验计算圆周率数值。The Comte de Buffon needle experiment, AD 1777dLp2dLNMdLp2MdNL2试验方案是：在平面内划一组相距试验方案是：在平面内划一组相距 d 的平行线，向此的平行线，向此平面随意的投掷长度平面随意的投掷长度 L 的细针，那末从针与平行线相的细针，那末从针与平行线相交的概率交的概率 p 可以得到可以得到 值。值。the Comte de Buffon needle experiment实验者实验者时间时间年份年份针长针长投针投针次数次数相交相交次数次数的的估计估计Wol

4、fWolf185018500.80.850005000253225323.15963.1596SmithSmith185518550.60.6320432041218.51218.53.15543.1554De Morgan CDe Morgan C186018601.01.0600600382.5382.53.1373.137FoxFox188418840.750.75103010304894893.15953.1595LezzeriniLezzerini190119010.830.8334083408180818083.14159293.1415929Reina19250.54192520

5、8593.1795历史上一些有名用投针试验计算历史上一些有名用投针试验计算值的结果值的结果Monte CarloMonte Carlo方法以概率统计为理论基础，方法以概率统计为理论基础，以随机抽样为主要手段。其基本思想是首先以随机抽样为主要手段。其基本思想是首先建立一个概率（或随机过程）模型，使它的建立一个概率（或随机过程）模型，使它的参数等于问题的解；然后通过对模型（或过参数等于问题的解；然后通过对模型（或过程）的抽样试验来获取有关参数的统计特征程）的抽样试验来获取有关参数的统计特征、解的近似值及精度估计。、解的近似值及精度估计。Monte Carlo的基本思想的基本思想NnnNN11设所要

6、求的量设所要求的量x是随机变量是随机变量的数学期望的数学期望E()，那么，那么用用Monte Carlo方法来近似确定方法来近似确定x的方法是对的方法是对进行进行N此重复抽样，产生相互独立的此重复抽样，产生相互独立的值的序列值的序列1，2，,N并计算其算术平均值：并计算其算术平均值：根据根据Kolmogorov的大数定理则有的大数定理则有即当即当N充分大时，充分大时， 成立的概率等于成立的概率等于1.亦即可以用亦即可以用 作为所求量作为所求量x的估计值。的估计值。1)lim(xPNNxEN)(NMeasuring ErrorN0Variance 方差方差: 2 = Var(X) = - 2St

7、andard deviation (standard error) of a random number is the square root of its variance2 = Var(X) = - 2Variance 实验者实验者时间时间年份年份针长针长投针投针次数次数相交相交次数次数的的估计估计WolfWolf185018500.80.850005000253225323.15963.1596SmithSmith185518550.60.6320432041218.51218.53.15543.1554De Morgan CDe Morgan C186018601.01.0600600

8、382.5382.53.1373.137FoxFox188418840.750.75103010304894893.15953.1595LezzeriniLezzerini190119010.830.8334083408180818083.14159293.1415929Reina19250.541925208593.1795历史上一些有名用投针试验计算历史上一些有名用投针试验计算值的结果值的结果N/1显然这一早期的显然这一早期的“古典古典Monte Carlo方法方法”及相应的及相应的抽样实践已寓示与这种数值计算方法相伴随的巨大抽样实践已寓示与这种数值计算方法相伴随的巨大工作量。所以可以理解

9、直到近半个世纪以来，工作量。所以可以理解直到近半个世纪以来，Monte Carlo方法的应用范围才不断扩展，形成计算方法的应用范围才不断扩展，形成计算数学的一个重要分支，是与电子计算机的相应发展数学的一个重要分支，是与电子计算机的相应发展不可分割的。不可分割的。 如果使如果使 值的精度达三位有效数字，则需数十万次。值的精度达三位有效数字，则需数十万次。实验者实验者时间时间年份年份针长针长投针投针次数次数相交相交次数次数的的估计估计WolfWolf185018500.80.850005000253225323.15963.1596SmithSmith185518550.60.6320432041

10、218.51218.53.15543.1554De Morgan CDe Morgan C186018601.01.0600600382.5382.53.1373.137FoxFox188418840.750.75103010304894893.15953.1595LezzeriniLezzerini190119010.830.8334083408180818083.14159293.1415929Reina19250.541925208593.1795历史上一些有名用投针试验计算历史上一些有名用投针试验计算值的结果值的结果N/1Monte Carlo模拟就是边产生随机数，边在计算机上进行随模

11、拟就是边产生随机数，边在计算机上进行随机过程模拟。机过程模拟。步骤如下步骤如下运行运行（1）在计算机上产生随机数，并使随机数）在计算机上产生随机数，并使随机数游动。游动。（2）进行条件判断，在满足条件的情况下）进行条件判断，在满足条件的情况下，按规则进行作业。，按规则进行作业。反复操作（反复操作（1）、（）、（2）。）。The Name of the GameMetropolis coined the name “Monte Carlo”, from its gambling Casino.Monte-Carlo, MonacoMonte Carlo 是位于欧洲地中海沿岸摩纳哥是位于欧洲地中海

12、沿岸摩纳哥（Monaco）国的一个重要城市。）国的一个重要城市。Stanislaw Ulam (1909-1984)S. Ulam is credited as the inventor of Monte Carlo method in 1940s, which solves mathematical problems using statistical sampling.大数定理是大数定理是Monte CaloMonte Calo模拟的理论基础，模拟的理论基础，我先做一下简单的介绍我先做一下简单的介绍大数定理大数定理 大数定律大数定律大数法则大数法则概率论历史上第一个极限定理属于伯努利，概率

13、论历史上第一个极限定理属于伯努利，后人称之为后人称之为“大数定律大数定律”。是概率论中讨论。是概率论中讨论随机变量序列的随机变量序列的算术平均值向常数收敛算术平均值向常数收敛的定的定律。是概率论与数理统计学的基本定律之一。律。是概率论与数理统计学的基本定律之一。又称弱大数理论。又称弱大数理论。数学家伯努利数学家伯努利 主要含义主要含义 有些随机事件无规律可循，但不少却是有规律有些随机事件无规律可循，但不少却是有规律的，这些的，这些“有规律的随机事件有规律的随机事件”中在大量重复中在大量重复出现的条件下，往往呈现几乎必然的统计特性，出现的条件下，往往呈现几乎必然的统计特性，这个规律就是大数定律。

14、这个规律就是大数定律。通俗地说，这个定理就是，在试验不变的条件通俗地说，这个定理就是，在试验不变的条件下，重复试验多次，随机事件的下，重复试验多次，随机事件的频率频率近似于它近似于它的的概率概率。比如，我们向上抛一枚硬币，硬币落下后哪一面朝上本比如，我们向上抛一枚硬币，硬币落下后哪一面朝上本来是偶然的，但当我们上抛硬币的次数足够多后，达到来是偶然的，但当我们上抛硬币的次数足够多后，达到上万次甚至几十万几百万次以后，我们就会发现，硬币上万次甚至几十万几百万次以后，我们就会发现，硬币每一面向上的次数约占总次数的二分之一。这种情况下，每一面向上的次数约占总次数的二分之一。这种情况下，偶然中包含着必然

15、。必然的规律与特性在大量的样本中偶然中包含着必然。必然的规律与特性在大量的样本中得以体现。得以体现。 又如称量某一物体的重量，假如衡器不存在系统偏差，又如称量某一物体的重量，假如衡器不存在系统偏差，由于衡器的精度等各种因素的影响，对同一由于衡器的精度等各种因素的影响，对同一物体重复称物体重复称量多次量多次，可能得到多个不同的重量数值，但它们的算术，可能得到多个不同的重量数值，但它们的算术平均值一般来说将随称量次数的增加而逐渐接近于物体平均值一般来说将随称量次数的增加而逐渐接近于物体的真实重量。的真实重量。 伯努利大数定律：伯努利大数定律：设设 n 为为 n 重伯努利实验中事件重伯努利实验中事件

16、 A 发生的次数，发生的次数，p 为每次实验中为每次实验中 A 出现的概率，则对任意的出现的概率，则对任意的 0 有（有（2）成立。）成立。 limn1pnPn(2)这一定理指出，不论随机变量的分布如何，只这一定理指出，不论随机变量的分布如何，只要要n n足够大，事件足够大，事件A A出现的频率趋向一个稳定的出现的频率趋向一个稳定的值，这个值就是值，这个值就是A A的概率。的概率。What is a Random Number? Follow a definite distribution, usually uniform distribution Uncorrelated Unpredict

17、ableThe Roulette and DiceMechanical random number generators运用运用Monte Carlo方法首先要产生一组随机数列。随机数列方法首先要产生一组随机数列。随机数列的产生方法，一种是找一个实际的随机物理过程，譬如说的产生方法，一种是找一个实际的随机物理过程，譬如说掷骰子，纪录放射辐射源中两个衰变之间的时间值等等。掷骰子，纪录放射辐射源中两个衰变之间的时间值等等。另一种方法是由计算机自己产生随机数列，由于这些随机另一种方法是由计算机自己产生随机数列，由于这些随机数不是从实际过程中得来的，故称之为数不是从实际过程中得来的，故称之为“伪随机伪

18、随机数数” (Pseudo-random numbers) 。产生随机数列的程序称。产生随机数列的程序称为随机数产生器为随机数产生器(random number generator)Pseudo-Random Numbers & random number generator伪随机伪随机数数 Truly random numbers can not be generated on a computer Pseudo-random numbers follows a well-defined algorithm, thus predictable and repeatable Pseudo-ra

19、ndom numbers have nearly all the properties of true random numbersPseudo-Random Numbers随机数可以有各种分布形式，它们都可以从随机数可以有各种分布形式，它们都可以从均匀分布均匀分布随随机数列中得出。而产生均匀分布随机数列的数学方法又有机数列中得出。而产生均匀分布随机数列的数学方法又有很多种，下面介绍几种。很多种，下面介绍几种。随机性和统计独立性要好随机性和统计独立性要好容易在计算机上实现容易在计算机上实现效率高效率高省时、占内存小省时、占内存小周期足够大周期足够大（1）平方取中）平方取中给定第一个数以后，则序

20、列中每个数可由前一个数按上给定第一个数以后，则序列中每个数可由前一个数按上式求得。式求得。0.6406, 0.0368, 0.135464061x4103683621xSnnx21110/把把2S位的数自乘，去头截尾只保留中间的位的数自乘，去头截尾只保留中间的2S位。位。13543x03682x0013542422x0183331623x)10),10/mod(int(221SSnnxx135403686406321xxx018333160013542441036836232221xxxSnnx21110/mod (a, p)a - INT(a / p) * pmod (7,3)3681000

21、0*)0368.41int(410368)10000,410368mod()10),10/mod(int(221SSnnxx)2),2/mod(int(221SSnnxxSnnx2112/平方取中平方取中)10),10/mod(int(221SSnnxxSnnx21110/)2),2/mod(int(212SSnnnxxxSnnx2222/乘积取中乘积取中二进位二进位)10),10/mod(int(221SSnnxxSnnx21110/)10),10/mod(int(212SSnnnxxxSnnx22210/乘积取中乘积取中平方取中平方取中十进位十进位（2）乘同余法）乘同余法(linear c

22、ongruential method)的叠代的叠代公式为公式为),mod(1cbaxxnn给定第一个数以后，则序列中每个数可由前一个数按上给定第一个数以后，则序列中每个数可由前一个数按上式求得。式求得。其中其中a, b and c are magic numbers: the choice of these numbers determines the quality of the generator.214748364712, 0,168077315cbacxranfn/()1MOD (a, p)a - INT(a / p) * pChoice of ParametersNameca (mu

23、ltiplier)bperiodANSIC rand()231110351524512345231Park-MillerNR ran0()231-1168070231-2drand48() 2482521490391711248Hayes 64-bit26463641362238467930051264(a x + b) mod cShort-Coming of LCGxnWhen (xn, xn+1) pairs are plotted for all n, a lattice structure is shown.xn+1),mod(1cbaxxnn1,1282, 0, 507xcbacx

24、ranfn/()1iiryrx212nnxx168071cxranfn/()1)2,2(313121474836471231c111cxxnn1*20iseedxIA/ICIQ )IA,IC(MODIR 12, 1 Iseed31Although popular, by virtue of the ease with which it can be programmed the linear congruential method does not satisfy all of the requirements that are now regarded as important in a r

25、andom number generator.随机性和统计独立性要好随机性和统计独立性要好容易在计算机上实现容易在计算机上实现效率高效率高省时、占内存小省时、占内存小周期足够大周期足够大mnnryrxxnxn+1（4） The Marsaglia random number generatorCombination generator：A period 2144.I) Lagged Fibonacci generator1yxyxyxyx snrnnxxx ii) Arithmetic sequence method16777216/16777213dcdcdcdc 16777216/7654

26、321cc1nn nnncxU Numerical integrationMonte Carlo integration1 定积分问题定积分问题 。dxxI10)exp(这个积分的值这个积分的值 I = - (e-1-1)0.63212为了用随机抽样方法来求解该积分为了用随机抽样方法来求解该积分, 我们先要构造一我们先要构造一个概率模型个概率模型. 设所要求的量设所要求的量x是随机变量是随机变量的数学期望的数学期望E()，那么用，那么用Monte Carlo方法来近似确定方法来近似确定x的方法是的方法是对对进行进行N次重复抽样，产生相互独立的次重复抽样，产生相互独立的值的序列值的序列1，2，,

27、N并计算其算术平均值：并计算其算术平均值：NnnNN11dxxI10)exp(图中方框的总面积为图中方框的总面积为1，而我们所求的积分值，而我们所求的积分值I即为图中即为图中阴影部分的面积。在正方形平面中均匀随机的投点，则阴影部分的面积。在正方形平面中均匀随机的投点，则落入阴影区的概率则为积分值落入阴影区的概率则为积分值I。NvI 掷点法掷点法Hit or miss method)exp( x00IENvI)(因此对于这里考虑的问题，我们拟构造这样的概率模型，即因此对于这里考虑的问题，我们拟构造这样的概率模型，即在在11的正方形平面中均匀随机的投点，则落入阴影区的概率的正方形平面中均匀随机的投

28、点，则落入阴影区的概率则为积分值则为积分值I。设在。设在N次投点试验中，落入阴影区的点为次投点试验中，落入阴影区的点为v次次,那那么观察频数么观察频数v也是一个随机变量，其数学期望也是一个随机变量，其数学期望 E(v)=NI。令。令表示观察频率，那么按照大数定理，当表示观察频率，那么按照大数定理，当N充分大时，频率收敛充分大时，频率收敛于概率，即以概率于概率，即以概率1成立。因此可由上述概率模型在成立。因此可由上述概率模型在N很大时很大时所得到的所得到的v/N等于所求的积分值等于所求的积分值I。图中方框的总面积为图中方框的总面积为1，而我们所求的积分值即为而我们所求的积分值即为图中阴影部分的面

29、积。图中阴影部分的面积。（1） 产生在产生在0，1区间上均匀分布独立的随机数区间上均匀分布独立的随机数r, r;NvI 计算机上的具体计算步骤为：计算机上的具体计算步骤为：（2） 令令r 和和r分别为所投点的分别为所投点的x, y坐标值，若坐标值，若re-r,则则表示所投点落在阴影区内，因此表示所投点落在阴影区内，因此v加上加上1, N也加上也加上1;若若re-r,则表示所投点落在阴影区外，因此则表示所投点落在阴影区外，因此v加上加上0,N加上加上1（3） 重复（重复（1），（），（2）直至）直至N足够大足够大.（4） 计算计算 以上称为一维定积分计算的掷点法。以上称为一维定积分计算的掷点法。

30、Nf1只要选取足够多的随机点，即当只要选取足够多的随机点，即当N充分大时充分大时)(iif)(NiifNI11（2）一维定积分计算的平均值法）一维定积分计算的平均值法若在若在x的定义域的定义域0,1上均匀随机的取点，即上均匀随机的取点，即选取随机数选取随机数，定义一个随机变量，定义一个随机变量1为为以概率以概率1成立。因此可由上述概率模型在成立。因此可由上述概率模型在N很大时得到所很大时得到所求的积分值求的积分值I。平均值法平均值法: Sample mean method)(NiifNI11)(11NiifNIiefi)(1) 产产 生在生在0，1区间上均匀分布独立的随机数区间上均匀分布独立的

31、随机数i。2) 令，令，, N加上加上1。3) 重复（重复（1），（），（2）直至）直至N足够大足够大.4) 计算计算一维定积分计算的平均值法具体计算步骤为一维定积分计算的平均值法具体计算步骤为dxxA1021Nf141412rMNA如果一维积分中的长度不是如果一维积分中的长度不是1，而是一般的情况（，而是一般的情况（a, b)， 那么，那么，f 在在(a, b)中中N个随机点上的平均值是对个随机点上的平均值是对dxxfba)(dyyfab10)()(作变换：作变换：x=a+(b-a)ydx=(b-a)dydxxf50)(NiiyfN1)(5dxxfba)(dyyfab10)()(作变换：作变

32、换：x=a+(b-a)yx=5yNiiyfNdxxf150)(5)(随机点在区域（随机点在区域（0, 1)中是均匀分布的。中是均匀分布的。dxxf52)(随机点在区域（随机点在区域（0,1)中是均匀分布的。中是均匀分布的。NiiyfN1)(3. 23 yxNf1x=5yNumerical Integration of Multi-dimensions 6152),(dxdyyxf NiiibafNdxdyyxf16152),(15),(AfdA用用Monte Carlo 方法计算积分方法计算积分, 最简单的办法最简单的办法就是在积分区域内随机地取一系列点就是在积分区域内随机地取一系列点, 计算

33、被计算被积函数在这些点上的数值并取平均积函数在这些点上的数值并取平均, 然后乘以然后乘以积分区域的积分区域的. 这种方法称为简单抽样法。这种方法称为简单抽样法。2/11 NNfSuppose we need to integrate from x0 to x1. We shall subdivide this interval into n steps of size x=(x1x0)/n as shown in figurex)f(xf(x)dxI(f)xxx000 xxxf(x)dxI(f)00 xxxdxxx)(xfxx)(xf)f(x00)( 21)(200000302006121x)

34、(xfx)(xfx)f(xError Nx/ 12N2/11 NNf1N Trapezium rulexxf(x)f(xf(x)dxI(f)xxx)210000 xxxxxxdxxx)(xf)f(xf(x)dxI(f)0000)(000)()21) 6121212130020000 xOxxf(x)f(xxx)(xfx)(xf)f(x)f(x（Error 复化梯型公式复化梯型公式 Compound Trapezium rulexxf(x)f(xf(x)dxI(f)xxx)210000The Compound Trapezium Rule approximation to the integra

35、l is thereforeWhile the error for each step is O(x3), the cumulative error is N times this or O(x2) O(N-2)2N Simpsons rulex0 x0+ xx0+2 x)()(4)(6bfcfafabS4NNumerical Integration of Multi-dimensionsIn general if the error goes as O(Na) in one dimension, then the error in d dimensions goes as O(Na/d).多

36、重积分的计算多重积分的计算物理量的平均值的计算归结为一个多重积分的计算物理量的平均值的计算归结为一个多重积分的计算, 设系统由设系统由100个粒子构成个粒子构成, 每个粒子有每个粒子有6个自由度个自由度, 所所以需要计算以需要计算300重积分重积分. 现在考虑在每一维取现在考虑在每一维取10个点个点, 总共有总共有10300个点个点. 假设计算机每秒可以计算假设计算机每秒可以计算1012个点个点, 计算计算 这个积分需要这个积分需要10288秒秒!问题比这个更严重问题比这个更严重! 如果如上述方式取点如果如上述方式取点, 则积分区则积分区域的内点数为域的内点数为 8300, 在总的点中所占比例

37、为在总的点中所占比例为 (8/10)300300, 也就是说也就是说, 取的点基本上都在表面上取的点基本上都在表面上!定积分问题定积分问题Monte Carlo一维定积分计算一维定积分计算 X 多重定积分计算多重定积分计算 goodRandom variables 随机变量随机变量Probability distribution functions (PDF).Expectation ValueStatistic of a Random VariableStatistic of random variables0,1,1,0,0,0,1,0,0,1,1,0,1Random variables

38、随机变量随机变量Random variables 随机变量随机变量532463451）在相同的确定实验条件下，对的观测无法给）在相同的确定实验条件下，对的观测无法给出单一固定值出单一固定值; ;）必须依据遍举测量原则，对所有可能取值给出）必须依据遍举测量原则，对所有可能取值给出发生概率发生概率Random variables 随机变量随机变量 X Random variables are hence characterized by a domain which contains all possible valuesthat the random value may take. This d

39、omain has a corresponding probability distribution functions (PDF).Random variables 随机变量随机变量Discrete variable 离散变量举例离散变量举例: :DomainTo this domain we have the corresponding probabilitiesiif there is a continuous range of values, such as an angle between 0 and 2 Continuous variable 连续型连续型随机变量随机变量: :在连

40、续区间取值，其取某确定值的概率由分在连续区间取值，其取某确定值的概率由分布密度函数给出布密度函数给出f(x)dx gives probability that X falls between x and x +dx.Random Variable A variable X that takes “random” values. We assume that it follows a probability distribution, p(x). Discrete variable: p1, p2, Continuous variable: p(x)dx gives probability th

41、at X falls between x and x +dx.11iip( )1p t dtProbability distribution functions (PDF).Expectation Value If the probability distribution is known, the expectation value (average value) can be computed as (for continuous variable)f()(f()f( )p( )XEXxx dxExpectation Value 期望值期望值3.5Expectation Value7Sta

42、tistic of a Random VariableVariance 方差方差: 2 = - 2Correlation 相关相关: C= - Standard deviation (standard error) of a random number is the square root of its variance随机变量随机变量x分布对期望值的离散程度分布对期望值的离散程度Variance 2 = - 2Variance 2 = - 22 =(1+4+9+16+25+36)/6-3.52=91/6-12.25=15.168-12.25=2.91VarianceVariance 2 = - 25.83Continuous Random VariablesFirst of all, if there is a continuous range of values, such as an angle between 0 and 2 , then the probability of getting exactly a speci