What is the semantic error in my code

[Bishop556]

Note that the code I am using is Python. Currently, I am attempting to numerically solve the equations of state that give rise to the Chandrasekar limit for a white dwarf star. My code works for a simple harmonic oscillator given the correct equations of motion, but does not compute the correct curve for my two coupled ODEs that I need for my project. My code is below. Note that the equations of state are defined in the comments.

import numpy as np
import matplotlib.pyplot as plt
from math import sqrt

#Variable Definitions

N = 500 #number of entries
dur = 25000 #radius taken up to in kilometers
dr = dur / float(N-1) #step size

#creating the array

y = np.zeros([N,2])

y[0,0] = 20 #initial central momentum
y[0,1] = 0  #initial mass

#runge kutta algorithm

def rk4(y, r, dr, deriv):
    k1 = dr * deriv(y,r)
    k2 = dr * deriv(y + 0.5*k1,r + 0.5*dr)
    k3 = dr * deriv(y + 0.5*k2, r + 0.5*dr)
    k4 = dr * deriv(y + k3, r + dr)
    y_next = y + (k1 + 2*(k2+k3)+k4)/6
    return y_next#Equations of state (Note that x[0] is my fermi momentum x and x[1] is my mass. r is my radius

def Eqns(x,r):
    y0 = (-5*x[1]*sqrt(1+x[0]**2))/(3*x[0]*r**2)  #dx/dr = -5/3(M/r^2)*sqrt(1+x^2)/x
    y1 = 3*(r**2)*(x[0]**3)                       #dM/dr = 3*r^2*x^3
    return np.array([y0,y1])

#forming data points

for i in range(N-1):
    y[i+1] = rk4(y[i],10**(-23), dr, Eqns) #started at some arbitrary smalll value to avoid
                                           #singularity
radius = np.linspace(10**(-23),25000, N)

#plotting
plt.plot(y[:,1],radius)

plt.show()

For some reason when I run this all I receive is a straight linear line which is nothing like the curve I expected that is attached to this post. Also, the graph indicates my mass is negative which is ridiculous. I have a feeling the error is a semantic error somewhere in my code, but would I not expect a curve instead of a straight line for ODEs that are obviously not linear. Thanks again.

EDIT: I forgot to mention that what the curve should resemble is the curve in the relativistic fermi gas line.

 

[D H]

Your bug is here:

y[i+1]= rk4(y[i],10**(-23), dr, Eqns)

You have r frozen at 10^-23. Change that 10**(-23) to i*r if i else 1e-23

Alternatively, change your loop to

r = 1e-23 # To avoid singularity at r=0
for i in range(N-1):
    y[i+1]= rk4(y[i],r, dr, Eqns)
    r = r+dr

 

[Bishop556]

After applying this change, I am receiving a straight horizontal line which is odd.

 

[SteamKing]

After applying this change, I am receiving a straight horizontal line which is odd.

What's the y location of this straight horizontal line?

 

[Bishop556]

import numpy as np
import matplotlib.pyplot as plt
from math import sqrt

#Variable Definitions

N = 500
dur = 2.0
 dr = dur/float(N-1)

#creating the array

y = np.zeros([N,2])

y[0,0] = 15
y[0,1] = 0#runge kutta algorithm

def rk4(y, r, dr, deriv):
    k1 = dr * deriv(y,r)
    k2 = dr * deriv(y + 0.5*k1,r + 0.5*dr)
    k3 = dr * deriv(y + 0.5*k2, r + 0.5*dr)
    k4 = dr * deriv(y + k3, r + dr)
    y_next = y + (k1 + 2*(k2+k3)+k4)/6
    return y_next#Equations of state

def Eqns(x,r):
    y0 = (-5*x[1]*sqrt(1+x[0]**2))/(3*x[0]*r**2)  #dx/dr = -5/3(M/r^2)*sqrt(1+x^2)/x
    y1 = 3*(r**2)*(x[0]**3)                       #dM/dr = 3*r^2*x^3
    return np.array([y0,y1])

#forming data pointsr = 1e-23 #starting value to avoid singularity
for i in range(N-1):
    if y[i,0]>0:
        y[i+1] = rk4(y[i],r,dr,Eqns)
        r = r + dr
      
    else:
        break #meant to stop negative values of x to be in the array

  
  
  
  #started at some arbitrary smalll value to avoid
                                           #singularity
radius = np.linspace(1e-23,r,N)#plotting
plt.xlim([0,2])
plt.ylim([0,.535])
plt.plot(radius, y[:,1])

plt.show()

 

[Bishop556]

My new code is listed above. I added a portion that would stop the array from receiving an negative values of x. The rest of the array will be zeros due to this effect which is expected. When computing the arrays, the data points seem to be somewhat correct, but the graph is wonky. My current abomination is linked here.

 

[Bishop556]

What's the y location of this straight horizontal line?

I am plotting radius in the x-axis and Mass in y axis.

 

[D H]

Do you have the right equations of state? Do you need to take smaller steps?

 

[Bishop556]

Do you have the right equations of state? Do you need to take smaller steps?

I believe I do. The equations of state are mentioned as comments in my code. Do you notice anything I am doing wrong to code them correctly?