第八章 Mathematica在量子力学中的应用举例 §8.1粒子在中心力场中的运动问题 设电子与原子核的约化质量为=mM,m(m)=-2e,哈密顿量为 +(r) +F(r) (8.1.1) 其中r为粒子所处的空间位置到中心势原点的距离。利用中心势的球对称性,我们将薛定格 方程写为在球坐标中的表示 (,9.o)=(E-().o 2 ar/(sins af sine ae/sin'e ag (8.1.2) 角动量算符的定义为:L=x×p。可以证明[L,H=0,所以角动量L是守恒量,即在中心 力场中运动粒子的一个重要特征是角动量守恒。由此可以得到L2(角动量的平方)也是守 恒量在求解中心力场作用下粒子的能量本征方程时,(B,2,构成对易算符的一个完全 集 (8.1.3) 其中在球坐标中的角动量平方算符可以表示为: n080(a0)sin20 a9 薛定格方程(8.1.2)则可以写为 2urorlor (8.1.5) 波函数v/(,)与极角(-x/2≤0≤丌/2)和方位角(0≤q≤x)的关联是由算符D和 L决定的。假定满足薛定格方程的本征波函数v(,O,)可以分离变量表示为 w(r, 0, p)=R()y(0, p)=r(Do(gp(o) (8.1.6)
第八章 Mathematica 在量子力学中的应用举例 §8.1 粒子在中心力场中的运动问题 设电子与原子核的约化质量为 m M m M e e + µ = , r Ze r 2 V( ) = − ，哈密顿量为 ( ) 2 (ˆ) 2 ˆ ˆ 2 2 2 2 V r V r p + ∇ = + = − µ µ = G = H . (8.1.1) 其中 r 为粒子所处的空间位置到中心势原点的距离。利用中心势的球对称性，我们将薛定格 方程写为在球坐标中的表示 ( ) = ∂ ∂ + ∂ ∂ ∂ ∂ + ∂ ∂ ∂ ∂ − ψ ϑ ϕ θ θ ϕ θ µ θ θ , , sin 1 sin sin 1 2 2 2 2 2 2 2 r r r r r = (E −V(r)) ( ψ r,ϑ,ϕ). (8.1.2) 角动量算符的定义为：L x p G G G = × 。可以证明[ ，所以角动量 是守恒量，即在中心 力场中运动粒子的一个重要特征是角动量守恒。由此可以得到 ] 0 ˆ ,ˆ L H = Lˆ 2 Lˆ （角动量的平方）也是守 恒量。在求解中心力场作用下粒子的能量本征方程时，( ) Lz H L ˆ , ˆ , ˆ 2 构成对易算符的一个完全 集。 − ∂ ∂ ∂ ∂ ≡ ∇ = 2 2 2 2 2 ˆ 1 = L r r r r ∆ . (8.1.3) 其中在球坐标中的角动量平方算符可以表示为： ∂ ∂ + ∂ ∂ ∂ ∂ = − 2 2 2 2 2 sin 1 sin sin 1 ˆ θ θ ϕ θ θ θ L = . (8.1.4) 薛定格方程(8.1.2)则可以写为 ψ( ) θ ϕ ( )ψ( θ ϕ µ , , , , ˆ 2 2 2 2 2 2 r E V r L r r r r = − − ∂ ∂ ∂ ∂ = = ) ) − . (8.1.5) 波函数ψ(r,θ,ϕ 与极角θ（−π / 2 ≤ θ ≤ π / 2 ）和方位角ϕ (0 ≤ ϕ ≤ π ) 的关联是由算符 和 决定的。假定满足薛定格方程的本征波函数 2 Lˆ Lz ˆ ψ(r,θ,ϕ)可以分离变量表示为 ψ( ) r,θ,ϕ ≡ R(r) ( Y θ,ϕ) ≡ R(r)Θ(ϑ)Φ(ϕ). (8.1.6)
L.在球坐标系中可以表示为:L=一i0,该算符的本征值由求解本征方程 i。Φ()=L(p), (8.1.7) 来得到。方程(8.1.7)的解为 Φ(q)=Ael2/h (8.18) 由于(8.1.8)式所示波函数解必须唯一确定,因而它也必定满足条件:Φ()=(2n+q), 并且角动量算符L的本征值应当是离散的,其本征值表示为:L:=mh,(m=0,±1,±2,) 由本征波函数的归一化条件,方程(8.1.7)归一化的解可以写为 d(p)= (8.1.9 类似地,对另一个守恒量角动量平方,我们有本征方程: i(.o)=-h2 H(O,q)=L2y(0,) (8.1.10) sinb0日 a0) sin ag 方程(8.110)的解是球谐函数Ym如果本征值满足L2=(+1)h2,方程(8.1.10)写为 +1(+1)}() (8.1.11) sin 80 角动量算符一作用在球谐函数m上的本征值由角量子数l=012决定。对应于确定的 角量子数/,算符L2的本征值则为l(+1)h2,此时磁量子数m则描写该角动量在二轴上的 投影,它的取值范围为:m=0,±1,±2,,+。这就是说:对确定的角动量量子数l,应当 有2+1个本征函数Y。对磁量子数m为正时的情况,球谐函数的完整表达式为 Yn(O,q)=(-) +m)4r P(cose le imp (8.1.12) 其中P(x)为阶的第m个伴随勒让德函数。如果磁量子数为负时(-m),其球谐函数满 足如下关系式
Lz ˆ 在球坐标系中可以表示为： ∂ϕ ∂ Lz = −i= ˆ . 该算符的本征值由求解本征方程 ( ) ϕ (ϕ) ϕ Φ = Φ ∂ ∂ − Lz i= ， (8.1.7) 来得到。方程（8.1.7）的解为 ( ) ϕ / = ϕ z iL Φ = Ae . (8.1.8) 由于（8.1.8）式所示波函数解必须唯一确定，因而它也必定满足条件：Φ( ) ϕ = Φ(2π +ϕ)， 并且角动量算符 的本征值应当是离散的，其本征值表示为：L Lz ˆ z = m= ，( . 由本征波函数的归一化条件，方程（8.1.7）归一化的解可以写为 m = 0,±1,±2,...) ( ) ϕ π ϕ im e 2 1 Φ = . (8.1.9) 类似地，对另一个守恒量-角动量平方，我们有本征方程： ( ) ( ) θ ϕ (θ ϕ θ θ ϕ θ θ θ θ ϕ , , sin 1 sin sin 1 , ˆ 2 2 2 2 2 2 L Y Y = L Y ∂ ∂ + ∂ ∂ ∂ ∂ = −= ). (8.1.10) 方程(8.1.10)的解是球谐函数Y 。如果本征值满足 L ，方程(8.1.10)写为 l,m 2 2 = l(l +1)= ( 1) ( ) , 0 sin 1 sin sin 1 2 , 2 2 = + + ∂ ∂ + ∂ ∂ ∂ ∂ θ ϕ θ θ ϕ θ θ θ Yl m l l . (8.1.11) 角动量算符 Lˆ2 作用在球谐函数Yl,m 上的本征值由角量子数l = 0,1,2,... m 决定。对应于确定的 角量子数l ，算符 的本征值则为l ，此时磁量子数 则描写该角动量在 轴上的 投影，它的取值范围为： 2 Lˆ 2 (l +1)= l z m = 0,±1,±2,...,± m 。这就是说：对确定的角动量量子数 ，应当 有 个本征函数Y 。对磁量子数 为正时的情况，球谐函数的完整表达式为 l 2l +1 l,m ( ) ( ) ( ) ( ) ( ) ( ) ϕ θ π θ ϕ m im l m l m P e l l m l m cos 4 2 1 ! ! , , 1 + + − Y = − . (8.1.12) 其中 P (x)为l 阶的第 个伴随勒让德函数。如果磁量子数为负时（ m l m − m ），其球谐函数满 足如下关系式
(,)=(-y 20) (8.1.13) 显然,球谐函数m也是算符L的本征函数。容易证明类似(⑧.1.7)式,球谱函数Y满足: (8.1.14) 因而球谐函数Y既是角动量算符平方L2,也是角动量算符的二分量L的本征函数。在 Mathematica中球谐函数表示为 SphericalHarmonicY]。勒让德多项式表示为 Legendre[] 将(8.1.6)式代入薛定格方程(8.1.2),再应用上面推导出的角动量部分波函数所满 足的薛定格方程,可以得到本征波函数v(r,,)表示中的径向部分R()应当满足的方程 dr 2 dR (8.1.15) 下面我们以氢原子为例进行分析。定义波尔半径a=≈529×10-1m为长度单位,即 p=r/a;以氢原子的电离能量E0=2a0h 4-≈13.5e为能量单位,即E=E/E0;定 义径向函数R(p)=l(p)/p。这时方程(8.1.15)写为 dap)++22 (p)=0 (8.1.16) 能量E的值是由方程(8.1.16)的本征值和本征函数决定的。 我们考虑稳定状态(束缚态),即E∞时的极限行 为,发现由波函数的幺正性条件要求上述两种表达方式下都可以推出 u(p)=pef,(p). (8.1.17 将(8.1.17)式代入(81.16)后,求解得到超几何函数(F)形式的解。 ()=c1F{1+1-2,1+22m (8.1.18)
( ) ( ) ( ) ( ) θ ϕ (θ ,ϕ ! ! , 1 * , l, m m l m Y l m l m + − = − − Y ) . (8.1.13) 显然，球谐函数Y 也是算符 的本征函数。容易证明类似(8.1.7)式，球谐函数Y 满足： l,m Lz ˆ l,m z l m Yl m m Yl m L Y i , , , ˆ = = = ∂ ∂ ≡ − ϕ . (8.1.14) 因而球谐函数Y 既是角动量算符平方 ，也是角动量算符的 分量 的本征函数。在 Mathematica 中 球谐函 数 表 示 为 SphericalHarmonicY[] 。 勒 让 德 多 项 式 表 示 为 LegendreP[]。 l,m 2 Lˆ z Lz ˆ 将（8.1.6）式代入薛定格方程(8.1.2), 再应用上面推导出的角动量部分波函数所满 足的薛定格方程，可以得到本征波函数ψ(r,θ ,ϕ)表示中的径向部分 R(r)应当满足的方程。 0 2 2 ( 1) 2 2 2 2 2 = + − + + + R r l l r Ze E dr dR dr r d R = µ . (8.1.15) 下面我们以氢原子为例进行分析。定义波尔半径 m m e a e 11 2 2 0 5.29 10− = ≈ × = 为长度单位，即 0 ρ = r / a ; 以氢原子的电离能量 eV m e a e E e 13.5 2 4 4 0 2 0 = = ≈ = 为能量单位，即 = E E0 ε ； 定 义径向函数 R(ρ) = u(ρ)/ ρ 。这时方程(8.1.15)写为 ( ) 0 ( ) 2 ( 1) 2 2 2 = + + + − ρ ρ ρ ε ρ ρ u Z l l d d u . (8.1.16) 能量ε 的值是由方程(8.1.16)的本征值和本征函数决定的。 我们考虑稳定状态（束缚态），即ε < 0的状态。分析表明函数u(ρ) 可以表示为多项式 或者指数形式。为了找出u(ρ) 的近似式，我们通过考察它在 r → 0 和 r → ∞ 时的极限行 为，发现由波函数的幺正性条件要求上述两种表达方式下都可以推出 u( ) ( ) . (8.1.17) 1 ρ ρ ρ γρ l l e f + − = 将(8.1.17)式代入(8.1.16)后，求解得到超几何函数( ) 1F1 形式的解。 ( ) = + − + γρ γ ρ 1 1 1 ,2l 2;2 Z f c F l l . (8.1.18)
其中y≡√-E。现在我们由式(81.17)得到电子在库仑势中的波函数的径向部分为 R(p)=MM/pc1F|+1-m2+23 (8.2.19) 由于归一化条件的要求,(8.1.18)的级数表示必须只有有限项。这个限制就给出了能量的 值 (n2=012 (8.1.20) 由此我们得到 (8.1.21) 由y和E的定义,则 E=E2=-E2. (8.1.22) (n1 其中n为主量子数(m=12…)。它是由径向量子数n,(n1=012,)和轨道角动量量子数 (=0.2,)决定的。在这里我们引入一组称为“拉盖尔( Laguerre,)多项式”的特殊正交多 项式[2],拉盖尔多项式式由级数定义为 e(x)=∑(-) 相应的归一化为 jdxr'expe-x()(t、 kl 超几何函数与拉盖尔多项式间有如下关系式 L(x)= nr(+a) na+1. x 这样电子在库仑势中的波函数的径向部分的解也可以写为 R(P=N /eln(2+1)(22 径向部分波函数(8.1.23)中的拉盖尔〔 Laguerre)多项式的性质见文献[2]。相应的波函数为 Wmlm(e, e, )=Nupe Z//7+1-n,21+2 (0,q) (8.1.24)
其中γ ≡ − ε 。现在我们由式(8.1.17)得到电子在库仑势中的波函数的径向部分为 = + − + − ρ ρ ρ ρ n Z R N e F l n l l Z n n l 2 ( ) 1 ,2 2; 1 1 / , . (8.2.19) 由于归一化条件的要求，(8.1.18)的级数表示必须只有有限项。这个限制就给出了能量的 值 = − + − γ Z n l r 1 ， ( = 0,1,2,...). (8.1.20) nr 由此我们得到 + +1 = n l Z r γ . (8.1.21) 由γ 和ε 的定义，则 ( ) 2 2 0 2 2 0 1 n E Z n l E Z E r = − + + = − . (8.1.22) 其中 为主量子数 。它是由径向量子数 n 和轨道角动量量子数 决定的。在这里我们引入一组称为“拉盖尔(Laguerre)多项式”的特殊正交多 项式 [2]，拉盖尔多项式式由级数定义为 n = 0, (γ Lk (n = 1,2,...) r ( = 0,1,2,...) r n l (l 1,2,...) ) ( ) ( ) ( ) ! 1 0 j x k j k L x k j j j k − + = ∑ − = γ γ . 相应的归一化为 ( ) ( ) ( ) ( ) ( ) ( ) k k kk k k dxx x L x L ′ x ′ ∞ Γ + + − = δ γ γ γ γ ! 1 exp 0 ∫ . 超几何函数与拉盖尔多项式间有如下关系式 ( ) ( ) F ( ) n x n n L x n , 1; ! 1 1 ( ) 1 1 ( ) − + Γ + Γ + + = α α α α . 这样电子在库仑势中的波函数的径向部分的解也可以写为 = ′ + + − ρ ρ ρ ρ n Z R N e L l n l l Z n n l 2 ( ) / (2 1) , . （8.1.23） 径向部分波函数(8.1.23)中的拉盖尔(Laguerre)多项式的性质见文献[2]。相应的波函数为 ψ ( ) ρ θ ϕ ρ ρ (θ ϕ ρ , 2 , , 1 ,2 2; 1 1 , / , , , l m l Z n n l m n l Y n Z N e F l n l = + − + − ) . (8.1.24)
在(8.1.19)和(8.1.24)式中的归一化常数为 yn-(21+1)V2n(n-1-1) (8.1.24) Mathematica Package file Coulomb. BeginPackage["CoulombPotential"] Clear [WaveF, WaveR, WaveA] WaveF:; usage=" WaveF[Z,r, theta,phi,n,1,mJ计算电子在库仑 势中本征波函数的表示。Z为原子核的电荷数.r为电子到中心势原点的距离 theta和phi为球坐标中的角度.n,1和m为能量和角动量算符的量子数 WayeR: usage=" Waver[Z,r,n,1]计算电子在库仑势中的本征波函数 径向部分的表示。Z为原子核的电荷数.r为电子到中心势原点的距离.n和1 为能量和角动量算符的量子数。 WaveA: usage=" WaveA[ theta,phi_,1,m]计算电子在库仑势中本征波函 数的角度关联部分表示。 theta和phi为球坐标中的角度.1和m表示角动 量算符的量子数。 定义公共变量 rnlm theta: usage phi: : usage Begin["’ Private’"] (*一产生库仑势中波函数的径向部分一* WaveR[Z,r,n, 1] Module[ lunit, tmpl 归一化常数 unit=(Sqrt[血n+1)!/(2n(n-1-1)!)]((2Z/n)^(1+3/2))/(2 1+1)! (*一产生波函数径向部分的定义一*) tmp unit r"1 Exp[-((z r)/n)] HypergeometriclF1[1+1-n, 21+2 (2zr)/n]
在(8.1.19)和(8.1.24)式中的归一化常数为 ( ) ( ) ( ) 3/ 2 , 2 2 1 ! ! 2 1 ! 1 + − − + + = l n l n Z n n l n l 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]
(*一产生库仑势中本征波函数的角度相关部分一* WaveALtheta, phi, 1 m]: Modulelitr tmp SphericalHarmonicY[l, m, theta, phi] (*一产生电子在库仑势中的本征波函数一*) WaveF[Z, r, theta_, phi,n, 1,m]: Modulelitmp tmp=WaveR[Z, r, n, 1] WaveA[theta, phi, 1 Endl EndPackage [ 当我们需要对电子在原子核的库仑势中的本征波函数习性进行分析时,我们 可以首先调入程序包 Coulomb. m,然后调用程序包中定义的函数。例如通过运 行下面的指令: r,u, Prolog->Thickness [o 0011] Plot [Abs [WaveAltheta, Pi/2, 2, 11]2 Itheta, 0, Pi], AxesLabel1->"theta,"y", Prolog->Thickness [0. 0011] Plot3D [Abs [WaveF[1, r, theta, Pi/2, 3, 2, 2]]2, r, 0, 15, theta, 0, Pi], Ligh ting->True] 0.025 -0.075 -0.1
] (* --- 产生库仑势中本征波函数的角度相关部分 --- *) 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，然后调用程序包中定义的函数。例如通过运 行下面的指令： (*----------------------------------------------------------------*) "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] (*---------------------------------------------------------------*)
0.14 0.08 0.06 0.04 0.02 th吧ta 0.5 1.5 图(8.1.1)本征波函数径向部分的四条曲线(当z=1,1=0和n=1,2,3,4时)。 图(8.1.2)本征波函数角度关联部分绝对值平方随极角变化的曲线(当=x/2,=2和 m=1时)。 0.0001 图(813)本征波函数绝对值平方随r和极角O变化的三维曲线(当=x/2,1=2和m=1 时)
图（8.1.1）本征波函数径向部分的四条曲线（当Z = 1，l = 0和n = 1,2,3,4 时）。 图（8.1.2）本征波函数角度关联部分绝对值平方随极角θ 变化的曲线（当ϕ = π 2，l 和 1.3）本 = 2 m = 1时）。 图（8. 征波函数绝对值平方随 r 和极角θ 变化的三维曲线（当ϕ = π 2， 和 l = 2 m = 1 时）
8.2求非相对论性薛定格方程本征能量限 用 Mathematica V4.0系统的指令,对应的计算过程可表述为 MATHEMATICA V4. 0 (*积分*) In[1]:=E-r"dr Out[1]= /[Reln]>-1&& Re[a]>0, t-"Gamma l +n]. erdr] 这一输出结果的含义是:如果Re4]>0,且Ren>-1,则以上积分的结果为 A-"T(1+n),否则将输出 dr 这意味着 Mathematica无法求解该问题。由此可以得归一化因子 归一化的“试验”波函数为 23 为保险起见,我们可以检验一下关系式两边的量纲。根据以前的讨论,我们知道关系式左 边的量纲为 DimE/2]。为使指数运算exp{-r}有意义,乘积必须是无量纲的量,即 Dm=1,由此有oa-Dme,甲 Dm平]=DmE32]=Dm312]。 很显然,在以上推导中至少量纲是正确的。下面我们演示一下如何运用 Mathematica语言 作以上定义和计算 采用 Mathematica v4.0的对应计算为:
8.2 求非相对论性薛定格方程本征能量限 用 Mathematica V4.0 系统的指令，对应的计算过程可表述为： MATHEMATICA V4.0 （* 积分 *） In[1]:= E ∫ r dr r n ∞ − 0 λ Out[1]= [Re[ ] 1&& Re[ ] 0, [1 ], ] 0 1 If n Gamma n e r dr n r n ∫ ∞ − − − > − > + λ λ λ 这一输出结果的含义是：如果 Re[λ] > 0 , 且 Re[n] > −1 则 ( ) 1 1 n λ n − − Γ + ， 以上积分的结果为 ，否则将输出 dr r r n ∫ ∞ 0 exp{λ } , 这意味着Mathematica无法求解该问题。由此可以得归一化因子 N = ， π λ 3/ 2 归一化的“试验”波函数为 Ψ(r, ) r e λ π λ λ − = 3/ 2 为保险起见，我们可以检验一 . 下关系式两边的量纲。根据以前的讨论，我们知道关系式左 边的量纲为 Dim[E3/ 2 ]。为使指数运算 exp{−λr}有意义，乘积 λr 必须是无量纲的量，即 Dim[λr] = 1。由此有 im[E]，即 [ ] [ ] D Dim r Dim λ = = 1 。 很显然，在以上推导中至少量纲是正确的。下面我们演示一下如何运用 Mathematica 语言 作以上定义和计算。 [ ] [ ] [ ] 3 / 2 3 / 2 Dim Ψ = Dim E = Dim λ 采用 Mathematica V4.0 的对应计算为：
MATHEMATICA V4. 0 In[2]:=4xr2Ed(*积分*) out[2]=4丌∥Re[4]>0, fe-2ir2dr] 输出的含义是:当Re>0时,计算结果为x/23,否则, Mathematica无法求解,将返回输 入形式4zc-2br2dh。完整的结果应是 4=1 下一步,我们将借助引入的“试验”波函数求动能项的期望值。由于我们只讨论基态的能 量本征值,而对基态量子数l=0,此时在径向中心力场势情况下可采用拉普拉斯算子形式 为 △≡V2d2,2d 其期望值为 ∫dx,,2)=Jdx (d22d ∫(2x-2 4x26)2、m(2 (2)(2x)2 =-2 我们可以看到这里的量纲检验仍然是正确的。我们在下式中省略了Dim[…]符号: y平→E-E32E2E32=E2→>2 动能项的期望值为 ∫d2(-2 8.2 相应的 Mathematica v4.0计算过程为 MATHEMATICA V4. 0 In【3]:=w-x1=PE(*定义“试验”波函数*)
MATHEMATICA V4.0 In[2]:= r E dr （* 积分 *） r ∫ ∞ − 0 2 2 4 λ π Out[2]= , ] 4 1 4 [Re[ ] 0, 0 2 2 3 If e r dr r ∫ ∞ − > λ λ π λ 输出的含义是：当Reλ > 0 时，计算结果为 3 π λ ,否则，Mathematica 无法求解，将返回输 入形式 dr 。完整的结果应是： 2 e r r 0 2 4 ∫ ∞ − λ π 1 4 4 3 2 = λ π N ⇒ π λ 3 N = . 下一步，我们将借助引入的“试验”波函数求动能项的期望值。由于我们只讨论基态的能 量本征值，而对基态量子数l ，此时在径向中心力场势情况下可采用拉普拉斯算子形式 为 = 0 dr d dr r d 2 2 2 2 ∆ ≡ ∇ = + . 其期望值为 ( ) ( ) r r e dr d dr r d d x r r d xe λ λ π λ λ λ − − Ψ ∆Ψ = + ∫ 2 , , 2 2 3 3 3 * ∫ ( ) r dr r r e λ π λ λ π λ 2 0 2 2 3 4 2 − ∞ ∫ = − ( ) ( ) ( ) ( ) Γ − Γ = 3 2 3 2 2 2 2 2 3 4 λ λ λ λ λ 3 2 4 1 4 λ λ λ = − = − . 我们可以看到这里的量纲检验仍然是正确的。我们在下式中省略了 Dim[…]符号: 3 3 / 2 2 3 / 2 2 2 2 3 1 Ψ Ψ → = → λ − E E E E E x x . 动能项的期望值为 ( ) ( ) µ λ λ µ λ µ 2 , 2 , 2 2 3 * 2 Ψ = ∆ = Ψ − ∫ d x r r p G . (8.2.12) 相应的 Mathematica V4.0 计算过程为 MATHEMATICA V4.0 In[3]:= [ ] r r E λ π λ ψ λ − = 3 / 2 _, _ : （* 定义“试验”波函数 *）
In[4]:=g]=D,192+2Dm *轨道角动量为零的有效拉普拉斯算符*) InL5]: =4rfrwfr,21glr,2]dr 0ut5=4xx2∥Rel>0- √ Mathematica表达式Dpsi[r,1 ambda],tr,n]]功能为,以r为变量对,)求n次偏导。幂 指数势函数()=ar的期望值为 Jd're'()()( 2)=4mar m2e-lMdr= 24/ma Tn+3) 即 G)-a)r,)=42+3 (8.2.13) 量纲分析要求 即耦合常数a的量纲为 则方程(8.2.13)的量纲也是正确的,即 对应的 Mathematica v4.0的指令为: MATHEMATICA V4. 0 In[6]:= (*积分* Out[6]= /[a]>0&& Re(n)>-3, 2-x-"Gamma(3+nl fe""drl 结合方程(82.12)和(8.2.13),可得能量表示为 2a(n+3) (8.2.14) 对于任何λ>0的值,这一能量解始终是能量值“真解”的上界:En≤E()通过求可使E() 取最小值的变分参数λ值,即解出λ-,就可以很容易地改进这一能量上限,使其逼近真解 显然,λ由下式给出
In[4]:= [ ] [ ] [ ] [ , ],{ ,1} 2 _, _ : [ , ],{ ,2} D r r r g r λ = Dψ r λ r + ψ λ （* 轨道角动量为零的有效拉普拉斯算符 *） In[5]:= 4 r [r, ]g[r, ]dr 0 2 π ψ λ λ ∫ ∞ Out[5]= ] 2 , 4 [Re[ ] 0, 0 5 / 2 7 / 2 3 / 2 2 dr e r e If e r r r r ∫ ∞ − − − > − − + π λ π λ π λ π λ λ λ λ λ 4 Mathematica 表达式 D[psi[r,lambda],{r,n}]功能为，以 r 为变量对Ψ( )求 n 次偏导。幂 指数势函数V 的期望值为 r,λ ( ) n r = ar ( ) ( ) ( ) ( ) ( ) 3 3 0 2 2 3 3 * 2 3 , , 4 4 + ∞ + − Γ + Ψ Ψ = = ∫ ∫ n n r n d x r V r r a r e dr a λ π π λ π π λ λ λ λ . 即 ( ) ( ) ( ) ( ) ( ) ( ) 3 3 * 3 2 3 , , 4 + Γ + ≡ Ψ Ψ = ∫ n n V r d x r V r r a λ λ λ λ . (8.2.13) 量纲分析要求 V → n = ar n aE− → E 即耦合常数 a 的量纲为 a → n+1 E . 则方程(8.2.13)的量纲也是正确的, 即 a −n λ → E E E n n = +1 − . 对应的 Mathematica V4.0 的指令为： ------------------------------------------------------------------------------ MATHEMATICA V4.0 ------------------------------------------------------------------------------ In[6]:= r E dr (* 积分 *) n 2λr 0 2 − ∞ + ∫ Out[6]= [Re[ ] 0& & Re[ ] 3,2 [3 ], ] 2 0 3 3 2 If n Gamma n e r dr n n r +n ∞ − − − − − ∫ > > − + λ λ λ ------------------------------------------------------------------------------ 结合方程(8.2.12)和(8.2.13)，可得能量表示为 ( ) ( ( ) ) n a n µ λ λ λ 2 3 2 2 2 Γ + E = + . （8.2.14） 对于任何λ > 0 的值，这一能量解始终是能量值“真解”的上界：Etrue ≤ E(λ)。通过求可使E(λ) 取最小值的变分参数λ 值，即解出λ min，就可以很容易地改进这一能量上限，使其逼近真解。 显然，λ min由下式给出