计算物理专题：有限差分法解决本征值问题

River Chandler

已于 2023-06-24 00:29:46 修改

阅读量666

点赞数

分类专栏： # 科学计算原理文章标签：算法 python 开发语言

于 2023-06-23 22:15:53 首次发布

本文链接：https://blog.csdn.net/Chandler_river/article/details/131354993

版权

科学计算原理专栏收录该内容

24 篇文章 2 订阅

订阅专栏

计算物理专题：有限差分法解决本征值问题

定态薛定谔方程差分形式

一维定态薛定谔方程

$-\alpha\frac{d^2}{dx^2}\psi(x)+V(x)\psi(x)=E\psi(x);\alpha=\frac{\hbar^2}{2\mu}$

$\left\{\begin{matrix} \psi(kh)=\psi_k\\ R=\alpha/h^2\\ a_k=2R+V(kh) \end{matrix}\right.$

$\begin{bmatrix} a_1& -R & & & \\ -R&a_2 &-R & & \\ & -R& a_3 &-R & \\ & & & \ddots& \\ & & & -R & a_n \end{bmatrix}\begin{bmatrix} \psi_1\\ \psi_2\\ \psi_3\\ \cdots\\ \psi_n \end{bmatrix}=E\begin{bmatrix} \psi_1\\ \psi_2\\ \psi_3\\ \cdots\\ \psi_n \end{bmatrix}$

谐振子

$V(x)=\frac{1}{2}\hbar\omega^2x^2$

$\hbar=1;\omega=1;\mu=1$

$x_0=0;x_1=20;h=0.01;$

解法代码

import numpy as np
def householder(symmetric_matrix):
    M = symmetric_matrix
    assert np.allclose(M,M.T),"matrix is not symmetric"
    N = len(M)
    for i in range(1,N-1):
        r = 0
        for j in range(i,N):
            r += M[i-1][j]**2
        r = r**0.5
        if r * M[i-1][i] > 0:
            r *= -1
        Ki = -1/(r**2 - r*M[i-1][i])
        Pi = np.zeros((N,1))
        Pi[i][0] = M[i-1][i] - r
        for j in range(i+2,N+1):
            Pi[j-1][0] = M[i-1][j-1]
        Fi = np.dot(Pi,Pi.T)*Ki
        for j in range(1,N+1):
            Fi[j-1][j-1] += 1        
        M = np.dot(np.dot(Fi,M),Fi)
    return M

def qr_decomposition(matrix):
    m,n = matrix.shape
    Q = np.zeros((m,n))
    R = np.zeros((n,n))
    for j in range(n):
        v = matrix[:,j]
        for i in range(j):
            R[i,j] = np.dot(Q[:,i],matrix[:,j])
            v = v - R[i,j] * Q[:,i]
        R[j,j] = np.linalg.norm(v)
        Q[:,j] = v/R[j,j]
    return Q,R

def QLmethod(matrix,tol=1e-6,maxiter=100):
    assert np.allclose(matrix,matrix.T),"matrix is not symmetric"
    matrix = householder(matrix)
    n = matrix.shape[0]
    eigenvalues = np.zeros(n)
    iterations = 0
    while np.max(np.abs(matrix.diagonal(offset=1))) > tol and iterations < maxiter:
        q,r = qr_decomposition(matrix)
        matrix = np.dot(r,q)
        iterations += 1
    eigenvalues = matrix.diagonal()
    return eigenvalues


hbar = 1
omega = 1
mu = 1
x0 = 0
x1 = 10
h = 0.1

V = lambda x:1/2*hbar*omega**2*x**2
R = hbar**2/2/mu/h**2
N = int((x1-x0)/h)
List = [2*R + V(k*h) for k in range(1,1+N)]

symmetric_matrix = np.diag(List)
for i in range(N-1):
    symmetric_matrix[i][i+1] = -R
    symmetric_matrix[i+1][i] = -R

#Three Symmetric matrix
eigenvalue = QLmethod(symmetric_matrix)

运行结果

array([199.05965221, 197.85154038, 196.69647447, 195.59251125,
       194.53779897, 193.53057225, 192.56914726, 191.65191722,
       190.7773482 , 189.94397515, 189.15039823, 188.39527932,
       187.6773388 , 186.99535243, 186.34814856, 185.73460537,
       185.15364836, 184.60424796, 184.08541726, 183.59620995,
       183.13571823, 182.70307101, 182.29743212, 181.91799858,
       181.56399908, 181.23469247, 180.92936632, 180.64733559,
       180.38794141, 180.15054982, 179.93455068, 179.73935657,
       179.56440178, 179.40914085, 179.27304046, 179.15548747,
       179.0549045 , 178.96231954, 178.82813537, 178.45299393,
       177.32375035, 174.73750414, 170.45930442, 165.03198466,
       159.18519843, 153.34588073, 147.67451019, 142.21724325,
       136.98169564, 131.96220481, 127.1481637 , 122.52724056,
       118.08680316, 113.81459325, 109.69905122, 105.72946646,
       101.89603735,  98.18988381,  94.60303514,  91.12840579,
        87.7597663 ,  84.49171355,  81.31964293,  78.2397245 ,
        75.24888443,  72.34479386,  69.52586706,  66.791272  ,
        64.14095638,  61.57569128,  59.0971279 ,  56.70783771,
        54.41122812,  52.21100812,  50.10939318,  48.1026652 ,
        46.17390347,  44.28931575,  42.41015044,  40.51421357,
        38.60050745,  36.67629929,  34.74758658,  32.81748476,
        30.88705093,  28.95607702,  27.02366776,  25.08868837,
        23.15008972,  21.20709216,  19.25922869,  17.30628391,
        15.34819203,  13.38495142,  11.41657896,   9.44309464,
         7.46451849,   5.4808703 ,   3.49216962,   1.49843574])

存疑

谐振子势确实是等距的，但是为什么计算出的结果与E_n=(n+\frac{1}{2})\hbar\omega 相差了一个1，这是为什么呢？没想明白。

向Fortran的翻译

有限差分法解决本征值问题的特点

引用有限差分法解决本征值问题时，步长的大小和边值条件会极大得影响计算的精度，有限差分法一般可以提供很好的近似解

量子隧穿效应的调研

\subsubsection{量子隧穿与核聚变}
   恒星原子核之间存在很高的排斥静电力。对于核聚变，需要大量的能量使原子核彼此足够接近以克服库仑力：
   原子核的平均动能可以被表示为
   \begin{equation}
       E_{kin}=\frac{3}{2}kT
   \end{equation}
   需要被克服的能量$E_{cb}$为：
   \begin{equation}
       E_{cd}=\frac{Z_1Z_2e^2}{r_c}
   \end{equation}
   由于恒星内部的原子核在静止燃烧过程中无法克服库伦势垒，因此除了热运动以外，必须要存在其他因素才能使得热核反应成为可能。尽管恒星内部的热运动低到无法直接引起热核反应，但它能在库仑势垒后面产生显著的概率密度使得量子隧穿能以一定的概率发生。这一事件所需要的能量远远低于克服库仑势垒所需的动能。对于运动速度为$\nu$的核子发生量子隧穿的概率$P_t(\nu)$正比于：
   \begin{equation}
       \exp{\frac{-4\pi ^2Z_1 Z_2 e^2}{h}\frac{1}{\nu}}
   \end{equation}

\subsubsection{分子生物学中量子隧穿}
   比氢原子重的化学元素合成与量子隧穿现象密切相关。对氘代细菌孢子的实验表明，与水基培养基中的自发突变率相比，在含氘氧化物的培养基中生长的细胞中的自发点突变率较低。这意味着突变偏差对所涉及的粒子质量的依赖性。而量子隧穿速率强烈依赖于粒子质量。$Löwdin DNA $突变模型描述了互补碱基对之间氢键所涉及的质子的量子性质。$Löwdin$的模型基于稳定的沃森-克里克碱基配对和稳定的遗传信息需要质子不改变其在碱基对中的位置这一前提，将质子在碱基对中的位置用所涉及的核酸碱基分子的规范形式来确定。$Löwdin$模型指出，通过质子对转移势垒的量子隧穿，当质子在DNA复制后立即处于量子力学非稳态时，质子在核酸碱基内的位置可能会以很小的概率因量子隧穿发生变化。\par
   Lödwin模型描述了DNA中质子的量子隧穿导致的普遍点突变偏差。更多的研究表明量子隧穿不仅能起到对DNA的破坏，而且还涉及到DNA修复问题。紫外线照射可以通过诱导DNA链中凸起的形成来损伤DNA。这些凸起是由于DNA链上相邻的嘧啶在紫外线照射下发生二聚而引起的。黄蛋白光解酶能够通过电子转移分裂共价连接的嘧啶来修复这种变形的DNA，从而使它们恢复到正常的单体形式，不再形成二聚体。长距离的超交换介导隧穿使得蛋白质中电子的量子隧穿可以在高达3纳米的距离上发生，这种势垒宽度远大于真空中的可能宽度。蛋白质通过超能变化介导的机制提供的低电子态提高了电子转移速率，即使在如此长的距离上它们仍然相对较高。蛋白质对损伤DNA的修复基本上基于长程电子隧穿。\par

\subsubsection{量子隧穿与语音识别$^{[3]}$}
   语音信号通常服从两种分布：较短语音信号服从Gauss分布，较长语音信号服从Laplace分布。对于同一说话人，由于自身身体条件，对于同样的语言内容，发音频率总在误差允许的范围内相似。对于不同的说话人，一方面发出的声音信号存在差异，这种差异表现为频率特征上的差异，可由处于不同稳定的量子态来描述。\par
   另一方面，对不同的语音信号分帧处理后，由于没帧的时间很短，基本上服从Gauss分布。因此，每一分帧语音信号可以被视作一个包含一组频率特征的量子波函数。不同的频率，通过相同的势垒，其隧穿系数不同，隧穿后的频率不同。\par
   因此，通过设置一组势垒，让每一个势垒有唯一的一个频率透射即可构成一组特征向量。这些特征都是非负的，再构成一个随机向量，可以用高斯向量表征。根据向量中元素拟合，降维后成的二维概率密度函数，通过最大似然估计就能实现对说话人的识别。\par