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

贡献者： addis

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

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

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

　　 把波函数取离散值，令 $\psi_{i,n} = \psi(x_i,t_n)$。用有限差分表示二阶导数（式 5 ），得

　　 失败的差分法程序如下

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

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