A new numerical strategy for solving nonlinear singular Emden-Fowler delay differential models with variable order

The present study is related to the numerical solutions of new mathematical models based on the variable order Emden-Fowler delay differential equations. The shifted fractional Gegenbauer, CS,j(α,μ)(t),\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C_{S,j}^{(\alpha ,\mu )}(t),$$\end{document} operational matrices (OMs) of VO differentiation, in conjunction with the spectral collocation method are used to solve aforementioned models numerically. The VO operator of differentiation will be used in the Caputo sense. The proposed technique simplifies these models by reducing them to systems of algebraic equations that are easy to solve. To determine the effectiveness and accuracy of the sugested technique, the absolute errors and maximum absolute errors for four realistic models are studied and illustrated by several tables and graphs at different values of the VO and the SFG parameters; α\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha$$\end{document} and μ.\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu.$$\end{document} Also numerical comparisons between the suggested technique with other numerical methods in the existing literature are held. The numerical results confirm that the suggested technique is accurate, computationally efficient and easy to implement.


Introduction
In recent years, many researchers have been interested in the study of initial value problems of some special second-order nonlinear singular ordinary differential equations (ODEs).One of these most important equations is the Emden-Fowler (EF) equation, which has a variety of applications in various fields.Robert Emden [1] and Fowler [2] have studied the EF equation in detail.The general mathematical model of EF is expressed as with the initial conditions (ICs) where g(t) and f(y(t)) are functions of t and y(t), respectively, r, , 1 , and 2 are constants.Solving differential equations with singularity behavior is a defy.In particular, the current EF equation problem, which involves a singularity at t = 0 , is crucial in practical applications.Analytically, solv- ing this equation is difficult, so to solve this problem various numerical techniques were used [3,4].For g(t) = 1, and = 1, Eq. (1.1) moderates to the Lane-Emden (LE) equation in its usual form.This equation has appeared in the work of Homer Lane and Robert Emden work.The inner construction of polytrophic stars, the gas cloud model, cluster galaxies, and radiative cooling are all represented by these models.No one can deny the usefulness and importance of these models, which have a wide range of applications in mathematical physics, astrophysics, and celestial mechanics such as isotropic continuous media, isothermal gas sphere models [5], dusty fluid models [6], oscillating magnetic systems [7], catalytic diffusion reactions [8], stellar structure (1.1) y �� (t) + r t y � (t) + g(t)f (y(t)) = 0, r ≥ 1, models [9], electromagnetic theory [10], gaseous star density [11], and classical quantum mechanics [12].
The nonlinear singular ODE (1.1) was first developed by Jonathan Homer Lane [13], a U.S. astrophysicist who was interested in estimating both the temperature and the density of mass on the surface and after 37 years was investigated in greater depth by Emden [1].This ODE in astrophysics is a dimensionless version of Poisson's equation for a simple stellar model's gravitational potential [14].Go back to [15,16] to get a sufficient explanation of deriving of the classic LE equation for modeling the thermal dynamics of a spherical cloud, which takes the following form: with the ICs In quantum mechanics and astrophysics, the values of m physically lie in the interval [0,5].In different literatures, the exact solutions of problem (1.3-1.4) have only been found at m = 0, 1, 5 [16], and for other values of m, several methods have been successfully used to approximate the solution of this problem [17,18].
One of the most essential significant classes of the differential equations is the delay differential (DD) equations.Due to their contributions significantly to solve the real-life problems.These equations have numerous applications in many biological models, as well as scientific phenomena such as communication system models, dynamical population models, economic systems, engineering systems, and transport models, so they have attracted a lot of attention from the research community [19,20], and various numerical/analytical techniques are introduced to solve them [21,22].
Fractional calculus is a new topic that appears as a generalization of integer calculus, in which the order will be non-integer (arbitrary).The dynamic behavior of various phenomena in engineering sciences, physics, and other disciplines of research can be explored more precisely using noninteger order differential equations.So, there are a growing number of applications of fractional-order calculus in different fields as e.g.Long electrical lines, electrochemical processes, dielectric polarization, colored noise, viscoelastic materials, chaos, control theory, and in many other areas.Recently, the VO differential equation in which the order is a function of space and/or time have been demonstrated their ability to accurately describe real world phenomena [23].In the VO model, the next state is influenced by all of its former states, not only upon its current state as in the case of integer-order derivatives.So, the VO derivative is more realistic and generalized than the standard and fractional ones.VO differential equations have gotten a lot of attention lately as a result of their suitability modeling a wide range of phenomena across many fields of science and engineering, such as anomalous diffusion [24], viscoelastic mechanics [25], control systems [26], petroleum engineering [27], and many other branches of physics and engineering [28].
With this respect, the importance of the VO-EF-DDEs are quite clear, due to their massive applications in mathematical physics and engineering problems.For this significance, obtaining the numerical solutions of this type of equations are required since the exact solutions are usually not easy to claim, because of the kernel of the VO operators have a changing exponent.So, the main objective of the present work is to develop an efficient and accurate numerical technique for solving the following new classes of VO-EF-DDEs: • A class of variable order Emden-Fowler delay differential equations (VO-EF-DDEs) of the following form: with the ICs where ≥ 1, (t), g(t), and F(t) are given functions, the parameters a, b, c, d, e, f , 1 , and 2 are known constants also D (t) t denotes the Caputo fractional derivative of variable order, 1 < (t) ≤ 2 and y(t) is an unknown function.By taking g(t) = 1 the EF model reduces to the Lane- Emden model.The description of problem (1.5) in the VO operator (t) has many advantages one of them; it describes the behavior of the system and its properties in a more accurate manner.Also, it's important to mention that the model (1.5-1.6)covers the classical and the constant fractional models when (t) is constant.
• A system of variable order Emden-Fowler delay of differential equations (VO-EF-DDEs) of the following form with the following ICs where 1 , 2 ≥ 1, a, b, c, d, e, f , 1 , 2 , 3 and 4 are known constants, the functions (t), g 1 (t), g 2 (t), F 1 (t) and F 2 (t) are assumed to be known functions, and also D (t) t is the Caputo fractional derivative of order 1 < (x) ≤ 2.
(1.5) D (t)  t y(at (1.7) ) By taking g 1 (t) = g 2 (t) = 1 , the EF model reduces to the Lane-Emden model.The suggested technique is based on the spectral collocation method with the OMs of variable/integer order differentiation of the SFGPs.To the best of our knowledge, the numerical treatments of VO-EF-DDEs and the system of VO-EF-DDEs have not been established by using SFGPs yet.
Spectral collocation methods are known to be one of the powerful numerical techniques for solving different kinds of differential equations due to their applicability in bounded and unbounded domains.The convergence speed is one of the major advantages of the spectral method.In general, spectral methods are promising candidates for solving fractional differential equations since their global nature fits well with the non-local definition of fractional operators [29,30].Spectral collocation methods have been widely used for solving problems with smooth solutions because of their efficiency, exponential rate of convergence, low computational cost, and high order of accuracy achieved using a minimal number of grid points.It is known that the accuracy of spectral methods increases with an increase in grid points but beyond a certain number of grid points, the accuracy rapidly deteriorates.Also, the limitation of spectral methods is that their accuracy deteriorates for complicated domains.OMs were restructured for solving numerous types of fractional differential equations (FDEs).The procedure of numerical methods in combining with OMs of some orthogonal polynomials for solving FDEs created extremely precise solutions for such equations.The essential step in the proposed technique is to represent the numerical solution of the problem as a finite sum using SFGPs as basis functions, then the OMs corresponding to VO/ integer order derivatives are derived.By using these OM in Eq. (1.5), and applying the collocation method which requires that the residual of VO-EF-DDEs is vanishing at N points that correspond to the Gauss quadrature nodes.When the obtained system of algebraic equations is completed with the initial conditions, (N + 1) algebraic equations are produced, which can be solved by utilizing any appropriate iterative technique.
The SFGPs have a number of characteristics that make them an ideal candidate for our suggested technique, like: • Due to the unique behavior of fractional differential and/ or integral equations, the direct use of spectral methods using traditional orthogonal polynomials such as Legendre, Chebyshev, and Jacobi has low convergence rates.To remedy this difficulty, we employ SFGPs (see [31]), in which the variable t in the SCPs is replaced by t  , 0 <  ≤ 1, t ∈ [0, L].One of the benefits of frac- tional-order polynomials is in dealing with problems with smooth solutions, where any choice of the fractional factor yields a highly accurate approximate solution.The fractional-order polynomials can lessen the loss of the order of convergence of the non-smoothness problems by suitably chosen of the factor μ. • For = 1, the SFGPS reduced to the SGPs.Numerous studies have shown that the Gegenbauer polynomials (GPs) are very effective in solving a wide range of issues [32][33][34].• The parameter  > −0.5 distributes the Gegenbauer poly- nomials (GPs).Every time this parameter was altered, a new polynomial was formed.As a result, they involve a limitless number of orthogonal polynomials like the Chebyshev polynomials (CPs) of the first kind with the parameter = 0, the CPs of the second kind with the parameter = 1, and the shifted Legendre polynomials (LPs) with the parameter = 1 2 .The same results are assumed for the SFGPs.
The main advantages of the proposed scheme are the lowcost computing, small CPU time, and the simple implementation.Also the present method has the ability to convert the given problem into a system of algebraic equations, which is easy to solve.Many numerical experiments are presented to measure the quality of the presented scheme.The obtained numerical results are compared with other methods in the existing literature.
The following is the contents of this paper: The preliminaries of VO-fractional calculus and SFGPs are presented in Sect. 2. The various types of shifted fractional Gegenbauer operational matrices (SFGOMs) are deduced in Sect.3. In Sect.4, the methodology of the proposed technique is explained.In Sect.5, some applications of the proposed algorithm to some numerical problems are offered to demonstrate the effectiveness of the proposed algorithm.At last, a conclusion is given in Sect.6.

Variable-order fractional derivative definition [35]
The Caputo fractional derivative of order (t) for a continu- ously differentiable function y ∶ [0, ∞] → R is defined as where (t) is a positive continuous bounded function and Because of the ICs for fractional differential equations with Caputo derivatives are the same as for integer order differential equations; the Caputo type definition is particularly useful in a wide range of applications.The results stated in definition (2.1) have the following characteristics: (2.1) • Linearity property:

Gegenbauer polynomials
The Gegenbauer polynomials C ( ) j (x), of degree j ∈ Z + , and associated with the parameter  > −1 2 are a sequence of real polynomials in the finite domain [−1, 1].They are a family of orthogonal polynomials which has many applications [32].
• A suitable standardization of the Gegenbauer polynomials C ( ) j (x) of degree j and associated with the parameter  > −1 2 dates back to Doha [34], where the Gegenbauer polynomials can be represented by or equivalently where P ( − 1 2 , − 1 2 ) j (x) is the Jacobi polynomial of degree j.
• As a result of this standardized C (0) j (x) becomes identical with the Chebyshev polynomials (CPs) of the first kind is the Legendre polynomials (LPs) L j (x).C (1)  j (x) is equal to 1 (j+1) U j (x), where U j (x) is the CPs of the second type [34].• The analytical form of the GPs is given by, • The Gegenbauer polynomials can be generated using the following useful recurrence equation • The Gegenbauer polynomials satisfy the following orthogonality relation (2.3) where ( ) (t) is the weight function, and it is an even function given by the relation and is the normalization factor, and i,j is the Kronecker delta function.• We denote the zeros of the Gegenbauer polynomials (Gegenbauer-Gauss nodes) C ( ) N+1 (x) by x ( ) N,j , j = 0, ⋯ , N.

Shifted Gegenbauer polynomials
• For using Gegenbauer polynomials in the interval [0, L], the shifted Gegenbauer polynomials (SGPs) are formed by replacing the variable x with z = 2x L − 1, 0 ≤ z ≤ L. So, we can write SGPs as • The explicit analytical form of the SGPs is given as • These polynomials recover the shifted CPs of the first kind T S,j (z) ≡ C (0) S,j (z), the shifted LPs L S,j (z) ≡ C ( S,j (z), and the shifted CPs of the second kind C (1)  S,j (z) ≡ 1 j+1 U S,j (z).• The orthogonal relation of SGPs is getting from where ( ) S (z) is the weight function, it is an even func- tion given by the relation and • We denote the zeros of the shifted Gegenbauer polynomial (shifted Gegenbauer-Gauss polynomial) C ( ) S,N+1 (z) by z ( ) S,N,j , j = 0, … , N. (2.4) ( ) (2.7) • The q-th derivative of C ( ) S,i (z), is obtained by

T h e o r e m 2 . 1 S u p p o s e t h a t t h e f u n ct i o n
is the best approximation to y(t) from ( , )  S,N , then the error bound is presented as where ‖D (N+1) y(t)‖ ≤ H , t ∈ [0, 1].
Proof Consider the generalized Taylor formula [36] (2.15) (2.17) Then, Since y N (t) is the best square approximation function of y(t), then By taking the square roots, the proof is completed.◻

SFGOMs of VO fractional differentiation
Theorem 3.1 The first order derivative of the SFG vector Φ (at + b) is approximated as where the OM Ψ is a square matrix of order (N + 1) × (N + 1) that takes the following form, where A is defined in (2.16), Υ N and Θ N are (N + 1) × (N + 1) matrices, and their elements ij , and ij , 0 ≤ i, j ≤ N are given, respectively, as follows (3.3) By substituting Eq. (3.5) into Eq.(3.2), the following relation is obtained Since A is invertible, then where Ψ is a square matrix, and this finishes the proof.

Theorem 3.2
The VO Caputo fractional derivative of the SFG vector Φ (at + b) can be approximated as where the OM Ξ is a square matrix of order (N + 1) × (N + 1) takes the following form, where A is defined in (2.16), Υ N and Λ N are (N + 1) × (N + 1) matrices, and their elements ij , and ij , 0 ≤ i, j ≤ N , respectively, are given as follows where Υ N, , Θ N are (N + 1) × (N + 1) matrices, and their ele- ments have the following forms, respectively, , for i = j, j = n, n + 1, … , N, 0, o.w. .

Since
By using Newton ′ s generalized binomial theorem, we can write this series converges for ∈ Z + or | at b | < 1.By approximating the series (3.10) by using the first (N+1) terms, the relation (3.9) can be written as By using the definition (2.1) of the VO fractional derivative in the Caputo sense, then Eq. (3.11) can be obtained as the below form: where Υ N, , Θ N are (N + 1) × (N + 1) matrices, and their ele- ments have the following forms, respectively, (3.8) (3.9) (3.11) , for i = j, j = n, n + 1, … , N, 0, o.w. .
By substituting from Eq. (3.12) into Eq.(3.8), the following relation is obtained where Ξ is a square matrix, and this finishes the proof.◻

The methodology
In this section, we will discuss how to use the SFG-OMs to solve the problem Eqs.(1.5-1.6)numerically.As a result, the problem will be converted into an algebraic system of equations that can be solved easily by any iteration method.By using the properties of the SFGPs which are mentioned in the definition (2.2), the solution can be approximated as by the same way, the terms y(ct + d), y(et + f ) are approximated.Moreover, the Caputo fractional derivative of order (t) of y N (at + b) is estimated as By using Theorem (3.2) and by using Theorem (3.13) Suppose the nodes t k which are the roots of the SFGPs.By substituting these nodes in Eqs.(4.4)-(4.5);therefore the collocation scheme can be written as with the ICs This yields a set of (N+1) algebraic equations with the required SFG coefficients ỹS,j , j = 0, 1, … , N, which may be solved using an appropriate iterative procedure.As a result, the approximate solution y N (t) of the system (1.7-1.8)can be obtained.

Numerical applications
In this section, we'll introduce four applications to evaluate the applicability of our approach.By using a Dell laptop with an Intel(R) Core(TM) i3 CPU M 370@ 2.40 GHz and 3.00GB RAM configuration, the software was written in Mathematica version 12.The results obtained by our approach are compared to those obtained in literature using other numerical methods.The results constructed in this paper are measured by means of: • The absolute errors (AEs) given by (4.5) where y(t i ) and y N (t i ) are the exact solution (ES) and the approximate solution (AS), respectively.
• The maximum absolute errors (MAEs) • The order of convergence (OC) where error(N) denotes the error corresponding to polynomial degree N.

Application (1)
Consider the following VO-EF-DDEs [37], with the ICs where the ES at (t) = 2 is y(t) = e t .In Figure 1a, the AE function with N = 14, = 1, = 0.5 displays highly accuracy, and the numerical outcomes at different values of VOs are shown in Fig. 1b.
The numerical results of application ( 2) are given in Tables 1 and 2 and are plotted in Fig. 2a-b.Table 1 introduces the computational time at = 1, the MAEs at different values of N, , and the order of convergence at = 1.Addi- tionally, Table 2 presents the AEs at N = 3, = 1, and different values of .In Fig. 2a the AEs graph is displayed, while Fig. 2b depicts the agreement between the ES and the AS graphs of y(t) for N = 3, = 1, = 0.5 at (t) = 1 + t 2 .
1 2  The numerical results of application (3) are given in Tables 3, 4, and 5 and are plotted in Fig. 3a-b.Table 3 introduces the computational time at = 1, the MAEs at = 0.5, and different values of N, , and the order of convergence at = 1.Additionally, Table 4 presents the AEs at N = 4, = 1, and different values of .Table 5 displays the efficiency and effectiveness of our method by comparing it to the method in Ref. [38] at (t) = 2 .In Fig. 3a we describe the AEs graph, while Fig. 3b plots the ES and the AS graphs of y(t) for N = 4, = 1, = 0.5 at (t) = 1 + sin (t) 2 .

Application (4)
Consider the following system of EF-DDEs-VO [39], with the initial conditions, where The numerical results of application (4) are listed in Tables 6, 7, and 8 and graphically illustrated in Fig. 4a-d.We tabulate the computational time at = 1, and the MAEs in Table 6 using the SFGOMs method with different choices of at different values of N. In addition, Table 7 lists the AEs for various values of .The author of Ref. [39] esti- mated the AS of y 1 (t) and y 2 (t) , and at Table 8, we calculated the AEs from these results and compared them with the AEs of our technique to demonstrate the effective and precise of our method.To illustrate the convergence of the proposed technique, we plot the AEs between the ES and the AS of  y 1 (t) and y 2 (t) with N = 3, = 1 and = 0.5 in Fig. 4a and  c, respectively.The agreement between the ES and AS of y 1 (t) and y 2 (t) is depicted in Fig. 4b and d.

Conclusion
In this article, an effective and accurate numerical technique was used to find the approximate solutions of the VO-EF-DDEs and the system of VO-EF-DDEs.The numerical technique depended on the OMs of variable and integerorder derivatives of SFGPs in conjunction with the spectral collocation method.The used technique has the ability to transform the aforementioned problems into systems of algebraic equations which are easy to solve.Numerical studies for four applications are provided to demonstrate the accuracy, efficiency, and applicability of the used technique.The obtained numerical results indicated that satisfactory results are obtained by using a few terms of the SFGPs, also the efficiency of the estimated technique is increased by using more terms of SFGPs.The SFGPs are simple to compute, as well as they have a fast convergence rate.Another benefit of SFGPs is that for every choice of the SFG parameters a good numerical solution is obtained.Also the computational results of Sect. 5 showed that the value of the parameter affected on the numerical solutions as follows; when = 1 , good solutions are obtained by using few terms of SFGPs while for non-integer values of more terms of SFGPs are needed to obtain good results.So, we can conclude that the SGPs which are considered as a special case of SFGPs are more flexible and effective than SFGPs in the solved applications.Also, the obtained solutions do not require much CPU time.Finally, the material studied in this paper is novel, and the proposed technique can be extended to solve other fractional delay-differential equations.
Acknowledgements The authors are very grateful to the referees, for their careful reading of the manuscript and insightful comments, which help to improve the quality of the paper.

•
The explicit analytic form of SFGPs is given as • Also, the orthogonal relation of SFGPs where and ( ) S,j are defined by Eq. (2.7).• The SFGPs comprise unlimited number of orthogonal polynomials, among them the shifted fractional-order CPs of the first kind T ( ) S,j (t) ≡ C (0, ) S,j (t), the shifted frac- t i o n a l -o r d e r C P s o f t h e s e c o n d k i n d U ( ) 1, … , N, 0, o.w. .Proof From Eq. (2.15), we can write Since By using Newton's generalized binomial theorem, we can write this series converges for i ∈ Z + or | at b | < 1. ◻ By approximating the series (3.4) by using the first (N+1) terms, the relation (3.3) can be written as (3.2) d dt Φ (at + b) = d dt C ( , ) S,0 (at + b), C ( , ) S,1 (at + b), … , C ( , ) S,N (at + b) (at + b).

Fig. 1
Fig. 1 The numerical results of SFGOM technique on application (1) at N = 14, = 1 and = 0.5.a The AEs graph of y(t).b The graphs of the ASs at different values of (t).

Fig. 4
Fig. 4 The numerical behavior of the proposed method on application (4).a shows the AEs of y 1 (t) at N = 3 and = 0.5, (t) = 1 + t 2 .b shows the graph of the ES and AS y 1 (t) on [0,1].c shows the AEs of

Table 7
The AEs of application (4) at N = 3 and different values of

Table 8
Comparison of the AEs for N = 3,