蒙特卡罗方法课件

蒙特卡罗方法湘潭大学物理学专业基础课计算物理及其应用

材料与光电物理学院孟利军2010年4月计算机模拟随机性模拟方法确定性模拟方法直接蒙特卡罗模拟蒙特卡罗积分Metropolis蒙特卡罗模拟蒙特卡罗通过不断产生随机数序列来模拟过程通过数值求解一个个粒子的运动方程来模拟整个系统的行为计算机模拟和蒙特卡罗方法间接蒙特卡罗模拟目录第一节蒙特卡罗方法概述第二节随机数与伪随机数第三节随机变量的抽样第四节蒙特卡罗方法的应用实例§1蒙特卡罗方法概述---基本思想基本思想：

针对待求问题，根据物理现象本身的统计规律，或人为构造一合适的依赖随机变量的概率模型，使某些随机变量的统计量为待求问题的解，进行大统计量（N→∞）的统计实验方法或计算机随机模拟方法。

理论依据：

大数定理：均匀分布的算术平均收敛于真值中心极限定理：置信水平下的统计误差两个例子：

Buffen投针实验求π射击问题（打靶游戏）Buffon投针实验(1777年）求π:３．各向同性随机投针，则夹角α在[０,π]均匀分布，所以：４．设投针N次，相交次数为M，则相交概率的期望值：N→∞大数定理§1蒙特卡罗方法概述---基本思想２．针与平行线垂直方向夹角为α，则相交概率为：１．平行线间距＝针长＝s一些人进行了实验，其结果列于下表：实验者年份投计次数π的实验值沃尔弗(Wolf)185050003.1596史密思(Smith)185532043.1553福克斯(Fox)189411203.1419拉查里尼(Lazzarini)190134083.1415929§1蒙特卡洛方法的基本思想§1蒙特卡罗方法概述---基本思想1.设r表示射击运动员弹着点到靶心的距离，ｇ(r)表示击中r处相应的得分数（环数），f(r)为该运动员弹着点的分布密度函数，它们反映运动员的射击水平。该运动员的射击成绩为：２.用概率语言来说，<g>是随机变量ｇ(r)的数学期望，即３.假设该运动员进行了Ｎ次射击，每次射击的弹着点依次为r1，r2，…，rN，则Ｎ次得分g(r1)，g(r2)，…，g(rN)的算术平均值代表了该运动员的成绩射击问题（打靶游戏）§1蒙特卡罗方法概述---基本思想４.用Ｎ次试验所得成绩的算术平均值作为数学期望<g>的估计值。例如，设射击运动员的弹着点分布为环数78910概率0.10.10.30.5用计算机作随机试验（射击）的方法为，选取一个随机数ξ，按右边所列方法判断得到成绩。这样，就进行了一次随机试验（射击），得到了一次成绩ｇ(r)，作Ｎ次试验后，得到该运动员射击成绩的近似值

§1蒙特卡罗方法概述---基本思想

收敛性:大数定理作为所求解的近似值。由大数定律可知，如果X1，X2，…，XN独立同分布，且具有有限期望值（E(X)<∞），则

即随机变量X的简单子样的算术平均值，当子样数Ｎ充分大时，以概率1收敛于它的期望值E(X)。由前面介绍可知，蒙特卡罗方法是由随机变量X的简单子样X1，X2，…，XN的算术平均值：§1蒙特卡洛方法概述---大数定理f(x)是X的分布密度函数。则当N充分大时，有如下的近似式蒙特卡罗方法的近似值与真值的误差问题，概率论的中心极限定理给出了答案。该定理指出，如果随机变量序列X1，X2，…，XN独立同分布，且具有有限非零的方差σ2，即§1蒙特卡洛方法概述---中心极限定理其中α称为置信度，1－α称为置信水平。这表明，不等式近似地以概率1－α成立，且误差收敛速度的阶为：O(N-1/2)上式中λα与置信度α是一一对应的，根据问题的要求确定出置信水平后，查标准正态分布表，就可以确定出λα。通常，蒙特卡罗方法的误差ε定义为给出几个常用的α与λα

的数值：α0.50.050.003λα

0.67451.963两点说明：(1)MC方法的误差为概率误差，这与其他数值计算方法是有区别的。(2)误差中的均方差σ是未知的，必须使用其估计值来代替，在计算所求量的同时，可计算出估计值。§1蒙特卡洛方法概述---中心极限定理（2）减小估计的均方差σ，比如降低一半，那误差就减小一半，这相当于N增大四倍的效果。

减小方差的各种技巧:（1）增大试验次数N。在σ固定的情况下，要把精度提高一个数量级，试验次数N需增加两个数量级。因此，单纯增大N不是一个有效的办法。显然当给定置信度α（λα）后，误差ε由σ和N决定。要减小ε：§1蒙特卡洛方法该概述---减小误差

效率:一般来说，降低方差的技巧，往往会使观察一个子样的时间增加。在固定时间内，使观察的样本数减少。所以，一种方法的优劣，需要由方差和观察一个子样的费用（使用计算机的时间）两者来衡量。这就是MC方法中效率的概念。它定义为σ2c，其中c是观察一个子样的平均费用。显然σ2c越小，方法越有效。（1）能够比较逼真地描述具有随机性质的事物的特点及物理实验过程§1蒙特卡罗方法概述---MC优点从这个意义上讲，蒙特卡罗方法可以部分代替物理实验，甚至可以得到物理实验难以得到的结果。用蒙特卡罗方法解决实际问题，可以直接从实际问题本身出发，而不从方程或数学表达式出发。它具有直观、形象的特点。（2）受几何条件限制小计算s维空间中的任一区域Ds上的积分：无论区域Ds的形状多么特殊，只要能给出描述Ds的几何特征的条件，就可以从Ds中均匀产生N个点：因此，在具有随机性质的问题中，如考虑的系统形状很复杂，难以用一般数值方法求解，而使用蒙特卡罗方法，不会有原则上的困难。其中Ds为区域Ds的体积。这是数值方法难以作到的。得到积分的近似值：§1蒙特卡罗方法概述---MC优点（3）收敛速度与问题的维数无关由误差定义可知，在给定置信水平情况下，MC方法的误差为O(N-1/2)

，与问题本身的维数无关。维数的变化，只引起抽样时间及估计量计算时间的变化，不影响误差。这一特点，决定了蒙特卡罗方法对多维问题的适应性。而一般数值方法，比如计算定积分时，计算时间随维数的幂次方而增加，而且由于分点数与维数的幂次方成正比，需占用相当数量的计算机内存，这些都是一般数值方法计算高维积分时难以克服的问题。（4）具有同时计算多个方案与多个未知量的能力（2）使用蒙特卡罗方法还可以同时得到若干个所求量。（1）对于那些需要计算多个方案的问题，使用蒙特卡罗方法有时不需要像常规方法那样逐个计算，而可以同时计算所有的方案，其全部计算量几乎与计算一个方案的计算量相当。例如对于屏蔽层为均匀介质的几何平板，要计算若干种厚度的穿透概率时，只需计算最厚的一种情况，其他厚度的穿透概率在计算最厚一种情况时稍加处理便可同时得到。例如在模拟粒子过程中，可以同时得到不同区域的通量、能谱、角分布等，而不像常规方法那样，需要逐一计算所求量。§1蒙特卡罗方法概述---MC优点（5）误差容易确定根据蒙特卡罗方法的误差公式，可以在计算所求量的同时计算出误差（6）程序结构简单，易于实现在计算机上进行蒙特卡罗方法计算时，程序结构简单，分块性强，易于实现。（1）收敛速度慢（2）误差具有概率性蒙特卡罗方法的收敛速度为O(N-1/2)

，一般不容易得到精确度较高的近似结果。对于维数少（三维以下）的问题，不如其他方法好。由于蒙特卡罗方法的误差是在一定置信水平下估计的，所以它的误差具有概率性，而不是一般意义下的误差。§1蒙特卡罗方法概述---MC缺点（3）计算结果与系统大小有关例如在粒子输运问题中:经验表明，只有当系统的大小与粒子的平均自由程可以相比较时（一般在十个平均自由程左右），蒙特卡罗方法计算的结果较为满意。但对于大系统或小概率事件的计算问题，计算结果往往比真值偏低。在使用蒙特卡罗方法时，可以考虑把蒙特卡罗方法与解析（或数值）方法相结合，取长补短，既能解决解析（或数值）方法难以解决的问题，也可以解决单纯使用蒙特卡罗方法难以解决的问题。真随机数：不可重复，物理方法产生。随机数是实现由已知分布抽样的基本量，在由已知分布的抽样过程中，将随机数作为已知量，用适当的数学方法可以由它产生具有任意已知分布的简单子样。由具有已知分布的总体中抽取简单子样，在蒙特卡罗方法中占有非常重要的地位。总体和子样的关系，属于一般和个别的关系，或者说属于共性和个性的关系。§2随机数的产生和检验蒙特卡罗模拟就是边产生随机数边进行随机模拟的方法准随机数：不具随机性质，只要处理问题能得到正确结果。如放射性衰变，电子设备的热噪音，宇宙射线的触发时间等。伪随机数：可重复，数学方法产生，必须通过统计检验区分：数列的随机性和随机数的分布一个完美的随机数序列可能具有某种分布（如均匀分布，高斯分布等），但具有某种分布的数列却可能完全不是随机的。良好统计分布，容易实现，效率高，周期长，可移植性好等§2随机数的产生和检验分布函数为：最简单、最基本，也是最重要的随机数是在单一区间[0,1]上的均匀分布的随机数，其分布密度函数为由于随机数在蒙特卡罗方法中占有极其重要的位置，我们用专门的符号ξ表示。用ξ1，ξ2

，…代表相互独立且具有相同单位均匀分布的随机数序列。独立性、均匀性是随机数必备的两个特点。如：掷筛子游戏，投掷硬币

用来作为随机数发生器的物理源主要有两种：一种是根据放射性物质的放射性，另一种是利用计算机的固有噪声。用物理方法产生的随机数序列无法重复实现，不能进行程序复算，给验证结果带来很大困难。而且，需要增加随机数发生器和电路联系等附加设备，费用昂贵。因此，该方法也不适合在计算机上使用。§2伪随机数的产生和检验---物理方法在计算机上用物理方法产生随机数的基本原理是：利用某些物理现象，在计算机上增加某些特殊设备，可以在计算机上直接产生随机数。这些特殊设备称为随机数发生器。（1）冯·诺伊曼平方取中法递推公式：以十进制数为例，平方取中法是把一个2S位的十进制自平方后，去头截尾只保留中间2S个数字，然后用102S来除，这样就可以得到在[0,1]上均匀分布的伪随机数序列。例如，设十进制数的2S=4，并取x1=6406，则有：相应的伪随机数序列是0.6406，0.3680，0.1354，0.8333，0.4388等§2伪随机数的产生和检验---数学方法具有周期性，有些数甚至会紧接着重复出现，很少使用。§2伪随机数的产生和检验---数学方法由Lehmer在1951年提出来的，它的一般形式是：对于任一初始值x0，伪随机数序列由下面递推公式确定：（2）Lehmer线性同余法例如乘同余法x0称为种子，改变它的值就得到基本序列的不同区段随机数。a---乘子，c---增量，m---模乘同余法具有在计算机上容易实现、快速等特点，所以乘同余法已被广泛采用。伪随机数的均匀性伪随机数的独立性均匀性是指在[0,1]区间内等长度子区间中随机数的数量是一样的。按先后顺序出现的随机数中，每个随机数的取值与其相距一定间隔的随机数取值之间无关。§2伪随机数序列的统计检验判断伪随机数序列是否满足均匀和相互独立的要求，要靠统计检验的方法实现。对于伪随机数的统计检验，一般包括两大类：均匀性检验和独立性检验伪随机数的均匀性将区间[0,1]分为K个子区间，统计随机数落在第k

个子区间的实际频数nk，它应当趋近于理论频数mk计算统计量如果χ2

值很大，表示远远偏离理想值，因此要求χ2值尽可能小，但如果它趋于0则有可能N已进入循环。通常求和中的每一项的大小约为1，因此χ2的值约为K

。§2随机数的产生和检验---统计检验K的取值不能太大也不能太小，太大反映不出“小区间”的均匀性，太小反映不出“大区间”的均匀性。限制条件（1）概率论中的Pearson定理说明，上式的极限概率分布是χ2分布给出了χ2≤x的概率。整数ν是系统自由度表示独立测量的次数，由于存在一个限制条件，故ν=K-1给出了χ2>x的概率余函数§2随机数的产生和检验---统计检验因此，当给定显著水平α(或置信度1−α)后，由方程Q(χ2|υ)=α或P(χ2|υ）=1−α解出χα值，或从已有的χ2表中查得χα值，如果由（1）式计算出来的χ小于χα，则认为在此置信度下，伪随机数序列在[0,1]中是均匀分布的。（1）顺序相关法它用相邻两个随机数的自相关函数（或相关系数）来标识伪随机数序列的独立性情况，间距为l的自相关函数是相关系数越小，独立性越好。§2随机数产生及检验---独立性统计检验（2）多维频率检验（1）将伪随机数序列用任意一种办法进行组合，每S个随机数作为S维空间中的一个点的坐标值，于是可以构成一个点序列。（2）把S维空间中的单位方体分成为K个子方体，方体边长（3）统计落在第k个子方体中的实际频数nk，它应当趋近于理论频数：例如将2N

个随机数序列分为两组：{ξ1

,

ξ

3,…}和{ξ

2,ξ

4,…}，分别作为平面中N

个点的x和y坐标值。在xy

平面中作K0×K0个小正方形网格区域，落在第(i,j)个网格区域中的实际频数为nij

，则§2随机数产生及检验---独立性统计检验用连续的两个随机数作为点(x,y)的坐标作图可以直观看出随机数之间的关联性。显然左边是不好的随机数。左上显示出条带结构，左下则是规则网格结构。我们也可以把伪随机数分为三列、四列等，用相似的方法进行多维独立性检验§3随机变量的抽样---离散型随机变量X：｛x1,x2,x3,···,xN｝例如：x可取3个值x1,x2,x3，它们出现的几率分别为2/8,5/8,1/8，则随机数小于2/8时实现x1

，在区间2/8,7/8中时实现x2

，大于7/8时实现x3.概率密度f：｛p1,p2,p3,···,pN｝如果从[0,1]区间中均匀抽样得到的随机数ξ满足下式时则物理量x取值为xn

。实际需要的大多数随机变量并不是[0,1]区间均匀分布的，而是有各种不同形式分布密度函数的随机变量。因此，对不均匀的随机变量抽样的关键问题是如何从均匀分布的伪随机变量样本中，抽取符合所要求的分布密度函数的简单子样。§3随机变量的抽样---离散型3MeV光子入射屏蔽铅板的全吸收反射过程反应类型X：光电效应康普顿散射电子对产生反应截面σ：σ1

σ2

σ3

反应几率f：ε1=σ1/σε2=σ2/σε3=σ3/σ归一化：例如Poisson分布是离散型分布对此分布进行抽样得到第n个事件发生的条件为电子对产生光电效应康普顿效应NNYY设连续型变量x在区间[a,b]中取值，可视为将上述的离散情形取连续极限:显然ξ(a)=0,ξ(b)=1且是单调增要求变量x，可由上式解析反解出x(ξ)的函数表达式，即求反函数。这对一些简单的几率密度函数解析表达式是很容易做到的。如粒子随机运动的自由程分布为指数分布：则求反函数后得§3随机变量的抽样---连续型累积函数注意（1-ξ）和ξ同样服从[0，1]的均匀分布§3随机变量的抽样---变换抽样法一维：变换抽样法的基本思想是将一个比较复杂的分布p(x)的抽样，变换为已知的简单分布g(y)的抽样我们希望找到x↔y

之间的对应关系，使得几率密度守恒：例如：黑体辐射的谱密度按频率ω表示时为当希望将谱密度用波长λ=2πcω表示时，有§3随机变量的抽样---变换抽样法显然，当g(y)取[0,1]均匀分布时，问题即化为：寻找y(x)，使其导数为p(x)，然后在[0,1]区间中对变量y

抽样得到均匀分布的随机数，再由x(y)关系得到对应几率密度函数p(x)的随机抽样x

。二维：有两个变量x和y的联合分布密度函数为p(x,y)，欲变换至变量u和v，它们的联合分布密度函数为g(u,v)取联合分布密度函数g(u,v)为均匀分布：§3随机变量的抽样---变换抽样法则我们的任务就是寻找变换式x=x(u,v),y=y(u,v)，以使p(x,y)=|∂(u,v)/∂(x,y)|

对均匀随机变量(u,v)进行抽样，代入变换式得x

和y的抽样。对于Gauss正态几率分布的抽样通过代换可以只考虑简单形式的分布令极角坐标系下的角度为2πv

，半径为现在我们试图通过一个两维联合分布的抽样获得该一维分布的抽样。u和v都是[0,1]区间中的均匀分布的随机抽样，则变换关系式为§3随机变量的抽样---变换抽样法Jacobi行列式即两维分布正为两个独立分布之积。显然抽样x

或y都满足正态分布。可见，为了得到满足一个复杂分布的随机抽样，这里用了两个满足简单分布的随机数。可得反变换§3随机变量的抽样---变换抽样法§4蒙特卡罗方法在积分计算中的应用一维积分平均值法：作变换：得标准积分：§4-1蒙特卡洛方法在积分计算中的应用直接抽样法：在x的定义域[0,1]上均匀随机取点，该均匀分布的随机变量记为ξ，定义随机变量η为：则有因此，只要抽取足够多的随机点，即当n足够大时，In就是积分I的一个无偏估计值。相应的方差为：可见，当f(x)在其定义域内变化较大时，方差较大。§4-1蒙特卡洛方法在积分计算中的应用重要抽样法：当f(x)在其定义域内有显著的起伏变化时，可采用重要抽样法。可将积分稍作变换：适当选取偏倚分布密度函数，使得f*(x)在定义域内变化比较平坦。然后产生[0,1]区间分布密度函数为g(x)的随机变量ξ’，定义：则有：偏倚分布密度函数§4-1蒙特卡洛方法在积分计算中的应用相应的方差为：蒙特卡洛计算结果的方差为：直接抽样重要抽样§4-1蒙特卡洛方法在积分计算中的应用高维积分平均值法标准形式：在实际物理问题中，被积函数在超立方体区域内可能强烈变化。若在积分区域内均匀抽样，积分贡献可能主要来自少数仅仅只有几个蒙特卡洛投点的小区域，从而导致很大的统计误差。所以采用重要抽样法，使得随机点更多地投在取值大的区间§4-1蒙特卡洛方法在积分计算中的应用选取偏倚分布密度函数，并定义（要求方差较小）有：按照偏倚分布密度函数在区域抽样N个子样则：积分的近似值§4-1蒙特卡洛方法在积分计算中的应用一维积分的掷点法一维标准积分：定义：则有：在单位正方形内投N个点，落在曲线f(x)下的有M个，则由于对y的积分可以解析计算，故此法的误差较平均值法大。§4-2蒙特卡洛方法在量子力学中的应用量子力学回顾薛定谔方程定义：费曼传播子量子力学的基本理论告诉我们：系统的所有信息，包括基态、激发态的能量、波函数等论可由费曼传播子给出。特别是和基态有关的信息，可以很方便的得到。§4-2蒙特卡洛方法在量子力学中的应用例如：基态波函数的模方可表示为费曼传播子可以表示为路径积分的形式，对于简单系统，即有其中，作用量§4-2蒙特卡洛方法在量子力学中的应用路径积分量子蒙特卡洛方法实际数值计算中，费曼传播子表达式为

实际计算中，Ｎ足够大即可。积分常数时§4-2蒙特卡洛方法在量子力学中的应用取：，有时§4-2蒙特卡洛方法在量子力学中的应用由此得到基态波函数：其中插入δ函数，得到§4-2蒙特卡洛方法在量子力学中的应用采用Ｍetropolis方法计算基态波函数：首先，选择任意的、连续Ｎ个时间间隔、且的一条路径，计算相应的能量。然后，再接着选一系列路径，每条路径与前一条路径最多只有一个时刻（如）有不同的空间点。采用Ｍetropolis方法来确定满足上述要求的新路径。其中，将随机定下的坐标改变到的过渡概率为。其中，为两条路径的能量差。对于每条路径，利用前述公式计算被积函数的估计值，并累加到求和之中。最终该求和所得值与抽样路径的总数相除得平均值，就得的数值结果。按上述方法，游走足够多的步数后，我们就得到的值。§4-2蒙特卡洛方法在量子力学中的应用变分量子蒙特卡洛方法对于任意的试探函数Ψ，其能量期望值满足对于的计算，采用重要抽样法。当给定试探函数后，由Ｍetropolis方法产生分布的Ｎ个位形。对于每个位形，计算出相应的局域能量，“局域能量”ε基态能量§4-2蒙特卡洛方法在量子力学中的应用则：变分法步骤：（１）选择一个物理上相对合理的基态试探波函数；（２）利用前述方法计算与之相应的能量期望值；（３）改变试探函数中的变分参数值，使试探函数改变一小量，记改变后的试探函数为，并计算相应的能量期望值；（４）计算能量改变值，若改变量小于０，则接受试探函数，否则拒绝，并回到第三步；（５）反复循环，直至能量期望值不再有明显变化为止。§4-2蒙特卡洛方法在量子力学中的应用若上述循环至第Ｍ次终止，则格林函数量子蒙特卡洛方法一维扩散方程：相应的格林函数：§4-2蒙特卡洛方法在量子力学中的应用

格林函数归一，且与ｘ，ｔ无关。扩散方程马尔科夫过程（格林函数）（单步游走的概率分布）

分布§4-2蒙特卡洛方法在量子力学中的应用Fokker-Planck方程：相应的格林函数：

（Δｔ的一阶近似）

构造马尔科夫链

分布力§4-2蒙特卡洛方法在量子力学中的应用３Ｎ维多体定态薛定谔方程的基态解扩散方程的定态解虚时薛定谔方程（取）：

（具有势函数的扩散方程）扩散项分支项格林函数：

（借用Ｄｉｒａｃ记号）§4-2蒙特卡洛方法在量子力学中的应用上述格林函数正是量子力学中的时间演化算符在坐标表象下的矩阵元，即费曼传播子。易知道当，或者说足够大时，上述右边的求和中只有基态才有贡献，这个算符（时间演化算符）行为就如同作用在基态波函数上。解析计算格林函数，得到§4-2蒙特卡洛方法在量子力学中的应用上面给出的是格林函数短时间近似结果。根据此结果，我们在格林函数蒙特卡洛模拟中，就必须进行大量的短时间间隔的游走，最终使其分布近似满足基态波函数。扩散步游走分支步由游走到的权重需乘因子模拟效率不高。Ｒeynoids方法§4-3蒙特卡洛方法在统计物理中的应用物理量的观测值微观粒子的某物理量在相空间的分布平均值高维积分相空间态矢热力学平衡状态下（恒温Ｔ）：观测量Ｈamilton量分布密度函数配分函数§4-3蒙特卡洛方法在统计物理中的应用上述公式涉及的是高维积分问题。只有理想气体、谐振子系统、Ising模型等极少数类型的问题可以解决求解。大多数情况下，只能借助近似方法求出。蒙特卡洛方法正则系综，

除掉动量以外的其它的相空间坐标当粒子间的相互作用与动量无关时，动量项的贡献可以被积分，这相当于将Hamilton量中的动量项去掉。则，平衡态下的概率分布为Boltzmann分布。§4-3蒙特卡洛方法在统计物理中的应用Boltzmann分布密度函数为：可见，所有对应于大能量的状态对应观测量积分的贡献都很小。随机选择（均匀抽样）ｎ各状态，则有§4-3蒙特卡洛方法在统计物理中的应用由于上式中大部分状态对求和的贡献很小，故抽样的效率比较低。为有效的进行计算，应采用重要抽样法。用Metropolis方法产生Boltzmann分布的ｎ个状态，则湘潭大学物理学专业基础课计算物理及其应用

材料与光电物理学院变分量子蒙特卡洛方法求解简单谐振子的基态能量

1.引言2.变分量子蒙特卡洛方法

3.计算步骤

4.结果及分析

目录1.引言经典力学中，质点在平衡位置附近的最基本最简单的运动是简谐振动。在量子物理中与此对应的微观粒子的运动就是谐振子。简单谐振子的理论在应用上有很大价值，因为经典力学告诉我们只要选择适当的坐标，任意粒子体系的微小振动都可以认为是一些相互独立无关的振子的集合的运动。普朗克在他的辐射理论中将辐射物质的中心当作一些谐振子。从而得到和实验相符合的结果。在分子光谱中，可以把分子的振动近似地当作谐振子的波函数。另外在量子场论中电磁场的问题也能归结成谐振子的形式。因此，谐振子问题的地位在物理学中非常重要。2.变分量子蒙特卡洛方法设量子体系的波函数为，则量子体系的Schrödinger方程为简单谐振子体系的哈密顿量为V(rD)为谐振子势，取其中，D表示的是维数，可取1、2和3维。取简单谐振子含调节参数的试探波函数形式2.变分量子蒙特卡洛方法使用量子单位，可得在计算过程中定义局部能量：波函数下的平均值为2.变分量子蒙特卡洛方法取N个含M个随机状态的系综，则：标准方差：其中3.计算步骤

计算步骤：取试探波函数形式（在实际的模拟过程还取到了试探波函数的形式，用以来对比优劣）；定义调节参数α的变化范围，用循环的语句来控制变化；取N个系综，第个系综随机取M个随机数；N个系综同步随机演化，每个系综变换按照Metropolis方法来产生马尔柯夫链，判定接受的依据为3.计算步骤计算不同参数α时，EL平均值：和标准方差：4.结果及分析

计算结果：当取一维谐振子基态试探波函数为时的结果如下：

图2横坐标为试探波函数的调节参数，纵坐标为变化调节参数后对应能量的标准方差图1横坐标为试探波函数的调节参数，纵坐标为变化调节参数后对应的能量值4.结果及分析当取一维谐振子基态试探波函数为时的结果如下：

图4横坐标为试探波函数的调节参数，纵坐标为变化调节参数后对应能量的标准方差图3横坐标为试探波函数的调节参数，纵坐标为变化调节参数后对应的能量值4.结果及分析结果分析：谐振子系统具有精确的解析波函数，用变分量子蒙特卡洛方法计算出的结果与解析结果很好的吻合。用变分量子蒙特卡洛方法求解量子体系的基态能和基态波函数时，可以加入对局部能量的标准方差的统计，在结果中取标准方差最小的试探波函数形式，再由能量取最小来确定的α值，这样，能够准确的从一系列试探波函数形式中选出最优基态波函数，并同时确定比较精确的基态能。证实了变分量子蒙特卡洛方法方法在解决量子体系问题方面的可行性和可靠性。统计力学：系综时间平均统计力学：各态历经假设系综平均统计力学：平均值随机过程动力学变量转移概率分布函数归一化条件Markov过程主方程Fokker-PlanckEquation布朗运动爱因斯坦模型概率演化的主方程差分形式微分形式平衡态简单抽样方法BoltzmannStatisticalAverage重要抽样方法BoltzmannStatisticalAverage经典粒子系统经典自旋系伊辛模型XY模型海森堡模型经典自旋系量子蒙特卡罗方法湘潭大学物理学专业基础课计算物理及其应用

材料与光电物理学院蒙特卡洛方法的简单应用

——求圆周率和模拟氢原子电子云

1.1蒙特卡洛处理的两类问题1.2求圆周率的数学模型1.3计算结果1.4氢原子电子的分布密度1.5模拟结果目录确定性问题1.1蒙特卡洛处理的两类问题随机性问题原子核物理问题、运筹学中的库存问题、随机服务系统中的排队问题，动物的生态竞争和传染病的蔓延计算多重积分、求逆矩阵、解线性代数方程组、解积分方程、解某些偏微分方程边值问题和计算微分算子的本征值基本思想：针对待求问题，根据物理现象本身的统计规律，或人为构造一合适的依赖随机变量的概率模型，使某些随机变量的统计量为待求问题的解，进行大统计量N→∞的统计实验方法或计算机随机模拟方法。1.2求圆周率的数学模型y2rr2rOrx在正方形区域内产生N个随机点{(x1,y1),(x2,y2),…,(xN,yN)}记录落在圆内的点数为M，如果取半径r=1，那么圆周率为：1.3计算结果π值随模拟次数的变化1.4氢原子电子的分布密度由原子物理理论和量子理论可知，氢原子态的波函数只是半径的函数，与和无关。而氢原子中电子半径的分布密度，即电子在半径处单位厚度球壳内出现的几率习惯上把这种分布形象地称为电子云。氢原子的基态即１s态（n=1,l=0,m=0），有，是D的最大值处的r值，其值与其中玻尔半径相同。r0是D收敛处的r值，即D的收敛点。1.4氢原子电子的分布密度氢原子的2s态(n=2,l=0,m=0),有

氢原子的3s态(n=3,l=0,m=0),有

1.5模拟结果氢原子态电子云模拟图氢原子态电子云模拟图氢原子态电子云模拟图氢原子组合电子云模拟图附录：程序具体FORTRAN程序如下：programmain

implicitnone

real(kind=8),parameter::m=1.0e8

integer::n,counter=0

real(kind=8)::x(1:m),y(1:m),Pi,r

open(unit=10,file='random_num.txt')

callrandom_seed()

CALLrandom_number(x)

callrandom_number(y)

don=1,m

r=sqrt(x(n)**2+y(n)**2)

if(r<=1)then

counter=counter+1

endif

enddo

pi=(counter*1.0)*4.0/(m*1.0)

writ