Лишняя работа, выполняемая при вызове (возможно, неправильный тип Dfun) с помощью scipy.integrate.odeint ⇐ Python

[Anonymous]

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

Я использую следующую замену уравнения Шрёдингера:

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, я найду энергию n=2.

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

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