http://www.paper.edu.cn
基于 MATLAB 的能级、波函数及几率密度图形
的处理
张树林 1,厉树忠 2
1 唐山市经济贸易学校,河北唐山(063700)
2 西北师范大学,兰州(730070)
E-mail:zhangshulin@yeah.net
摘 要:Matlab 语言是由美国的 Clever Moler 博士于 1980 年开发的,它集科学计算、图象
处理、声音处理于一身,并提供了丰富的 Windows 图形界面设计方法。Matlab 语言提供了
一套功能强大的绘图命令,这些命令可以根据输入的数据自动完成图形的绘制,为计算过程
和结果的可视化提供了极佳的手段。本文着重介绍运用 Matlab 编程语言完成量子力学中一
维无限深势阱和线性谐振子的能级、波函数及几率密度图形的绘制问题。
关键词:Matlab 语言;能级;波函数;几率密度
Matlab 语言是由美国的 Clever Moler 博士于 1980 年开发的,其设计初衷是为解决“线性
代数”课程的矩阵运算问题。目前,Matlab 已经不仅仅是 “矩阵实验室”了,它集科学计算、
图象处理、声音处理于一身,并提供了丰富的 Windows 图形界面设计方法。Matlab 语言提
供了一套功能强大的绘图命令,这些命令可以根据输入的数据自动完成图形的绘制,为计算
可以完成对图形加标号、加标题或画上网格标线,可设置不同颜色、线型、视角等功能
过程和结果的可视化提供了极佳的手段,在图形的绘制过程中,用户可选择多种类型的绘图,
。
在量子力学中,用薛定谔方程处理了一维无限深势阱和线性谐振子等简单的定态问题,
并通过研究这些问题,得到了一些严密的结果。本文运用 Matlab 编程语言完成一维无限深
势阱和线性谐振子的能级、波函数及几率密度等结果的图形绘制。
1. 一维无限深势阱
1.1 量子力学问题
[1][2][3]
在许多情况中,如金属中的电子、原子中的电子、原子核中的质子和中子等粒子的运动
有一个共同点,即粒子的运动都被限制在有限的空间范围内,或者说,粒子处于束缚态。为
了分析束缚态粒子的共同特点,我们可以将上述情况简单化、理想化,建立无限深势阱模型
[4]
。
无限深势阱的势能在一定区域内(
粒子的势能为:
xU
)(
=
如图(1)所示:
,0
⎧
⎨
∞
⎩
0
x
x <<0
x
<<
x
,0
≤
a
≥
a
)为零,而在此区域外势能为无限大,即
a
图 1 一维无限深势阱
- 1 -
1.2 无限深势阱问题的具体解法
ax
≥
=xψ
在势阱外:
0)(
≤ ,0
[
x
http://www.paper.edu.cn
] 在势阱内:
=xU
0)(
,故定态薛定谔方程为:
−
2
m
2
d
2
ψ E
=
dx
2
ψ
x <<0
a
,令
k =
mE
2
2
,方程可化为标准形式:
k
+
=
+
kx
0,0
2
ψ
<<
ax
d
2
ψ
dx
2
式中 A ,δ为两个待定常数,根据波函数连续性要求,
)
δ
。令波函数不能恒为零,而 A 不能为零,
0)0( =ψ
得:
A
x
kx
,所以 ka 必
0
sin
n 取负数给不出新的波函数。这告诉我们 k 只能取下列
......
sin =δA
)(
,据
0
=aψ
aψ
)(
sin
式得
ka
A
=
=
0
)( =ψ
3,2,1=n
=
A
ψ
x
)(
其通解为:
sin(
在势阱边界上,由通解和
所以必须 0=δ ,于是
须满足:
值:
ka = ,
πn
n
k π
=
a
3,2,1=n
......
由
k =
mE
2
2
式可知,粒子的能量只能取下列值:
En
=
2
π
n
2
2
2ma
2
3,2,1=n
......
即只有当能量值取
En
=
2
π
n
2
2
2ma
2
所给出那些 nE 值时对应的波函数才有满足边值条件,
这样我们就能很自然地得到,对应于量子数 n 的全部可能值,有无限多个能量值,它们组成
体系的分立能级。而被束缚在阱中的粒子的能量也只能取一系列离散的数值,即能量是量子
化的。
代入到
=ψ
x
)(
A
sin
kx
中,并把势阱外的波函数也包括在内,便得到能量为
最后得到能量为 nE 的归一化波函数为:
ψ
x
)(
=
1.3 能量、波函数和几率密度的 MATLAB 实现
1.3.1 能量的 MATLAB 实现
一维无限深势阱能量公式
En
=
2
π
n
2
2
ma
2
2
=
hn
2
2
ma
8
2
- 2 -
,式中 h 是普朗克常数,m 为粒子
将
n
k π
=
a
En
=
2
π
n
2
2
2ma
2
n
≠
,0
=
ψn
,0
∫∞
∞−
的 波 函 数 :
ψ
xn
)(
=
⎧
0
⎪⎪
⎨
⎪
⎪
⎩
A
sin
n
π
a
x
x
≤
,0
x
≥
a
0
<<
x
a
3,2,1=n
......
≡
0
,波函数无意义。波函数中 A 可由归一化条件确定
|
ψ
n
x
|)(
2
dx
=
a
∫
0
|
ψ
x
|)(
2
dx
=
A
a
0
2
sin
n
π
a
dxx
=
1
知:
A
=
2
a
2
∫
⎧
0
⎪⎪
⎨
⎪
⎪
⎩
x
≤
,0
x
≥
a
。
0
<<
x
a
sin2
a
n
π
a
x
http://www.paper.edu.cn
5,4,3,2,1=n
的质量,a 是势阱的直径,我们取一个质量单位 u=1.66e-27 千克, a=1 个长度单位。下面我
时的图形。在 MATLAB 命令窗口中输入如下程序:
们应用 MATLAB 软件作出
>>clear,close all
E=zeros(1,5);
for n=1:5
a=1;
E(n,a)=(6.626e-34).^2*n.^2/(8*(1.66e-27)*a.^2);
end;stem(E,'.')
>> ylabel('能量 E');xlabel('N')
按回车键后,在图形窗口中出现如下图形(2):
E
量
能
9
8
7
6
5
4
3
2
1
0
x 10-40
1
1.5
2
2.5
3
N
3.5
4
4.5
5
图 2 一维无限深势阱的分立能级
我们根据图形可以容易直观地得出讨论:1、能量量子化
π
n
2
2
ma
2
2
hn
2
2
ma
8
=
2
n=1,2,3... 。一维势阱中粒子的能量必须是量子化的,只
2
En
=
能取分立值。2、n=1 的状态称粒子的基态,是能量最低的状态,这个最低能量
E =
1
h
2
8ma
2
称作零点能,它表明粒子永远在运动着,零点能的存在是不确定关系的必然结果。
1.3.2 波函数和节点的 MATLAB 实现
一维势阱中粒子能量为 nE 的归一化波函数为:
ψ
x
)(
=
⎧
0
⎪⎪
⎨
⎪
⎪
⎩
sin2
a
n
π
a
x
x
≤
,0
x
≥
a
,
0
<<
x
a
在势阱内是正弦函数,下图描述 n=1,2,3 时三个波函数。在 MATLAB 命令窗口中输入如
下程序:
>>clear,close all
x=0:1/100:1;a=1;
n=1;y1=sqrt(2/a).*sin(n.*pi.*x./a);
n=2;y2=sqrt(2/a).*sin(n.*pi.*x./a);
n=3;y3=sqrt(2/a).*sin(n.*pi.*x./a);
n=4;y4=sqrt(2/a).*sin(n.*pi.*x./a);
- 3 -
http://www.paper.edu.cn
subplot(221),plot(x,y1);xlabel('x'),ylabel('\psi1');grid on;
subplot(222),plot(x,y2); xlabel('x'),ylabel('\psi2');grid on;
subplot(223),plot(x,y3); xlabel('x'),ylabel('\psi3');grid on;
subplot(224),plot(x,y4); xlabel('x'),ylabel('\psi4');grid on
按回车键后,在图形窗口中出现如下图形(3):
1.5
1
0.5
1
ψ
0
0
2
1
0
-1
-2
0
3
ψ
2
1
0
-1
-2
0
2
1
0
-1
-2
0
2
ψ
4
ψ
1
1
0.5
x
0.5
x
图 3 一维无限深势阱的能量本征函数
0.5
x
0.5
x
1
1
根据图形可以容易直观地得出结论:①、节点数为量子数 n-1 个,即 n=1,2,3...节点数为
0,1,2...
②、波长 λ 也随量子数的改变而改变,当 n=1 时,λ=2a, 当 n=2 时,λ=a,当 n=3 时,λ=2a/3 ...
1.3.3 几率密度分布的 MATLAB 实现
一维无限深势阱的几率密度公式为 =nJ
2)(xnψ ,是波函数的平方,因此我们在波函
数程序的基础上很方便地得到其 MATLAB 实现。在 MATLAB 指令窗口中输入下列程序:
>> clear,close all
x=0:1/100:1;a=1;
n=1;Y1=(sqrt(2/a).*sin(n.*pi.*x./a)).^2;
n=2;Y2=(sqrt(2/a).*sin(n.*pi.*x./a)).^2;
n=3;Y3=(sqrt(2/a).*sin(n.*pi.*x./a)).^2;
n=4;Y4=(sqrt(2/a).*sin(n.*pi.*x./a)).^2;
subplot(221),plot(x,Y1);xlabel('x'),ylabel('|\psi1|^2');grid on;
subplot(222),plot(x,Y2); xlabel('x'),ylabel('|\psi2|^2');grid on;
subplot(223),plot(x,Y3); xlabel('x'),ylabel('|\psi3|^2');grid on;
subplot(224),plot(x,Y4); xlabel('x'),ylabel('|\psi4|^2');grid on
>>按回车键后,在图形窗口中得到图形(4):
- 4 -
|ψ1|2
2
ψ3|
|
3
2
1
0
0
3
2
1
0
0
0.5
x
0.5
x
1
1
|ψ2|2
2
ψ4|
|
3
2
1
0
0
2
1.5
1
0.5
0
0
图 3 一维无限深势阱中粒子位置几率密度分布
http://www.paper.edu.cn
0.5
x
0.5
x
1
1
根据图形可以直观地得出一维无限深势阱的几率密度,节点数与量子数 n 的关系为 n-1
个,即 n=1,2,3,4 节点数为 0,1,2,3。
2. 线性谐振子
2.1 量子力学问题
如果在一维空间内运动的粒子的势能为
xU
)(
为线性谐振子
。该体系的定态薛定谔方程为
[4]
−
2
2
µ
+
1
2
2
ψψµω
=
E
x
2
。
2
µω=
2
x
,ω 是常数,则这种体系就称
1
2
d
2
ψ
dx
2
2.2 线性谐振子问题的理论解法
先将定态薛定谔方程
−
2
2
µ
d
2
ψ
dx
2
+
1
2
引进无量纲的参数
ααξ
= x
=
此时 α=3.9663e+003),则上方程变成:
µω
d
2
ψ
d
2
ξ
2
ψψµω
=
E
x
2
简化。
和
λ
=
E2
ω
(因为 ω 是常数,不妨设为 ω=1,
+
[
ψξλ
−
]
2
=
0
d
2
ψ
d
2
ξ
2/2ξψ ±= e
2/2ξψ e=
,其中
不满足边
当ξ很大时,λ与ξ相比可以忽略,方程进一步近似为:
2
− ψξ
=
0
(1)
可证,当
±∞→ξ
值条件,故只能取:
时,方程(1)的渐近解为:
2/2ξψ −= e
在渐进解形式的启发下,令方程
将它代入方程。得:
Hd
2
d
2
ξ
−
2
ξ
dH
d
ξ
d
2
ψ
d
2
ξ
+
[
ψξλ
−
]
2
=
0
的精解形式为:
ψ ξ H
e−=
2/2
)(
ξ
,
+
[
−
λ
]1
H
=
0
是厄密方程,解为
)(ξH ,从而得ψ。
- 5 -
http://www.paper.edu.cn
分析:只有当
λ
21
=−
n
n
3210
、、、=
......
时,方程才能满足要求,此时,方程的解
为厄密多项式,通常认为:
H
n
)(
ξ
−=
)1(
n
e
2
ξ
n
d
2
d
ξ
2
ξ
−
e
它是ξ的 n 次多项式,如前四项:
=
;1
H
H
=
;2
ξ
1
0
H
2
=
4
2
ξ
进一步可得一维谐振子的能量可能取值为:
En
=
对应于能量 nE 的波函数是:
ψ
n
x
)(
=
eN
n
2
α−
x
2/2
H
3
=
8
2
ξ
−
12
ξ
n
3210
、、、=
......
ω
−
;2
1
2
xH
(
α
+
n
)
(
)
n
2.3 能量、波函数和几率密度的 MATLAB 实现
2.3.1 能量的 MATLAB 实现
在能量式
En
=
(
n
+
1
2
)
ω
n
3210
、、、=
......
中,
=
h
π2
,取 1=ω ,h 是普朗克常数,
在 MATLAB 中 , n 不 能 取 0 值 , 为 得 到 正 确 的 能 量 图 形 , 我 们 做 修 改 如
下:
En
=
[(
n
1)1
+−
2
]
h
π2
ω
。在 MATLAB 指令窗口中输入下列程序:
>> clear,close all
E=zeros(1,5);
for n=1:5
E(n)=((n-1)+1/2)*(6.626e-34)/(2*pi);
end;stem(E,'.')
ylabel('能量 E');xlabel('N')
按回车键后,在图形窗口中得到图形(4):
x 10-34
5
4.5
4
3.5
3
2.5
2
1.5
1
0.5
E
量
能
0
1
1.5
2
2.5
3
N
3.5
4
4.5
5
图 4 线性谐振子的能级分布
通过对图形的观察,很容易得出结论:线性谐振子的能量只能取分立值且振子的基态
(n=0)能量
E =
0
。
ω
h
π4
- 6 -
2.3.2 波函数的 MATLAB 实现
在波函数
ψ
n
x
)(
=
eN
n
α−
x
2
2/2
http://www.paper.edu.cn
H
n
x
(
α
)
中,归一化因子
N
n
=
α
n
n
2!
π
针对
不同的 n 值有不同的取值,在 ω=1,α=3.9663e+003 的情况下,我们取得 N=0、1、2、3 时
的 nN 值分别为:N0=47.3048,N1=33.4496;N2=16.7248;N3=6.8279。同时厄密多项式也不
相同,因此我们做分别编写程序处理。其 MATLAB 实现为:
>>clear, close all
x=-0.002:0.002/100:0.002;
h0=47.3048./exp(((3.9663e+003)*x).^2./2);
subplot(221),plot(x,h0);xlabel('\xi'),ylabel('\psi0');
legend('n=0');grid on
x=-0.05:0.05/100:0.05;
h1=33.4496*2*((3.9663e+003)*x)./exp((3.9663e+003)*x.^2./2);
subplot(222),plot(x,h1);xlabel('\xi'),ylabel('\psi1');
legend('n=1');grid on
x=-0.002:0.002/100:0.002;
h2=16.7248*(4*((3.9663e+003)*x).^2-2)./exp(((3.9663e+003)*x).^2./2);
subplot(223),plot(x,h2);xlabel('\xi'),ylabel('\psi2');
legend('n=2');grid on
x=-0.001:0.001/100:0.001;
h3=(6.8279*(8*((3.9663e+003)*x).^3-12*(3.9663e+003)*x))./exp(((3.9663e+003)*x).^2./2);
subplot(224),plot(x,h3);xlabel('\xi'),ylabel('\psi3');
legend('n=3');grid on
按回车键后,在图形窗口中得到图形(5):
60
40
20
ψ0
0
-2
-1
ψ2
40
20
0
-20
-40
-2
-1
ψ1
ψ3
4000
2000
0
-2000
-4000
-0.05
40
20
0
-20
-40
-1
n=0
1
2
x 10-3
n=2
1
2
x 10-3
0
ξ
0
ξ
n=1
0
ξ
0.05
n=3
-0.5
0
ξ
0.5
1
x 10-3
图 5 线性谐振子的前四个本征函数
通过对图形的观察,很容易得出结论:线性谐振子的波函数
)(ξψn
在有限范围内与 ξ
轴相交 n 次,即
)(ξψn
=0 有 n 个根,或者说
)(ξψn
有 n 个节点。
- 7 -
2.3.3 几率密度的 MATLAB 实现
http://www.paper.edu.cn
2
2
α−
x
2
n
|
=
2/
|)
H
eN
n
,是波函数的
2)(ξψn
线性谐振子的几率密度公式为 =nJ
x
(
α
平方,因此我们在波函数程序的基础上很方便地得到其 MATLAB 实现:
>> clear, close all
x=-0.002:0.002/100:0.002;
h0=(47.3048./exp(((3.9663e+003)*x).^2./2)).^2;
subplot(221),plot(x,h0);xlabel('\xi'),ylabel('|\psi|^2');
legend('n=0');grid on;
x=-0.05:0.05/100:0.05;
h1=(33.4496*2*((3.9663e+003)*x)./exp((3.9663e+003)*x.^2./2)).^2;
subplot(222),plot(x,h1);xlabel('\xi'),ylabel('|\psi|^2');
legend('n=1');grid on;
x=-0.002:0.002/100:0.002;
h2=(16.7248*(4*((3.9663e+003)*x).^2-2)./exp(((3.9663e+003)*x).^2./2)).^2;
subplot(223),plot(x,h2);xlabel('\xi'),ylabel('|\psi|^2');
legend('n=2');grid on;
x=-0.001:0.001/100:0.001;
h3=((6.8279*(8*((3.9663e+003)*x).^3-12*(3.9663e+003)*x))./exp(((3.9663e+003)*x).^2./2)).^2;
subplot(224),plot(x,h3);xlabel('\xi'),ylabel('|\psi|^2');
legend('n=3');grid on
按回车键后,在图形窗口中得到图形(6):
2
|
ψ
|
2
|
ψ
|
3000
2000
1000
0
-2
-1
1500
1000
500
0
-2
-1
0
ξ
0
ξ
n=0
1
2
x 10-3
n=2
x 106
8
6
4
2
2
|
ψ
|
0
-0.05
1500
1000
500
2
|
ψ
|
n=1
0
ξ
0.05
n=3
1
2
x 10-3
0
-1
-0.5
图 6 线性谐振子的前四个几率密度
0
ξ
0.5
1
x 10-3
观察图形得出结论:线性谐振子的几率密度 =nJ
2)(ξψn
在有限范围内与 ξ 轴相交也
为 n 次,即
2)(ξψn
=0 有 n 个根,或者说也有 n 个节点。
- 8 -