Das Lösen von ODE mit SciPy ergibt einen ungültigen Wert in double_scalars Fehler

Question

Ich versuche, die Lane-Emden-Gleichung für beliebige Werte von n zu lösen (polytropischer Index). Um SciPy zu verwenden, habe ich die ODE zweiter Ordnung als eine Menge von zwei gekoppelten ODEs erster Ordnung dargestellt. Ich habe den folgenden Code:

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 

Hier habe ich definiert phi = theta‘

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() 

jedoch in dem folgenden Fehler dieser Code Ergebnisse ausgeführt wird:

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

Zuerst habe ich dachte, dies sei ein Ergebnis der Singularität der Gleichung bei xi = 0, so änderte ich mein Integrationsintervall:

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

Dies behebt das Problem bei n = 0., aber nicht bei anderen Werten von n. Die Plots, die ich bekomme, ergeben keinen Sinn und sind einfach falsch.

Warum bekomme ich diese Fehlermeldung, und wie kann ich meinen Code beheben?

Answer 1

Das folgende funktioniert für mich (das ähnelt der Wikipedia entry). Die Ableitung wurde falsch geschrieben (negativ für das falsche Element) und die Eingabe in die Ableitung wird einfach in n geändert.

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() 

---

Die ursprüngliche Frage Bruchwerte des polytropischer Index verwendet, so etwas wie Folge verwendet werden könnten, um diese Werte ...

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]) 

zu untersuchen