2維勢能中的波包演化相關(guān)的matlab代碼。

渣代碼，輕噴。

哈密頓量和時(shí)間演化相關(guān)部分

L = 20; % Interval Length

N = 500; % 設(shè)置 N*N 格點(diǎn). 有大動(dòng)量時(shí)格點(diǎn)不能太稀疏

% x = linspace(-L/2,L/2-L/(N),N)'; % Coordinate vector， 周期性邊界時(shí)

x = linspace(-L/2,L/2,N)'; % Coordinate vector，硬邊界時(shí)

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

% 設(shè)置勢能

% 可按需調(diào)整

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)... %雖然記號(hào)是x，實(shí)際上應(yīng)當(dāng)被理解為 g(y)

); % 得到 V*f(x)*g(y)形式的勢能。也可以多取幾個(gè)這樣的勢能在相加

% figure

% surf(x,x,reshape(diag(U),N,N)); % 查看勢能

% 邊界附近的虛勢能，使得運(yùn)動(dòng)到邊界的波衰減

% 不一定必要

% 不一定對任意動(dòng)量都適用，有時(shí)需要調(diào)整參數(shù)

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 % 可調(diào)整展開階數(shù)

? ? 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; % 可手動(dòng)調(diào)整

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; % 依據(jù)波函數(shù)模平方的大致范圍而定

xlim manual

ylim manual

zlim manual

% 波函數(shù)模平方的顏色初始設(shè)置

CO=zeros(N,N,3);

CO(:,:,1) = 0.25*ones(N,N);

CO(:,:,3) = 0.8*ones(N,N);

Nb = 1000; % 給相位插值時(shí)使用（讓相位看起來更平滑）

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

% 計(jì)算psi模平方

? ? rho = abs(psi0).^2;

? ? rho1 = reshape(rho,[N,N]);

%設(shè)置著色

? ? CO(:,:,2) = 15*ones(N,N).*(rho1)/(h);

? ? CO(:,:,3) = 0.8*ones(N,N)+5*ones(N,N).*atan(rho1)/h;

% 繪制波函數(shù)模平方的曲面，按照數(shù)組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); %時(shí)間

? ? text(ax,-7,7,1*h,...

? ? ? ? ['$|\langle\mathbf{k}_0\rangle|=',num2str(abs(k0),'%.3f'),'$'],...

? ? ? ? 'Interpreter','latex','FontSize',20); % 初態(tài)的平均動(dòng)量的大小

? ? % 固定坐標(biāo)軸范圍

? ? 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平面上用顏色顯示勢能的大致形貌

? ? % 細(xì)節(jié)參數(shù)可按需調(diào)整

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

? ? %計(jì)算相位，以?0.001*h?為波函數(shù)模平方的閾值，低于閾值的部分不保留結(jié)果

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

? ? % 在相應(yīng)平面用顏色繪制相位

? ? Ph = surf(ax,Xq,Yq,Z0,Phase,'FaceAlpha',0.3);

? ? shading interp

? ? Ph.EdgeColor = 'none';

? ??

? ? hold off

? ? % 顏色設(shè)置

? ? colormap(hsv)

? ? caxis([-pi*h pi*h])

? ? % 3維繪圖的參數(shù)設(shè)置

? ? pbaspect([1 1 0.5]);

? ? campos([-18 -35 2.55*h])

? ? camtarget([0 0 0*h])

? ? camva(20)

? ??

? ? drawnow

end