无限深势井的电子波包模拟
假设我们有一个电子,处于 (0, L) 的一维无限深势井中,初始状态为中心位于 x0 = L2 处的波包:
其中 L = 10^(−8)m, σ = 10^(−10)m, k = 5×10^10m^(−1)。
另外物理常数 ̄h = 1.0546 × 10^(−34)J · s,Me = 9.109 × 10^(−31)kg。 设定时间步长 ht = 10^(−18)s,坐标分割数 N = 1000,坐标步长 hx = L/N。
以下模拟是在时间 t = 0, 300ht, 500ht, 800ht, 2000ht 时刻,波函数以及概率密度的分布和变化
代码:
import numpy as np
import matplotlib.pyplot as plt
#设定参数
L = 10.0**(-8)
sgm = 10.0**(-10)
k = 5.0*10**10
hbar = 1.0546*10**(-34)
Me = 9.109*10**(-31)
x0 = L/2
N = 1000
ht = 10.0**(-18)
hx = L/N
a1 = 1 + (1j*hbar*ht)/(2*Me*hx**2)
a2 = -(1j*hbar*ht)/(4*Me*hx**2)
b1 = 1 - (1j*hbar*ht)/(2*Me*hx**2)
b2 = (1j*hbar*ht)/(4*Me*hx**2)
#定义初始波函数
x = np.arange(0,L+hx,hx)
t = np.arange(0,2000*ht+ht,ht)
psi0 = np.exp(-(x-x0)**2/(2*sgm**2))*np.exp(1j*k*x)
#定义中间格点波函数所构成的矩阵
phi = np.zeros((len(psi0), len(t)),dtype=complex)
phi[:,0]=psi0
#定义传递矩阵
A = np.zeros((len(psi0), len(psi0)),dtype=complex)
B = np.zeros((len(psi0), len(psi0)),dtype=complex)
for n in range(len(psi0)):
if n==0:
A[n,0]=a1
A[n,1]=a2
B[n,0]=b1
B[n,1]=b2
elif n==len(psi0)-1:
A[n,n]=a1
A[n,n-1]=a2
B[n,n]=b1
B[n,n-1]=b2
else:
A[n,n-1]=a2
A[n,n]=a1
A[n,n+1]=a2
B[n,n-1]=b2
B[n,n]=b1
B[n,n+1]=b2
#计算中间格点波函数
for n in range(1,len(t)):
phi[:,n]=np.linalg.inv(A) @ B @ phi[:,n-1]
#画图
for m in [0, 300, 500, 800, 2000]:
prob = phi[:,m]*np.conjugate(phi[:,m]) #计算概率分布
lbl='t='+str(m)
plt.figure(1,figsize=(15,3),dpi=300)
plt.plot(x,np.real(phi[:,m]),label=lbl)
plt.xlabel('x')
plt.ylabel('real')
plt.legend()
plt.grid()
plt.figure(2,figsize=(15,3),dpi=300)
plt.plot(x,np.imag(phi[:,m]),label=lbl)
plt.xlabel('x')
plt.ylabel('image')
plt.legend()
plt.grid()
plt.figure(3,figsize=(15,3),dpi=300)
plt.plot(x,np.sqrt(prob),label=lbl)
plt.xlabel('x')
plt.ylabel('probability')
plt.legend()
plt.grid()
输出:
最后需要坦白,跑代码的时间需要大概10分钟,但我还没研究更优化的算法和代码,嗯,因为懒。所以大家谨慎参考。
Like my work? Don't forget to support and clap, let me know that you are with me on the road of creation. Keep this enthusiasm together!
- Author