♥初体験♥数值求解一维谐振子本征问题

Posted on 2021-10-31 In 笔记 Views: Symbols count in article: 4.1k Reading time ≈ 4 mins.

介绍

Scientific Computation?

这是我第一次干和“科学计算”相关的事。之前虽然一直惦记着Scientific Computation，但实际上连Import Numpy都没有干过。于是在这个作业的驱动(?)下，我写了第一个有用的Python程序。虽然结果的代码量很少，但过程实属是历尽坎坷，有一些记录的价值和趣味在。

这也成为了这个博客的第一篇文章。

1 2	import numpy as np import matplotlib.pyplot as plt

一些背景

最近实在太忙，不然还想梳理一下到目前学习的内容。总的来说，老师想用一种不同的方法讲课，把一些理论力学、量子力学、统计力学的知识用“和化学相关”的方式串下来。我感觉这种课还是挺少见的。当然我也逐渐明白《化学中的数学》这个课名中的纠结与妥协——真的是很难说这是个什么课。

量子力学的部分是按着 Sakurai 的书来引入的，即直接给出“量子态存于希尔伯特空间”、“可观测量对应厄米算符”的公理化体系。接下来，我们给出动量算符的形式，解出了动量本征态（波函数）。之后进入了“一维无穷深势阱”，要做的是解哈密顿量的本征体系。一维无穷深势阱的势能是这个： \[ V(x)= 0 \text{ }(-L/2 \leq x \leq L/2) \] \[ V(x)=\infty,\text{otherwise} \]

在势阱的内部，哈密顿量只有动能项：\(\hat{H}=\hat{p}^2\big/2m\) .

利用相容算符有共同本征态的性质，将哈密顿量的本征态写成（一对正负的）动量本征态的线性组合，由边值条件解出了系数。然后就得到了一维无穷深势阱的本征波函数： \[ \braket{x|\psi_n}=\sqrt{\frac{2}{L} }\sin\left(\frac{n\pi(x+L/2)}{L}\right) \] 以及对应的能量本征值： \[ E=\frac{\hbar^2}{2m}\left(\frac{n\pi}{L}\right)^2 \] 量子数n在这里取正整数，因为对于sin函数而言，负的n只有相位上的区别，由波函数的相位自由度可以弥补。

至此我们拿到了这个本征问题的解。可以证明，这一组解构成了一个正交完备基组，可用来表出希尔伯特空间里所有的（可观测）态（其实就是傅里叶级数）。接下来，我们用这一个模型导出了配分函数等。在研究了“分子的平动”后，我们开始研究振动。于是就用到了谐振子模型。

任意一维体系的数值求解

既然我们知道无穷深势阱的解是一组完备基，那么可以用来表出任意我们想要的本征态。求解本征问题会变为求级数的系数，进一步地会变成一个矩阵对角化问题。为什么要用复杂的级数形式来求解呢？这是因为大多实际体系是无法给出解析解的。对于化学而言关注的总是数值解。

既然要做计算，我们就不能让n趋向无穷了，而是进行一个截断。得到N个“近似完备基组”，然后通过结果的收敛性来判断取值的合理（调参捏）。到了有限维，算符就便乘了我们熟悉的矩阵（厄米矩阵hso）。波函数自然也是要离散化的，也成为了一个很长的向量。

具体来讲，我们要求解的本征问题长成这样： \[ \hat{H}\ket{\psi_n}=\varepsilon_n\ket{\psi_n} \] 我们令解组成的矩阵\((\psi_1,\psi_2,\dots,\psi_n)={A}\) ，令上面提到的已知基组\((\phi_1,\phi_2,\dots,\phi_n)=B\) ，线性表出则可以用一个系数矩阵\(C\)写成： \[ A=CB \] 求解这个系数矩阵的方式就是：确定\(\hat{H}\)在基组\(B\)下的矩阵表示，然后对角化 \[ \hat{H}=C\lambda C^\dagger \] 然后就有 \[ \hat{H}(CB)=\lambda (CB) \] \(\lambda\)的对角元就是本征值，厄米矩阵的性质保证了它们都是实数。\(CB\)就包含了对应的各个本征矢。[写的时候有点仓促草]

所以要做的关键就是求哈密顿算符的矩阵元。等一下，\(\hat{H}=\hat{T}+\hat{V}\)，动能项我们已经知道了，是一个对角元是\(E_n\)的对角矩阵（和之前的并无区别）。需要求的是势能算符的矩阵元，方法是在位置空间展开单位算符： \[ V_{mn}=\bra{\phi_m}\hat{V}\ket{\phi_n}=\int_{-L/2}^{L/2}\phi_m^*(x)V(x)\phi_n(x) \mathrm{d}x \] 积分可以用梯形积分的方法进行数值计算。

我们要算的一维谐振子体系的势能为： \[ V(x)=\frac{1}{2}m\omega^2x^2 \] 为简便，令\(\hbar=1,m=1,\omega=2\), 而N，L，以及算积分时的划分点数量S是要调的参数。

结果与讨论

解决Python算得慢的问题

众所周知，Python是一种解释型语言，比C/C++、Fortran之类速度慢几个数量级。虽然不是用不了，但显然搞快点会更爽。办法就是”向量化“，然后交给Numpy库。Numpy背后是C代码，且带有了一些高级的线性代数算法，所以比Python自带的循环之类快得很。这是助教教给我们的策略，我在看懂以后就直接用了。而且这使得代码看起来很简洁。

ns=np.arange(1,N+1)
#build T matrix
t_ev=(1/2)*(np.pi/L)**2 * ns**2 #eigen values
Tmat=np.diag(t_ev)

这是动能的对角矩阵，构建方法就像是直接对n维向量作用一个函数。

1
2
3

#base vector v(x,n)
nx=np.outer(xs,ns)
bases = np.sin((np.pi/L) * nx)

基向量的构建方式是我很难想到的。事实上写成一个矩阵是我在写这篇文章时才梳理清楚的，虽然说本来就该这样……但线性代数的知识属实是忘了一些。在这里直接由x和n”张成“了一个矩阵，然后用我们的基组的表达式进行赋值，同样是避免了循环的出现。

#potential funcion
def fp(x):
    return 2*x**2
#build V matrix
xs=np.linspace(0,L,S)
dx=xs[1]-xs[0]
diag_fp = np.diag(fp(xs-L/2))
Vmat = np.dot(bases.T,np.dot(diag_fp,bases)) * 2/L *dx

定义一个势能。注意到为了计算的方便，令x的取值为0到L，这样势能函数也应当作相应的平移。我忘了这样的事因而卡了一段时间。

而当我们的基组写成了矩阵，求势能矩阵元就不用像定义种那般，用一个个基分别算积分。在这里就是三个矩阵相乘的形式（x已经离散化过了捏）。

那么事情基本上就做完了，然后就是对角化：

1
2
3

Hmat = Tmat + Vmat
H_e = np.linalg.eig(Hmat)
print(sorted(H_e[0]))

用到Numpy库里的ligalg.eig()，它得到的第一个array是特征值。但是……

解决我脑子的问题

run出来的结果不符合设想。本征值应当是 \[ \varepsilon_n = (n+\frac{1}{2})\hbar\omega \] 但是输出了一堆奇怪的数……

我坐在理教二层，调试了近三个小时，心态大崩。一直以为是势能矩阵写错了，后来又怀疑是参数问题，然而每调一回参数结果都会向意料外的方向发展。求助于同学才发现是，动能矩阵少写了n上的平方……

脑子的问题仍需要解决。

[1.0000000001548224, 3.000000000003902, 5.000000000046899, 7.000000000002921, 9.000000000025032, 11.000000000003851, 13.000000000021144, 15.00000000000944, 17.000000000017277, 19.00000000005756, 21.00000000053884, 23.00000000387303, 25.000000029432897, 27.000000170822677, 29.000001044443803, 31.000004965387035, 33.00002472004807, 35.000097260956586, 37.000396429511035, 39.0012947512046] #20 eigenvalues

本征态画图

首先，我们要对本征态（向量）按本征值进行排序。这部分由 Nipporita 提供了帮助，原因是我对于 Python 的掌握实属不彳亍。

H_eva, H_eve = np.linalg.eig(Hmat)
H_e = [(H_eva[i], H_eve.T[i]) for i in range(N)]
H_e_sort = sorted(H_e, key = lambda x: x[0])
print([i[0] for i in H_e_sort[:20]]) #new method

但是我接下来画出的是一坨奇怪的东西。和想象中的完全不同，甚至处处不连续。

在修改时破坏了 Nipporita 留下的代码，不会改回了，于是又求助。顺便发现了我的问题：我试图对系数矩阵画图而没有乘上基组（什么脑子）。本质上是在没有想明白原理的时候就开干了，但这没有带来任何的方便。

最终有了画图的代码:

1
2
3

xs_=np.linspace(-L/2,L/2,S)
for i in range(20): #plot
    plt.scatter(xs_,(np.dot(bases, H_e_sort[i][1])), s = 0.5)

以及图：

以下是换了一种方法，把波函数在纵向排开：

完工

好！

结语

没什么可说的，这个文章没有价值。该干的事情还有很多。写这个只是为了测试一下 Hexo-NexT 的 MathJax 渲染以及代码高亮效果。

我被更多的作业狂暴鸿儒。等有空了再来维护博客。