✅作者简介：热爱科研的Matlab仿真开发者，修心和技术同步精进，matlab项目合作可私信。
🍎个人主页：Matlab科研工作室
🍊个人信条：格物致知。
⛄ 内容介绍
在量子力学中，薛定谔方程是研究粒子行为的基本方程之一。它描述了粒子的波函数随时间的演化，并提供了关于粒子在不同位置和动量上的概率分布的信息。在本文中，我们将讨论如何求解具有恒定质量的二维薛定谔方程。
二维薛定谔方程的一般形式如下：
iħ∂ψ/∂t = -ħ²/2m(∂²ψ/∂x² + ∂²ψ/∂y²) + V(x, y)ψ
其中，i是虚数单位，ħ是约化普朗克常数，t是时间，m是粒子的质量，(x, y)是粒子的位置，V(x, y)是势能函数，ψ是波函数。
为了求解这个方程，我们可以采用分离变量的方法。假设波函数可以表示为两个单变量函数的乘积形式，即ψ(x, y) = X(x)Y(y)。将这个形式代入方程中，并将方程两边除以ψ(x, y)，我们可以得到：
(iħ/X)dX/dx + (iħ/Y)dY/dy = -ħ²/2m(X''/X + Y''/Y) + V(x, y)
由于等式两边只依赖于不同的变量，所以它们必须等于一个常数，我们将其记为E。这样我们就可以将方程分解为两个单变量方程：
(iħ/X)dX/dx + V(x)X = EX
(iħ/Y)dY/dy + V(y)Y = EY
解这两个方程，我们可以得到X(x)和Y(y)的解。然后，我们可以将它们乘在一起，得到波函数的解ψ(x, y)。
值得注意的是，势能函数V(x, y)在不同的问题中会有不同的形式。在一些简单的情况下，它可能是一个常数，表示一个均匀势场。在其他情况下，它可能是一个依赖于位置的函数，表示一个非均匀势场。根据具体问题的不同，我们需要选择适当的势能函数来求解方程。
一旦我们找到了波函数的解，我们就可以计算粒子在不同位置上的概率分布。根据量子力学的原理，波函数的模的平方给出了粒子在给定位置上的概率密度。通过对波函数的模的平方进行归一化，我们可以得到粒子在整个空间中的总概率为1。
总结起来，求解具有恒定质量的二维薛定谔方程是一个重要的问题，在量子力学研究中具有广泛的应用。通过采用分离变量的方法，我们可以将方程分解为两个单变量方程，并求解它们以得到波函数的解。这样我们就可以计算粒子在不同位置上的概率分布，从而了解粒子的行为。
希望本文对读者理解求解具有恒定质量的二维薛定谔方程有所帮助。在实际应用中，我们可能会遇到更加复杂的情况，需要采用其他方法来求解方程。然而，分离变量法是一个重要的起点，它为我们理解量子力学中的基本概念和现象提供了一个坚实的基础。
⛄ 核心代码
function[E,psi]=Schroed2D_PWE_f(x,y,V0,Mass,n,Nx,Ny,NGx,NGy)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Constants %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
h=6.62606896E-34; %% Planck constant [J.s]
hbar=h/(2*pi);
e=1.602176487E-19; %% electron charge [C]
m0=9.10938188E-31; %% electron mass [kg]
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%% Interpolation on a grid that have 2^N points %%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
NGx = 2*floor(NGx/2); %% round to lower even number
NGy = 2*floor(NGy/2); %% round to lower even number
[X,Y] = meshgrid(x,y);
xx=linspace(x(1),x(end),Nx);
yy=linspace(y(1),y(end),Ny);
[XX,YY] = meshgrid(xx,yy);
V=interp2(X,Y,V0,XX,YY);
dx=x(2)-x(1);
dxx=xx(2)-xx(1);
dy=y(2)-y(1);
dyy=yy(2)-yy(1);
Ltotx=xx(end)-xx(1);
Ltoty=yy(end)-yy(1);
[XX,YY] = meshgrid(xx,yy);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%% Building of the potential in Fourier space %%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Vk = fftshift(fft2(V))*dxx*dyy/Ltotx/Ltoty;
Vk =Vk(Ny/2-NGy+1:Ny/2+NGy+1 , Nx/2-NGx+1:Nx/2+NGx+1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%% Reciprocal lattice vectors %%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Gx = (-NGx/2:NGx/2)'*2*pi/Ltotx;
Gy = (-NGy/2:NGy/2)'*2*pi/Ltoty;
NGx=length(Gx);
NGy=length(Gy);
NG=NGx*NGy;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%% Building Hamiltonien %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
idx_x = repmat((1:NGx), [NGy 1 ]);
idx_x = idx_x(:);
idx_y = repmat((1:NGy)', [1 NGx]);
idx_y = idx_y(:);
%idx_X = (idx_x-idx_x') + NGx; %% work only in Octave
%idx_Y = (idx_y-idx_y') + NGy; %% work only in Octave
idx_X = (repmat(idx_x,[1 NG])-repmat(idx_x',[NG 1])) + NGx; %% work in Octave and Matlab
idx_Y = (repmat(idx_y,[1 NG])-repmat(idx_y',[NG 1])) + NGy; %% work in Octave and Matlab
idx = sub2ind(size(Vk), idx_Y(:), idx_X(:));
idx = reshape(idx, [NG NG]);
GX = diag(Gx(idx_x));
GY = diag(Gy(idx_y));
D2 = GX.^2 + GY.^2 ;
H = hbar^2/(2*m0*Mass)*D2 + Vk(idx)*e ;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%% Solving Hamiltonien %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
H = sparse(H);
[psik, Ek] = eigs(H,n,'SM');
E = diag(Ek) / e;
%E=abs(E);
E=real(E);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%% Transforming & Scaling the waves functions %%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for i=1:n
PSI = reshape(psik(:,i),[NGy,NGx]);
PSI = invFFT2D(PSI,Ny,Nx)/(dxx*dyy) ;
psi_temp = interp2(XX,YY,PSI,X,Y);
psi(:,:,i) = psi_temp / sqrt( trapz( y' , trapz(x,abs(psi_temp).^2 ,2) , 1 ) ); % normalisation of the wave function psi
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% here is a small patch due to differences between Octave and Matlab
% Matlab order the eigen values while Octave reverse it
if E(1)>E(2)
psi=psi(:,:,end:-1:1);
E=E(end:-1:1);
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [Vxy] = invFFT2D(Vk2D,Ny,Nx)
Nkx=length(Vk2D(1,:));
Nky=length(Vk2D(:,1));
Nx1=Nx/2-floor(Nkx/2);
Nx2=Nx/2+ceil(Nkx/2);
Ny1=Ny/2-floor(Nky/2);
Ny2=Ny/2+ceil(Nky/2);
Vk2D00=zeros(Ny,Nx);
Vk2D00( Ny1+1:Ny2 , Nx1+1:Nx2)=Vk2D;
Vxy=ifft2(ifftshift(Vk2D00));
end⛄ 运行结果
⛄ 参考文献
[1]刘晓军.有限差分法解薛定谔方程与MATLAB实现[J].高师理科学刊, 2010(3):3.DOI:10.3969/j.issn.1007-9831.2010.03.022.
