数值计算方法讲座2

《数值计算方法讲座2》由会员分享，可在线阅读，更多相关《数值计算方法讲座2（17页珍藏版）》请在装配图网上搜索。

1、kkfk= f (x, f (x)axbx第二章 数值积分和Monte Carlo方法第一节数值积分S = Ib f G )dxa令 h = x - x , x = a,x = b , 则k+1 k 0nS = IB J f (x)dxk =0 %f (x)= f +(x- x )f +k=f( xk)零阶近似 f 1 )= f +0(h) kS = hB f +0(h)k=0=hB+0(h)k=0一阶近似f )= fk +1 - x)fk+1 - fk +0(h 2 ) hJ xk+i (x - x )dx = x+1 - M - x (x- x )=J_ h2xk22 k k+1k 2kS

2、 = Eff h +1 (f - f )h + oC2)，k 2 k+1 kk=0 v/-h 习 2 (f +f )+0 C 2)k=01 /、/ 从直观看，用2(fk + fk+近似f (，)比只用f或、好。这方法也称Trapezoid 方法。这样的数值积分方法的优点： 简单直观，误差可以控制缺点： “平均主义”，在I f G)B0的区域，f (x)Ax对S贡献很小，但消耗 1k同等的机时。在多自由度系统这弱点尤为特出。问题：直观地看，零级近似和一级近似的差别在哪？习题：编程序数值计算高斯积分。第二节 Monte Carlo方法如何用随机方法求积分？例如，可用抛石子方法。但这方法不比简单的数

3、值积分有效。1 .简单抽样的Monte Carlo方法均匀地随机地选取a M中M个点x ，显然，kS - 1 尤f (x ) + 0I/JM)M kk=1当M足够大，当然可以得到足够好的积分值。问题为什么误差是O（/vM）?o(/,M)比较误差:Monte Carlo 方法数值积分Trapezoid方法AS x 1/ M 2(h = 1/ M)对单自由度而言，数值积分方法要有效得多。对多自由度，例如d自由度，Monte Carlo 方法AS x 1/M！数值积分Trapezoid方法AS x 1/ M 2/ d(hd = 1/ M)当d非常大，数值积分方法根本没法和Monte Carlo方法比

4、较。我们当然可以再改进数值积分方法的精度，但这种改进的量级没 法和d的大小比拟。在多体系统的数值模拟中，d通常至少是10 4！Monte Carlo方法真的是完美了吗?当然不是。平均主义的弱点其实还没改进卜面我们引入所谓的重要抽样Monte Carlo方法当引入重要抽样方法后，每次抽样的样品可能不独立如何取得独立的抽样，是Monte Carlo方法的重点所在!2.重要抽样的Monte Carlo方法如果被积函数KM是均匀的函数，则简单抽样方法已经可以得到 相当准确的积分值。如果被积函数KM不是均匀的函数-这在高维积分十分常见，则 必须引入重要抽样。我们希望在对积分贡献大的区域多取样品，在贡 献

5、小的区域少取样品。设积分为S =Ib f G)W (尤)dxa其中fx)是均匀函数，W(x)是不均匀函数，而且W(x) 0，可以归一 化，给予概率分布的含义假设我们可以按照分布W(x)得到M个点XJ，则S = M 义 f (xz)+O(/、M)l=1注意：现在Wx)不出现在求和式子中，而是体现在X1的分布里。证明：如第一节方法，把a ,b分为n个小区间，x,落在，x +力区间的数目约为W (x ) h - Mk对 x ex , x + h f (x ). f (x )lk1 Mf (x ). 1 云 f (x )W (x ) hMMi Mk kk=0显然limM siLf (x )= lim

6、f (x)W (x )h = * f (x机=SM i /， kkk关键：如何按分布w (x),产生M个七(即上面的再)问题：可以用产生随机数的方法产生W (x)吗？答，：多自由度，有相互作用时不可能3. Markov 过程产生x 的方法：k构造一个Markov过程，即给出一个动力学规则，由随机地产生+1。那么，给一个x0，可产生x , x ,x ,x 0 1M 0M 0 + M如果随机过程满足一定条件，则当M0足够大时，即达到“平衡态”时，x ,x按W (x)分布。M0 +1M0 +M两个条件|：各态历经从概率上说，在有限时间内x.可走遍a , b :细致平衡Markov过程由从x.到x.+

7、1的转移概率定义。设T (x T x，)为过程从x转移(或跃迁)到x，的概率，T (x，T x)为从x，到x的概率。则细致平衡条件为T (x xf) _ W (xf)T (xf Tx) = Wtx)记忆：从转移到r的概率正比与w (x)证明：设W (x, t)为x的“非平衡态”分布W(x 其)=T (x xf)w (x, t)+ T (x x)W (x , t) dtxxt 8W (x , t) W (x, 8)dW (x, t)0 dt T (x xW (x, 8)= T (xf x )W (x, 8 )x，x，细致平衡条件显然保证W (x )满足上述方程，所以，W (x, 8) = W (

8、x)不过，细致平衡条件是充分条件。第三节 Metropolis算法和Heatbath算法Markov 过程的全部信息包含在转移概率T (x x，)中细致平衡条件是平衡态W (x )的要求是否各态历经常常可以直观判断T (x xf)必须给出从任意k到任意*的概率，但这并不一定要求从x到任意x的概率都是零。各态历经只要求在有限时间内，即T (x x)的有限次作用，能到达 x。1. Metropolis 算法Metropolis (1915-1999)The Paper was cited 7500 times from 1988 to 2003令T (x x，)_W (x )/W (x)W (/)

9、v W (x )一 1W (x 贝 W (x)=Min (1, W(x )/W(x)设 W (xf) W (x)T (x T x，)_1_W (x * )T (x，T x) W(x)/W(x，) W(x)即T(x T x)满足细致平衡条件。但是，注意到T (x T xf)是转移概率，所以应有归一化条件z T (x T xf)= 1 x上面的T (x T x，)还不满足归一化条件。 m所以，完整的Metropolis算法的转移概率是T (x T x，)= S(x T x，) T (x T x，)其中S(x T x)是从x选中x，的概率，而且S(x T x) = S(x T x)。这样的转移概率仍

10、然满足细致平衡条件。当然，归一化条件也能保证。例如，可选1.S(x T x)=常数，即均匀选取x，S (x T x，) = J常数x e x-s,x + s 0x史x-s,x + s归纳起来，Metropolis算法包括如下步骤：1 ) 设x广x，按S(x T x)选取尝试x，2)以概率T (x x)取x+广x以概率 1 -T (x X)取= x (!)第二步要记住，初学者容易误解，以为一旦x不被采纳，重新选取 尝试x。不难理解，无论是1或2方案的S(xtx)，Metropolis算法都满 足各态历经条件，因为只要即皿)，0，在有限时间取到x的概率不 为零。Metropolis算法非常普遍。一

11、个重要的例子，在物理学中，常取W (x)为正则分布，W ( x) x e -H ( x )/ kT其中H(x)代表能量，T是温度，k是Boltzmann常数。所以，跃迁概率为T (x X，)_e-(H (x，)-H(x) / kT H (x，) H (x )m X |1H (xr) r，取 X= x，否则 X= xe -(H (x，)-H (x)/ kT 如何取M0和M ？这与具体系统有关。通常通过尝试改变M以判断M是 o0否够大，即随机过程是否已到达“平衡”态。一般 M 10M0。不过，对单自由度问题，M0不重要。练习：H (x)= Lb x22计算，和解析结果比较H (x)=1 b x2

12、+J人 x4， 计算 24讨论的作用 对没有解析解的系统，如何判断数值结果的可靠性和正 确性?有兴趣的同学还可以和数值积分比较，会发现MonteCarlo方法在单自由度情形没有优势2. Heat-bath 算法Heat-bath算法没有Metropolis算法那么普遍但也是多体系统研究中的典型方法之一Heat-bath算法直接取转移概率为T ( 19= W x)这算法自然是各态历经并且满足细致平衡条件的。 Heat-bath算法的特点是T、x T X)与x无关在某些情形会显示出优越性当W (x)较复杂，在计算机上不一定能够找到简单的实现事实上，对单自由度情形，已经是直接产生平衡态分布匝(x)

13、例如，如果x只取x1和x2,可取W (x1)x = x1W (x 2)x = x 2T (x T x，) =计算机上的实现：取随机数r e 0,1 如果 W( x1 r 取 x, = x1，否则 x, = x2 这一算法可简单推广到x取多个分立值情形。当x是连续变量时，不一定能实现Heat-bath算法。但我们可以结合 Metropolis算法和Heat-bath算法。例如，取T (x T x，)= S(x T x ) T (x T x，)W(x)/(W(x) + W(x)x = xT (x T x=W (x )/(W (x) + W (x )x = xh 1 1 10otherwise这算法

14、Metropolis的特征多些。归纳起来，Metrpolis算法或Heat-bath算法可以相当普遍地实现 重要抽样的Monte Carlo模拟方法。但是，对多体系统，抽样得到的样品可能不独立，即有关联。设关联时间为T ,即每T个样品才有一个是独立的。那么，误差为S x O(OM)当T很大时，Monte Carlo方法遇到极大困难。对物理学家而言， 这是Monte Carlo算法需要解决的重要问题。第四节 Langevin 方程问题：系统的Hamiltonian量为S，平衡态分布为exp(-S)，(这里温度已吸收到S)。假设系统t = 0时处于一初始状态，系统如何演化至平衡态?单自由度S =

15、S G)尤：实数Langevin 方程dx (t)dS (x)f )=+ 门 ) dtdxn:高斯随机数(n (t)； = 0；n (t)n (t); = 2 5 (tt，)对固定tp () e bn 2n (t) = j+a)dn *n - e bn2n Z 8Z = J dn e-bn2如果没有随机力，平衡态为国也=0，即能量取极小值。 dtdx如果存在随机力，体系会被推离能量极小，处于某种能量较高的平衡态。由于随机力的存在，Langevin方程存在复杂性，因为我们必须考 虑对随机力带来的奇异性。为了简单起见，我们对 时间分立化在数值模拟中应用较直观AtLangevin 方程A (t)=

16、x (t + At)-x (t)= At+At Mt) dx:.nAx (t )=A?nn. 2(t)n( t -)门=5 td S (x (t)a()一Ft At+2Atn U方程的解x*是随机变量，在数值模拟中给定初始值，xn(t)还 不确定，与随机力有关。也就是说，在t时刻，x遵从一个分布物理量。G)的平均值 =j dx P (x ; t )。( x ) 问题： + d x (t)n 2 d 2 x (t)nnii(t )r dS (x 0) _z(r- -r *+E。),、/d2o(x (t)r dS(x (t)_+圣*+Atn (t),.xn Q与n Q无关，只与n (t -At以及

17、更早的随机力有关a3)Ax(t)二11t如3)a s 3) axX)-nAta x vt)又=j dx P G; t)-)ao as a20、k a x a x a x 2 /=j dx O (x)(a2a a s)ax2 ax ax/个P (x ; t)分步积分还作用于P （x ； ）这里做分步积分时，假设P （8；t） = 0另一方面A 叫，* =j dx oQapW Ata tFokker-Planck 方程a p (x ； t) =H atH 5FPP (x ; t)a s Q)-+ax axax J当,竺公t 0a tHpP (x ; 8)= 0FP(a显然P (x ; 3)s (x )