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 the semantic error in my code

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

    Attached Files:

    Last edited by a moderator: May 11, 2015
  2. jcsd
  3. May 11, 2015 #2

    D H

    User Avatar
    Staff Emeritus
    Science Advisor

    Your bug is here:
    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

    Alternatively, change your loop to
    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
  4. May 11, 2015 #3
    After applying this change, I am receiving a straight horizontal line which is odd.
     
  5. May 11, 2015 #4

    SteamKing

    User Avatar
    Staff Emeritus
    Science Advisor
    Homework Helper

    What's the y location of this straight horizontal line?
     
  6. May 12, 2015 #5
    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
    radius = np.linspace(1e-23,r,N)


    #plotting
    plt.xlim([0,2])
    plt.ylim([0,.535])
    plt.plot(radius, y[:,1])

    plt.show()
     
     
  7. May 12, 2015 #6
    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:

  8. May 12, 2015 #7
    I am plotting radius in the x axis and Mass in y axis.
     
  9. May 12, 2015 #8

    D H

    User Avatar
    Staff Emeritus
    Science Advisor

    Do you have the right equations of state? Do you need to take smaller steps?
     
  10. May 12, 2015 #9
    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?
     
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook




Similar Discussions: What is the semantic error in my code
  1. Error in code (Replies: 4)

Loading...