Journal of Applied Mathematics

Volume 2013 (2013), Article ID 850170, 9 pages

http://dx.doi.org/10.1155/2013/850170

Research Article

A New Tau Method for Solving Nonlinear Lane-Emden Type Equations via Bernoulli Operational Matrix of Differentiation

E. Tohidi,¹ Kh. Erfani,¹ M. Gachpazan,² and S. Shateyi³

¹Department of Mathematics, Islamic Azad University, Zahedan Branch, Zahedan, Iran

²Department of Applied Mathematics, School of Mathematical Sciences, Ferdowsi University of Mashhad, Mashhad, Iran

³Department of Mathematics, University of Venda, Private Bag X5050, Thohoyandou 0950, South Africa

Received 18 February 2013; Accepted 24 April 2013

Academic Editor: Mehmet Sezer

Copyright © 2013 E. Tohidi et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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)

Figure 1: (a) Comparison of the exact and approximated solutions and (b) history of the errors.

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.

Table 1: The exact and approximated solution of Example 3 together with the error comparisons for .

Solution of PM
Exact Solution [7]
Error of PM
Error of [4]
Error of [1]

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].

Table 2: The exact and approximated solution of Example 4 together with the error comparisons for .

Solution of PM
Exact Solution [7]
Error of PM
Error of [4]
Error of [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.

Table 3: The exact and approximated solution of Example 5 together with the error comparisons.

Solution of PM
Exact Solution [35]
Error of PM
Error of [1]
Error of [36]

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].

Table 4: The exact and approximated solution of Example 6 together with the error comparisons.

Solution of PM
Exact Solution [35]
Error at
Error of [1]
Error [36]

Table 5: Vanishing of the Bernoulli coefficients as the number of index is increased.

Example 2
Example 3
Example 4
Example 5
Example 6

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.

References

K. Parand, M. Dehghan, A. R. Rezaei, and S. M. Ghaderi, “An approximation algorithm for the solution of the nonlinear Lane-Emden type equations arising in astrophysics using Hermite functions collocation method,” Computer Physics Communications, vol. 181, no. 6, pp. 1096–1108, 2010.

N. Kumar, R. K. Pandey, and C. Cattani, “Solution of the Lane-Emden equation using the bernstein operational matrix of integration,” ISRN Astronomy and Astrophysics, vol. 2011, Article ID 351747, 7 pages, 2011.

R. K. Pandey and N. Kumar, “Solution of Lane-Emden type equations using Bernstein operational matrix of differentiation,” New Astronomy, vol. 17, pp. 303–308, 2012.

R. K. Pandey, N. Kumar, A. Bhardwaj, and G. Dutta, “Solution of Lane-Emden type equations using Legendre operational matrix of differentiation,” Applied Mathematics and Computation, vol. 218, no. 14, pp. 7629–7637, 2012.

H. T. Davis, Introduction to Nonlinear Differential and Integral Equations, Dover Publications Inc., New York, NY, USA, 1962.

J. H. Lane, “On theoretical temperature of the sun under the hypothesis of a gaseous mass maintaining its internal heat and depending on the laws of gases known to terrestrial experiment,” American Journal of Science and Arts, 2nd Series, vol. 50, pp. 57–74, 1870.

A.-M. Wazwaz, “A new algorithm for solving differential equations of Lane-Emden type,” Applied Mathematics and Computation, vol. 118, no. 2-3, pp. 287–310, 2001.

A.-M. Wazwaz, “The modified decomposition method for analytic treatment of differential equations,” Applied Mathematics and Computation, vol. 173, no. 1, pp. 165–176, 2006.

K. Parand, M. Shahini, and M. Dehghan, “Rational Legendre pseudospectral approach for solving nonlinear differential equations of Lane-Emden type,” Journal of Computational Physics, vol. 228, no. 23, pp. 8830–8840, 2009.

J. I. Ramos, “Linearization techniques for singular initial-value problems of ordinary differential equations,” Applied Mathematics and Computation, vol. 161, no. 2, pp. 525–542, 2005.

J. I. Ramos, “Piecewise-adaptive decomposition methods,” Chaos, Solitons & Fractals, vol. 40, no. 4, pp. 1623–1636, 2009.

J. I. Ramos, “Series approach to the Lane-Emden equation and comparison with the homotopy perturbation method,” Chaos, Solitons & Fractals, vol. 38, no. 2, pp. 400–408, 2008.

J. I. Ramos, “Linearization methods in classical and quantum mechanics,” Computer Physics Communications, vol. 153, no. 2, pp. 199–208, 2003.

S. A. Yousefi, “Legendre wavelets method for solving differential equations of Lane-Emden type,” Applied Mathematics and Computation, vol. 181, no. 2, pp. 1417–1422, 2006.

A. Yıldırım and T. Öziş, “Solutions of singular IVPs of Lane-Emden type by the variational iteration method,” Nonlinear Analysis: Theory, Methods & Applications, vol. 70, no. 6, pp. 2480–2484, 2009.

A. H. Bhrawy, E. Tohidi, and F. Soleymani, “A new Bernoulli matrix method for solving high-order linear and nonlinear Fredholm integro-differential equations with piecewise intervals,” Applied Mathematics and Computation, vol. 219, no. 2, pp. 482–497, 2012.

K. Erdem and S. Yalçinbaş, “Bernoulli polynomial approach to high-order linear differential-difference equations,” AIP Conference Proceedings, vol. 1479, pp. 360–364, 2012.

K. Erdem and S. Yalçinbaş, “A Bernoulli polynomial approach with residual correction for solving mixed linear Fredholm integro-differential-difference equations,” Journal of Difference Equations and Applications, 2013.

M. Sezer, S. Yalçinbaş, and N. Şahin, “Approximate solution of multi-pantograph equation with variable coefficients,” Journal of Computational and Applied Mathematics, vol. 214, no. 2, pp. 406–416, 2008.

E. Tohidi, A. H. Bhrawy, and K. Erfani, “A collocation method based on Bernoulli operational matrix for numerical solution of generalized pantograph equation,” Applied Mathematical Modelling, vol. 37, no. 6, pp. 4283–4294, 2013.

E. Tohidi, “Legendre approximation for solving linear HPDEs and comparison with Taylor and Bernoulli matrix methods,” Applied Mathematics, vol. 3, no. 5, pp. 410–416, 2012.

E. Tohidi, F. Soleymani, and A. Kilicman, “Robustness of operational matrices of differentiation for solving state-space analysis and optimal control problems,” Abstract and Applied Analysis, vol. 2013, Article ID 535979, 9 pages, 2013.

E. Tohidi and A. Kilicman, “A collocation method based on the Bernoulli operational matrix for solving nonlinear bvps which arise from the problems in calculus of variation,” Mathematical Problems in Engineering, vol. 2013, Article ID 757206, 9 pages, 2013.

F. Toutounian, E. Tohidi, and A. Kilicman, “Fourier operational matrices of dierentiation and transmission: introduction and applications,” Abstract and Applied Analysis, vol. 2013, Article ID 198926, 11 pages, 2013.

F. Toutounian, E. Tohidi, and S. Shateyi, “A collocation method based on Bernoulli operational matrix for solving high order linear complex differential equations in a rectangular domain,” Abstract and Applied Analysis, vol. 2013, Article ID 823098, 12 pages, 2013.

S. Yüzbaşı, “A numerical approach for solving the high-order linear singular differential-difference equations,” Computers & Mathematics with Applications, vol. 62, no. 5, pp. 2289–2303, 2011.

O. R. N. Samadi and E. Tohidi, “The spectral method for solving systems of Volterra integral equations,” Journal of Applied Mathematics and Computing, vol. 40, no. 1-2, pp. 477–497, 2012.

E. Tohidi and O. R. N. Samadi, “Optimal control of nonlinear Volterra integral equations via Legendre polynomials,” IMA Journal of Mathematical Control and Information, vol. 30, no. 1, pp. 67–83, 2013.

N. Koblitz, p-Adic Numbers, p-Adic Analysis and Zeta-Functions, vol. 58 of Graduate Texts in Mathematics, Springer, New York, NY, USA, 2nd edition, 1984.

S. Lang, Introduction to Modular Forms, Springer, Berlin, Germany, 1976.

R. P. Boas and R. C. Buck, Polynomial Expansions of Analytic Functions, Springer, New York, NY, USA, 1964.

D. B. Fairlie and A. P. Veselov, “Faulhaber and Bernoulli polynomials and solitons,” Physica D, vol. 152-153, pp. 47–50, 2001.

M.-P. Grosset and A. P. Veselov, “Elliptic Faulhaber polynomials and Lamé densities of states,” International Mathematics Research Notices, vol. 2006, Article ID 62120, 31 pages, 2006.

B. Doyon, J. Lepowsky, and A. Milas, “Twisted vertex operators and Bernoulli polynomials,” Communications in Contemporary Mathematics, vol. 8, no. 2, pp. 247–307, 2006.

G. P. Horedt, Polytropes: Applications in Astrophysics and Related Fields, Kluwer Academic Publishers, Dordrecht, The Netherlands, 2004.

A. M. Rismani and H. Monfared, “Numerical solution of singular IVPs of Lane-Emden type using a modified Legendre-spectral method,” Applied Mathematical Modelling, vol. 36, no. 10, pp. 4830–4836, 2012.

M. Abramowitz and I. A. Stegun, Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables, National Bureau of Standards, Wiley, New York, NY, USA, 1972.