贡献者: addis
我们用式 8 计算波函数的含时演化。另外我们用离散的 $k_i$ 代替连续的 $k$,对 $k$ 积分改为对 $i$ 求和并乘以 $\Delta{k}$。另一种可能更常见的数值解 TDSE 算法见 “Crank-Nicolson 算法解一维含时薛定谔方程(Matlab)”。
代码中使用了 “方势垒定态波函数程序” 中的 FSB.m
和 FSB2.m
。
% 方势垒散射高斯波包
clear; close all;
% ==== 参数 ====
% 高斯波包
x0 = -17; t0 = 0; m = 1;
sig_x = 2; k0 = 4;
% 方势垒 [-L,L], 高 V0
L = 4; V0 = 4;
mode = 1; % 1: 使用类 sin,cos 基底; 2: 使用类 exp 基底
% 位置、动量、时间网格
xmin = -50; xmax = 50; Nx = 2000;
int_xmin = x0-6*sig_x; int_xmax = x0+6*sig_x;
Nk = 500; kmax = k0+4/sig_x; kmin = max(kmax/Nk, k0-4/sig_x);
tmin = 0; tmax = 20; Nt = 200;
ax = [xmin, xmax, -0.5, 0.5];
% ==============
psi = @(x) 1/(2*pi*sig_x^2)^0.25/sqrt(1 + 1i*t0/(2*m*sig_x^2))...
*exp(-(x-x0-k0*t0/m).^2/(2*sig_x)^2/(1 + 1i*t0/(2*m*sig_x^2)))...
.*exp(1i*k0*(x-k0*t0/(2*m)));
% 画 FFT 谱
x = linspace(int_xmin, int_xmax, Nx); dx = (int_xmax-int_xmin)/(Nx-1);
k = linspace(kmin, kmax, Nk); dk = (kmax-kmin)/(Nk-1);
[g1, k1] = FFT(fftresize(psi(x),Nx*2), int_xmin, dx);
g = interp1(k1, g1, k);
figure; plot(k, abs(g), '.');
axis([kmin, kmax, 0, max(abs(g))*1.1]);
xlabel k; title 'FFT of \psi(x)';
E = k.^2/(2*m);
x = linspace(xmin, xmax, Nx);
t = linspace(tmin, tmax, Nt);
figure; set(gcf, 'Unit', 'Normalized', 'Position', [0.1, 0.1, 0.4, 0.3]);
for it = 1:Nt
psi1 = zeros(size(x));
if mode == 1 % 使用类 sin,cos 基底
A = zeros(1,Nk); B = A;
for i = 1:Nk
A(i) = integral(@(x)psi(x).*FSB(x,E(i),L,m,V0,0), ...
int_xmin, int_xmax, 'RelTol', 1e-16);
B(i) = integral(@(x)psi(x).*FSB(x,E(i),L,m,V0,1), ...
int_xmin, int_xmax, 'RelTol', 1e-16);
psi1 = psi1 + dk*(A(i)*FSB(x,E(i),L,m,V0,0)...
+ B(i)*FSB(x,E(i),L,m,V0,1))*exp(-1i*E(i)*t(it));
end
else % 使用类 exp 基底
for i = 1:Nk
C1 = zeros(1,Nk); C2 = C1;
C1(i) = integral(@(x)psi(x).*conj(FSB2(x,E(i),L,m,V0,1)), ...
xmin, xmax, 'RelTol', 1e-16);
C2(i) = integral(@(x)psi(x).*conj(FSB2(x,E(i),L,m,V0,-1)), ...
xmin, xmax, 'RelTol', 1e-16);
psi1 = psi1 + dk*(C1(i)*FSB2(x,E(i),L,m,V0,1) ...
+ C2(i)*FSB2(x,E(i),L,m,V0,-1))*exp(-1i*E(i)*t(it));
end
end
clf;
plot(x, real(psi1)); hold on; plot(x, imag(psi1));
plot([xmin, -L, -L, L, L, xmax], [0, 0, 0.4, 0.4, 0, 0], 'k');
axis(ax); xlabel x; title(['t = ' num2str(t(it), '%.2f')]);
saveas(gcf, [num2str(it) '.png']);
end