《计算物理学》课程教学资源（习题）第二至八章 Computational Physics

第二章习题 (1)采用线性同余法(参见公式(2.2.3)产生伪随机数。取a=5,C=1,m=16和 x=1.记录下产生出的前20数,它产生数列的周期是多少? (2)取a=137,c=187,m=256和x0=1,用线性同余法产生出三维数组 (5n,5n,5+2}和二维数组{5n,5n},然后分别绘出其三维和二维分布图形 (3)用“投针法”计算出圆周率的数值,画出程序流程框图,并编写程序 (4)已知电子在物质中的作用截面总=0mpm+0光电+O电子对,试写出电子在物质 层中相互作用的抽样程序框图和程序 (5)编写一个程序按照门=-h5产生随机数序列{},并绘图表明其分布满足分布 密度函数 f(x)=2,x>0,2>0 0.其它 (6)了轻子的平均寿命为34×10-13s,试写出N个了轻子在实验室系中以速度v运动 的飞行距离的抽样程序框图和程序 (7)写出各向同性分布的角度θ,q抽样程序(d2= sin e de) (8)如分布密度函数为f(x,y)=c,(其中,x21,y≥0.n为整数,试写出抽样 程序框图和程序 (9)证明 Breit- Wigner分布 f(x) 可以通过x1=x0-TcO(z5)抽样得到 (10)归一化黑体辐射频谱为 15 f(x)dx re-1a(其中hv 证明如下抽样步骤得到的抽样分布满足上面的分布,求出它的抽样效率。 抽样步骤:让L等于满足下面不等式的整数l的最小值 ∑25x2/0 然后置r≈、1 l(253545),其中5为[0,1区间均匀分布的伪随机数

第二章 习题 （1） 采用线性同余法（参见公式(2.2.3)）产生伪随机数。取 a = 5 ，c =1，m =16 和 0 x =1.记录下产生出的前 20 数，它产生数列的周期是多少？ （2） 取 a =137 ， c =187 ， m = 256 和 0 x =1 ， 用 线 性同 余法 产 生出 三维 数 组    n n n , , + + 1 2 和二维数组   n n , +1 ，然后分别绘出其三维和二维分布图形。 （3） 用“投针法”计算出圆周率的数值，画出程序流程框图，并编写程序。 （4） 已知电子在物质中的作用截面  总 =  compton + 光电 + 电子对 ，试写出电子在物质 层中相互作用的抽样程序框图和程序。 （5） 编写一个程序按照   ln  −1 = − 产生随机数序列 i ，并绘图表明其分布满足分布 密度函数      = − 0, 其它 , 0, 0 ( )    e x f x x 。 （6）  轻子的平均寿命为 s 13 3.4 10−  ，试写出 N 个  轻子在实验室系中以速度 v 运动 的飞行距离的抽样程序框图和程序。 （7） 写出各向同性分布的角度  , 抽样程序（ d = sin  d d ）。 （8） 如分布密度函数为 n x y x n e f x y − ( , ) = ，（其中， x 1, y  0, n 为整数），试写出抽样 程序框图和程序。 （9） 证明 Breit-Wigner 分布 2 2 0 ( ) 1 ( ) − +   = x x f x  可以通过 0 ( ) i i x x cot = −    抽样得到 。 （10） 归一化黑体辐射频谱为 ) ( ) 1 ( 15 ( ) 4 4 k T h dx x e x f x dx x   = − = 其中 证明如下抽样步骤得到的抽样分布满足上面的分布，求出它的抽样效率。 抽样步骤：让 L 等于满足下面不等式的整数 l 的最小值， 4 4 1 1 1 90 l j j   =   然后置 ln( ) 1  2  3  4  5 L x = − ，其中  i 为 [0，1] 区间均匀分布的伪随机数

(11)对正则高斯分布抽样 p(x)dk=1x7)、(x-) (12) Gamma函数的一般形式为 f(x)h=a,x"eh,x≥0 (n 抽样证明其抽样方法可以为 7=--l(5152….5n15n) (13)x2分布的一般形式为 f(x)dx c edx x>0 22T(n/2) 抽样证明其抽样方法可以为 7=∑x2,其中x1,x2,…xn为标准正态分布的n个独立抽样值 (14)选择偏倚分布密度函数g(x)=eˉx,用蒙特卡洛重要抽样法求积分 d x (15)编写一个程序,采用 Metropolis随机游走的方法产生按高斯分布 f(x)=Aep[-x2/(2a2)(o2=) 的随机点。抽样中常数A的值需要知道吗?试决定接受点与试探步数之比,到达平 衡分布的时间与最大试探步长O的关系。(提示:判断“平衡”的标准是 ≈a2)。δ选多大较合理? (16)用 Metropolis随机游走的方法计算积分 xe,0o≤x≤4 (17) Laplace方程及其边界条件为 0 o(x,0)=q(x)=0,o(0,y)=(y)=1 用随机游走的蒙特卡洛方法数值求解正方形场域(0≤x≤1,1≤y≤1)的势函数。 第三章习题 (1)利用蒙特卡洛方法计算三维、四维、五维和六维空间的单位半径球的体积

（11） 对正则高斯分布抽样: ( ) dx x p x dx       − = − 2 2 2 exp 2 1 ( )     . （12） Gamma 函数的一般形式为 , 0 ( 1)! ( ) 1  − = − − x e dx x n a f x dx n ax n 抽样证明其抽样方法可以为 ln( ..... ) 1 1 2 n 1 n a  = −    −  . （13） 2  分布的一般形式为 , 0 2 ( / 2) 1 ( ) / 2 1 / 2 / 2   = − − x e dx x n f x dx n x n 抽样证明其抽样方法可以为 = = n i i x 1 2  , 其中 n x , x ,....x 1 2 为标准正态分布的 n 个独立抽样值. （14） 选择偏倚分布密度函数 x g x e − ( ) = ，用蒙特卡洛重要抽样法求积分   − 0 3/ 2 x e dx x . （15） 编写一个程序，采用 Metropolis 随机游走的方法产生按高斯分布  ( ) 2 2 f (x) = Aexp − x / 2 ， ( ) 2  =1 的随机点。抽样中常数 A 的值需要知道吗？试决定接受点与试探步数之比，到达平 衡分布的时间与最大试探步长  的关系。（提示：判断“平衡”的标准是 2 2  x   ）。  选多大较合理? （16） 用 Metropolis 随机游走的方法计算积分 , (0 4) 4 0 2    − x e dx x x 。 （17） Laplace 方程及其边界条件为 ( )  ( ) ( ) ( ) ( )   = = = =  = ,0 ,1 0, 0, 1, 1 , 0 2 x x y y x y      ， 用随机游走的蒙特卡洛方法数值求解正方形场域 (0  x 1, 1 y 1) 的势函数。 第三章 习题 （1）利用蒙特卡洛方法计算三维、四维、五维和六维空间的单位半径球的体积

(2)利用分布密度函数f(x)=Ae做重要抽样来求积分,并分析误差与投点数的关系。 dx (3)用事例证明蒙特卡洛求积分的标准误差为 其中A为物理观测量,N为蒙特卡洛投点个数 (4)采用 Metropolis方法产生一维分子速度分布密度函数为 f(v)=Cv'e 的游走样本点,并将其分布和上述分布函数曲线进行比较(上式中C,a为常数)。 (5)写出采用 Metropolis方法对高斯分布f(x)=Aexp(-x2/2a)的抽样框图和程序 (A和G为常数)。 (6)编写采用 Metropolis算法计算一维积分 4-x2/2x2 的程序。用该程序计算三维积分 xdx[y dy[=2exp (-r2/2)d23 然后将得到的结果与解析计算得到的精确结果进行比较,并分析模拟游走点数与误差 的关系。 (7)对以下一维扩散方程 U(x,0)=f(x) a-0 aU C(0,1)=U 0(1, 6=0 可以通过在h×τ的矩形格点上的随机游走来求解(其中h,τ分别为x和t划分的格 点长度在x方向向前和向后游走的几率为on1=on3=(2+h2/x),而在t方向 向前和向后游走的几率为n2=o4=(+2x/h2)。试编程予以计算 (8)编写程序,采用路径积分量子力学蒙特卡洛方法求液态He4的基态能量。 (9)编写程序,采用变分量子力学蒙特卡洛方法求氢分子的基态能量。其中两质子和两电 子应当按四体系统来处理。 (10)修改习题(8)的程序,采用格林函数量子蒙特卡洛方法求氢分子基态能量,并与习 题(8)的结果进行比较

（2）利用分布密度函数 x f x Ae− ( ) = 做重要抽样来求积分，并分析误差与投点数的关系。  + − = 0 5 / 2 I x e dx x （3）用事例证明蒙特卡洛求积分的标准误差为 N A A 2 2 2 1  = −  ， 其中 A 为物理观测量，N 为蒙特卡洛投点个数。 （4）采用 Metropolis 方法产生一维分子速度分布密度函数为 2 2 ( ) v f v Cv e − = 的游走样本点，并将其分布和上述分布函数曲线进行比较（上式中 C, 为常数）。 （5）写出采用 Metropolis 方法对高斯分布 ( ) exp( / 2 ) 2 f x = A −x  的抽样框图和程序 （ A 和  为常数）。 （6）编写采用 Metropolis 算法计算一维积分 ( )  + − − x x dx 2 2 exp / 2 的程序。用该程序计算三维积分 ( )    + − + − + − x dx y dy z exp − r / 2 dz 2 2 2 2 . 然后将得到的结果与解析计算得到的精确结果进行比较，并分析模拟游走点数与误差 的关系。 （7）对以下一维扩散方程      = = =   =   1 2 0 2 (1, ) (0, ) ( ,0) ( ) , U t U U t U U x f x t U x U  ， 可以通过在 h 的矩形格点上的随机游走来求解（其中 h ， 分别为 x 和 t 划分的格 点长度）。在 x 方向向前和向后游走的几率为 ( ) 1 2 01 03 2 / −  =  = + h  ，而在 t 方向 向前和向后游走的几率为 ( ) 1 2 02 04 1 2 / −  =  = +  h 。试编程予以计算。 （8）编写程序，采用路径积分量子力学蒙特卡洛方法求液态 4 He 的基态能量。 （9）编写程序，采用变分量子力学蒙特卡洛方法求氢分子的基态能量。其中两质子和两电 子应当按四体系统来处理。 （10）修改习题（8）的程序，采用格林函数量子蒙特卡洛方法求氢分子基态能量，并与习 题（8）的结果进行比较

第四章习题 (1)用有限差分法发展一个程序,数值求解正方形场域(0≤x≤1,0≤y≤1)的拉普拉斯 方程, o(x,0)=o(x,)=0,o0,y)=o1,y)=1 (2)用有限差分法发展一个程序,数值求解极坐标下的泊松方程 Vo(,O)=-4(,O) o(a,0)=1,(b,0)=V2 然后,选择p(,)=0,边界条件a=1,b=2,V=0,V2=1时,两个圆圈 中间的势分布。 (3)第(2)题中若 p(r,0)={10 exp[(r-a)/al a<r<b 0, 其它 其余取值相同,数值求解两个圆圈中间的势分布 (4)在一个二维L×L的反应堆中,中子的扩散方程为 L八叫(L 用hn=h,=h的正方形网格离散化后,证明它的有限差分方程满足 TTh m = SIn L L h 2L 第五章习题 (1)公式(5.2.20)是更一般的公式 NN Nm dxdy=An△ 的特殊情况。试给出证明 (2)用有限元素法发展一个程序,数值求解正方形场域(0≤x≤10≤y≤1)的拉普 拉斯方程, o(x,0)=o(x)=0,g(0,y)=o(L,y)=1

第四章 习题 (1) 用有限差分法发展一个程序，数值求解正方形场域 (0 1, 0 1     x y ) 的拉普拉斯 方程， ( )  ( ) ( ) ( ) ( )   = = = =  = ,0 ,1 0, 0, 1, 1 , 0 2 x x y y x y      。 （2）用有限差分法发展一个程序，数值求解极坐标下的泊松方程 ( ) ( ) ( ) ( )    = =  = − 1 2 2 , , , , 4 , a V b V r r          。 然后，选择 (r,) = 0 ，边界条件 a =1, b = 2, V1 = 0, V2 =1 时，两个圆圈 中间的势分布。 （3）第（2）题中若 ( )  ( )      − −   = 0， 其它 exp / , 10 1 , r a a a r b  r  ， 其余取值相同，数值求解两个圆圈中间的势分布。 （4）在一个二维 L  L 的反应堆中，中子的扩散方程为 ( , ) sin sin 0 2 2 2 2 2  =            +         −   +   L y L x a x y x y    ， 用 hx = hy = h 的正方形网格离散化后，证明它的有限差分方程满足 1 2 2 2 2 sin 8 sin sin −              +            = L h h a L j h L i h ij     。 第五章 习题 （1） 公式(5.2.20)是更一般的公式 ( )  + + +  = e n m l j k i k l n k l n N N N dxdy 2 ! ! ! ! 的特殊情况。试给出证明。 （2） 用有限元素法发展一个程序，数值求解正方形场域 (0 1, 0 1     x y ) 的拉普 拉斯方程， ( )  ( ) ( ) ( ) ( )   = = = =  = ,0 ,1 0, 0, 1, 1 , 0 2 x x y y x y     

(3)修改(2)题中的程序,仍然采用有限元素法数值求解三角形场域 D∈(0≤x≤1,y≤1-x)的拉普拉斯方程 V2(x,y)=0 p(x:0)=0,o(0,y)=o(x,y=1-x)=1 第六章习题 (1)写用 Verlet速度算法求解三维分子运动方程的程序。 (2)编写一个三维,元胞尺寸为L的周期边界条件计算程序。 (3)试做总能量固定的单原子系统的分子动力学模拟。元胞为Lx=L,=L2=10,划分 为10×10×10的正方形网格。元胞内原子数N=64。原子质量m=1。位势为 Lenard- Jones势,其中ε=σ=1,边界条件为周期性边界条件,初始位置是随机 分布在正则节点上,初始速度为按[-1,1随机分布。分子动力学模拟步长取为 M=002,模拟100-200步后原子的速度分布和位置分布如何? (4)试做二维单原子系统的分子动力学模拟。系统温度T=085保持固定,模拟参数及 其他条件同上题。 第七章习题 (1)用 Mathematica语言定义求解一元二次方程ax2+bx+c=0的函数,该函数还要求 能处理各种常数a,b,c的情况 (2)用 Mathematica语言定义一个能画出任意给定n值的正n边形的函数 (3)用 Mathematica语言定义一个操作函数,它可以在给定二维矩形xy平面区间,按 给定步长h,b划分矩形网格,并列出节点的坐标表{x, (4)接着上题,用 Mathematica语言定义一个操作函数,对确定步数n,生成一个在x-y 平面格点上n步服从均匀分布的随机游走的图 (5)用 Mathematica语言实现一个产生任意阶的勒让德多项式的 Mathematica程序包。勒 让德多项式的递推公式为: 0(x)=1,P(x)=x,Pn1(x)=[(2n+1)xP(x)-mP1(x)n+1) (6)用 Mathematica语言实现一个产生任意阶的埃米尔特多项式的 Mathematica程序包。 埃米尔特多项式的递推公式为: Ho(x)=1, H,(x)=x, H+(x)=2xH, (x)-2nHn-(x) (7) Mathematica语言编写一个从某点出发求多元函数的局部极小或极大值的程序包。 (8)用 Mathematica语言编写一个程序包,它能实现平面图形的(a)平移,(b)旋转,(c) 对x座标轴的反射

（3） 修改（ 2 ） 题 中 的 程 序 ， 仍 然 采 用 有 限 元 素 法 数 值 求 解 三 角 形 场 域 D(0  x 1, y 1− x) 的拉普拉斯方程， ( )  ( ) ( ) ( )   = = = − =  = ,0 0, 0, , 1 1 , 0 2 x y x y x x y     第六章 习题 （1） 写用 Verlet 速度算法求解三维分子运动方程的程序。 （2） 编写一个三维，元胞尺寸为 3 L 的周期边界条件计算程序。 （3） 试做总能量固定的单原子系统的分子动力学模拟。元胞为 Lx = Ly = Lz = 10 ，划分 为 101010 的正方形网格。元胞内原子数 N = 64 。原子质量 m =1 。位势为 Lenard-Jones 势，其中  = =1 ，边界条件为周期性边界条件，初始位置是随机 分布在正则节点上，初始速度为按[-1，1]随机分布。分子动力学模拟步长取为 t = 0.02 ，模拟 100-200 步后原子的速度分布和位置分布如何？ （4） 试做二维单原子系统的分子动力学模拟。系统温度 T = 0.85 保持固定,模拟参数及 其他条件同上题。 第七章 习题 （1） 用 Mathematica 语言定义求解一元二次方程 0 2 ax + bx + c = 的函数，该函数还要求 能处理各种常数 a,b,c 的情况。 （2） 用 Mathematica 语言定义一个能画出任意给定 n 值的正 n 边形的函数 （3） 用 Mathematica 语言定义一个操作函数，它可以在给定二维矩形 x-y 平面区间，按 给定步长 hx hy , 划分矩形网格，并列出节点的坐标表 xi , yi。 （4） 接着上题，用 Mathematica 语言定义一个操作函数，对确定步数 n，生成一个在 x-y 平面格点上 n 步服从均匀分布的随机游走的图。 （5）用 Mathematica 语言实现一个产生任意阶的勒让德多项式的 Mathematica 程序包。勒 让德多项式的递推公式为： P0 (x) =1， P (x) = x 1 ， ( ) (2 1) ( ) ( ) ( 1) Pn+1 x = n + xPn x − nPn−1 x n + . （6）用 Mathematica 语言实现一个产生任意阶的埃米尔特多项式的 Mathematica 程序包。 埃米尔特多项式的递推公式为： H0 (x) =1， H (x) = x 1 ， H x xH (x) nH (x) n 1 2 n 2 n 1 ( ) + = − − （7） Mathematica 语言编写一个从某点出发求多元函数的局部极小或极大值的程序包。 （8）用 Mathematica 语言编写一个程序包，它能实现平面图形的（a）平移，（b）旋转，(c) 对 x 座标轴的反射

第八章习题 (1)画出不同主量子数、轨道量子数n,下,氢原子径向部分波函数随r的变化图形,并 讨论原子序数Z变化的作用 (2)画出不同轨道量子数、磁量子数、l,m下,氢原子O,q部分波函数随b和q的变化的 维图形。 (3)采用诺曼诺夫( Numerov)法(参见附录C),编写一个程序求解电子的一维薛定格方程 的最低的两个能量本征值和波函数。该电子所在势阱的势函数为(我们选择原子单位 方 l0x,0<x<5 V(x) 50.其它 (4)编写求解不同势阱高度和宽度的一维势阱薛定格方程波函数和能量谱的程序包 势阱的势函数为 Vo x< (x) 0其它 (5)试用 Schroedinger.m程序计算()=r2时薛定格方程的基态,第一激发态的能量值, 并与变分法求得的结果进行比较

第八章 习题 （1）画出不同主量子数、轨道量子数 n,l 下，氢原子径向部分波函数随 r 的变化图形，并 讨论原子序数 Z 变化的作用。 （2）画出不同轨道量子数、磁量子数、 l,m 下，氢原子 , 部分波函数随  和  的变化的 三维图形。 （3）采用诺曼诺夫(Numerov)法（参见附录 C），编写一个程序求解电子的一维薛定格方程 的最低的两个能量本征值和波函数。该电子所在势阱的势函数为(我们选择原子单位 a.u.，即 me = e =  = c =1 )      = 50, 其它 10 , 0 5 ( ) x x V x 。 （4）编写求解不同势阱高度和宽度的一维势阱薛定格方程波函数和能量谱的程序包。 势阱的势函数为   −  = 0 其它 ( ) V0 x a V x 。 （5）试用 Schroedinger.m 程序计算 2 V(r) = r 时薛定格方程的基态，第一激发态的能量值， 并与变分法求得的结果进行比较