计算物理(二)Schrodinger方程
fengxiaot Lv4

本篇文章是西安交大《计算物理与程序设计》课程的笔记,程序语言为 MATLAB。内容主要包括用打靶法求解定态Schrodinger方程的能级,以及用Thomas算法求解化为三对角线性系统的含时薛定谔方程。

定态Schrodinger方程

求解方法

一维定态Schrodinger方程为

Hφ=EφH \varphi = E \varphi

22md2φdx2+Vφ=Eφ-\frac{\hbar^2}{2m} \frac{\mathrm{d}^{2}\varphi}{\mathrm{d}x^2} + V\varphi=E\varphi

这实质上是二阶常微分方程的本征值问题,其中 EE 是待求本征值,或者用物理的术语「能级」。束缚态要求波函数在边界处趋于 00 ,不妨假设 ±a\pm a 处有一个无限深势阱壁,阱外已经超出了我们研究的区域

φ(±a)=0\varphi(\pm a)=0

这样一来便可以使用打靶法配合 Numerov 算法求解此边值问题。此处初值为 φ(a)=0\varphi(- a)=0φ(a)=K\varphi^\prime (-a) = KKK 为可以任意选取的常数,靶心为 φ(a)=0\varphi(a)=0,变化的子弹为 EE

/* 注:子弹不是 φ(a)\varphi^\prime (-a) !只要找到本征能级 EE,即便波函数导数初值 φ(a)\varphi^\prime (-a) 任取,Schrodinger 方程也能使之很快收敛到正确的本征波函数上。 */

斜势阱壁

处理斜势阱壁需要使用双向积分方法,以保证经典禁戒区的波函数不包含指数增长成分。

如图所示,选取能级与势阱的交点为接合点 xcx_c ,此时靶心为

1R[φ<(xc)φ>(xc+)]=0\frac{1}{R}\left[ \varphi_<^\prime (x_c^-) -\varphi_>^\prime (x_c^+) \right] = 0

其中 RR 为合适的缩放因子,使得接合点 xcx_c 处波函数左右导数之差与二分法误差限 δ\delta 相当,从而跳过间断点。RR 一般取 11 即可。


含时Schrodinger方程

差分形式

一维含时 Schrodinger 方程为

iψt=22m2ψx2+Vψ\mathrm{i} \hbar\frac{\partial \psi}{\partial t} = -\frac{\hbar^2}{2m} \frac{\partial^{2}{\psi}}{\partial{x}^{2}} + V\psi

按照 Crank-Nicolson 算法,等式右边的每个 ψ\psi 应当替换为当前时刻和下一时刻波函数的算术平均 12(ψj+ψj+1)\frac{1}{2}\left(\psi^j + \psi^{j+1}\right) 再作差分

iψij+1ψijΔt=22m12[ψi+1j2ψij+ψi1j(Δx)2+ψi+1j+12ψij+1+ψi1j+1(Δx)2]+12Vi(ψij+ψij+1)\mathrm{i} \hbar\frac{\psi^{j+1}_i-\psi^{j}_i}{\Delta t} = -\frac{\hbar^2}{2m} \frac{1}{2} \left[\frac{\psi_{i+1}^j-2\psi_i^j+\psi_{i-1}^j}{(\Delta x)^2} + \frac{\psi_{i+1}^{j+1}-2\psi_i^{j+1}+\psi_{i-1}^{j+1}}{(\Delta x)^2} \right] + \frac{1}{2} V_i\left(\psi^j_i + \psi^{j+1}_i\right)

定义 r=iΔt2m(Δx)2r=\frac{\mathrm{i} \hbar \Delta t }{2m(\Delta x)^2} ,则有差分方程

r2ψi1j+1+(1+r+iΔt2Vi)ψij+1r2ψi+1j+1=r2ψi+1j+(1riΔt2Vi)ψij+r2ψi1j-\frac{r}{2} \psi_{i-1}^{j+1}+\left(1+r+\frac{\mathrm{i} \Delta t}{2 \hbar}V_i\right) \psi_{i}^{j+1} - \frac{r}{2}\psi_{i+1}^{j+1} = \frac{r}{2} \psi_{i+1}^{j}+\left(1-r-\frac{\mathrm{i} \Delta t}{2 \hbar}V_i\right)\psi_{i}^{j} + \frac{r}{2}\psi_{i-1}^{j}

仍然为三对角方程,使用附录中Thomas算法即可求解。

算符形式

差分方程亦可通过算符的形式写出。定义哈密顿量算符

H=22mδ2+VH = -\frac{\hbar^2}{2m} \delta^2 +V

系统的演化可以写为

ψ(t+Δt)=eiHΔt/ψ(t)=eiHΔt/2eiHΔt/2ψ(t)    eiHΔt/2ψ(t+Δt)=eiHΔt/2ψ(t)    (1+iΔt2H)ψ(t+Δt)=(1iΔt2H)ψ(t)\begin{aligned} &|\psi(t+\Delta t)\rangle = \mathrm{e}^{-\mathrm{i}H\Delta t / \hbar}|\psi(t)\rangle=\frac{\mathrm{e}^{-\mathrm{i}H\Delta t/2\hbar}}{\mathrm{e}^{\mathrm{i}H\Delta t/2\hbar}} |\psi(t)\rangle \\ \implies \, & \mathrm{e}^{\mathrm{i}H\Delta t/2\hbar} |\psi(t+\Delta t)\rangle = \mathrm{e}^{-\mathrm{i}H\Delta t/2\hbar} |\psi(t)\rangle \\ \implies \, & \left( 1+ \frac{\mathrm{i}\Delta t}{2\hbar}H\right) \psi(t+\Delta t) = \left( 1- \frac{\mathrm{i}\Delta t}{2\hbar}H\right) \psi(t) \end{aligned}

可以验证,结果与上述差分方程完全相同。实际计算中仍然需要使用差分方程。但算符形式好处在于可以推广到哈密顿量 HH 也含时的情形,可以参考这篇来自 Computer Physics Communications 的文献[1]

MATLAB代码

实际计算中,可以进行一个小变换,令 χij=ψij+1+ψij\chi_i^j = \psi_{i}^{j+1} + \psi_{i}^{j},于是差分方程化为

r2χi1j+(1+r+iΔt2Vi)χijr2χi+1j=2ψij-\frac{r}{2} \chi_{i-1}^j+\left(1+r+\frac{\mathrm{i} \Delta t}{2 \hbar}V_i\right) \chi_{i}^j - \frac{r}{2}\chi_{i+1}^j = 2\psi_{i}^{j}

然后便可求解此三对角方程,注意在等式右边首项和末项加上 r2χ1j\frac{r}{2} \chi_1^jr2χNj\frac{r}{2} \chi^j_N 。每一步求出 χij\chi_i^j 后减去 ψij\psi_{i}^{j} 即可得到 ψij+1\psi_i^{j+1},再进行下一步求解。此变换意义不明,但许多教材上都保留了下来。


附录

三对角线性系统

在线性代数中,三对角矩阵是矩阵的一种。一个三对角矩阵的非零系数在如下的三条对角线上:主对角线、低对角线、高对角线。在许多物理问题中,三对角矩阵常常作为原始数据出现,因此它们本身是很重要的,这种矩阵仅有 2n12n-1 个独立的元素。

三对角矩阵的生成

设主对角线上的元素为 nn 维向量 b\boldsymbol{b},低对角线和高对角线上的元素分别为 n1n-1 维向量 a\boldsymbol{a}c\boldsymbol{c} 。则 MATLAB 中可用 diag 命令生成

1
A = diag(b,0) + diag(a,-1) + diag(c,1);

特殊地,若对角线上各个元素均相同,则可以用 ones 生成

1
A = diag(b*ones(n,1),0) + diag(a*ones(n-1,1),-1) + diag(c*ones(n-1,1),1);

三对角线性方程组的求解

三对角线性方程组指形如

aixi1+bixi+cixi+1=dia_{i}x_{i-1}+b_{i}x_{i}+c_{i}x_{i+1}=d_{i}

其中 a1=0a_1=0cn=0c_n=0 。利用三对角矩阵 AA ,可以记为

Ax=dA \boldsymbol{x} = \boldsymbol{d}

[b1c1a1b2c2a2cn2an2bn1cn1an1bn][x1x2x3xn1xn]=[d1d2d3dn1dn]\begin{bmatrix} b_1 & c_1 & & & & \\ a_1 & b_2 & c_2 & & & \\ & a_2 & \ddots & \ddots & & \\ & & \ddots & \ddots & c_{n-2} & \\ & & & a_{n-2} & b_{n-1} & c_{n-1} \\ & & & & a_{n-1} & b_n \end{bmatrix} \begin{bmatrix} x_1 \\ x_2 \\ x_3 \\ \vdots \\ x_{n-1} \\ x_n \end{bmatrix}= \begin{bmatrix} d_1 \\ d_2 \\ d_3 \\ \vdots \\ d_{n-1} \\ d_n \end{bmatrix}

对于三对角线性方程组,不必使用时间复杂度为 O(n3)O(n^3) 的高斯消元法,可以使用Thomas算法达到仅 O(n)O(n) 的时间复杂度。

Thomas算法的MATLAB代码为

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
function x = tridiag(A,d)
% Solves the tridiagonal linear system Ax = d for x
% using the matrix implementation of the tridiagonal matrix algorithm
% tridiag.m

% determines n
n = length(d);

% preallocates x
x = zeros(n,1);

% forward elimination
for i = 2:n
w = A(i,i-1)/A(i-1,i-1);
A(i,i) = A(i,i)-w*A(i-1,i);
d(i) = d(i)-w*d(i-1);
end

% backward substitution
x(n) = d(n)/A(n,n);
for i = (n-1):(-1):1
x(i) = (d(i)-A(i,i+1)*x(i+1))/A(i,i);
end

end

三对角矩阵的特征值

待更

数值积分

求解Schrodinger方程过程中常常需要对波函数进行归一化,而波函数往往是离散点列而非连续函数,不方便调用 integral 进行积分。

这里提供一个 numint(X,Y) 函数,根据 XX 指定的坐标与标量间距对 YY 进行积分。当 XX 包含偶数个区间(奇数个散点)时,使用Simpson算法;否则使用MATLAB自带的梯形算法 trapz(X,Y)

1
2
3
4
5
6
7
8
9
10
11
12
function integration = numint(X,Y)
% numint computes the integral of Y with respect to X
% Use the Simpson rule or Trapz automatically based on dim of X
% numint.m
N = length(X);
if mod(N,2) == 0 % if N is even, use trapz
integration = trapz(X,Y);
else % if N is odd, use simpson
dx = X(2) - X(1);
integration = dx/3*(Y(1)+Y(N)+4*sum(Y(2:2:N))+2*sum(Y(3:2:N-1)));
end
end

间断函数的零点

面对间断函数时,即便是预搜索二分法也难以找到函数零点。

如图所示,实心点为周期间断函数真正的零点,而空心点为间断点。预搜索二分法固然能判断并找到实心点作为零点,但在空心点两侧函数异号,预搜索认为此处也存在一个零点。无论二分区间 (a,b)(a,b) 的长度 ab|a-b| 怎样逼近给定的 ε\varepsilon,区间两侧总是异号,陷入死循环。

处理方法是,引入二分法偏差返回值 deviation 和误差限 δ\delta 。偏差返回值 deviation 定义为 ab|a-b| 缩小到 ε\varepsilon 时,函数值在区间端点的差值 f(a)f(b)|f(a)-f(b)| 。对于真正的零点,deviation 应当甚小,能够小于给定的误差限 δ\delta 。对于间断点,无论 (a,b)(a,b) 如何缩小都始终分布在间断点异侧,f(a)f(b)|f(a)-f(b)| 甚大,因为函数在此有阶跃。

因而可将二分法搜寻完成后的 deviation 是否小于 δ\delta 作为判据进行优化。若 deviation 大于 δ\delta ,则说明遭遇了间断点,从 root+step 开始重新搜索。以下为优化后的预搜索二分查找零点函数 prebisection(f,x0,step,epsilon,delta)

1
2
3
4
5
6
7
8
9
10
11
12
function root = prebisection(f,x0,step,epsilon,delta)
% prebisection.m
while f(x0+step)*f(x0)>0
x0 = x0 + step;
end
[root,deviation] = bisection(f,x0,x0+step,epsilon);
% If deviation > delta even though |a-b| < epislon, it indicates a leap here
% When finding a leap, skip it, start new searching from root+step
if deviation > delta
root = prebisection(f,root+step,step,epsilon,delta);
end
end

二分查找零点参见计算物理第一篇文章中的 bisection.m



  1. H. Gharibnejad, B.I. Schneider, M. Leadingham, H.J. Schmale, A comparison of numerical approaches to the solution of the time-dependent Schrödinger equation in one dimension, Computer Physics Communications, Volume 252, 2020, 106808, ISSN 0010-4655, https://doi.org/10.1016/j.cpc.2019.05.019. ↩︎