用 Numerov 方法求解一维薛定谔方程 (python)
solving 1D Schrödinger equation with Numerov method (python)
晚上好。
我目前正在尝试求解一维薛定谔方程。 (与时间无关)与 Numerov 方法。该方法的推导对我来说很清楚,但我在实现时遇到了一些问题。我试图在 google 上寻找解决方案,并且有一些 (like this one or this one),但我真的不明白他们在他们的代码中做了什么......
问题:
通过一些数学运算,您可以得到以下形式的等式:
其中 。首先,我想看看潜力 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()
由于情节看起来根本不像波函数,所以肯定有问题,但我很难找出它是什么。如果有人能提供一点帮助,那就太好了。
谢谢司徒
不幸的是我不太记得量子物理所以我不明白一些细节。我仍然在您的代码中看到一些错误:
为什么在 numerov_step
里面你平方 k1
、k2
和 k3
?
在你的主周期
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)
- 这是我不太明白的部分。在你的 first link 上查看动画,看起来这个等式只对
V(x)
和 E
的某些组合有很好的解决方案,在其他情况下它很快就会变得疯狂。看起来你的 V(x)
和 E
与 hbar
和 V(x)
的比例与引用的文章完全不同,这可能是解决方案变得疯狂的另一个原因。
晚上好。
我目前正在尝试求解一维薛定谔方程。 (与时间无关)与 Numerov 方法。该方法的推导对我来说很清楚,但我在实现时遇到了一些问题。我试图在 google 上寻找解决方案,并且有一些 (like this one or this one),但我真的不明白他们在他们的代码中做了什么......
问题:
通过一些数学运算,您可以得到以下形式的等式:
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()
由于情节看起来根本不像波函数,所以肯定有问题,但我很难找出它是什么。如果有人能提供一点帮助,那就太好了。
谢谢司徒
不幸的是我不太记得量子物理所以我不明白一些细节。我仍然在您的代码中看到一些错误:
为什么在
numerov_step
里面你平方k1
、k2
和k3
?在你的主周期
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)
- 这是我不太明白的部分。在你的 first link 上查看动画,看起来这个等式只对
V(x)
和E
的某些组合有很好的解决方案,在其他情况下它很快就会变得疯狂。看起来你的V(x)
和E
与hbar
和V(x)
的比例与引用的文章完全不同,这可能是解决方案变得疯狂的另一个原因。