1
$\begingroup$

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.

$\endgroup$
  • 1
    $\begingroup$ What is your initial condition? Do the various coefficients become singular at one end of your integration interval? $\endgroup$ – Wolfgang Bangerth Apr 17 '18 at 3:03
  • $\begingroup$ 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. $\endgroup$ – Java Newbie Apr 17 '18 at 11:28
  • $\begingroup$ Is this a boundary value problem then? Have you thought about using units that are more numerically friendly (or dimensionless variables)? $\endgroup$ – nicoguaro Apr 17 '18 at 13:12
  • 6
    $\begingroup$ Besides numerical issues, are you sure that the stated equations are correct? Subtracting pressure from density or mass from radial coordinate seems fishy. $\endgroup$ – Bort Apr 17 '18 at 15:01
1
$\begingroup$

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

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.