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:(adsbygoogle = window.adsbygoogle || []).push({});

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 (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()

**Physics Forums | Science Articles, Homework Help, Discussion**

The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

# What is my malfunction?

**Physics Forums | Science Articles, Homework Help, Discussion**