模拟一维波包散射的部分代码
部分 matlab 代码。
% 设置格点化的空间坐标
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
%***************************************************************
%设置 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); %边界附近的虚势能，把跑到边界的粒子移出系统
% 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
%***************************************************************
%***************************************************************
%***************************************************************
% 下面部分可以和上面部分分开运行。
% 仅修改波函数以及画图的相关参数时，不必重复运行上面的部分。
% 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
% 离散化的空间坐标对最大动量有限制，不然会出现Bloch振荡这类现象
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
% 设置波函数为 平面波*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)
% 画图的乱七八糟的参数
% y用于设置fft图的横坐标
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;
% 右侧（透射）区域的动量密度分布
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);
% 左侧（入/反射）区域的动量密度分布
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
% 对 fft 的结果重排，把后半部分挪到 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
