CODEWITHSUNDEEP
pythonscipyode
Grim Reaper
Grim Reaper

Reputation: 89

Solving ODE with SciPy gives invalid value encountered in double_scalars error

I'm trying to solve the Lane-Emden equation for arbitrary values of n (polytropic index). In order to use SciPy, I have represented the second-order ODE as a set of two coupled first-order ODEs. I have the following 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

Here, I defined 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()

However, running this code results in the following error:

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

At first, I thought this was a result of the singularity of the equation at xi = 0, so I changed my integration interval:

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

This fixes the problem at n = 0., but not at other values of n. The plots I get don't make any sense and are just plain wrong.

Why do I get this error message, and how can I fix my code?

Upvotes: 2

Views: 1648

Answers (1)

p-robot
p-robot

Reputation: 4904

The following works for me (which looks similar to the Wikipedia entry). The derivative was written incorrectly (negative on the wrong element) and input to the derivative is changed simply to 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()

enter image description here

The original question used fractional values of the polytropic index, so something like follow could be used to investigate those values ... 

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

fractional polytropic index values

Upvotes: 1

Related Questions