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

Question

我正在尝试求解区间 [0,L] 上无限势阱 V=0 的 TISE。练习告诉我们，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_bohr 是这个问题中的 L 值。我们正在尝试使用一些 scipy 方法从数值上找到能量特征值

Answer 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