粒子在2维势能中的波函数演化模拟
2维势能中的波包演化相关的matlab代码。
渣代码，轻喷。
哈密顿量和时间演化相关部分
L = 20; % Interval Length
N = 500; % 设置 N*N 格点. 有大动量时格点不能太稀疏
% x = linspace(-L/2,L/2-L/(N),N)'; % Coordinate vector， 周期性边界时
x = linspace(-L/2,L/2,N)'; % Coordinate vector，硬边界时
dx = x(2) - x(1); % Coordinate step
% operators in 1D:
e = ones(N,1);
% Dif1 = spdiags([-e 0*e e -e e],[-1 0 1 N-1 -N+1],N,N)/(2*dx);% 有周期性边界
% Lap1 = spdiags([e -2*e e e e],[-1 0 1 N-1 -N+1],N,N)/dx^2;
Dif1 = spdiags([-e 0*e e],[-1 0 1],N,N)/(2*dx);% 硬边界
Lap1 = spdiags([e -2*e e],[-1 0 1],N,N)/dx^2;
I1 = speye(N,N);
% operators in 2D:
Lap2 = kron(Lap1,I1)+kron(I1,Lap1); % 2D Laplacian
% PDx = kron(Dif1,I1); % partial diff of x
% PDy = kron(I1,Dif1);% partial diff of y
% xy = spdiags(x,0,N,N);
% X = kron(xy,I1); % x position operator
% Y = kron(I1,xy); % y position operator
% R2 = kron(xy.^2,I1)+kron(I1,xy.^2); % R^2=x^2+y^2
% 设置势能
% 可按需调整
V=100000;
a = 0.3;
U = V*kron(...
spdiags((x.^2).*exp(-(x).^2/a^2),0,N,N),... % f(x)
spdiags((x.^2).*exp(-(x).^2/a^2),0,N,N)... %虽然记号是x，实际上应当被理解为 g(y)
); % 得到 V*f(x)*g(y)形式的势能。也可以多取几个这样的势能在相加
% figure
% surf(x,x,reshape(diag(U),N,N)); % 查看势能
% 边界附近的虚势能，使得运动到边界的波衰减
% 不一定必要
% 不一定对任意动量都适用，有时需要调整参数
G = 20;
l = 1.2;
Ui = -1i*G*(...
kron(...
spdiags(exp(-(x+0.5*L)/l)+exp(-(-x+0.5*L)/l),0,N,N),I1)+...
kron(...
I1,spdiags(exp(-(x+0.5*L)/l)+exp(-(-x+0.5*L)/l),0,N,N))...
);
% figure
% surf(x,x,reshape(diag(imag(Ui)),N,N)); % 查看势能虚部
H = 0.5*(-Lap2)+U+Ui; % Hamiltonian
NT = 50000; % No. of time steps
TF = 50; T = linspace(0,TF,NT); % Time vector
dT = T(2)-T(1); % Time step
% E = expm(-1i*H*dT); % Evolution
I = speye(size(H));
A = -1i*H*dT;
E = I;
for m=1:10 % 可调整展开阶数
E = E+(A^m/factorial(m));
end
% spy(E)
（沿x向传播的）高斯波包：
e1 = e;
xf = kron(x,e1);
yf = kron(e1,x);
xo = -4.2; yo = 0; % center position
k0 = 10; a1 = 2; % 可手动调整
psi0 = exp(-((xf-xo).^2+(yf-yo).^2)/(a1^2)).*exp(1i*xf*k0);
psi0 = psi0/sqrt(psi0'*psi0); %归一化
画图：
fig=figure;
fig.Position = [0 0 2000 1000];
ax = axes(fig);
h=0.2e-3; % 依据波函数模平方的大致范围而定
xlim manual
ylim manual
zlim manual
% 波函数模平方的颜色初始设置
CO=zeros(N,N,3);
CO(:,:,1) = 0.25*ones(N,N);
CO(:,:,3) = 0.8*ones(N,N);
Nb = 1000; % 给相位插值时使用（让相位看起来更平滑）
Z0 = -1.5*h*ones(Nb,Nb); % 显示相位的平面高度
for t = 1:NT % time index for loop
% calculate probability density rho(x,T)
psi0 = E*psi0; % calculate new psi from old psi
% if rem(t,4) %跳过部分帧
% continue
% end
% 计算psi模平方
rho = abs(psi0).^2;
rho1 = reshape(rho,[N,N]);
%设置着色
CO(:,:,2) = 15*ones(N,N).*(rho1)/(h);
CO(:,:,3) = 0.8*ones(N,N)+5*ones(N,N).*atan(rho1)/h;
% 绘制波函数模平方的曲面，按照数组CO着色
s0 = surf(ax,x,x,rho1+0.1*h,CO,'FaceAlpha',0.4);
shading interp
s0.EdgeColor = 'none';
text(ax,-7,7,1.2*h,...
['$T=',num2str(dT*t,'%.2f'),'$'],...
'Interpreter','latex','FontSize',20); %时间
text(ax,-7,7,1*h,...
['$|\langle\mathbf{k}_0\rangle|=',num2str(abs(k0),'%.3f'),'$'],...
'Interpreter','latex','FontSize',20); % 初态的平均动量的大小
% 固定坐标轴范围
zlim(ax,[-1.5*h 1.5*h]);
xlim(ax,[-0.5*L 0.5*L]);
ylim(ax,[-0.5*L 0.5*L]);
hold on
% 在xy平面上用颜色显示势能的大致形貌
% 细节参数可按需调整
imagesc(ax,[-0.5*L,0.5*L],[-0.5*L,0.5*L],...
reshape(diag(real(10*U/(V))),N,N),'AlphaData',0.7)
%计算相位，以 0.001*h 为波函数模平方的阈值，低于阈值的部分不保留结果
phase = reshape(...
h*angle(psi0).*heaviside(abs(psi0).^2-0.001*h)...
,[N,N]);
%给相位插值，不一定必要
[Xq,Yq] = meshgrid(linspace(x(1),x(end),Nb));
Phase = interp2(kron(x',e1),kron(e1',x),phase,Xq,Yq,'makima');
% 在相应平面用颜色绘制相位
Ph = surf(ax,Xq,Yq,Z0,Phase,'FaceAlpha',0.3);
shading interp
Ph.EdgeColor = 'none';
hold off
% 颜色设置
colormap(hsv)
caxis([-pi*h pi*h])
% 3维绘图的参数设置
pbaspect([1 1 0.5]);
campos([-18 -35 2.55*h])
camtarget([0 0 0*h])
camva(20)
drawnow
end
