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?

  • $\begingroup$ Note that everything works fine if I start integrating from $t_0 > 0$. I just want to make sure there is no other way. $\endgroup$ – astrojuanlu Jun 3 '12 at 18:09
  • 1
    $\begingroup$ 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$? $\endgroup$ – Bill Barth Jun 3 '12 at 18:38
  • $\begingroup$ @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. $\endgroup$ – 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.


\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.

| cite | improve this answer | |
  • $\begingroup$ 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. $\endgroup$ – astrojuanlu Jun 4 '12 at 5:49
  • $\begingroup$ @Juanlu001: Good call; I corrected the mistake. $\endgroup$ – Geoff Oxberry Jun 4 '12 at 5:54
  • $\begingroup$ You were right! It was as simple as a simple if clause. Thank you! $\endgroup$ – 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).

| cite | improve this answer | |
  • $\begingroup$ Very valuable package, didn't know it! Thank you very much. $\endgroup$ – astrojuanlu Jun 4 '12 at 14:35
  • $\begingroup$ @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$? $\endgroup$ – Geoff Oxberry Jun 4 '12 at 15:48
  • 1
    $\begingroup$ @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. $\endgroup$ – GertVdE Jun 4 '12 at 20:22

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.