MATLAB在量子力学学习中的应用

周群益，莫云飞，周丽丽，侯兆阳

(1.广州理工学院，广东 广州 510540；2.长沙学院 电子信息与电气工程学院，湖南 长沙 410022；3.赣南医学院 信息工程学院，江西 赣州 341000；4.长安大学 理学院应用物理系，陕西 西安 710064)

量子力学中应用了大量的数学公式，朱涵睿等作者总结了一些常用数学工具.这是一项有意义的工作.可惜，公式中出现了一些明显的差错[1]，有必要更正.MATLAB是一个强大的数学工具，可以简单地解决许多计算问题.

1 常用极限的差错

(1)

因此，sinx与x是等价无穷小量，表示为

sinx～x，(x→0)

(2)

其他8个函数的极限和等价无穷小量分别为：

(3)

(4)

(5)

(6)

(7)

(8)

(9)

(10)

检验极限正确与否的最简单的方法就是用MATLAB计算.例如，对于比正弦函数，在命令窗口用下列指令即可计算极限：

>>syms x,y=sin(x)/x,L=limit(y,x,0)

其中，“>>”是MATLAB的提示符.指令syms x定义符号变量；y表示符号表达式；limit函数求极限，第1个参数y是表达式，第2个参数x是自变量，第3个参数0表示自变量的极限；L表示函数的极限.这些指令和函数可通过help指令了解其功能.例如，在命令窗口输入：

>>help limit

回车就能说明函数limit的功能和格式.用如下指令还能画出曲线：

>>figure,ezplot(y),grid

>>hold on,stem(0,L)

其中，指令figure开创图形窗口；指令ezplot是简易绘图指令，第1个参数是表达式，自变量的默认范围是[-2π,2π]；指令grid加网格；指令hold on保留已画的曲线；指令stem画杆线，第1是参数是横坐标，第2个参数是纵坐标.如图1所示，表头显示比正弦函数，此函数是偶对称曲线，当x→0时的极限是1.

图1 比正弦函数的曲线和x趋于0时的极限

2 常微分方法的求法

文献[1]列举了各类常微分方程的求法公式，这都是高等数学教材的内容.其实，对于有显函数解的常微分方程，可以用dsolve函数直接求解.例如，一阶线性微分方程：

(11)

可化为标准形式：

(12)

根据公式可解得y=(ex+C)(x+ 1)2.

利用MATLAB的指令：

>>y=dsolve(′(x+1)*Dy=2*y+exp(x)*(x+1)^3′,′x′)

回车立即可得

y =

exp(x)*(x + 1)^2 + C1*(x + 1)^2

函数中的第1个参数是微分方程，第2个参数是自变量，两项都放在单引号中.系统默认变量为t，所以第二项′x′是必不可少的.

不论是伯努利方程还是常系数线性微分方程，不论是齐次还是非齐次微分方程，有显函数解的常微分方程都能用dsolve函数求解.

3 欧拉公式的验证

MATLAB可对复数进行计算，很容易验证欧拉公式.例如：

eix=cosx+ isinx

(13)

利用MATLAB的指令：

>>y=cos(x)+i*sin(x)-exp(i*x)

>>simplify(y)

化简的结果为0，这就验证了欧拉公式.i是MATLAB默认的虚数单位.对于著名公式：

eiπ+1=0

(14)

输入MATLAB的函数：

>>exp(i*sym(pi))+1

结果就是0.注意：pi是MATLAB的是默认的变量，表示圆周率，本来是一个近似的数值，sym(pi)将数值转化为符号，可以得到精确的符号结果.直接用如下函数：

>>exp(i*pi)+1

只能得到很小的复数值.

4 常用不定积分的计算

>>int(sin(x))

回车立即可得

-cos(x)

>>int(cot(x))

回车立即可得

log(sin(x))

5 常用定积分的计算

MATLAB的函数int很容易计算常用的定积分.例如，欲检验无穷积分公式：

(15)

在命令窗口输入指令：

>>int(x^2/(1+x^2)^4,0,Inf)

回车立即可得

pi/32

就验证了式(15).用同样的方法可以验证无穷积分公式：

(16)

用int函数还能验证三角函数的正交性.

6 常用的级数求和公式的计算

MATLAB的级数求和函数为symsum，如果要验证公式：

(17)

可以在命令窗口用如下指令：

>>syms n,symsum(1/n^2,n,1,Inf)

结果为

pi^2/6

用同样的方法可以验证公式：

(18)

如果要验证公式：

(19)

可以在命令窗口用如下指令：

>>symsum(1/(2*n)^2,1,Inf)

>>symsum(1/(2*n-1)^2,1,Inf)

>>syms n p positive,symsum(1/n^p,n,1,Inf)

其中，参数positive表示变量n和p是正数.结果是

piecewise(1

这是条件函数或分段函数，意思是当1

>>figure,ezplot(zeta(p),[1,10]),grid

如图2所示，表头显示函数ζ(p)，p不限于偶数，也不限于正整数，ζ(p)随着p的增加而减小，当p→+∞时，ζ(p)→1.

图2 ζ(p)函数的曲线和极限

如果在命令窗口用如下指令：

>>syms i,symsum(2*i+1,0,n-1)

结果为

n + n*(n -1)

7 梯度算符和拉普拉斯算符

>>syms x y z,r=sqrt(x^2+y^2+z^2)

>>gradient(r)

结果为

x/(x^2 + y^2 + z^2)^(1/2)

y/(x^2 + y^2 + z^2)^(1/2)

z/(x^2 + y^2 + z^2)^(1/2)

gradient函数还能对数值矩阵求数值梯度.拉普拉斯算符为del2，只能对数值矩阵求值.

8 狄拉克函数的计算

MATLAB的狄拉克函数为dirac，在命令窗口输入指令：

>>syms x a,int(dirac(x-a),-Inf,Inf)

回车可得1，这就验证了公式：

(20)

在命令窗口输入指令

>>int(dirac(x-a),a-eps,a+eps)

回车可得1.指令中，eps是MATLAB的默认变量，表示数值之间的最小间隔，其值为2.2204e-16，即2.204×10-16.这就验证了如下公式：

(21)

在命令窗口输入指令：

>>syms f(x),int(f(x)*dirac(x-a),x,-Inf,Inf)

回车可得

f(a)

这就验证了如下公式：

(22)

在命令窗口输入指令：

>>dirac(-x)

回车可得

dirac(x)

说明狄拉克函数是偶函数:

δ(-x) =δ(x)

(23)

在命令窗口输入指令：

>>d=f(x)*dirac(x-a)-f(a)*dirac(x-a)

>>simplify(d)

回车可得0，这就验证了如下公式：

f(x)δ(x-a)=f(a)δ(x-a)

(24)

不过，用MATLAB不能验证下式：

(25)

也不能根据极限验证下式：

(26)

这是因为δ(x)函数是一种广义函数.

9 伽玛函数和贝塔函数的计算

1)伽玛函数定义为

(27)

在命令窗口输入指令：

>>syms x t positive

>>int(exp(-t)*t^(x-1),t,0,Inf)

回车可得

gamma(x)

输入指令：

>>figure,ezplot(gamma(x)),grid

可画出函数曲线.如图3所示，Γ函数在区间(0,+∞)是先降后升的曲线，通过解析延拓，Γ函数可以扩展到小于零的区域.

图3 Γ函数的曲线和解析延拓曲线

输入指令

>>gamma(sym(1/2))

回车可得

pi^(1/2)

>>syms p,g=gamma(p)*gamma(1-p)

>>simplify(g)

回车可得

pi/sin(pi*p)

这就验证了公式：

(28)

为了验证定积分：

(29)

在命令窗口输入指令：

>>syms x n a positive

>>int(x^(2*n)*exp(-a*x^2),0,Inf)

回车可得

gamma(n + 1/2)/(2*a^(n + 1/2))

这个结果比式(29)的右边的公式更简单.

2) B函数定义为

(30)

B函数与Γ函数的关系为

(31)

B函数是一个二元函数，用如下指令可以画出B函数的曲面：

>>syms x y,figure,ezsurf(beta(x,y),[0,2]),alpha(0.5)

其中，ezsurf指令画二元函数的曲面，第1个参数是二元函数，第2个参数是两个自变量的范围；alpha(0.5)指令使曲面半透明.如图4所示，B函数随着x或y的增加而减小，当x或y→+∞时，B→0.

图4 B函数的曲面

为了检验正弦幂函数和余弦幂函数的定积分公式：

(32)

在命令窗口输入指令：

>>syms x n positive,In=int(sin(x)^n,x,0,pi/2)

或者

>>In=int(cos(x)^n,x,0,pi/2)

回车可得

beta(1/2,n/2 + 1/2)/2

这是用B函数表示的结果，n不分奇数和偶数，因而更简单.用如下指令可以画出曲线：

>>figure,ezplot(In,[0,10]),grid

如图5所示，n不限于整数，In随着n的增加而减小，当n→+∞时，In→0.

图5 正弦和余弦幂函数的定积分曲线In

10 特殊函数的微分方程的求解

在量子力学的计算中要解一些微分方程并用到一些特殊函数[4-6].

1)厄密微分方程为

y″-2xy′+2ny=0

(33)

用下列函数可以求方程的解：

>>y=dsolve(′D2y-2*x*Dy+2*n*y′,′x′)

结果很长.用函数hermiteH可以计算厄密多项式的值.

2) 贝塞尔微分方程为

x2y″ +xy′ +(x2-n2)y=0

(34)

用下列函数可以求方程的解：

>>y=dsolve(′x^2*D2y+x*Dy+(x^2-n^2)*y′,′x′)

结果为

C1*besselj(n,x) + C2*bessely(n,x)

其中，besselj(n,x)表示第n阶第一类贝塞尔函数，bessely(n,x)表示第n阶第二类贝塞尔函数.贝塞尔多项式的值可用函数besselj和bessely计算.

3) 连带勒让德微分方程为

(35)

用下列函数可以求方程的解：

>>y=dsolve(′(1-x^2)*D2y-2*x*Dy+

(l*(1+l)-m^2/(1-x^2))*y′,′x′)

结果也很长.用函数legendre可以计算连带勒让德多项式的值.

4) 拉盖尔微分方程为

xy″+(1-x)y′+ny=0

(36)

用下列函数可以求方程的解：

>>y=dsolve(′x*D2y+(1-x)*Dy+n*y′,′x′)

结果也很长.用函数laguerreL可以计算拉盖尔多项式的值.

11 MATLAB在量子力学中的应用举例

以线性谐振子为例[4]，其薛定谔方程为

(37)

其中，x是位移，ψ是定态波函数，ћ=h/2π是约化普朗克常数，m是振子的质量，E是能量.U是其势能：

(38)

其中，ω是振子的圆频率.设一个由能量决定的无量纲参数：

(39)

再设

(40)

α的倒数具有长度的量纲.还设一个无量纲的位移ξ=αx，则薛定谔方程可简化为

(41)

要使方程有解，能量必须取分立值：

(42)

其中，n是自然数0,1,2,….线性谐振子的波函数为

(43)

其中，Nn是归一化常数：

(44)

而Hn是n阶厄密多项式.概率密度为

(45)

利用MATLAB编程如下：

n=8;syms x%量子数,定义符号变量

Hn=hermiteH(n,x)%厄密多项式

Nn=sqrt(1/sqrt(sym(pi))/2^n/factorial(n));%归一化常数

Pn=Nn*exp(-x^2/2)*Hn%波函数

figure,ezplot(Pn)%创建图形窗口,画波函数曲线

hold on,ezplot(Pn^2)%保持图像,画概率密度曲线

grid on,axis tight%加网格,贴框

当量子数n取8时，如图6所示，波函数是波动的，有8个零点，其值大约在-0.5到0.6之间；概率密度也是波动的，其值不小于零.当n取不同的正整数时，可得类似的曲线(图略).增加一些指令，还能将量子化的概率密度与经典粒子的概率密度进行比较.

图6 一维线性谐振子的波函数和概率密度

勒让德多项式和拉盖尔多项式在氢原子角向概率密度和径向概率密度的计算中具有重要应用.

量子力学的计算比较复杂，我们精选了10个范例[7]，从黑体辐射到氢原子光谱的规律，从德布罗意波长到隧道效应，全部做到了可视化.特别是氢原子的量子力学处理，对于不同的量子数，用MATLAB绘制的概率密度的等值面图，形象描述其分布规律.本书有200多个MATLAB程序，其中有10多个关于量子力学的程序，可供读者参考.

12 结论

在量子力学的学习中要用大量数学工具，手工计算是必不可少的，更重要的是充分利用计算机语言这种工具.MATLAB是一种计算和绘图功能都十分强大的工具，不但能够检错，还可以解决许多复杂的计算问题，被誉为“工科神器”.我们建议在大一就开设MATLAB课程，并且将各门学科与MATLAB编程方法结合起来，提高学生提出问题和解决问题的能力.