# Solving non-linear singular ODE with SciPy odeint / ODEPACK

I want to solve the Lane-Emden isothermal equation [PDF, eq. 15.2.9]

$$\frac{d^2 \!\psi}{d \xi^2} + \frac{2}{\xi} \frac{d \psi}{d \xi} = e^{-\psi}$$

with the initial conditions

$$\psi(\xi = 0) = 0 \quad \left.\frac{d\psi}{d \xi}\right|_{\xi = 0} = 0$$

using SciPy odeint() but, as it can be seen, the equation is singular at the origin. The documentation states it uses ODEPACK.

I already know the power series of the solution in a neighbourhood of $\xi = 0$ (ref):

$$\psi(\xi) \simeq \frac{\xi^2}{6} - \frac{\xi^4}{120} + \frac{\xi^6}{1890}$$

I tried setting tcrit to np.array([0.0]), but didn't work: I get a warning about invalid values and then my solution is all NaN. Should I integrate starting from 0.01 maybe? Or is there any other solution?

• Note that everything works fine if I start integrating from $t_0 > 0$. I just want to make sure there is no other way. – astrojuanlu Jun 3 '12 at 18:09
• There seem to be some problems with your formulation. I don't know this equation, so I Googled it and only came up with the "Lane"-Emden equation which is somewhat different. Additionally, did you mean your derivatives to be with respect to $\xi$ or $t$? – Bill Barth Jun 3 '12 at 18:38
• @BillBarth, you were right there were errors, I meant $\xi$. I also added a reference to the equation, taken from a Stellar Structure and Evolution course (www2.astro.psu.edu/users/rbc/astro534.html) at the Pennsylvania State University (the Polytropes lesson). It's equation 15.2.9. – astrojuanlu Jun 3 '12 at 19:50

All right, this answer is a shot in the dark, but here goes.

First, transform the second-order ODE into a system of two ODEs. Let

\begin{align} \varphi_{1} &= \psi, \\ \varphi_{2} &= \dot{\psi}, \end{align}

where the dots atop functions correspond to differentiation with respect to the independent variable (in this case, $\xi$).

Then the second-order implicit ODE

\begin{align} \ddot{\psi}(\xi) + 2\xi^{-1}\dot{\psi}(\xi) &= e^{-\psi(\xi)} \\ \psi(0) &= 0 \\ \dot{\psi}(0) &= 0 \end{align}

can be expressed as the first-order explicit ODE

\begin{align} \dot{\varphi}_{1}(\xi) &= \varphi_{2}(\xi) \\ \dot{\varphi}_{2}(\xi) &= -2\xi^{-1}\varphi_{2}(\xi) + e^{-\varphi_{1}(\xi)} \\ \varphi_{1}(0) &= 0 \\ \varphi_{2}(0) &= 0. \end{align}

At first, it would appear that we cannot evaluate the right-hand side of this explicit ODE system at $\xi = 0$, like a numerical integrator requires. If a solution to this system exists, then it must be differentiable. On the assumption that a solution exists, take the limit of the right-hand side as $\xi \rightarrow 0$.

First, we know that

\begin{align} \lim_{\xi \rightarrow 0} \varphi_{2}(\xi) = 0, \end{align}

because we've assumed that a solution exists, so $\varphi_{2}$ is differentiable, which means it must be continuous. The limit of a continuous function at a point is its value at that point, and we know the value of $\varphi_2(0)$ because it is an initial condition.

We also know that

\begin{align} \lim_{\xi \rightarrow 0} e^{-\varphi_{1}(\xi)} = 1 \end{align}

for similar reasons; we've assumed that $\varphi_{1}$ is differentiable, so it is continuous, and $\varphi_{1}(0) = 0$ because it is an initial condition.

Finally,

\begin{align} \lim_{\xi \rightarrow 0} \frac{-2\varphi_{2}(\xi)}{\xi} = \lim_{\xi \rightarrow 0} -2\dot{\varphi}_{2}(\xi), \end{align}

by using L'Hôpital's rule on the indeterminate form $0/0$.

To proceed further, we have to make another assumption: $\dot{\varphi}_{2}$ is continuous at $\xi = 0$. Then it follows that

\begin{align} \lim_{\xi \rightarrow 0} -2\dot{\varphi}_{2}(\xi) = -2\dot{\varphi_{2}}(0). \end{align}

Revisiting the first-order ODE, and evaluating the right-hand side at $\xi = 0$, we can see that we have:

\begin{align} \dot{\varphi}_{1}(0) = 0 \\ \dot{\varphi}_{2}(0) = -2\dot{\varphi}_{2}(0) + 1, \end{align}

from which it follows that $\dot{\varphi}_{2}(0) = 1/3$.

Using this analysis, you could plug in an if statement that returns these values of the right-hand side function at $\xi = 0$, which should get you past the singularity. That said, this analysis requires a couple assumptions about continuity that may or may not hold, so take the resulting solution with a grain of salt.

• Geoff, maybe you meant $\varphi_1 = \psi$ and $\varphi_2 = \dot{\psi}$? I added to my answer that I already know the power series of the solution at the origin (maybe I should have done that before). Anyway, you showed that $\ddot{\psi}(0) = 1 / 3$, which is the case. I'll try the if statement thing. – astrojuanlu Jun 4 '12 at 5:49
• @Juanlu001: Good call; I corrected the mistake. – Geoff Oxberry Jun 4 '12 at 5:54
• You were right! It was as simple as a simple if clause. Thank you! – astrojuanlu Jun 4 '12 at 10:24

If you want more options for your ODE solver, have a look at the Assimulo package that implements bindings to the CVODE package (and RADAU and some simple integrators for that matter).

• Very valuable package, didn't know it! Thank you very much. – astrojuanlu Jun 4 '12 at 14:35
• @GertVdE: Would any of these integrators be able to handle the singularity without resorting to using an if statement to fill in the proper limit of the right-hand side as $\xi \rightarrow 0$? – Geoff Oxberry Jun 4 '12 at 15:48
• @GeoffOxberry: I believe the IDA solver of the Sundials suite (which Assimulo wraps for Python) allows to search for a consistent initial value starting from a user's guess. This would allow Juanlu001 to start from his series expansion as initial guess and let IDA solve for the correct (numerically, that is) IV. – GertVdE Jun 4 '12 at 20:22