- #1
Bishop556
- 37
- 4
Hello, I am trying to now use my Runge Kutta code to numerically solve the equations of state for an ordinary star and then transcribe that to a neutron star. My code is below:
What I am trying to do is solve the two equations of state while using the polytropic relation between density and pressure.Now, I am slightly confused on how I can accomplish this as I have three unknowns: mass, pressure, and density which requires three separate equations are needed to develop a unique solution. I, however, am only receiving a straight line that diverges to infinity on the second iteration. Does anyone have any ideas why this is occurring?
Code:
import numpy as np
import matplotlib.pyplot as plt
from math import pi
#Constants
G = 6.674*1e-11
N = 5000
Gamma = 1.5
K = 5000
dur = 5
dr = dur/float(N-1)
state = np.zeros([N,2])
#initial conditions
state[0,0] = 0
state[0,1] = 5
def rk4(y,r, dr, den, deriv):
k1 = dr * deriv(y, den, r)
k2 = dr * deriv(y + 0.5*k1,den, r + 0.5*dr)
k3 = dr * deriv(y + 0.5*k2,den, r + 0.5*dr)
k4 = dr * deriv(y + k3, den, r + dr)
y_next = y + (k1 + 2*(k2+k3)+k4)/6
return y_next
def density(state):
return ((1/K)*state[1])**(-Gamma) # row = (P/K)^(-Gamma)
def EoS(state,den, r):
y0 = 4*pi*(r**2)*den #dM/dr = 4*pi*(r^2)*row
y1 = (-G*state[1]*den)/(r**2) #dP/dr = -G*M*row / r^2
return np.array([y0,y1])
r = 1e-23 #starting value to avoid singularity
for i in range(N-1):
den = density(state[i])
state[i+1]= rk4(state[i],r,dr,den, EoS)
r = r + drplt. plot(state[:,1], state[:,0])
plt. show()
What I am trying to do is solve the two equations of state while using the polytropic relation between density and pressure.Now, I am slightly confused on how I can accomplish this as I have three unknowns: mass, pressure, and density which requires three separate equations are needed to develop a unique solution. I, however, am only receiving a straight line that diverges to infinity on the second iteration. Does anyone have any ideas why this is occurring?
Last edited: