一维薛定谔方程不稳定的差分法数值解（Matlab）

贡献者： addis

预备知识　薛定谔方程（单粒子一维），一维波动方程的简单数值解（Matlab）

　　本文使用原子单位制。本文演示如何用差分法解一维薛定谔方程。

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

　　乍看之下，类比 “一维波动方程的简单数值解（Matlab）”，我们似乎也能用简单的差分法求解式 1 ，但事实证明这是行不通的，迭代几步以后，数值误差就会盖过波函数本身。这说明对薛定谔方程来说，差分法是一个不稳定的算法。但为了教学，还是把公式和代码给出来。

　　一些稳定的算法如直接对角化（未完成）、Crank-Nicolson、以及 Lanczos 算法（通常配合 FEDVR 网格一起使用）。

　　把波函数取离散值，令 $ψ_{i, n} = ψ (x_{i}, t_{n})$ 。用有限差分表示二阶导数（式 5 ），得

\begin{matrix} (2) & - \frac{1}{2 m} \frac{ψ_{i - 1, n} - 2 ψ_{i, n} + ψ_{i + 1, n}}{Δ x^{2}} + V_{i, n} ψ_{i, n} = i \frac{ψ_{i, n + 1} - ψ_{i, n}}{Δ t} . \end{matrix}

于是

\begin{matrix} (3) & ψ_{i, n + 1} = \frac{i Δ t}{2 m Δ x^{2}} (ψ_{i - 1, n} - 2 ψ_{i, n} + Δ t ψ_{i + 1, n}) + (1 - i Δ t V_{i, n}) ψ_{i, n} . \end{matrix}

　　失败的差分法程序如下

图 1：程序运行的第 24 个循环，噪音已经使波函数扭曲。即使把时间步长和空间步长设置得很小也不会有多大改善。

代码 1：TDSE_1d_failed.m

% 差分法解一维薛定谔方程

function TDSE_1d_failed
% ==== 参数设置 ======
m = 1; % 质量，角频率
xmin = -10; xmax = 30; Nx = 1000; % x 网格
tmin = 0; tmax = 10; Nt = 3000; % 时间网格
Nplot = 1; % 画图步数
t0 = 0; % 高斯波包的初始时间
p0 = 5; % 初始动量
x0 = 0; sig_x = 2; % 初始位置， 波包宽度
V = @(x,t) zeros(size(x));
% ===================
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);
figure; plot(x, real(psi)); hold on;
plot(x, imag(psi));
for it = 1:Nt
    d2psi = [0, psi(1:end-2) - 2*psi(2:end-1) + psi(3:end), 0];
    psi = 1i*dt/(2*m*dx^2)*d2psi + (1 - 1i*dt*V(x,t(it))).*psi;
    psi(1) = 0; psi(end) = 0;
    if mod(it, Nplot) == 0
        clf;
        plot(x, real(psi)); hold on;
        plot(x, imag(psi)); axis([xmin,xmax,-0.5,0.5]);
    end
end
end

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