如何在 python 上解决 TISE 的简单边值问题答案

【问题标题】：How to solve a simple boundary value problem for TISE on python如何在 python 上解决 TISE 的简单边值问题
【发布时间】：2020-05-17 10:26:11
【问题描述】：

我正在尝试在区间 [0,L] 上求解 TISE 以获得无限势阱 V=0。练习告诉我们，0 处的波函数及其导数的值分别为0,1。这允许我们使用scipy.integrate.odeint 函数来解决给定能量值的问题。

现在的任务是在给定L 处的波函数为0 的进一步边界条件下找到能量特征值，使用 python 上的求根函数。我做了一些研究，只能找到一种叫做“拍摄方法”的东西，我不知道如何实施。另外，我遇到了求解 BVP scipy 函数，但是我似乎无法理解该函数的第二个输入中到底发生了什么（边界条件残差）

m_el   = 9.1094e-31      # mass of electron in [kg]
hbar   = 1.0546e-34      # Planck's constant over 2 pi [Js]
e_el   = 1.6022e-19      # electron charge in [C]
L_bohr = 5.2918e-11      # Bohr radius [m]

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

def eqn(y, x, energy):              #array of first order ODE's     
    y0 = y[1]
    y1 = -2*m_el*energy*y[0]/hbar**2

    return np.array([y0,y1])

def solve(energy, func):           #use of odeint
    p0 = 0
    dp0 = 1
    x = np.linspace(0,L_bohr,1000)
    init = np.array([p0,dp0])
    ysolve = odeint(func, init, x, args=(energy,))
    return ysolve[-1,0]

这里的方法是在solve(energy,func)中输入eqn作为func。 L_bohr 是这个问题中的 L 值。我们正在尝试使用一些 scipy 方法在数值上找到能量特征值

【问题讨论】：

那么当你在下面添加一行 print(solve(0.1,eqn)) 或用 0.1 代替其他能量时会发生什么？您可能会在 odeint 文档 docs.scipy.org/doc/scipy/reference/generated/… 中找到一些线索
对于 1000 个 x 值直至 L，您得到了能量 = 0.1 的 ode 求解？但是这与最终的边界条件不匹配
那么您能更准确地陈述您的问题吗？您是否在某个地方遇到错误？
我需要使用一些 scipy 方法找到能量特征值
您是否需要以数值方式执行此操作，或者您可以利用无限井中薛定谔方程的解析解这一事实吗？

标签： python physics ode

【解决方案1】：

对于 scipy 中的所有其他求解器，参数顺序为 x,y，即使在 odeint 中，也可以通过提供选项 tfirst=True 来使用此顺序。于是改成

def eqn(x, y, energy):              #array of first order ODE's     
    y0, y1 = y
    y2 = -2*m_el*energy*y0/hbar**2

    return [y1,y2]

对于 BVP 求解器，您必须将能量参数视为具有零导数的额外状态分量，因此添加了第三个插槽在边界条件下。 Scipy 的solve_bvp 允许将其保留为参数，这样您就可以在边界条件中获得 3 个插槽，从而允许将一阶导数固定在 x=0 以从特征空间中选择一个非平凡解。

def bc(y0, yL, E):
    return [ y0[0], y0[1]-1, yL[0] ]

接下来构造一个接近可疑基态的初始状态并调用求解器

x0 = np.linspace(0,L_bohr,6);
y0 = [ x0*(1-x0/L_bohr), 1-2*x0/L_bohr ]
E0 = 134*e_el

sol = solve_bvp(eqn, bc, x0, y0, p=[E0])
print(sol.message, "  E=", sol.p[0]/e_el," eV")

然后产生情节

x = np.linspace(0,L_bohr,1000)
plt.plot(x/L_bohr, sol.sol(x)[0]/L_bohr,'-+', ms=1)
plt.grid()

The algorithm converged to the desired accuracy. E= 134.29310361903723 eV

【讨论】：