1 Introduction

The Lane-Emden equation (LE) is an initial value nonlinear ordinary differential equation describing the equilibrium density distribution of the self-gravitating sphere. LE is often employed to model many problems in physics and astrophysics, like the thermal background of a spherical cloud of gas, stellar structure, and particle currents (^{Chandrasekhar 1967}, ^{Horedt 2004}). Polytropic and isothermal gas spheres have been the subject of many studies given their importance as fundamental models for stellar structure, (^{Chandrasekhar 1967}, Shore 2007, Kippenhahn et al. 2012, ^{Maciel 2016}). The Lane-Emden equation for a polytropic gas sphere can be written as (^{Chandrasekhar 1967})

where θ is called the Emden function and n is the polytropic index. The solution of equation (1) holds under the initial conditions

For equation (1), there are exact solutions only for the polytropic indexes n = 0, 1, and 5.

The LE equation for the isothermal gas sphere is given by (^{Chandrasekhar 1967})

subject to the conditions

The main application of isothermal gas spheres is in the context of molecular clouds, stellar cluster dynamics, and nonlinear analysis in applied mathematics (^{Chandrasekhar 1942}, ^{Kurth 1957}, ^{Davis 1962}, ^{Ito et al. 2018}). In addition to numerical approaches (i.e., neural networks, ^{Ahmad et al. 2017}, ^{Nouh et al. 2021}), genetic algorithms (^{Ahmad et al. 2016}), and the pattern search optimization technique (^{Lewis et al. 2000}), various analytical methods exist to solve Lane-Emden-type equations (i.e., ^{Chowdhury 2009}, ^{Ibrahim and Darus 2008}, Podlubny 1999, ^{Momani and Ibrahim 2008}, ^{Nouh 2004}, ^{Nouh and Saad 2013}).

Monte Carlo (MC) calculations first appeared in astronomy in the late 1960s (^{Auer 1968}, ^{Avery and House 1968}, and ^{Magnan 1968}, ¹⁹⁷⁰). The Markov Chains method (MCMC) is one of the most powerful Monte Carlo algorithms. The MCMC is used extensively for astrometric orbits (^{Tuomi & Kotiranta 2009}, ^{Otor et al. 2016}), visual binary orbits (^{Mendez et al. 2017}), and other problems (see, for instance, ^{Hestroffer 2012}, ^{Mede & Brandt 2014}).

Compared to traditional approaches to solving differential equations, the MC method has been adopted as an independent method as science, technology, and computers have improved over the last century. MC approaches were implemented to solve the initial value problems of ordinary differential equations by Zhong and Tian (2011) and ^{Akhtar et al. (2015)}. ^{Akhtar et al. (2015)} effectively refined the complexities of Zhong and Tian’s algorithm. ^{Uslu and Sari (2020)} presented a Monte Carlo-based stochastic technique to solve various Lotka-Volterra equation systems. ^{El-Essawy et al. (2023}) introduced a novel numerical solution to LE equations (i.e., positive and negative index polytropes, isothermal gas spheres, and the white dwarf equation) using the MC technique.

In the present paper, we solve the LE equations of the polytropic and isothermal gas spheres using two different MC algorithms (MC₁ and MC₂). The results from the two MC algorithms are compared to exact and numerical solutions. As an application, we compute the physical parameters of the polytropic and isothermal gas spheres. The structure of the paper is as follows: we introduce the two MC algorithms in § 2, §3 deals with the results, and § 4 summarises the results.

2 The MC Algorithms

In the following subsections, we shall describe the two algorithms implemented to calculate the Emden function (θ) of the polytropes and its derivative (θ'). The first algorithm (MC1) is proposed by ^{El-Essawy et al. (2023)}, and the second one is based on the algorithm proposed by ^{Akhtar et al. (2015}).

3 Results

To solve the LE differential equations of the polytropic and isothermal gas spheres, we coded MC ₁ and MC ₂ using R language. A total of 1000000 random samples were used. For the numerical solution, we used the R package function rkMethod that implements the Runge-Kutta method (RK). We used an oddly spaced spacing with 0.001 increments. We solved the LE equation (equation 1) for the polytropic indexes n = 0, 1, 1.5, 3, and 5, and the LE equation of the isothermal gas sphere (equation 4); the results are listed in Tables A1-A6 and plotted in Figure 1. The designations of the columns of Tables A1-A6 are as follows: Column 1 is the dimensionless parameter x, Column 2 is the Emden function calculated by the exact solutions (using equation (3) for n = 0,1 , and 5) or from the numerical integration for n = 1.5, 3, and the isothermal gas sphere; Column 3 is the MC solution using algorithm 1 (MC ₁); Column 4 is the MC solution using algorithm 2 (MC ₂); Column 5 is the absolute error such that E1=(exact/numerical-MC1); Column 6 represents the absolute error such that E2=(exact/numerical-MC2). The maximum absolute errors between the two MC algorithms and the exact/numerical solutions are listed in Table 1. As shown in Figure 1 and Table 1, the Emden functions (θ) calculated using the second MC algorithm (MC ₂) have a smaller error than the first MC algorithm (MC ₁) except in the two situations of the polytrope with n = 5 and the isothermal gas sphere.

Fig. 1 The distribution of the Emden function for the polytropes with n = 0, 1, 1.5, 3, 5, and the isothermal gas sphere. The color figure can be viewed online. 

Similar to the scheme we applied in a previous paper (MC₁, ^{El-Essawy et al., 2023}), we determined the absolute errors between the MC₂ solution and several numerical techniques to assess the effect of the MC calculations on the accuracy of the computed Emden function (θ). For the polytropes with n = 1, we used the active-set algorithm-based neural networks (AST-NN, ^{Ahmad et al. 2017}), Chebyshev neural networks (Ch-NN, ^{Mall and Chakraverty 2014}), the Pattern Search Optimization Technique (PS, ^{Lewis et al. 2000}), and genetic algorithms (GA, ^{Ahmad et al. 2016}).

For n = 3, we compared our results with the Homotopy Analysis Method (HAM), Padé approximants (PA) of an order [4; 4], and the numerical solution with the Simpson rule (SIMP); ^{Al-Hayani et al. (2017)}. Tables 2 and 3 compare θ computed by the MC₂ and by different numerical methods. Figures 2 and 3 plot the absolute errors between the MC₂ values and the above mentioned numerical methods. For the polytrope with n = 1, the maximum absolute error between the exact solution and MC₂ is the smallest, while for the polytrope with n = 3, the maximum absolute error between the Simpson rule and MC₂ is the smallest.

Fig. 2 Comparative study based on values of absolute errors for the polytrope with n = 1. The color figure can be viewed online. 

Fig. 3: Comparative study based on values of absolute errors for the polytrope with n = 3. The color figure can be viewed online. 

Polytropic stellar models may be constructed using the polytropic index n = 0-1 for neutron stars, n = 1.5 for white dwarfs, and n = 3 for solar-type stars. The mass contained in a radius R0 is given by

The radius of the polytrope is given by

The pressure is provided by

The central density is computed from the equation

and finally, the ratio of central density to mean density is

Table 4 compares the density ratio ρc/ρm (equation 12) computed from the three methods (ex/num,MC1,MC2). The table shows that the ratio calculated from the two MC algorithms agrees with the exact/numerical values.

We construct mass-radius relation and pressure profiles for the model with a polytropic index n = 3; the zero of the Emden function (x1) and the first derivative of the Emden function (θ'). The results are listed in Tables 5 and 6 for the mass-radius relation and the pressure profile, respectively. Based on the two tables with absolute errors in Columns 5 and 6, one can demonstrate an excellent agreement between the numerical and MC calculations.

4 Conclusion

This paper employed two Monte Carlo (MC) approaches to solve Lane-Emden equations of the polytropic and isothermal gas spheres. We wrote R codes to implement the two MC algorithms, where 10⁶ random samples were used in each integration procedure. The numerical solution of the differential equations is performed using the R package function rkMethod, which is based on the Runge-Kutta technique (RK) for solving ordinary differential equations. An oddly spaced interval with 0.001 increments was chosen.

The MC results were compared to the exact and numerical solutions to determine the accuracy and efficiency of the proposed approach. When the MC and numerical/analytical models are compared, they agree well for the polytropic indexes n = 0, 1, 1.5, 3, while there are significant discrepancies from the models with n = 5 and the isothermal gas sphere.

For the polytropic indexes n = 1 and 3, we compared MC solutions and other numerical approaches like active-set algorithm-based neural networks, Chebyshev neural networks, the Pattern Search Optimization Technique, and genetic algorithms. We may infer that the MC algorithms have the advantage of needing no data training like neural network methods. Based on absolute error, the MC solutions are consistent with solutions from different numerical methods. Moreover, we found that the second algorithm (MC₂) is faster than the first one (MC₁) and more accurate. We computed the mass-radius relation and the pressure profile for the polytrope with n = 3, and we found good agreement between the MC models and the numerical one.