摘要:此文介绍了一种使用MATLAB求解一维定态薛定谔方程的方法。利用充分格式进行离散化,得出相应的矩阵方程,用MATLAB求解本征值和本征函数。此方法简单可靠,可以处理各种时间无关的束缚态问题。所用的程序附于文后。
宏观物体遵循的是牛顿运动方程,给定初始条件以及受力条件,我们就可以求出任意时刻粒子的状态。而在原子尺度上,所有粒子都表现出波的行为,粒子的状态用波函数
![]()
来描述,描述微观粒子运动的方程是由薛定谔于1925年提出薛定谔方程。微观粒子的运动与其所处的势场相关,当势场不随时间变化时,称为定态,一维情况下的定态薛定谔方程为
其中,
![]()
表示粒子所处的势场。由定态波函数可以得出总波函数
1926年, Born提出,波函数模的平方
![]()
代表在位置
![]()
,时刻
![]()
寻找到电子的概率,这就是波函数的条统计解释。可以由遵循微分方程的波函数表示,薛定谔波方程涉及空间坐标和时间。对于一维情况,在时刻
![]()
时,在
![]()
和
![]()
之间的某处找到电子的概率由下式给出
时间无关薛定谔方程(1)是系统的能量本征方程。该特征值方程通常由一组特定的函数
![]()
和相应的一组常数
![]()
满足解的条件,被称为是哈密顿算符的本征函数和相应的本征值。当测量处于
![]()
状态的系统的能量时,结果将始终为
![]()
。对于束缚态情况,必须有
一维定态薛定谔方程是二阶微分方程,但是,能够解析求解的情况屈指可数,如氢原子,谐振子,无限深势阱等。随着计算机技术的发展,我们可以数值上进行求解。本文利用MATLAB软件,使用矩阵方法求解束缚态的本征值问题。
对于原子系统,以nm和eV的能量测量长度更方便。我们可以使用缩放因子
因此我们可以将等式(1)写成
为了求解上述方程式,首先进行离散化处理。将坐标
![]()
离散化为
![]()
个点,用
![]()
来表示,若
![]()
,则有
![]()
另外,还需要对方程式进行差分格式处理,二阶导数处理如下
因此,
![]()
的二阶导数矩阵可以写成
![]()
矩阵
![]()
矩阵大小为
![]()
而不是
![]()
,因为函数的二阶导数无法在终点进行计算,即
![]()
和
![]()
。动能矩阵为
![]()
,势能矩阵为对角矩阵,即
![]()
。则我们现在可以将哈密顿矩阵定义为
![]()
。用于生成哈密顿矩阵的代码是
% Make Second Derivative Matrix ---------------
off = ones(num-3,1);
SD_matrix = (-2*eye(num-2) + diag(off,1) + diag(off,-1))/dx2;
% Make KE Matrix
K_matrix = Cse * SD_matrix;
% Make Hamiltonian Matrix
H_matrix = K_matrix + U_matrix;
因此,矩阵形式的薛定谔方程是
![]()
。MATLAB有内置函数可以求解本征值问题,其代码为
[e_funct, e_values] = eig(H)
其中e_funct是具有对应于第n个本征函数的第n列的
![]()
矩阵,并且e_values是按递增顺序的N个本征值的列向量。一般可以求解出N-2个本征值和本征函数,但只有e_values的负值才有意义。为了获得完整的特征向量,我们需要包括端点,其中
![]()
。
接下来举一个例子,计算有限深方势阱问题。如图所示是求解得到的有限深方势阱的本征能量谱。
有限深方势阱的本征能量谱同时还可以求出本征函数以及几率分布,如图所示,
本征波函数以及几率分布我们还可以改变势函数去计算各种各样的定态束缚态问题,也可以去计算已知解的问题,以验证此方法的可靠性。
程序:
clear;
clc;
tic;
num = 2001; % Number of data points (odd number)
% Constants -----------------
hbar = 1.055e-34;
e = 1.602e-19;
m = 9.109e-31;
eps0 = 8.854e-12;
Ese = 1.6e-19; % Energy scaling factor
Lse = 1e-9; % Length scaling factor
Cse = -hbar^2/(2*m) / (Lse^2*Ese); % Schrodinger Eq constant
% Potential well parameters
U = zeros(num,1);
U_matrix = zeros(num-2);
% Potential Wells square well
% Enter energies in eV and distances in nm
xMin = -0.1;
xMax = +0.1;
x1 = 0.05; % 1/2 well width
U1 = -400; % Depth of well (eV)
x = linspace(xMin,xMax, num);
for cn = 1 : num
if abs(x(cn)) <= x1
U(cn) = U1;
end
end
s = sprintf('Potential Well: SQUARE');
% Graphics -----------------------
figure(1);
set(gcf,'Name','Potential Energy','NumberTitle','off')
plot(x,U,'LineWidth',3);
axis([xMin-eps xMax min(U)-50 max(U)+50]);
title(s);
xlabel('x (nm)');
ylabel('energy (eV)');
grid on
% Make potential energy matrix
dx = (x(2)-x(1));
dx2 = dx^2;
for cn =1:(num-2)
U_matrix(cn,cn) = U(cn+1);
end
% Make Second Derivative Matrix
off = ones(num-3,1);
SD_matrix = (-2*eye(num-2) + diag(off,1) + diag(off,-1))/dx2;
% Make KE Matrix
K_matrix = Cse * SD_matrix;
% Make Hamiltonian Matrix
H_matrix = K_matrix + U_matrix;
% Find Eignevalues E_n and Eigenfunctions psi_N
[e_funct, e_values] = eig(H_matrix);
% All Eigenvalues 1, 2 , ... n where E_N < 0
flag = 0;
n = 1;
while flag == 0
E(n) = e_values(n,n);
if E(n) > 0
flag = 1;
end % if
n = n + 1;
end % while
E(n-1) = [];
n = n-2;
% Corresponding Eigenfunctions 1, 2, ... ,n: Normalizing the wavefunction
for cn = 1 : n
psi(:,cn) = [0; e_funct(:,cn); 0];
area = simpson1d((psi(:,cn) .* psi(:,cn))',xMin,xMax);
psi(:,cn) = psi(:,cn)/sqrt(area); % normalize
prob(:,cn) = psi(:,cn) .* psi(:,cn);
if psi(5,cn) < 0
psi(:,cn) = -psi(:,cn);
end % curve starts positive
end % for
% Display eigenvalues in Command Window
disp(' ');
disp('========================= ');
disp(' ');
fprintf('No. bound states found = %0.0g n',n);
disp(' ');
disp('Quantum State / Eigenvalues En (eV)');
for cn = 1 : n
fprintf(' %0.0f ',cn);
fprintf(' %0.5g n',E(cn));
end
disp(' ')
disp(' ');
% Plot energy spectrum
xs(1) = xMin;
xs(2) = xMax;
figure(2);
set(gcf,'Units','Normalized');
set(gcf,'Position',[0.5 0.1 0.4 0.6]);
set(gcf,'Name','Energy Spectrum','NumberTitle','off')
set(gcf,'color',[1 1 1]);
set(gca,'fontSize',12);
plot(x,U,'b','LineWidth',2);
xlabel('position x (nm)','FontSize',12);
ylabel('energy U, E_n (eV)','FontSize',12);
h_title = title(s);
set(h_title,'FontSize',12);
hold on
cnmax = length(E);
for cn = 1 : cnmax
ys(1) = E(cn);
ys(2) = ys(1);
plot(xs,ys,'r','LineWidth',2);
end %for
axis([xMin-eps xMax min(U)-50 max(U)+50]);
% Plots first 5 wavefunctions & probability density functions
if n < 6
nMax = n;
else
nMax = 5;
end
figure(3)
clf
set(gcf,'Units','Normalized');
set(gcf,'Position',[0.05 0.1 0.4 0.6]);
set(gcf,'NumberTitle','off');
set(gcf,'Name','Eigenvectors & Prob. densities');
set(gcf,'Color',[1 1 1]);
%nMax = 8;
for cn = 1:nMax
subplot(nMax,2,2*cn-1);
y1 = psi(:,cn) ./ (max(psi(:,cn)-min(psi(:,cn))));
y2 = 1 + 2 * U ./ (max(U) - min(U));
plot(x,y1,'lineWidth',2)
hold on
plot(x,y2,'r','lineWidth',1)
%plotyy(x,psi(:,cn),x,U);
axis off
%title('psi cn);
title_m = ['psi n = ', num2str(cn)] ;
title(title_m,'Fontsize',10);
subplot(nMax,2,2*cn);
y1 = prob(:,cn) ./ max(prob(:,cn));
y2 = 1 + 2 * U ./ (max(U) - min(U));
plot(x,y1,'lineWidth',2)
hold on
plot(x,y2,'r','lineWidth',1)
title_m = ['psi^2 n = ', num2str(cn)] ;
title(title_m,'Fontsize',10);
axis off
end
toc
m文件simpson1d.m
function integral = simpson1d(f,a,b)
% [1D] integration - Simpson's 1/3 rule
% f function a = lower bound b = upper bound
% Must have odd number of data points
% Simpson's coefficients 1 4 2 4 ... 2 4 1
numS = length(f); % number of data points
if mod(numS,2) == 1
sc = 2*ones(numS,1);
sc(2:2:numS-1) = 4;
sc(1) = 1;
sc(numS) = 1;
h = (b-a)/(numS-1);
integral = (h/3) * f * sc;
else
integral = 'Length of function must be an ODD number'
end
参考资料:
http://www.physics.usyd.edu.au/teach_res/mp/mphome.htm