计算物理学练习题及参考解答

N( , 2 )
1
2
exp[( x ) 2 / 2 2பைடு நூலகம்] 1 2
比如期望值为 1，方差为 1 的高斯分布表达式为 N(1,1)
exp[( x 1) 2 / 2]
16．对物理问题的计算机模拟所需要的伪随机数应当满足什么样的标准？有哪些统计检验方法？ 答：良好的统计分布特性；高效率；循环周期长；产生程序可以移植性好；可以重复产生。 统计检验有：均匀性检验；独立性检验；组合规律检验；无连贯性检验；参数检验等等。 17．在蒙特卡洛方法应用中减小方差的基本技术：重要抽样法，分层抽样法，控制变量法和对偶变量法。 然而，单独使用这四种减小方差的技巧仍然有其局限性。人们发展了一些复合蒙特卡洛计算技术，如 适应性蒙特卡洛方法和多道蒙特卡洛抽样方法等。这些蒙特卡洛技巧对于被积函数在积分范围内具有 多个尖峰的情况，特别具有实用价值。 18．真随机数列 答：真随机数列是不可预测的，因而也是不能重复产生的数据序列。 19．伪随机数列 答：通过某些数学公式计算而产生的随机数列 20．同余产生器及程序代码 答：一般通过如下的线性同余关系式来产生数列
2
x n 1 (ax n c)(mod m)
n xn / m
float Random(int n, int m, float seed, float a, float b) { int i; float r; r = seed; for (i = 1; i <= n; i++) r = (a * r + b) % m; return r / m; } 21．均匀性检验 答：均匀性检验又称频率检验，它用来检验用随机数（样本值）确定的经验频率和均匀分布频率是否有显 著性差异。常用的统计检验方法有 x 2 检验和累积频率检验（K-S 检验） 。 22．随机变量抽样 答：指的就是由给定分布函数产生随机数的方法。首先，在[0,1]区间抽取均匀分布的伪随机数列，再从中 抽取满足给定分布密度函数的简单子样，并且各个伪随机数相互独立。 23．连续分布的随机变量抽样一般有哪些方法？ 答：直接抽样法；变换抽样法；舍选抽样法；复合抽样法；近似抽样方法 24．试述离散型分布的随机变量的直接抽样 答：对于离散序列数 x1 , x2 , xi , 给定每个数的取值概率 p1 , p2 , pi , ，则我们可定义其分布函数 F(x) 如下：
λe -λλ 27．对指数分布 f(x) 0
答：第一步：求分布函数 F(x) 第二步：抽样 ~ [0,1] 第三步：计算分布函数的反函数
x 0,λ 0
其他
直接抽样

x
-
f(x)dx λe λx dx 1 e λx
＝F( ) 1 - e -
g (u )du P[u u u du ] g (u ) 称为 u 的概率密度分布函数，它表示随机变量 u 取 u 到 u du 之间值的概率。 而 G (u ) g ( x) dx
u
则称为 u 的分布函数。 15．高斯分布可以由给定的期望值 和方差 2 完全确定下来，通常用 N( , 2 ) 来表示
3．计算物理学是物理学、数学、计算机科学三者结合的产物，它也是物理学的一个分支，与理论物理、 实验物理有着密切的联系。 4．计算机在物理学中有哪些应用？ 答：计算机数值分析、计算机符号处理、计算机模拟、计算机实时控制 5．计算机技术有各种各样的算法，可以概括为最基本的两类：串行计算和并行计算。 6．理论物理在实际计算中遇到许多困难：非线性问题求解和非对称问题的求解；自变量较多问题求解； 非规则界面问题求解等。 7．计算物理的优点有：省时省钱；具有更大的自由度和灵活性；能够模拟极端条件下的实验。 8．第一原理方法是基于量子力学基本原理建立起来的；分子动力学方法是基于经典力学基本原理建立起 来的；蒙特卡罗方法是基于统计力学基本原理建立来的。 9．计算机模拟一般有哪两种类型？ 答：随机模拟和确定性模拟，比如蒙特卡罗模拟和分子动力学模拟。 10．什么是蒙特卡罗模拟？它的应用一般有哪三种形式？ 答：通过不断产生随机数序列来模拟过程。直接蒙特卡罗模拟、蒙特卡罗积分、Metropolis 蒙特卡罗模拟。 11．蒙特卡洛方法的理论依据 答： （1）大数法则：人们发现，在一个随机事件中，随着试验次数的增加，事件发生的频率趋于一个稳 定值；人们同时也发现，在对物理量的测量实践中，测定值的算术平均也具有稳定性。大数法则反映了大 量随机数之和的性质。 （2）中心极限定理：中心极限定理，是概率论中讨论随机变量和的分布以正态分布为极限的一组定 理。这组定理是数理统计学和误差分析的理论基础，指出了大量随机变量近似服从正态分布的条件。中心 极限定理告诉我们：在有足够大，但又有限的抽样数 n 的情况下，蒙特卡洛估计值是如何分布的。 12．请简述著名的巴夫昂（Buffon）投针实验。并写出用 Matlab 实现的代码。 （给出方程、算法框图、程 序） 答：在平滑桌面上画一组相距为 s 的平行线，向此桌面随意地投掷长度 l s 的细针，那么从针与平行线相
如图第一项限中单位正方形内投点在圆内的概率即为单位圆面积的四分之一。
2 数学方程： 4 dx1 dx2 (1 x12 x2 )

1
0

1
0
算法框图： 产生随机点 （ξ， η） M 个； 统计其中满足条件 2 2 1 的点的个数 N； 计算π值 4 N / M 。 Matlab 程序：P=4/100000*length(find(sum(rand(2,100000).^2)<1))
y h ( x) ， x g ( y )
选择 ( y) ，使得
f ( x) ( y)
dy ( h( x ))h( x) dx
则可对 ( y) 抽样得到δ，通过变换 g ( ) ，得到满足分布密度函数 f ( x ) 的随机抽样。
3
更为一般的情况，设连续型随机抽样 ( , ) 的分布密度函数为 f ( x , y ) 。考虑变换
1 ln(1 )
28．梅氏游走法计算氦原子中两电子间库伦作用的平均值。（给出动力学方程、数值解方程、算法框图、 程序） 答：! 梅氏游走法计算氦原子中两电子间库伦作用的平均值 program main dimension x(6) !输入抽样参数 Nt,Ng,Nf,Ns,dx Nt=1000 !热化步数
4
10
20
30
40
Ng=100 !样本组数 Nf=10 !抽样间隔 Ns=10000 !每组抽样个数 dx=0.1 !游走步长 !初始化：随机设置初始值 x do 10 i=1,6 call random(RND) x(i)=0.01*(rnd-0.5) !两电子在核附近 !热化：消除初始化影响，趋于平衡分布 do 20 j=1,Nt call walk(RND,dx,x) !分组间隔抽样，计算平均值和误差 su=0 su2=0 sdu=0 do 40 ig=1,Ng !样本分成 Ng 个组 ug=0 ug2=0 do 30 k=1,Ns !Ns 间隔 Nf 抽样，Ns 个样本为一组 call walk(RND,dx,x) if(mod(k,Nf).ne.0) goto 30 x12=x(1)-x(4) y12=x(2)-x(5) z12=x(3)-x(6) r12=dist(x12,y12,z12) u=1/r12 ug=ug+u ug2=ug2+u*u !组内求和 continue ug=ug/Ns sigmag2=ug2/Ns-ug*ug dug=sqrt(sigmag2/Ns) !组内平均、方差和误差 su=su+ug su2=su2+ug*ug sdu=sdu+dug !组间求和 avu=su/Ng sigma2=su2/Ng-avu*avu du1=sqrt(sigma2/Ng) du2=sdu/Ng del=du1-du2 !组间平均、方差和误差
1
交的概率就可以得到 的数值。 p
2N 2 M ， M p N
clear S=1; %平行线间距 L=1; %针长 N=1000000; %总投针次数 M=0; for i=1:N x=rand*S/2; %针到最近平行线的距离 a=rand*pi/2; %偏角 if(x<sin(a)*L/2) M=M+1; %统计相交次数 end end testpi=2*N/M %pi 的实验值 13． 在考虑蒙特卡罗模拟的精确度时， 不能只是简单地减少方差和增加模拟次数， 还要同时兼顾计算费用， 即机时耗费。 14．假定我们研究连续的随机变量，由随机变量的分布可以得到它取某给定值的概率，即
F ( x ) pi 。
xi x
在区间[0,1]上取均匀分布的随机数ξ，判断满足下式的 j 值：
F ( x j 1 ) F ( x j )
则抽样值η为 x j ，η分布符合分布函数 F(x)的要求为。 25、试述连续分布的随机变量的变换抽样法。 答：设连续型随机变量η的分布密度函数为 f ( x ) 。要对满足分布密度函数 f(x)的随机变量η 抽样较难时 可考虑通过其它已知函数的抽样来得到。考虑变换
!输出 avu,du1,du2,del 100 open(12,file='out.dat') write(12,1000) Nt,Ng,Nf,Ns,dx,avu,du1,du2,del close(12)
5
1000 format(4i10,5f15.4) end 计算距离的函数子程序 function dist(x,y,z) dist=sqrt(x*x+y*y+z*z) return end ! 计算权重的函数子程序 subroutine weight(x,f) dimension x(6) r1=dist(x(1),x(2),x(3)) r2=dist(x(4),x(5),x(6)) f=exp(-3.375*(r1+r2)) return end ! 梅氏游动一步的子程序 subroutine walk(RND,dx,x) dimension x(6),x0(6) call weight(x,f0) do 10 i=1,6 x0(i)=x(i) call random(RND) ! 存旧 10 x(i)=x(i)+dx*(RND-0.5) ! 生新 call weight(x,f) call random(RND) if(f.ge.f0*RND) goto 30 !游动 do 20 i=1,6 20 x(i)=x0(i) !不动 30 return End 29．有限差分法 答：微分方程和积分微分方程数值解的方法。基本思想是把连续的定解区域用有限个离散点构成的网格来 代替，这些离散点称作网格的节点；把连续定解区域上的连续变量的函数用在网格上定义的离散变量函数 来近似；把原方程和定解条件中的微商用差商来近似，积分用积分和来近似，于是原微分方程和定解条件 就近似地代之以代数方程组，即有限差分方程组 ，解此方程组就可以得到原问题在离散点上的近似解。 然后再利用插值方法便可以从离散解得到定解问题在整个区域上的近似解。 30．采用有限差分法求解微分方程时可以用直接法、随机游走法和迭代求解法。其中迭代法被广泛采用， 有直接迭代法、高斯－赛德尔迭代法和超松弛迭代法。 !
u h1 ( x, y ) ， v h2 ( x, y ) ，
f ( x, y ) g (u , v) J g ( h1 ( x, y ), h2 ( x, y )) J ， J
u / x u / y ， v / x v / y
这样就可以通过抽取满足分布密度函数 g (u , v ) 随机抽样 ( , ) 得到待求的满足分布密度函数 f ( x , y ) 的 随机抽样 ( , ) ，其中 h1 ( , ) ， h2 ( , ) 。 26．试给出一个用随机数计算π的 Matlab 程序。（10 分） 解：物理模型：
计算物理学练习题及参考解答
1．计算物理学的英文表示：computatioal physics 或者 computer physics 2．什么是计算物理学？它与理论物理、实验物理有什么区别和联系？ 答：计算物理是指以计算机及计算机技术为工具和手段，运用计算数学的方法解决复杂物理问题的一门应 用科学。 计算物理方法是除理论方法和实验方法之外的第三种研究手段，计算物理现已成为物理学研究的三大支柱 之一，它与实验物理和理论物理的关系如下图：