Решение ODE с помощью SciPy дает недопустимое значение, обнаруженное в ошибке double_scalars

Question

Я пытаюсь решить уравнение Лейн-Эмдена для произвольных значений n (политропный индекс). Чтобы использовать SciPy, я представил ODE второго порядка как набор из двух связанных ODE первого порядка. У меня есть следующий код:

import matplotlib.pyplot as plt 
import numpy as np 
from scipy.integrate import odeint 

def poly(y,xi,n): 
    theta, phi = y 
    dydt = [phi/(xi**2), -(xi**2)*theta**n] 
    return dydt 

Здесь я определил Фи = тета»

y0 = [1.,0.] 
xi = np.linspace(0., 16., 201)  
for n in range(0,11): 
    sol = odeint(poly, y0, xi, args=(n/2.,)) 
    plt.plot(xi, sol[:, 0], label=str(n/2.)) 
plt.legend(loc='best') 
plt.xlabel('xi') 
plt.grid() 
plt.show() 

Однако выполнение этого кода приводит к следующей ошибке:

RuntimeWarning: invalid value encountered in double_scalars
app.launch_new_instance()

На первый, я что это было связано с особенностью уравнения при xi = 0, поэтому я изменил интервал интеграции:

xi = np.linspace(1e-10, 16., 201) 

Это устраняет проблему при n = 0., но не при других значениях n. Графики, которые я получаю, не имеют никакого смысла и просто ошибаются.

Зачем мне это сообщение об ошибке и как я могу исправить свой код?

Answer 1

Следующие работы для меня (что похоже на Wikipedia entry). Производная была написана неправильно (отрицательная на неправильном элементе), и ввод в производную изменен просто на n.

import matplotlib.pyplot as plt 
import numpy as np 
from scipy.integrate import odeint 

def poly(y,xi,n): 
    theta, phi = y 
    dydxi = [-phi/(xi**2), (xi**2)*theta**n] 
    return dydxi 

fig, ax = plt.subplots() 

y0 = [1.,0.] 
xi = np.linspace(10e-4, 16., 201) 

for n in range(6): 
    sol = odeint(poly, y0, xi, args=(n,)) 
    ax.plot(xi, sol[:, 0]) 

ax.axhline(y = 0.0, linestyle = '--', color = 'k') 
ax.set_xlim([0, 10]) 
ax.set_ylim([-0.5, 1.0]) 
ax.set_xlabel(r"$\xi$") 
ax.set_ylabel(r"$\theta$") 
plt.show() 

---

Оригинальный вопрос используется дробные значения показателя политропы, так что-то вроде последующих могут быть использованы для изучения этих значений ...

for n in range(12): 
    sol = odeint(poly, y0, xi, args=(n/2.,)) 

    if np.isnan(sol[:,1]).any(): 
     stopper = np.where(np.isnan(sol[:,1]))[0][0] 
     ax.plot(xi[:stopper], sol[:stopper, 0]) 
    else: 
     ax.plot(xi, sol[:, 0])