用 Numerov 方法求解一维薛定谔方程 (python)

问题描述 投票:0回答:1

我目前正在尝试求解一维薛定谔方程。 (与时间无关)使用 Numerov 方法。该方法的推导对我来说很清楚,但我在实现时遇到了一些问题。我尝试在谷歌上寻找解决方案,有一些(像这个这个),但我不太明白他们在代码中做了什么......

问题:

通过一些数学计算,您可以得到以下形式的方程: 哪里。首先,我想看看潜力

V(x)=1 if -a<x<a

由于我没有能量值或 Psi 的第一个值(启动算法所需的),我只是猜测一些......

代码如下所示:

import numpy as np
import matplotlib.pyplot as plt
from scipy.constants import hbar

m= 1e-27
E= 0.5

def numerov_step(psi_1,psi_2,k1,k2,k3,h):
    #k1=k_(n-1), k2=k_n, k3=k_(n+1)
    #psi_1 = psi_(n-1) and psi_2=psi_n
    m = 2*(1-5/12. * h**2 * k2**2)*psi_2
    n = (1+1/12.*h**2*k1**2)*psi_1
    o = 1 + 1/12. *h**2 *k3**2
    return (m-n)/o

def numerov(N,x0,xE,a):
    x,dx = np.linspace(x0,xE,N+1,retstep=True)

    def V(x,a):
        if (np.abs(x)<a):
            return 1
        else:
            return 0
    
    k = np.zeros(N+1)
    for i in range(len(k)):
        k[i] = 2*m*(E-V(x[i],a))/hbar**2

    psi= np.zeros(N+1)
    psi[0]=0
    psi[1]=0.1    

    for j in np.arange(2,N):
        psi[j+1]= numerov_step(psi[j],psi[j+1],k[j-1],k[j],k[j+1],dx)

    return psi

x0 =-10
xE = 10
N =1000

psi=numerov(N,x0,xE,3)

x = np.linspace(x0,xE,N+1)

plt.figure()
plt.plot(x,psi)
plt.show()

由于情节看起来根本不像波函数,肯定有什么问题,但我很难找出它是什么..如果有人能帮忙一点就好了。

python math physics numerical-methods differential-equations
1个回答
2
投票

不幸的是我不太记得量子物理学所以我不明白一些细节。我仍然在您的代码中看到一些错误:

  1. 为什么在

    numerov_step
    里面你方形
    k1
    k2
    k3

  2. 在你的主周期中

   for j in np.arange(2,N):
      psi[j+1]= numerov_step(psi[j],psi[j+1],k[j-1],k[j],k[j+1],dx)

你搞乱了索引。看起来这条线应该是

   for j in np.arange(2, N):
     psi[j] = numerov_step(psi[j - 2], psi[j - 1], k[j - 2], k[j - 1], k[j], dx)
  1. 这是我不太明白的部分。在您的第一个链接上查看动画,看起来这个方程仅对于
    V(x)
    E
    的某些组合有很好的解决方案,在其他情况下它很快就会变得疯狂。看来您的
    V(x)
    以及
    E
    hbar
    V(x)
    的比例与参考文章有很大不同,这可能是解决方案变得疯狂的又一个原因。
© www.soinside.com 2019 - 2024. All rights reserved.