Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

What is my malfunction?

  1. May 25, 2015 #1
    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


    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. jcsd
  3. May 25, 2015 #2
    Also, note that state[0] is mass, state[1] is pressure, and den is density.
  4. May 25, 2015 #3


    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.
Share this great discussion with others via Reddit, Google+, Twitter, or Facebook