- Email: [email protected]
A numerical approach for solving singular nonlinear Lane–Emden type equations arising in astrophysics
New Astronomy 34 (2015) 178–186
Contents lists available at ScienceDirect
New Astronomy journal homepage: www.elsevier.com/locate/newast
A numerical approach for solving singular nonlinear Lane–Emden type equations arising in astrophysics A. Kazemi Nasab, A. Kılıçman ⇑, Z. Pashazadeh Atabakan, W.J. Leong Department of Mathematics, University Putra Malaysia (UPM), 43400 Serdang, Selangor, Malaysia
h i g h l i g h t s To present an efficient method to solve Lane–Emden type equations. To represent smooth and especially piecewise smooth functions properly. To improve the accuracy by properly increasing either the number of subintervals or the number of collocation points. To compare present solutions for more accuracy and efficiency.
a r t i c l e
i n f o
a b s t r a c t
Article history: Received 16 May 2013 Received in revised form 28 March 2014 Accepted 17 June 2014 Available online 11 July 2014 Communicated by W. Soon Keywords: Wavelet analysis method Chebyshev finite difference method Lane–Emden type equations Singular IVPs Isothermal gas spheres
In this paper, we suggest a numerical method based upon hybrid of Chebyshev wavelets and finite difference methods for solving well-known nonlinear initial-value problems of Lane–Emden type. The useful properties of the Chebyshev wavelets and finite difference method are utilized to reduce the computation of the problem to a set of nonlinear algebraic equations. Making a comparison among the obtained results using the present method with those ones reported in literature by some other well-known methods confirms the accuracy and computational efficiency of the present technique. Ó 2014 Elsevier B.V. All rights reserved.
1. Introduction In recent years, solving singular initial value problems formulated by second order nonlinear differential equations has gained a great deal of attention from many scientists in mathematics and physics. Among equations of this type, the Lane–Emden type equations are well-known and are used as the benchmark to verify new methods. The Lane–Emden differential equation of index m was first introduced by Lane in 1870 (Lane, 1870) and then studied more in detail by Emden (1907), describing the equilibria of nonrotating polytropic fluids in a self-gravitating star and formulated as following:
(
y00 ðtÞ þ 2t y0 ðtÞ þ yðtÞm ¼ 0; yð0Þ ¼ 1;
tP0
y0 ð0Þ ¼ 0:
⇑ Corresponding author. E-mail address: [email protected] (A. Kılıçman). http://dx.doi.org/10.1016/j.newast.2014.06.008 1384-1076/Ó 2014 Elsevier B.V. All rights reserved.
ð1:1Þ
The solutions of the Lane–Emden equation for a given index m are known as polytropes of index m. The closed form solutions are known only for certain integer values of m, the exact solutions for m ¼ 0; 1 and 5 are given in Bender et al. (1989), respectively as follows,
yðtÞ ¼ 1
t2 ; 6
yðtÞ ¼
sinðtÞ ; t
1=2 t2 yðtÞ ¼ 1 þ : 3
In addition, the zeros of yðtÞ were found asymptotically by Bender et al. (1989) as following:
n ¼ p þ 0:885273956d þ 0:24222d2 ; where d ¼ 0:5; 0; 0:5; 1:0 and 1.5 correspond to m ¼ 0; 1; 2; 3 and 4, respectively. The application of the Lane–Emden equation in astrophysics and also its importance in the kinetics of combustion and the Landau–Ginzburg critical phenomena motivates physicists to pay considerable attention to solve it (Dixon and Tuszynski, 1990; Fermi, 1927; Fowler, 1930; Frank-Kamenetskii, 1955; Chandrasekhar,
179
A. Kazemi Nasab et al. / New Astronomy 34 (2015) 178–186
1939; Eddington, 1926; Spitzer, 1942). On the other hand, its nonlinearity and singular behaviour at the origin makes it fascinating for mathematicians to consider it as a prototype for testing new methods for solving nonlinear differential equations. A lot of analytical and numerical approaches have been proposed in literature to solve Lane–Emden equation. Most methods are based on either series solution or perturbation methods. A new perturbation method called d-method was proposed by Bender et al. (1989). Mandelzweig and Tabakin (2001) applied quasilinearization method for solving standard Lane–Emden equation that approximates the solution of a nonlinear differential equation by considering the nonlinear terms as a perturbation about the linear ones. The Adomian decomposition method was employed in Shawagfeh (1993) and Wazwaz (2001). Wazwaz (2006) employed the modified decomposition method for solving the analytical treatment of nonlinear differential equations like the Lane–Emden equation. In Liao (2003), Singh et al. (2009) and Van Gorder and Vajravelu (2008), the homotopy analysis method was proposed for solving Lane–Emden type equations which provides us with a convenient way to adjust convergence regions. Ritz’s method was employed in He (2003) to obtain an analytical solution of the problem. The author of Ramos (2003, 2005, 2008, 2009) proposed different techniques for solving the Lane–Emden equation such as linearization technique, series approach and piecewise-adaptive decomposition methods. The solution of the Lane–Emden equation is numerically challenging because of the singularity behavior at the origin. Yousefi (2006) presented a numerical method based on converting Lane–Emden equations to integral equations and then solve the obtained integral equations using Legendre wavelet approximations. Rational Legendre pseudospectral approach (Parand et al., 2009), Hermite functions collocation method (Parand et al., 2010a), Sinc-Collocation method (Parand and Pirkhedri, 2010), Lagrangian method (Parand et al., 2010b), hybrid functions collocation method (Marzban et al., 2008) and a modified Legendre-spectral method (Adibi and Rismani, 2010) are some other numerical techniques. Furthermore, Karimi Vanani and Aminataei (2010) proposed a numerical method based on converting the Lane–Emden equations into integral equations. They converted the acquired integral equation into a power series and then transformed the power series into Pade series to obtain an approximate polynomial of arbitrary order for solving Lane–Emden equations. Recently, Pandey and Kumar (2012) and Pandey et al. (2012) used Brenstein and Legendre operational matrix of differentiation to solve the problem. A Jacobi–Gauss collocation method proposed in Bhrawy and Alofi (2012). Some more articles concerning Lane–Emden equations can be found in Van Gorder (2011a,b), Boubaker and Van Gorder (2012), Herbst and Momoniat (2012), Doha et al. (2013) and Pandey and Kumar (2012). In recent years, wavelets have received considerable attention by researchers in different fields of science and engineering. One advantage of wavelet analysis is the ability to perform local analysis (Misiti et al., 2000). Wavelet analysis is able to reveal signal aspects that other analysis method miss, such as trends, breakdown points, discontinuities, etc. In comparison with other orthogonal functions, multiresolution analysis aspect of wavelets permit the accurate representation of a variety of functions and operators. In other words, we can change M and k simultaneously to get more accurate solution. Another benefit of wavelet method for solving equations is that after discreting the coefficients matrix of algebraic equations is sparse. So it is computationally efficient the use of wavelet methods for solving equations. In addition, the solution is convergent. In this paper we consider the general Lane–Emden type equations which are formulated as
a
y00 ðtÞ þ y0 ðtÞ þ pðtÞqðyÞ ¼ rðtÞ; t
at P 0
ð1:2Þ
subject to the initial conditions
yð0Þ ¼ A;
y0 ð0Þ ¼ B;
ð1:3Þ
where a; A and B are real constants and pðtÞ; qðyÞ and rðtÞ are some given functions. The well known Lane–Emden type equations with some other particular forms of gðyÞ can be modeled several phenomena in mathematical physics and astrophysics such as the theory of stellar structure, the thermal behavior of a spherical cloud of gas, isothermal gas spheres and the theory of thermionic currents (Chandrasekhar, 1939; Richardson et al., 1921). The paper is organized as follows. Section 2 included some necessary definitions and mathematical preliminaries of Chebyshev polynomials and Chebyshev wavelets. In Section 3, we first introduce the Chebyshev finite difference method and then the Chebyshev wavelet finite difference method (CWFD) is presented. Section 4 included convergence analysis of the proposed method. In Section 5, the proposed approach is used to approximate the Lane–Emden type equations. As a result the Lane–Emden type equation is converted to a system of algebraic equations which is simply solved. Some illustrative examples of different types are given to demonstrate the efficiency and accuracy of the method in Section 6. In Section 7, concluding remarks are given. 2. Preliminaries and notations In this Section, we present some notations, definitions and preliminary facts that will be used further in this work. 2.1. Chebyshev polynomials Chebyshev polynomials of the first kind of degree m, can be defined as follows,
T m ðtÞ ¼ cos mb;
b ¼ arccos t;
ð2:1Þ
which are orthogonal with respect to the weight function pffiffiffiffiffiffiffiffiffiffiffiffiffi wðtÞ ¼ 1= 1 t 2 ,
Z
1
T m ðtÞT n ðtÞwðtÞdt ¼
1
p 2
cm dnm ;
ð2:2Þ
where dnm is the Kronecker delta function and
cm ¼
2;
m ¼ 0;
1;
m=1:
ð2:3Þ
Chebyshev polynomials also satisfy the following recursive formula:
T 0 ðtÞ ¼ 1;
T 1 ðtÞ ¼ t;
T mþ1 ðtÞ ¼ 2tT m ðtÞ T m1 ðtÞ;
m ¼ 1; 2; . . . ð2:4Þ
The set of Chebyshev polynomials is a complete orthogonal set in the Hilbert space ½1; 1. A function f 2 ½1; 1 can be written in terms of Chebyshev polynomials as
f ðtÞ ¼
1 X ^f m T m ðtÞ; m¼0
^f m ¼ 2 pc m
Z
1
f ðtÞT m ðtÞwðtÞdt:
ð2:5Þ
1
2.2. Wavelets and Chebyshev wavelets Wavelets are mathematical tools used to cut up a given function or data into different scale components, and then study each component with a resolution correspond with its scale. They are applicable in a variety of areas, such as engineering, biology, medicine, physics, and mathematics. A wavelet or precisely a mother wavelet is a function w 2 L2 ðRÞ such that the class of functions f2r=2 wð2r t sÞ : r; s 2 Zg is an orthonormal basis for L2 ðRÞ. we have the following family of continuous wavelets as:
180
A. Kazemi Nasab et al. / New Astronomy 34 (2015) 178–186
1 tb ; wa;b ðtÞ ¼ jaj 2 w a
a; b 2 R; a – 0:
ð2:6Þ
Chebyshev wavelets wn;m ¼ wðk; n; m; tÞ, have four arguments, n ¼ 1; . . . ; 2k1 ; k can assume any positive integer, m is degree of Chebyshev polynomials of the first kind and t denotes the time.
(
wn;m ðtÞ ¼
k
22 pm T m ð2k t 2n þ 1Þ;
n1 2k1
0;
n 6 t < 2k1 ;
otherwise;
where the summation symbol with double primes denotes a sum with both the first and last terms halved. Moreover, t m are the extrema of the Mth-order Chebyshev polynomial T M ðtÞ and defined as
tm ¼ cos
mp ; M
tM ¼ 1 < tM1 < < t1 < t0 ¼ 1 are the zeros of ð1
8 1 < pffiffipffi ; m ¼ 0; pm ¼ qffiffiffi : 2; m P 1 p
ð2:8Þ
wn ðtÞ ¼ wð2k t 2n þ 1Þ:
ð2:9Þ
family
of
Chebyshev
wavelets
fm ¼
w1;2 w2;0 w2;1 w2;2
9 ¼ p2ffiffipffi > > = pffiffi 1 2pffiffi2 ffi ¼ p ð4t 1Þ 06t< ; > 2 pffiffi > ; ¼ 2pffiffip2ffi ½2ð4t 1Þ2 1 9 ¼ p2ffiffipffi > > > pffiffi = ¼ 2pffiffip2ffi ð4x 3Þ i> pffiffi h > > ¼ 2pffiffip2ffi 2ð4x 3Þ2 1 ;
The first two derivatives of the function f ðtÞ at the points tm ; m ¼ 0; 1; . . . ; M are given by Elbarbary and El-Kady (2003) as
1 6 t < 1: 2
1 X 1 X cn;m wn;m ðtÞ;
ð2:10Þ
where cn;m ¼ ðf ðtÞ; wn;m ðtÞÞwn , in which ð:; :Þ denotes the inner product in L2 ½0; 1. If we consider truncated series in (2.10), we obtain
cn;m wn;m ðtÞ ¼ C T WðtÞ;
ð1Þ
M k1 X 4hj X khk T k ðt j ÞT l ðtm Þ; M k¼1 l¼0;ðkþlÞodd cl M k1 X 4hj X khk kjp lmp cos cos ¼ M k¼1 l¼0;ðkþlÞodd cl M M
dm;j ¼
2 2 M k2 X 2hj X kðk l Þhk T k ðtj ÞT l ðt m Þ; M k¼2 l¼0;ðkþlÞev en cl 2 2 M k2 X 2hj X kðk l Þhk kjp lmp cos cos ¼ M k¼2 l¼0;ðkþlÞev en cl M M
C ¼ ½c1;0 ; . . . ; c1;M ; c2;0 ; . . . ; c2;M ; . . . ; c2k1 ;1 ; . . . ; c2k1 ;M T ;
3. Chebyshev wavelet finite difference method In this Section, we present the Chebyshev wavelet finite difference method (CWFD). It makes sense to give first a brief introduction to Chebyshev finite difference method. Clenshaw (1960) introduced the following approximation of the function f ðtÞ denoted by ðPM Þf ðtÞ as
m¼0
M
2 X00 f ðt k ÞT m ðt k Þ; M k¼0
1 2k
ðtm þ 2n 1Þ;
ðPM Þf ðtÞ ¼
WðtÞ ¼ ½w1;0 ; . . . ; w1;M ; w2;0 ; . . . ; w2;M ; . . . ; w2k1 ;0 ; . . . ; w2k1 ;M T : ð2:12Þ
fm ¼
ð3:8Þ
ð3:9Þ
then a function f ðtÞ can be written in terms of Chebyshev wavelet basis functions as follows
where C and WðtÞ are 2 ðM þ 1Þ 1 matrices given by
fm T m ðtÞ;
ð3:7Þ
with h0 ¼ hM ¼ 12 ; hj ¼ 1 for j ¼ 1; 2; . . . M 1. As can be seen from Eq. (3.6), the first two derivatives of the function f ðtÞ at any point of the Chebyshev–Gauss–Lobatto points is expanded as a linear combination of the values of the function as these points. Now consider t nm ; n ¼ 1; 2; . . . ; 2k1 T, k is a finite integer number, m ¼ 0; 1; . . . ; M, as the corresponding Chebyshev–Gauss–Lobh n such that atto collocation points at the nthsubinterval 2n1 k1 ; k1 2
ð2:11Þ
k
M X 00
ð3:6Þ
where
tnm ¼
n¼1 m¼0
ðPM Þf ðtÞ ¼
n ¼ 1; 2;
j¼0
ð2Þ
n¼1 m¼0
f ðtÞ ¼
M X ðnÞ dm;j f ðtj Þ;
dm;j ¼
A function f ðtÞ 2 L2 ½0; 1 may be expanded as
2k1 X M X
ð3:5Þ
and
2.3. Function approximation
f ðtÞ ¼
ð3:4Þ
M 2 X00 mkp : f ðtk Þ cos M k¼0 M
fwn;m ðtÞj
n ¼ 1; 2; . . . ; 2k1 ; m 2 N [ f0gg forms an orthonormal basis for L2 ½0; 1 with respect to the weight function wn ðtÞ, see Derili and Sohrabi (2008). For k ¼ 2 and M ¼ 2, Chebyshev wavelets are as following:
w1;1
ð3:3Þ
Using Eq. (3.2) we have,
mkp ; T m ðtk Þ ¼ cos M
f ðnÞ ðt m Þ ¼ 2.1. The
M ðtÞ t 2 Þ dTdt .
so fm can be rewritten as
and m ¼ 0; 1; . . . ; M and n ¼ 1; . . . ; 2k1 . In Eq. (2.8) the coefficients are used for orthonormality. Note that by dealing with the Chebyshev wavelets, to get orthogonal wavelets, the weight func~ tion wðtÞ ¼ wð2t 1Þ need to be dilated and translated as follows
Lemma
ð3:2Þ
The well known Chebyshev–Gauss–Lobatto interpolated points,
ð2:7Þ
where
w1;0
m ¼ 0; 1; 2; . . . ; M:
ð3:1Þ
M 2k1 X X 00
cnm wnm ðtÞ;
ð3:10Þ
n¼1 m¼0
where cnm ; n ¼ 1; 2; . . . ; 2k1 T; m ¼ 0; 1; . . . ; M, are the expansion h n coefficients of the function f ðtÞ at the subinterval 2n1 and k1 ; k1 2 wn;m ðtÞ; n ¼ 1; 2; . . . ; 2k1 T; m ¼ 0; 1; . . . ; M, are defined on ½0; TÞ as in Eq. (2.7) and (2.8). In view of Eqs. (2.5) and (3.1), we can obtain the coefficients cnm as M
2 X00 f ðt np Þwnm ðt np Þ ð2k=2 pm Þ M p¼0 M mpp 1 2 X00 : f ðt np Þ cos ¼ k=2 M M 2 pm p¼0
cnm ¼
1
2
ð3:11Þ
181
A. Kazemi Nasab et al. / New Astronomy 34 (2015) 178–186
Using Eq. (3.7) and (3.8), the first two derivatives of the function f (t) at the points t nm ; n ¼ 1; 2; . . . ; 2k1 T; m ¼ 0; 1; . . . ; M, can be obtained as
f ðrÞ ðtnm Þ ¼
M X ðrÞ dn;m;j f ðt nj Þ;
r ¼ 1; 2;
ð3:12Þ
j¼0
where ð1Þ
M k1 X 4hj X khk w ðt nj Þwnl ðtnm Þ; M k¼1 l¼0;ðkþlÞodd cl pk pl nk M k1 X 4hj X 2k khk kjp lmp cos ; ¼ cos M k¼1 l¼0;ðkþlÞodd cl M M
dn;m;j ¼
Table 2 Comparison of yðtÞ values of standard Lane–Emden equation, for the present method and exact numerical values (Horedt, 1986), for m ¼ 2 with M ¼ 8; k ¼ 4. t
Our method
Exact values
0.5 1 1.5 2 2.5 3 3.5 4
0.9593527 0.8486541 0.6953671 0.5298364 0.3747393 0.2418241 0.133969 0.0488401
0.9593527 0.8486541 0.6953671 0.5298364 0.3747393 0.2418241 0.133969 0.0488402
11=2 #2 M X cnm wn;m ðtÞ wn ðtÞdtA ;
and
ð3:13Þ where
0
and ð2Þ
2 2 M k2 X 2hj X 2k kðk l Þhk wnk ðtnj Þwnl ðt nm Þ; M k¼2 l¼0;ðkþlÞev en cl pk pl 2 2 M k2 X 2hj X 4k kðk l Þhk kjp lmp cos : ð3:14Þ ¼ cos M k¼2 l¼0;ðkþlÞev en cl M M
¼@
rk;M
dn;m;j ¼
Z
1
0
" f ðtÞ
k1 2X T
n¼1 m¼0
pffiffiffiffiffiffiffi 2pB Rp3 þ : 5 12 ð2nÞ2
C¼
5. Discretization of problem 4. Convergence analysis In this section, we investigate the convergence analysis of the proposed method. Interested reader can find a proof of the theorems in Appendix A. Lemma 4.1. If the Chebyshev wavelet expansion of a continuous function f ðtÞ converges uniformly, then the Chebyshev wavelet expansion converges to the function f ðtÞ, see Adibi and Assari (2010) and Folland (1999). Theorem 4.2. A function f ðtÞ 2 L2wn ½0; 1Þ, with bounded second derivative, say jf 00 ðtÞj 6 B, can be expanded as an infinite sum of Chebyshev wavelets, and the series converges uniformly to f ðtÞ, that is, 2k1 X 1 X ^cnm wnm ðtÞ: f ðtÞ ¼
ð4:1Þ
n¼1 m¼0
L2wn ½0; 1Þ
Theorem 4.3. Suppose f ðtÞ 2 with bounded second derivative, say jf 00 ðtÞj 6 B, then its Chebyshev wavelet finite difference expansion converges uniformly to f ðtÞ, that is, 1 2k1 X X 0
cnm wnm ðtÞ ¼ f ðtÞ;
ð4:2Þ
n¼1 m¼0
where the summation symbol with prime denotes a sum with the first term halved.
In this Section, the Chebyshev wavelet finite difference method (CWFD) is employed for solving the Lane–Emden type equations. Consider the following general form:
a
y00 ðtÞ þ y0 ðtÞ þ pðtÞqðyðtÞÞ ¼ rðtÞ; t
a P 0; t 2 ½0; 1
ð5:1Þ
subject to the initial conditions
y0 ð0Þ ¼ B;
yð0Þ ¼ A;
ð5:2Þ
where a; A and B are real constants and pðtÞ; qðyðtÞÞ and rðtÞ are some given functions. Now suppose the interval ½0; 1 is divided into h k1 n 2k1 subintervals In ¼ 2n1 . We also consider k1 ; k1 ; n ¼ 1; 2; . . . ; 2 2 the Chebyshev–Gauss–Lobatto collocation points
n1 2k1
¼ t n0 < t n1 < < tn;M1 < tn;M ¼
n 2k1
on the nth subinterval In ; n ¼ 1; 2; . . . ; 2k1 , where tns is defined as following:
tns ¼
1 2k
ðt s þ 2n 1Þ;
s ¼ 1; 2; . . . ; M 1:
ð5:3Þ
In order to obtain the solution yðtÞ in Eq. (5.1), we first consider yn ðtÞ as restriction of yðtÞ on In and then approximate yn ðtÞ in terms of CWFD basis functions as follows
yn ðtÞ ¼
M X 00
cnm wnm ðtÞ;
ð5:4Þ
m¼0
Theorem 4.4 (Accuracy estimation). Suppose f ðtÞ 2 L2wn ½0; 1Þ with bounded second derivative, say jf 00 ðtÞj 6 B, then we have the following accuracy estimation:
rk;M 6
2k1 X 1 X
!1=2
C2
2 n¼1 m¼Mþ1 ðm 1Þ
2
ð4:3Þ
;
where cnm ; n ¼ 1; 2; . . . ; 2k1 ; m ¼ 0; 1; . . . ; M, are defined in (3.11). We now collocate Eq. (5.1) at the Chebyshev–Gauss–Lobatto points tns ; s ¼ 1; . . . ; M 1, as
y00n ðt ns Þ þ
a tns
y0n ðtns Þ þ pðt ns Þqðyn ðt ns ÞÞ ¼ rðt ns Þ:
ð5:5Þ
In view of Eq. (3.12), we obtain
Table 1 Comparison of the first zeros of standard Lane–Emden equations, for the present method, other methods and exact numerical values (Horedt, 2004). m
Our method
Ramos (2005)
0 1 2 3 4
2.4494897427832 3.1415926535898 4.3528746 6.89684862 14.9715464
2.44899 3.14048 4.35086 6.89312 14.96518
Parand et al. (2010a)
4.3528746 6.89684862 14.9715463
Parand et al. (2010b)
Exact values
4.352875 6.8968486 14.971546
2.4494897427832 3.1415926535898 4.3528746 6.89684862 14.9715463
182
A. Kazemi Nasab et al. / New Astronomy 34 (2015) 178–186
Table 3 Comparison of yðtÞ values of standard Lane–Emden equation, for the present method, other methods and exact numerical values (Horedt, 2004), for m ¼ 3 with M ¼ 10; k ¼ 5. t
Our method
0 0.1 0.5 1 5 6 6.8 6.896
Error
Parand et al. (2010a) 0
1 0.998335830 0.959839070 0.855057569 0.110819835 0.043737984 0.004167789 3.60111E05
0.00 10 2.96 108 3.01 108 8.60 109 3.51 108 3.89 109 3.65 1010 5.45 1011
0
0.00 10 1.40 106 2.99 106 1.99 106 3.90 107 1.14 106 1.05 105 8.88 108
Table 4 Comparison of yðtÞ values of standard Lane–Emden equation, for the present method, method in Parand et al. (2010a) and exact numerical values (Horedt, 2004), for m ¼ 4 with M ¼ 10; k ¼ 5.
Parand and Pirkhedri (2010) 0
Parand et al. (2010b) 0
0.00 10 1.82 105 1.12 104 1.07 104 2.02 104 5.00 105 2.93 104 2.10 105
0.00 10 1.28 105 1.81 105 6.00 107 2.00 107 2.00 105 2.79 106 1.01 106
Exact values 1 0.9983358 0.9598391 0.8550576 0.1108198 0.04373798 0.004167789 3.60112E05
Table 6 Comparison of yðtÞ values of isothermal gas sphere, for the present method, method in Parand et al. (2010a) and series solution given by Wazwaz (2001), for m ¼ 4 with M ¼ 12; k ¼ 2.
t
Our method
Error
Parand et al. (2010a)
Exact values
t
Our method
Error
Parand et al. (2010a)
Wazwaz
0 0.1 0.2 0.5 1 5 10 14 14.9
1 0.99833666 0.993386214 0.960310902 0.860813812 0.235922731 0.059672742 0.008330527 0.000576419
0.00 100 4.05 108 1.35 108 2.36 109 1.22 108 3.09 108 1.59 109 2.79 1010 3.46 1012
0.00 100 2.51 104 2.48 104 2.05 104 1.93 104 8.59 105 6.22 105 2.47 105 5.00 107
1 0.9983367 0.9933862 0.9603109 0.8608138 0.2359227 0.05967274 0.008330527 0.000576419
0 0.1 0.2 0.5 1 1.5 2 2.5
0 0.001666419 0.006653971 0.041154515 0.158828174 0.338019831 0.559823312 0.806341085
0.00 100 3.81 1011 3.98 1013 4.92 1010 3.24 107 6.31 106 1.40 104 3.67 103
0.00 100 5.85 107 6.04 107 5.58 107 8.20 107 6.72 106 1.39 104 3.68 103
0 0.001665834 0.006653367 0.041153957 0.158827354 0.33801311 0.55996266 0.810019671
Table 5 ~T for m ¼ 5. The maximum absolute errors ET and E L 50 30 10 5
Et
~ Et 7
1.4 10 2.0 1010 2.0 1014 1.8 1017
3.0 107 1.6 109 1.0 1013 1.0 1016
Fig. 2. Graph of solutions to the Lane–Emden equation modelling an isothermal gas sphere in comparing the current method and Wazwaz solution (Wazwaz, 2001). Table 7 Comparison of yðtÞ values in Example 6.3, for the present method, method in Parand et al. (2010a) and exact solution. t
Fig. 1. Graph of solutions to the standard Lane–Emden equation for m ¼ 0; 1; 2; 3 and 4 respectively.
0 0.01 0.1 0.5 1 2 3 4 5 6 7 8 9 10
Our method 1 0.00019999 0.019900662 0.446287103 1.386294361 3.218875825 4.605170186 5.666426688 6.516193076 7.221835825 7.824046011 8.34877454 8.813438495 9.230241034
Error
Parand et al. (2010a) 0
0.00 10 2.80 1015 2.30 1015 2.55 1015 1.01 1013 1.07 1013 5.79 1014 1.96 1014 6.06 1015 2.28 1014 3.36 1014 4.05 1014 4.47 1014 3.96 1011
0
0.00 10 2.93 106 3.93 106 3.02 106 9.31 107 5.00 107 8.10 107 7.69 107 6.64 107 5.48 107 1.70 107 1.09 106 1.21 105 3.83 105
Exact values 1 0.00019999 0.019900662 0.446287103 1.386294361 3.218875825 4.605170186 5.666426688 6.516193076 7.221835825 7.824046011 8.34877454 8.813438495 9.230241034
183
A. Kazemi Nasab et al. / New Astronomy 34 (2015) 178–186
Fig. 3. (a) The absolute error between approximate and exact solution, (b) the absolute error between the first derivative of approximate and exact solution for Example 6.3.
yðrÞ n ðt ns Þ ¼
M X ðrÞ dn;s;j yn ðtnj Þ;
r ¼ 1; 2;
ð5:6Þ
~T ¼ maxfjy0 ðtÞ y0 ðtÞj : 0 6 t 6 Lg; E e a
j¼0 ðrÞ dn;s;j ;
where r ¼ 1; 2 are defined in (3.13) and (3.14). In addition, substituting Eqs. (3.10) and (3.12) into the initial conditions (5.2), we get other two equations. Besides, it is also necessary that the approximate solution and it’s first derivative be continuous at the interface of subintervals, then
yðrÞ n
n
2k1
¼
ðrÞ
lim yn1 ðtÞ:
t!
ð5:7Þ
n 2k1
Therefore, we have a system of ðM þ 1Þ algebraic equations, which can be solved for the unknowns yn ðt ns Þ. Consequently, we obtain the solution yn ðtÞ to the given Eqs. (5.1) and (5.2). Note that for solving the given equation at the interval I1 , we make use of the given initial conditions and for other subintervals, we use Eq. (5.7) to get initial conditions. After that we obtain yn ðtÞ and the summation will be the required solution yt . 6. Illustrative examples The Lane–Emden equation can be analytically solved only for a few special, integer values of certain index. For other values of index, the numerical solutions can be obtained. In this Section, we consider some numerical examples to demonstrate the validity of the method (CWFD) that we propose to solve the Lane–Emden type equations. The reason that we consider these examples there are closed form solutions that available in the literature. This allows one to compare the results which were obtained by using this method and the analytical solution as well as with the solutions obtained by using other methods. Example 6.1. As the first example, we consider the standard Lane– Emden equation that was used to model the thermal behavior of a spherical cloud of gas acting under the mutual attraction of its molecules. If we set pðtÞ ¼ 1; qðyÞ ¼ ym ; rðtÞ ¼ 0; a ¼ 2; A ¼ 1 and B ¼ 0 in Eq. (1.2), we get
2 y00 ðtÞ þ y0 ðtÞ þ yðtÞm ¼ 0; t
tP0
subject to the boundary conditions
yð0Þ ¼ 1;
and
y0 ð0Þ ¼ 0:
Consider
ET ¼ maxfjye ðtÞ ya ðtÞj : 0 6 t 6 Lg
ð6:1Þ
where ya ðtÞ and ye ðtÞ denote the approximate and the exact solutions, respectively. We applied the method introduced in Section 6 and solved this problem for M ¼ 10; k ¼ 5. In Table 1 we compare the first zeros of standard Lane–Emden equations obtained using the present method with M ¼ 10; k ¼ 5, those reported in literature resulted from other well-known methods and exact values given by Horedt (2004) for m ¼ 0; 1; 2; 3 and 4, respectively. As can be seen from Table 1, the present method provides very accurate predictions of the zeros of yðtÞ even in large intervals. Tables 2–4 show the approximation of yðtÞ for the standard Lane–Emden equation for m ¼ 2; 3 and 4 respectively obtained by the proposed method and those reported in Horedt (2004, 1986). In addition, Table 5 shows the maximum absolute errors for m ¼ 5 and different values of L, for which the exact solution exists. Our obtained results confirm the accuracy of the present method even for large values of L. The resulting graph for m ¼ 0; 1; 2; 3 and 4 is shown in Fig. 1. Next is a display of the solutions for several values of m. Example 6.2. As the second example, we consider the isothermal gas spheres equation with the following formula:
2 y00 ðtÞ þ y0 ðtÞ þ eyðtÞ ¼ 0; t
tP0
ð6:2Þ
subject to the boundary conditions
yð0Þ ¼ 0;
y0 ð0Þ ¼ 0:
Wazwaz (2001), Liao (2003), Singh et al. (2009) and Ramos (2008) obtained a series solution for this equation using ADM, MHAM and series expansion respectively:
1 1 4 8 6 122 8 61:67 10 yðtÞ ffi x2 þ x x þ x x : 6 5:4! 21:6! 81:8! 495:10! This equation has been also solved using Hermite functions collocation, HPM, series solutions, ADM and HAM methods respectively in Parand et al. (2010a), Chowdhury and Hashim (2009), Aslanov (2008), Aslanov (2008) and Bataineh et al. (2009). We apply Chebyshev wavelet finite difference method to solve the equation with M ¼ 12; k ¼ 2. In Table 6 we make a comparison among the values of yðtÞ obtained by the present method and those obtained in Wazwaz (2001) and Parand et al. (2010a) which
184
A. Kazemi Nasab et al. / New Astronomy 34 (2015) 178–186
indicate that our results are more accurate than those reported by Parand et al. (2010a). Fig. 2 shows the resulting graph of the isothermal gas spheres equation in comparison to the current method and that obtained by Wazwaz (2001). Example 6.3. If we set pðtÞ ¼ 1; qðyðtÞÞ ¼ 4 2eyðtÞ þ eyðtÞ=2 ; rðtÞ ¼ 0; A ¼ 0 and B ¼ 0 in Eq. (5.1), we will get one of the Lane–Emden type equations with the following formula:
2 y00 ðtÞ þ y0 ðtÞ þ 4 2eyðtÞ þ eyðtÞ=2 ¼ 0; t
tP0
ð6:3Þ
subject to the boundary conditions
yð0Þ ¼ 0;
y0 ð0Þ ¼ 0;
which has the analytical solution yðtÞ ¼ 2 lnð1 þ t 2 Þ. This equation has been solved in Parand et al. (2010a), Chowdhury and Hashim (2009) and Yildirim and Ozis (2009) with Hermite functions collocation, HPM and VIM method. We apply the introduced method with M ¼ 15; k ¼ 5 for solving Eq. (6.3). In Table 7, we compare the values of yðtÞ obtained by the present method, those reported by Parand et al. (2010a) and exact values. The absolute error between approximate and exact solution, also the absolute error between the first derivative of approximate and exact solution is shown in Fig. 3.
^cnm ¼
rffiffiffiffi Z p cos a þ 2n 1 2 cos ma da f k k p 2 2 2 0 1
using integration by part, we get
pffiffiffi pffiffiffi #p 2 2 cos a þ 2n 1 sin ma f þ pffiffiffiffi k pffiffiffiffi 3k m 2m 2k p 22 p 2 0 Z p cos a þ 2n 1 sin ma sin a da; f0 2k 0
^cnm ¼
the first part is zero, therefore,
^cnm ¼
pffiffiffi Z p 2 0 cos a þ 2n 1 sin ma sin a da: f pffiffiffiffi 3k 2k 22 m p 0
Using integration by part again, it yields
#p 1 sinðm 1Þa sinðm þ 1Þa 0 cos a þ 2n 1 f pffiffiffiffiffiffiffi 3k m1 mþ1 2k 2 2 m 2p 0 Z p 1 00 cos a þ 2n 1 r m ðaÞ da; þ 5k pffiffiffiffiffiffiffi f 2k 2 2 m 2p 0
^cnm ¼
where
sinðm 1Þa sinðm þ 1Þa : rm ðaÞ ¼ sin a m1 mþ1 Thus, we get
7. Conclusion An efficient and accurate method based on hybrid of Chebyshev wavelets and finite difference methods was introduced. The useful properties of Chebyshev wavelets and finite difference method make it a computationally efficient method to approximate the solution of nonlinear Lane–Emden type equations in a semi-infinite interval. We converted the given problem to a system of algebraic equations using collocation points. The main advantage of the present method is the ability to represent smooth and especially piecewise smooth functions properly. It was also shown that the accuracy can be improved either by increasing the number of subintervals or by increasing the number of collocation points in subintervals. Moreover, this method is valid for large-domain calculations. Several examples have been provided to demonstrate the powerfulness of the proposed method. A comparison was made among the present method, some other well-known approaches and exact solution which confirms that the introduced method are more accurate and efficient.
Z p
1
00 cos a þ 2n 1 ^ r m ðaÞ da
jcnm j ¼ 5k pffiffiffiffiffiffiffi f k
2 2 m 2p 0 2 !Z
p
1
f 00 cos a þ 2n 1 rm ðaÞ da 6 pffiffiffiffiffiffiffi 5k
k 2 2 2 m 2p 0 Z p B 6 5k pffiffiffiffiffiffiffi jrm ðaÞj da: 2 2 m 2p 0
Z p
sin a sinðm 1Þa sinðm þ 1Þa da
m1 mþ1 0
Z p
sin a sinðm 1Þa sin a sinðm þ 1Þa
þ
da 6
m1 mþ1 0
0
6
2mp ; m2 1
m > 1:
Since n 6 2k1 , then we obtain
pffiffiffiffiffiffiffi 2pB
The authors are very grateful to the referees for their valuable suggestions and comments that improved the paper.
Now, in the case of m ¼ 1, apply Eq. (7.1), and we have
j^cn1 j <
Here, we include proofs of the theorems, see Kazemi Nasab et al. (2013). Proof of Theorem 4.2. We have
¼
Z
n 2k1
n1 2k1
n;k
¼
Z
1
f ðtÞwnm ðtÞwn ðtÞdt 0
k
22 pm f ðtÞT m ð2k t 2n þ 1Þwð2k t 2n þ 1Þdt;
if m > 1, by substituting 2k t 2n þ 1 ¼ cos a, it yields
ð7:3Þ
jrm ðaÞj da ¼
Acknowledgments
Appendix A
ð7:2Þ
However
Z p
j^cnm j 6
^cn;m ¼ ðf ðtÞ; wn;m ðtÞÞ w
ð7:1Þ
5
ð2nÞ2 ðm2 1Þ
:
pffiffiffi pp ffiffi 2
max jf 0 ðtÞj: 3 ð2nÞ2 06t61
ð7:4Þ
ð7:5Þ
It is also worth to mention that the set fwn0 g1 n¼1 form an orthogonal system which was constructed by Haar scaling function with P 2kffiffiffi ^ p respect to the weight function wðtÞ, so 1 n¼1 p c n0 wn0 ðtÞ is convergent, see Daubechies et al. (1992). Hence, we will have
k1
k1
X
X
X 2 1 2 2k1 X 1
^
X^ c c w ðtÞ 6 w ðtÞ þ j^c jjw ðtÞj
nm nm n0 n0
n¼1
n¼1 m¼1 nm nm
n¼1 m¼0
k1
X
X 2 2k1 X 1
j^c j < 1 6 ^cn0 wn0 ðtÞ þ
n¼1 m¼1 nm
n¼1 P P1 ^ Therefore, with the aid of Lemma 4.1, the series 1 n¼1 m¼1 c nm wnm ðtÞ converges to f ðtÞ uniformly, see Adibi and Assari (2010). h
185
A. Kazemi Nasab et al. / New Astronomy 34 (2015) 178–186
Proof of Theorem 4.3. From Theorem 4.2, we have
Proof of Theorem 4.4. Since
k1
2 1 X X ^cnm wnm ðtÞ;
f ðtÞ ¼
2 k;M
r
n¼1 m¼0
¼
¼
n;k
¼
Z
n 2k1
n1 2k1
¼
Z
f ðtÞwnm ðtÞwn;k ðtÞdt;
¼
0
k 2
2 pm f ðtÞT m ð2k t 2n þ 1Þwð2k t 2n þ 1Þdt;
ð7:6Þ
Z p cos a þ 2n 1 cos ma da: f k 2k 22 0
aj ¼
jp : M
f ðt nj Þ cos
j¼0
Rp3 12M 2
mjp ¼ cnm : M
where
(
) 2 d cosðaÞ þ 2n 1 cosðmaÞ ; 0 6 a 6 p : R ¼ max f da2 2k In view of Eq. (7.1) and by triangle inequality, we get
6
Rp3
12M pffiffiffiffiffiffiffi 2pB
þ 2
pffiffiffiffiffiffiffi 2pB 5
ð2nÞ2 ðm2 1Þ
Rp3 C 6 ; þ 12ðm2 1Þ ð2nÞ52 ðm2 1Þ ðm2 1Þ
ð7:7Þ
where
C¼
pffiffiffiffiffiffiffi Rp3 2p B þ : 5 12 ð2nÞ2
Since jwnm ðtÞj 6 1, then it follows that
k1
k1
k1 1
X
X
X 2 2 2 1
X0 c w ðtÞ 6
c w ðtÞ þ c w ðtÞ
n¼1 m¼0 nm nm n¼1 2 n0 n0 n¼1 n1 n1
þ
2k1 X 1 X
k1
jcnm jjwnm ðtÞj 6
n¼1 m¼2
2 1X jcn0 j 2 n¼1
2k1 2k1 X 1 X X þ jcn1 j þ jcnm j < 1: n¼1
n¼1 m¼0
1 X
c2nm ðwn;m ðtÞÞ2 wn ðtÞdt
n¼1 m¼Mþ1
2k1 X 1 X
c2nm
Z
1
ðwn;m ðtÞÞ2 wn ðtÞdt:
0
2k1 X 1 X
c2nm 6
2k1 X 1 X n¼1 m¼Mþ1
C2 ðm2
2
1Þ
;
where
pffiffiffiffiffiffiffi Rp3 2pB þ : 5 12 ð2nÞ2
References
;
jcnm j 6 jcnm ^cnm j þ j^cnm j 6
n¼1 m¼0 2k1 1X
That completes the proof. h
According to approximation error for the trapezoidal rule, we have
j^cnm cnm j 6
Z
n¼1 m¼Mþ1
C¼
From Eqs. (3.9) and (3.11), for 1 < m < M, we have
k 22 p m M
n¼1 m¼0
" k1 #2 2 1 2k1 X M 1 X X X cnm wn;m ðtÞ cnm wn;m ðtÞ wn;k ðtÞdt
We know that the family fwn;m ðtÞj n 2 N; m 2 N [ f0gg forms orthoR1 normal basis for L2wn ½0; 1Þ, so 0 ðwn;m ðtÞÞ2 wn ðtÞdt ¼ 1. Therefore, in view of Eq. (7.7), we will have
r2k;M ¼
M p p X00 cosðaj Þ þ 2n 1 cos maj ; ’ mk f k M 2 2 2 j¼0
M 00 p 2X
#2 2k1 X M X f ðtÞ cnm wn;m ðtÞ wn ðtÞdt
n¼1 m¼Mþ1
pm
Using the trapezoidal rule for integration with M equidistant subintervals gives
^cnm ’
¼
1
by substituting 2 t 2n þ 1 ¼ cos a, it yields
^cnm
Z
0
k
^cnm ¼
"
0
We prove that this series converges to f ðtÞ uniformly. First of all we show that cnm converges to ^cnm . We have
^cn;m ¼ ðf ðtÞ; wn;m ðtÞÞ w
1
0
where
^cn;m ¼ ðf ðtÞ; wn;m ðtÞÞ : wn
Z
ð7:8Þ
n¼1 m¼2
Therefore, in view of Lemma 4.1, series (4.2) is uniformly convergent to f ðtÞ. h
Adibi, H., Assari, P., 2010. Math. Probl. Eng. 2010. Article ID 138408. Adibi, H., Rismani, A.M., 2010. Comput. Math. Appl. 60, 2126. Aslanov, A., 2008. Phys. Lett. A 372, 3555. Aslanov, A., 2008. Int. J. Comput. Math. 85, 1709. Bataineh, A.S., Noorani, M.S.M., Hashim, I., 2009. Commun. Nonlinear Sci. Numer. Simul. 14, 1121. Bender, C.M., Milton, K.A., Pinsky, S.S., Simmons Jr., L.M., 1989. J. Math. Phys. 30, 1447. Bhrawy, A.H., Alofi, A.S., 2012. Commun. Nonlinear Sci. Numer. Simulat. 17, 62. Boubaker, K., Van Gorder, Robert A., 2012. New Astron. 17, 565. Chandrasekhar, S., 1939. An Introduction to The Study of Stellar Structure. Dover, New York. Chowdhury, M.S.H., Hashim, I., 2009. Nonlinear Anal. Real World Appl. 10, 104. Clenshaw, C.W., Curtis, A.R., 1960. Numer. Math. 2, 197. Daubechies, I., 1992. Ten Lectures on Wavelets. CBMS-NSF Regional Conference Series in Applied Mathematics, vol. 61. SIAM, Philadelphia, Pa, USA. Derili, H., Sohrabi, S., 2008. Math. Sci. (QJMS) 2 (3), 261. Dixon, J.M., Tuszynski, J.A., 1990. Phys. Rev. A 41, 4166. Doha, E.H., Abd- Elhameed, W.M., Youssri, Y.H., 2013. New Astron. 23–24, 113. Eddington, A.S., 1926. The Internal Constitution of The Stars. Cambridge University Press, London. Elbarbary, Elsayed M.E., El-Kady, M., 2003. Appl. Math. Comput. 139, 513. Emden, R., 1907. Gaskugeln Anwendungen der Mechan. Warmtheorie, Druck und Verlag Von B.G. Teubner, Leipzig and Berlin. Fermi, E., 1927. Rend. Accad. Naz. Lincei 6, 602. Folland, G.B., 1999. Real Analysis: Modern Techniques and Their Applications, Pure and Applied Mathematics, second ed. John Wiley & Sons, New York, NY, USA. Fowler, R.H., 1930. MNRAS 91, 63. Frank-Kamenetskii, D.A., 1955. Diffusion and Heat Exchange in Chemical Kinetics. Princeton University Press, Princeton. He, J.H., 2003. Appl. Math. Comput. 143, 539. Herbst, R.S., Momoniat, E., 2012. New Astron. 17, 388. Horedt, G.P., 1986. Astrophys. Space Sci. 126 (2), 357. Horedt, G.P., 2004. Polytropes: Applications in Astrophysics and Related Fields. Kluwer Academic Publishers, Dordrecht. Karimi Vanani, S., Aminataei, A., 2010. Comput. Math. Appl. 59, 2815. Kazemi Nasab, A., Kilicman, A., Pashazadeh Atabakan, Z., Abbasbandy, S., 2013. Abstr. Appl. Anal. [Article ID: 916456]. Lane, J.H., 1870. Am. J. Sci. Arts 50, 57. Liao, S., 2003. Appl. Math. Comput. 142, 1. Mandelzweig, V.B., Tabakin, F., 2001. Comput. Phys. Commun. 141, 268. Marzban, H.R., Tabrizidooz, H.R., Razzaghi, M., 2008. Phys. Lett. A 372 (37), 5883. Misiti, M., Misiti, Y., Oppenheim, G., Poggi, J.-M., 2000. Wavelets Toolbox Users Guide. The MathWorks, Wavelet Toolbox, for use with Matlab. Pandey, Rajesh K., Kumar, Narayan, 2012. New Astron. 17, 303. Pandey, Rajesh K., Kumar, Narayan, 2012. New Astron. 17, 303.
186
A. Kazemi Nasab et al. / New Astronomy 34 (2015) 178–186
Pandey, Rajesh K., Kumar, Narayan, Bhardwaj, Abhinav, Dutta, Goutam, 2012. Appl. Math. Comput. 218, 7629. Parand, K., Pirkhedri, A., 2010. New Astron. 15, 533. Parand, K., Shahini, M., Dehghan, M., 2009. J. Comput. Phys. 228, 8830. Parand, K., Dehghan, M., Rezaeia, A.R., Ghaderi, S.M., 2010a. Comput. Phys. Commun. 181, 1096. Parand, K., Rezaei, A.R., Taghavi, A., 2010b. Acta Astronaut. 67, 673. Ramos, J.I., 2003. Comput. Phys. Commun. 153, 199. Ramos, J.I., 2005. Appl. Math. Comput. 161, 525. Ramos, J.I., 2008. Chaos Solitons Fract. 38, 400. Ramos, J.I., 2009. Chaos Solitons Fract. 40, 1623. Richardson, O.W., 1921. The Emission of Electricity from Hot Bodies. Longmans, Green and co., London, New York.
Shawagfeh, N.T., 1993. J. Math. Phys. 34, 4364. Singh, O.P., Pandey, R.K., Singh, V.K., 2009. Comput. Phys. Commun. 180, 1116. Spitzer, L., 1942. Ap. J. 95, 329. Van Gorder, Robert A., 2011a. New Astron. 16, 65. Van Gorder, Robert A., 2011b. New Astron. 16, 492. Van Gorder, R.A., Vajravelu, K., 2008. Phys. Lett. A 372, 6060. Wazwaz, A., 2001. Appl. Math. Comput. 118, 287. Wazwaz, A., 2006. Appl. Math. Comput. 173, 165. Yildirim, A., Ozis, T., 2009. Ser. A Theor. Methods Appl. 70, 2480. Yousefi, S.A., 2006. Appl. Math. Comput. 181, 1417.
Contents lists available at ScienceDirect
New Astronomy journal homepage: www.elsevier.com/locate/newast
A numerical approach for solving singular nonlinear Lane–Emden type equations arising in astrophysics A. Kazemi Nasab, A. Kılıçman ⇑, Z. Pashazadeh Atabakan, W.J. Leong Department of Mathematics, University Putra Malaysia (UPM), 43400 Serdang, Selangor, Malaysia
h i g h l i g h t s To present an efficient method to solve Lane–Emden type equations. To represent smooth and especially piecewise smooth functions properly. To improve the accuracy by properly increasing either the number of subintervals or the number of collocation points. To compare present solutions for more accuracy and efficiency.
a r t i c l e
i n f o
a b s t r a c t
Article history: Received 16 May 2013 Received in revised form 28 March 2014 Accepted 17 June 2014 Available online 11 July 2014 Communicated by W. Soon Keywords: Wavelet analysis method Chebyshev finite difference method Lane–Emden type equations Singular IVPs Isothermal gas spheres
In this paper, we suggest a numerical method based upon hybrid of Chebyshev wavelets and finite difference methods for solving well-known nonlinear initial-value problems of Lane–Emden type. The useful properties of the Chebyshev wavelets and finite difference method are utilized to reduce the computation of the problem to a set of nonlinear algebraic equations. Making a comparison among the obtained results using the present method with those ones reported in literature by some other well-known methods confirms the accuracy and computational efficiency of the present technique. Ó 2014 Elsevier B.V. All rights reserved.
1. Introduction In recent years, solving singular initial value problems formulated by second order nonlinear differential equations has gained a great deal of attention from many scientists in mathematics and physics. Among equations of this type, the Lane–Emden type equations are well-known and are used as the benchmark to verify new methods. The Lane–Emden differential equation of index m was first introduced by Lane in 1870 (Lane, 1870) and then studied more in detail by Emden (1907), describing the equilibria of nonrotating polytropic fluids in a self-gravitating star and formulated as following:
(
y00 ðtÞ þ 2t y0 ðtÞ þ yðtÞm ¼ 0; yð0Þ ¼ 1;
tP0
y0 ð0Þ ¼ 0:
⇑ Corresponding author. E-mail address: [email protected] (A. Kılıçman). http://dx.doi.org/10.1016/j.newast.2014.06.008 1384-1076/Ó 2014 Elsevier B.V. All rights reserved.
ð1:1Þ
The solutions of the Lane–Emden equation for a given index m are known as polytropes of index m. The closed form solutions are known only for certain integer values of m, the exact solutions for m ¼ 0; 1 and 5 are given in Bender et al. (1989), respectively as follows,
yðtÞ ¼ 1
t2 ; 6
yðtÞ ¼
sinðtÞ ; t
1=2 t2 yðtÞ ¼ 1 þ : 3
In addition, the zeros of yðtÞ were found asymptotically by Bender et al. (1989) as following:
n ¼ p þ 0:885273956d þ 0:24222d2 ; where d ¼ 0:5; 0; 0:5; 1:0 and 1.5 correspond to m ¼ 0; 1; 2; 3 and 4, respectively. The application of the Lane–Emden equation in astrophysics and also its importance in the kinetics of combustion and the Landau–Ginzburg critical phenomena motivates physicists to pay considerable attention to solve it (Dixon and Tuszynski, 1990; Fermi, 1927; Fowler, 1930; Frank-Kamenetskii, 1955; Chandrasekhar,
179
A. Kazemi Nasab et al. / New Astronomy 34 (2015) 178–186
1939; Eddington, 1926; Spitzer, 1942). On the other hand, its nonlinearity and singular behaviour at the origin makes it fascinating for mathematicians to consider it as a prototype for testing new methods for solving nonlinear differential equations. A lot of analytical and numerical approaches have been proposed in literature to solve Lane–Emden equation. Most methods are based on either series solution or perturbation methods. A new perturbation method called d-method was proposed by Bender et al. (1989). Mandelzweig and Tabakin (2001) applied quasilinearization method for solving standard Lane–Emden equation that approximates the solution of a nonlinear differential equation by considering the nonlinear terms as a perturbation about the linear ones. The Adomian decomposition method was employed in Shawagfeh (1993) and Wazwaz (2001). Wazwaz (2006) employed the modified decomposition method for solving the analytical treatment of nonlinear differential equations like the Lane–Emden equation. In Liao (2003), Singh et al. (2009) and Van Gorder and Vajravelu (2008), the homotopy analysis method was proposed for solving Lane–Emden type equations which provides us with a convenient way to adjust convergence regions. Ritz’s method was employed in He (2003) to obtain an analytical solution of the problem. The author of Ramos (2003, 2005, 2008, 2009) proposed different techniques for solving the Lane–Emden equation such as linearization technique, series approach and piecewise-adaptive decomposition methods. The solution of the Lane–Emden equation is numerically challenging because of the singularity behavior at the origin. Yousefi (2006) presented a numerical method based on converting Lane–Emden equations to integral equations and then solve the obtained integral equations using Legendre wavelet approximations. Rational Legendre pseudospectral approach (Parand et al., 2009), Hermite functions collocation method (Parand et al., 2010a), Sinc-Collocation method (Parand and Pirkhedri, 2010), Lagrangian method (Parand et al., 2010b), hybrid functions collocation method (Marzban et al., 2008) and a modified Legendre-spectral method (Adibi and Rismani, 2010) are some other numerical techniques. Furthermore, Karimi Vanani and Aminataei (2010) proposed a numerical method based on converting the Lane–Emden equations into integral equations. They converted the acquired integral equation into a power series and then transformed the power series into Pade series to obtain an approximate polynomial of arbitrary order for solving Lane–Emden equations. Recently, Pandey and Kumar (2012) and Pandey et al. (2012) used Brenstein and Legendre operational matrix of differentiation to solve the problem. A Jacobi–Gauss collocation method proposed in Bhrawy and Alofi (2012). Some more articles concerning Lane–Emden equations can be found in Van Gorder (2011a,b), Boubaker and Van Gorder (2012), Herbst and Momoniat (2012), Doha et al. (2013) and Pandey and Kumar (2012). In recent years, wavelets have received considerable attention by researchers in different fields of science and engineering. One advantage of wavelet analysis is the ability to perform local analysis (Misiti et al., 2000). Wavelet analysis is able to reveal signal aspects that other analysis method miss, such as trends, breakdown points, discontinuities, etc. In comparison with other orthogonal functions, multiresolution analysis aspect of wavelets permit the accurate representation of a variety of functions and operators. In other words, we can change M and k simultaneously to get more accurate solution. Another benefit of wavelet method for solving equations is that after discreting the coefficients matrix of algebraic equations is sparse. So it is computationally efficient the use of wavelet methods for solving equations. In addition, the solution is convergent. In this paper we consider the general Lane–Emden type equations which are formulated as
a
y00 ðtÞ þ y0 ðtÞ þ pðtÞqðyÞ ¼ rðtÞ; t
at P 0
ð1:2Þ
subject to the initial conditions
yð0Þ ¼ A;
y0 ð0Þ ¼ B;
ð1:3Þ
where a; A and B are real constants and pðtÞ; qðyÞ and rðtÞ are some given functions. The well known Lane–Emden type equations with some other particular forms of gðyÞ can be modeled several phenomena in mathematical physics and astrophysics such as the theory of stellar structure, the thermal behavior of a spherical cloud of gas, isothermal gas spheres and the theory of thermionic currents (Chandrasekhar, 1939; Richardson et al., 1921). The paper is organized as follows. Section 2 included some necessary definitions and mathematical preliminaries of Chebyshev polynomials and Chebyshev wavelets. In Section 3, we first introduce the Chebyshev finite difference method and then the Chebyshev wavelet finite difference method (CWFD) is presented. Section 4 included convergence analysis of the proposed method. In Section 5, the proposed approach is used to approximate the Lane–Emden type equations. As a result the Lane–Emden type equation is converted to a system of algebraic equations which is simply solved. Some illustrative examples of different types are given to demonstrate the efficiency and accuracy of the method in Section 6. In Section 7, concluding remarks are given. 2. Preliminaries and notations In this Section, we present some notations, definitions and preliminary facts that will be used further in this work. 2.1. Chebyshev polynomials Chebyshev polynomials of the first kind of degree m, can be defined as follows,
T m ðtÞ ¼ cos mb;
b ¼ arccos t;
ð2:1Þ
which are orthogonal with respect to the weight function pffiffiffiffiffiffiffiffiffiffiffiffiffi wðtÞ ¼ 1= 1 t 2 ,
Z
1
T m ðtÞT n ðtÞwðtÞdt ¼
1
p 2
cm dnm ;
ð2:2Þ
where dnm is the Kronecker delta function and
cm ¼
2;
m ¼ 0;
1;
m=1:
ð2:3Þ
Chebyshev polynomials also satisfy the following recursive formula:
T 0 ðtÞ ¼ 1;
T 1 ðtÞ ¼ t;
T mþ1 ðtÞ ¼ 2tT m ðtÞ T m1 ðtÞ;
m ¼ 1; 2; . . . ð2:4Þ
The set of Chebyshev polynomials is a complete orthogonal set in the Hilbert space ½1; 1. A function f 2 ½1; 1 can be written in terms of Chebyshev polynomials as
f ðtÞ ¼
1 X ^f m T m ðtÞ; m¼0
^f m ¼ 2 pc m
Z
1
f ðtÞT m ðtÞwðtÞdt:
ð2:5Þ
1
2.2. Wavelets and Chebyshev wavelets Wavelets are mathematical tools used to cut up a given function or data into different scale components, and then study each component with a resolution correspond with its scale. They are applicable in a variety of areas, such as engineering, biology, medicine, physics, and mathematics. A wavelet or precisely a mother wavelet is a function w 2 L2 ðRÞ such that the class of functions f2r=2 wð2r t sÞ : r; s 2 Zg is an orthonormal basis for L2 ðRÞ. we have the following family of continuous wavelets as:
180
A. Kazemi Nasab et al. / New Astronomy 34 (2015) 178–186
1 tb ; wa;b ðtÞ ¼ jaj 2 w a
a; b 2 R; a – 0:
ð2:6Þ
Chebyshev wavelets wn;m ¼ wðk; n; m; tÞ, have four arguments, n ¼ 1; . . . ; 2k1 ; k can assume any positive integer, m is degree of Chebyshev polynomials of the first kind and t denotes the time.
(
wn;m ðtÞ ¼
k
22 pm T m ð2k t 2n þ 1Þ;
n1 2k1
0;
n 6 t < 2k1 ;
otherwise;
where the summation symbol with double primes denotes a sum with both the first and last terms halved. Moreover, t m are the extrema of the Mth-order Chebyshev polynomial T M ðtÞ and defined as
tm ¼ cos
mp ; M
tM ¼ 1 < tM1 < < t1 < t0 ¼ 1 are the zeros of ð1
8 1 < pffiffipffi ; m ¼ 0; pm ¼ qffiffiffi : 2; m P 1 p
ð2:8Þ
wn ðtÞ ¼ wð2k t 2n þ 1Þ:
ð2:9Þ
family
of
Chebyshev
wavelets
fm ¼
w1;2 w2;0 w2;1 w2;2
9 ¼ p2ffiffipffi > > = pffiffi 1 2pffiffi2 ffi ¼ p ð4t 1Þ 06t< ; > 2 pffiffi > ; ¼ 2pffiffip2ffi ½2ð4t 1Þ2 1 9 ¼ p2ffiffipffi > > > pffiffi = ¼ 2pffiffip2ffi ð4x 3Þ i> pffiffi h > > ¼ 2pffiffip2ffi 2ð4x 3Þ2 1 ;
The first two derivatives of the function f ðtÞ at the points tm ; m ¼ 0; 1; . . . ; M are given by Elbarbary and El-Kady (2003) as
1 6 t < 1: 2
1 X 1 X cn;m wn;m ðtÞ;
ð2:10Þ
where cn;m ¼ ðf ðtÞ; wn;m ðtÞÞwn , in which ð:; :Þ denotes the inner product in L2 ½0; 1. If we consider truncated series in (2.10), we obtain
cn;m wn;m ðtÞ ¼ C T WðtÞ;
ð1Þ
M k1 X 4hj X khk T k ðt j ÞT l ðtm Þ; M k¼1 l¼0;ðkþlÞodd cl M k1 X 4hj X khk kjp lmp cos cos ¼ M k¼1 l¼0;ðkþlÞodd cl M M
dm;j ¼
2 2 M k2 X 2hj X kðk l Þhk T k ðtj ÞT l ðt m Þ; M k¼2 l¼0;ðkþlÞev en cl 2 2 M k2 X 2hj X kðk l Þhk kjp lmp cos cos ¼ M k¼2 l¼0;ðkþlÞev en cl M M
C ¼ ½c1;0 ; . . . ; c1;M ; c2;0 ; . . . ; c2;M ; . . . ; c2k1 ;1 ; . . . ; c2k1 ;M T ;
3. Chebyshev wavelet finite difference method In this Section, we present the Chebyshev wavelet finite difference method (CWFD). It makes sense to give first a brief introduction to Chebyshev finite difference method. Clenshaw (1960) introduced the following approximation of the function f ðtÞ denoted by ðPM Þf ðtÞ as
m¼0
M
2 X00 f ðt k ÞT m ðt k Þ; M k¼0
1 2k
ðtm þ 2n 1Þ;
ðPM Þf ðtÞ ¼
WðtÞ ¼ ½w1;0 ; . . . ; w1;M ; w2;0 ; . . . ; w2;M ; . . . ; w2k1 ;0 ; . . . ; w2k1 ;M T : ð2:12Þ
fm ¼
ð3:8Þ
ð3:9Þ
then a function f ðtÞ can be written in terms of Chebyshev wavelet basis functions as follows
where C and WðtÞ are 2 ðM þ 1Þ 1 matrices given by
fm T m ðtÞ;
ð3:7Þ
with h0 ¼ hM ¼ 12 ; hj ¼ 1 for j ¼ 1; 2; . . . M 1. As can be seen from Eq. (3.6), the first two derivatives of the function f ðtÞ at any point of the Chebyshev–Gauss–Lobatto points is expanded as a linear combination of the values of the function as these points. Now consider t nm ; n ¼ 1; 2; . . . ; 2k1 T, k is a finite integer number, m ¼ 0; 1; . . . ; M, as the corresponding Chebyshev–Gauss–Lobh n such that atto collocation points at the nthsubinterval 2n1 k1 ; k1 2
ð2:11Þ
k
M X 00
ð3:6Þ
where
tnm ¼
n¼1 m¼0
ðPM Þf ðtÞ ¼
n ¼ 1; 2;
j¼0
ð2Þ
n¼1 m¼0
f ðtÞ ¼
M X ðnÞ dm;j f ðtj Þ;
dm;j ¼
A function f ðtÞ 2 L2 ½0; 1 may be expanded as
2k1 X M X
ð3:5Þ
and
2.3. Function approximation
f ðtÞ ¼
ð3:4Þ
M 2 X00 mkp : f ðtk Þ cos M k¼0 M
fwn;m ðtÞj
n ¼ 1; 2; . . . ; 2k1 ; m 2 N [ f0gg forms an orthonormal basis for L2 ½0; 1 with respect to the weight function wn ðtÞ, see Derili and Sohrabi (2008). For k ¼ 2 and M ¼ 2, Chebyshev wavelets are as following:
w1;1
ð3:3Þ
Using Eq. (3.2) we have,
mkp ; T m ðtk Þ ¼ cos M
f ðnÞ ðt m Þ ¼ 2.1. The
M ðtÞ t 2 Þ dTdt .
so fm can be rewritten as
and m ¼ 0; 1; . . . ; M and n ¼ 1; . . . ; 2k1 . In Eq. (2.8) the coefficients are used for orthonormality. Note that by dealing with the Chebyshev wavelets, to get orthogonal wavelets, the weight func~ tion wðtÞ ¼ wð2t 1Þ need to be dilated and translated as follows
Lemma
ð3:2Þ
The well known Chebyshev–Gauss–Lobatto interpolated points,
ð2:7Þ
where
w1;0
m ¼ 0; 1; 2; . . . ; M:
ð3:1Þ
M 2k1 X X 00
cnm wnm ðtÞ;
ð3:10Þ
n¼1 m¼0
where cnm ; n ¼ 1; 2; . . . ; 2k1 T; m ¼ 0; 1; . . . ; M, are the expansion h n coefficients of the function f ðtÞ at the subinterval 2n1 and k1 ; k1 2 wn;m ðtÞ; n ¼ 1; 2; . . . ; 2k1 T; m ¼ 0; 1; . . . ; M, are defined on ½0; TÞ as in Eq. (2.7) and (2.8). In view of Eqs. (2.5) and (3.1), we can obtain the coefficients cnm as M
2 X00 f ðt np Þwnm ðt np Þ ð2k=2 pm Þ M p¼0 M mpp 1 2 X00 : f ðt np Þ cos ¼ k=2 M M 2 pm p¼0
cnm ¼
1
2
ð3:11Þ
181
A. Kazemi Nasab et al. / New Astronomy 34 (2015) 178–186
Using Eq. (3.7) and (3.8), the first two derivatives of the function f (t) at the points t nm ; n ¼ 1; 2; . . . ; 2k1 T; m ¼ 0; 1; . . . ; M, can be obtained as
f ðrÞ ðtnm Þ ¼
M X ðrÞ dn;m;j f ðt nj Þ;
r ¼ 1; 2;
ð3:12Þ
j¼0
where ð1Þ
M k1 X 4hj X khk w ðt nj Þwnl ðtnm Þ; M k¼1 l¼0;ðkþlÞodd cl pk pl nk M k1 X 4hj X 2k khk kjp lmp cos ; ¼ cos M k¼1 l¼0;ðkþlÞodd cl M M
dn;m;j ¼
Table 2 Comparison of yðtÞ values of standard Lane–Emden equation, for the present method and exact numerical values (Horedt, 1986), for m ¼ 2 with M ¼ 8; k ¼ 4. t
Our method
Exact values
0.5 1 1.5 2 2.5 3 3.5 4
0.9593527 0.8486541 0.6953671 0.5298364 0.3747393 0.2418241 0.133969 0.0488401
0.9593527 0.8486541 0.6953671 0.5298364 0.3747393 0.2418241 0.133969 0.0488402
11=2 #2 M X cnm wn;m ðtÞ wn ðtÞdtA ;
and
ð3:13Þ where
0
and ð2Þ
2 2 M k2 X 2hj X 2k kðk l Þhk wnk ðtnj Þwnl ðt nm Þ; M k¼2 l¼0;ðkþlÞev en cl pk pl 2 2 M k2 X 2hj X 4k kðk l Þhk kjp lmp cos : ð3:14Þ ¼ cos M k¼2 l¼0;ðkþlÞev en cl M M
¼@
rk;M
dn;m;j ¼
Z
1
0
" f ðtÞ
k1 2X T
n¼1 m¼0
pffiffiffiffiffiffiffi 2pB Rp3 þ : 5 12 ð2nÞ2
C¼
5. Discretization of problem 4. Convergence analysis In this section, we investigate the convergence analysis of the proposed method. Interested reader can find a proof of the theorems in Appendix A. Lemma 4.1. If the Chebyshev wavelet expansion of a continuous function f ðtÞ converges uniformly, then the Chebyshev wavelet expansion converges to the function f ðtÞ, see Adibi and Assari (2010) and Folland (1999). Theorem 4.2. A function f ðtÞ 2 L2wn ½0; 1Þ, with bounded second derivative, say jf 00 ðtÞj 6 B, can be expanded as an infinite sum of Chebyshev wavelets, and the series converges uniformly to f ðtÞ, that is, 2k1 X 1 X ^cnm wnm ðtÞ: f ðtÞ ¼
ð4:1Þ
n¼1 m¼0
L2wn ½0; 1Þ
Theorem 4.3. Suppose f ðtÞ 2 with bounded second derivative, say jf 00 ðtÞj 6 B, then its Chebyshev wavelet finite difference expansion converges uniformly to f ðtÞ, that is, 1 2k1 X X 0
cnm wnm ðtÞ ¼ f ðtÞ;
ð4:2Þ
n¼1 m¼0
where the summation symbol with prime denotes a sum with the first term halved.
In this Section, the Chebyshev wavelet finite difference method (CWFD) is employed for solving the Lane–Emden type equations. Consider the following general form:
a
y00 ðtÞ þ y0 ðtÞ þ pðtÞqðyðtÞÞ ¼ rðtÞ; t
a P 0; t 2 ½0; 1
ð5:1Þ
subject to the initial conditions
y0 ð0Þ ¼ B;
yð0Þ ¼ A;
ð5:2Þ
where a; A and B are real constants and pðtÞ; qðyðtÞÞ and rðtÞ are some given functions. Now suppose the interval ½0; 1 is divided into h k1 n 2k1 subintervals In ¼ 2n1 . We also consider k1 ; k1 ; n ¼ 1; 2; . . . ; 2 2 the Chebyshev–Gauss–Lobatto collocation points
n1 2k1
¼ t n0 < t n1 < < tn;M1 < tn;M ¼
n 2k1
on the nth subinterval In ; n ¼ 1; 2; . . . ; 2k1 , where tns is defined as following:
tns ¼
1 2k
ðt s þ 2n 1Þ;
s ¼ 1; 2; . . . ; M 1:
ð5:3Þ
In order to obtain the solution yðtÞ in Eq. (5.1), we first consider yn ðtÞ as restriction of yðtÞ on In and then approximate yn ðtÞ in terms of CWFD basis functions as follows
yn ðtÞ ¼
M X 00
cnm wnm ðtÞ;
ð5:4Þ
m¼0
Theorem 4.4 (Accuracy estimation). Suppose f ðtÞ 2 L2wn ½0; 1Þ with bounded second derivative, say jf 00 ðtÞj 6 B, then we have the following accuracy estimation:
rk;M 6
2k1 X 1 X
!1=2
C2
2 n¼1 m¼Mþ1 ðm 1Þ
2
ð4:3Þ
;
where cnm ; n ¼ 1; 2; . . . ; 2k1 ; m ¼ 0; 1; . . . ; M, are defined in (3.11). We now collocate Eq. (5.1) at the Chebyshev–Gauss–Lobatto points tns ; s ¼ 1; . . . ; M 1, as
y00n ðt ns Þ þ
a tns
y0n ðtns Þ þ pðt ns Þqðyn ðt ns ÞÞ ¼ rðt ns Þ:
ð5:5Þ
In view of Eq. (3.12), we obtain
Table 1 Comparison of the first zeros of standard Lane–Emden equations, for the present method, other methods and exact numerical values (Horedt, 2004). m
Our method
Ramos (2005)
0 1 2 3 4
2.4494897427832 3.1415926535898 4.3528746 6.89684862 14.9715464
2.44899 3.14048 4.35086 6.89312 14.96518
Parand et al. (2010a)
4.3528746 6.89684862 14.9715463
Parand et al. (2010b)
Exact values
4.352875 6.8968486 14.971546
2.4494897427832 3.1415926535898 4.3528746 6.89684862 14.9715463
182
A. Kazemi Nasab et al. / New Astronomy 34 (2015) 178–186
Table 3 Comparison of yðtÞ values of standard Lane–Emden equation, for the present method, other methods and exact numerical values (Horedt, 2004), for m ¼ 3 with M ¼ 10; k ¼ 5. t
Our method
0 0.1 0.5 1 5 6 6.8 6.896
Error
Parand et al. (2010a) 0
1 0.998335830 0.959839070 0.855057569 0.110819835 0.043737984 0.004167789 3.60111E05
0.00 10 2.96 108 3.01 108 8.60 109 3.51 108 3.89 109 3.65 1010 5.45 1011
0
0.00 10 1.40 106 2.99 106 1.99 106 3.90 107 1.14 106 1.05 105 8.88 108
Table 4 Comparison of yðtÞ values of standard Lane–Emden equation, for the present method, method in Parand et al. (2010a) and exact numerical values (Horedt, 2004), for m ¼ 4 with M ¼ 10; k ¼ 5.
Parand and Pirkhedri (2010) 0
Parand et al. (2010b) 0
0.00 10 1.82 105 1.12 104 1.07 104 2.02 104 5.00 105 2.93 104 2.10 105
0.00 10 1.28 105 1.81 105 6.00 107 2.00 107 2.00 105 2.79 106 1.01 106
Exact values 1 0.9983358 0.9598391 0.8550576 0.1108198 0.04373798 0.004167789 3.60112E05
Table 6 Comparison of yðtÞ values of isothermal gas sphere, for the present method, method in Parand et al. (2010a) and series solution given by Wazwaz (2001), for m ¼ 4 with M ¼ 12; k ¼ 2.
t
Our method
Error
Parand et al. (2010a)
Exact values
t
Our method
Error
Parand et al. (2010a)
Wazwaz
0 0.1 0.2 0.5 1 5 10 14 14.9
1 0.99833666 0.993386214 0.960310902 0.860813812 0.235922731 0.059672742 0.008330527 0.000576419
0.00 100 4.05 108 1.35 108 2.36 109 1.22 108 3.09 108 1.59 109 2.79 1010 3.46 1012
0.00 100 2.51 104 2.48 104 2.05 104 1.93 104 8.59 105 6.22 105 2.47 105 5.00 107
1 0.9983367 0.9933862 0.9603109 0.8608138 0.2359227 0.05967274 0.008330527 0.000576419
0 0.1 0.2 0.5 1 1.5 2 2.5
0 0.001666419 0.006653971 0.041154515 0.158828174 0.338019831 0.559823312 0.806341085
0.00 100 3.81 1011 3.98 1013 4.92 1010 3.24 107 6.31 106 1.40 104 3.67 103
0.00 100 5.85 107 6.04 107 5.58 107 8.20 107 6.72 106 1.39 104 3.68 103
0 0.001665834 0.006653367 0.041153957 0.158827354 0.33801311 0.55996266 0.810019671
Table 5 ~T for m ¼ 5. The maximum absolute errors ET and E L 50 30 10 5
Et
~ Et 7
1.4 10 2.0 1010 2.0 1014 1.8 1017
3.0 107 1.6 109 1.0 1013 1.0 1016
Fig. 2. Graph of solutions to the Lane–Emden equation modelling an isothermal gas sphere in comparing the current method and Wazwaz solution (Wazwaz, 2001). Table 7 Comparison of yðtÞ values in Example 6.3, for the present method, method in Parand et al. (2010a) and exact solution. t
Fig. 1. Graph of solutions to the standard Lane–Emden equation for m ¼ 0; 1; 2; 3 and 4 respectively.
0 0.01 0.1 0.5 1 2 3 4 5 6 7 8 9 10
Our method 1 0.00019999 0.019900662 0.446287103 1.386294361 3.218875825 4.605170186 5.666426688 6.516193076 7.221835825 7.824046011 8.34877454 8.813438495 9.230241034
Error
Parand et al. (2010a) 0
0.00 10 2.80 1015 2.30 1015 2.55 1015 1.01 1013 1.07 1013 5.79 1014 1.96 1014 6.06 1015 2.28 1014 3.36 1014 4.05 1014 4.47 1014 3.96 1011
0
0.00 10 2.93 106 3.93 106 3.02 106 9.31 107 5.00 107 8.10 107 7.69 107 6.64 107 5.48 107 1.70 107 1.09 106 1.21 105 3.83 105
Exact values 1 0.00019999 0.019900662 0.446287103 1.386294361 3.218875825 4.605170186 5.666426688 6.516193076 7.221835825 7.824046011 8.34877454 8.813438495 9.230241034
183
A. Kazemi Nasab et al. / New Astronomy 34 (2015) 178–186
Fig. 3. (a) The absolute error between approximate and exact solution, (b) the absolute error between the first derivative of approximate and exact solution for Example 6.3.
yðrÞ n ðt ns Þ ¼
M X ðrÞ dn;s;j yn ðtnj Þ;
r ¼ 1; 2;
ð5:6Þ
~T ¼ maxfjy0 ðtÞ y0 ðtÞj : 0 6 t 6 Lg; E e a
j¼0 ðrÞ dn;s;j ;
where r ¼ 1; 2 are defined in (3.13) and (3.14). In addition, substituting Eqs. (3.10) and (3.12) into the initial conditions (5.2), we get other two equations. Besides, it is also necessary that the approximate solution and it’s first derivative be continuous at the interface of subintervals, then
yðrÞ n
n
2k1
¼
ðrÞ
lim yn1 ðtÞ:
t!
ð5:7Þ
n 2k1
Therefore, we have a system of ðM þ 1Þ algebraic equations, which can be solved for the unknowns yn ðt ns Þ. Consequently, we obtain the solution yn ðtÞ to the given Eqs. (5.1) and (5.2). Note that for solving the given equation at the interval I1 , we make use of the given initial conditions and for other subintervals, we use Eq. (5.7) to get initial conditions. After that we obtain yn ðtÞ and the summation will be the required solution yt . 6. Illustrative examples The Lane–Emden equation can be analytically solved only for a few special, integer values of certain index. For other values of index, the numerical solutions can be obtained. In this Section, we consider some numerical examples to demonstrate the validity of the method (CWFD) that we propose to solve the Lane–Emden type equations. The reason that we consider these examples there are closed form solutions that available in the literature. This allows one to compare the results which were obtained by using this method and the analytical solution as well as with the solutions obtained by using other methods. Example 6.1. As the first example, we consider the standard Lane– Emden equation that was used to model the thermal behavior of a spherical cloud of gas acting under the mutual attraction of its molecules. If we set pðtÞ ¼ 1; qðyÞ ¼ ym ; rðtÞ ¼ 0; a ¼ 2; A ¼ 1 and B ¼ 0 in Eq. (1.2), we get
2 y00 ðtÞ þ y0 ðtÞ þ yðtÞm ¼ 0; t
tP0
subject to the boundary conditions
yð0Þ ¼ 1;
and
y0 ð0Þ ¼ 0:
Consider
ET ¼ maxfjye ðtÞ ya ðtÞj : 0 6 t 6 Lg
ð6:1Þ
where ya ðtÞ and ye ðtÞ denote the approximate and the exact solutions, respectively. We applied the method introduced in Section 6 and solved this problem for M ¼ 10; k ¼ 5. In Table 1 we compare the first zeros of standard Lane–Emden equations obtained using the present method with M ¼ 10; k ¼ 5, those reported in literature resulted from other well-known methods and exact values given by Horedt (2004) for m ¼ 0; 1; 2; 3 and 4, respectively. As can be seen from Table 1, the present method provides very accurate predictions of the zeros of yðtÞ even in large intervals. Tables 2–4 show the approximation of yðtÞ for the standard Lane–Emden equation for m ¼ 2; 3 and 4 respectively obtained by the proposed method and those reported in Horedt (2004, 1986). In addition, Table 5 shows the maximum absolute errors for m ¼ 5 and different values of L, for which the exact solution exists. Our obtained results confirm the accuracy of the present method even for large values of L. The resulting graph for m ¼ 0; 1; 2; 3 and 4 is shown in Fig. 1. Next is a display of the solutions for several values of m. Example 6.2. As the second example, we consider the isothermal gas spheres equation with the following formula:
2 y00 ðtÞ þ y0 ðtÞ þ eyðtÞ ¼ 0; t
tP0
ð6:2Þ
subject to the boundary conditions
yð0Þ ¼ 0;
y0 ð0Þ ¼ 0:
Wazwaz (2001), Liao (2003), Singh et al. (2009) and Ramos (2008) obtained a series solution for this equation using ADM, MHAM and series expansion respectively:
1 1 4 8 6 122 8 61:67 10 yðtÞ ffi x2 þ x x þ x x : 6 5:4! 21:6! 81:8! 495:10! This equation has been also solved using Hermite functions collocation, HPM, series solutions, ADM and HAM methods respectively in Parand et al. (2010a), Chowdhury and Hashim (2009), Aslanov (2008), Aslanov (2008) and Bataineh et al. (2009). We apply Chebyshev wavelet finite difference method to solve the equation with M ¼ 12; k ¼ 2. In Table 6 we make a comparison among the values of yðtÞ obtained by the present method and those obtained in Wazwaz (2001) and Parand et al. (2010a) which
184
A. Kazemi Nasab et al. / New Astronomy 34 (2015) 178–186
indicate that our results are more accurate than those reported by Parand et al. (2010a). Fig. 2 shows the resulting graph of the isothermal gas spheres equation in comparison to the current method and that obtained by Wazwaz (2001). Example 6.3. If we set pðtÞ ¼ 1; qðyðtÞÞ ¼ 4 2eyðtÞ þ eyðtÞ=2 ; rðtÞ ¼ 0; A ¼ 0 and B ¼ 0 in Eq. (5.1), we will get one of the Lane–Emden type equations with the following formula:
2 y00 ðtÞ þ y0 ðtÞ þ 4 2eyðtÞ þ eyðtÞ=2 ¼ 0; t
tP0
ð6:3Þ
subject to the boundary conditions
yð0Þ ¼ 0;
y0 ð0Þ ¼ 0;
which has the analytical solution yðtÞ ¼ 2 lnð1 þ t 2 Þ. This equation has been solved in Parand et al. (2010a), Chowdhury and Hashim (2009) and Yildirim and Ozis (2009) with Hermite functions collocation, HPM and VIM method. We apply the introduced method with M ¼ 15; k ¼ 5 for solving Eq. (6.3). In Table 7, we compare the values of yðtÞ obtained by the present method, those reported by Parand et al. (2010a) and exact values. The absolute error between approximate and exact solution, also the absolute error between the first derivative of approximate and exact solution is shown in Fig. 3.
^cnm ¼
rffiffiffiffi Z p cos a þ 2n 1 2 cos ma da f k k p 2 2 2 0 1
using integration by part, we get
pffiffiffi pffiffiffi #p 2 2 cos a þ 2n 1 sin ma f þ pffiffiffiffi k pffiffiffiffi 3k m 2m 2k p 22 p 2 0 Z p cos a þ 2n 1 sin ma sin a da; f0 2k 0
^cnm ¼
the first part is zero, therefore,
^cnm ¼
pffiffiffi Z p 2 0 cos a þ 2n 1 sin ma sin a da: f pffiffiffiffi 3k 2k 22 m p 0
Using integration by part again, it yields
#p 1 sinðm 1Þa sinðm þ 1Þa 0 cos a þ 2n 1 f pffiffiffiffiffiffiffi 3k m1 mþ1 2k 2 2 m 2p 0 Z p 1 00 cos a þ 2n 1 r m ðaÞ da; þ 5k pffiffiffiffiffiffiffi f 2k 2 2 m 2p 0
^cnm ¼
where
sinðm 1Þa sinðm þ 1Þa : rm ðaÞ ¼ sin a m1 mþ1 Thus, we get
7. Conclusion An efficient and accurate method based on hybrid of Chebyshev wavelets and finite difference methods was introduced. The useful properties of Chebyshev wavelets and finite difference method make it a computationally efficient method to approximate the solution of nonlinear Lane–Emden type equations in a semi-infinite interval. We converted the given problem to a system of algebraic equations using collocation points. The main advantage of the present method is the ability to represent smooth and especially piecewise smooth functions properly. It was also shown that the accuracy can be improved either by increasing the number of subintervals or by increasing the number of collocation points in subintervals. Moreover, this method is valid for large-domain calculations. Several examples have been provided to demonstrate the powerfulness of the proposed method. A comparison was made among the present method, some other well-known approaches and exact solution which confirms that the introduced method are more accurate and efficient.
Z p
1
00 cos a þ 2n 1 ^ r m ðaÞ da
jcnm j ¼ 5k pffiffiffiffiffiffiffi f k
2 2 m 2p 0 2 !Z
p
1
f 00 cos a þ 2n 1 rm ðaÞ da 6 pffiffiffiffiffiffiffi 5k
k 2 2 2 m 2p 0 Z p B 6 5k pffiffiffiffiffiffiffi jrm ðaÞj da: 2 2 m 2p 0
Z p
sin a sinðm 1Þa sinðm þ 1Þa da
m1 mþ1 0
Z p
sin a sinðm 1Þa sin a sinðm þ 1Þa
þ
da 6
m1 mþ1 0
0
6
2mp ; m2 1
m > 1:
Since n 6 2k1 , then we obtain
pffiffiffiffiffiffiffi 2pB
The authors are very grateful to the referees for their valuable suggestions and comments that improved the paper.
Now, in the case of m ¼ 1, apply Eq. (7.1), and we have
j^cn1 j <
Here, we include proofs of the theorems, see Kazemi Nasab et al. (2013). Proof of Theorem 4.2. We have
¼
Z
n 2k1
n1 2k1
n;k
¼
Z
1
f ðtÞwnm ðtÞwn ðtÞdt 0
k
22 pm f ðtÞT m ð2k t 2n þ 1Þwð2k t 2n þ 1Þdt;
if m > 1, by substituting 2k t 2n þ 1 ¼ cos a, it yields
ð7:3Þ
jrm ðaÞj da ¼
Acknowledgments
Appendix A
ð7:2Þ
However
Z p
j^cnm j 6
^cn;m ¼ ðf ðtÞ; wn;m ðtÞÞ w
ð7:1Þ
5
ð2nÞ2 ðm2 1Þ
:
pffiffiffi pp ffiffi 2
max jf 0 ðtÞj: 3 ð2nÞ2 06t61
ð7:4Þ
ð7:5Þ
It is also worth to mention that the set fwn0 g1 n¼1 form an orthogonal system which was constructed by Haar scaling function with P 2kffiffiffi ^ p respect to the weight function wðtÞ, so 1 n¼1 p c n0 wn0 ðtÞ is convergent, see Daubechies et al. (1992). Hence, we will have
k1
k1
X
X
X 2 1 2 2k1 X 1
^
X^ c c w ðtÞ 6 w ðtÞ þ j^c jjw ðtÞj
nm nm n0 n0
n¼1
n¼1 m¼1 nm nm
n¼1 m¼0
k1
X
X 2 2k1 X 1
j^c j < 1 6 ^cn0 wn0 ðtÞ þ
n¼1 m¼1 nm
n¼1 P P1 ^ Therefore, with the aid of Lemma 4.1, the series 1 n¼1 m¼1 c nm wnm ðtÞ converges to f ðtÞ uniformly, see Adibi and Assari (2010). h
185
A. Kazemi Nasab et al. / New Astronomy 34 (2015) 178–186
Proof of Theorem 4.3. From Theorem 4.2, we have
Proof of Theorem 4.4. Since
k1
2 1 X X ^cnm wnm ðtÞ;
f ðtÞ ¼
2 k;M
r
n¼1 m¼0
¼
¼
n;k
¼
Z
n 2k1
n1 2k1
¼
Z
f ðtÞwnm ðtÞwn;k ðtÞdt;
¼
0
k 2
2 pm f ðtÞT m ð2k t 2n þ 1Þwð2k t 2n þ 1Þdt;
ð7:6Þ
Z p cos a þ 2n 1 cos ma da: f k 2k 22 0
aj ¼
jp : M
f ðt nj Þ cos
j¼0
Rp3 12M 2
mjp ¼ cnm : M
where
(
) 2 d cosðaÞ þ 2n 1 cosðmaÞ ; 0 6 a 6 p : R ¼ max f da2 2k In view of Eq. (7.1) and by triangle inequality, we get
6
Rp3
12M pffiffiffiffiffiffiffi 2pB
þ 2
pffiffiffiffiffiffiffi 2pB 5
ð2nÞ2 ðm2 1Þ
Rp3 C 6 ; þ 12ðm2 1Þ ð2nÞ52 ðm2 1Þ ðm2 1Þ
ð7:7Þ
where
C¼
pffiffiffiffiffiffiffi Rp3 2p B þ : 5 12 ð2nÞ2
Since jwnm ðtÞj 6 1, then it follows that
k1
k1
k1 1
X
X
X 2 2 2 1
X0 c w ðtÞ 6
c w ðtÞ þ c w ðtÞ
n¼1 m¼0 nm nm n¼1 2 n0 n0 n¼1 n1 n1
þ
2k1 X 1 X
k1
jcnm jjwnm ðtÞj 6
n¼1 m¼2
2 1X jcn0 j 2 n¼1
2k1 2k1 X 1 X X þ jcn1 j þ jcnm j < 1: n¼1
n¼1 m¼0
1 X
c2nm ðwn;m ðtÞÞ2 wn ðtÞdt
n¼1 m¼Mþ1
2k1 X 1 X
c2nm
Z
1
ðwn;m ðtÞÞ2 wn ðtÞdt:
0
2k1 X 1 X
c2nm 6
2k1 X 1 X n¼1 m¼Mþ1
C2 ðm2
2
1Þ
;
where
pffiffiffiffiffiffiffi Rp3 2pB þ : 5 12 ð2nÞ2
References
;
jcnm j 6 jcnm ^cnm j þ j^cnm j 6
n¼1 m¼0 2k1 1X
That completes the proof. h
According to approximation error for the trapezoidal rule, we have
j^cnm cnm j 6
Z
n¼1 m¼Mþ1
C¼
From Eqs. (3.9) and (3.11), for 1 < m < M, we have
k 22 p m M
n¼1 m¼0
" k1 #2 2 1 2k1 X M 1 X X X cnm wn;m ðtÞ cnm wn;m ðtÞ wn;k ðtÞdt
We know that the family fwn;m ðtÞj n 2 N; m 2 N [ f0gg forms orthoR1 normal basis for L2wn ½0; 1Þ, so 0 ðwn;m ðtÞÞ2 wn ðtÞdt ¼ 1. Therefore, in view of Eq. (7.7), we will have
r2k;M ¼
M p p X00 cosðaj Þ þ 2n 1 cos maj ; ’ mk f k M 2 2 2 j¼0
M 00 p 2X
#2 2k1 X M X f ðtÞ cnm wn;m ðtÞ wn ðtÞdt
n¼1 m¼Mþ1
pm
Using the trapezoidal rule for integration with M equidistant subintervals gives
^cnm ’
¼
1
by substituting 2 t 2n þ 1 ¼ cos a, it yields
^cnm
Z
0
k
^cnm ¼
"
0
We prove that this series converges to f ðtÞ uniformly. First of all we show that cnm converges to ^cnm . We have
^cn;m ¼ ðf ðtÞ; wn;m ðtÞÞ w
1
0
where
^cn;m ¼ ðf ðtÞ; wn;m ðtÞÞ : wn
Z
ð7:8Þ
n¼1 m¼2
Therefore, in view of Lemma 4.1, series (4.2) is uniformly convergent to f ðtÞ. h
Adibi, H., Assari, P., 2010. Math. Probl. Eng. 2010. Article ID 138408. Adibi, H., Rismani, A.M., 2010. Comput. Math. Appl. 60, 2126. Aslanov, A., 2008. Phys. Lett. A 372, 3555. Aslanov, A., 2008. Int. J. Comput. Math. 85, 1709. Bataineh, A.S., Noorani, M.S.M., Hashim, I., 2009. Commun. Nonlinear Sci. Numer. Simul. 14, 1121. Bender, C.M., Milton, K.A., Pinsky, S.S., Simmons Jr., L.M., 1989. J. Math. Phys. 30, 1447. Bhrawy, A.H., Alofi, A.S., 2012. Commun. Nonlinear Sci. Numer. Simulat. 17, 62. Boubaker, K., Van Gorder, Robert A., 2012. New Astron. 17, 565. Chandrasekhar, S., 1939. An Introduction to The Study of Stellar Structure. Dover, New York. Chowdhury, M.S.H., Hashim, I., 2009. Nonlinear Anal. Real World Appl. 10, 104. Clenshaw, C.W., Curtis, A.R., 1960. Numer. Math. 2, 197. Daubechies, I., 1992. Ten Lectures on Wavelets. CBMS-NSF Regional Conference Series in Applied Mathematics, vol. 61. SIAM, Philadelphia, Pa, USA. Derili, H., Sohrabi, S., 2008. Math. Sci. (QJMS) 2 (3), 261. Dixon, J.M., Tuszynski, J.A., 1990. Phys. Rev. A 41, 4166. Doha, E.H., Abd- Elhameed, W.M., Youssri, Y.H., 2013. New Astron. 23–24, 113. Eddington, A.S., 1926. The Internal Constitution of The Stars. Cambridge University Press, London. Elbarbary, Elsayed M.E., El-Kady, M., 2003. Appl. Math. Comput. 139, 513. Emden, R., 1907. Gaskugeln Anwendungen der Mechan. Warmtheorie, Druck und Verlag Von B.G. Teubner, Leipzig and Berlin. Fermi, E., 1927. Rend. Accad. Naz. Lincei 6, 602. Folland, G.B., 1999. Real Analysis: Modern Techniques and Their Applications, Pure and Applied Mathematics, second ed. John Wiley & Sons, New York, NY, USA. Fowler, R.H., 1930. MNRAS 91, 63. Frank-Kamenetskii, D.A., 1955. Diffusion and Heat Exchange in Chemical Kinetics. Princeton University Press, Princeton. He, J.H., 2003. Appl. Math. Comput. 143, 539. Herbst, R.S., Momoniat, E., 2012. New Astron. 17, 388. Horedt, G.P., 1986. Astrophys. Space Sci. 126 (2), 357. Horedt, G.P., 2004. Polytropes: Applications in Astrophysics and Related Fields. Kluwer Academic Publishers, Dordrecht. Karimi Vanani, S., Aminataei, A., 2010. Comput. Math. Appl. 59, 2815. Kazemi Nasab, A., Kilicman, A., Pashazadeh Atabakan, Z., Abbasbandy, S., 2013. Abstr. Appl. Anal. [Article ID: 916456]. Lane, J.H., 1870. Am. J. Sci. Arts 50, 57. Liao, S., 2003. Appl. Math. Comput. 142, 1. Mandelzweig, V.B., Tabakin, F., 2001. Comput. Phys. Commun. 141, 268. Marzban, H.R., Tabrizidooz, H.R., Razzaghi, M., 2008. Phys. Lett. A 372 (37), 5883. Misiti, M., Misiti, Y., Oppenheim, G., Poggi, J.-M., 2000. Wavelets Toolbox Users Guide. The MathWorks, Wavelet Toolbox, for use with Matlab. Pandey, Rajesh K., Kumar, Narayan, 2012. New Astron. 17, 303. Pandey, Rajesh K., Kumar, Narayan, 2012. New Astron. 17, 303.
186
A. Kazemi Nasab et al. / New Astronomy 34 (2015) 178–186
Pandey, Rajesh K., Kumar, Narayan, Bhardwaj, Abhinav, Dutta, Goutam, 2012. Appl. Math. Comput. 218, 7629. Parand, K., Pirkhedri, A., 2010. New Astron. 15, 533. Parand, K., Shahini, M., Dehghan, M., 2009. J. Comput. Phys. 228, 8830. Parand, K., Dehghan, M., Rezaeia, A.R., Ghaderi, S.M., 2010a. Comput. Phys. Commun. 181, 1096. Parand, K., Rezaei, A.R., Taghavi, A., 2010b. Acta Astronaut. 67, 673. Ramos, J.I., 2003. Comput. Phys. Commun. 153, 199. Ramos, J.I., 2005. Appl. Math. Comput. 161, 525. Ramos, J.I., 2008. Chaos Solitons Fract. 38, 400. Ramos, J.I., 2009. Chaos Solitons Fract. 40, 1623. Richardson, O.W., 1921. The Emission of Electricity from Hot Bodies. Longmans, Green and co., London, New York.
Shawagfeh, N.T., 1993. J. Math. Phys. 34, 4364. Singh, O.P., Pandey, R.K., Singh, V.K., 2009. Comput. Phys. Commun. 180, 1116. Spitzer, L., 1942. Ap. J. 95, 329. Van Gorder, Robert A., 2011a. New Astron. 16, 65. Van Gorder, Robert A., 2011b. New Astron. 16, 492. Van Gorder, R.A., Vajravelu, K., 2008. Phys. Lett. A 372, 6060. Wazwaz, A., 2001. Appl. Math. Comput. 118, 287. Wazwaz, A., 2006. Appl. Math. Comput. 173, 165. Yildirim, A., Ozis, T., 2009. Ser. A Theor. Methods Appl. 70, 2480. Yousefi, S.A., 2006. Appl. Math. Comput. 181, 1417.