python - Lane-emden solutions using RK4 (hard-coded)

Question

I'm working on finding an plotting solutions to the Lane-Emden equation for values n=[0,6], in intervals of 1/2. I'm new to Python, and can't seem to figure out how to use RK4 to make this work. Please help!

Current progress.

TypeError: unsupported operand type(s) for Pow: 'int' and 'list' on line 37 in main.py

The error just appeared after I added in the equations defined as r2, r3, r4 and k2, k3, k4.

import numpy as np
import matplotlib.pyplot as plt

n = [0,1,2,3,4,5,6,7,8,9,10,11,12,13]
theta0 = 1
phi0 = 0
step = 0.01
xi0 = 0
xi_max = 100

theta = theta0
phi = phi0
xi = xi0 + step

Theta = [[],[],[],[],[],[],[],[],[],[],[],[],[],[]]
Phi = [[],[],[],[],[],[],[],[],[],[],[],[],[],[],[]]
Xi = [[],[],[],[],[],[],[],[],[],[],[],[],[],[],[]]

for i in n:
    Theta[i].append(theta)
    Phi[i].append(phi)
    Xi[i].append(xi)
 
def dTheta_dXi(phi,xi): #r1
    return -phi/xi**2
    
def r2(phi,xi):
    return dTheta_dXi(phi+step,xi+step*dTheta_dXi(phi,xi))

def r3(phi,xi):
    return dTheta_dXi(phi+step,xi+step*r2(phi,xi))
    
def r4(phi,xi):
    return dTheta_dXi(phi+step,xi+step*r3(phi,xi))

def dPhi_dXi(theta,xi,n): #k1
    return theta**(n)*xi**2
    
def k2(theta,xi,n):
    return dPhi_dXi(theta+step,xi+step*dPhi_dXi(theta,xi,n),n)

def k3(theta,xi,n):
    return dPhi_dXi(theta+step,xi+step*k2(theta,xi,n),n)
    
def k4(theta,xi,n):
    return dPhi_dXi(theta+step,xi+step*k3(theta,xi,n),n)
    
for i in n: 
    while xi < xi_max:
        if theta < 0:
            break
        dTheta = (step/6)*(dTheta_dXi(phi,xi)+2*r2(phi,xi)+2*r3(phi,xi)+r4(phi,xi))
        dPhi = (step/6)*(dPhi_dXi(theta,xi,i/2.)+2*k2(theta,xi,n)+2*k3(theta,xi,n)+k4(theta,xi,n))
        theta += dTheta
        phi += dPhi
        xi += step
        Theta[i].append(theta)
        Phi[i].append(phi)
        Xi[i].append(xi)
    print i/2., round(xi,4), round(dTheta_dXi(phi,xi),4), round(xi/3./dTheta_dXi(phi,xi),4), round(1./(4*np.pi*(i/2.+1))/dTheta_dXi(phi,xi)**2,4)
    theta = theta0
    phi = phi0
    xi = xi0 + step

Answer 1

Using RK4 for coupled systems

If one understands a first-order system as a vector-valued system working with vector states, then the scalar version of RK4

k1 = f(x,y)
k2 = f(x+0.5*h, y+0.5*h*k1)
k3 = f(x+0.5*h, y+0.5*h*k2)
k4 = f(x+h, y+h*k3)
x,y = x+h, y+h/6*(k1+2*k2+2*k3+k4)

can also directly be used for the vector case. Sometimes it appears educational to implement this component-wise. While in mathematical texts it is preferred to use one-letter variable names, possibly with sub- or superscripts, the variables in the code of programs usually are multi-lettered. So instead of r2 and k2 it would be more descriptive to use k2_Theta and k2_Phi.

Then it becomes rather intuitive that the state used to evaluate the k3 components has the arguments theta+0.5*step*k2_Theta and phi+0.5*step*k2_Phi.

k2_Xi etc. is always 1 for the independent variable, so the value for the 3rd stage is simply xi+0.5*step.

Implementation details RK4

The values of k1 etc. are fixed inside the step and result from the evaluation of the derivatives function. It makes absolutely no sense to declare them as functions themselves. That is, the RK4 step specialized to this situation becomes just

def RK4_update(theta, phi, xi, step, n):
    k1_Theta = dTheta_dXi(phi, xi)
    k1_Phi = dPhi_dXi(theta, xi, n)
    k2_Theta = dTheta_dXi(phi+0.5*step*k1_Phi, xi+0.5*step)
    k2_Phi = dPhi_dXi(theta+0.5*step*k1_Theta, xi+0.5*step, n)
    k3_Theta = dTheta_dXi(phi+0.5*step*k2_Phi, xi+0.5*step)
    k3_Phi = dPhi_dXi(theta+0.5*step*k2_Theta, xi+0.5*step, n)
    k4_Theta = dTheta_dXi(phi+step*k3_Phi, xi+step)
    k4_Phi = dPhi_dXi(theta+step*k3_Theta, xi+step, n)
    dTheta = (step/6)*(k1_Theta+2*k2_Theta+2*k3_Theta+k4_Theta)
    dPhi = (step/6)*(k1_Phi+2*k2_Phi+2*k3_Phi+k4_Phi)
    return dTheta, dPhi

On the singularity of the Lane-Emden equation

For the solution to exist at xi=0 one needs at least that phi ~ xi^k with k>=2. This gives that theta is almost constant, which in turn leads to an integration phi = theta0^n*xi^3/3 which then in the other equation gives theta = theta0 - theta0^n*xi^2/6. This allows to take the first step away from the singularity without using the numerical method.

xi = step 
theta, phi = theta0 - theta0**n*xi**2/6, theta0**n*xi**3/3
Xi[i] = [0, xi]
Theta[i] = [theta0, theta]
Phi[i] = [0, phi]

Then the main loop can be written as

for i in range(N): 
    n = i/2
    xi = step 
    theta, phi = theta0 - theta0**n*xi**2/6, theta0**n*xi**3/3
    Xi[i] = [0, xi]
    Theta[i] = [theta0, theta]
    Phi[i] = [0, phi]
    while xi < xi_max:
        if theta < 0:
            break
        dTheta, dPhi = RK4_update(theta,phi,xi,step,n)
        theta += dTheta
        phi += dPhi
        xi += step
        Theta[i].append(theta)
        Phi[i].append(phi)
        Xi[i].append(xi)

Then plotting with

for i in range(N):
    plt.plot(Xi[i],Theta[i], label=f"n={i/2}")
plt.grid(); plt.legend(); plt.show()

results in

Trick used: To avoid rational powers of negative values, replace theta**n with theta*abs(theta)**(n-1) or similar continuations.

---

Old contents

You should once again explore what update goes where. xi is the independent variable and thus only gets updates 0.5*step and step, the updates of theta use the derivatives dTheta_dXi and similarly phi is updated using the slopes dPhi_dXi

def r2(phi,xi):
    return dTheta_dXi(phi+0.5*step*dPhi_dXi(theta,xi,n),xi+0.5*step)

def k2(theta,xi,n):
    return dPhi_dXi(theta+0.5*step*dTheta_dXi(phi,xi),xi+0.5*step,n)

def r3(phi,xi):
    return dTheta_dXi(phi+0.5*step*k2(theta,xi,n),xi+0.5*step)

etc.

Now one can see that due to the coupled nature of the equation you need both theta and phi as arguments everywhere. Further, even if that works in the end you end up computing many of the values multiple times where assembling everything in one loop only requires one computation.