Я пытаюсь численно решить уравнение Лейна-Эмдена в Python, используя scipy.integrate.ode.
По какой-то причине мой код работает с целыми значениями n (индекс политропы), такими как 3, но не работает с нецелыми значениями, такими как 2,9 и 3,1.
Я получаю предупреждение во время выполнения: «7: RuntimeWarning: недопустимое значение встречается в double_scalars», а оператор печати «print f0, f1» показывает, что несколько значений являются nan, когда n не является целым числом.
from scipy.integrate import ode
import numpy as np
def rhs(zet, u, n):
x = u[0] # theta
y = u[1] # phi
f0 = -u[1]/(zet**2)
f1 = (u[0]**(n))*(zet**2)
print f0, f1
return np.array([f0,f1])
r = ode(rhs).set_integrator("vode", method="adams", atol=1.e-13, rtol=1.e-13, nsteps=15000, order=12)
y0 = np.array([1.,0.])
zeta0 = 0.000001
r.set_initial_value(y0, zeta0)
n=3.1
r.set_f_params(n)
dzeta = 0.000001
zeta = [zeta0]
theta = [y0[0]]
while r.successful() and r.y[0] > 0.0:
r.integrate(r.t + dzeta)
zeta.append(r.t)
theta.append(r.y[0])
Что может быть причиной этой проблемы?
