Abstract
A new and efficient numerical approach is developed for solving nonlinear Lane-Emden type equations via Bernoulli operational matrix of differentiation. The fundamental structure of the presented method is based on the Tau method together with the Bernoulli polynomial approximations in which a new operational matrix is introduced. After implementation of our scheme, the main problem would be transformed into a system of algebraic equations such that its solutions are the unknown Bernoulli coefficients. Also, under several mild conditions the error analysis of the proposed method is provided. Several examples are included to illustrate the efficiency and accuracy of the proposed technique and also the results are compared with the different methods. All calculations are done in Maple 13.
1. Introduction
In recent years, the researches on the singular initial value problems (SIVPs) in several special second-order ordinary differential equations (ODEs) have received considerable attention among mathematicians and physicists. One of the most well-known classes of such equations are the Lane-Emden type equations which model many phenomena in mathematical physics and astrophysics. They are nonlinear ordinary differential equations which describe the equilibrium density distribution in self-gravitating sphere of polytrophic isothermal gas and have a singularity at the origin [1]. It must be noted that these equations have fundamental importance in the field of radiative cooling and modeling of clusters of galaxies. Moreover, it has been recently observed that the density profiles of dark matter halos are often modeled by the isothermal Lane-Emden equation with suitable boundary conditions at the origin [2]. Since getting the analytic solution of these equations is a difficult task in many cases, robust numerical schemes must be constructed for obtaining the approximated solutions. In this paper, we will present an efficient method for computing the numerical solution of the Lane-Emden type equations [3, 4] with the initial conditions where the prime denotes the differentiation with respect to , is a constant, and are nonlinear continuous functions. Selecting , , and yields [5, 6] which has another form subject to the initial conditions Similarly, isothermal gas spheres equation (in the case of ) is modeled by [5] with the zero values of and .
Recently, many analytical and numerical methods have been used to solve Lane-Emden type equations (1), (3), and (6). We note that the main difficulty arises in the singularity of the equations at . In [7, 8], Wazwaz employed the Adomian decomposition method with an alternate framework designed to overcome the difficulty of the singular point. The authors of [9] also applied pseudospectral method based on rational Legendre functions. Ramos [10–13] and Yousefi [14] solved Lane-Emden equation through several powerful numerical methods. Also, Yıldırım and Öziş [15] presented approximate exact solutions of a class of Lane-Emden type singular IVPs problems, by the VIM.
Since the beginning of 1994, the Bernoulli, Bernstein, Legendre, Taylor, Fourier, and Bessel matrix methods have been used in the works [2, 14, 16–26] to solve high-order linear and nonlinear differential (including hyperbolic partial differential equations) Fredholm Volterra integro-differential difference delay equations and their systems. The main characteristic of these approaches is based on the operational matrices of differentiation instead of integration. However, we can use high order Gauss quadrature rules such as [27, 28] for problems in different integration forms. To the best of our knowledge, this is the first work concerning the Bernoulli matrix method for solving nonlinear SIVPs such as Lane-Emden type equations (1), (3), and (6). This partially motivated our interest in such method. Therefore, in the light of the above-mentioned methods and by means of the matrix relations between the Bernoulli polynomials and their derivatives together with the idea of the Tau scheme, we develop a new method for solving nonlinear Lane-Emden type equations.
The Bernoulli polynomials play an important role in several branches of mathematical analysis, such as the theory of distributions in p-adic analysis [29], the theory of modular forms [30], the study of polynomial expansions of analytic functions [31], and so forth. Recently, new applications of the Bernoulli polynomials have also been found in mathematical physics, in connection with the theory of the Korteweg-de Vries equation [32] and Lamé equation [33], and in the study of vertex algebras [34].
In this paper, we are concerned with the direct solution technique for solving (1), (3), and (6) (with the associated initial conditions) using the Tau method based on Bernoulli operational matrix, such that it can be implemented efficiently and at the same time has a good convergence property. It should be noted that, the Tau method originate from vanishing the inner product of the residual and each of the test functions. By the aid of Tau scheme, the mentioned equations would be transformed into the systems of algebraic equations.
The remainder of the paper is organized as follows. In Section 2 we introduce some mathematical preliminaries of Bernoulli polynomials and also their operational matrix of differentiation. Section 3 is devoted to applying the Tau method for solving Lane-Emden type equations using the Bernoulli operational matrix. In Section 4 we provide the error analysis of the proposed method. In Section 5 the proposed method is applied to several examples together with a full comparison to other recent methods. Also conclusions are given in Section 6.
2. Bernoulli Polynomials and Their Operational Matrix of Differentiation
Bernoulli polynomials and their applications can be found in number theory and classical analysis initially. They also appear in the integral representation of the differentiable periodic functions, since they are employed for approximating such functions in terms of polynomials. They are also used for representing the remainder term of the composite Euler-Maclaurin quadrature rule [16].
In this section, we recall some properties of the Bernoulli polynomials which will be of fundamental importance in the sequel.
Property 1 (differentiation: see [16]). , .
Property 2 (integral means conditions: see [16]). , .
Property 3 (differences: see [16]). , .
Property 4 (monomials representation: see [16]). , .
If we introduce the Bernoulli vector in the form , then the derivative of the , with the aid of the first property, can be expressed in the matrix form by where is the operational matrix of differentiation.
Accordingly, the th derivative of can be given by where is defined in (7).
3. Implementation of the Tau Scheme Using Bernoulli Operational Matrix
In this section, we apply the Tau method for solving the Lane-Emden type equations numerically. Therefore, we first assume that the solution of the Lane-Emden type equations can be written in terms of linear combination of Bernoulli polynomials (with unknown coefficients). Then, inner product of the residual with each of the test functions should be assumed to zero. Thus, the considered problem would be transformed into a system of nonlinear algebraic equations in which its solutions are the associated Bernoulli coefficients. Again consider the Lane-Emden equation with the initial conditions Now we approximate , and by Bernoulli polynomials in the following forms where the unknowns are meanwhile are known and are the Bernoulli polynomials.
Using operational matrix of differentiation of Bernoulli polynomials, (1) can be written approximately as follows: The residual for (12) can be written as Applying the typical Tau method, which is used in the sense of a particular form of the Petrov-Galerkin method [4], (12) can be converted in nonlinear equations by applying The initial conditions are given by Equations (14) and (15) generate sets of nonlinear equations in terms of unknown Bernoulli coefficients , . These nonlinear equations can be solved using Newton’s iterative method that was implemented in command of MAPLE software for obtaining the unknown coefficients of the row vector and hence the solution can be approximated easily by .
4. Error Analysis and Accuracy of the Solution
In this section, we will illustrate convergence of the proposed method assuming the known functions and also the unknown solution are in the space (where is the degree of smoothness of the problem) with bounded derivatives. But some lemmas should be recalled from the literature and then the main theorem of this section would be provided.
Lemma 1. Assume that is an enough smooth function and also is approximated by the Bernoulli series , then the coefficients for all can be calculated from the following relation:
Proof. See [20].
Lemma 2. Assume that one approximates the function on the interval by Bernoulli polynomials as argued in Lemma 1. Then the coefficients decays as follows: where .
Proof. Since it is trivial we omit the proof.
The previous lemma implies that Bernoulli coefficients are decayed rapidly as increasing of under the condition of boundedness of all derivatives of with respect to . We will illustrate this fact experimentally in our numerical examples. It should be noted that, there are some examples such as which has bounded derivatives, but does not go to zero rapidly.
Theorem 3 (see [37]). Suppose that is in the space and is approximated by Bernoulli polynomials as done in Lemma 1. With more details assume that is the approximate polynomial of in terms of Bernoulli polynomials and is the remainder term. Then, the associated formulas are stated as follows: where and denotes the largest integer not greater than .
Proof. See [37].
Lemma 4. Suppose (with bounded derivatives) and is the approximated polynomial using Bernoulli polynomials. Then the error bound would be obtained as follows: where denotes a bound for all the derivatives of function (i.e., , for ) and is a positive constant.
Proof. See [25].
Next, we will provide the main theorem of this section. We show that, if both of the and are approximated by the Bernoulli polynomials in the Lane-Emden equation (1), then the error of the approximation of depends directly on the approximation of . Therefore, using enough values of may gives us high order approximation of the desired solutions.
Theorem 5. Assume that and , where is a linear integral operator. If we approximate and by and , respectively, using Bernoulli polynomials, then where is the Lipschitz constant of the function with respect to its second variable and also .
Proof. Consider the integral operator . Applying to both sides of (1) yields Considering the following assumptions: transforms (21) into the following Volterra integral equations: Now approximate both functions and by and using Bernoulli polynomials. Therefore Thus In other words and this completes the proof.
5. Numerical Illustrations
In this section, several numerical examples are given to illustrate the accuracy and effectiveness of the proposed method. All calculations are done on a 64 bits PC laptop. As we claimed in Section 2, our method which is based on the Bernoulli polynomials has more efficiency with respect to the Legendre method [4] in isothermal gas spheres equation (i.e., Example 3) and also has the same accuracy with respect to methods which use high-order Gauss quadrature rules [36]. We must recall that, the results of our method and the method of [36] are the same in Examples 5 and 6, but our results were obtained in less CPU time; meanwhile the results of [36] were obtained in more CPU time by using the same PC and equipment. Moreover in Table 5, we show the vanishing of the Bernoulli coefficients (as shown theoretically in Lemma 2) for all the examples of this section experimentally.
Example 1 (see [4]). At first we consider the equation with the initial conditions , . Trivially the exact solution of this equation is .
Note that , and . We apply the method that was explained in Section 3 for . Thus assume that
Our aim is to determine the unknown Bernoulli coefficients , , and by using Tau method. Also the operational matrix of Bernoulli polynomials and its square are as follows:
By using (14) we have
where
Therefore
By applying the initial conditions we have
Solving (31)-(32) yields , , . Thus
which is the exact solution.
Example 2 (see [4]). We now consider the following Lane-Emden equation:
with initial conditions , which has the exact solution .
Again, we apply our proposed method for solving the above equation for different values of such as and . In other words we assume that and . The numerical results of our method are depicted in Figure 1. In the Figure 1(a), the exact solution together with the numerical solutions and is illustrated in a larger computational interval. Moreover the error history of these approximated solutions is shown in the Figure 1(b). From this figure one can conclude that our method obtained highly accurate solutions even in large computational intervals.
(a)
(b)
Example 3 (see [1, 4, 7]). As the third example we consider the isothermal gas spheres equation as follows:
with the initial conditions and .
In this case we have and . We approximate by using the five terms of its Maclaurin expansion (i.e., ). Again we use instead of in the procedure of approximation. In other words
We apply our method for solving this problem by using several values of . We provide the numerical solutions at the points , and in the case of in Table 1 together with the exact solution that was reported in [7]. Moreover we make an interesting comparison with the methods that are based on Legendre operational matrix of differentiation [4] and the Hermite collocation method [1] in the associated values of the errors by the assumption of . From these comparisons we see that the presented method (PM) has more efficiency with respect to the above-mentioned methods.
Example 4 (see [1, 4, 7]). As the fourth example we consider the following Lane-Emden equation:
with the initial conditions and .
Similar to the third example we have and . Approximating by using the three terms of its Maclaurin expansion yields and also by applying the Bernoulli polynomials in the procedure of approximations we have
Again the method that is proposed in Section 3 will be used. The numerical solutions at the points , and in the case of are written in Table 2 together with the exact solution that was reported in [7]. Moreover the error of the presented method (PM) together with the errors of [1, 4] by the assumption of is provided. Evidently the results of our methods are very close to [4] and are superior with respect to the results that were obtained in [1].
Example 5 (see [1, 35, 36]). We now consider the Lane-Emden equation with the initial conditions , .
The exact solution of this equation was reported in [35]. The numerical results of our scheme together with two other methods [1, 36] are provided in Table 3. Not only does our method needs to lower values of Bernoulli polynomials (as the test functions) with respect to the methods that are based on Gauss quadrature rules [36] for obtaining high accurate solutions, but also the CPU time of our method for solving the associated system of nonlinear algebraic is mush less than others such as [36]. The basic reason for this claim is based upon using the operational matrices of differentiation and this fact leads to sparse equations.
Example 6 (see [1, 35, 36]). As the final example we consider the Lane-Emden equation with initial conditions , .
Again we recall that the exact solution of this equation was reported in [35]. The numerical results of this Example are given in Table 4. Similar to the previous example we can see that the results of our approach are the same [36] by using lower values of Bernoulli polynomials with respect to the shifted Legendre polynomials and also are superior with respect to the results of [1].
Hint
As we claimed in Lemma 2 the Bernoulli coefficients in the procedure of approximating any arbitrary continuous function must tend to zero as the number of the test functions (i.e., the Bernoulli polynomials) tend to infinity. This fact is shown experimentally in Table 5.
6. Conclusions
The Lane-Emden type equations describe a variety of phenomena in theoretical physics and astrophysics, including the aspects of stellar structure, the thermal history of a spherical cloud of gas, isothermal gas spheres, and thermionic currents. Lane-Emden type equations have been considered by many mathematicians as mentioned before [4]. The fundamental goal of this paper has been to construct an approximation to the solution of nonlinear Lane-Emden type equations in the computational interval . A set of Bernoulli polynomials is proposed to provide an effective but simple way to improve the convergence of the solution by the Tau method. The validity of the method is based on the assumption that it converges by increasing the number of Bernoulli polynomials. A comparison is made among the numerical (and exact) solutions of [35] and the series solutions of [7], Legendre operational matrix [4], Hermite collocation method [1] and high-order Gauss quadrature rule [36], and the current work. It has been shown that the present work provides acceptable approach for Lane-Emden type equations.
Conflict of Interests
The authors declare that they do not have any conflict of interests in their submitted paper.
Acknowledgments
The authors thank the editor and referees for their valuable comments and helpful suggestions which led to the improved version of the paper.