4

I am trying to numerically solve the Lane-Emden equation in Python using the scipy.integrate.ode class.

For some reason, my code works for integer values of n (the polytropic index) such as 3, but not for non-integer values such as 2.9 and 3.1.

I am getting a Runtime Warning, "7: RuntimeWarning: invalid value encountered in double_scalars" and the print statement "print f0, f1" shows that several values are nan's when n is not an integer.

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

What could be causing this issue?

user112829
  • 491
  • 1
  • 4
  • 7
  • 3
    Are the numbers you are attempting to exponentiate ever negative? If so, are you aware that raising a negative number to a non-integral power is undefined (or at least ventures into the realm of complex numbers)? – twalberg Jul 01 '15 at 18:35
  • Looks to be a numpy error but I don't use numpy or scipy. See [this answer](http://stackoverflow.com/a/3767475/2588818) for some possibly helpful advice in how to debug this. – Two-Bit Alchemist Jul 01 '15 at 18:40
  • 1
    @twalberg I just printed the bases and it looks like when they become very close to 0 they go negative. I don't understand why this would be, though, as I test whether r.y[0] is greater than 0 at the top of the while loop. Do you have any ideas? – user112829 Jul 01 '15 at 18:57
  • 3
    When you look at x,y you will see that, in particular, x (u[0]) does indeed go negative - which is not surprising since the algorithm does not know that it cannot go there. In most ODE algorithms, one jumps forwards a ways, looks back, and adjusts. If it jumps to far 'forward' it will try to evaluate at a negative u[0] and then all bets are off. You will have to be smarter about dealing with the cliff you have constructed. – Jon Custer Jul 01 '15 at 22:20

0 Answers0