Scipy odeint Unexpected Results

Question

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
R = 10              #radius

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.xlabel('Radius (km)')
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.

Answer 1

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
R = 10              #radius
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')
plt.xlabel('Radius (km)')
plt.title('Neutron Star')
plt.legend()
plt.show()

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$.