Crank-Nicolson 算法解一维含时薛定谔方程（Matlab）

贡献者： addis

预备知识 1　薛定谔方程（单粒子一维）

　　¹本文使用原子单位制。薛定谔方程为

\begin{matrix} (1) & - \frac{1}{2} \frac{\partial^{2} ψ}{\partial x^{2}} + V ψ = i \frac{\partial ψ}{\partial t} . \end{matrix}

传播子作用于波函数为

\begin{matrix} (2) & ψ (x, t + Δ t) = \exp (- i H Δ t) ψ (x, t) . \end{matrix}

用 Crank-Nicolson 或 Caley scheme² 得到的结果是

\begin{matrix} (3) & (1 + \frac{i}{2} H_{n + 1} Δ t) ψ_{n + 1} = (1 - \frac{i}{2} H_{n} Δ t) ψ_{n} . \end{matrix}

其中

ψ_{n}

是时刻

t_{n}

的波函数列矢量（已知），

ψ_{n + 1}

为时刻

t_{n + 1}

的波函数列矢量（未知），

H_{n}

是

t_{n}

时刻的哈密顿矩阵。

　　但事实上，还可以继续减少计算量。若近似认为³ $H_{n + 1} \approx H_{n}$ ，将式 3 整理后得

\begin{matrix} (4) & (\frac{1}{2} + \frac{i}{4} H_{n} Δ t) (ψ_{n + 1} + ψ_{n}) = ψ_{n} . \end{matrix}

解这个方程，再减去

ψ_{n}

即可。

1. 等间距网格

　　对于等间距坐标网格 $x_{1}, x_{2}, \dots$ ，可以用差分法（式 5 ）计算二阶导数，表示为矩阵有

\begin{matrix} (5) & D_{2} = \frac{1}{Δ x^{2}} (\begin{matrix} - 2 & 1 & 0 & 0 & \dots \\ 1 & - 2 & 1 & 0 & \dots \\ 0 & 1 & - 2 & 1 & \dots \\ 0 & 0 & 1 & - 2 & \dots \\ ⋮ & ⋮ & ⋮ & ⋮ & ⋱ \end{matrix}) . \end{matrix}

那么

H_{n} = - D_{2} / (2 m) + V_{n}

，其中

V_{n}

是对角矩阵，第

i

个对角元为

V (x_{i}, t_{n})

。现在，就可以代入式 3 或式 4 求解了。

2. 虚时间法

预备知识 2　虚时间法求基态波函数

　　虚时间法用于求解势阱中的基态，使用虚时间后，式 3 和式 4 分别变为

\begin{matrix} (6) & (1 + \frac{1}{2} H Δ t) ψ_{n + 1} = (1 - \frac{1}{2} H Δ t) ψ_{n}, \end{matrix}

\begin{matrix} (7) & (\frac{1}{2} + \frac{1}{4} H Δ t) (ψ_{n + 1} + ψ_{n}) = ψ_{n} . \end{matrix}

3. 代码

预备知识 3　Matlab 的稀疏矩阵

　　 Matlab 代码如下，使用式 4 ，以及 Matlab 的稀疏矩阵。势能函数可以在 V_fun 中设置，我们以方势垒为例，所有参数和 “高斯波包的方势垒散射数值计算（Matlab）” 相同。不同的是，由于我们使用迪利克雷边界条件，波函数到达边界后会发生全反射。要避免反射，可以用开放边界条件以及吸收势能。

图 1：运行结果

代码 1：TDSE_cn1d.m

% Crank-Nicolson 法解一维薛定谔方程
% 等间距网格，稀疏矩阵

function TDSE_cn1d
% ==== 参数设置 ======
xmin = -80; xmax = 80; Nx = 1000; % x 网格
tmin = 0; tmax = 20; Nt = 400; % 时间网格
Nplot = 10; % 画图步数
ax = [xmin, xmax, -0.5, 0.5];
% 高斯波包
x0 = -17; t0 = 0; m = 1; % 高斯波包的初始时间
p0 = 4; % 初始动量
sig_x = 2; % 初始位置， 波包宽度
% 方势垒 [-L,L], 高 V0
L = 4; V0 = 4;
V = @(x,t) V_fun(x,t,L,V0); % 势能函数
% ===================
close all;
psi_gs = @(x) 1/(2*pi*sig_x^2)^0.25/sqrt(1 + 1i*t0/(2*m*sig_x^2))...
      *exp(-(x-x0-p0*t0/m).^2/(2*sig_x)^2/(1 + 1i*t0/(2*m*sig_x^2)))...
      .*exp(1i*p0*(x-p0*t0/(2*m)));
x = linspace(xmin, xmax, Nx); dx = (xmax-xmin)/(Nx-1);
t = linspace(tmin, tmax, Nt); dt = (tmax-tmin)/(Nt-1);
psi = psi_gs(x).';
% 准备稀疏矩阵 [对角线, 上对角线, 下对角线]
ind1 = [1:Nx, 1:Nx-1, 2:Nx]; % 行标
ind2 = [1:Nx, 2:Nx, 1:Nx-1]; % 列标
% 动能矩阵非零元
T = -1/(2*m)*[-2*ones(1,Nx), ones(1,2*Nx-2)]/dx^2;
A = (1i*dt/4)*T; % 对角元稍后更新
figure; plot(x, real(psi)); hold on;
plot(x, imag(psi));
for it = 1:Nt
    % 更新对角元
    A(1:Nx) = 0.5 + (1i*dt/4)*(T(1:Nx) + V(x, t(it)+dt/2));
    Asp = sparse(ind1, ind2, A, Nx, Nx);
    tmp = Asp \ psi;
    % err = norm(Asp*tmp - psi);
    psi = tmp - psi;
    if mod(it, Nplot) == 0
        clf;
        plot(x, real(psi)); hold on;
        plot(x, imag(psi));
        plot([xmin, -L, -L, L, L, xmax], [0, 0, 0.4, 0.4, 0, 0], 'k');
        axis(ax);
        pause(0.2);
    end
end
end

% 势能函数：方势垒
function ret = V_fun(x,~,L,V0)
ret = zeros(size(x));
ret(x > -L & x < L) = V0;
end

1. ^ 参考 [1]。
2. ^ 二者是一回事，见 [1] 19.2 节。
3. ^ 或者更公平地，把式 4 中 $H_{n}$ 改为 $H_{n + 1 / 2}$ ，即 $(t_{n} + t_{n + 1}) / 2$ 时的哈密顿矩阵。

[1] ^ W. H. Press, et al. Numerical Recipes 3rd edition.

致读者：小时百科一直以来坚持所有内容免费无广告，这导致我们处于严重的亏损状态。长此以往很可能会最终导致我们不得不选择大量广告以及内容付费等。因此，我们请求广大读者热心打赏 ，使网站得以健康发展。如果看到这条信息的每位读者能慷慨打赏 20 元，我们一周就能脱离亏损，并在接下来的一年里向所有读者继续免费提供优质内容。但遗憾的是只有不到 1% 的读者愿意捐款，他们的付出帮助了 99% 的读者免费获取知识，我们在此表示感谢。