A quintic B-spline technique for a system of Lane-Emden equations arising in theoretical physical applications

1.   Introduction

The Lane-Emden equation represents a dimensionless form of Poisson's equation that arises in astrophysics for the spherically symmetric, polytrophic fluids, and the gravitational potential of Newtonian self-gravitating ^[1,2,3]. Modeling diverse phenomena in astrophysics, physical, and mathematical physics, such as the stellar structure theory, isothermal gas spheres, the thermal dynamic of a spherical gas cloud, the thermionic current theory, chemical reactions, population evolution, and pattern formation, results in scalar and systems of Lane-Emden equations, see ^[4,5] and the references therein. While there has been little research on Lane-Emden equation systems, recent attention to studying this type of systems has increased considerably ^[5].

In this study, we consider the following Lane-Emden system of the form

subject to

where ${\delta }_{1}, {\delta }_{2}, {\varepsilon }_{1}, \;\mathrm{a}\mathrm{n}\mathrm{d}\;{\varepsilon }_{2}$ are real constants, and ${\hslash }_{i}$ and ${\aleph }_{i}\left(\tau \right), i = 1, 2$ are given continuous functions.

Numerous approaches have been established for solving scalar and systems of Lane-Emden equations, including the Haar wavelet collocation method ^[6], Laplace transform and residual error function ^[7], Bernoulli wavelets functional matrix technique ^[5], B-spline methods ^[8,9,10], Adomian decomposition method ^[9], Chebyshev operational matrix method ^[4], variational iteration method ^[11,12], discontinuous finite element method ^[13], Bernstein collocation method ^[14], Bessel-collocation procedure ^[15], and Legendre Polynomials ^[16,17,18].

The literature survey reveals that collocation methods are an important tool in obtaining approximate solutions for different types of differential equations, including different classes of initial and boundary value problems, Singular differential equations, partial and fractional partial differential equations, system of partial differential equations, fractional Volterra integro-differential equations, and Abel's integral equations, ^{[19,20,21,22,23,24,25,26,27,28,29,30,31,32]}, among others. One well-known established method among collocation methods is the so-called B-spline method, where the letter "B" represents "basis". This method was originally introduced by Schoenberg in 1946. The primary motivation for introducing the B-splines is the creation of a stable interpolating function across finite number of points, which maintain the smoothness and the shape of the data ^[33,34]. Recently, B-spline methods have been demonstrated to be useful in approximation theory, image processing, and numerical computation due to their valuable properties such as numerical computation stability, local effects of coefficient changes, and built-in smoothness between adjacent polynomial pieces.

The spline methods, as is known, provide inaccurate solutions with the presence of singularity. To defeat the drawback of these methods, we, in this work, develop an effective method based on quintic B-spline functions, known as the quintic B-spline method (QBSM), to approximate the solution of (1). To construct the QBSM, the approximate solution is forced to fulfill the considered system at the grid points, converting it into a set of algebraic equations with unknown coefficients. Solving the set of algebraic equations determines the values of these coefficients. Note that the considered problem has a singularity at $\tau = 0.$ When addressing the singularity of (1) numerically, it is important to efficiently deal with the singularity via certain means. In our case, we employ the L'Hôpital rule to its second term. To the best of our knowledge, the results presented in this work are new and have not been previously presented in the literature. The method is illustrated with several test problems. It is demonstrated that the accuracy of the method is of fourth-order convergence, superior to the convergence of the cubic B-spline method, which is proven to be of second-order convergence, derived in our prior work ^[8]. Outcomes are compared with some other numerical solutions to demonstrate the advantage of the method.

The structure of this paper is as follows: Section 2 provides the preliminaries of the quintic B-spine functions and their properties. Section 3 is dedicated to the construction of QBSM for obtaining the solution of the considered system. Section 4 discusses the convergence of the method. Section 5 provides the numerical illustration, and, finally, Section 6 summarizes and concludes our work.

2.   Preliminaries

In this section, we define quintic B-spline functions and their main properties to be utilized in constructing QBSM. We construct this method upon a uniform mesh. To do this, we partition the solution domain $\mathrm{\Gamma } = [\alpha, \beta]$ into $k$ subintervals ${\mathrm{\Gamma }}_{i} = [{\tau }_{i}, {\tau }_{i+1}$] by the grid points ${\tau }_{i} = \alpha +i\mathrm{\Lambda }$ $(i = 0, 1, ..., k)$, where $\mathrm{\Lambda } = (\beta -\alpha)/k.$ Let $\mathrm{\Omega }$ be the set of these grid points of the solution domain $\mathrm{\Gamma }$, referred to as the partition of $\mathrm{\Gamma }$, and is defined as $\mathrm{\Omega } = \left\{{\tau }_{0}, {\tau }_{1}, \dots, {\tau }_{k}\right\}.$ To provide proper support for the quintic B-spline functions, it is essential to introduce an additional five grid points on each side of the solution domain $\mathrm{\Gamma }$. Consequently, the solution domain $\mathrm{\Gamma }$ is extended to $\stackrel{-}{\mathrm{\Gamma }} = [\alpha -5\mathrm{\Lambda }, \beta +5\mathrm{\Lambda }]$ with ${\tau }_{i} = \alpha +i\mathrm{\Lambda }$ $\left(i = -5, ..., k+5\right).$ The linear space of quintic splines over this defined partition is expressed as

where ${\mu \left(\tau \right)|}_{{\mathrm{\Gamma }}_{i}}$ indicates the restriction of $\mu \left(\tau \right)$ over ${\mathrm{\Gamma }}_{i}$ and ${P}_{5}$ designates the set of one-variable quintic polynomials. The dimension of the linear space ${\mathrm{{\rm M}}}_{5}\left(\mathrm{\Gamma }\right)$ is $(k+5).$ According to ^[30], the quintic B-spline ${\mathrm{{\rm K}}}_{r}\left(\tau \right)(r = -2, ..., k+2)$ is defined as

The basis functions ${\mathrm{{\rm K}}}_{r}$, $r = -2, ..., k+2$, are nonnegative and linearly independent on the domain $\left[\alpha, \beta \right].$ The values of ${\mathrm{{\rm K}}}_{r}\left(\tau \right)$ and their derivatives up to the third order at the grid points are recorded in Table 1.

For an appropriately smooth function $\omega \left(\tau \right)$, one can uniquely define a quintic spline

that fulfill the interpolation conditions $\mu \left({\tau }_{i}\right) = \omega \left({\tau }_{i}\right), i = 0, ..., k,$ and ${\mu }^{\text{'}}\left(\alpha \right) = {\omega }^{\text{'}}\left(\alpha \right).$ From Table 1, for the discretization knots ${\tau }_{j}(j = 0, ..., k)$, we get

Equations (3)–(6) serve as the fundamental relations in the construction of the QBSM.

3.   Construction and convergence analysis of the method

This part of the study discusses the method and the convergence analysis.

3.1.   Construction of QPSM

In this section, we present the development of a collocation method based on quintic B-spline functions for (1) and (2). Let ${\mu }_{i}\left(\tau \right) = \sum _{r = -2}^{k+2}{{\mathrm{\lambda }}_{i, r}\mathrm{{\rm K}}}_{r}\left(\tau \right), i = 1, 2,$ represents the quintic B-spline approximate solution of the exact solution ${\omega }_{i}\left(\tau \right)$ to (1). To overcome the singularity behavior of (1), we employ the L'Hôpital rule on the second term at $\tau = 0,$ to obtain

Discretizing (7), we get

where $j = 1, \cdots, k.$ Using (3)–(5), we have

where $i = 1, 2$ and $j = 1, \cdots, k.$ Additionally, from (2), we derive the following four equations

where $i = 1, 2.$ Four equations are still needed. Therefore, by differentiating (1), we obtain:

Applying the L'Hôpital rule and using (2)–(4), (6) and (12) becomes

Similarly, at $\tau = 1$, we obtain

where $i = 1, 2$.

Expressing (9)–(11), (13), and (14) in matrix form as

where A represents a coefficient matrix of dimension $2(k+5)\times 2(k+5)$, $\mathrm{\Phi }$ is a column vector defined as

and $\mathrm{\Psi }$ is a column vector with $2(k+5)$ entries. Solving this system yields the coefficients of the approximate solution ${\mu }_{i}\left(\tau \right)$ for (1).

3.2.   Convergence analysis

In this section, we demonstrate the convergence analysis of QBSM. To facilitate this analysis, we assume that ${\omega }_{i}\left(\tau \right)\in {C}^{5}\left[\mathrm{0, 1}\right], i = \mathrm{1, 2}.$ From (3)–(6), we have ^[35,36]

With the operator notations, $\mathrm{\Xi }{\omega }_{i}\left({\tau }_{j}\right) = {\omega }_{i}\left({\tau }_{j+1}\right), D{\omega }_{i}\left({\tau }_{j}\right) = {\omega }_{i}^{\text{'}}\left({\tau }_{j}\right)$, and $I{\omega }_{i}\left({\tau }_{j}\right) = {\omega }_{i}\left({\tau }_{j}\right)$, Eq (16) can be expressed as

$i = \mathrm{1, 2}.$ Setting $\mathrm{\Xi } = {e}^{\mathrm{\Lambda }D}$ in (17) gives

Expanding the exponential functions in (18) in powers of $\mathrm{\Lambda }D$, we obtain

$i = \mathrm{1, 2}.$ Next, let's define truncation error as follows

As ${\mu }_{i}\left({\tau }_{j}\right) = {\omega }_{i}\left({\tau }_{j}\right), i = 1, 2$ and $j = 1, ..., k,$ Eq (20) can be simplified as

Hence, by using (19) in (21), we can conclude that

and for j = 0, we have

As ${\mu }_{i}\left({\tau }_{0}\right) = {\omega }_{i}\left({\tau }_{0}\right), i = 1, 2,$ Eq (23) can be simplified as

Hence, by using (19) in (24), we find

In light of (22) and (25), it can be deduced that the truncation error for the Lane-Emden system is of the order $O\left({\mathrm{\Lambda }}^{4}\right).$

4.   Numerical results

In this section, five test problems are considered to demonstrate the accuracy and applicability of QBSM. Additionally, the obtained numerical results corresponding to the considered system have been compared with those achieved previously ^[4,8,37]. Note that, in our calculations, "$\mathrm{E}-n$" means ${10}^{-n}.$

The absolute error (${Abs}_{i}$) and ${L}_{\infty }$ error are defined by

where ${\omega }_{i}\left(\tau \right)$ and ${\mu }_{i}\left(\tau \right)$ represent the exact and QBSM solutions at the grid point ${\tau }_{j}$, respectively. Moreover, the order of convergence (OC) of the method is computed by applying the following formula:

Problem 1. Consider the following system

subject to

The exact solution for this system is ${\omega }_{1}\left(\tau \right) = 1+{\tau }^{2}, {\omega }_{2}\left(\tau \right) = 1-{\tau }^{2}.$

We apply the proposed QBSM to solve this problem for $\mathrm{\Lambda } = 0.1$. Table 2 presents the exact and approximate solutions at the grid points. It is worth mentioning that, for this problem, the outcomes are exact and the errors are only incurred caused by round-off errors in computational processes.

Problem 2. Consider the following system

subject to

The exact solution of this system is

We apply the proposed QBSM to solve this problem for $\mathrm{\Lambda } = 0.1, 0.01$. The logarithmic plots of absolute errors for various values of $k$ are depicted in Figure 1, which exhibits that if the value of $k$ is increased, the error decreases. The absolute errors obtained by QBSM are given in Tables 3 and 4 along with those obtained by CBSM ^[8] and Chebyshev operational matrix method (COMM) ^[4]. Comparison reveals that QBSM yields more accurate solutions than the methods in ^[4,8]. The outcomes of ${L}_{\infty }^{i}\left(k\right)$ errors are listed using $k = 16, 32, 64,$ and $128$. In addition, the ${OC}^{i}, i = \mathrm{1, 2},$ are computed and the results are tabulated in Table 5. It can be observed that the achieved ${OC}^{i}, i = \mathrm{1, 2},$ is four. The method's computational time (CPU time) is reported in the same Table, which confirms that the QBSM is computationally effective.

Problem 3. Consider the following system

subject to

where the exact solution is

We apply the proposed QBSM to obtain the approximate solutions to this problem for $\mathrm{\Lambda } = 0.1, 0.01$. Absolute errors of QBSM for $\mathrm{\Lambda } = 0.1, 0.01$ are listed in Tables 5 and 6, respectively, along with those obtained by the CBSM ^[8]. From these tables, it can be observed that QBSM provides lesser error than CBSM. The logarithmic plots of absolute errors for various values of $k$ are depicted in Figure 2. The outcomes of ${L}_{\infty }^{i}\left(k\right)$ errors are listed using $k = 16, 32, 64,$ and $128$. In addition, the ${OC}^{i}, i = \mathrm{1, 2},$ are computed and the results are tabulated in Tables 7 and 8. The table show that the achieved ${OC}^{i}, i = \mathrm{1, 2},$ is four. The method's CPU time is reported in the same table, which confirms that the QBSM is computationally effective.

Problem 4. Consider the following system of LEE

subject to

where the exact solution is given by ${\omega }_{1}\left(\tau \right) = \sqrt{{1+\tau }^{2}}, {\omega }_{2}\left(\tau \right) = \frac{1}{\sqrt{{1+\tau }^{2}}}.$ We solve this system using the proposed QBSM for $\mathrm{\Lambda } = 0.1, 0.01$. Absolute errors obtained by QBSM for $\mathrm{\Lambda } = 0.1, 0.01$ are given in Tables 9 and 10, along with the errors obtained by the CBSM ^[8] and COMM ^[4]. From these tables, it seems that the errors of QBSM are less than the errors of CBSM and COMM. The logarithmic graphs of absolute errors for different values of n are presented in Figure 3. The outcomes of ${L}_{\infty }^{i}\left(k\right)$ errors are listed using $k = 16, 32, 64,$ and $128$. In addition, the ${OC}^{i}, i = \mathrm{1, 2},$ are computed and the results are tabulated in Table 11. The table show that the achieved ${OC}^{i}, i = \mathrm{1, 2},$ is four. The method's CPU time is reported in the same table, which confirms that the QBSM is computationally effective.

Problem 5. Consider the following system of LEE

subject to

The exact solution for this system is

We solve this system using the proposed QBSM for $\mathrm{\Lambda } = 0.1, 0.01$. Absolute errors obtained by QBSM for $\mathrm{\Lambda } = 0.1, 0.01$ are tabulated in Tables 12 and 13, along with the errors reported in CBSM ^[8] and Dickson operational matrix (DOM) ^[37]. We note that QBSM yields results more accurate than those obtained in ^[8,37]. The logarithmic graphs of absolute errors for different values of $k$ are displayed in Figure 4. The outcomes of ${L}_{\infty }^{i}\left(k\right)$ errors are listed using $k = 16, 32, 64,$ and $128$. In addition, the ${OC}^{i}, i = \mathrm{1, 2},$ are computed and the results are tabulated in Table 14. The table show that the achieved ${OC}^{i}, i = \mathrm{1, 2},$ is four. The method's CPU time is reported in the same Table, which confirms that the QBSM is computationally effective.

As can be observed from the above tables, the proposed QBSM is fourth-order accurate and the practical convergence order aligns consistently with the theoretical convergence order obtained in the previous section.

5.   Conclusions

In this study, we have established a numerical method for solving systems of Lane-Emden equations. The QBSM has been constructed using quintic B-spline functions on the uniform mesh. We investigate the convergence analysis of the QBSM and found it exhibited fourth-order convergence. To strengthen the significance of the QBSM method and validate theoretical results, we examined five test problems. We have presented tabular and graphical exhibitions to confirm the effectiveness of QBSM. Notably, the numerical solutions of QBSM are in good agreement with the exact ones, and their accuracy improves as the step sizes decrease. Moreover, we compared the QBSM with other numerical methods such as CBSM, DOM, and COMM, and the comparison exposed that the QBSM produces more accurate numerical results than the other methods. In conclusion, the method is computationally efficient, accurate, robust, easy to address the singularity, and, therefore, it can be employed to solve different classes of nonlinear singular differential equations.

Use of AI tools declaration

The authors declare that they have not used Artificial Intelligence (AI) tools in the creation of this article.

Conflict of interest

The authors declare no conflicts of interest.