项目场景:
本文介绍使用龙格-库塔法(Runge-Kutta method)和弦截法(Secant method)求解一维无限深势阱的定态薛定谔方程的本征值(eigenvalue)和本征函数( eigenfunction ),并进行可视化。
模拟环境
- Jupyter Notebook
- python3(numpy,matplotlib)
理论基础:
1.薛定谔方程
1.力场中粒子的薛定谔方程( Schrödinger equation)
−
ℏ
2
2
m
∇
2
Ψ
+
V
(
r
⃗
)
Ψ
=
E
Ψ
(
1
)
\bf{-\hbar^2 \over 2m}\nabla^2\Psi+V(\vec{r})\Psi=E\Psi (1)
2m−ℏ2∇2Ψ+V(r)Ψ=EΨ(1)
2.根据(1)式化简得一维无限深势阱中粒子的薛定谔方程(a Particle in an Infinite Potential Well)
V
(
x
)
=
{
0
if
−
a
<
x
<
a
∞
if
x
>
a
,
x
<
−
a
\bf V(x) = \begin{cases} 0 &\text{if } -a<x <a\\ \infin &\text{if } x >a , x<-a \end{cases}
V(x)={0∞if −a<x<aif x>a,x<−a
− ℏ 2 2 m d Ψ ( x ) 2 d x 2 + V ( x ) Ψ ( x ) = E Ψ ( x ) ( 2 ) \bf {-\hbar^2 \over 2m}\dfrac{d\Psi(x)^2}{dx^2}+V(x)\Psi(x)=E\Psi(x)(2) 2m−ℏ2dx2dΨ(x)2+V(x)Ψ(x)=EΨ(x)(2)
2.龙格-库塔法
1.龙格-库塔法(Runge-Kutta method)
四阶龙格-库塔法是常用的常微分方程数值计算方法,局部截断误差为五阶小量,计算精度相当高。
- 对于一阶微分方程,若在给定区间进行划分,第(i+1)个点的函数值
y
i
+
1
y_{i+1}
yi+1与第i个点的函数值
y
i
y_{i}
yi的有下面的关系
(h为划分的间距,f(x,y)为所求函数的一阶导)
{ d y d x = f ( x , y ) k 1 = f ( x i , y i ) k 2 = f ( x i + 1 2 h , y i + k 1 1 2 h ) k 3 = f ( x i + 1 2 h , y i + k 2 1 2 h ) k 3 = f ( x i + h , y i + k 3 h ) y i + 1 = y i + 1 6 h ( k 1 + 2 k 2 + 2 k 3 + k 4 ) \\ \textcolor{blue}{ \begin{cases} \dfrac{dy}{dx}=f(x,y) \\ k_{1} = f(x_i,y_i) \\ k_{2} = f(x_i+{1 \over 2}h,y_i+k_1{1 \over 2}h) \\ k_{3} = f(x_i+{1 \over 2}h,y_i+k_2{1 \over 2}h) \\ k_{3} = f(x_i+h,y_i+k_3h) \\ y_{i+1} = y_{i} + {1 \over 6}h(k_{1}+2k_{2}+2k_{3}+k_{4}) \end{cases} } ⎩⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎧dxdy=f(x,y)k1=f(xi,yi)k2=f(xi+21h,yi+k121h)k3=f(xi+21h,yi+k221h)k3=f(xi+h,yi+k3h)yi+1=yi+61h(k1+2k2+2k3+k4) - 对于二阶微分方程,需要引入一个中间量(φ),将二阶微分方程化为两个一阶微分方程组(如下所示),这样,Ψ的增量可以用φ表示,φ的增量可以用题目中的二阶微分方程表示,给定φ和Ψ的初始值φ(0),Ψ(0),代入下面的微分方程组得到各自的导数(变化率)得到相应k1,k2,k3,k4,然后乘以h(分割的间距),便可以得到φ(1),Ψ(1),以此类推便可以得到所求区间的所有φ,Ψ的值
{ d Ψ d x = ϕ d ϕ d x = 2 m ℏ 2 ( V ( X ) − E ) Ψ ( x ) \begin{cases} \dfrac{d\Psi}{dx}=\phi \\ \\ \dfrac{d\phi}{dx}={2m \over \hbar^2}(V(X)-E)\Psi(x) \end{cases} ⎩⎪⎪⎪⎨⎪⎪⎪⎧dxdΨ=ϕdxdϕ=ℏ22m(V(X)−E)Ψ(x)
程序实现:
import numpy as np
import matplotlib.pyplot as plt
#预先准备好所需物理量
electron_mass = 9.109383702e-31 #electron mass (kg)
h_bar = 1.054571817e-34 #h bar (J.s)
electron_charge = 1.602176634e-19 #electron charge (C)
#定义模拟的区域宽度
a = 5e-11
xstart = -a #(m)
xend = a #(m)
N = 1000
h = (2*a) / N
x_points=np.arange(xstart,xend,h) #划分区间,获得离散的点
r = np.array([0,1])
#定义势能函数
def V(x):
return 0.0
#此函数以上面的微分方程组相匹配
#用于给定φ,Ψ获得相应的导数值,返回一个2*1数组
def function(r,x,E,Potential):
psi = r[0]
phi = r[1]
fpsi = phi
fphi = 2*electron_mass*(1/h_bar**2)*(Potential(x)-E)*psi
return np.array([fpsi,fphi])
#龙格-库塔法函数实现
def RungeKutta2d(r,x_points,function,E,Potential):
xpoints = []
ypoints = []
for t in x_points:
xpoints.append(r[0])
ypoints.append(r[1])
k1 = h*function(r,t,E,Potential)
k2 = h*function(r+0.5*k1, t+0.5*h,E,Potential)
k3 = h*function(r+0.5*k2, t+0.5*h,E,Potential)
k4 = h*function(r+k3, t+h,E,Potential)
r = r + (k1 + 2*k2 + 2*k3 + k4)/6
xpoints.append(r[0])
ypoints.append(r[1])
return np.array([xpoints, ypoints])
#此函数用于获得理论解用于和数值解对比
def get_Analytical(n):
E_n = (np.pi**2 * h_bar**2 * n**2) / (2 * electron_mass * (2*a)**2)
return E_n
#此函数用于可视化
def draw_Image(WaveFunction,ProbablityDensity):
x_points2 = np.arange(xstart,xend+h,h)
plt.figure(figsize=(10, 4))
plt.subplot(1,2,1)
#波函数图像绘制
plt.title('Wavefuction')
plt.plot(x_points2,WaveFunction,'r')
plt.subplot(1,2,2)
#概率密度图像绘制
plt.title('ProbabilityDensity(Ψ^2)')
plt.plot(x_points2,(1/a)*np.square(ProbablityDensity),'g')
def get_Calculated(E1,E2,n):
wave1 = RungeKutta2d(r,x_points,function,E1,V)[0,N]
wave2 = RungeKutta2d(r,x_points,function,E2,V)[0,N]
tolerance = electron_charge / 1000
#这里使用的是弦割法,见下面说明
while abs(E2-E1) > tolerance:
E3 = E2 - wave2*(E2-E1)/(wave2-wave1)
E1 = E2
E2 = E3
wave1 = RungeKutta2d(r,x_points,function,E1,V)[0,N]
wave2 = RungeKutta2d(r,x_points,function,E2,V)[0,N]
solutionE = RungeKutta2d(r,x_points,function,E3,V)
E_n = get_Analytical(n)
print("理论解{0:0.9e} J".format(E_n))
print("数值解 {0:0.9e} J".format(E3))
draw_Image(solutionE[0],solutionE[0])
我们要明确一点,对这个物理问题的解是唯一的(给定初始条件求解波函数,粒子的状态,本征值,本征函数必然是唯一确定的),但我们使用龙格库塔法获得波函数的数值解并不是唯一的,这是因为这是一个线性微分方程,即使我们给定初始条件,也未必是我们想要的真实粒子的波函数(因为我们的给定的初始条件不一定是粒子实际运动的初始状态),所以不妨我们先通过理论计算获得解析解,然后给出一个大致范围(如E1,E2,理论值要处于E1和E2之间),利用弦割法的原理,在这个E1和E2之间寻找符合真实粒子状态的本征值和本征函数。
#以下是各个能级态的能量本征值对应的区间,不唯一,仅作为参考
n = 1
E1 = 3e-18
E2 = 8e-18
# n = 2
# E1 = 2e-17
# E2 = 3e-17
# n = 3
# E1 = 5e-17
# E2 = 6e-17
# n = 4
# E1 = 9e-17
# E2 = 10e-17
# n = 5
# E1 = 1e-16
# E2 = 2e-16
get_Calculated(E1,E2,n)
#输出结果如下