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.