含時勢能的部分代碼。

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; % “平滑化”的階梯函數(shù)

U0 = V0*(Hat(x-1.5)+Hat(-x-1.5)+Hat(x+0.3).*Hat(-x+0.3));% 靜態(tài)的勢能

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); % 動態(tài)的勢能函數(shù)

% 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'); % 求基態(tài)波函數(shù)

psi0 = V(:,1); %初態(tài)

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