如何在python上解决TISE的简单边值问题 [英] How to solve a simple boundary value problem for TISE on python

问题描述

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

I am trying to solve the TISE for an infinite potential well V=0 on the interval [0,L]. The exercise gives us that the value of the wavefunction and its derivative at 0 is 0,1 respectively. This allows us to using the scipy.integrate.odeint function in order to solve the problem for a given energy value.

现在的任务是使用python上的根查找函数，在给定进一步边界条件(在L处的波动函数为0)的情况下找到能量特征值.我进行了一些研究，只能找到一种叫做射击方法"的东西，但我无法弄清楚如何实现.另外，我遇​​到了solve BVP scipy函数，但是我似乎无法理解该函数的第二个输入(边界条件残差)到底是什么

The task is to now find the energy eigenvalues given the further boundary condition that the wavefunction at L is 0, using a root finding function on python. I have done some research and could only find something called the 'shooting method' which I cannot figure out how to implement. Also, I have come across the solve BVP scipy function, however I can't seem to understand what exactly goes in the second input for this function (boundary condition residuals)

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]

此处的方法是将eqn作为func输入到solve(energy，func)中. L_bohr是此问题中的L值.我们正在尝试使用某种scipy方法在数值上找到能量本征值

The method here is to input eqn as func in solve(energy,func). L_bohr is the L value in this problem. We are trying to numerically find the energy eigenvalues using some scipy method

推荐答案

对于scipy中所有其他求解器，参数顺序x,y，甚至在odeint中，也可以通过提供选项tfirst=True来使用此顺序.因此更改为

For all the other solvers in scipy the argument order x,y, and even in odeint one can use this order by giving the option tfirst=True. Thus change to

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处，从而从特征空间中选择一个非平凡解.

For the BVP solver you have to think of the energy parameter as an extra state component with zero derivative, thus adding a third slot in the boundary conditions. Scipy's solve_bvp allows to keep it as parameter, so that you get 3 slots in the boundary conditions, allowing to fix the first derivative at x=0 to select one non-trivial solution from the eigenspace.

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

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

Next construct an initial state that is close to the suspected ground state and call the solver

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

这篇关于如何在python上解决TISE的简单边值问题的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！