Numerov算法解一维无限深势阱的问题 (含量子力学导论)

写在前面

• 原子物理期末论文要写一篇和量子隧穿的问题,然后,啊,写不出来.
• 请有志之人给点帮助!!!球球了!
• 之中可能会遇到这样的微分方程

• 所以需要学学Numerov算法来解决一些问题.

均匀格点上的Numerov算法

• 推导
• 事实上你可以自己推导出一种算法,然后坚称这是你自己发明的

• 整理：

 

一维无限深势阱

波函数

边界条件:

显然,在无限深势阱中的能量是分立的

为了说明这个问题,我们首先要回顾一下量子力学导论的部分

量子力学导论回顾

• 禁闭的波必然导出量子化的条件

习题选(杨福家 第五版)

程序设计

• 该问中,能量的解析解为

• 我们将不断调整的值,
• 代入上式,得:

• 为了能够找到合适的k值,我们应当设置一个tol值,当满足边界条件时,记录当前的k值

完整代码

import numpy as np
import matplotlib.pyplot as plt
 
def numerov(f,y0,dy0):
    h = 0.01
    y = np.zeros(len(f))
    y[0] = y0
    y[1] = y0 + h * dy0
    for i in range(1,len(f)-1):
        k1 = (1 + h**2/12. * f[i+1])
        k2 = 2 * (1 + h**2/12. * f[i])
        k3 = (1 + h**2/12. * f[i])
        y[i+1] = (k2*y[i] - k3*y[i-1] - h**2 * f[i] * y[i])/k1
    return y
 
tmin,tmax = 0,1
T = np.linspace(tmin,tmax,int((tmax-tmin)/0.01))
 
#归一化之前波函数启动点处的导数值可以随意取
 
y0  = 0
dy0 = 1
 
tol = 10 ** -6
k = 0
dk = 0.0001
u = 1
K,Kx = [],[]
 
for n in range(1000000):
    F = lambda r:k**2
    f = np.array([F(i) for i in T])
    Y = numerov(f,y0,dy0)
 
    if abs(Y[-1]) <= tol:
        K.append(Y[-1])
        Kx.append(k)
        
    if not n%5000:
        print(n/1000000)
    k += dk
    
 
import joblib
joblib.dump(K,"model-Numerov-infinity-K.pkl")
joblib.dump(Kx,"model-Numerov-infinity-Kx.pkl")
 
#plt.plot(Kx,K)
plt.scatter(Kx,K,s=2,c="red")
plt.plot(Kx,K,linewidth=2,c="blue")
plt.scatter([i*np.pi for i in range(1,30)],[0 for i in range(29)],s=2,c="green")
#plt.scatter([np.pi*i for i in range(10)],[0.03 for i in range(10)],s=1)
plt.show()

图像

 

•  它的误差在1e-6的量级上,至少说明还是和理论符合的差不多的