含时势能情形的波函数时间演化示例

含时势能的部分代码。

L = 6.25; % Interval Length

N = 1000; % 设置 N*N 格点

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

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

NT = 50000; % No. of time steps

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

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

I1 = speye(N,N);

% 动能部分

% 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];

UMmt = exp(-0.5*1i*dT*((K').^2)/2); % exp(-i*(-0.5*Laplacian)*dT/2), 处于动量空间

% 势能部分

% 势能不可变化太剧烈，除非修改后面计算时间演化的方法

V0 = 20;

Hat = @(y) ((atan(y*10)+pi/2)/pi).^2; % “平滑化”的阶梯函数

U0 = V0*(Hat(x-1.5)+Hat(-x-1.5)+Hat(x+0.3).*Hat(-x+0.3));% 静态的势能

e = ones(N,1); 

Lap1 = spdiags([e -2*e e e e],[-1 0 1 N-1 -N+1],N,N)/dx^2; % 二阶导（差分）算子

r0=1;

Ut = @(r)(r0*r)*Hat(-x); % 动态的势能函数

% UVh = exp(-1i*(U0+Ut(t*dT))*dT); 

H0 = -0.5*Lap1+spdiags(U0+Ut(-0.5),0,N,N); % 离散化的哈密顿算子

[V,D]=eigs(full(H0),1,'smallestabs'); % 求基态波函数

psi0 = V(:,1); %初态

psi0 = psi0/sqrt(psi0'*psi0); 

rate = 0.05; %速率

h = 1.0e-2;

for t = 1:NT*10

    %当前势能

    DelTl = min(-0.5+rate*t*dT,0.5);

    Unow = U0+Ut(DelTl);

    

    UVh = exp(-1i*Unow*dT); 

    

    % 时间演化，Strang方法

    psik = full(UMmt.*fftshift(fft(psi0)));

    psi0 = ifft(ifftshift(psik));

    psi0 = full(UVh.*psi0);

    

    psik = full(UMmt.*fftshift(fft(psi0)));

    psi0 = ifft(ifftshift(psik));    

    

    %计算密度、相位等

    rho = abs(psi0).^2;

    phase = angle(psi0).*heaviside(rho-0.001*h);

    % 按需加上其他东西

    

end