如何在python上求解TISE的一个简单边值问题

Question

我试图求解区间[0,L]上无限位势阱V=0的问题。这个练习给出了波函数及其导数在0处的值分别是0,1。这允许我们使用scipy.integrate.odeint函数来解决给定能量值的问题。在

现在的任务是在进一步的边界条件下，即L处的波函数是0的进一步边界条件下，使用python上的寻根函数来找到能量本征值。我做了一些研究，只找到了一种叫做“射击法”的东西，我不知道该如何实施。另外，我遇到了solve-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值。我们试图用一些scipy方法数值求解能量本征值

Answer 1

对于scipy中的所有其他解算器，参数顺序x,y，甚至在{}中，可以通过给出选项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的bvp_solve允许将其作为参数保存， 这样你就可以

^{pr2}$

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

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