第八章 Mathematica在量子力学中的应用举例 § 8.1 粒子在中心力场中的运动问题 设电子与原子核的约化质量为 Mm Mm e e + =μ , r Ze r 2 )( ?=V ，哈密顿量为 )( 2 )?( 2 ? ? 2222 rVrV p + ? ?=+= μμ = G = H . (8.1.1) 其中 r 为粒子 所处的空 间 位置到中 心 势原点的 距 离 。 利 用 中 心 势的球对 称 性 ， 我 们 将薛定 格 方程写为在球坐标中的表示 ()= ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? +? ? ? ? ? ? ? ? ? ? +? ? ? ? ? ? ? ? ? ? ? ??ψ ?θθ θ θθμ ,, sin 1 sin sin 1 2 2 2 2 2 2 2 r r r rr = ( )( )( )??ψ ,,rrVE? . (8.1.2) 角动量算 符 的定义为 ： L px GG G ×= 。 可以证 明 [ ， 所 以 角动量 是 守 恒量， 即 在中 心 力场中运动 粒子的一个 重要特征是 角动量守恒 。由此可以 得到 0] ? , ? =HL L ? 2 ? L （角动 量的平方）也 是 守 恒量。 在求 解中心力 场 作用下粒 子 的能量本 征 方程时， ( ) z L ? LH , ? , ? 2 构成对易算 符 的一个 完全 集 。 ? ? ? ? ? ? ?? ? ? ? ? ? ? ? ? ? =?≡ 2 2 2 2 2 ? 1 = L r r rr ? . (8.1.3) 其中在球坐标中的角动量平方算符可以表示为： ? ? ? ? ? ? ? ? + ? ? ? ? ? ? ? ? ? ? ?= 2 2 2 22 sin 1 sin sin 1 ? ?θθ θ θθ =L . (8.1.4) 薛定格方程 (8.1.2)则可以写为 ()()(?θψ?θψ μ ,,,, ? 2 2 2 2 2 2 rVEr L r r rr ?= ? ? ? ? ? ? ?? ? ? ? ? ? ? ? ? ? = = ) ) ? . (8.1.5) 波函数 (?θψ ,,r 与极 角 θ（ 2/2/ πθπ ≤≤? ） 和方位 角 ? )0( π? ≤≤ 的关联 是 由算符 和 决定的。假定满足薛定格方程的本征波函数 2 ? L z L ? ( )?θψ ,,r 可以 分离变量 表示为 ()( )( ) ()( ) ( )???θ?θψ ΦΘ≡≡ rRYrRr ,,, . (8.1.6) z L ? 在球坐标系中可以表示为： ?? ? ?= =iL z ? . 该算符的本征值由求解本征方程 () ()?? ? Φ=Φ ? ? ? z Li= ， (8.1.7) 来得到。方程（ 8.1.7）的解为 ( ) =/? ? z iL Ae=Φ . (8.1.8) 由于（ 8.1.8）式所示 波 函数解必 须 唯一确定 ， 因而它也 必 定满足条 件 ： () ( )?π? +Φ= 2Φ ， 并 且 角 动量算 符 的本征 值应当是 离 散的， 其本 征 值 表 示 为 ： L z L ? =m z = ， ( . 由本征波函数的归一化条件，方程（ 8.1.7）归一化的解可以写为 ,...)2,1,0 ±±=m () ? π ? im e 2 1 =Φ . (8.1.9) 类似地，对另一个守恒量 -角动量平方，我们有本征方程： () () (?θ?θ ?θθ θ θθ ?θ ,, sin 1 sin sin 1 , ? 2 2 2 2 22 YLYYL = ? ? ? ? ? ? ? ? +? ? ? ? ? ? ? ? ? ? ?= = ) . (8.1.10) 方程 (8.1.10)的解是球谐函数 Y 。如果本征值满足 L ，方程 (8.1.10)写为 ml, 22 )1( =+= ll ()0,)1( sin 1 sin sin 1 , 2 2 2 = ? ? ? ? ++ ? ? +? ? ? ? ? ? ? ? ? ? ?θ ?θθ θ θθ ml Yll ? ? . (8.1.11) 角动 量算符 作用 在 球谐函数 Y 上 的本征值 由角量子数 2 ? L ml, ,...2,1,0=l m 决 定。对应 于确定的 角量子数 l，算 符 的本征 值则为 l ， 此 时磁量子 数 则描写该 角 动量在 轴上 的 投影，它的取值范围为： 2 ? L 2 )1( =+l l z m ±±±1,0= ,...,2, m 。这就是说：对确定的角动量量子数 ，应当 有 个本征函数 Y 。对磁量子数 为正时的情况，球谐函数的完整表达式为 l 12 +l ml, ()() () () ( ) () ? θ π ?θ imm l m ml eP l ml ml cos 4 12 ! ! 1, , + + ? ?=Y . (8.1.12) 其中 为 l阶的 第 个 伴随 勒 让德函数 。 如果磁量 子 数为负时 （()xP m l m m? ） ， 其 球 谐函数 满 足如下关系式 ()() ( ) () (?θ?θ , ! ! 1, * ,, ml m ml Y ml ml + ? ?= ? )Y . (8.1.13) 显然， 球谐 函 数 Y 也是算 符 的本征函 数。 容易证明 类 似 (8.1.7)式 ， 球谐函数 Y 满足 ： ml, z L ? ml, mlmlmlz YmYiYL ,,, ? == = ? ? ?≡ ? . (8.1.14) 因而球 谐函 数 Y 既是 角动 量算符 平方 ，也是 角动 量算符 的 分量 的 本 征 函数。在 Mathematica 中 球谐函 数表 示为 SphericalHarmonicY[] 。勒 让德 多项 式表 示为 LegendreP[]。 ml, 2 ? L z z L ? 将（ 8.1.6） 式代入薛定 格方程 (8.1.2), 再应用 上面推导出 的角动量部 分波函数所 满 足的薛定 格 方程，可 以 得到本征 波 函 数 ( )?θψ ,,r 表示中 的径向部 分 ( )rR 应当满足 的 方 程 。 0 )1(22 2 2 22 2 = ? ? ? ? ? ? + ? ? ? ? ? ? ? +++ R r ll r Ze E dr dR rdr Rd = μ . (8.1.15) 下面我们以 氢原子为例 进行分析。 定义波尔半 径 m em a e 11 2 2 0 1029.5 ? ×≈= = 为长度单 位，即 0 /ar=ρ ; 以 氢原子 的电离能 量 eV em a e E e 5.13 2 4 4 0 2 0 ≈== = 为 能 量 单位， 即 0 EE=ε ； 定 义径向函数 ρρρ /)()( uR = 。这时方程 (8.1.15)写为 0)( )1(2)( 22 2 = ? ? ? ? ? ? + ?++ ρ ρρ ε ρ ρ u llZ d ud . (8.1.16) 能量 ε 的值是由方程 (8.1.16)的本征值和本征函数决定的。 我们 考 虑稳定状 态 （束缚态） ，即 0<ε 的状态。 分析表明 函 数 )(ρu 可以表 示 为多项式 或者指数形 式。为了找 出 )(ρu 的近似式 ，我们通过 考察它在 r 和0→ ∞→r 时的极限行 为，发现由波函数的幺正性条件要求上述两种表达方式下都可以推出 u . (8.1.17) )()( 1 ρρρ γρ l l fe ?+ = 将 (8.1.17)式代入 (8.1.16)后，求解得到 超几何函数 ( ) 11 F 形式的解。 () ? ? ? ? ? ? ? ? +?+= γρ γ ρ 2;22,1 11 l Z lFcf l . (8.1.18) 其中 εγ ?≡ 。现在我们由式 (8.1.17)得到电子在库仑势中的波函数的径向部分为 ? ? ? ? ? ? +?+= ? ρρρ ρ n Z lnlFeNR nZl ln 2 ;22,1)( 11 / , . (8.2.19) 由于归 一化 条件的 要求 ， (8.1.18)的级数 表示 必须只 有有 限项。 这个 限制就 给出 了能量 的 值 ? ? ? ? ? ? ? ? ?+?= γ Z ln r 1 ， ( . (8.1.20) ,...)2,1,0= r n 由此我们得到 1++ = ln Z r γ . (8.1.21) 由 γ 和 ε 的定义，则 () 2 2 0 2 2 0 1 n ZE ln ZE E r ?= ++ ?= . (8.1.22) 其中 为主量 子数 。它是 由径向量子 数 n 和轨道角 动量量子 数 决定的。在 这里我们引 入一组称为 “拉盖尔 (Laguerre)多项式 ”的特殊正交 多 项式 [2]，拉盖尔多项式式由级数定义为 n ,0= (γ k L ,...)2,1( =n r ,...)2,1,0( = r n l ,...)2,1(l ) () () ( ) ! 1 0 j x jk k xL jk j j k ? ? ? ? ? ? ? ? ? + ?= ∑ = γ γ . 相应的归一化为 () () () () () ( ) kkkk k k xLxLxxdx ′′ ∞ ++Γ =? δ γ γγγ ! 1 exp 0 ∫ . 超几何函数与拉盖尔多项式间有如下关系式 () () ()xnF n n xL n ;1, 1! 1 )( 11 )( +? +Γ ++Γ = α α α α . 这样电子在库仑势中的波函数的径向部分的解也可以写为 ? ? ? ? ? ? ′= + + ? ρρρ ρ n Z LeNR l ln nZl ln 2 )( )12(/ , . （ 8.1.23） 径 向 部 分波函 数 (8.1.23)中的 拉盖尔 (Laguerre)多项式的 性 质见文献 [2]。 相应的 波 函数为 () (?θρρ?θρψ ρ , 2 ;22,1,, ,11 / ,,, ml nZl lnmln Y n Z lnlFeN ? ? ? ? ? ? +?+= ? ) . (8.1.24) 在 (8.1.19)和 (8.1.24)式中的归一化常数为 () () () 2/3 , 2 !12 ! !12 1 + ? ? ? ? ? ? ?? + + = l ln n Z lnn ln l N . (8.1.24) (*------------------------------------------------------------------ Mathematica Package file Coulombp.m ------------------------------------------------------------------*) BeginPackage["CoulombPotential`"] Clear[WaveF,WaveR,WaveA]; WaveF::usage = "WaveF[Z_, r_, theta_, phi_, n_, l_, m_]计算电子在库仑 势中本征波函数的表示。 Z 为原子核的电荷数 . r为电子到中心势原点的距离 . theta 和 phi 为球坐标中的角度 . n, l和 m为能量和角动量算符的量子数。 " WaveR::usage = "WaveR[Z_, r_, n_, l_] 计算电子在库仑势中的本征波函数 径向部分的表示。 Z 为原子核的电荷数 . r为电子到中心势原点的距离 . n和 l 为能量和角动量算符的量子数。 " WaveA::usage = "WaveA[theta_, phi_, l_, m_]计算电子在库仑势中本征波函 数的角度关联部分表示。 theta 和 phi 为球坐标中的角度 . l 和 m 表示角动 量算符的量子数。 " (* --- 定义公共变量 --- *) r::usage n::usage l::usage m::usage theta::usage phi::usage Begin["’ Private’ "] (* --- 产生库仑势中波函数的径向部分 --- *) WaveR[Z_, r_, n_, l_] := Module[{unit, tmp}, (* --- 归一化常数 --- *) unit = (Sqrt[(n + l)!/(2 n (n - l - 1)!)] ((2 Z)/n)^(l + 3/2)) /(2 l + 1)!; (* --- 产生波函数径向部分的定义 --- *) tmp = unit r^l Exp[-((Z r)/n)] Hypergeometric1F1[l + 1 - n, 2 l + 2, (2 Z r)/n] ] (* --- 产生库仑势中本征波函数的角度相关部分 --- *) WaveA[theta_, phi_, l_, m_] := Module[{tmp}, tmp = SphericalHarmonicY[l, m, theta, phi] ] (* -- 产生电子在库仑势中的本征波函数 --- *) WaveF[Z_, r_, theta_, phi_, n_, l_, m_] := Module[{tmp}, tmp = WaveR[Z, r, n, l] WaveA[theta, phi, l, m] ] End[] EndPackage[] 当我们需要对电子在原子核的库仑势中的本征 波函数习性进行分析时， 我们 可以首先调入程序包 Coulombp.m，然后调用程序包中定义的函数。例如通过运 行下面的指令： (*----------------------------------------------------------------*) << Coulombp.m Plot[WaveR[1,r,1,0],WaveR[1,r,2,0],WaveR[1,r,3,0],WaveR[1,r,4,0], {r,0,35},AxesLabel->"r","u",Prolog->Thickness[0.001]] Plot[Abs[WaveA[theta,Pi/2,2,1]]^2, {theta,0,Pi},AxesLabel->"theta","Y",Prolog->Thickness[0.001]] Plot3D[Abs[WaveF[1,r,theta,Pi/2,3,2,2]]^2,{r,0,15},{theta,0,Pi},Ligh ting->True] (*---------------------------------------------------------------*) 图（ 8.1.1）本征波函数径向部分的四条曲线（当 1=Z ， 0=l 和 4,3,2,1=n 时） 。 图（ 8.1.2） 本征波函 数 角度关联 部 分绝对值 平 方随极角 θ 变 化的曲线 （ 当 2π? = ， l 和 1.3） 本 2= 1=m 时） 。 图（ 8. 征 波 函 数 绝对值平 方 随 r 和极角 θ 变 化的三维 曲 线 （当 2π? = ， 和 2=l 1=m 时） 。 8.2 求非相对论性薛定格方程本征能量限 用 Mathematica V4.0 系统的指令，对应的计算过程可表述为： MATHEMATICA V4.0 （ * 积分 *） In[1]:= E ∫ drr nr ∞ ? 0 λ Out[1]= & ]],1[,0]Re[&1][Re[ 0 1 drrenGammanIf nrn ∫ ∞ ??? +>?> λ λλ 这 一 则输 出 结果 的含义 是： 如果 0]Re[ >λ , 且 1]Re[ ?>n ( ) 1 1 n nλ ?? Γ+ ， 以 上 积分 的结果 为 ，否则将输出 dr r r n ∫ ∞ 0 }exp{λ , 这意味着 Mathematica无法求解该问题。由此可以得归一化因子 =N ， π λ 2/3 归 一 化 的 “ 试 验 ” 波 函 数 为 (rΨ , ) r e λ π λ λ ? = 2/3 为保险 起见 ，我们 可以 检验一 . 下关 系式两 边的 量纲。 根据 以前的 讨论 ，我们 知道 关系式 左 边的量 纲为 。为使 指数 运算][ 2/3 EDim }exp{ rλ? 有意义 ，乘积 rλ 必须 是无量 纲的 量，即 1][ =rDim λ 。由此有 ][Eim ，即 ][ ][ D rDim Dim ==λ 1 。 很显然， 在 以上推导 中 至少量纲 是 正确的。 下 面我们演 示 一下如何 运 用 Mathematica 语言 作以上定义和计算。 ][][][ 2/32/3 λDimEDimDim ==Ψ 采用 Mathematica V4.0 的对应计算为： MATHEMATICA V4.0 In[2]:= （ * 积分 *） drEr r ∫ ∞ ? 0 22 4 λ π Out[2]= ], 4 1 ,0][Re[4 0 22 3 drreIf r ∫ ∞ ? > λ λ λπ 输出的含 义 是：当 Re 0>λ 时， 计算结果 为 3 λπ ,否则， Mathematica 无法 求解，将 返 回 输 入形式 。完整的结果应是： dr 2 re r 0 2 4 ∫ ∞ ? λ π 1 4 4 3 2 = λ π N ? π λ 3 =N . 下一步 ，我 们将借 助引 入的“ 试验 ”波函 数求 动能项 的期 望值。 由于 我们只 讨论 基态的 能 量本征值 ， 而对基态 量 子 数 l ， 此 时 在径向中 心 力场势情 况 下可采用 拉 普拉斯算子形 式 为 0= dr d rdr d 2 2 2 2 +=?≡? . 其期望值为 ()() rr e dr d rdr d xedrrxd λλ π λ λλ ?? ? ? ? ? ? ? ? ? +=?ΨΨ ∫ 2 ,, 2 2 3 3 *3 ∫ () r errdr λ λλπ π λ 2 0 22 3 24 ? ∞ ∫ ?= ( ) () ( ) () ? ? ? ? ? ? Γ ? Γ = 23 23 2 2 2 2 3 4 λ λ λ λλ 23 4 1 4 λ λ λ ?=? ? ? ? ? ? ?= . 我们可以看到这里的量纲检验仍然是正确的。我们在下式中省略了 Dim[… ]符号 : 222/322/33 2 3 1 λ→=→ΨΨ ? EEEEE x x . 动能项的期望值为 () () μ λ λ μ λ μ 2 , 2 , 2 2 *3 2 =Ψ ? ? ? ? ? ? ? ? ? ?Ψ= ∫ rrxd p G . (8.2.12) 相应的 Mathematica V4.0 计算过程为 MATHEMATICA V4.0 In[3]:= [] r Er λ π λ λψ ? = 2/3 :__, （ * 定义“试验”波函数 *） In[4]:= [][ ][]}1,{],,[ 2 }2,{],,[:__, rrD r rrDrg λψλψλ += （ * 轨道角动量为零的有效拉普拉斯算符 *） In[5]:= drrgrr ],[],[4 0 2 λλψπ ∫ ∞ Out[5]= ] 2 , 4 ,0][Re[ 0 2/72/5 22/3 dr e r e reIf rr r ∫ ∞ ?? ? ? ? ? ? ? ? ? ? +??> π λ π λ π λ λλπ λλ λ 4 Mathematica 表达式 D[psi[r,lambda],{r,n}]功能为 ， 以 r 为 变量对 ( )Ψ 求 n 次偏导。 幂 指数势函数 V 的期望值为 λ,r () n arr = ()()() ( ) () 3 3 0 22 3 *3 2 3 44,, + ∞ ?+ +Γ ==ΨΨ ∫∫ n rn n adrerarrVrxd λ π π λ π π λ λλ λ . 即 () ()()() ( ) () 3 3*3 2 3 4,, + +Γ =ΨΨ≡ ∫ n n arrVrxdrV λ λλλ . (8.2.13) 量纲分析要求 V → n ar= n aE ? → E 即耦合常数 的量纲为 a a → 1+n E . 则方程 (8.2.13)的量纲也是正确的 , 即 a n? λ → EEE nn = ?+1 . 对应的 Mathematica V4.0 的指令为： ------------------------------------------------------------------------------ MATHEMATICA V4.0 ------------------------------------------------------------------------------ In[6]:= (* 积分 *) drEr rn λ2 0 2 ? ∞ + ∫ Out[6]= ]],3[2,3]Re[&&0][Re[ 2 0 233 drrenGammanIf nrnn + ∞ ????? ∫ +?>> λ λλ ------------------------------------------------------------------------------ 结合方程 (8.2.12)和 (8.2.13)，可得能量表示为 () ( () ) n na λμ λ λ 2 3 22 2 +Γ +=E . （ 8.2.14） 对于任何 0>λ 的值 ， 这 一能量 解始终是 能 量值“真 解 ”的上界 ： ( )E 。 通过求可 使 λE true ≤ ( )λE 取最小值 的 变分参数 λ 值， 即解出 min λ ， 就可 以很容易 地 改进这一 能 量上限， 使其 逼近真解 。 显然， min λ 由下式给出 ( ) 0= ? ? λ λE . 对公式（ 8.2.14）求偏导，可得 () () () 0 2 3 1 = +Γ ?= ? ? +n n an E λμ λ λ λ . 求解该方程得到 () )2/(1 1 min 2 3 + + ? ? ? ? ? ? +Γ = n n nanμ λ . (8.2.15) 将这一结果反代回关系式 (8.2.14)，即可得改进后的能量本征值上限 () () ? ? ? ? ? ? + ? ? ? ? ? ? +Γ ? ? ? ? ? ? ? ? == + + + n nan EE n n nn 2 1 2 31 2 1 )2/(2 1 )2/( minvar μ λ . (8.2.16) 对应的 Mathematica V4.0 的程序为： ------------------------------------------------------------------------------ MATHEMATICA V4.0 ------------------------------------------------------------------------------ In[7]:= e (* 定义函数 n nGammaa )2/(]3[2/)2/(:_][ 2 λμλλ ++= ( )λE *) In[8]:= D[e[λ],λ] （ * 对参数 λ 作微分（注意 :D[e[λ],{λ,1}]等效于 D[e[λ],λ]） ): *) Out[8]= λ μ ?2 ?1?n anλ ?1?n Gamma@3+ nD (* 解方程求 min λ *)} In[9]:= ],0]3[2/[ 11 λλμλ ==+? ???? nGammaanSolve nn Out[9]= ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? +→ + ?? n n nGammaan 2 1 1 )]3[2( μλ In[10]:= e （ * 计算 ]])3[2[( )2/(11 nn nGammaan +?? +μ ( ) min λE *) Out[10]= ]3[])3[2(2 2 ])3[2( 2 1 11 2 2 1 nGammanGammaana nGammaan n n nn n n + ? ? ? ? ? ? ? ? ++ + ? + ???? + ?? μ μ μ （ * 注意：指令 PowerExpand[expr] 的功能为将所有乘积和指数作幂次展开。 % 代表 Mathematica 输出的最 后 的一个表 达 式，在此 即 为上面最 后 一个表达 式 ， 即 Out[11]}。 *） In[11]:= PowerExpand[%] Out[11]= () ( ) nn n n n nn n nn n nnn n nGammananGammana ++ ? + ? ++ ?? ++ ? +++ ?? +++ 2 2 222 2 2 22 2 2 22 2 2 2 2 34 ]3[2]3[ μμ2 (*--------------------------------------------------------------------------*) 通过 仔 细的比较 ， Mathematica V3.0 输出 Out[12]和 Mathematica V4.0 输出 Out[11] 给出的结 果 与公式（ 8.2.16）是 一致 的。下面 我 们 将 (8.2.16)式应 用 于一些 典 型 的势函 数 上。 库仑势： 将 α?=a 和 n 代入，可得 1?= 2 2 μα ?≤ true E . 可以看 出， 表达式 的右 边正好 是氢 原子基 态能 量，即 等号 是严格 成立 的。显 然， 这是由 于 我们恰好 选 取氢原子 基 态波函数 作 为“试验 ” 波函数引 来 的 。 对 于 1=α 和 1=μ 情况， 能量 数值解为 。这一处理的 Mathematica V3.0 表述式为 5.0?=E 与上面指令对应的 Mathematica V4.0 版本如下： ------------------------------------------------------------------------------ MATHEMATICA V4.0 ------------------------------------------------------------------------------ In[13]:= e[λ_,n_,a_,μ _] := n nGammaa )2/(]3[2/)2/( 2 λμλ ++ (* 定义需要最小化的函数 *) In[14]:= FindMinimum[e[λ,-1,-1,1], {λ,0.5}] (* 以变分参数 5.0=λ 为起始点，求最小值 *) Out[14]= {-0.5, {λ 1.} → （ * 这可以理解为 1 min =λ 时，能量最小值 *） 5.0 min ?=E 线性势： V 。将 代入关系式（ 8.2.16） ，得到 () arr = 1=n 3/1 2 3/5 var 2 3 ? ? ? ? ? ? ? ? ? ? ? ? ? ? =≤ μ a E true E . 将这一结果按关系式 (8.2.6)的格式重写出来，有 3/1 2 3/1 2 3/4 3/5 2 4764.2 22 3 ? ? ? ? ? ? ? ? = ? ? ? ? ? ? ? ? ≤ μμ aa E true . (8.2.17) 而基态能量的 “真解 ” 可由 Airy 函数 [8][9]的第一个零点给出， true E 3/1 2 2 3381.2 ? ? ? ? ? ? ? ? = μ a true E . 比较两种结果，可以看到我们求得的原始上限 与真值的相对误差非常小 . var E %6 var ? ? true true E EE . 由此我 们在 不具体 求解 薛定格 方程 的情况 下， 解析地 推导 出线性 势基 态本征 能量 ，其数 值 结果与 “真解 ”有很好的近似。 对应的 Mathematica V4.0 版本为： MATHEMATICAV4.0 In[15]:= e[λ_,n_,a_,μ _] := n nGammaa )2/(]3[2/ )2/(2 λλ μ + (* 定义能量函数 *) In[16]:= FindMinimum[e[λ,1,1,1],{ λ,0.5}] （ * 以变分参数 5.0=λ 为起始点，求最小值。 *） Out[16]= {1.96556,{ λ -> 1.14471}} 将 代入 公式 (8.2.15),(8.2.16)和 (8.2.17)，可对以上 结果进行检 验。通常将 计算结果绘 图显示，有 助于进行比 较。令 1== an 12 =μ ，我 们分别绘制 函数 和 ，然后借助 Mathematica V3.0 的 Show 指令将两图形合在一起进行比较。 3/2 3381.2 aE true = 3/2 4764.2 a var E = MATHEMATICA V3.0 定义函数 In[17]:= etrue[a_] := 2.3381 a^(2/3) In[18]:= eupper[a_] := 2.4764 a^(2/3) In[19]:=plot1=Plot[etrue[a],{a,0,3}, AxesLabel->{”a”,”E”}, TextStyle-> {FontSlant->”Italic”,FontSize->14}] (* 绘制 (plot1) *) true E （ * AxesLabel: 定义坐标轴的表述； TextStyle 定义字型和字体大小。 *） 图 8.2.2 线性势基态能量真值 。 In[20]:=plot2=Plot[eupper[a],{a,0,3},AxesLabel->{”a”,”E”},TextStyle-> 2) *) true E {FontSlant ->”Italic”,FontSize->14}] (* 绘制 var E (plot 图 8.2.3 线性势基态的变分上限能量 。 var E In[21]:= Show[plot1,plot2] (* 将 {\ tt plot1}与 {\ tt plot2}合并绘制 *) 8.2.4 线性势基态能量真值和变分上限。 向 激 发 态 图 径 论 了 基态 的能量 本征 值问题 。下 面，我 们将 讨论如 何求 解径向 激发 态的能量 上一 节 讨 上限。 这一 求解过程 的 详细介绍 可 以参看文 献 [6]。 在这里 我们将求 解 一个线性 势 束缚系统 的基态 和第 一激发 态能 量，以 演示 求解这 类问 题的一 般方 法。首 先， 需要选 取一 组正交 基 （ “试验 ”波函数） 。对于基态，我们再次选择氢原子本征函数 λ 2/3 () r er λ π λ ? =Ψ 0 , ， 而对于第一激发态，我们选用 () () r err λ λ π λ λ ? ?=Ψ 23 3 , 2/3 1 ， (8.2.18) 激 发 态 的波函 数是 有节点 的， 我们选 取的 “试验 ”函数 必须 反映 这一 特性。 第一 激发态 有 一 个节点 ，因 此 “试验 ”函 数必须 有一 个零点 ，如 等式 (8.2.18)所示 。第 一激 发 态 的 “试验 ”波 函数要求与基态本征函数相互正交 0 10 =ΨΨ . 我们可以运用 Mathematica V3.0 语言系统直接检验其正交归一性。 MATHEMATICA In[22]:= psi0[lambda_,r_] := Sqrt[lambda^3/Pi] Exp[-lambda r] (* 定义基态 “实验 ”函数 *) 0 Ψ In[23]:= psi1[lambda_,r_] := Sqrt[lambda^3/(3 Pi)] (3-2 lambda r) Exp[-lambda r] (* 定义第一激发态 “实验 ”波函数 1 Ψ *) In[24]:= 4 Pi Integrate[r^2 psi0[lambda,r]^2, {r,0,Infinity}] (* 检验 的归一性 *) 0 Ψ Out[24]:= Pi Pi 4 1 4 In[25]:= 4 Pi Integrate[r^2 psi1[lambda,r]^2, {r,0,Infinity}] (* 检验 的归一性 *) 1 Ψ Out[25]:= Pi Pi 4 1 4 In[26]:= 4 Pi Integrate[r^2 psi0[lambda,r] psi1[lambda,r], {r,0,Infinity}] (* 检验波函数 0 Ψ 和 的正交性 *) 1 Ψ Out[26]= 4 Pi 0 In[27]:= 4 Pi Integrate[r^2 psi0[lambda,r] a r psi0[lambda,r], {r,0,Infinity}] (* 线性势期望值 0000 ΨΨ= arV *) Out[27]= Pilambda a Pi 8 3 4 In[28]:= 4 Pi Integrate[r^2 psi0[lambda,r] a r psi1[lambda,r], {r,0,Infinity}] (* 线性势矩阵元 101001 ΨΨ== arVV *) Out[28]= ? ? ? ? ? ? ? ? ? Pilambda aSqrt Pi 8 ]3[ 4 In[29]:= 4 Pi Integrate[r^2 psi1[lambda,r] a r psi1[lambda,r], {r,0,Infinity}] (* 线性势期望值 1111 ΨΨ= arV *) Out[29]= Pilambda a Pi 8 5 4 In[30]:=laplacepsi0[lambda_,r_]:=D[psi0[lambda,r],{r,2}]+2/r D[psi0[lambda,r],r] (* 定义基态波函数 0 Ψ 的拉普拉斯量 *) In[31]:= laplacepsi1[lambda_,r_] := D[psi1[lambda,r],{r,2}]+2/r D[psi1[lambda,r],r] (* 定义第一激发态 1 Ψ 的拉普拉斯量 *) In[32]:= 4 Pi Integrate[r^2 psi0[lambda,r](-laplacepsi0[lambda,r]/(2 mu)), {r,0,Infinity}] (* 动能项期望值 0 2 000 2 Ψ ? ?Ψ= μ T *) Out[32]= Pimu lambda Pi 8 2 4 In[33]:= 4 Pi Integrate[r^2 psi0[lambda,r] (-laplacepsi1[lambda,r]/(2 mu)), {r,0,Infinity}] (* 动能项矩阵元 1 2 01001 2 Ψ ? ?Ψ== μ TT *) Out[33]= PimuSqrt lambda Pi ]3[4 2 4 In[34]:= 4 Pi Integrate[r^2 psi1[lambda,r] (-laplacepsi1[lambda,r]/(2 mu)), {r,0,Infinity}] (* 动能项期望值 1 2 111 2 Ψ ? ?Ψ= μ T *) Out[34]= Pimu lambda Pi 24 7 2 4 根据以上 矩 阵元定义 ， 即可求解 方 程 (8.2.7)。 我们将运 用 Mathematica 语言， 分 别 以 解 析 和数值两种方式求能量本征值。 ------------------------------------------------------------------------------ MATHEMATICA ------------------------------------------------------------------------------ (* 能量矩阵元 : *) In[35]:= e00[lambda_,mu_,a_] := lambda^2/(2 mu)+3 a/(2 lambda) In[36]:= e11[lambda_,mu_,a_] := 7 lambda^2/(6 mu)+ 5 a/(2 lambda) In[37]:= e10[lambda_,mu_,a_] := lambda^2/(Sqrt[3] mu)- Sqrt[3] a/(2 lambda) In[38]:= ematrix[lambda_,mu_,a_] := {{e00[lambda,mu,a],e10[lambda,mu,a]}, {e10[lambda,mu,a],e11[lambda,mu,a] }} (* 定义以参数 λ 、 μ 和 为自变量的能量矩阵元a ( )λ ij E *) In[39]:= eeigen[lambda_,mu_,a_] := Eigenvalues[ematrix[lambda,mu,a]] (* 定义以参数 λ 、 μ 和 为自变量的能量本征值a ( )λE *) In[40]:= eeigen[lambda,mu,a] （ * 解析推导本征值 *） Out[40]= {(5*lambda^4*mu + 12*a*lambda*mu^ - 2*lambda*mu*Sqrt[4*lambda^6 – 6*a*lambda^3*mu + 9*a^2*mu^2])/(6*lambda^2*mu^2), (5*lambda^4*mu + 12*a*lambda*mu^2 + 2*lambda*mu*Sqrt[4*lambda^6 - 6*a*lambda^*mu + 9*a^2*mu^2])/(6*lambda^2*mu^2) } In[41]:= FullSimplify[%] (* 化简上式 *) Out[41]= {(5*lambda^3 + 12*a*mu - 2*Sqrt[4*lambda^6 - 6*a*lambda^3*mu + 9*a^2*mu^2])/(6*lambda*mu), (5*lambda^3 + 2*(6*a*mu + Sqrt[4*lambda^6 - 6*a*lambda^3*mu + 9*a^2*mu^2]))/(6*lambda*mu)} In[42]:= eeigen[1,1/2,1] (* 在 λ=1 (GeV)、 μ =0.5 (GeV)及 a 条件下，数值求解本征值 2 )(1 GeV= E *) Out[42]= {(11 - Sqrt[13])/3, (11 + Sqrt[13])/3} In[43]:= N[%] (* 数值化 Mathematica 指令 N[expr]可求表达式 expr 的数值 *） Out[43]= {2.46482, 4.86852} ------------------------------------------------------------------------------ - 由上式可知，能量本征值的变分上限分别为 GeV2.46482)1( =S upper E , E GeVS 86852.4)2( upper = , 其相应的真值解为（由 Airy 函数的头两个零点给出） GeV2.33811)1( =SE true GeV4.08795)2( =SE true 对于基态（以 1S 表示） ，我们所得的变分上限值 的相对误差为 upper E %4.5 )1( )1()1( = ? SE SESE true trueupper 这一结果 比 前面给出 的 6% 相 对误差 稍好一些。 由此可以 看 出 ， 通过增 加矩阵的 大 小 （此处 是由 1 变成 2 ） ， 可以提高 上 限值的近 似 程度。对 于 第一激发 态 （ 以 2S 表示 ） ，变 分 上限值 的相对误差为 1× E 2× upper %1.19 )2( )2()2( = ? SE SESE true trueupper 这一误差 值 稍大了一 些 ，我们可 以 运用方程 (8.2.8)及 (8.2.9)，来 改进 这一结果 。 这 样 我 们就能找到覆盖在所选的这组 “试验 ”波函数上的能量本征值最小值。 但是， 改变 变分参数 λ 往往 会 破 坏来自 于 不 同激发 态 “试验 ”波函数的 正 交性。 原则 上， 特征方 程并 非由方 程 (8.2.7)给出 ， 关于这 一困 难的详 细讨 论参见 文献 [6]。 考虑到 这一因 素，下面我们将仅限于基态 1S 的讨论，因为基态不涉及正交性破坏的问题。 MATHEMATICA In[44]:= e00eigen[lambda_] := Part[Eigenvalues[ematrix[lambda,1/2,1]],1] (* 在 及 条件下，定义以 GeV5.0=μ 2 1 GeVa = λ 为自变量的基态本征值 *) In[45]:= e00eigen[lambda] (* 解析计算基态能量本征值 *) Out[45]= (6*lambda + 5*lambda^4 - lambda*Sqrt[9 - 12*lambda^3 + 16*lambda^6])/(3*lambda^2) In[46]:= FindMinimum[%,{lambda,0.5}] (* 寻找 1S-态能量的最小值（以 GeV5.0=λ 为起始点） *) (* 虽然我们希望能够用 Mathematica 指令 FindMinimum[e00eigen[lambda],{lambda,0.5}]求 得正确的最 小值，但实 际运行显示 Mathematica 无法同步 实 现计算本 征 值 、 抽取矩 阵元相关 部 分并求解 最 小值。 因此 ， 在 实 际 运用中，我们是将计算过程仔细地划分成几部分来进行的。 *) Out[46]= {2.4322, {lambda -> 0.665633}} 通过对变分参数最小化后，得到改进的新结果为 GeV2.43220)1( var =SE 此时参数为 GeV665633.0 min =λ . 显然，我们已经成功地减小了相对误差。 %9.3 )1( )1()1( var = ? SE SESE true true . 正 如 前 面所说 ， 最 小化过 程一般将 导 致基态和 激 发态具有 不 同的变分 参 数值： ji λλ ≠ 。 由此 ， 所有这些态一般将不再正交。即 ( ) () ijjjii δλλ ≠ΨΨ . 此时，本征值问题的特征方程 (8.2.7)将变为 () ( ) ( ) ( )[ ] 0 ? =ΨΨ?ΨΨ jjii upper jjii EH λλλλdet . 在下面的 内 容中，我 们 将讨论如 何 通过扩大 矩 阵大小， 来 进一步提 高 薛定格能 级上 限 的精确度。 5. 勒盖尔 (Laguerre)上限 从前 面 的讨论 可以 看出， 精确 地求解 能量 上限的 关键 步骤是 选择 一组合 适的 “试探 ”波 函数。 我们 这里引入 一 组勒盖尔 多 项 式 ( ) ，以 提 高 “试探 ”波 函数的性 能 。 并进 一步 引入两 个变分参 数 ： 含质量 量 纲的参数 γ k L λ 及无 量 纲 的参数 β ； 以 及 选 择可以得 到 含角动量 及其 投l 影 的 “试探 ”波函数m ()x lmk G , ψ ，其相应的坐标空间表述为 () () () ( ) ( )??= ?+ lmk l lmk YxLxxNx GGGG λλψ γ β 2exp 1 , . (8.2.19) 波函数的 归 一化条件 要 求变分参 数 λ 数值为正 ： 即 0>λ 。 这里 ( )? lm Y 标识在 与立体角 内 的角动量为 ，其投影为 的球谐函数，其正交关系约定为 ? l m . (8.2.20) () () mmllmllm YYd ′′′′ =??? ∫ δδ * (8.2.19)式 所定义的 波 函数正交 归 一，不但 确 定了归一 化 常 数 ，还对 参数N γ 的取值作 出 限制： βγ 22 += l 。由此可得 () () () () ()()?? +++Γ = + ?+ ++ lm l k l l lmk YxLxx kl k x GGGG λλ β λ ψ β β β 2exp 122 !2 22 1 122 , . 它满足的正交归一化条件为 () () mmllkkmlklmk xxxd ′′′′′′ = ∫ δδδψψ GG , * , 3 . 很显然， 正 交归一化 同 时也对第 二 个变分参 数 β 作出限制 2 1?>β ，即 其 取 值域为 2 1 ?>β 。 为便于讨 论 ，我们做 以 下化简： 质 量标度 GeVm 1 1 m 2 == ， 变 分参数 GeV1=λ 及 1=β 。 为表述求 解 的一般过 程 以及便于 比 较 ， 我们 仍将 讨论线性 势 的情况 V ar= ， 并取 。 这里的计算 至少在矩阵 大小扩大到 4 2 1GeVa = 4× 阶时，所有 计算仍能够 解析求解（ 手工推导） 。这 一工作留给读者作为一个练习。下面将演示运用 Mathematica 语言进行计算的过程： ? 定义所选取的 “试探 ”波函数 ( )x lmk G , ψ ； ? 计算拉普拉斯算符（动能项）的矩阵元； ? 计算径向坐标 r （势能项）的矩阵元； ? 确定所得的总能量矩阵 ()λ ij E 的本征值 ( )λE ； ? 比较不同矩阵大小对应的能量本征值，以期观测到收敛行为。 ------------------------------------------------------------------------------ - MATHEMATICA ------------------------------------------------------------------------------ - In[47]:= psix[k_,l_,m_,r_] := Sqrt[2^(2 l+3) k!/Gamma[2 l+3+k]] r^l Exp[-r]*LaguerreL[k,2 l+2,2 r]*SphericalHarmonicY[l,m,theta,phi] (* 定义 “试探 ”波函数 ( )x lmk G , ψ *) In[48]:= psi[k_,r_] := psix[k,0,0,r] （ * 由于我们仅讨论 S-波情况，可使 “ 试探 ”波函数有 0== ml *） In[49]:= delta[k_,r_] := D[psi[k,r],{r,2}]+2/r D[psi[k,r],{r,1}] （ * 定义拉普拉斯算符作用在 0=l （ S-波）态上的 ( )r k ψ? *） In[50]:= intks[k_,s_,r_] := psi[s,r] delta[k,r] （ * 定义积分 () () rr ks ψψ ? *） In[51]:= kinen[k_,s_] := -4 Pi Integrate[r^2 intks[k,s,r],{r,0,Infinity}] （ * 动能算符 ??=T 的矩阵元 （注：由于已取两粒子质量 () ()(rrrdr ks ψψ ?? ∫ ∞ 2 0 为 ，约化质量 GeVmm 1 21 == μ 也将等于 1 ） *） GeV2/ In[52]:= poten[k_,s_] := 4 Pi Integrate[r^3 psi[s,r] psi[k,r],{r,0,Infinity}] (* 势能算符 ( ) rr =V 的矩阵元 *) () ()rrrrdr ks ψψ 2 0 ∫ ∞ In[53]:= toten[k_,s_] := kinen[k,s]+poten[k,s] (* 总能量矩阵元 *) （ * 下面我 们将用 Table 指 令来构 造矩阵。 因 为矩阵指 标 是 从 0 开始 计数的， 我 们需要重 新定义矩阵，以使 时给出 1 矩阵等等。 *） 1=x 1× In[54]:= totenmat[x_] := Table[toten[k,s],{k,0,x-1},{s,0,x-1}] （ * 借助指 令 Eigenvalues[M]， 定 义函数 eeigen[x]。 这 一 指令将给 出 xx× 阶矩阵 M 的 本征值，亦即对任意的 xx× 阶矩阵对角化。 *） In[55]:= eeigen[x_] := Eigenvalues[totenmat[x]] In[56]:= eeigen[1] （ * 1 能量矩阵的本征值 *） 1× Out[56]= ? ? ? ? ? 2 5? In[57]:= eeigen[2] （ * 2 能量矩阵的本征值 *） 2× Out[57]= ()() ? ? ? ? +? 1311 3 1 ,1311 3 1 ? ? In[58]:= N[%] （ * 运用 N[%] 指令对上式输出进行数值化处理 *） Out[58]= {2.46482,4.86852} （ * 数值化的基态和第一激发态能量 *） In[59]:= eeigen[3] （ * 3 能量矩阵的本征值 *） 3× Out[59]= ? ? ? ? ? +? 75,75, 2 9? In[60]:= N[%] （ * 运用 N[%] 指令对上式输出进行近似数值计算 *） Out[60]= {2.35425,4.5,7.64575} （ * 数值化的基态和第一、第二激发态能量 *） （ * 通常，一个 5 矩阵的本征值是无法解析求解的。我们来进行数值求解 *） 5× In[61]:= N[eeigen[5]] （ * 求基态和前四个激发态的能量的近似数值 *） Out[61]= {2.34136,4.13334,5.72535,8.11424,15.519} In[62]:= N[eeigen[10]] （ * 注： 10 10× 能量矩阵本征值的计算将耗费一定机时 *） Out[62]= {2.33812,4.08858,5.53209,6.83859,8.14892,9.91409,12.195,14.096,17.146, 49.7026} （ * 基态和前九个径向激发态的能量 *） 在表 8.2.1 中，我们 对 基态、激 发 态能级上 限 精确度与 能 量矩阵大 小 的关系进 行了 一 个对照比较。 表 8.2.1 不同大小矩阵对应的能量本征值相对误差的比较 矩阵大小 1S 态 2S 态 1 1× %6 -- 22× %5 19 % 33× 0 %7. 10 % 55× 0 %1. 1 % 10 10× %104 4? × 2 %10 2? × 可见，这 样 的处理对 精 度的提高 极 为显著。 需 要指出的 是 ，这种处 理 方法应用 范围 很 广， 对 不 同 于幂指数 V 的 势函数 n r= ( )V ， 以 及不同于 非 相对论r mp 2/ 2 G 的微 分算符都 适 用 。 因此， 我们 可以用 这一 方法来 求解 更为复 杂的 哈密顿 算符 的能量 本征 值问题 ，如 处理包 含 平方根相 对 论动能算 符 22 mp + G 的 准 相 对论过 程 [10]。 另 一方面， 由 于 直 到 4 阶能量矩 阵 4× E 是能够 解 析对角化 的 ， 这一特 性 使我们能 够 对纯数值 计 算给出的 结 果进行控 制 。 但 是 ， 我们需要 牢 记此处计 算 得出的数 值 结果， 仅是 表示能量 本 征值真实 解 的上限： ()λE≤E 。 true 总结 ： 在本节中 ， 我们讨论 了 如何运用 Mathematica 语言和一 些 基本物理 原 理简便 地 处理束 缚态 问题， 并且 由此计 算给 出的数 值结 果具有 较好 的精度 。我 们在计 算中 并未真 正 涉及相 关波 函数的 确定 问题。 然而 ，我们 的这 种处理 和计 算方法 说明 [10]： 只要矩 阵尺 度 足够大 ，即 可实现 对于 波函数 任意 精度的 逼近 。另一 方面 ，对于 小尺 度矩阵 而言 ，一般 来 说即使能级计算给出的结果相当准确，我们也完全不能相信所选取的波函数。 作为一个例 子，下面将 绘制 0=k 和 5=k 情况 下 “试探 ”波函数的三维 图形。以便 直 观了解一下这些波函数到底是什么样的。 ------------------------------------------------------------------------------ - MATHEMATICA ------------------------------------------------------------------------------ - In[63]:= psix[k_,l_,m_,x_,y_] := Sqrt[2^(2 l+3) k!/Gamma[2 l+3+k]] Sqrt[x^2 + y^2]^l Exp[-Sqrt[x^2+y^2]] LaguerreL[k,2 l+2,2 Sqrt[x^2+y^2]] （ * 将 “试探 ”波函数 ( )x lmk G , ψ 化为自变量 x 、 的图形函数 *） y (yx lmk , , ψ ) In[64]:= Plot3D[psix[0,0,0,x,y],{x,-4,4},{y,-4,4},TextStyle -> {FontSlant-> ”Italic”,FontSize->12}] （ * 绘制 0=== mlk 情况下 “试探 ”函数 ( )yx lmk , , ψ 的三维图形。 这时的波函数本征基态 *） 图 8.2.5 0=== mlk 时的 “试探 ”波 函 数 ( )yx lmk , , ψ 。 （ * 注： 对 ，函 数0≠k ()yx lmk , , ψ 并 不是对应 于 特定的激 发 能级上的 试 探波函数 。 合 适 的 求 阵 的本 征 *） In[65]:= Plot3D[psix[5,0,0,x,y],{x,-4,4},{y,-4,4}, TextStyle-> （ * 绘制 且 试探函 数要 通 过 解 能 量 矩 矢 量 来决 定，试 探波 函数则 是各 种试探 函数 的叠加 。 {FontSlant->”Italic”,FontSize->12}] 5=k 0== ml 时的 ψ ( )yx lmk , , “试探 ”函数 的三维图形。 *） 图 8.2.6 k 且5= 0== ml 时的“试探”波函数 ( )yx lmk , , ψ 。