What is the semantic error in my code

Tags:
1. May 11, 2015

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.

Code (Python):

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

#plotting

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.

Attached Files:

File size:
13.3 KB
Views:
78
• graph from code.png
File size:
3.4 KB
Views:
66
Last edited by a moderator: May 11, 2015
2. May 11, 2015

D H

Staff Emeritus
Code (Text):
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

Code (Text):

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

Last edited: May 11, 2015
3. May 11, 2015

Bishop556

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

4. May 11, 2015

SteamKing

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

5. May 12, 2015

Bishop556

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

r = 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

#plotting
plt.xlim([0,2])
plt.ylim([0,.535])

plt.show()

6. May 12, 2015

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.

Attached Files:

• graph from code.png
File size:
2.4 KB
Views:
74
7. May 12, 2015

Bishop556

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

8. May 12, 2015

D H

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

9. May 12, 2015

Bishop556

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?