Избыточная работа, выполненная по вызову (возможно, неправильный тип dfun) с scipy.integrate.odeint ⇐ Python

[Anonymous]

Я пишу код в Python, который предсказывает уровень энергии водорода, который я буду использовать в качестве шаблона для исследования уровней энергии Quarkonium. Я использую функцию scipy.integrate.odeint () для решения уравнения Shroedinger, и она работает нормально для более низких уровней энергии до n = 6 . Я не ожидаю, что мне нужно выйти за рамки этого, но Odeint возвращает избыточную работу, выполненную по этому вызову (возможно, неправильный тип dfun). , что только побуждает меня расширить то, что я могу предсказать.

Замена уравнения Шродингера, которую я использую: < /p>

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
< /code>

Затем я использую scipy.integrate.odeint () < /code> на нем и итерации через энергии и используя другие функции, которые я определил, что оценивает поворотные точки и узлы в результате. То, как я нахожу уровни энергии, - это найти самое низкое значение E, где количество поворотных точек и узлов соответствует тому, что должно быть; Затем увеличение l на 1 и поиск новой энергии земли, например, Если l = 0 я найду n = 1 энергия и если l = 3 , я найду n = 2 энергия.

Как только код увеличивает до L = 7 < /code>, он не возвращает ничего полезного. Диапазон r был расширен, но я пытался сохранить его таким же, чтобы уменьшить количество шагов, но безрезультатно. Код самоучка, поэтому в своем исследовании я читал о якобианских. Боюсь, я еще не выработал, что они есть или если мне это нужно. Есть идеи? < /P>

def v_emag(r):
v = -alpha/r
return v

def s_e(y,r,l,E): #Shroedinger equation for electromagntism
x = numpy.zeros_like(y)
x[0] = y[1]
x[1] = ((l*(l+1))/(r**2) - 2.0*mu_e*(E - V_emag(r)))*y[0]
return x

def i_s_e(l,E,start=0.001,stop=None,step=(0.005*bohr)):
if stop is None:
stop = ((l+1)*30-10)*bohr
r = numpy.arange(start,stop,step)
y = odeint(s_e,y0,r,args=(l,E))
return y

def inormalise_e(l,E,start=0.001,stop=None,step=(0.005*bohr)):
if stop is None:
stop = ((l+1)*30-10)*bohr
r = numpy.arange(start,stop,step)
f = i_s_e(l,E,start,stop,step)[:,0]
f2 = f**2
area = numpy.trapz(f2,x=r)
return f/(numpy.sqrt(area))

def inodes_e(l,E,start=0.001,stop=None,step=(0.005*bohr)):
if stop is None:
stop = ((l+1)*30-10)*bohr
x = i_s_e(l,E,start,stop,step)
r = numpy.arange(start,stop,step)
k=0
for i in range(len(r)-1):
if x[i,0]*x[i+1,0] < 0: #If u value times next u value

Подробнее здесь: https://stackoverflow.com/questions/275 ... ate-odeint