# Scipy odeint Unexpected Results

I am attempting to numerically integrate the equation $$\frac{\mathrm{dP} }{\mathrm{d} r}=-\left ( P+\rho\left ( r \right ) \right )\frac{m\left ( r \right )+4\pi r^{3}P}{r\left [ r-2m\left ( r \right ) \right ]}$$ to find $P\left ( r \right )$. I have defined the following: $$\rho\left ( r \right )=\rho_{c}\left [ 1-\left ( \frac{r}{R} \right )^{2} \right ]$$ $$m\left ( r \right )=\frac{4\pi \rho_{c}\left ( 5r^{3}R^{2}-3r^{5} \right )}{15R^{2}}.$$ My code is as follows:

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

def dP_dr(x, r):
return -(x+density(r, p_c, R))*(mass(r, p_c, R)+4*np.pi*r**3*x)/(r*(r-
2*mass(r, p_c, R)))

def mass(r, p_c, R):
return 4*np.pi*p_c*(5*r**3*R**2-3*r**5)/(15*R**2)

def density(r, p_c, R):
return p_c*(1-(r/R)**2)

p_c = 7*10**17      #central density
P_c = 1.6*10**34    #central pressure

deltaR = 0.01
r = np.linspace(deltaR, R, int(R/deltaR))

P = odeint(dP_dr, P_c, r)

plt.figure()
plt.plot(r, P, label='Pressure')
plt.title('Neutron Star')
plt.legend(loc=0)
plt.show()


I'm unsure what is wrong, but I am getting an overflow error. The graphed result also is garbage. The pressure, $P\left ( r \right )$ should also decrease with $r$, but it jumps in very sharp spikes. Any help with what I am doing wrong would be appreciated.

• What is your initial condition? Do the various coefficients become singular at one end of your integration interval? – Wolfgang Bangerth Apr 17 '18 at 3:03
• My initial condition is P_c (central pressure). I am essentially using a model for a neutron star. The center of the star has pressure P_c and density p_c. I integrate this outward until i get to the radius of the star, R. At this point, the density and pressure should be zero. – Java Newbie Apr 17 '18 at 11:28
• Is this a boundary value problem then? Have you thought about using units that are more numerically friendly (or dimensionless variables)? – nicoguaro Apr 17 '18 at 13:12
• Besides numerical issues, are you sure that the stated equations are correct? Subtracting pressure from density or mass from radial coordinate seems fishy. – Bort Apr 17 '18 at 15:01

I have changed your code to this

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

def dP_dr(P, r, p_c, R, G, c):
return - G * (P + density(r, p_c, R) / c**2) * (mass(r, p_c, R)
+ 4 * np.pi * r**3 * P / c**2) / (r * (r - 2 * G* mass(r, p_c, R) / c**2))

def mass(r, p_c, R):
return 4 * np.pi * p_c * (5 * r**3 * R**2 - 3 * r**5) / (15 * R**2)

def density(r, p_c, R):
return p_c * (1 - (r/R)**2)

p_c = 7 * 10**10     #central density
P_c = 1.6 * 10**15    #central pressure
c = 3 * 10**8
G = 6.674 * 10**(-11)

deltaR = 0.0001
r = np.linspace(deltaR, R, int(R/deltaR))

P, dic = odeint(dP_dr, P_c, r, args=(p_c, R, G, c), full_output=True, printmessg=True)

num = 10000
plt.figure()
plt.plot(r[:num], P[:num], label='Pressure')

You have to play with num, p_c, P_c. With the numbers you show, $$P(r)$$ decays very fast. By the way, I have added two constants to the function dP_dr, $$G$$ and $$c$$.