3維波函數(shù)相關(guān)

% 粒子質(zhì)量 m 、普朗克常數(shù) hbar都取為 1

% 位矢空間格點(diǎn)

L = 10; % 立方形位矢空間邊長(zhǎng)

N = 200; % 設(shè)置 N*N*N 格點(diǎn)

if N>300

? ? error('huge N');

end

x = linspace(-L/2,L/2-L/N,N)'; % 1 維格點(diǎn)

dx = x(2) - x(1); % 格點(diǎn)最小間距

[X,Y,Z]=meshgrid(x); % 3 維位矢空間格點(diǎn)

% 時(shí)間數(shù)組

dT = 0.01;

TF = 50;?

T = 0:dT:TF; % 離散化時(shí)間

NT = length(T);

%*****************************************

% 選擇方向建立球坐標(biāo)

% 定義旋轉(zhuǎn)矩陣 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; % 完整的旋轉(zhuǎn)矩陣

% 將球坐標(biāo) theta 為 0 的方向選擇為 Rt^(-1)*[0;0;1]

% 用自定義函數(shù)得到3維位矢格點(diǎn)中各個(gè)點(diǎn)的 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);

% 設(shè)置勢(shì)能

% 至少在波函數(shù)振幅顯著非零的區(qū)域內(nèi)，

% 讓勢(shì)能函數(shù)盡量平緩

% 下面僅僅是一個(gè)例子

Am = 50;

al = 2;

SphP =@(x) -Am./(exp(-al*x)+exp(al*x)); % 一個(gè)勢(shì)阱函數(shù)

Up = SphP(R); % 一個(gè)球?qū)ΨQ勢(shì)能

G = 2000;

d = 0.28*L;

Hat = @(y) ((atan(y*6)+pi/2)/pi).^4; % 一個(gè)連續(xù)化的“階梯”函數(shù)

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)); % 邊長(zhǎng)為 2*d 的立方形勢(shì)阱

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); % 將原本定義的勢(shì)阱旋轉(zhuǎn) Rt^(-1)

Ut = Up+Ub;% 總勢(shì)能

% 時(shí)間平移算符的勢(shì)能部分

UVh = exp(-1i*Ut*dT/2); % exp(-i*U*dT/2),處于位矢表象

% 時(shí)間平移算符的動(dòng)能部分

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維動(dòng)量空間格點(diǎn)

UMmt = exp(-1i*dT*(Kx.^2+Ky.^2+Kz.^2)/2); % exp(-i*(-0.5*Laplacian)*dT), 處于動(dòng)量表象

%*****************************************

% 設(shè)置初態(tài)波函數(shù)

% 有相當(dāng)多種選擇，下面給出兩種

% 球?qū)ΨQ勢(shì)能的能量本征態(tài)

% 構(gòu)造比較精細(xì)的1維位矢格點(diǎn)

N1 = 2000;?

x1 = linspace(1e-2,L/2,N1)'; % 最好避開(kāi) 0，具體范圍和格點(diǎn)數(shù)可以視情況調(diào)整

l = 1; % 選擇角動(dòng)量量子數(shù)

m = 0; % 沿 Rt^(-1)*[0;0;1] 方向的角動(dòng)量分量量子數(shù)

nl = 1;% 該角動(dòng)量（l）下的第n個(gè)能級(jí)（按照能量從低到高），取1~5

% 計(jì)算徑向波函數(shù)

[En,Vl] = RadEig(x1,SphP(x1),l); %可以用En觀察能量本征值

% 注意 nl 對(duì)應(yīng)于束縛態(tài)還是非束縛態(tài)

% 完整的波函數(shù)

psi3 = spline(x1,Vl(:,nl).*(x1).^(-1),R).*...% 插值得到徑向波函數(shù)

? ? YLM(R,Tt,P,l,m);% 球諧函數(shù)

% %高斯波包

% % 球坐標(biāo)下的動(dòng)量中心，

%% 注意極軸方向仍是[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);

% % 球坐標(biāo)下的位矢中心

%% 注意極軸仍是[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);

% % 高斯波包波函數(shù)

% 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)))); %歸一化

%*****************************************

% 時(shí)間演化并繪圖，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; % 控制顏色深度（其實(shí)是不透明度）的參數(shù)

scale = 1;% 控制球坐標(biāo)輪廓的大小比例的參數(shù)

% 密度分布

rho3 = abs(psi3).^2;

rho2 = sum(rho3,3); % 沿 [0;0;1] 方向投影為二維

% 顯示初態(tài)波函數(shù)

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')

% 球坐標(biāo)輪廓

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);

% 繪圖的一些參數(shù)

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);

% 時(shí)間演化

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);% 顯示時(shí)間

? ??

? ??

? ? 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')

%*****************************************

%自定義函數(shù)

function Rho = xyz2rho(x,y,z) % Cartesian 坐標(biāo) to 球坐標(biāo)的 r

Rho = sqrt(x.^2+y.^2+z.^2);

end

function th = xyz2theta(x,y,z)% 球坐標(biāo)的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) % 球坐標(biāo)的phi

ph = atan2(y,x);

end

function [En,V] = RadEig(x1,U1,l)?

% 計(jì)算球?qū)ΨQ情形的等效 1 維本征函數(shù)，l 為角動(dòng)量量子數(shù)

% 要求格點(diǎn) x1 間隔均勻，間隔不能過(guò)小，但數(shù)組長(zhǎng)度不宜過(guò)大

% x1 中元素應(yīng)當(dāng)都是正實(shí)數(shù)，最小值接近 0

N1 = length(x1);

dx1 = x1(2)-x1(1);

e = ones(N1,1);?

% 離散化二階導(dǎo)數(shù)，硬邊界

Lap1 = spdiags([e -2*e e],[-1 0 1],N1,N1)/dx1^2;?

% 有效勢(shì)能

xp2 = x1.^(-2);

% xp2(isinf(xp2)) = dx1^(-2);% 截?cái)酂o(wú)窮大值

U1 = U1+0.5*l*(l+1).*xp2;

% 1維 Hamiltonian

H0 = -0.5*Lap1+spdiags(U1+100,0,N1,N1);

% 計(jì)算 |En+100| 最小的 5 個(gè)本征矢和相應(yīng)的 本征值+100

% 因此基態(tài)能量最好大于 -100?？勺孕行薷木唧w數(shù)值

[V,D]=eigs(full(H0),5,'smallestabs');

% 本征能量

En = diag(D)-100;

end

function Ylm = YLM(R,T,P,l,m) % 球諧函數(shù) 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