Polytropic Dynamical Systems with Time Singularity

In this paper we consider a class of second order singular homogeneous differential equations called the Lane-Emden-type with time singularity in the drift coefficient. Lane-Emden equations are singular initial value problems that model phenomena in astrophysics such as stellar structure and are governed by polytropics with applications in isothermal gas spheres. A hybrid method is proposed to approximate the solution of this type of dynamic equations.


Introduction
Laplace's equation and Poisson's equation are important examples of elliptic partial differential equations which used broadly in applied mathematics and theoretical physics, see, e.g., [22].For instance, Poisson's equation used to calculate gravitational field in potential theory and can be seen as generalization of Laplace's equation.By removing or reducing dimensions from Poisson's equation, we obtain a second-order nonlinear differential equation called Lane-Emden-type equation (LE, for short).The Lane-Emden equation (a.k.a.polytropic dynamic equation) is one of the well studied classical dynamical systems that has many applications in nonlinear mathematical physics and non-Newtonian fluid mechanics (see, for instance, [2,3,6,10,11,21,26]).A preliminary study on the LE equations (polytropic and isothermal) was undertaken by astrophysicists Lane (1870) and Emden (1907), such that the interest of the LE derived from its nonlinearity and singular behavior at the origin.The point x 0 is called ordinary point (or regular point) of the dynamic equation (2) if the coefficients of x, x ′ are analytic in an interval about x 0 .Otherwise, it is called singular point.In solving singular boundary value problems (BVPs) some numerical techniques are based on the idea of replacing a two-point BVP by two suitable initial value problem [14,21,25].In this paper we adopt such idea (called the shooting method) to study dynamical models that play an essential role in the theory of star structure and evolutions, thermodynamics, and astrophysics (see, e.g., [9]).Equation (1) describes and models the mechanical structure of a spherical body of gas such as a self-gravitating star and also appeared in the study of stellar dynamics (see, for instance [8,11] and the references therein).The solutions to the LE, which are known as polytropes, are functions of density versus the radius expressed by x(t) in (2).The index n determines the order of that solution.Nonlinear singular LE equations can be formulated as or, subject to x(0) = 1 , x ′ (0) = 0.
The dynamical system model ( 2) along with initial conditions form a special type of initial value problems (IVP) for which it has several applications in the fields of celestial mechanics, quantum physics and astrophysics [6,10,14,26].The following figure is a motivation example shows finite solutions of Lane-Emden equation for the value of n in equation ( 1) or (2) given by n = 0, 1, 2, 3, 4, 5, 6.
For some special cases when n = 0, 1, 5 exact analytical solutions were obtained by Chandrasekhar [8], while for all other values of n approximate analytical methods were obtained such as: the Adomian decomposition method [20,27], homotopy analysis method [5], power series expansions [16], variational method [13], and linearization techniques [23] (provide accurate closed-form solutions around the singularity.).Numerical discretization for equation (1) has been the object of several studies in the last decades (see, e.g., [?, [1][2][3]6,10,21,25,26] and the references therein).In [16], the authors presented numerical method for solving singular IVPs by converting Lane-Emden-type equation ( 1) to an integral operator form then rewriting the acquired Voltera integral equation in terms of a power series.Ramos [23] applied linearization method for the numerical solution of singular initial value problems of linear and nonlinear, homogeneous and nonhomogeneous second-order dynamic equations.Russell and Shampine in [25] discussed the solution of the singular nonlinear BVP for certain dynamical systems in the context of analytical geometry and symmetry as follows and with boundary conditions x ′ (0) = 0 (or equivalently x(0) is finite), x(b) = λ, for some scalar λ, and the convergence is uniform over the interval [0, 1].Biles et al. in [6], have considered an initial value problem for Lane-Emden type of the form where a ∈ R and p(t) may be singular at t = 0.They introduced the following definition and theorem, respectively; where the theorem gives the conditions of existence and uniqueness of solution of second-order linear BVPs.
x is a solution of the above equation (4) if and only if there exist some T > 0, such that x, x ′ are absolutely continuous on [0,T].
Our paper is organized in the following fashion.In section 2, we provide some necessary notations and essential background.In section 3 we present the second-order dynamical system of Lane-Emden type, and the BVP is transformed to IVP by shooting method.Then applying Euler's method on the resulted initial value problem to get approximations for the solution of the LE.The convergence results and error estimation are analyzed in section 4. Finally, numerical examples are provided to demonstrate the validity and efficiency of the proposed technique.

Preliminaries
In this section we introduce some basic definitions and conventional notations.Let C 1 (I) be the space of all continuously differentiable functions defined on an interval I.A set D in the Euclidean space R n is compact set if and only if it is closed and bounded set.The basic space used throughout this paper is the space of continuous functions C[0, 1] on the compact set [0, 1] with the associated norm (distance) function defined by, In particular, all C 1 functions are locally Lipschitz.The following two theorems address existence and uniquness of solutions to any IVP.
has a unique solution x = x(t).
To better understand the theorem we illustrate it by giving an example on the interval [1,2] instead of [0, 1]: Consider the BVP, x ′′ (t) + sin x ′ + e −tx = 0 with x(1) = x(2) = 0 and t ∈ [1,2].Now apply the theorem to then the condition is satisfied and the BVP has a unique solution.Now reader might ask how can we apply this Theorem to Lane-Emden equation.Theorem 2.2 can be simplified by taking into account that the functions sin x ′ x ′ and e −tx are continuous on the interval (0, ∞) to assure the differential equation is linear.

Computational Mehtods For Dynamical systems
In this section, we start by presenting the methods (shooting to transform from BVP to IVP, and Euler's for regular singularity in the drift term) and apply them on the second order singular dynamical system.

Shooting method
The shooting method treats the two-point BVP as an IVP.The idea basically, is to write the BVP in a vector form and begin the solution at one end of the BVP, and then "shooting" to the other end with any IVP solver, such as; Runge-Kutta method or multistep method for linear case and Secant method or Newton's method for nonlinear case, until the boundary condition at the other end converges to its correct value.To be precise, the ordinary differential equation of second order, associated with its initial conditions must normally be written as a system of first order equation before it can be solved by standard numerical methods.Next figure shows graphically the mechanism of the shooting.
Roughly speaking, we 'shoot' out trajectories in different directions until we find a trajectory that has the desired boundary value.The drawback of the method is that it is not as robust as those used to solve BVPs such as finite difference or collocation methods presented in [1,25], and there is no guarantee of convergence.Shooting method can be used widely for solving a BVP by reducing it to an associated IVP, and is valid for both linear (also called chasing method) and non linear BVPs, by [18], Next theorem provides existence and uniqueness to the BVP's solution.
and assume f is continuous function on D such that it satisfies the BVP: Suppose that f x and f x ′ are continuous on the same set D.

and
(ii) There exists M > 0 such that then the BVP (8) has a unique solution.
A special case of this theorem is the following corollary, i.e., when the right hand side of (8) is linear.For linear Lane-Emden equations, one can use Frobenius method to determine the analytical solutions of (1) near the singularity, see, for instance, [23].
Corollary 3.2.Consider (8) given by and the time-dependent coefficients p(t) , q(t) , r(t) are continuous functions on the domain [a, b] and further q(t) > 0, then the BVP (8) has a unique solution.
Proof.We need to consider two cases: (i) When equation ( 9) given with boundary conditions x(a) = α , x ′ (a) = 0, has a unique solution x 1 (t).(ii) When equation ( 9) with r(t) = 0, x(a) = 0 , x ′ (a) = 1, has a unique solution x 2 (t).Therefore, one can easily check that the linear combination x 2 (t) is the unique solution to (9), and hence to (8) due to the existence and uniqueness guaranteed by Picard-Lindelof theorem (2.1).

Euler's Method
Euler's method is a numerical approach for solving (iteratively) initial value problems, as follows: We divide the time interval [t 0 , T ] into N equal subintervals, each of length h = ∆t = t n+1 −t n , for n ≥ 0, and start by initial value x(0) then move forward using the step size towards x(T ), that is, given the second-order ordinary differential equation (7), converting it into two first-order dynamic equations (i.e., dynamical system).Discretize the interval [t 0 , T ] into subintervals, and by assuming y n the approximation to x(t n ) and v n the approximation to u(t n ).Euler's method is then can be expanded, as a twoterms truncated Taylor series, by the following Euler's method for solving a second-order differential equation is given by: Forward Euler's Algorithm.
Step 1. (Forward step): Given t n , y n , v n define The local error at every step is proportional to the square of the step size h and the global error at a given time is proportional to h.Moreover, the order of the global error can be calculated from the order of the local error ( i.e. by summing up the local error).We can understand Euler's method by appealing the idea that some differential equations provide us with the slope at all points of the function , while an initial value provides a point on the function.Using this information we can approximate the function with a tangent line at the initial point.It is known that the tangent line is only a good approximation over a small interval.When moving to a new point, we can construct an approximate tangent line, using the actual slope of the function, and an approximation to the value of the function at the tangency point.Repeating this manner, we eventually construct a piecewise-linear approximation to the solution of the differential equation.Moreover, this approximation can be seen as a discrete function and to make it a continuous function, we interpolate (linearly) between each pair of these points.
In the following, we study and analyse the Lane-Emden-type equation with an endpoint singularity in terms of the independent variable which has the form where â(t, x(t), x ′ (t)) : [0, 1) × R × R → R, and the Lipschitz functions a(t, x), g(t, x) ∈ C 1 ([0, 1) × R), for all 0 ≤ t < 1.At t = 1, the −a(t, x) 1 − t term is singular, but symmetry implies the boundary condition x ′ (0) = 0.With this boundary condition, the term −a(t, x) 1 − t dx dt is well defined as t → 1.The solution of ( 10) can be given by the system: Define x t := x(t), x ′ t := x ′ (t).By the fundamental theorem of calculus and provided that all integrals are exist (finite), we notice that equation ( 11) is equivalent to the nonlinear system of integral equations: Where 0 = t 0 < t 1 < t 2 < ... < 1.
Expanding the integrands in ( 12) so we have: Or in the equivalent form, For simplicity we assume t tn s tn â(u, x u , x ′ u ) du ds.
Thus the system becomes, where h n+1 = t n+1 − t n .
In order to estimate the error, we need to find a bound for the integrands in L n and L (2) n .The double integrals in both L (1) , L (2) yield the local truncation error, if we define the numerical value by: where h n+1 = t n+1 − t n .

Discretization and Convergence Analysis
Consider a sequences of times 0 = t 0 < t 1 < t 2 < ... < 1, and the corresponding step sizes h n = t n − t n−1 .Define x n = x(t n ) and x ′ n = x ′ (t n ) where (x(t), x ′ (t)) is a solution of (5).Writing (8) in the form: Use y n as defined in ( 9) and let n By using the inequality (x + y) 2 ≤ 2x 2 + 2y 2 , the error can be estimated as, Next, we introduce some assumptions on the functions a(t, x(t)), g(t, x(t)) and their partial derivatives for t ∈ [0, 1), x ∈ R .But before that we remind ourselves of the value of â from section 3, â(t, x(t), x ′ (t)) = −a(t, x(t)) Also, for any T 1 , T 2 ∈ [0, 1) the Lipschitz conditions are: Our required bounds explicitly are: The partial derivatives bounds are: This final bound applies along the path Taking the difference between the computed and the exact values of â, By adding and subtracting the required terms, we have Thus, the difference 18 becomes, Note that, We now apply a very well known result from functional analysis, Cauchy-Schwarz inequality twice on L (1) andL (2) : for some Constant D 1 , which does not depend on h n+1 and n.
where D 2 is independent of n and h n+1 .
To avoid the singularity and produce a better estimation to test the efficiency of the algorithm, we introduce a variable step size by fixing ĥ > 0 and then defining step size h n and node points t n using ĥ: ĥ In the process of estimating the global error, we need to use the following two fundamental lemmas: Lemma 4.1.For all x ≥ −1, and any m > 0, we have 0 The proof of this result follows by applying Taylor's theorem with f (x) = e x , x 0 = 0, and n = 1.Lemma 4.2.if M 1 ≥ −1 and M 2 ≥ 0 are real numbers and {a n } N n=0 is a sequence with a 0 ≥ 0 such that then, Proof.Fix a positive integer n, then (20) can be written as Now if we add the two inequalities in (11) together, we will have Using the definition of the norm ∥ϵ n ∥ = (ϵ ′ n ) 2 + (ϵ n ) 2 , then system (13) can be simplified as ĥ) 3 where m 1 and m 2 are independent constants of h n+1 and t n+1 .Now we apply Lemma 4.2 for a n = ∥ϵ n ∥ 2 , followed by a foundation for the step size order, with M 1 = 1 + m 1 ( ĥ) and M 2 = m 2 ( ĥ) 3 such that if The following theorem can assur the variable step size and the uniform convergence for solutions of the method.

Simulation and Numerical Experiments
In this section we run the algorithm over some examples to show the validity of the method.We used MATLAB with bulit-in functions such as; ode45 and EulerSolver Example 5.1.Consider the second order differential equation (10) with a(t, x) = sin x, and g(t, x) = x 5 , where the step size is 0.05 and time interval [0, 1] along with initial conditions x(0) = 0, x ′ (0) = 2; i.e., Table 1 compares the two dependent solutions x(t) and x ′ (t) for equation (10) given the above numerical values, and figures below draw the relationships between trajectories of the differential equation and the time.The analytical solution to this problem is somewhat lower than our approximation.By shrinking the size of the interval ∆t, we could calculate a more accurate estimate.Example 5.5.In this example we consider the non-autonomous inhomogeneous second order system with the right-hand side being t 3 e 2t , a(t, x) = 4, and g(t, x) = 4x, where the step size is 0.01 and along with initial conditions x(0) = 0, x ′ (0) = 0; with the absence of singularity.The graphs shown below and the tables as well.

Conclusion and Extensions
In this paper our primary goal was to investigate the second-order singular Lane Emden type equations and we have successfully arrived at the solutions by the forward Euler's algorithm combined with the shooting method, which in turn, reduces the boundary value problem into initial value problem, so the method showed that it is a precise and timesaving method.The Lane Emden equations are solved for the values of the polytropic indices varies from 1, 2, 3 and 5 with having constants, linear functions and periodic functions in the drift term.The numerical solution of the problem for these values of indices replaces the unsolvable version of equation and any closed form solution that we wish to find.For the case of n = 2 the solution is obtained as an infinite power series.Graphical representations of these results give us information about polytropes for different values of polytropic indices which may be helful in the study of the behavior of stellar structures in astrophysics.One good extension for this work is through implementing backward Euler formula for a second-order differential equations where the recursion formula is the same, except that the dependent variable is a vector.Another possible modification for the work is by using the reliable Runge-Kutta method which promises accurate results in deriving the solutions of the Lane Emden equations.It is also significant in handling highly nonlinear differential equations with less computations and a larger interval of convergence.For thinking globally, finite difference methods may be used to replace the shooting method to treat the boundary value problem.Finally, we may think of adding the additive noice to the second order differential equation (it will be called stochastic differential equation) and in this case, Euler's method will be replaced by Euler-Maruyama Algorithm, see, for instance, [12,15].

Figure 1 :
Figure 1: Comparison between approximated solution by Euler's method and the actual solution for the equation x ′′ + 4x ′ + 4x = t 3 e 2t .
Define a continuous function f : D → R n where D is an open subset of R n+1 , and consider Let D be a nonempty set.Suppose there is a function f from D to itself, and 0 ≤ L < 1, where L is free of x and y.If for any two points x, y ∈ D we have|f (x) − f (y)| ≤ L|x − y| , ∀ x, y ∈ D,then f is called a contraction.The smallest such value of L is called the Lipschitz constant of f , and f is then called a Lipschitz function.
Definition 2.2.A function f : D ⊂ R n+1 → R n is said to be locally Lipschitz in x if for each compact set contained in D, and each x, y ∈ D, there exists L > 0 such that