High-Order Shifted Gegenbauer Integral Pseudospectral Method for Solving Differential Equations of Lane-Emden Type

TL;DR

This paper introduces a high-order, exponentially convergent shifted Gegenbauer integral pseudospectral method for efficiently solving Lane-Emden differential equations with mixed boundary conditions, demonstrating superior accuracy and stability.

Contribution

The paper develops a novel integral pseudospectral method using shifted Gegenbauer polynomials, including new algorithms and rigorous error analysis for solving Lane-Emden equations.

Findings

Exhibits exponential convergence and high accuracy.

Demonstrates effectiveness through five numerical examples.

Outperforms existing methods in stability and efficiency.

Abstract

We present a novel, high-order, efficient, and exponentially convergent shifted Gegenbauer integral pseudospectral method (SGIPSM) to solve numerically Lane-Emden equations provided with some mixed Neumann and Robin boundary conditions. The framework of the proposed method includes: (i) recasting the problem into its integral formulation, (ii) collocating the latter at the shifted flipped-Gegenbauer-Gauss-Radau (SFGGR) points, and (iii) replacing the integrals with accurate and well-conditioned numerical quadratures constructed via SFGGR-based shifted Gegenbauer integration matrices. The integral formulation is eventually discretized into linear/nonlinear system of equations that can be solved easily using standard direct system solvers. The implementation of the proposed method is further illustrated through four efficient computational algorithms. Moreover, we furnish rigorous error and convergence analyses of the SGIPSM. Five numerical test examples are presented to verify the effectiveness, accuracy, exponential convergence, and numerical stability of the proposed method. The numerical simulations are associated with extensive numerical comparisons with other rival methods in the literature to demonstrate further the power of the proposed method. The SGIPSM is broadly applicable and represents a strong addition to common numerical methods for solving linear/nonlinear differential equations when high-order approximations are required using a relatively small number of collocation points.

Tables

Table 1. Table 1: The REs of the B-spline method [ Caglar and Caglar (2006) ] and the SGIPSM.

$x$	B-spline method [Caglar and Caglar (2006)]	SGIPSM
	$n=20$	$n=5,\alpha =0.1$
0	9.8160e-005	1.8405e-006
0.05	9.8756e-005	1.8512e-006
0.1	1.0122e-004	3.2213e-006
0.2	1.0231e-004	9.0436e-006
0.3	1.0119e-004	3.6570e-006
0.4	1.0425e-004	7.1760e-006
0.5	1.0616e-004	1.3045e-005
0.6	1.0911e-004	6.5036e-006
0.7	1.0925e-004	1.1139e-005
0.8	1.1398e-004	2.5024e-005
0.9	1.1117e-004	2.4411e-006
AE1	0	0

Table 2. Table 2: The REs of Bessel’s method [ Yüzbaşı and Sezer (2013) ] and the SGIPSM using eight collocation points.

$x$	Bessel’s method [Yüzbaşı and Sezer (2013)]	SGIPSM
		$\alpha =1.1$
0	5.9957e-006	7.9499e-008
0.2	6.0946e-006	8.6474e-009
0.4	6.8724e-006	5.4461e-010
0.6	8.8719e-006	2.8767e-008
0.8	1.5473e-005	4.0944e-008
AE1	1.6198e-016	0

Table 3. Table 3: The REs of the Cubic spline method [ Kanth (2007) ] and the SGIPSM

$x$	Cubic spline method [Kanth (2007)]	SGIPSM
	$n=50$	$n=7,\alpha =1.4$
0	2.2341e-006	1.1021e-007
0.1	2.1616e-006	3.7718e-009
0.2	1.9619e-006	1.7350e-008
0.3	1.6583e-006	8.4623e-008
0.4	1.2894e-006	3.6471e-008
0.5	9.0037e-007	1.6719e-008
0.6	5.3675e-007	3.4350e-008
0.7	2.3891e-007	7.9305e-008
0.8	3.8253e-008	4.9866e-009
0.9	4.4620e-008	1.9892e-009
AE1	1.3323e-015	4.4409e-016

Table 4. Table 4: The MAEs of [ Turkyilmazoglu (2013) ] ’s method and the SGIPSM.

$n$	[Turkyilmazoglu (2013)]’s method	SGIPSM
		MAE/($\alpha$)
3	1.1100e-003	3.9711e-004/($\alpha =0.8$)
6	5.5622e-006	1.7118e-006/($\alpha =-0.1$)
8	5.2440e-008	2.6347e-008/($\alpha =0.8$)

Table 5. Table 5: The REs of the cubic-spline method [ Kanth and Reddy (2005) ] and the SGIPSM.

$x$	Cubic spline method [Kanth and Reddy (2005)]	SGIPSM
	$n=20$	$n=6,\alpha =-0.2$
0.25	-	1.4056e-008
0.05	8.9610e-005	2.2797e-007
0.075	-	4.4184e-007
0.1	8.9087e-005	5.7829e-007
0.2	8.6627e-005	1.9466e-007
0.3	8.2151e-005	7.0896e-007
0.4	7.6256e-005	7.5662e-007
0.5	6.8275e-005	7.4427e-008
0.6	5.8753e-005	7.2706e-007
0.7	4.7165e-005	4.2395e-007
0.8	3.3711e-005	3.3767e-007
0.9	1.7861e-005	3.0010e-007
AE1	0	0

Table 6. Table 6: The REs of the SGIPSM using n = 5 𝑛 5 n=5 .

$x$	SGIPSM
	$\alpha =0.9$
0	2.9203e-005
0.3	5.3820e-005
0.6	9.1514e-005
0.9	7.3363e-005
1.2	5.2870e-005
AE1.5	4.996e-016

Table 7. Table 7: The MAEs of [ Khuri and Sayfy (2010) ] ’s method and the SGIPSM using n = 5 ​ ( 5 ) ​ 20 𝑛 5 5 20 n=5(5)20 .

$n$	[Khuri and Sayfy (2010)]’s method	SGIPSM
		MAE/($\alpha$)
5	2.37e-005	1.8012e-005/(0.9)
10	6.18e-006	1.3886e-009/(0.5)
15	3.03e-006	1.1013e-013 /(-0.1)
20	1.56e-006	5.8287e-015 /(2)

Table 8. Table 8: The REs of the SGIPSM using n = 7 𝑛 7 n=7 .

$x$	SGIPSM
	$\alpha =0.9$
0	2.8272e-016
0.1	1.6607e-015
0.2	2.0410e-014
0.3	5.2419e-015
0.4	1.2707e-014
0.5	4.1473e-015
0.6	1.1436e-015
0.7	2.8049e-015
0.8	5.0413e-015
0.9	3.3102e-015
AE1	4.4409e-016

Equations

$$y^{\prime \prime }(x)+\frac{{\alpha }_2}{x}{y}^{\prime }(x)+f(x,y)=0,x\in (0,b],$$

$$y^{\prime \prime }(x)+\frac{{\alpha }_2}{x}{y}^{\prime }(x)+p(x)y(x)=g(x),x\in (0,b],$$

$$y^{\prime }(0)={\alpha }_1,$$

$$\beta y(b)+\gamma y^{\prime }(b)=\delta ,$$

$$(n+2\alpha )G_{n+1}^{(\alpha )}(x)=2(n+\alpha )x{G}_n^{(\alpha )}(x)-n{G}_{n-1}^{(\alpha )}(x),n=1,2,\dots ,$$

$$K_n^{(\alpha )}=2^{n-1}\frac{\Gamma (n+\alpha )\Gamma (2\alpha +1)}{\Gamma (n+2\alpha )\Gamma (\alpha +1)},n=0,1,\dots .$$

$${\int }_{-1}^1{G}_m^{(\alpha )}(x)G_n^{(\alpha )}(x)w^{(\alpha )}(x)dx={\lambda }_n^{(\alpha )}{\delta }_{m,n},$$

$${\lambda }_n^{(\alpha )}=\frac{2^{2\alpha -1}n!{\Gamma }^2(\alpha +\frac{1}{2})}{(n+\alpha )\Gamma (n+2\alpha )},$$

$${\varphi }_j^{(\alpha )}(x)=\frac{G_j^{(\alpha )}(x)}{\sqrt{{\lambda }_j^{(\alpha )}}},j=0,\dots ,n,$$

$$\sum _{j=0}^n{\varpi }_{n,j}^{(\alpha )}{\varphi }_s^{(\alpha )}\left(x_{n,j}^{(\alpha )}\right){\varphi }_k^{(\alpha )}\left(x_{n,j}^{(\alpha )}\right)={\delta }_{sk},s,k=0,\dots ,n,$$

$${\varpi }_{n,0}^{(\alpha )}$$

$${\varpi }_{n,j}^{(\alpha )}$$

$${\aleph }_{n,j}^{\left(\alpha \right)}=2^{2\alpha -1}\frac{{\Gamma }^2\left(\alpha +\frac{1}{2}\right)n!}{\left(n+\alpha +\frac{1}{2}\right)\Gamma (n+2\alpha +1)}\left(1+x_{n,j}^{(\alpha )}\right){\left(G_n^{(\alpha )}\left(x_{n,j}^{(\alpha )}\right)\right)}^{-2},j=0,1,\dots ,n.$$

$$f(x)\approx P_{n}f(x)=\sum _{i=0}^n{a}_i{G}_i^{(\alpha )}(x),$$

$${\mathbf{I}}_n^{(\alpha )}={\left({\int }_{-1}^{x_{n,0}^{(\alpha )}}{P}_{n}f(x)dx,{\int }_{-1}^{x_{n,1}^{(\alpha )}}{P}_{n}f(x)dx,\dots ,{\int }_{-1}^{x_{n,n}^{(\alpha )}}{P}_{n}f(x)dx\right)}^T,$$

$${\mathbf{I}}_n^{(\alpha )}={\mathbf{Q}}^{(1)}\mathbf{F}.$$

$${\int }_{-1}^{x_{n,i}^{(\alpha )}}f(x)dx=\sum _{k=0}^n{Q}_{i,k}^{(1)}{f}_{n,k}^{(\alpha )}+{}_f{E}_n^{(1)}\left(x_{n,i}^{(\alpha )},{\xi }_{n,i}^{(1)}\right)\forall i,$$

$$Q_{i,k}^{(1)}=\sum _{j=0}^n{\left({\lambda }_j^{(\alpha )}\right)}^{-1}{\varpi }_{n,k}^{(\alpha )}{G}_j^{(\alpha )}\left(x_{n,k}^{(\alpha )}\right){\int }_{-1}^{x_{n,i}^{(\alpha )}}{G}_j^{(\alpha )}(x)dx\forall i,k,$$

$${}_f{E}_n^{(1)}\left(x_{n,i}^{(\alpha )},{\xi }_{n,i}^{(1)}\right)=\frac{f^{(n+1)}\left({\xi }_{n,i}^{(1)}\right)}{(n+1)!K_{n+1}^{(\alpha )}}{\int }_{-1}^{x_{n,i}^{(\alpha )}}{q}_n^{(\alpha )}(x)dx\forall i,$$

$$\sum _{j=0}^n{\varpi }_{n,j}^{(\alpha )}{G}_s^{(\alpha )}\left(x_{n,j}^{(\alpha )}\right)f_{n,j}^{(\alpha )}$$

$$=\sum _{k=0}^n{a}_k\sum _{j=0}^n{\varpi }_{n,j}^{(\alpha )}\sqrt{{\lambda }_s^{(\alpha )}{\lambda }_k^{(\alpha )}}{\varphi }_s^{(\alpha )}\left(x_{n,j}^{(\alpha )}\right){\varphi }_k^{(\alpha )}\left(x_{n,j}^{(\alpha )}\right),$$

$$\stackrel{\text{Eq. }\text{eqref}eq:13}{=}\sum _{k=0}^n{a}_k\sqrt{{\lambda }_s^{(\alpha )}{\lambda }_k^{(\alpha )}}{\delta }_{sk}=a_s{\lambda }_s^{(\alpha )},s=0,\dots ,n.$$

$$\Rightarrow a_s=\frac{1}{{\lambda }_s^{(\alpha )}}\sum _{j=0}^n{\varpi }_{n,j}^{(\alpha )}{G}_s^{(\alpha )}\left(x_{n,j}^{(\alpha )}\right)f_{n,j}^{(\alpha )}\forall s.$$

$$P_{n}f(x)=\sum _{k=0}^n{f}_{n,k}^{(\alpha )}{L}_{n,k}^{(\alpha )}(x),$$

$$L_{n,k}^{(\alpha )}(x)={\varpi }_{n,k}^{(\alpha )}\sum _{j=0}^n\frac{1}{{\lambda }_j^{(\alpha )}}{G}_j^{(\alpha )}\left(x_{n,k}^{(\alpha )}\right)G_j^{(\alpha )}(x)\forall k.$$

$$f(x)=\sum _{k=0}^n{f}_{n,k}^{(\alpha )}{L}_{n,k}^{(\alpha )}(x)+{}_f{E}_n\left(x,\xi \right)\forall x\in [-1,1],$$

$${}_f{E}_n\left(x,\xi \right)=\frac{f^{(n+1)}\left(\xi \right)}{(n+1)!}\prod _{k=0}^n\left(x-x_{n,k}^{(\alpha )}\right).$$

$${\int }_{-1}^{x_{n,i}^{(\alpha )}}{\int }_{-1}^{x}f(t)dtdx=\sum _{k=0}^n{Q}_{i,k}^{(2)}{f}_{n,k}^{(\alpha )}+{}_f{E}_n^{(2)}\left(x_{n,i}^{(\alpha )},{\xi }_{n,i}^{(2)}\right)\forall i,t\in [-1,1],$$

$$Q_{i,k}^{(2)}=\left(x_{n,i}^{(\alpha )}-x_{n,k}^{(\alpha )}\right)Q_{i,k}^{(1)}\forall i,k,$$

$${}_f{E}_n^{(2)}\left(x_{n,i}^{(\alpha )},{\xi }_{n,i}^{(2)}\right)=\frac{\left(x_{n,i}^{(\alpha )}-{\xi }_{n,i}^{(2)}\right)f^{(n+1)}\left({\xi }_{n,i}^{(2)}\right)-(n+1)f^{(n)}\left({\xi }_{n,i}^{(2)}\right)}{(n+1)!K_{n+1}^{(\alpha )}}{\int }_{-1}^{x_{n,i}^{(\alpha )}}{q}_n^{(\alpha )}(x)dx\forall i.$$

$$x_{b,n,k}^{(\alpha )}=\frac{b}{2}\left(x_{n,k}^{(\alpha )}+1\right)\forall k.$$

$$f(x)\approx \sum _{i=0}^n{a}_i{G}_{b,i}^{(\alpha )}(x),$$

Full text

1 Introduction

Lane-Emden equations describe various physical phenomena in mathematical physics and astrophysics such as thermal explosions, stellar structure, thermal behavior of a spherical cloud of gas, thermionic currents, radiative cooling, etc.; cf. [Chandrasekhar (1958), Harley and Momoniat (2007), Singh et al. (2009), Kanth and Aruna (2010), Yüzbaşı and Sezer (2013), Tohidi et al. (2013), Abd-Elhameed (2015)]. They arise also in chemistry, chemical engineering, and other branches to govern the concentrations of oxygen and the carbon substrate [Wazwaz et al. (2016)], deriving analytical expressions for effectiveness factors $\eta$ in non-isothermal, spherical catalysts [Scott et al. (1984)], etc.

In the present work, we are concerned with the numerical solution of the following two forms of nonlinear and linear Lane-Emden equations

[TABLE]

[TABLE]

respectively, where $y$ is the unknown solution function, $f$ is a given nonlinear smooth function, $p(x)$ and $g(x)$ are some given sufficiently smooth functions, $\alpha_{2}\in\mathbb{R}$ and $b\in\mathbb{R}^{+}$ are some given constants. Both equations are provided with the following mixed Neumann and Robin boundary conditions

[TABLE]

for some given real numbers $\alpha_{1},\beta,\gamma$, and $\delta$. The real parameters $\beta$ and $\gamma$ in the boundary condition (1.3b) do not vanish simultaneously. We shall refer to Eqs. (1.1), (1.3a), and (1.3b) by Problem 1. Also, we refer to Eqs. (1.2), (1.3a), and (1.3b) by Problem 2.

The theoretical picture for second-order initial-value problems of Lane-Emden type equations is well-known. Local existence, global existence, and uniqueness of solutions were proven by [Biles et al. (2002)]. For problem 1 with $f(x,y)=y^{m}$ and $\alpha_{2}=2$, the exact solution is available only for $m=0,1$, and $5$. For other values of $m$, the problem can only be solved numerically [Mohan and Al-Bayaty (1980), Parks (1984)]; therefore, the problem of finding solutions for other values of $m$ is very important from both theoretical and practical viewpoints. This topic captivated the interest of many researchers through many decades in attempt to derive new efficient computational algorithms to find accurate numerical solutions for these types of singular differential equations. In the following, we mention a few: [Seidov and Kuzakhmedov (1977), Mohan and Al-Bayaty (1980), Parks (1984)] approached the numerical solution using power series methods. [Horedt (1986)] published an article later listing seven-digit tables of the numerical solutions of the Lane-Emden equation for the plane-parallel, cylindrical, and spherical case using a Runge-Kutta method– a work considered a succession to the works of [Sadler and Miller (1932)] and [Chandrasekhar and Wares (1949)] without the aid of electronic computers. [Shawagfeh (1993)] used a coalition of the Adomian decomposition method and the Padé approximants method to accelerate the convergence of the power series. A variational approach for solving Lane-Emden equations was presented by [He (2003)]. [Ramos (2004)] handled the Lane-Emden equation using a piecewise quasilinearization method that produces piecewise analytical solutions calculated through the analytic integration of the reduced linear constant-coefficients ordinary differential equation. A rational Legendre tau method was employed in the same year by [Parand and Razzaghi (2004)] for solving special initial-value problems of Lane-Emden type equations on a semi-infinite interval. [Momoniat and Harley (2006)] applied Lie group analysis to achieve a larger radius of convergence than the power series solution of nonlinear Lane-Emden equations with $f(x,y)=e^{y}$ and $\alpha_{2}=2$. Differential transformation method based on Taylor series expansion for solving singular initial-value problems of Lane-Emden type equations was introduced by [Ertürk (2007)]. [Dehghan and Shakeri (2008)] studied the Lane-Emden equation using variational iteration method. [Marzban et al. (2008)Marzban, Tabrizidooz, and Razzaghi] applied a hybridization of block-pulse functions and Lagrange interpolating polynomials to reduce the computation of nonlinear initial-value problems to a system of non-algebraic equations. [Yang and Hou (2010)] used a spectral method based on a truncated series of Chebyshev expansions to solve singular initial value problems of Lane-Emden type equations. Bernstein polynomial approximations were employed by [Pandey and Kumar (2012)] to solve linear and nonlinear Lane-Emden equations. [Bhrawy and Alofi (2012)] treated singular initial-value problems of nonlinear Lane-Emden type equations defined on a semi-infinite domain using shifted Jacobi polynomials for the spatial approximation integrated with shifted Jacobi-Gauss points as collocation nodes. [Yüzbaşı and Sezer (2013)] presented a predictor-corrector algorithm to solve linear Lane-Emden equations provided with certain mixed conditions. In particular, the predicted solution is computed using a Bessel collocation method, and then correctified by solving an additional error problem that is constructed using a residual error function. [Öztürk and Gülsu (2014)] applied a collocation method based on Hermite polynomial approximations. [Šmarda and Khan (2015)] presented a computational method integrating both the differential transformation method and a modified general formula for the Adomian polynomials for solving singular initial-value problems of nonlinear Lane-Emden type equations. Recently, [Calvert et al. (2016)] used modified rational Bernoulli function approximations to solve nonlinear Lane-Emden equations on the semi-infinite domain.

To avoid the ill-conditioning of differential operators, an alternative direction to the aforementioned methods is to recast the differential equation into its integral formulation to take advantage of the well-conditioning of integral operators, then discretize the latter using various discretization techniques. Perhaps the first approach in this direction for solving Lane-Emden equations was presented by [Yousefi (2006)] in 2006, using Legendre wavelet approximations and the Gaussian integration method. A similar approach was later presented by [Karimi Vanani and Aminataei (2010)] in which the transformed integral equation was eventually converted into Padé series form. [Wazwaz et al. (2013)] applied the Adomian decomposition method to handle the integral form of the Lane-Emden equations. The Adomian decomposition method was later applied by [Rach et al. (2014)] on some alternate derivations for the Volterra integral forms of the Lane-Emden equation.

Motivated by the abundant advantages and large success of pseudospectral methods in broad scientific areas and applications [Fornberg (1998), Boyd (2001), Mason and Handscomb (2002), Elgindy (2017), Suh (2009), Ng and Yan (2003), Xu et al. (2016), Elgindy (2016), Xiao et al. (2016)], the exceedingly accurate and stable approximations obtained via integral reformulations [Du (2016), Elgindy (2016), Elgindy and Smith-Miles (2013a), Elgindy and Smith-Miles (2013b), Greengard (1991)], and the useful fundamental properties of Gegenbauer basis polynomials [Elgindy (2013), Keiner (2009)], we present the shifted Gegenbauer integral pseudospectral method (SGIPSM): a high-order, efficient, and exponentially convergent pseudospectral method to numerically solve Problems 1 and 2. The proposed method employs shifted Gegenbauer polynomial expansions to directly discretize the integral forms of Problems 1 and 2 on $(0,b]$, and adopts collocations in the nodal space at the shifted flipped-Gegenbauer-Gauss-Radau (SFGGR) points to overcome the singularity of the Lane-Emden equation at the origin that reflects the main difficulty, and to accurately impose the boundary condition at $x=b$. Moreover, instead of directly integrating both sides of Eqs. (1.1) and (1.2), the proposed method follows closely the method of [Elgindy and Smith-Miles (2013b)] by introducing a useful substitution to the second-order solution derivative, solving Problems 1 and 2 in terms of the new substituting function, and then stably recovering the original solution function via successive integration. All necessary integral evaluations are accomplished using some newly developed SFGGR-based shifted Gegenbauer integration matrices. This approach entails the approximation of the 1- and 2-fold integrals of the solution and substituting functions, whereas the former approach may require the evaluation of the 1- and 2-fold integrals of highly nonlinear and intricate integrands.

The rest of the article is organized as follows: in Section 2, we begin by briefly reviewing some basic preliminaries concerned with Gegenbauer polynomials. In Section 3, we show how to construct the SFGGR-based shifted Gegenbauer integration matrices from the standard Gegenbauer integration matrices in addition to some error analysis results. In Section 4, we introduce the SGIPSM for solving Problems 1 and 2. Section 5 is devoted for the study of error bounds and convergence analysis of the SGIPSM. Section 6 is designed for numerical simulations and comparisons in which five numerical test examples are studied followed by some final concluding remarks and a discussion in Section 7. Appendix A is reserved for the presentation of four efficient computational algorithms to implement the SGIPSM.

2 Preliminaries

In this section, we present some prior properties of the Gegenbauer polynomials. We rely heavily on the useful standardization given by [Elgindy and Smith-Miles (2013), Eq. (A.2)].

The Gegenbauer polynomial, $G_{n}^{(\alpha)}(x)$, of degree $n\in\mathbb{Z}^{+}_{0}=\mathbb{Z}^{+}\cup\{0\}$, and associated with the parameter $\alpha>-1/2$, is a real-valued function that appears as an eigensolution to the singular Sturm-Liouville problem in the interval $[-1,1]$ [Elgindy (2016)]. Chebyshev polynomials of the first kind and Legendre polynomials are special cases of this family of orthogonal polynomials for $\alpha=0$ and $0.5$, respectively. Gegenbauer polynomials can be easily generated by the following recurrence relation

[TABLE]

starting with $G_{0}^{(\alpha)}(x)=1$ and $G_{1}^{(\alpha)}(x)=x$. The leading coefficient, $K_{n}^{(\alpha)}$, of $G_{n}^{(\alpha)}(x)$ is defined by

[TABLE]

Gegenbauer polynomials are orthogonal with respect to the weight function $w^{(\alpha)}(x)=\left(1-x^{2}\right)^{\alpha-\frac{1}{2}}$, and their orthogonality relation is defined by,

[TABLE]

where

[TABLE]

is the normalization factor, and $\delta_{m,n}$ is the Kronecker delta function. We denote the zeros of the $(n+1)$th-degree polynomial $q_{n}^{(\alpha)}(x)=G_{n+1}^{(\alpha)}(x)-G_{n}^{(\alpha)}(x)$ (aka FGGR nodes) sorted in descending order by $x_{n,k}^{(\alpha)},k=0,1,\ldots,n$, and denote their set by $\mathbb{S}_{n}^{(\alpha)}$. The orthonormal Gegenbauer basis polynomials are defined by

[TABLE]

and they satisfy the discrete orthonormality relation

[TABLE]

where $\varpi_{n,j}^{(\alpha)},j=0,1,\ldots,n$, are the corresponding Christoffel numbers of the FGGR quadrature formula on the interval $[-1,1]$ that are defined by

[TABLE]

[TABLE]

3 The SFGGR-based shifted Gegenbauer integration matrix

Suppose that we approximate a real-valued function $f(x)\in C^{n+1}[-1,1]$ by the following truncated series expansion of Gegenbauer basis polynomials

[TABLE]

where $a_{i},i=0,\ldots,n,$ are the Gegenbauer spectral coefficients. Denote the function values evaluated at the FGGR nodes, $f\left(x_{n,i}^{(\alpha)}\right)$, by $f_{n,i}^{(\alpha)}$, for all $i$. The first-order Gegenbauer integration matrix calculated at the FGGR nodes is a linear map, $\mathbf{Q}^{(1)}$, which takes a vector of $(n+1)$ function values, ${\bm{F}}={\left({{f_{n,0}^{(\alpha)}},{f_{n,1}^{(\alpha)}},\ldots,{f_{n,n}^{(\alpha)}}}\right)^{T}}$, to a vector of $(n+1)$ integral values

[TABLE]

such that

[TABLE]

Following the method presented by [Elgindy and Smith-Miles (2013)], one can construct the elements of $\mathbf{Q}^{(1)},0\leq i,k\leq n,$ by the following theorem.

Theorem 3.1**.**

Let $f(x)\in C^{n+1}[-1,1]$ be approximated by the Gegenbauer interpolant $P_{n}f$ given by Eq. (3.1), where the coefficients $a_{i},\;i=0,\ldots,n,$ are obtained by interpolating $f$ at the FGGR points. Then there exist a first-order integration matrix $\mathbf{Q}^{(1)}=\left(Q_{i,j}^{(1)}\right),0\leq i,j\leq n$, and some numbers $\xi_{n,i}^{(1)}\in(-1,1)$, $i=0,\ldots,n,$ such that

[TABLE]

where

[TABLE]

[TABLE]

and the weights $\varpi_{n,k}^{(\alpha)}$ are as defined by Eqs. (2.7).

Proof.

Since $f_{n,j}^{(\alpha)}=\sum\limits_{k=0}^{n}a_{k}\,G_{k}^{(\alpha)}\left(x_{n,j}^{(\alpha)}\right),\;j=0,\ldots,n$, then

[TABLE]

[TABLE]

Substituting Eq. (3.6) into Eq. (3.1) yields the Gegenbauer interpolant in Lagrange form as follows:

[TABLE]

where

[TABLE]

Therefore, we can write

[TABLE]

for some $\xi\in(-1,1)$, where ${}_{f}E_{n}$ is the interpolation truncation error at the FGGR points defined by

[TABLE]

The proof is established by realizing that $q_{n}^{(\alpha)}(x)=K_{n+1}^{(\alpha)}\prod_{k=0}^{n}\left(x-x_{n,k}^{(\alpha)}\right)$, and integrating Eq. (3.9) on $\left[-1,x_{n,i}^{(\alpha)}\right]$, for all $i$. ∎

Let us denote the $q$th-order Gegenbauer integration matrix by $\mathbf{Q}^{(q)}$, its $i$th row by $\mathbf{Q}_{i}^{(q)}$, for all $i$, and its entries by $Q_{i,k}^{(q)}$, for all $i,k$. The following theorem provides a useful means to calculate the 2-fold integral of $f$ through the second-order Gegenbauer integration matrix generated directly from the first-order Gegenbauer integration matrix.

Theorem 3.2**.**

Given the assumptions of Theorem 3.1, there exist a second-order integration matrix $\mathbf{Q}^{(2)}=\left(Q_{i,j}^{(2)}\right),0\leq i,j\leq n$, and some numbers $\xi_{n,i}^{(2)}\in(-1,1)$, $i=0,\ldots,n,$ such that

[TABLE]

where

[TABLE]

[TABLE]

Proof.

The proof can be easily derived using Cauchy’s formula for repeated integration. ∎

To discretize the nonlinear Lane-Emden equation (1.1) on the shifted domain $(0,b]$, the shifted forms of Gegenbauer basis polynomials and their associated integration matrices are needed. To this end, let $G_{b,n}^{(\alpha)}(x)$ denote the shifted Gegenbauer polynomial $G_{n}^{(\alpha)}(2x/b-1)$, for all $x\in[0,b]$. We denote the SFGGR nodes by $x_{b,n,k}^{(\alpha)},k=0,1,\ldots,n$, and denote their set by $\mathbb{S}_{b,n}^{(\alpha)}$. Clearly

[TABLE]

Following the above convention, let $\mathbf{Q}_{b}^{(q)},\mathbf{Q}_{b,i}^{(q)}$, and $Q_{b,i,k}^{(q)}$, for all $i$ and $k$ denote the $q$th-order shifted integration matrix, its $i$th row, and its entries, respectively, and denote $f\left(x_{b,n,i}^{(\alpha)}\right)$ by $f_{b,n,i}^{(\alpha)}$, for all $i$. The following two theorems highlight the generation of the first-and second-order shifted Gegenbauer integration matrices $\mathbf{Q}_{b}^{(1)}$ and $\mathbf{Q}_{b}^{(2)}$, and mark the truncation errors of their associated quadratures.

Theorem 3.3**.**

Let $f(x)\in C^{n+1}[0,b]$ be approximated by a truncated series expansion of shifted Gegenbauer basis polynomials of the following form:

[TABLE]

where $a_{i},i=0,\ldots,n,$ are the shifted Gegenbauer spectral coefficients obtained by interpolating $f$ at the SFGGR points. Then there exist a first-order integration matrix $\mathbf{Q}_{b}^{(1)}=\left(Q_{b,i,j}^{(1)}\right),\;0\leq i,j\leq n$, and some numbers $\xi_{b,n,i}^{(1)}\in(0,b),i=0,\ldots,n,$ such that

[TABLE]

where

[TABLE]

[TABLE]

$q_{b,n}^{(\alpha)}(x)=q_{n}^{(\alpha)}\left(2x/b-1\right)$, for all $x\in(0,b)$, and $K_{b,n}^{(\alpha)}={\left({2/b}\right)^{n}}K_{n}^{(\alpha)}$ is the leading coefficient of the shifted Gegenbauer polynomial $G_{b,n}^{(\alpha)}(x)$.

Proof.

Using the change of variable $x=b\left(y+1\right)/2$, and Theorem 3.1, we find that

[TABLE]

for some $\xi_{n,i}^{(1)}\in(-1,1)$, $i=0,\ldots,n$, where $\bar{f}(y)=f(b(y+1)/2)$, for all $y\in[-1,1]$, and $\xi_{b,n,i}^{(1)}=b\left(\xi_{n,i}^{(1)}+1\right)/2$, for all $i.$ ∎

Theorem 3.4**.**

Given the assumptions of Theorem 3.3, there exist a second-order integration matrix $\mathbf{Q}_{b}^{(2)}=\left(Q_{b,i,j}^{(2)}\right),\;0\leq i,j\leq n$, and some numbers $\xi_{b,n,i}^{(2)}\in(0,b)$,  $i=0,\ldots,n,$ such that

[TABLE]

where

[TABLE]

[TABLE]

Proof.

The proof follows similar to that of Theorem 3.2. ∎

In the rest of the article, by $f(\bm{x})$, we mean ${[f({x_{1}}),f({x_{2}}),\ldots,f({x_{n}})]^{T}}$, for any real-valued function $f:\bm{\Omega}_{1}\subseteq\mathbb{R}\to\mathbb{R}$, and $\bm{x}=[x_{1},x_{2},\ldots,x_{n}]^{T}\in\mathbb{R}^{n}$. Similarly, by $g(\bm{x},\bm{y})$, we mean ${[g({x_{1}},y_{1}),g({x_{2}},y_{2}),\ldots,g({x_{n}},y_{n})]^{T}}$, for any real-valued function $g:\bm{\Omega}_{2}\subseteq\mathbb{R}^{2}\to\mathbb{R}$, and $\bm{y}=[y_{1},y_{2},\ldots,y_{n}]^{T}\in\mathbb{R}^{n}$.

4 The SGIPSM

In this section, we present the SGIPSM to solve numerically Problems 1 and 2. The derivation of the proposed SGIPSM is divided into two cases according to whether $\beta\neq 0$ or $\beta=0$; therefore, we refer to these two cases by Assumptions 1 and 2, respectively.

We commence our numerical scheme by recasting Problem 1 into its integral formulation. By introducing the useful substitution

[TABLE]

for some unknown function $\phi(x)$, we can recover the solution function and its derivative in terms of $\phi(x)$ by recursive integration; in particular

[TABLE]

Substituting Eqs. (4.1), (4.2), and (4.3) into Eq. (1.1) yields the integral nonlinear Lane-Emden equation as follows:

[TABLE]

If Assumption 1 holds, then we can approximate $y(0)$ using Eqs. (1.3b), (4.2), and (4.3) as follows:

[TABLE]

where $\bm{\Phi}=\phi\left(\bm{x}_{b,n}^{(\alpha)}\right)$, and ${\bm{x}_{b,n}^{(\alpha)}}=\left(x_{b,n,0}^{(\alpha)},x_{b,n,1}^{(\alpha)},\ldots,x_{b,n,n}^{(\alpha)}\right)^{T}$. Let $\mathbf{D}_{\bm{v}}=\mathop{\mathrm{diag}}\left(v_{0},\ldots,v_{n}\right)$, for any vector $\bm{v}=\left(v_{i}\right)\in\mathbb{R}^{n+1},\;i=0,\ldots,n$, and denote the vector of reciprocals $\left(1/v_{0},\ldots,1/v_{n}\right)^{T}$ by $\bm{v}^{\div}$. Moreover, let $\bm{1}_{n}$ and $\mathbf{I}_{n}$ denote the all ones column vector and the identity matrix each of size $n$, respectively. Then collocating the integral nonlinear Lane-Emden Eq. (4.4) at the SFGGR nodes yields the following (n+1)th-order nonlinear system

[TABLE]

and “$\otimes$” denotes the Kronecker product. The solution of the nonlinear system (4.6a) produces the approximate solution vector $\bm{\Phi}$. The original solution vector, $y\left(\bm{x_{b,n}^{(\alpha)}}\right)$, can then be recovered using Eq. (4.3) as follows:

[TABLE]

We can find the approximate solution at any point $x\in(0,b]$ using interpolation at the SFGGR nodes as follows:

[TABLE]

where $y_{b,n,k}^{(\alpha)}=y\left(x_{b,n,k}^{(\alpha)}\right)$, for all $k$,

[TABLE]

is the $n$th-degree shifted Lagrange interpolating polynomial, $\lambda_{b,j}^{(\alpha)}={\left({b/2}\right)^{2\,\alpha}}\lambda_{j}^{(\alpha)}$, and $\varpi_{b,n,j}^{(\alpha)}={\left({b/2}\right)^{2\,\alpha}}\varpi_{n,j}^{(\alpha)}$, for all $j$. This technique can be carried out using Algorithm A.1 in A. For Problem 2, Eqs. (4.6a) can be reduced into the following linear algebraic system

[TABLE]

The original solution vector $y\left(\bm{x_{b,n}^{(\alpha)}}\right)$ can again be recovered using Eq. (4.7). This technique can be implemented using Algorithm A.2 in A.

Now consider the second scenario when Assumption 2 holds. In this case, we can readily impose the boundary condition (1.3b) via Eq. (4.2) by the following approximation:

[TABLE]

Collocating Eq. (4.4) at the SFGGR nodes and combining the result with Eq. (4.11) gives the following $(n+2)$th-order nonlinear system in the $(n+2)$ unknowns $\bm{\Phi}$ and $y(0)$:

[TABLE]

and “$[.;.]$” is the vertical matrix concatenation along columns defined by “${[{.^{T}},{.^{T}}]^{T}}$.” The solution of the nonlinear system (4.12a) produces the approximate solution vector $\bm{\Psi}=\left[\bm{\Phi};y(0)\right]$ of Problem 1. The original solution vector, $y\left(\bm{x_{b,n}^{(\alpha)}}\right)$, can then be recovered using Eq. (4.3) as follows:

[TABLE]

This technique can be achieved using Algorithm A.3 in A. For Problem 2, we can easily show after some mathematical manipulation that Eqs. (4.12a) can be reduced into the following linear algebraic system

[TABLE]

The original solution vector, $y\left(\bm{x_{b,n}^{(\alpha)}}\right)$, can again be recovered using Eq. (4.13). This technique can be accomplished using Algorithm A.4 in A.

5 Error bounds and convergence analysis

We begin this section by stating one technical lemma and three theorems that are needed for the subsequent derivation of the error estimates of the SGIPSM. Let $\mathbb{R}_{0}^{+}=\mathbb{R}^{+}\cup\{0\}$, and $\|f\|_{\infty,\mathbb{S}}=\sup\left\{\,\left|f(x)\right|:x\in\mathbb{S}\,\right\}$, be the uniform norm (or sup norm) of a real-valued bounded function $f$ defined on a set $\mathbb{S}\subset\mathbb{R}$. The following lemma is a modified version of [Elgindy (2016), Lemma 4.1].

Lemma 5.1**.**

The uniform norm of the $n$th-degree Gegenbauer polynomial, ${G_{n}^{(\alpha)}}$, and its associated shifted form, ${G_{b,n}^{(\alpha)}}$, is given by

[TABLE]

where $D^{(\alpha)}>1$ is a constant dependent on $\alpha$, but independent of $n$.

The following three theorems highlight the error bounds of the Gegenbauer quadratures associated with $\mathbf{Q}^{(1)},\mathbf{Q}_{b}^{(1)}$, and $\mathbf{Q}_{b}^{(2)}$, respectively.

Theorem 5.1**.**

Assume that $f(x)\in C^{n+1}[-1,1]$ and $\left\|f^{(n+1)}\right\|_{\infty,[-1,1]}=A\in\mathbb{R}^{+}$, for some constant $A$ independent of $n$. Moreover, let $\int_{-1}^{x_{n,i}^{(\alpha)}}f(x)\,dx$ be approximated by the Gegenbauer quadrature, $\sum\nolimits_{k=0}^{n}{Q_{i,k}^{(1)}\,f_{n,k}^{(\alpha)}},\;i=0,\ldots,n.$ Then there exist some positive constants ${B^{(\alpha)}}$ and ${C^{(\alpha)}}$ dependent on $\alpha$ and independent of n such that the quadrature truncation error is bounded by the following inequalities

[TABLE]

where

[TABLE]

Moreover, when $n\rightarrow\infty$, we have

[TABLE]

for $\alpha\geq 0$, and

[TABLE]

for $-1/2<\alpha<0$, where $\mathop{<}\limits_{\sim}$ means “less than or asymptotically equal to.”

Proof.

Inequality (5.1) is determined through Eqs. (2.2) and (3.5), and Lemma 5.1. For  $n\rightarrow\infty$, we find through [Elgindy (2016), Lemma 4.2] that

[TABLE]

where $B^{(\alpha)}$ is a constant depend on $\alpha$ and independent of $n$. Inequality (5.3) is established by combining this result with Lemma 5.1 such that

[TABLE]

where $D_{1}^{(\alpha)},D_{2}^{(\alpha)}>1$ are constants dependent on $\alpha$ and independent of $n$. The asymptotic inequality (5.4) is derived by realizing that $n+1\sim n$ as $n\to\infty$, and setting $C^{(\alpha)}=B^{(\alpha)}\left(D_{1}^{(\alpha)}+D_{2}^{(\alpha)}\right)$. ∎

Theorem 5.2**.**

Assume that $f(x)\in C^{n+1}[0,b]$ and $\left\|f^{(n+1)}\right\|_{\infty,[0,b]}=A\in\mathbb{R}^{+}$, for some constant $A$ independent of $n$. Moreover, let $\int_{0}^{x_{b,n,i}^{(\alpha)}}f(x)\,dx$ be approximated by the shifted Gegenbauer quadrature, $\sum\limits_{k=0}^{n}{Q_{b,i,k}^{(1)}{\kern 1.0pt}f_{b,n,k}^{(\alpha)}},\;i=0,\ldots,n$. Then there exist some positive constants $B^{(\alpha)}$ and $C^{(\alpha)}$ dependent on $\alpha$ and independent of n such that the quadrature truncation error is bounded by the following inequalities

[TABLE]

Moreover, when $n\rightarrow\infty$, we have

[TABLE]

for $\alpha\geq 0$, and

[TABLE]

for $-1/2<\alpha<0$.

Proof.

The proof is straightforward using Eqs. (2.2) and (3.18), Lemma 5.1, and [Elgindy (2016), Lemma 4.2]. ∎

Theorem 5.3**.**

Assume that $f(x)\in C^{n+1}[0,b]$ and ${\left\|{{f^{(n+k)}}}\right\|_{\infty,[0,b]}}={A_{k}}\in{\mathbb{R}^{+}},k=0,1,$ where the constants ${A_{0}}$ and ${A_{1}}$ are independent of $n$. Moreover, let $\int_{0}^{x_{b,n,i}^{(\alpha)}}\int_{0}^{x}f\left(t\right)\,dt\,dx$ be approximated by the shifted Gegenbauer quadrature, $\sum\limits_{k=0}^{n}{Q_{b,i,k}^{(2)}\,f_{b,n,k}^{(\alpha)}},\;i=0,\ldots,n.$ Then there exist some positive constants $C^{(\alpha)}$ and $D^{(\alpha)}$ dependent on $\alpha$ and independent of $n$ such that the quadrature truncation error is bounded by the following inequalities

[TABLE]

Moreover, when $n\rightarrow\infty$, we have

[TABLE]

for $\alpha\geq 0$, and

[TABLE]

for $-1/2<\alpha<0$.

Proof.

The proof is straightforward using Eqs. (2.2) and (3.22), Lemma 5.1, and [Elgindy (2016), Lemma 4.2]. ∎

Now suppose that we denote the exact and approximate solutions of each of Problems 1 and 2 by $y(x)$ and $y_{\text{app},n}(x)$, respectively. We are now ready to give the a priori and asymptotic a priori error estimates through the following key theorem, which proves the elegant exponential convergence rate of the SGIPSM for Problem 1.

Theorem 5.4** (A priori and asymptotic a priori error estimates).**

Assume that $\left\|y^{(n+k+2)}\right\|_{\infty,(0,b]}=A_{k}\in\mathbb{R}^{+},\;k=0,1$, and $y$ is approximated by the nodal truncated series expansion of shifted Gegenbauer basis polynomials (4.8), where the second-order derivative values $y^{\prime\prime}\left({x_{b,n,i}^{(\alpha)}}\right),i=0,\ldots,n,$ are obtained by collocating the integral nonlinear Lane-Emden equation (4.4) at the SFGGR points. Suppose also that the boundary conditions (1.3a) and (1.3b) and Assumption 1 hold. Then the solution error values, $y_{b,n,i}^{(\alpha)}-{y_{\text{app},n}}\left({x_{b,n,i}^{(\alpha)}}\right),\;i=0,\ldots,n$, are bounded by the following inequalities

[TABLE]

Moreover, when $n\rightarrow\infty$, the asymptotic solution error values are bounded by

[TABLE]

for $\alpha\geq 0$, and

[TABLE]

for $\frac{-1}{2}<\alpha<0$, where $\mu_{1,b}^{(\alpha)}$ and $\mu_{2,b}^{(\alpha)}$ are positive constants dependent on $b$ and $\alpha$, and independent of $n$.

Proof.

Using Eqs. (4.3) and (4.5), and Theorems 5.2 and 5.3, we can easily show that

[TABLE]

where ${}_{y^{\prime\prime}}E_{b,n,i}^{(j),\alpha}={}_{y^{\prime\prime}}E_{b,n}^{(j)}\left(x_{b,n,i}^{(\alpha)},\xi_{b,n,i}^{(j)}\right)$, for some $\xi_{b,n,i}^{(j)}\in(0,b),i=0,\ldots,n;j=1,2$. This completes the proof of the theorem. ∎

Theorem 5.4 shows that

[TABLE]

with an error upper bound given by

[TABLE]

The following theorem marks the priori and the asymptotic priori error estimates of the solution derivative when Assumption 2 holds.

Theorem 5.5** (A priori and asymptotic a priori error estimates).**

Assume that $\left\|y^{(n+3)}\right\|_{\infty,(0,b]}=A\in\mathbb{R}^{+}$, and $y$ is approximated by the nodal truncated series expansion of shifted Gegenbauer basis polynomials (4.8), where the second-order derivative values $y^{\prime\prime}\left({x_{b,n,i}^{(\alpha)}}\right),i=0,\ldots,n,$ are obtained by collocating the integral nonlinear Lane-Emden equation (4.4) at the SFGGR points. Suppose also that the boundary conditions (1.3a) and (1.3b) and Assumption 2 hold. Then

[TABLE]

and the solution derivative error values, $y^{\prime(\alpha)}_{b,n,i}-{y^{\prime}_{\text{app},n}}\left({x_{b,n,i}^{(\alpha)}}\right),\;i=1,\ldots,n$, are bounded by the following inequalities

[TABLE]

where $y^{\prime(\alpha)}_{b,n,k}=y^{\prime}\left(x_{b,n,k}^{(\alpha)}\right)$, for all $k$. Moreover, there exist some positive constants $B^{(\alpha)}$ and $C^{(\alpha)}$ dependent on $\alpha$ and independent of n such that the asymptotic solution derivative error values are bounded by

[TABLE]

when $n\rightarrow\infty$.

Proof.

Eqs. (5.18) result by imposing the boundary conditions (1.3a) and (1.3b). Inequalities (5.19), (5.20), and (5.21) are established using Eq. (4.2) and Theorem 5.2. ∎

Theorem 5.5 shows that

[TABLE]

with an error upper bound given by

[TABLE]

The following corollary highlights the uniform convergence of the solution function approximations $y_{\text{app},n},n=1,2,\ldots$ under Assumption 2.

Corollary 5.1** (Uniform convergence).**

Assume that $\left\|y^{(n+3)}\right\|_{\infty,(0,b]}<\infty$, and $y$ is approximated by the nodal truncated series expansion of shifted Gegenbauer basis polynomials (4.8), where the second-order derivative values $y^{\prime\prime}\left({x_{b,n,i}^{(\alpha)}}\right),i=0,\ldots,n,$ are obtained by collocating the integral nonlinear Lane-Emden equation (4.4) at the SFGGR points. Suppose also that the boundary conditions (1.3a) and (1.3b) and Assumption 2 hold. Then the sequence of solution function approximations $\{y_{\text{app},n}\}_{n=1}^{\infty}$ converges uniformly on $(0,b]$ to the solution function $y$ with exponential rate.

Proof.

Eq. (5.18b) implies that ${\lim_{n\to\infty}}{y_{\text{app},n}}(b)$ exists. We can also infer from inequalities (5.20) and (5.21) that there exists a positive integer $N$ such that

[TABLE]

for all $n\geq N$ and $x\in(0,b]$. Hence, the SGIPSM generates a sequence of derivative approximations $\{y^{\prime}_{\text{app},n}\}_{n=1}^{\infty}$ that converges uniformly to the solution derivative function $y^{\prime}$, for all $x\in(0,b]$ with exponential rate. This completes the required proof. ∎

Now, let us define the residual function $\mathcal{R}$ by

[TABLE]

Using Eqs. (4.1)-(4.3), we have

[TABLE]

Hence, the residual function can be estimated at the SFGGR points by

[TABLE]

We can further estimate the residual error bounds of the SGIPSM for Problem 1 at the SFGGR points through the following two useful theorems.

Theorem 5.6**.**

Assume that $\left\|y^{(n+k+2)}\right\|_{\infty,(0,b]}=A_{k}\in\mathbb{R}^{+},\;k=0,1$, and suppose also that the function $f(x,y)$ satisfies Lipschitz condition in the variable $y$ on the set $\mathbb{S}_{b,n}^{(\alpha)}\times\mathbb{R}$ with Lipschitz constant $\lambda>0$. If the boundary conditions (1.3a) and (1.3b) and Assumption 1 hold, then the SGIPSM discretizes the integral nonlinear Lane-Emden equation (4.4) at the SFGGR nodes $x_{b,n,i}^{(\alpha)}\in\mathbb{S}_{b,n}^{(\alpha)},\;i=0,\ldots,n$, such that the residual function values $\mathcal{R}\left(x_{b,n,i}^{(\alpha)}\right),\;i=0,\ldots,n$ are bounded by the inequalities

[TABLE]

Moreover, when $n\rightarrow\infty$, the asymptotic residual function values are bounded by

[TABLE]

for $\alpha\geq 0$, and

[TABLE]

for $\frac{-1}{2}<\alpha<0$, where $\mathcal{L}_{1,\lambda,b}^{(\alpha)}$ and $\mathcal{L}_{2,\lambda,b}^{(\alpha)}$ are positive constants dependent on $\lambda,b,\alpha$, and independent of $n$.

Proof.

Subtracting Eq. (1.1) from Eq. (5.26), and using Lipschitz condition together with Eqs. (1.3b), (4.3), (4.5), Theorems 5.2 and 5.3 yield

[TABLE]

from which the proof is established. ∎

Theorem 5.7**.**

Assume that $\left\|y^{(n+3)}\right\|_{\infty,(0,b]}<\infty$, and suppose also that the function $f(x,y)$ satisfies Lipschitz condition in the variable $y$ on the set $\mathbb{S}_{b,n}^{(\alpha)}\times\mathbb{R}$ with Lipschitz constant $\lambda>0$. If the boundary conditions (1.3a) and (1.3b) and Assumption 2 hold, then there exists a relatively small positive number $0<\varepsilon\ll 1$ and a positive integer $N$ such that the SGIPSM discretizes the integral nonlinear Lane-Emden equation (4.4) at the SFGGR nodes $x_{b,n,i}^{(\alpha)}\in\mathbb{S}_{b,n}^{(\alpha)},\;i=0,\ldots,n$ with asymptotic residual function values $\mathcal{R}\left(x_{b,n,i}^{(\alpha)}\right),\;i=0,\ldots,n$ bounded by the inequalities

[TABLE]

for all $n\geq N$, where ${B^{(\alpha)}}$ and ${C^{(\alpha)}}$ are positive constants dependent on $\alpha$, and independent of $n$.

Proof.

By the uniform convergence of Corollary 5.1, we find that

[TABLE]

The proof is established by combining this result with inequality (5.34). ∎

Similarly, to furnish the exponential convergence rate of the SGIPSM for Problem 2, notice that the residual term in this case takes the following special form:

[TABLE]

where the differential operator $L=\displaystyle{\frac{d^{2}}{dx^{2}}+\frac{\alpha_{2}}{x}\frac{d}{dx}}+p(x)$. The following two corollaries are direct results of Theorems 5.4 and 5.7.

Corollary 5.2**.**

Assume that $\left\|y^{(n+k+2)}\right\|_{\infty,(0,b]}=A_{k}\in\mathbb{R}^{+},\;k=0,1$, and $\left\|p\right\|_{\infty,(0,b]}=M\in\mathbb{R}_{0}^{+}$. If the boundary conditions (1.3a) and (1.3b) and Assumption 1 hold, then the SGIPSM discretizes the integral nonlinear Lane-Emden equation (4.4) at the SFGGR nodes $x_{b,n,i}^{(\alpha)}\in\mathbb{S}_{b,n}^{(\alpha)},\;i=0,\ldots,n$, such that the residual function values $\mathcal{R}\left(x_{b,n,i}^{(\alpha)}\right),\;i=0,\ldots,n$ are bounded by the inequalities

[TABLE]

Moreover, when $n\rightarrow\infty$, the asymptotic residual function values are bounded by

[TABLE]

for $\alpha\geq 0$, and

[TABLE]

for $\frac{-1}{2}<\alpha<0$, where $\mathcal{L}_{1,M,b}^{(\alpha)}$ and $\mathcal{L}_{2,M,b}^{(\alpha)}$ are positive constants dependent on $M,b,\alpha$, and independent of $n$.

Corollary 5.3**.**

Assume that $\left\|p\right\|_{\infty,(0,b]}=M\in\mathbb{R}_{0}^{+}$. If the boundary conditions (1.3a) and (1.3b) and Assumption 2 hold, then there exists a relatively small positive number $0<\varepsilon\ll 1$ and a positive integer $N$ such that the SGIPSM discretizes the integral nonlinear Lane-Emden equation (4.4) at the SFGGR nodes $x_{b,n,i}^{(\alpha)}\in\mathbb{S}_{b,n}^{(\alpha)},\;i=0,\ldots,n$ with asymptotic residual function values $\mathcal{R}\left(x_{b,n,i}^{(\alpha)}\right),\;i=0,\ldots,n$ bounded by the inequalities

[TABLE]

for all $n\geq N$, where ${B^{(\alpha)}}$ and ${C^{(\alpha)}}$ are positive constants dependent on $\alpha$, and independent of $n$.

Proof.

The proof is similar to that of Theorem 5.7. ∎

6 Numerical experiments

In this section, we apply the proposed SGIPSM on five test examples of Problems 1 and 2. All numerical experiments were carried out using MATLAB R2011b (7.13.0.564) software installed on a personal laptop equipped with an Intel(R) Core(TM) i3-4500U CPU with 1.70 GHz speed running on Windows 7 Ultimate 32-bit operating system. The exact solutions were calculated using MATHEMATICA 7 with $18$ digits of precision maintained in internal computations. In all numerical tests, by AE and RE, we mean the absolute and relative errors, respectively. MAE and AEb denote the maximum AE and the AE in the computation of the solution function $y$ at $x=b$, respectively. Furthermore, by $\kappa_{\infty}$, we mean the infinity norm condition number of the resulting linear algebraic system associated with Problem 2. The resulting linear algebraic systems of equations were solved using MATLAB “mldivide” algorithm. Nonlinear algebraic systems associated with Problem 1 were solved by MATLAB “fsolve” solver with termination tolerance “TolFun” set at $10^{-15}$.

Example 1. Consider Problem 2 with $\alpha_{1}=0,\alpha_{2}=1,\beta=1,\gamma=0,\delta=0,b=1,p(x)=0$, and $g(x)=\left(8/\left(8-x^{2}\right)\right)^{2}$. The exact solution is $y(x)=2\log\left(7/\left(8-x^{2}\right)\right)$. Figure 1 shows the plots of the exact solution, approximate solution, AE, and MAE on $[0,1]$; in addition, the figure shows $\kappa_{\infty}$ obtained by the SGIPSM for several values of $n$ and $\alpha$.

Example 1 was previously solved by [Caglar and Caglar (2006)] and later by [Yüzbaşı and Sezer (2013)] using a B-spline method and an improved Bessel’s method, respectively. Table 1 shows the power of the state-of-the-art SGIPSM when compared with the B-spline method of [Caglar and Caglar (2006)], as illustrated by gaining extra digits of precision using far lesser number of collocation points. The SGIPSM outperforms also the method of [Yüzbaşı and Sezer (2013)] as clearly observed from Table 2, where higher-order approximations are achieved using the same number of collocation points. Moreover, the SGIPSM performs a single discretization to the integral form of the problem against two discretizations for the method of [Yüzbaşı and Sezer (2013)] applied to the problem and its associated error problem; thus, significantly reducing the necessary computational cost to establish the same level of accuracy. In addition, thanks to the discretization process employing the SFGGR collocation points, the boundary condition, $y(1)=0$, is satisfied exactly by the SGIPSM.

Example 2. Consider Problem 1 with $\alpha_{1}=0,\alpha_{2}=2,\beta=1,\gamma=0,\delta=\sqrt{3}/2,b=1$, and $f(x,y)=y^{5}$. The exact solution is $y(x)=1/\sqrt{1+x^{2}/{3}}$. Figure 2 shows the plots of the exact solution, approximate solution, AE, and MAE on $[0,1]$ for several values of $n$ and $\alpha$. [Kanth (2007)] solved the problem using a cubic spline method combined with a quesilinearization technique that reduces the nonlinear problem into a sequence of linear problems. [Turkyilmazoglu (2013)] later approached the solution of the problem using a continuous functions series expansion technique in which the problem is discretized into a system of nonlinear algebraic equations by expanding the solution in a truncated Maclaurin series. The expansion coefficients are then calculated through a Galerkin-like method. Table 3 shows the rapid convergence rates and the memory minimizing feature of the SGIPSM against the cubic spline method of [Kanth (2007)], where the former method produces higher accuracy using relatively much smaller number of collocation points. On the other hand, the latter method requires the solution of a nonlinear system of equations of order $53$, versus only $8$ for the SGIPSM to finally yield approximations of lower accuracy. The SGIPSM gives also greater accuracy than that of [Turkyilmazoglu (2013)] for several values of $n$ as clearly observed from Table 4. The small approximation error at $x=1$ is due to the inexactness in the computer representation of the real number $\sqrt{3}/2$.

Example 3. Consider Problem 2 with $\alpha_{1}=0,\alpha_{2}=2,\beta=1,\gamma=0,\delta=5.5,b=1,p(x)=-4$, and $g(x)=-2$. The exact solution is $y(x)=0.5+5\sinh(2x)/(x\sinh(2))$. Figure 3 shows the plots of the exact solution, approximate solution, AE, MAE, and $\kappa_{\infty}$ obtained by the SGIPSM for several values of $n$ and $\alpha$. [Kanth and Reddy (2005)] solved the problem by modifying the original differential equation (1.2) at the singular point $x=0$ via L’Hôpital’s rule, and then collocating the modified differential equation at equispaced grid points using a cubic spline method. Such a numerical scheme discretized the problem into a tridiagonal system of equations that was solved by Thomas algorithm. Table 5 shows the preference of the proposed SGIPSM against the method of [Kanth and Reddy (2005)] in terms of higher accuracy and reduced computational cost.

Example 4. Consider Problem 1 with $\alpha_{1}=0,\alpha_{2}=1,\beta=1,\gamma=0,\delta=2\ln\left(\left(4-2\sqrt{2}\right)/\left(7.75-4.5\sqrt{2}\right)\right),b=1.5$, and $f(x,y)=e^{y}$. The exact solution is $y(x)=2\ln\left((c+1)/\left(cx^{2}+1\right)\right)$ with $c=3-2\sqrt{2}$. Figure 4 shows the plots of the exact solution, approximate solution, and AE on $[0,1.5]$ using $n=5$ and $\alpha=0.9$, in addition to the MAE for $n=2^{k},k=2(1)7$ and $\alpha=-0.4(0.1)2$. The REs of the SGIPSM using $n=5$ and $\alpha=0.9$ is shown in Table 6. This problem was previously solved by [Khuri and Sayfy (2010)] using a decomposition method in combination with the cubic B-spline collocation technique. In particular, the former method is employed in the vicinity of the singularity $x=0$, where the nonlinear term $f(x,y)$ is decomposed in terms of the Adomian polynomials, while the latter method is implemented outside this range. Table 7 shows however higher-order approximations in favor of the SGIPSM using several values of $n$.

Example 5. Consider Problem 1 with $\alpha_{1}=-1,\alpha_{2}=2,\beta=0,\gamma=1,\delta=-1,b=1$, and $f(x,y)=\sin(y)-\cos(x)+2/x$. The exact solution is $y(x)=\pi/2-x$. Figure 5 shows the plots of the exact solution, approximate solution, AE, and MAE on $[0,1]$ for several values of $n$ and $\alpha$. Table 8 verifies the strength of the proposed method as shown by the exceedingly accurate approximations obtained by the method using as small as $8$ collocation points.

7 Conclusion and discussion

In this work, we introduced the SGIPSM: a novel, highly accurate, exponentially convergent, and efficient numerical method for solving singular linear and nonlinear Lane-Emden equations provided with certain mixed Neumann and Robin boundary conditions. The proposed method exploits the well conditioning of numerical integral operators by working on the integral formulation of the problem. The reduced problem is then collocated at the SFGGR nodes, and the definite integrals are approximated using some novel SFGGR-based shifted Gegenbauer integration matrices. Such numerical integration operators are constant operators that can be stored for certain sets of SFGGR points, and invoked later when implementing the proposed computational algorithms. The numerical scheme eventually discretizes the problem into a system of linear/nonlinear equations that can be solved using standard linear/nonlinear system solvers. The exponential convergence of the SGIPSM is verified theoretically through Theorem 5.4 and Corollary 5.1. Moreover, estimates of the residual error bounds were provided by Theorems 5.6 and 5.7 and Corollaries 5.2 and 5.3. Five test examples were solved to assess the power of the SGIPSM. In all test examples, the exponential convergence of the SGIPSM is clearly demonstrated, where higher-order approximations are achieved using relatively small numbers of collocation nodes. Moreover, the MAE of the method generally approaches the machine epsilon at either $n=16$ or $n=32$, and then degrades slightly for increasing values of $n$ due to the presence of round-off errors. Adopting the SFGGR nodes for collocating Problems 1 and 2 permits for approximating the solution function $y$ at $x=b$ very accurately. The numerical simulations of Examples 1 and 3 manifest that $\kappa_{\infty}$ increases monotonically for increasing values of $n$ with a swelling rate of increase for increasing positive values of $\alpha$. For instance, the value of $\kappa_{\infty}$ of Example 1 jumps from about $4$ at $\alpha=-0.4$ and $n=128$ into more than $10^{5}$ at $\alpha=2$ for the same value of $n$. A similar pattern is also observed in Example 3, where the value of $\kappa_{\infty}$ surges from about $10$ at $\alpha=-0.4$ and $n=128$ into more than $10^{6}$ at $\alpha=2$ for the same value of $n$. Therefore, for very large values of $n$, we would expect a significant drop in the precision of the approximate solutions for increasing positive values of $\alpha$– a result that is reflected by the theoretical aftermath of Eqs. (5.17) and (5.23). We draw the attention of the reader also to the significant observation pointed out by [Elgindy and Smith-Miles (2013)] that the Gegenbauer quadratures employed in the discretization process ‘may become sensitive to round-off errors for positive and large values of the parameter $\alpha$ due to the narrowing effect of the Gegenbauer weight function,’ which drives the quadratures to become more extrapolatory and possibly give rise to poor integral approximations. On the other hand, we can infer from both examples that the resulting linear algebraic systems remain well conditioned as $n$ tends to infinity for $\alpha\in[-0.4,0]$ with a general monotonic increase of $\kappa_{\infty}$ for increasing values of $n$ and $\alpha$. That is, Gegenbauer basis polynomials associated with negative values of $\alpha$ generally lead to more well conditioned linear systems than their rival Chebyshev and Legendre bases polynomials as $n$ tends to infinity. We observe though that $\kappa_{\infty}$ remains plausible for $\alpha\in[-0.4,1]$ in the range $n\leq 16$. In fact, if we denote the interval $[-1/2+\varepsilon,r]$ by ${I_{\varepsilon,r}^{G}}$, for some relatively small positive number $\varepsilon$ and $r\in[1,2]$, then we closely follow the rule of thumb recommended by [Elgindy (2016)] in the sense of generally performing collocations for values of $\alpha\in{I_{\varepsilon,r}^{G}}$, for small/medium numbers of collocation points and shifted Gegenbauer expansion terms; however, collocations at the shifted flipped-Chebyshev-Gauss-Radau points should be put into effect for large numbers of collocation points and shifted Gegenbauer expansion terms if the approximations are sought in the infinity norm (Chebyshev norm). Numerical comparisons with other traditional methods shown in Tables 1–5, and 7 prove the preference of the proposed method over traditional methods in the literature with respect to the accuracy of the calculated approximations, convergence rate, and the employed number of collocation points. The presented SGIPSM is easily programmed, and can be extended to solve various problems in many areas of science and applications.

Appendix A Computational algorithms