项目场景：

本文介绍使用龙格－库塔法（Runge-Kutta method）和弦截法（Secant method）求解一维无限深势阱的定态薛定谔方程的本征值（eigenvalue）和本征函数（ eigenfunction ），并进行可视化。

模拟环境

Jupyter Notebookpython3(numpy,matplotlib)

理论基础：

1.薛定谔方程

1.力场中粒子的薛定谔方程（ Schrödinger equation）

−ℏ22m∇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)={0if−a<x<a∞ifx>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<aifx>a,x<−a​

−ℏ22mdΨ(x)2dx2+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−ℏ2​dx2dΨ(x)2​+V(x)Ψ(x)=EΨ(x)（2）

2.龙格－库塔法

1.龙格－库塔法（Runge-Kutta method）

四阶龙格－库塔法是常用的常微分方程数值计算方法，局部截断误差为五阶小量，计算精度相当高。

对于一阶微分方程，若在给定区间进行划分，第（i+1）个点的函数值yi+1y_{i+1}yi+1​与第i个点的函数值yiy_{i}yi​的有下面的关系

(h为划分的间距，f（x,y）为所求函数的一阶导)

{dydx=f(x,y)k1=f(xi,yi)k2=f(xi+12h,yi+k112h)k3=f(xi+12h,yi+k212h)k3=f(xi+h,yi+k3h)yi+1=yi+16h(k1+2k2+2k3+k4)\\ \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​+21​h,yi​+k1​21​h)k3​=f(xi​+21​h,yi​+k2​21​h)k3​=f(xi​+h,yi​+k3​h)yi+1​=yi​+61​h(k1​+2k2​+2k3​+k4​)​对于二阶微分方程，需要引入一个中间量（φ），将二阶微分方程化为两个一阶微分方程组(如下所示)，这样，Ψ的增量可以用φ表示，φ的增量可以用题目中的二阶微分方程表示，给定φ和Ψ的初始值φ（0），Ψ（0），代入下面的微分方程组得到各自的导数（变化率）得到相应k1,k2,k3,k4，然后乘以h（分割的间距），便可以得到φ（1），Ψ（1），以此类推便可以得到所求区间的所有φ，Ψ的值

{dΨdx=ϕdϕdx=2mℏ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 npimport 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 = 1000h = (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 = phifphi = 2*electron_mass*(1/h_bar**2)*(Potential(x)-E)*psireturn 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)/6xpoints.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 = E2E2 = E3wave1 = 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 = 1E1 = 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-16get_Calculated(E1,E2,n)#输出结果如下