Singular Lane–Emden–Fowler Equations

• Singular Lane–Emden–Fowler equations are nonlinear differential equations characterized by singular coefficients or boundary conditions that induce unbounded solution behavior.
• The topic integrates asymptotic analysis, boundary regularity techniques, and a variety of numerical methodologies, including collocation, spectral, and neural network-based approaches.
• Applications span astrophysics and thermal explosion models, demonstrating the equations’ significance in modeling stellar structures and complex boundary phenomena.

Searching arXiv for papers on singular Lane–Emden–Fowler equations and related numerical/analytic methods. The singular Lane–Emden–Fowler equation denotes a family of nonlinear differential equations in which singular behavior is generated either by a coefficient that blows up at a distinguished point, typically $x=0$, or by boundary data that force the unknown to vanish where the nonlinearity contains a negative power of the solution. In the material surveyed here, two principal settings appear. The first is the semilinear elliptic PDE

$-\Delta u=f(X)\,u^{-\gamma},$

posed in a bounded Lipschitz domain with Dirichlet data, where the singularity is induced by the boundary condition $u=0$ on part of $\partial\Omega$ (Guo et al., 20 Mar 2025). The second is the classical ODE family

$y''+\frac{\alpha}{x}y'+Q(x)\,y=f(x,y),$

or closely related Emden–Fowler forms, where the term $\alpha/x$ is singular at the origin (Patel et al., 2024). Across both settings, the subject combines local singular analysis, existence and uniqueness theory, asymptotics, boundary regularity, and a wide range of numerical methodologies.

1. Problem classes and sources of singularity

In the elliptic PDE setting studied by Guo, Li, and Zhang, one works in a bounded Lipschitz domain $\Omega\subset\mathbb R^n$, $n\ge2$, with a positive Hölder-continuous weight $f$ satisfying $0<\lambda\le f\le\Lambda$, and seeks a classical solution

$-\Delta u=f(X)\,u^{-\gamma},$0

together with arbitrary continuous data on $-\Delta u=f(X)\,u^{-\gamma},$1. The singularity arises precisely because $-\Delta u=f(X)\,u^{-\gamma},$2 on $-\Delta u=f(X)\,u^{-\gamma},$3, making the right-hand side unbounded (Guo et al., 20 Mar 2025).

In the ODE literature, the singular Lane–Emden–Fowler equation typically refers to a second-order equation with a radial term $-\Delta u=f(X)\,u^{-\gamma},$4. A representative form is

$-\Delta u=f(X)\,u^{-\gamma},$5

with singular initial conditions at $-\Delta u=f(X)\,u^{-\gamma},$6 (Patel et al., 2024). Closely related formulations include

$-\Delta u=f(X)\,u^{-\gamma},$7

for singular initial-value or boundary-value problems (Seiler et al., 2023), and

$-\Delta u=f(X)\,u^{-\gamma},$8

for Emden–Fowler type IVPs (Bildik et al., 2016). In each case, the term $-\Delta u=f(X)\,u^{-\gamma},$9 or $u=0$0 becomes unbounded as $u=0$1, so regularity conditions such as $u=0$2 or $u=0$3 are imposed to ensure bounded behavior.

A separate class of singularity appears in perturbed boundary-value problems with singular endpoints at both $u=0$4 and $u=0$5. One example is

$u=0$6

where $u=0$7 and $u=0$8 contain poles at $u=0$9 and $\partial\Omega$0 vanishes at $\partial\Omega$1 (Kycia, 2018). Another variant places the singular coefficient at the upper endpoint $\partial\Omega$2, in the time-singular form

$\partial\Omega$3

so that $\partial\Omega$4 becomes singular as $\partial\Omega$5 (Hazaimah, 2024).

These formulations are related by the common feature that a formally lower-order term, or the nonlinearity itself, becomes dominant near a singular point or boundary. This suggests that “singular Lane–Emden–Fowler equation” is best viewed as a structural class rather than a single canonical equation.

2. Boundary regularity and cone-frequency classification for the PDE

For the semilinear elliptic problem in a Lipschitz domain, the principal geometric device is the limiting cone at a boundary point. If $\partial\Omega$6 is a spherical open set with Lipschitz boundary, its cone is

$\partial\Omega$7

Let $\partial\Omega$8 be the first Dirichlet eigenvalue on $\partial\Omega$9, and define the minimal frequency $y''+\frac{\alpha}{x}y'+Q(x)\,y=f(x,y),$0 by

$y''+\frac{\alpha}{x}y'+Q(x)\,y=f(x,y),$1

The pair $y''+\frac{\alpha}{x}y'+Q(x)\,y=f(x,y),$2 is called sub-critical if $y''+\frac{\alpha}{x}y'+Q(x)\,y=f(x,y),$3, critical if $y''+\frac{\alpha}{x}y'+Q(x)\,y=f(x,y),$4, and super-critical if $y''+\frac{\alpha}{x}y'+Q(x)\,y=f(x,y),$5 (Guo et al., 20 Mar 2025).

This classification governs the boundary growth rate. If $y''+\frac{\alpha}{x}y'+Q(x)\,y=f(x,y),$6 solves $y''+\frac{\alpha}{x}y'+Q(x)\,y=f(x,y),$7 in a Lipschitz cylinder $y''+\frac{\alpha}{x}y'+Q(x)\,y=f(x,y),$8, then in $y''+\frac{\alpha}{x}y'+Q(x)\,y=f(x,y),$9 there is a constant $\alpha/x$0 such that:

Regime	Upper growth rate
Sub-critical	$\alpha/x$1
Critical	$\alpha/x$2
Super-critical	$\alpha/x$3

Moreover, if $\alpha/x$4 near $\alpha/x$5, these rates are optimal (Guo et al., 20 Mar 2025).

Along the inward normal ray $\alpha/x$6, the corresponding lower bounds are also regime-dependent: $\alpha/x$7

$\alpha/x$8

and

$\alpha/x$9

(Guo et al., 20 Mar 2025). The critical logarithmic correction is a distinctive feature of the borderline frequency balance.

The significance of this framework is that boundary regularity is not described by a single universal exponent. Instead, the interaction between the singular power $\Omega\subset\mathbb R^n$0 and the local cone geometry determines the asymptotic profile. A plausible implication is that non-smooth boundary geometry is not merely a technical complication but part of the principal asymptotic mechanism.

3. Well-posedness, comparison principles, and analytic continuation

For the Lipschitz-domain PDE, the paper establishes well-posedness together with boundary growth and boundary Harnack estimates (Guo et al., 20 Mar 2025). Several auxiliary results support this theory. A maximum-principle lemma states that if $\Omega\subset\mathbb R^n$1 and $\Omega\subset\mathbb R^n$2 are respectively super- and sub-solutions of $\Omega\subset\mathbb R^n$3 and $\Omega\subset\mathbb R^n$4 on $\Omega\subset\mathbb R^n$5, then $\Omega\subset\mathbb R^n$6 in $\Omega\subset\mathbb R^n$7. A non-degeneracy lemma states that any super-solution $\Omega\subset\mathbb R^n$8 in $\Omega\subset\mathbb R^n$9 satisfies $n\ge2$0, preventing excessively flat behavior near the boundary (Guo et al., 20 Mar 2025).

For singular ODEs, existence and uniqueness are often obtained by reformulating the equation as an equivalent regular problem. Seiler and Seiß rewrite

$n\ge2$1

as a three-dimensional autonomous system in $n\ge2$2: $n\ge2$3 The point $n\ge2$4 is a stationary point, and the unique solution of the singular initial-value problem is exactly the integral curve lying on the one-dimensional unstable manifold of $n\ge2$5 (Seiler et al., 2023). Under the hypotheses $n\ge2$6 and $n\ge2$7 in the generalized formulation $n\ge2$8, they state that there is a unique smooth solution through the singular point whose prolonged graph is that unstable manifold (Seiler et al., 2023).

In the perturbed endpoint-singular problem on $n\ge2$9, local analytic solutions are constructed separately at $f$0 and $f$1 using convergent Frobenius-type series, and a Lyapunov function

$f$2

is used to show $f$3 on $f$4, ruling out movable blow-up in $f$5 (Kycia, 2018). Under the stated matching and Lyapunov-type conditions, the paper concludes that there exists at least one nontrivial solution analytic at both singular endpoints.

These approaches illustrate three distinct well-posedness paradigms: comparison and barrier methods for singular semilinear PDEs, dynamical-systems desingularization for ODEs at $f$6, and endpoint power-series matching plus Lyapunov control for problems singular at both ends.

4. Boundary Harnack theory and the failure of the classical linear picture

A central contribution in the Lipschitz-domain PDE theory is the boundary Harnack principle for the singular equation. If $f$7 both solve

$f$8

in $f$9 with zero boundary data on $0<\lambda\le f\le\Lambda$0, then in $0<\lambda\le f\le\Lambda$1,

$0<\lambda\le f\le\Lambda$2

with $0<\lambda\le f\le\Lambda$3 (Guo et al., 20 Mar 2025). The paper identifies this as a generalization of the classical Kemper estimate for harmonic functions and states that, to the authors’ knowledge, it is the first Kemper-type estimate for singular semilinear equations (Guo et al., 20 Mar 2025).

The semilinear setting differs sharply from the linear one. In particular, the boundedness of $0<\lambda\le f\le\Lambda$4 does not imply continuity up to $0<\lambda\le f\le\Lambda$5, and this failure is formalized in Theorem 1.6 of the paper (Guo et al., 20 Mar 2025). Under sub-critical or critical hypotheses, one can show $0<\lambda\le f\le\Lambda$6 at the boundary, while under additional smoothness or convexity of $0<\lambda\le f\le\Lambda$7 one recovers full continuity (Guo et al., 20 Mar 2025).

A key technical input is the lack of a suitable upper barrier in the singular semilinear case. To overcome this, the paper introduces iterative upper-barrier constructions and a Campanato-type iteration. One lemma states that if $0<\lambda\le f\le\Lambda$8 in $0<\lambda\le f\le\Lambda$9, with $-\Delta u=f(X)\,u^{-\gamma},$00 and $-\Delta u=f(X)\,u^{-\gamma},$01 for $-\Delta u=f(X)\,u^{-\gamma},$02, then $-\Delta u=f(X)\,u^{-\gamma},$03 in $-\Delta u=f(X)\,u^{-\gamma},$04; this underlies the iteration proving $-\Delta u=f(X)\,u^{-\gamma},$05 in the sub-critical and critical regimes (Guo et al., 20 Mar 2025).

The resulting theory corrects a common overextension of linear intuition: for singular semilinear equations, the ratio of two positive solutions need not inherit the boundary continuity properties familiar from harmonic analysis. This distinction is one of the main conceptual differences between the singular Lane–Emden–Fowler PDE and linear elliptic theory.

5. Asymptotic constructions, barriers, and auxiliary functions

The proof of sharp boundary growth in cones is built around a rescaled difference

$-\Delta u=f(X)\,u^{-\gamma},$06

where $-\Delta u=f(X)\,u^{-\gamma},$07 solves the constant-weight problem $-\Delta u=f(X)\,u^{-\gamma},$08. This function is nonnegative, subharmonic in the cone, and vanishes on its boundary (Guo et al., 20 Mar 2025). By comparison with the homogeneous harmonic profile

$-\Delta u=f(X)\,u^{-\gamma},$09

one obtains

$-\Delta u=f(X)\,u^{-\gamma},$10

(Guo et al., 20 Mar 2025).

The barrier iteration is then encoded in the dyadic quantities

$-\Delta u=f(X)\,u^{-\gamma},$11

for which

$-\Delta u=f(X)\,u^{-\gamma},$12

A geometric or log-geometric summation yields the three growth regimes (Guo et al., 20 Mar 2025). The auxiliary function $-\Delta u=f(X)\,u^{-\gamma},$13 thus measures how far $-\Delta u=f(X)\,u^{-\gamma},$14 is from exact homogeneity of degree $-\Delta u=f(X)\,u^{-\gamma},$15, and its decay controls the dyadic increments of the barrier sequence.

In endpoint-singular ODE problems, analogous asymptotic information is encoded in local power-series expansions. At $-\Delta u=f(X)\,u^{-\gamma},$16, one sets

$-\Delta u=f(X)\,u^{-\gamma},$17

while at $-\Delta u=f(X)\,u^{-\gamma},$18, with $-\Delta u=f(X)\,u^{-\gamma},$19,

$-\Delta u=f(X)\,u^{-\gamma},$20

and derives recurrences whose convergence is guaranteed by a normal-form reduction and an analytic existence theorem (Kycia, 2018). In the time-singular model, the classical Lane–Emden cases $-\Delta u=f(X)\,u^{-\gamma},$21 are given explicitly as

$-\Delta u=f(X)\,u^{-\gamma},$22

and for other integer indices a Frobenius-type power series about $-\Delta u=f(X)\,u^{-\gamma},$23 is used (Hazaimah, 2024).

These constructions indicate that singular Lane–Emden–Fowler problems are often controlled not by a direct closed-form solution but by carefully chosen comparison objects: harmonic profiles in cones, dyadic barrier sequences, unstable manifolds, or endpoint Frobenius expansions.

6. Numerical methods and computational formulations

The numerical literature represented here is diverse, but it is organized around a shared objective: regularize or bypass the singular structure while preserving the nonlinear dynamics.

One family uses basis expansions and collocation. The orthonormal polynomial wavelet method considers

$-\Delta u=f(X)\,u^{-\gamma},$24

with $-\Delta u=f(X)\,u^{-\gamma},$25 and $-\Delta u=f(X)\,u^{-\gamma},$26, expands the unknown in wavelet bases generated by orthonormal polynomials, and enforces the ODE at midpoint collocation nodes (Verma et al., 2019). Nonlinearity is handled either by Newton’s quasilinearization or by Newton–Raphson on the discrete system (Verma et al., 2019). The paper states that, as the resolution is increased, the computed solutions converge to exact or known solutions, and reports errors of order $-\Delta u=f(X)\,u^{-\gamma},$27–$-\Delta u=f(X)\,u^{-\gamma},$28 in a Chandrasekhar stellar-structure test and maximum errors $-\Delta u=f(X)\,u^{-\gamma},$29 in a thermal explosion problem (Verma et al., 2019).

A second family uses rational spectral or reproducing-kernel representations. The Rational Chebyshev of Second Kind collocation method solves

$-\Delta u=f(X)\,u^{-\gamma},$30

on $-\Delta u=f(X)\,u^{-\gamma},$31 by mapping the semi-infinite interval to $-\Delta u=f(X)\,u^{-\gamma},$32 and using the ansatz

$-\Delta u=f(X)\,u^{-\gamma},$33

so that the center conditions are built in (Parand et al., 2015). The reproducing-kernel Hilbert-space method rewrites the singular IVP as $-\Delta u=f(X)\,u^{-\gamma},$34, constructs a reproducing kernel $-\Delta u=f(X)\,u^{-\gamma},$35, and represents the solution by a convergent series in an orthonormal system derived from kernel functions (Abdullah et al., 2017). The latter paper states that the $-\Delta u=f(X)\,u^{-\gamma},$36-term approximation converges uniformly on $-\Delta u=f(X)\,u^{-\gamma},$37 and reports maximum absolute errors below $-\Delta u=f(X)\,u^{-\gamma},$38 in its examples (Abdullah et al., 2017).

A third family is based on dynamical systems and shooting. For

$-\Delta u=f(X)\,u^{-\gamma},$39

Seiler and Seiß compute the unstable manifold of the associated autonomous vector field and combine this with a shooting procedure for boundary-value problems (Seiler et al., 2023). The paper states that no special series expansions at $-\Delta u=f(X)\,u^{-\gamma},$40 are needed, and reports that default Maple RKF4/5 tolerances $-\Delta u=f(X)\,u^{-\gamma},$41, $-\Delta u=f(X)\,u^{-\gamma},$42 gave 6–7 digits of accuracy on test problems (Seiler et al., 2023).

A fourth family uses neural networks. The shifted Legendre neural network method builds the trial solution

$-\Delta u=f(X)\,u^{-\gamma},$43

so that $-\Delta u=f(X)\,u^{-\gamma},$44 and $-\Delta u=f(X)\,u^{-\gamma},$45 are satisfied automatically and the singular $-\Delta u=f(X)\,u^{-\gamma},$46 term does not cause numerical blow-up (Patel et al., 2024). For the standard $-\Delta u=f(X)\,u^{-\gamma},$47 Lane–Emden case, the paper reports max absolute error $-\Delta u=f(X)\,u^{-\gamma},$48, while for an inhomogeneous Lane–Emden–Fowler problem it reports max error $-\Delta u=f(X)\,u^{-\gamma},$49 (Patel et al., 2024). Physics-informed neural networks were also benchmarked for second-, third-, and fourth-order singular ODEs, using both soft constraints in the loss and hard-constraint trial solutions such as $-\Delta u=f(X)\,u^{-\gamma},$50 or $-\Delta u=f(X)\,u^{-\gamma},$51 (Baty, 2023). The hard-constraint formulation guarantees exact boundary satisfaction, and the paper reports errors ranging from $-\Delta u=f(X)\,u^{-\gamma},$52 to $-\Delta u=f(X)\,u^{-\gamma},$53 across the benchmark suite (Baty, 2023).

A fifth line of work emphasizes high-order iterative discretization. For the singular boundary-value problem

$-\Delta u=f(X)\,u^{-\gamma},$54

Dang et al. derive a continuous fixed-point iteration using the Green’s function and then discretize it by a corrected trapezoidal, Euler–Maclaurin-based rule (A et al., 27 Aug 2025). They present three discrete iterative schemes for three specific cases and state that the resulting methods achieve eighth-order accuracy and convergence (A et al., 27 Aug 2025). In one benchmark, the reported errors decrease from $-\Delta u=f(X)\,u^{-\gamma},$55 at $-\Delta u=f(X)\,u^{-\gamma},$56 to $-\Delta u=f(X)\,u^{-\gamma},$57 at $-\Delta u=f(X)\,u^{-\gamma},$58, with observed orders $-\Delta u=f(X)\,u^{-\gamma},$59, $-\Delta u=f(X)\,u^{-\gamma},$60, and $-\Delta u=f(X)\,u^{-\gamma},$61 on successive refinements (A et al., 27 Aug 2025).

Taken together, these methods show no single dominant numerical paradigm. Collocation, spectral rational bases, reproducing kernels, unstable-manifold shooting, perturbation iteration, and neural-network residual minimization each target a different structural aspect of the singularity.

7. Applications, interpretations, and recurrent themes

The classical Lane–Emden equation appears in astrophysics through stellar structure and polytropic gas spheres. The time-singular study states that the radial Lane–Emden equation models phenomena in astrophysics such as stellar structure and is governed by polytropics with applications in isothermal gas spheres (Hazaimah, 2024). In the same source, the dimensionless variable $-\Delta u=f(X)\,u^{-\gamma},$62 is interpreted as proportional to the density profile of a self-gravitating polytropic gas sphere, with the polytropic index

$-\Delta u=f(X)\,u^{-\gamma},$63

determined by the adiabatic exponent (Hazaimah, 2024).

Other model problems represented in the numerical literature include Chandrasekhar’s stellar structure equation, thermal explosion models, shallow membrane cap equations, human-head heat conduction, Michaelis–Menten kinetics, oxygen-uptake problems, and the Thomas–Fermi equation after the Majorana transform (Verma et al., 2019, Seiler et al., 2023). These examples show that the singular Lane–Emden–Fowler framework extends well beyond polytropes.

A recurrent theme across the analytic and computational work is that the singularity is often geometrically or structurally removable, but not negligible. In ODEs it may be “blown up” into a regular equilibrium of an autonomous system (Seiler et al., 2023); in PDEs it may be encoded by cone geometry and frequency (Guo et al., 20 Mar 2025); in numerical schemes it is absorbed into the ansatz, Green’s function, or choice of collocation points (Verma et al., 2019, Patel et al., 2024). This suggests that successful analysis of singular Lane–Emden–Fowler equations depends less on suppressing the singularity than on reformulating it in a structure-compatible way.

A further recurrent theme is that singular problems do not automatically inherit the regularity properties of their nonsingular analogues. The failure of continuity of $-\Delta u=f(X)\,u^{-\gamma},$64 in the semilinear boundary Harnack problem (Guo et al., 20 Mar 2025) is one example; the need for trial solutions that encode exact endpoint behavior in neural and collocation methods is another (Patel et al., 2024, Baty, 2023). In this sense, the singular Lane–Emden–Fowler equation functions as a meeting point of asymptotic analysis, nonlinear elliptic theory, singular dynamical systems, and high-accuracy numerical computation.