Solving The Lane-Emden Equation


Introduction

Most of theoretical astrophysics involves the numerical solution of various differential equations, both ordinary DEs (ODEs) and partial DEs (PDEs). Our aim here is to solve the Lane-Emden equation. This is a second order ODE.

The Lane Emden Equation

(NOTE: Sometimes the equations does not render upon page load. You might want to refresh the page for the equations to display properly.)

The dimensionless Lane-Emden equation is:

$$\dfrac{1}{{\xi}^2}\dfrac{d}{d\xi}{\xi}^2\dfrac{d\theta}{d\xi} = - {\theta}^n $$

Here \( \xi \) is the scaled radius, \( \theta \) is the scaled temperature, and \(n\) is the index.

Where the initial conditions at \(\xi = 0\) are, \( \theta = 1 \) and \( \dfrac{d\theta}{d\xi} = 0 \).

Now we will be looking at how to approach a solution to such a second order ODE.

We first need to write the 2nd order Lane-Emden equation as two first order equations. In the code we will use \(t\) for \(\xi\), \(y\) for \(\theta\) and \(z\) for \(\dfrac{d\theta}{d\xi}\). The Lane-Emden equation is then:

$$\dfrac{d}{dt}({t}^2\dfrac{dy}{dt}) = - {t}^2{y}^n $$

We write this as two first order equations by setting \(z = \dfrac{dy}{dt}\). The two equations thus become:

\( \dfrac{dy}{dt} = z \)
\( \dfrac{dz}{dt} = -\dfrac{2}{t}z - {y}^n \)

and the initial conditions become:

\( y(0) = 1 \) and \( z(0) = 0 \)

Solution using Euler Method

I will update the contents later.

Solution using Runge-Kutta Method

The Runge-Kutta Method

The Euler method is a particular a class of techniques known as "Runge-Kutta" method, which have the general form:
$$ y_{i+1} = y_{i} + h\phi $$
where, \(\phi\) is some approximation to the slope. In Euler Method \(\phi = k_{1} = f(x_{i},y_{i})\).

The most used Runge-Kutta method is the 4th order Runge-Kutta Method (hereafter RK4), which has a global error of order \(h^{4}\). In this method we take:

\( k_{1} = f(x_{i},y_{i}) \)
\( k_{2} = f(x_{i}+\dfrac{h}{2},y_{i}+\dfrac{hk_{1}}{2}) \)
\( k_{3} = f(x_{i}+\dfrac{h}{2},y_{i}+\dfrac{hk_{1}}{2}) \)
\( k_{4} = f(x_{i}+h,y_{i}+hk_{3}) \)

and

\( y_{i+1} = y_{i} + h\phi \)

where

\( \phi = \dfrac{1}{6}(k_{1}+2k_{2}+2k_{3}+k_{4}) \)

Using RK4 for 2nd order ODE

Any second order differential equation can be written as two coupled first order equations,

\( \dfrac{dy}{dt} = z \)   and   \( \dfrac{dz}{dt} = F(t, y, z) \)

with initial conditions $ y(0) = y_{0} $   and   $ z(0) = z_{0} $ and step size $ h $;

These coupled equations can be solved numerically using a fourth order Runge-Kutta method as follows:

$ k_{y1} = z(0) = z_{0} $     and     $ k_{z1} = F(t, y_{0}, z_{0}) $

$ k_{y2} = z_{0} + \dfrac{hk_{z1}}{2} $     and     $ k_{z2} = F(t+\dfrac{h}{2}, y_{0}+\dfrac{hk_{y1}}{2}, z_{0}+\dfrac{hk_{z1}}{2}) $

$ k_{y3} = z_{0} + \dfrac{hk_{z2}}{2} $     and     $ k_{z3} = F(t+\dfrac{h}{2}, y_{0}+\dfrac{hk_{y2}}{2}, z_{0}+\dfrac{hk_{z2}}{2}) $

$ k_{y4} = z_{0} + hk_{z3} $     and     $ k_{z4} = hF(t+h, y_{0}+hk_{y3}, z_{0}+hk_{z3}) $

$ dy = \dfrac{h}{6}(k_{y1}+2k_{y2}+2k_{y3}+k_{y4}) $   and   $ dz = \dfrac{h}{6}(k_{z1}+2k_{z2}+2k_{z3}+k_{z4}) $

So now the RK4 approximation of next values are:

$ y_{1} = y_{0}+dy $    and    $ z_{1} = z_{0}+dz $

and so on.

Pseudo code

Declare function F(t,y,z)
Declare Tinitial
Declare Tfinal
Declare N = Number of steps
Declare h = (Tinitial - Tfinal)/N
Declare Yinitial
Declare Zinitial
Y = Yinitial
Z = Zinitial
Print Yinitial, Zinitial
LOOP for T from Tinitial to Tfinal with increment of h
        ky1 = Z;
        kz1 = F(t, Y, Z);
        
        ky2 = (Z + h*kz1/2);
        kz2 = F(t + h/2, Y + h*ky1/2, Z + h*kz1/2);
        
        ky3 = (Z + h*kz2/2);
        kz3 = F(t + h/2, Y + h*ky2/2, Z + h*kz2/2);
        
        ky4 = (Z + h*kz3);
        kz4 = F(t + h, Y + h*ky3, Z + h*kz3);
        
        dy = h/6*(ky1 + ky2 + ky3 + ky4);
        dz = h/6*(kz1 + kz2 + kz3 + kz4);
        
        Y = Y + dy;
        Z = Z + dz;
        
        Print Y, Z
END LOOP

Define function F(t,y,z) {
    equation for F(t,y,z) goes here
}
Show Comments