pisces.math_utils.integration.compute_lame_emden_solution#

pisces.math_utils.integration.compute_lame_emden_solution(xi_max: float, n: float = 1, xi_min: float = 1e-05, end_at_first_zero: bool = False, **kwargs)[source]#

Compute a numerical solution to the Lane-Emden equation for a given polytropic index \(n\).

The Lane-Emden equation models the structure of polytropic spheres such as stars, and is given by:

\[\frac{1}{\xi^2} \frac{d}{d\xi} \left( \xi^2 \frac{d\theta}{d\xi} \right) + \theta^n = 0,\]

with boundary conditions:

\[\theta(0) = 1, \quad \theta'(0) = 0.\]

Because the equation is singular at \(\xi = 0\), we begin integration at a small nonzero value \(\xi_{\rm min}\) using a Taylor expansion to approximate the initial conditions.

Parameters:

  • xi_max (float) –

    Maximum value of the dimensionless radial coordinate \(\xi\) to integrate to. Depending on the index and the value of the end_at_first_zero flag, this may not be reached if the solution crosses zero. If the solution does not trigger a termination event, the integration will continue to this value.

    In general, it is best to set this well beyond the expected extent of the solution, such as 10 or 20.

  • n (float, optional) – Polytropic index. Default is 1.

  • xi_min (float, optional) – Starting value of \(\xi\) to avoid the singularity at zero. Default is 1e-5.

  • end_at_first_zero (bool, optional) – If True, halt integration when \(\theta\) first crosses below zero. This is overruled if n is a non-integer, in which case the integration will always truncate to avoid non-real solutions.

  • **kwargs – Additional keyword arguments passed to scipy.integrate.solve_ivp(). This can be used to control the integration method, tolerances, and other parameters.

Returns:

solution – Full solution object from scipy.integrate.solve_ivp(). Includes:

  • solution.t: array of \(\xi\) values

  • solution.y: 2D array of corresponding \(\theta(\xi)\) and auxiliary variables

  • solution.status: solver status code (0=success, 1=event-triggered)

  • solution.t_events: event times, if any

  • solution.message: diagnostic message

Return type:

scipy.integrate.OdeResult

Notes

For small \(\xi\), we approximate:

\[\theta(\xi) \approx 1 - \frac{1}{6} \xi^2, \quad \theta'(\xi) \approx -\frac{1}{3} \xi,\]

which are valid for small deviations from the origin. These values are used to initialize the solver at \(\xi = \xi_{\rm min}\).

Additionally, when end_at_first_zero is True, or when n is not an integer, we define an event function that triggers when \(\theta\) crosses below a small threshold (e.g., 0.01). This allows the integration to stop gracefully when the solution is no longer valid, avoiding numerical instability.

Examples

from pisces.math_utils.integration import (
    compute_lame_emden_solution,
)

xi, theta = compute_lame_emden_solution(
    xi_max=10.0,
    n=1.5,
    end_at_first_zero=True,
    max_step=0.1,
)