✅作者简介：热爱科研的Matlab仿真开发者，修心和技术同步精进，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 numberNGy = 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 Octaveidx_X = (repmat(idx_x,[1 NG])-repmat(idx_x',[NG 1])) + NGx; %% work in Octave and Matlabidx_Y = (repmat(idx_y,[1 NG])-repmat(idx_y',[NG 1])) + NGy; %% work in Octave and Matlabidx = 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 psiend%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% here is a small patch due to differences between Octave and Matlab% Matlab order the eigen values while Octave reverse itif E(1)>E(2) psi=psi(:,:,end:-1:1); E=E(end:-1:1);endend%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%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.
