§8.2求非相对论性薛定格方程本征能量限 用 Mathematica v4.0系统的指令,对应的计算过程可表述为 MATHEMATICA V4. 0 (*积分*) In[1]: =E-r"dr Out[1]=I[Re[n]>-1&& Re[a]>0,t1-Gammdl+n], Le-ir"dr 这一输出结果的含义是:如果Re[4]>0,且Re>-1,则以上积分的结果为 xr(+n),否则将输出 这意味着 Mathematica无法求解该问题。由此可以得归一化因子 A 对库仑势,归一化的“试验”波函数为 pr 为保险起见,我们可以检验一下关系式两边的量纲。根据以前的讨论,我们知道关系式左 边的量纲为DiE2]。为使指数运算exp{-r}有意义,乘积必须是无量纲的量,即 Dm]=1。由此有Dm]=mm例”即 Dm]=DimE312]=Dm32]。 很显然,在以上推导中至少量纲是正确的。下面我们演示一下如何运用 Mathematica语言 作以上定义和计算 采用 Mathematica v4.0的对应计算为: MATHEMATICA V4. 0 In2:47[rE-2b(*积分*
§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 ，则以上积分的结果为 Γ( + n ，否则将输出 − − 1 1 λ λ ) dr r r n ∫ ∞ 0 exp{λ } , 这意味着Mathematica无法求解该问题。由此可以得归一化因子 π λ 3/ 2 N = ， 对库仑势，归一化的“试验”波函数为 ( ) r r e λ π λ λ − Ψ = 3/ 2 , . 为保险起见，我们可以检验一下关系式两边的量纲。根据以前的讨论，我们知道关系式左 边的量纲为 [ ]。为使指数运算 exp{ 3/ 2 Dim E −λr}有意义，乘积 λr 必须是无量纲的量，即 Dim[λr] = 1。由此有 Dim[E] [ ] 1 [ ] Dim r Dim λ = = ，即 [ ] [ ] [ ]。 3 / 2 3 / 2 Dim Ψ = Dim E = Dim λ 很显然，在以上推导中至少量纲是正确的。下面我们演示一下如何运用 Mathematica 语言 作以上定义和计算。 采用 Mathematica V4.0 的对应计算为： MATHEMATICA V4.0 In[2]:= r E dr （* 积分 *） r ∫ ∞ − 0 2 2 4 λ π
out2]=4x∥Re[]>0 4235e2b 输出的含义是:当ReA>0时,计算结果为x/3,否则, Mathematica无法求解,将返回输 入形式4[e2br2b。完整的结果应是: 下一步,我们将借助引入的“试验”波函数求动能项的期望值。由于我们只讨论基态的能 量本征值,而对基态量子数l=0,此时在径向中心力场势情况下可采用拉普拉斯算子形式 为 △=V2d2,2d 其期望值为 =24xa(2x2-22)2 我们可以看到这里的量纲检验仍然是正确的。我们在下式中省略了Dm[…]符号: 平平→E-E32E2E3=E2→22 动能项的期望值为 p=(dNl, H(,) (82.12) 相应的 Mathematica v4.0计算过程为 MATHEMATICA V4. 0 In[3]:=vx= (*定义“试验”波函数*) In4]:=s,]=D,2]+2pr,l明 (*轨道角动量为零的有效拉普拉斯算符* 4[rwr, a][r, a]dr
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]:= [ ] [ ] [ ] [ , ],{ ,1} 2 _, _ : [ , ],{ ,2} D r r r g r λ = Dψ r λ r + ψ λ （* 轨道角动量为零的有效拉普拉斯算符 *） In[5]:= 4 r [r, ]g[r, ]dr 0 2 π ψ λ λ ∫ ∞
0t=42mR0-z 2e-r25/2 er27 Mathematica表达式Dpsi[r,1 ambda],r,n]]功能为,以r为变量对r,求n次偏导。幂 指数势函数()=ar的期望值为 Jd'rYGr, a)()Y(, )=" dr=24m r(n+3) (》=Jdr(,xy(Pv,4) (n+3) (8.2.13) (2A) 量纲分析要求 即耦合常数a的量纲为 则方程(82.13)的量纲也是正确的,即 E 对应的 Mathematica v4.0的指令为: MATHEMATICA V4. 0 积分*) 0ut6=∥Ra1>0&&Re>-32))am3+n12+"1 结合方程(8.2.12)和(8.2.13),可得能量表示为 E(2) 2 3) (8.2.14 对于任何λ>0的值,这一能量解始终是能量值“真解”的上界:E。≤E()。通过求可使E(a) 取最小值的变分参数λ值,即解出λ-就可以很容易地改进这一能量上限,使其逼近真解 显然,L由下式给出 对公式(8.2.14)求偏导,可得 r(n+3)
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)的量纲也是正确的, 即 n a − λ → 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 E µ λ λ λ 2 3 2 2 2 Γ + = + . （8.2.14） 对于任何λ > 0 的值，这一能量解始终是能量值“真解”的上界：Etrue ≤ E(λ)。通过求可使E(λ) 取最小值的变分参数λ 值，即解出λ min，就可以很容易地改进这一能量上限，使其逼近真解。 显然，λ min由下式给出 ( ) = 0 ∂ ∂ λ E λ . 对公式（8.2.14）求偏导，可得 ( ) ( ) ( ) 0 2 3 1 = Γ + = − ∂ ∂ n+ n an E µ λ λ λ λ
求解该方程得到 (8.2.15) 将这一结果反代回关系式(8.2.14),即可得改进后的能量本征值上限 E=E() anI(m+3)7(12 (8.2.16) 对应的 Mathematica v40的程序为: MATHEMATICA V4. 0 In[7]:=411=2/20)+a12 Gammon+3(2)(*定义函数E()*) In[8]:=D[e[A],] (*对参数λ作微分(注意:De[],{,1}等效于D[e[λ],λ])):* 2--nanx--n Gamma [3+nj 0ut[8] *解方程求*)} In[9]: Solve(a/u-2--"an 2-lm Gamma(3+n]=0,a OutL9-332>(2--anu Gamma[3+nD) In[101:=c(21 "an u Gan3+n])10m(*计算E(a)*) (*注意:指令 Powerexpand[expr]的功能为将所有乘积和指数作幂次展开。%代表 Mathematica输出的最后的一个表达式,在此即为上面最后一个表达式,即0ut[11]}。*) In[11]: Power Expand [% 0ut11=22+a2+n2+"42+ n Gamme(3+n21+22+a2+hn2+2+" Gamma3+n 库仑势:将a=-a和n=-1代入,可得 可以看出,表达式的右边正好是氢原子基态能量,即等号是严格成立的。显然,这是由于 我们恰好选取氢原子基态波函数作为“试验”波函数引来的。对于a=1和a=1情况,能量
求解该方程得到 ( ) 1/( 2) min 1 2 3 + + Γ + = n n anµ n λ . (8.2.15) 将这一结果反代回关系式(8.2.14)，即可得改进后的能量本征值上限 ( ) ( ) + Γ + = = + + + n an n E E n n n n 2 1 2 1 3 2 1 2 /( 2) 1 /( 2) var min µ λ . (8.2.16) 对应的 Mathematica V4.0 的程序为： ------------------------------------------------------------------------------ MATHEMATICA V4.0 ------------------------------------------------------------------------------ In[7]:= e[λ _]:= λ 2 /(2µ) + a / 2 Gamma[n + 3]/(2λ) n (* 定义函数 E(λ) *) In[8]:= D[e[ λ ], λ ] （* 对参数λ 作微分（注意:D[e[ λ ],{ λ ,1}]等效于 D[e[ λ ], λ ]）): *) Out[8]= λ µ − 2−1−n a n λ−1−n Gamma@3 + nD (* 解方程求λ min *)} In[9]:= [ / 2 [3 ] 0, ] 1 1 λ µ − λ + == λ − − − − Solve an Gamma n n n Out[9]= → + − −n +n an Gamma n 2 1 1 λ (2 µ [3 ]) In[10]:= e[(2−1−n an µ Gamma[3 + n])1/(2+n) ] （* 计算 ( ) E λ min *) Out[10]= 2 (2 [3 ]) [3 ] 2 (2 [3 ]) 2 1 1 1 2 2 1 a an Gamma n Gamma n an Gamma n n n n n n n + + + + − − − − − + − − + µ µ µ （* 注意：指令 PowerExpand[expr] 的功能为将所有乘积和指数作幂次展开。% 代表 Mathematica 输出的最后的一个表达式，在此即为上面最后一个表达式，即 Out[11]}。*） In[11]:= PowerExpand[%] Out[11]= ( ) ( ) n n n n n n n n n n n n n n n a n Gamma n a n Gamma n + + − + − + + − − + + − + + + − − + + + 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 4 3 2 µ [3 ] 2 µ [3 ] (*--------------------------------------------------------------------------*) 库仑势：将a = −α 和n = −1代入，可得 2 2 α µ Etrue ≤ − . 可以看出，表达式的右边正好是氢原子基态能量，即等号是严格成立的。显然，这是由于 我们恰好选取氢原子基态波函数作为“试验”波函数引来的。对于α =1和 µ = 1情况，能量
数值解为E=-0.5。这一处理的 MathematicaⅤ3.0表述式为 与上面指令对应的 Mathematica v4.0版本如下: MATHEMATICA V4. 0 In[13]:=e[λ_n_,B,_]:x212a)+a/2Gmdn+3]/(x) (*定义需要最小化的函数* In[14]:= Findminimum[e[λ,-1,-1,1],{A,0.5}] (*以变分参数λ=0.5为起始点,求最小值* 0ut[14]={-0.5,{λ→1.} *这可以理解为m=1时,能量最小值E==-0.5*) 线性势:()=ar。将n=1代入关系式(8.2.16),得到 E_≤E= 将这一结果按关系式(8.26)的格式重写出来,有 24764 (8.2.17) 而基态能量的“真解”E可由Airy函数[8][9]的第一个零点给出 E=2.338 比较两种结果,可以看到我们求得的原始上限E。与真值的相对误差非常小 E-E 由此我们在不具体求解薛定格方程的情况下,解析地推导出线性势基态本征能量,其数值 结果与“真解”有很好的近似。 对应的 Mathematica v4.0版本为: MATHEMAtICAV4 0 In15]:=e[Da_,]:=2120a12 Amman+3]/(2A) (*定义能量函数*) In[6]:= Findminimum[e[λ,1,1,1],{,0.5}] (*以变分参数=0.5为起始点,求最小值
数值解为 E = −0.5。这一处理的 Mathematica V3.0 表述式为 λ ( )r = E λ 与上面指令对应的 Mathematica V4.0 版本如下： ------------------------------------------------------------------------------ MATHEMATICA V4.0 ------------------------------------------------------------------------------ In[13]:= e[ _,n_,a_, µ _] := n /(2 ) a / 2 Gamma[n 3]/(2 ) 2 λ µ + + λ (* 定义需要最小化的函数 *) In[14]:= FindMinimum[e[ λ ,-1,-1,1], { λ ,0.5}] (* 以变分参数λ = 0.5为起始点，求最小值 *) Out[14]= {-0.5, { λ → 1.} （* 这可以理解为 1 λ min = 时，能量最小值 0.5 *） Emin = − 线性势：V ar 。将n = 1代入关系式（8.2.16），得到 1 / 3 2 5 / 3 var 2 3 ≤ = µ a Etrue E . 将这一结果按关系式(8.2.6)的格式重写出来，有 1 / 3 2 1 / 3 2 4 / 3 5 / 3 2 2.4764 2 2 3 = ≤ µ µ a a true . (8.2.17) 而基态能量的“真解” 可由 Airy 函数[8][9]的第一个零点给出， Etrue 1/ 3 2 2 2.3381 = µ a Etrue . 比较两种结果，可以看到我们求得的原始上限 与真值的相对误差非常小. Evar var ≅ 6% − true true E E E . 由此我们在不具体求解薛定格方程的情况下，解析地推导出线性势基态本征能量，其数值 结果与“真解”有很好的近似。 对应的 Mathematica V4.0 版本为： MATHEMATICAV4.0 In[15]:= e[ _,n_,a_, µ _] := n a / 2 Gamma[n 3]/(2 ) 2 /(2 ) λ λ µ + (* 定义能量函数 *) In[16]:= FindMinimum[e[ λ ,1,1,1],{ λ ,0.5}] （* 以变分参数λ = 0.5为起始点，求最小值。 *）
0ut[16]={1.96556,{A->1.14471} 将n=a=1代入公式(8.2.15),(8.2.16)和(8.2.17),可对以上结果进行检验。通常将 计算结果绘图显示,有助于进行比较。令2=1,我们分别绘制函数E=2.3812和 E=2.4764a23,然后借助 Mathematica v3.0的Show指令将两图形合在一起进行比较。 MATHEMATICA V3. 0 定义函数 In[17]:= etrue[a]:=2.3381a^(2/3) In[18]:= upper[a]:=2.4764a^(2/3) In[19]: =plot=Plot [etrue[a], fa, 0, 3), AxesLabe1->("a","E", Textstyle-> Fontslant-> Italic”, Fontsize->14}](*绘制E(plot1)*) (* AxesLabe1:定义坐标轴的表述; Textstyle定义字型和字体大小。* 0.511.522.53 图82.2线性势基态能量真值E In [20]: =plot2=Plot [eupper[a],(a, 0, 31, AxesLabel->("a""E"I, TextStyle-> FontSlant ->"Italic", FontSize->14] (*绘制E(plot2)*) 0.511.5 图823线性势基态的变分上限能量E
Out[16]= {1.96556,{ λ -> 1.14471}} 将 代入公式(8.2.15),(8.2.16)和(8.2.17)，可对以上结果进行检验。通常将 计算结果绘图显示，有助于进行比较。令 n = a = 1 2µ = 1，我们分别绘制函数 和 ，然后借助 Mathematica V3.0 的 Show 指令将两图形合在一起进行比较。 2 / 3 E 2.3381a true = 2 / 3 2.4764a Evar = 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) *) Etrue （* AxesLabel: 定义坐标轴的表述； TextStyle 定义字型和字体大小。*） 图8.2.2 线性势基态能量真值 。 In[20]:=plot2=Plot[eupper[a],{a,0,3},AxesLabel->{”a”,”E”},TextStyle-> 2) *) 图8.2.3 线性势基态的变分上限能量 。 true E {FontSlant ->”Italic”,FontSize->14}] (* 绘制 Evar (plot Evar
In[21]:=Show[plot1,plot2](*将 tt plot1}与[\ tt plot2合并绘制*) 2.5 图8.2.4线性勢基态能量真值和变分上限
In[21]:= Show[plot1,plot2] (* 将{\tt plot1} 与{\tt plot2}合并绘制 *) 图8.2.4 线性势基态能量真值和变分上限