Using the Runge Kutta Method to determine mass

[Bishop556]

I am confused on how to use the Runge Kutta method to solve for a relationship between the Chrandrasekhar Mass and radius on the following two equations of state:

dx/dr = (-5/3)*(M/r^2)*[sqrt(1+x^2) /x] where x(r=0) = x_c

dM/dr = +3*(r^2)*(x^3) where M(r=0) = 0

where M is the mass, r is the radius of the white dwarf star, and x is the Fermi momentum. Also, x_c is the central Fermi momentum which we can take as some arbitrary high value.

I understand that for both equations that I have two dependent variables on r, M(r) and x(r), but I do not know how I can relate them. I wish to find radius needed such that x(r) approaches zero as this would indicate the Chrandrasekhar Mass at the same given r. If both equations were not coupled, then this would be a simple problem, but they are. Note that both equations are dimensionless as they have been scaled appropriately.

How would one go about using the Runge Kutta Method to determine the mass/radius relation? I know I need to somehow find a relationship between x and M which is given in the second equation, but that seems difficult to code.

 

[DrClaude]

How would one go about using the Runge Kutta Method to determine the mass/radius relation? I know I need to somehow find a relationship between x and M which is given in the second equation, but that seems difficult to code.

You don't need anything more than a set of coupled differential equations, which you have. The coding will actually be trivial.

In what language are you going to program? C, Fortran, MATLAB, ...?

 

[Bishop556]

I'm programming in Python.

 

[Bishop556]

I just do not have much exposure to the Runge-Kutta method, so I am trying to educate myself online. What would be my thought process to code these equations?

 

[DrClaude]

I just do not have much exposure to the Runge-Kutta method, so I am trying to educate myself online. What would be my thought process to code these equations?

I don't know much python, so some of what I say may need to be adapted for the peculiarities of the language

I would separate the problem into two parts. First, you need to setup a function that, given an array y(n) of dependent variables and the value x of the independent variable, returns an array dydx(n) containing the derivatives of the dependent variables. Second, you setup another function that, given a vector y0(n) of initial conditions for the dependent variables at point x0, returns y(n) at a point x0+h, where the y(n) were calculated by the Runge-Kutta method, using the evaluation of dydx from the first function. Then, you have a main function that will loop N times over the RK integration, to take the system in its initial state to some final x = x0 + N h.

If you wish, you may not even need to implement the second function above, as I am pretty sure you can find ready-to-use RK integrators for python (maybe in sci.py?).

 

[fluidistic]

I just do not have much exposure to the Runge-Kutta method, so I am trying to educate myself online. What would be my thought process to code these equations?

A google search returns the cooked solution. Warning, spoiler: http://rosettacode.org/wiki/Runge-Kutta_method#Python.

 

[Bishop556]

xc = input("Enter the central fermi momentum")
print (xc)

def fa1(x):
    return -1*(5/3.0)*(x[1] / x[0]**2)*math.sqrt(1+x[2]**2)/x[2]                       ''' First ODE'''
    
def fb1(x):
    return 3*(x[0]**2)*(x[2]**3)       '''Second ODE'''

def fc1(x):
    return xc                   '''Fermi Momentum at r=0'''

def fd1(x):
    return 0                    '''Mass at r=0'''

def rKN(x,fx,n,h):
    k1 = []
    k2 = []
    k3 = []
    k4 = []
    xk = []
    while x[2]>0:   '''I want the code to run until the momentum is zero.'''
        for i in range(n):
            k1.append (fx[i](x)*h)
        for i in range(n):
            xk.append(x[i]+k1[i]*0.5)
        for i in range(n):
            k2. append(fx[i](xk)*h)
        for i in range(n):
            xk[i]=x[i]+k2[i]*0.5
        for i in range(n):
            k3.append(fx[i](xk)*h)
        for i in range(n):
            xk[i] = x[i]+k3[i]
        for i in range(n):
            k4.append(fx[i](xk)*h)
        for i in range(n):
            x[i] = x[i] + (k1[i]+2*(k2[i]+k3[i])+k4[i])/6
        return xdef ODES():
    fx = [fa1, fb1, fc1, fd1]
    x = [10**(-23),0,xc,0]
    for i in range (20000):
        x = rKN(x,fx,4,0.0025)
    print (x[0])

So, I took this code from another website, more or less. that was used for a specific type of oscillator. Note that x[0] = radius, x[1] = Mass, x[2] = fermi momentum. My problem is that I wish to allow the code to run until x[2] = 0. However, the code won't print any of the values of the individual elements of my array x after the code runs. I believe I am making a big semantics error as I am a little confused if incorporating the independent variable, that radius known as x[0], in the ODES such that I have was a good idea.

 

[DrClaude]

xc = input("Enter the central fermi momentum")
print (xc)

def fa1(x):
    return -1*(5/3.0)*(x[1] / x[0]**2)*math.sqrt(1+x[2]**2)/x[2]                       ''' First ODE'''
   
def fb1(x):
    return 3*(x[0]**2)*(x[2]**3)       '''Second ODE'''

def fc1(x):
    return xc                   '''Fermi Momentum at r=0'''

def fd1(x):
    return 0                    '''Mass at r=0'''

You have only two dependent variables and two first-order ODEs. You shouldn't define four equations here.

    while x[2]>0:   '''I want the code to run until the momentum is zero.'''

This should not be here. rKN should take one step of length h, that's it.

    fx = [fa1, fb1, fc1, fd1]
    x = [10**(-23),0,xc,0]

Again, fx and x should be of length 4.

    for i in range (20000):
        x = rKN(x,fx,4,0.0025)
    print (x[0])

it is in the loop that you should check if zero momentum is achieved.