# Numerical solution to the Lane-Emden equation

I am new to Julia programming, but for fun and practice, I wrote a program that numerically solves the Lane-Emden equation. I even wrote an interactive version that changes the values of n and log_delta_xi with sliders from PlutoUI (not included here).

I will leave the program here and encourage any suggestions as to how it can be improved.

using LaTeXStrings
using Plots
using Markdown
gr()

println("Lane-Emden Equation")
L"""\dfrac{\text{d}}{\text{d}\xi} \
\left( \xi^2 \dfrac{\text{d}\theta}{\text{d}\xi} \right) \
= \
-\xi^2\theta^n \ """

display(md"""
Separation of variables


\dfrac{\text{d}y}{\text{d}\xi}   =  \dfrac{z}{\xi^2}



\dfrac{\text{d}z}{\text{d}\xi}  = -\xi^2y^n

""")

function solveLaneEmden(log_delta_xi=-4, n=3)
delta_xi = 10.0^log_delta_xi

# Inner boundary condition
y0 = 1 - delta_xi^2/6
z0 = -delta_xi^3/3

ys  = [y0]
zs  = [z0]
xis = [delta_xi]
ycs = [y0]
zcs = [z0]

while true
y  =  last(ys)
z  =  last(zs)
xi =  last(xis)
yc =  last(ycs)
zc =  last(zcs)

## Primitive method
yi = y + delta_xi * z/xi^2
zi = z + delta_xi * -xi^2*y^n

## Predictor-corrector technique
xii = xi + delta_xi
yci = yc + 1/2 * delta_xi * (z/xi^2 + zi/xii^2)
zci = zc + 1/2 * delta_xi * (-xi^2*y^n - xi^2*yi^n)

# Outer boundary condition
if (yi < 1e-10 || yci < 1e-10)
break
end

push!(xis, xii)
push!(ys, yi)
push!(zs, zi)
push!(ycs, yci)
push!(zcs, zci)

end

return (xis, ys, ycs)

end

function plotLaneEmden(log_delta_xi=-4, n=3)
xis, ys, ycs = solveLaneEmden(log_delta_xi, n)

xi2 = range(0,sqrt(6),step=1e-3)
# @. will add . to every operator
Plots.plot(xi2, 1 .- xi2.^2/6, label="n = 0")

xi2 = range(0, pi, step=1e-3)
Plots.plot!(xi2, sin.(xi2)./xi2, linecolor = :orange, label="n = 1")

Plots.plot!(xis, ys,  linecolor = :green, label="n = \$n")
Plots.plot!(xis, ycs, linecolor = :black, linestyle = :dash, label="P-C")

Plots.xlabel!("ξ")
Plots.ylabel!("θ")

end

plotLaneEmden()

• Instead of separation of variables, which usually has a different meaning in this context, may I suggest a variation of equivalent first order system? Jan 10 at 9:41
• @unrelenting nosedive I am not familiar with this?
– Dila
Jan 10 at 15:10
• Is this a pedagogical exercise where you try to avoid DifferentialEquations.jl on purpose? Jan 13 at 14:21
• @phipsgabler I suppose, but I do not know how to implement DifferentialEquations.jl with this.
– Dila
Jan 13 at 16:25