部分 matlab 代碼。

% 設(shè)置格點(diǎn)化的空間坐標(biāo)

L = 100; % Interval Length

N = 5000; % No of points

x = linspace(-L/2,L/2,N)'; % Coordinate vector

dx = x(2) - x(1); % Coordinate step

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

%設(shè)置 Hamiltonian

V=83; %barrier height

U = V.*heaviside(x-1).*heaviside(-x+2)+ V.*heaviside(-x-1).*heaviside(x+2); % scattering off potential

G = 20;%控制邊界上的衰減

U = U-heaviside(x-45)*(1i*G).*(x-45).*(x-45)-heaviside(-x-45)*(1i*G).*(x+45).*(x+45); %邊界附近的虛勢(shì)能，把跑到邊界的粒子移出系統(tǒng)

% Finite-difference representation of Laplacian and Hamiltonian

e = ones(N,1);

Lap = spdiags([e -2*e e e e],[-1 0 1 N-1 -N+1],N,N)/dx^2; % 周期性邊界條件下的離散Laplacian?

H = -(1/2)*Lap + spdiags(U,0,N,N); % Hamiltonian

% Parameters for computing psi(x,T) at different times 0 < T < TF

NT = 4000; % No. of time steps

TF = 150; T = linspace(0,TF,NT); % Time vector

dT = T(2)-T(1); % Time step

hbar = 1;

% Time displacement operator E=exp(-iHdT/hbar)

E = expm(-1i*full(H)*dT/hbar); % time desplacement operator

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

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

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

% 下面部分可以和上面部分分開運(yùn)行。

% 僅修改波函數(shù)以及畫圖的相關(guān)參數(shù)時(shí)，不必重復(fù)運(yùn)行上面的部分。

% Parameters for making intial wavefunction psi

ko1 = 12.50;?

ko2 = 12.73;

ko3 = 12.96;

%ko = 12.7292206; % Peak momentum

a = 7; % package width

xo=-22;% peak center

% 離散化的空間坐標(biāo)對(duì)最大動(dòng)量有限制，不然會(huì)出現(xiàn)Bloch振蕩這類現(xiàn)象

if ko1>(pi*N/L)

? ? error('too large ko1');

end

if ko2>(pi*N/L)

? ? error('too large ko2');

end

if ko3>(pi*N/L)

? ? error('too large ko3');

end

% 設(shè)置波函數(shù)為 平面波*Gauss輪廓，及其歸一化

psi1 = exp(-(x-xo).^2/(2*a.^2)).*exp(1i*x*ko1);

psi1 = psi1/sqrt(psi1'*psi1*dx); % normalize the psi(x,0)

psi2 = exp(-(x-xo).^2/(2*a.^2)).*exp(1i*x*ko2);

psi2 = psi2/sqrt(psi2'*psi2*dx); % normalize the psi(x,0)

psi3 = exp(-(x-xo).^2/(2*a.^2)).*exp(1i*x*ko3);

psi3 = psi3/sqrt(psi3'*psi3*dx); % normalize the psi(x,0)

% 畫圖的亂七八糟的參數(shù)

% y用于設(shè)置fft圖的橫坐標(biāo)

D=floor(20*(L/2-2)/(2*pi));

y = 2*pi*(0:length(floor(N/2+N*(2/L)):N)-1)/(L/2-2);

N1=length(y);

H=0.1;

fig=figure;

fig.Position = [0 0 2000 1000];

dy = y(2)-y(1);

for t = 1:NT? % time index for loop

% calculate probability density rho(x,T)

psi1 = E*psi1; % calculate new psi from old psi

rho1 = conj(psi1).*psi1; % rho(x,T)

psi2 = E*psi2;

rho2 = conj(psi2).*psi2;?

psi3 = E*psi3;

rho3 = conj(psi3).*psi3;?

% 右側(cè)（透射）區(qū)域的動(dòng)量密度分布

f1t = fft(psi1(N-N1+1:N))/sqrt(N1);

f2t = fft(psi2(N-N1+1:N))/sqrt(N1);

f3t = fft(psi3(N-N1+1:N))/sqrt(N1);

% 左側(cè)（入/反射）區(qū)域的動(dòng)量密度分布

f1i = fft(psi1(1:N1))/sqrt(N1);

f2i = fft(psi2(1:N1))/sqrt(N1);

f3i = fft(psi3(1:N1))/sqrt(N1);

% 畫圖

subplot(223);

plot([ko1,ko1],[-9,5],'color',[0,0.4,0.5],'LineStyle','--');

hold on

% 對(duì) fft 的結(jié)果重排，把后半部分挪到 k<0 的位置（類似于把 0~2π 改為 -π）

plot(...

? ? [y(N1-D:N1)-y(N1)-dy,y],log([abs(f1t(N1-D:N1)).^2;abs(f1t).^2]),'c',...

? ? [y(N1-D:N1)-y(N1)-dy,y],log([abs(f2t(N1-D:N1)).^2;abs(f2t).^2]),'k',...

? ? [y(N1-D:N1)-y(N1)-dy,y],log([abs(f3t(N1-D:N1)).^2;abs(f3t).^2]),'r');

axis([-17 17 -9 5]);?

hold off

xlabel('k_t', 'FontSize', 24);

subplot(222);

plot([ko1,ko1],[-9,5],'color',[0,0.4,0.5],'LineStyle','--');

hold on

plot([-ko1,-ko1],[-9,5],'color',[0,0.4,0.5],'LineStyle','--');

plot(...

? ? [y(N1-D:N1)-y(N1)-dy,y],log([abs(f1i(N1-D:N1)).^2;abs(f1i).^2]),'c',...

? ? [y(N1-D:N1)-y(N1)-dy,y],log([abs(f2i(N1-D:N1)).^2;abs(f2i).^2]),'k',...

? ? [y(N1-D:N1)-y(N1)-dy,y],log([abs(f3i(N1-D:N1)).^2;abs(f3i).^2]),'r');

hold off

axis([-17 17 -9 5]);?

xlabel('k_{ir}', 'FontSize', 24);

subplot(221);

plot(x,rho1,'c',x,real(psi1).^2,'g--',...

? ? x,rho2,'k',x,real(psi2).^2,'b--',...

? ? x,rho3,'r',x,real(psi3).^2,'m--'); % plot rho(x,T) vs. x

hold on?

plot([1,1],[0,H],'color',[0,0.5,0.4],'LineStyle','--');

plot(-[1,1],[0,H],'color',[0,0.5,0.4],'LineStyle','--');

plot([2,2],[0,H],'color',[0,0.5,0.4],'LineStyle','--');

plot(-[2,2],[0,H],'color',[0,0.5,0.4],'LineStyle','--');

hold off

axis([-L/2 L/2 0 H]); % set x,y axis parameters for plotting

xlabel('x', 'FontSize', 24);

subplot(224);

plot(x,rho1,'c',x,real(psi1).^2,'g--',...

? ? x,rho2,'k',x,real(psi2).^2,'b--',...

? ? x,rho3,'r',x,real(psi3).^2,'m--'); % plot rho(x,T) vs. x

hold on?

plot([1,1],[0,H],'color',[0,0.5,0.4],'LineStyle','--');

plot(-[1,1],[0,H],'color',[0,0.5,0.4],'LineStyle','--');

plot([2,2],[0,H],'color',[0,0.5,0.4],'LineStyle','--');

plot(-[2,2],[0,H],'color',[0,0.5,0.4],'LineStyle','--');

hold off

axis([-3 3 0 H/5]); % set x,y axis parameters for plotting

xlabel('x', 'FontSize', 24);

drawnow

end