§ 8.2 求非相对论性薛定格方程本征能量限
用 Mathematica V4.0 系统的指令，对应的计算过程可表述为：
MATHEMATICA V4.0
（ * 积分 *）
In[1]:= drrE
nr
∫
∞
?
0
λ
Out[1]= ]],1[,0]Re[&&1][Re[
0
1
drrenGammanIf
nrn
∫
∞
???
+>?>
λ
λλ
这 一输出 结果 的含义 是： 如果 Re[ 0]>λ , 且
1]Re[ ?>n
， 则以上 积分 的结果 为
，否则将输出 (n+Γ
??
1
1 λ
λ )
dr
r
r
n
∫
∞
0
}exp{λ
,
这意味着 Mathematica无法求解该问题。由此可以得归一化因子
π
λ
2/3
=N ，
对 库仑势 ，归一化的 “试验”波函数 为
()
r
er
λ
π
λ
λ
?
=Ψ
2/3
, .
为保险 起见 ，我们 可以 检验一 下关 系式两 边的 量纲。 根据 以前的 讨论 ，我们 知道 关系式 左
边的量 纲为 。为使 指数 运算 exp{][
2/3
EDim
}rλ? 有意义 ，乘积 rλ 必须 是无量 纲的 量，即
1][ =rDim λ 。由此有
][EDim
][
1
][
rDim
==λDim
，即
。 ][][][
2/32/3
λDimEDimDim ==Ψ
很显然， 在 以上推导 中 至少量纲 是 正确的。 下 面我们演 示 一下如何 运 用 Mathematica 语言
作以上定义和计算。
采用 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)的量纲也是正确的 , 即
n
a
?
λ → 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
E
λμ
λ
λ
2
3
22
2
+Γ
+=
. （ 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
(*--------------------------------------------------------------------------*)
库仑势 ：将 α?=a 和 n 代入，可得 1?=
2
2
μα
?≤
true
E
.
可以看 出， 表达式 的右 边正好 是氢 原子基 态能 量，即 等号 是严格 成立 的 。显 然， 这是由 于
我们恰好 选 取氢原子 基 态波函数 作 为“试验 ” 波函数引 来 的 。 对 于 1=α 和 1=μ 情况， 能量
数值解为 。这一处理的 Mathematica V3.0 表述式为 5.0?=E
λ
()r =
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） ，得到
ar
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
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) *)
图 8.2.3 线性势基态的变分上限能量 。
true
E
{FontSlant ->”Italic”,FontSize->14}] (* 绘制
var
E (plot
var
E
In[21]:= Show[plot1,plot2] (* 将 {\ tt plot1}
与 {\ tt plot2}合并绘制 *)
图 8.2.4 线性势基态能量真值和变分上限。
