3维波函数相关代码
% 粒子质量 m 、普朗克常数 hbar都取为 1
% 位矢空间格点
L = 10; % 立方形位矢空间边长
N = 200; % 设置 N*N*N 格点
if N>300
error('huge N');
end
x = linspace(-L/2,L/2-L/N,N)'; % 1 维格点
dx = x(2) - x(1); % 格点最小间距
[X,Y,Z]=meshgrid(x); % 3 维位矢空间格点
% 时间数组
dT = 0.01;
TF = 50;
T = 0:dT:TF; % 离散化时间
NT = length(T);
%*****************************************
% 选择方向建立球坐标
% 定义旋转矩阵 in O(3)
% 使用 Euler 角
Alpha = pi/80;
Ry = [cos(Alpha),0,sin(Alpha);...
0,1,0;...
-sin(Alpha),0,cos(Alpha)];
Beta = pi/2-0.05*pi;
Rx = [1,0,0;...
0,cos(Beta),-sin(Beta);...
0,sin(Beta),cos(Beta)];
Gamma = -pi/80;
Rz = [cos(Gamma),-sin(Gamma),0;...
sin(Gamma),cos(Gamma),0;...
0,0,1];
Rt = Ry*Rx*Rz; % 完整的旋转矩阵
% 将球坐标 theta 为 0 的方向选择为 Rt^(-1)*[0;0;1]
% 用自定义函数得到3维位矢格点中各个点的 r,theta,phi
R = xyz2rho(Rt(1,1)*X+Rt(1,2)*Y+Rt(1,3)*Z,...
Rt(2,1)*X+Rt(2,2)*Y+Rt(2,3)*Z,...
Rt(3,1)*X+Rt(3,2)*Y+Rt(3,3)*Z);
Tt = xyz2theta(Rt(1,1)*X+Rt(1,2)*Y+Rt(1,3)*Z,...
Rt(2,1)*X+Rt(2,2)*Y+Rt(2,3)*Z,...
Rt(3,1)*X+Rt(3,2)*Y+Rt(3,3)*Z);
P = xyz2phi(Rt(1,1)*X+Rt(1,2)*Y+Rt(1,3)*Z,...
Rt(2,1)*X+Rt(2,2)*Y+Rt(2,3)*Z,...
Rt(3,1)*X+Rt(3,2)*Y+Rt(3,3)*Z);
% 设置势能
% 至少在波函数振幅显著非零的区域内，
% 让势能函数尽量平缓
% 下面仅仅是一个例子
Am = 50;
al = 2;
SphP =@(x) -Am./(exp(-al*x)+exp(al*x)); % 一个势阱函数
Up = SphP(R); % 一个球对称势能
G = 2000;
d = 0.28*L;
Hat = @(y) ((atan(y*6)+pi/2)/pi).^4; % 一个连续化的“阶梯”函数
BOX =@(x,y,z) G*(Hat(x-d)+ Hat(-x-d))+...
G*(Hat(y-d)+ Hat(-y-d))+...
G*(Hat(z-d)+ Hat(-z-d)); % 边长为 2*d 的立方形势阱
Ub = BOX(Rt(1,1)*X+Rt(1,2)*Y+Rt(1,3)*Z,...
Rt(2,1)*X+Rt(2,2)*Y+Rt(2,3)*Z,...
Rt(3,1)*X+Rt(3,2)*Y+Rt(3,3)*Z); % 将原本定义的势阱旋转 Rt^(-1)
Ut = Up+Ub;% 总势能
% 时间平移算符的势能部分
UVh = exp(-1i*Ut*dT/2); % exp(-i*U*dT/2),处于位矢表象
% 时间平移算符的动能部分
e1 = ones(N,1);
% for even N
K = [2*pi*(-N/2:-1)/L,2*pi*(0:N/2-1)/L]';
% % for odd N
% K = [2*pi*(-N/2+1/2:-1)/L,2*pi*(0:N/2-1/2)/L];
% Kx = zeros(N,N,N);Ky = zeros(N,N,N);Kz = zeros(N,N,N);
[Kx,Ky,Kz] = meshgrid(K);% 3维动量空间格点
UMmt = exp(-1i*dT*(Kx.^2+Ky.^2+Kz.^2)/2); % exp(-i*(-0.5*Laplacian)*dT), 处于动量表象
%*****************************************
% 设置初态波函数
% 有相当多种选择，下面给出两种
% 球对称势能的能量本征态
% 构造比较精细的1维位矢格点
N1 = 2000;
x1 = linspace(1e-2,L/2,N1)'; % 最好避开 0，具体范围和格点数可以视情况调整
l = 1; % 选择角动量量子数
m = 0; % 沿 Rt^(-1)*[0;0;1] 方向的角动量分量量子数
nl = 1;% 该角动量（l）下的第n个能级（按照能量从低到高），取1~5
% 计算径向波函数
[En,Vl] = RadEig(x1,SphP(x1),l); %可以用En观察能量本征值
% 注意 nl 对应于束缚态还是非束缚态
% 完整的波函数
psi3 = spline(x1,Vl(:,nl).*(x1).^(-1),R).*...% 插值得到径向波函数
YLM(R,Tt,P,l,m);% 球谐函数
% %高斯波包
% % 球坐标下的动量中心，
%% 注意极轴方向仍是[0;0;1]
% ko = 2;
% Th = 0.5*pi;
% Ph = 0.5*pi;
% kxo = ko*sin(Th)*cos(Ph);
% kyo = ko*sin(Th)*sin(Ph);
% kzo = ko*cos(Th);
% % 球坐标下的位矢中心
%% 注意极轴仍是[0;0;1]
% ao = 1;% 控制波包大小
% R0 = 3;
% Thr = 0.5*pi;pi/3;atan(sqrt(2));
% Phr = -pi/3;0;pi/4;
% xo = R0*sin(Thr)*cos(Phr);
% yo = R0*sin(Thr)*sin(Phr);
% zo = R0*cos(Thr);
% % 高斯波包波函数
% psi3 = exp(-((X-xo).^2+(Y-yo).^2+(Z-zo).^2)/(ao^2))...
% .*exp(1i*(X*kxo+Y*kyo+Z*kzo));
psi3 = psi3/sqrt(sum(sum(sum(abs(psi3).^2)))); %归一化
%*****************************************
% 时间演化并绘图，Strang 方法
% exp(-i*H*dT)=exp(-i*U*dT/2)*exp(-i*(-0.5*Laplacian)*dT)*exp(-i*U*dT/2)+O(dT^3)
%建立 figure 和 axes
fig=figure;
fig.Position = [0 0 2000 1000];
ax = axes(fig);
ax.Position = [0.2 0.05 0.7 0.9];
whitebg(fig,'black');
alphamap(ax,'default' );
Z0 = zeros(N,N);
Al = 5e4; % 控制颜色深度（其实是不透明度）的参数
scale = 1;% 控制球坐标轮廓的大小比例的参数
% 密度分布
rho3 = abs(psi3).^2;
rho2 = sum(rho3,3); % 沿 [0;0;1] 方向投影为二维
% 显示初态波函数
Pr = surf(ax,x,x,Z0,Z0,...
'FaceColor','interp','EdgeColor','none',...
'AlphaData',rho2*Al,'FaceAlpha','interp',...
'AlphaDataMapping','direct');
ax.XAxis.Visible = 'off';
ax.YAxis.Visible = 'off';
grid(ax,'off');
hold(ax,'on')
% 球坐标轮廓
theta = 0:0.01*pi:2*pi;
plot(ax,scale*4*(1/al)*cos(theta),...
scale*4*(1/al)*sin(theta),'y.','linewidth',0.01);
reqx = scale*4*(1/al)*cos(theta);
reqy = scale*4*(1/al)*sin(theta);
reqz = 0*reqx;
req = [reqx;reqy;reqz];
req = Rt'*req;
plot3(ax,req(1,:),req(2,:),req(3,:),'y--','linewidth',0.01);
va = Rt'*[0;0;1];
ls = -2*scale*(4/al):1e-1*scale*(4/al):2*scale*(4/al);
plot3(ax,va(1)*ls,va(2)*ls,va(3)*ls+2.1*scale*(4/al),...
'y--','linewidth',0.01);
plot3(ax,va(1)*[ls(1)/2,ls(end)/2],...
va(2)*[ls(1)/2,ls(end)/2],...
va(3)*[ls(1)/2,ls(end)/2]+2.1*scale*(4/al),...
'yo','linewidth',0.01);
plot3(ax,0,0,0+2.1*scale*(4/al),'y*','linewidth',0.01);
% 绘图的一些参数
pbaspect([1 1 1]);
campos([0 0 L])
camtarget([0 0 0])
axis(ax,[-0.5*L 0.5*L -0.5*L 0.5*L]);
drawnow
delete(Pr);
% 时间演化
for t = 1:NT
% exp(-i*U*dT/2)
psi3 = full(UVh.*psi3);
%exp(-i*(-0.5*Laplacian)*dT)
psik = full(UMmt.*fftshift(fftn(psi3)));
% exp(-i*U*dT/2)
psi3 = full(UVh.*ifftn(ifftshift(psik)));
% 密度分布
rho3 = (abs(psi3)).^2;
rho2 = sum(rho3,3);% 投影至二维
% 绘图
Pr = surf(ax,x,x,Z0,Z0,'FaceColor','interp','EdgeColor','none',...
'AlphaData',rho2*Al,'FaceAlpha','interp','AlphaDataMapping','direct');
TxtT = text(ax,-1.1*L,0.2*L,...
['$T=',num2str(dT*t,'%.3f'),'$'],...
'Interpreter','latex','FontSize',30);% 显示时间
axis(ax,[-0.5*L 0.5*L -0.5*L 0.5*L]);
pbaspect([1 1 1]);
campos([0 0 L])
camtarget([0 0 0])
drawnow
delete(Pr);
delete(TxtT);
end
hold(ax,'off')
%*****************************************
%自定义函数
function Rho = xyz2rho(x,y,z) % Cartesian 坐标 to 球坐标的 r
Rho = sqrt(x.^2+y.^2+z.^2);
end
function th = xyz2theta(x,y,z)% 球坐标的theta
r = sqrt(x.^2+y.^2+z.^2);
c = z./r;
c(isinf(c))=0;
c(isnan(c))=0;
th = acos(c);
end
function ph = xyz2phi(x,y,z) % 球坐标的phi
ph = atan2(y,x);
end
function [En,V] = RadEig(x1,U1,l)
% 计算球对称情形的等效 1 维本征函数，l 为角动量量子数
% 要求格点 x1 间隔均匀，间隔不能过小，但数组长度不宜过大
% x1 中元素应当都是正实数，最小值接近 0
N1 = length(x1);
dx1 = x1(2)-x1(1);
e = ones(N1,1);
% 离散化二阶导数，硬边界
Lap1 = spdiags([e -2*e e],[-1 0 1],N1,N1)/dx1^2;
% 有效势能
xp2 = x1.^(-2);
% xp2(isinf(xp2)) = dx1^(-2);% 截断无穷大值
U1 = U1+0.5*l*(l+1).*xp2;
% 1维 Hamiltonian
H0 = -0.5*Lap1+spdiags(U1+100,0,N1,N1);
% 计算 |En+100| 最小的 5 个本征矢和相应的 本征值+100
% 因此基态能量最好大于 -100。可自行修改具体数值
[V,D]=eigs(full(H0),5,'smallestabs');
% 本征能量
En = diag(D)-100;
end
function Ylm = YLM(R,T,P,l,m) % 球谐函数 Ylm(r,theta,phi)
Nx = size(R);
Nx = Nx(1);
C1 = sqrt(((2*l+1)*factorial(l-abs(m)))/(4*pi*factorial(l+abs(m))));
lgd = legendre(l,cos(T));
if l
Ylm = C1*reshape(lgd(abs(m)+1,:,:,:),Nx,Nx,Nx)...
.*exp(1i*m*P);
else
Ylm = C1*lgd...
.*exp(1i*m*P);
end
end
