项目场景：

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

模拟环境

Jupyter Notebook
python3(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∞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

−ℏ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−ℏ2dx2dΨ(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+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Ψ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 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 = 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 = 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)
#输出结果如下

一维无限深势阱定态薛定谔方程相关推荐

matlab二维势阱简谐振动程序,常规解法与MATLAB解决一维无限深势阱中的粒子问题...
龙源期刊网 http://www.doczj.com/doc/ddcba3222d60ddccda38376baf1ffc4ffe47e2cd.html 常规解法与MATLAB解决一维无限深势阱中的粒 ...
matlab求条件概率密度_Matlab对量子力学中的一维无限深势阱的模拟计算
小作品主要是借助 MATLAB 语言,对量子力学中一维无限深方势阱的波函数和概率密度分布做了可视化分析,通过对模拟图形的分析,形象的理解了波函数的特性,提高了自己使用计算机处理抽象物理问题的能力. 量 ...
matlab向量的模_基于MATLAB使用矩阵方法求解一维定态薛定谔方程
摘要:此文介绍了一种使用MATLAB求解一维定态薛定谔方程的方法.利用充分格式进行离散化,得出相应的矩阵方程,用MATLAB求解本征值和本征函数.此方法简单可靠,可以处理各种时间无关的束缚态问题.所用 ...
matlab求薛定谔方程,定态薛定谔方程的MATLAB求解(一)
ψλ:本征值为λ的本征函数,也有多个,甚至无穷多个,有时一个本征值对应多个不同的本征函数,这称为简并.若一个本征值对应的不同本征函数数目为N,则称N重简并. 1.4 定态情况下的薛定谔方程一般解 1. ...
影响堪比登月！谷歌等设计DL新方式让神经网络无限深无限窄
https://www.toutiao.com/a6656587016664252931/ 2019-02-11 11:57:19 [新智元导读]一个关于计算机如何学习的新理论的蓝图正在形成,其影响甚 ...
神经网络与定态薛定谔方程
定态薛定谔方程是如果是定态的自由粒子,这个方程的解是因为是定态的波函数与时间无关,这个粒子的能量E不随时间变化假设E=1,让t→0 所以波函数变成让A和都等于1 让神经网络里的节点都是在位形空 ...
com.sun.istack.SAXException2: 在对象图中检测到循环。这将产生无限深的 XML
错误如下所示: javax.xml.ws.soap.SOAPFaultException: Marshalling Error: 在对象图中检测到循环.这将产生无限深的 XML: org.entity ...
matlab解薛定谔方程,有限差分法解薛定谔方程与MATLAB实现
第30卷第3期高师理科学刊 Vol. 30 No.3 2010年 5月 Journal of Science of Teachers′ College and University Ma ...
物理化学笔记（1）量子化学基础
物理化学笔记(1)量子化学基础化学是人类关于原子和分子的知识和智慧的结晶.一个优秀的化学家需要适当了解其他科学分支的观测角度,在一定程度上听得懂其他学者的语言.但更重要和根本的是化学家一定要能流利地 ...
《高等统计物理学》Cookbook（持续更新）
在<高等统计物理学>系列文章中,有时候会用到一些重要的量子力学基础知识,比如定态薛定谔方程的求解.产生和消灭算符等等.本部分为读者建立一个即用即查的"工具箱",为方便正 ...

一维无限深势阱定态薛定谔方程