# What is my malfunction?

1. May 25, 2015

### Bishop556

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:
Code (Text):
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 + dr

plt. 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: May 25, 2015
2. May 25, 2015

### Bishop556

Also, note that state[0] is mass, state[1] is pressure, and den is density.

3. May 25, 2015

### Staff: Mentor

I'm not sure about this, but it seems to me that you defined your state array two be two-dimensional (your code: state = np.zeros([N,2]) ), but you are using it later as if it were a one-dimensional array, as in what you wrote above.