待命完成的多余工作（可能是错误的Dfun类型）scipy.integrate.odein公司

Question

我正在用Python编写一个预测氢能级的代码，我将用它作为研究夸克偶素能级的模板。我用scipy.integrate.odeint()函数来解Shroedinger方程，它对低能级直到n=6起作用。我不希望我有太多的需要去超越它，但是odeint返回Excess work done on this call (perhaps wrong Dfun type).，这只会鼓励我扩展我可以预测的范围。在

我使用的Shroedinger方程代换是：

u'' - (l*(l+1)/r**2 - 2mu_e(E-V_emag(r))) * u = 0
=>
u' = v
v' = ((l*(l+1))/(r**2) - 2.0*mu_e*(E - V_emag(r)))*u

然后我在上面使用scipy.integrate.odeint()，遍历能量，并使用我定义的其他函数来评估结果中的转折点和节点。我找到能级的方法是找到可能的最低E值，其中转折点和节点的数目与它应该匹配；然后将L递增1，找到新的基态能量，例如，如果L=0我将找到n=1能量，如果L=3，我将找到{}能量。在

一旦代码增加到L=7，它就不会返回任何有用的信息。r的范围已经扩展，但是我尝试保持不变以减少步骤数，但是没有效果。代码是自学的，所以在我的研究中我读到了关于雅各比的书。恐怕我还没有弄清楚它们是什么，或者我是否需要。有什么想法吗？在

^{pr2}$

Answer 1

我不知道你的代码到底出了什么问题，也不完全确定你的while i<40:循环在做什么，所以如果我错了，也许你可以纠正以下问题。在

如果你想要这个系统中某个n, l的波函数，你可以计算能量为E=RH/n^2，其中RH是里德堡常数，所以你不需要计算节点。如果您确实需要计算节点，那么(n,l)对应的数量是n-l-1，因此您可以改变E并观察固定l的节点数量的变化

在我看来，主要的问题是，你的r范围不够大，不能包含所有的节点（对于大的n~l），而且odeint不知道要远离另一个渐近解psi~exp（+cr），因此在某些情况下，对于大r，^{

如果这有帮助的话，这就是我找到SE方程的数值解的方法：你需要根据n,l改变r-范围来避免上述问题（例如，如果你要求n, l = 10, 9，看看会发生什么）。在

import numpy as np
import scipy as sp
from scipy.integrate import odeint

m_e, m_p, hbar = sp.constants.m_e, sp.constants.m_p, sp.constants.hbar
mu_e = m_e*m_p/(m_e + m_p)
bohr = sp.constants.physical_constants['Bohr radius'][0]
Rinfhc = sp.constants.physical_constants['Rydberg constant times hc in J'][0]
RHhc = Rinfhc * mu_e / m_e

fac = sp.constants.e**2/4/sp.pi/sp.constants.epsilon_0

def V(r):
    return -fac/r

def deriv(y, r, l, E):
    y1, y2 = y
    dy1dr = y2
    dy2dr = -2*y2/r - (2*mu_e/hbar**2*(E - V(r)) - l*(l+1)/r**2)*y1
    return dy1dr, dy2dr

def solveSE(l, E, y0):
    rstep = 0.001 * bohr
    rmin = rstep
    rmax = 200*l * bohr         # 
    r = np.arange(rmin, rmax, rstep)
    y, dydt = odeint(deriv, y0, r, args=(l,E)).T
    return r, y, dydt

n = 10
l = 2
y0 = (bohr, -bohr)
E = -RHhc / n**2
r, psi, dpsi_dr = solveSE(l, E, y0)

import pylab
pylab.plot(r, psi)
pylab.show()