- Email: [email protected]
On the rational second kind Chebyshev pseudospectral method for the solution of the Thomas–Fermi equation over an infinite interval
Journal of Computational and Applied Mathematics 257 (2014) 79–85
Contents lists available at ScienceDirect
Journal of Computational and Applied Mathematics journal homepage: www.elsevier.com/locate/cam
On the rational second kind Chebyshev pseudospectral method for the solution of the Thomas–Fermi equation over an infinite interval A. Kılıçman a,∗ , Ishak Hashim b , M. Tavassoli Kajani c , Mohammad Maleki b a
Department of Mathematics and Institute of Mathematical Research, Universiti Putra Malaysia (UPM), 43400 UPM, Serdang, Selangor, Malaysia b
School of Mathematical Sciences, Universiti Kebangsaan Malaysia, 43600 Bangi Selangor, Malaysia
c
Department of Mathematics, Khorasgan Branch, Islamic Azad University, Isfahan, Iran
article
info
Article history: Received 17 February 2012 Received in revised form 24 June 2013 Keywords: Thomas–Fermi equation Rational second kind Chebyshev functions Pseudospectral method Semi-infinite interval Singular ODE
abstract In this paper, we propose a pseudospectral method for solving the Thomas–Fermi equation which is a nonlinear singular ordinary differential equation on a semi-infinite interval. This approach is based on the rational second kind Chebyshev pseudospectral method that is indeed a combination of tau and collocation methods. This method reduces the solution of this problem to the solution of a system of algebraic equations. The slope at origin is provided with high accuracy. Comparison with some numerical solutions shows that the present solution is effective and highly accurate. © 2013 Elsevier B.V. All rights reserved.
1. Introduction A number of problems arising in science and engineering are set in semi-infinite domains. We can apply different spectral methods that are used to solve problems in semi-infinite domains. The first approach is using Laguerre or Hermite polynomials [1–5]. The second approach is replacing the semi-infinite domain with [0, L] interval by choosing L, sufficiently large. This method is named domain truncation [6,7]. The third approach is reformulating original problem in the semi-infinite domain to singular problem in the bounded domain by variable transformation and then using the Jacobi polynomials to approximate the resulting singular problem [8–12]. The fourth approach of the spectral method is based on rational orthogonal functions. Boyd [13,14] defined a new spectral basis, named rational Chebyshev functions on the semi-infinite interval, by mapping to the Chebyshev polynomials. Guo et al. [15] introduced a new set of rational Legendre functions which are mutually orthogonal in L2 (0, +∞). They applied a spectral scheme using the rational Legendre functions for solving the Korteweg–de Vries equation on the half line. Boyd et al. [16] applied pseudospectral methods on a semi-infinite interval and compared rational Chebyshev, Laguerre and mapped Fourier sine. The authors of [17–19] applied the spectral method to solve nonlinear ordinary differential equations on semi-infinite intervals. Their approach was based on a rational Tau method. They obtained the operational matrices of derivative and product of rational Chebyshev and Legendre functions and then they applied these matrices together with the Tau method to reduce the solution of these problems to the solution of the system of algebraic equations. Furthermore, the authors of [20,21] introduced the rational second and third kind Chebyshev tau method for solving the Lane–Emden equation and Volterra’s population model as nonlinear differential equations over the infinite interval. In [22], rational second kind Chebyshev bases
∗
Corresponding author. E-mail addresses: [email protected], [email protected] (A. Kılıçman).
0377-0427/$ – see front matter © 2013 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.cam.2013.07.050
80
A. Kılıçman et al. / Journal of Computational and Applied Mathematics 257 (2014) 79–85
and Galerkin method are used to obtain the approximate solution of a system of high-order integro-differential equations on the semi-infinite interval. One of the most important nonlinear singular ordinary differential equations that occurs in a semi-infinite interval is the Thomas–Fermi equation, as follows [23,24]: d2 y dx2
1
3
= √ y 2 (x), x
x > 0,
(1)
which appears in the problem of determining the effective nuclear charge in heavy atoms. Boundary conditions for this equation are as follows: y(0) = 1,
lim y(x) = 0.
(2)
x→∞
The Thomas–Fermi equation is useful for calculating form factors and for obtaining effective potentials which can be used as initial trial potentials in self-consistent field calculations. For an understanding of the early developments of the Thomas–Fermi equation, we refer the reader to see the review [25], which especially pointed out the considerable contributions to this field made by Fermi and Majorana. The Thomas–Fermi problem has been solved by different techniques. The authors of [26–28] used a perturbative approach to determine analytic solutions for the Thomas–Fermi equation. Bender et al. [26] replaced the right-hand side of this equation by one which contains the parameter δ , i.e., y′′ (x) = y(x) series in δ y = y0 + δ y1 + δ 2 y2 + δ 3 y3 + · · · .
y(x) x
δ
; the potential is then expanded in a power
(3)
This procedure reduced Eq. (1) into a set of linear equations with associated boundary conditions. Laurenzi [27] applied a perturbative method by combining it with an alternate choice of the nonlinear term of Eq. (1) to produce a rapidly converging analytic solution. Cedillo [28] wrote Eq. (1) in terms of density and then the δ -expansion was employed to obtain an absolute converging series of equations. Adomian [29] applied the decomposition method for solving the Thomas–Fermi equation and then Wazwaz [30] proposed a non-perturbative approximate solution to this equation by using the modified decomposition method in a direct manner without any need to a perturbative expansion or restrictive assumptions. He combined the series obtained with the Padé approximation which provided a promising tool to handle problems on an unbounded domain. Liao [31] solved the Thomas–Fermi equation by the homotopy analysis method. This method provided a convenient way to control the convergence of approximation series and adjusted convergence regions when necessary, which was a fundamental qualitative difference in analysis between the homotopy analysis method and all other reported analytic techniques. Khan [32] used the homotopy analysis method with a new and better transformation which improved the results in comparison with Liao’s work. In [33], the quasilinearization approach was applied for solving Eq. (1). This method approximated the solution of a nonlinear differential equation by treating the nonlinear terms as a perturbation about the linear ones, and unlike perturbation theories it is not based on the existence of some kind of a small parameter. Ramos [34] presented two piecewise quasilinearization methods for the numerical solution of Eq. (1). Both methods were based on the piecewise linearization of ordinary differential equations. The first method (C1-linearization) provided global smooth solutions, whereas the second one (C0-linearization) provided continuous solutions. Yao [34] solved the Thomas–Fermi equation by an analytic technique named the homotopy analysis method with a more generalized set of basis function, and consequential auxiliary linear operator was introduced to provide a series solution. Zhu et al. [35] approximated the original Thomas–Fermi equation by a nonlinear free boundary value problem. They then used an iterative method to solve the free boundary value problem. Parand et al. [36] investigated the Sinc-collocation method on the half-line for solving Eq. (1) by using Sinc functions. Similar, for a singular boundary-value problem we refer to. In this paper, we introduce a combination of tau and pseudospectral methods based on rational second kind Chebyshev (RSC) functions. The proposed method requires the definition of RSC functions, operational matrix of derivative and rational second kind Chebyshev–Gauss collocation points and weights. The application of the method to the Thomas–Fermi equation leads to a nonlinear algebraic system. We employed this method to the Thomas–Fermi equation because first, it is easy to apply and numerically achieve spectral convergence; second, because of singularity in this equation this method can handle this problem; and third, the limit of the RSC functions at infinity is computable and thus the boundary conditions at infinity can easily be handled. This paper is arranged as follows: in Section 2, we describe the formulation and some properties of rational second kind Chebyshev functions required for our subsequent development. Section 3 summarizes the application of this method for solving the Thomas–Fermi equation and a comparison is made with existing methods in the literature. The conclusions are described in the final section. 2. Properties of RSC functions In this section, we present some properties of rational second kind Chebyshev functions and introduce the rational second kind Chebyshev pseudospectral approach.
A. Kılıçman et al. / Journal of Computational and Applied Mathematics 257 (2014) 79–85
81
2.1. RSC functions
√ The second kind Chebyshev polynomials are orthogonal in the interval [−1, 1] with respect to the weight function 1 − x2 and we find that Un (x) satisfies the recurrence relation U0 (x) = 1,
U1 (x) = 2x,
Un (x) = 2xUn−1 (x) − Un−2 (x),
(4)
n ≥ 2.
The RSC functions are defined by
R n ( x) = U n
x−L
x+L
;
thus RSC functions satisfy R 0 ( x) = 1 , R n ( x) = 2
R1 (x) = 2
x−L
x+L
x−L
x+L
,
Rn−1 (x) − Rn−2 (x),
(5) n ≥ 2.
2.2. Function approximation √
Let w(x) = (4x+LLx)3 denote a non-negative, integrable, real-valued function over the interval I = [0, +∞). We define L2w (I ) = {y : I → R | y is measurable and ∥y∥w < ∞} ,
(6)
where ∞
|y(x)| w(x)dx 2
∥y∥w =
12
,
(7)
0
is the norm induced by the scalar product
⟨y, z ⟩w =
∞
y(x)z (x)w(x)dx.
(8)
0
Thus {Rn (x)}n≥0 denote a system which is mutually orthogonal under Eq. (8), i.e., ∞
Rn (x)Rm (x)w(x) dx = 0
π 2
δnm ,
(9)
where δnm is the Kronecker delta function. This system is complete in L2w (I ); as a result, any function y ∈ L2w (I ) can be expanded as follows: y(x) =
∞
ak Rk (x),
(10)
k=0
with ak =
2
π
⟨y, Rk ⟩w .
(11)
The ak ’s are the expansion coefficients associated with the family {Rk (x)}. If the infinite series in Eq. (10) is truncated, then it can be written as y(x) ≃
N
ak Rk (x) = AT R(x),
(12)
k=0
where A = [a0 , a1 , . . . , aN ]T , R(x) = [R0 (x), R1 (x), . . . , RN (x)]T . Moreover, from the recurrence relation in Eq. (5) we have R(0) = [1, −2, 3, −4, . . . , (−1)N (N + 1)]T = e1 , lim R(x) = [1, 2, 3, 4, . . . , N + 1]T = e2 .
x→∞
(13) (14)
82
A. Kılıçman et al. / Journal of Computational and Applied Mathematics 257 (2014) 79–85
2.3. Operational matrix of derivative The derivative of the vector R(x) defined in Eq. (12) can be approximated by R′ (x) ≃ DR(x),
(15)
where D is the n × n operational matrix for derivative. Differentiating Eq. (5) we get R′0 (x) = 0,
5
R′1 (x) =
4
R0 (x) − R1 (x) +
R′n (x) = (R1 (x)Rn−1 (x))′ − R′n−2 (x),
1 4
R2 (x)
n ≥ 2.
(16) (17)
By using Eqs. (16) and (17) the matrix D can be calculated. The matrix D is a lower Hessenberg matrix and can be expressed as D = D1 + D2 , where D1 is a tridiagonal matrix which is obtained from
−2 + 7i
D1 = diag
4
, −i ,
i
4
,
i = 0, . . . , n − 1
and the dij elements of matrix D2 are obtained from
dij =
0,
i≤j+1 i > j + 1.
(−1)i+j+1 (2j),
2.4. RSC collocation points and weights Abramowitz and Stegun [37] introduced rational second kind Chebyshev–Gauss points. Let
RN = span{R0 , R1 , . . . , RN },
(18)
and
tj = − cos
(j + 1)π N +2
,
j = 0, 1, . . . , N ,
(19)
are the N + 1 second kind Chebyshev–Gauss points; thus we define xj = L
1 + tj 1 − tj
,
j = 0, 1, . . . , N ,
(20)
which are named rational second kind Chebyshev–Gauss nodes. The constant L is a user-defined mapping parameter. Boyd [38,39] introduced some sophisticated ways to estimate the best choice of L; however, in practice the criterion for ‘‘optimum’’ is the rate of convergence. The relations between rational second kind Chebyshev orthogonal systems and rational second kind Gauss integration are as follows: ∞
y(x)w(x)dx =
1
y L
−1
0
=
N
1+t
1−t
y(xj )wj ,
ρ(t )dt
∀y ∈ R2N +1 ,
(21)
j =0
where ρ(t ) =
√
1 − t 2 and wj = Nπ+2 sin2
(j+1)π N +2
, 0 ≤ j ≤ N.
3. Numerical solution of the Thomas–Fermi equation In this phase at first we rewrite the Thomas–Fermi equation introduced in Eqs. (1) and (2) as
√
3
xy′′ (x) − y 2 (x) = 0,
y(0) = 1,
(22)
lim y(x) = 0.
x→∞
By applying Eqs. (12) and (15) on Eq. (22) we define Res(x) =
√
3
xAT D2 R(x) − (AT R(x)) 2 .
(23)
A. Kılıçman et al. / Journal of Computational and Applied Mathematics 257 (2014) 79–85
83
Table 1 Approximations of y′ (0) for the present method. N
L
y′ (0)
Relative error
5 7
0.1485860 0.0958885
−1.588071834 −1.588071347
8.34 × 10−7 3.47 × 10−7 y′ (0) = −1.588071
Kobayashi result [40]
Table 2 Comparison between methods in [31,32,41,35] and the present method for y′ (0). Padé
Liao [31]
Khan [32]
Yao [41]
N
Zhu [35]
[20, 20] −1.58281 −1.582901807 −1.585148733 1000 −1.58794357 [30, 30] −1.58606 −1.586494973 −1.588004950 2000 −1.58807411
N
Present method
5 7
−1.588071834 −1.588071347
Fig. 1. Thomas–Fermi graph obtained by the present method.
As in a typical tau method and using Eq. (8) we can write
⟨Res(x), Rk (x)⟩w =
∞
Res(x)Rk (x)w(x)dx = 0,
k = 0, 1, . . . , N − 2.
(24)
0
Now, a pseudospectral method is defined by applying Eq. (21) on (24) to generate (N − 1) algebraic equations as follows: N
Res(xj )Rk (xj )wj = 0,
k = 0, 1, . . . , N − 2.
(25)
j=0
In addition, using Eqs. (12)–(14) the boundary conditions in Eq. (22) can be approximated as AT e1 = 1,
AT e2 = 0.
(26)
Solving the system of N + 1 nonlinear equations in Eqs. (25) and (26) using Newton’s method for the unknown coefficients aj and substituting the obtained results in Eq. (12) the values of y(x) can be obtained. The initial slope y′ (0) is difficult to compute by any means and plays an important role in determining many physical properties of the Thomas–Fermi atom. It determines the energy of a neutral atom in the Thomas–Fermi approximation:
E=
6 7
4π 3
32
7
Z 3 y′ (0),
(27)
where Z is the nuclear charge. Kobayashi [40] provided a highly accurate numerical value of the initial slope as y′ (0) = −1.588071. The approximations of y′ (0) computed by the present method and their relative errors are shown in Table 1. Obviously, this method is convergent by increasing the number of collocation points and obtaining suitable L. The comparison of the initial slope y′ (0) calculated by the present work with values obtained by Liao [31], Khan [32], Yao [41] and Zhu et al. [35] is given in Table 2, which shows that the present solution is highly accurate. Table 3 shows the approximations of y(x) obtained by the method proposed in this paper for N = 7 and L = 0.0958885, and those obtained by Khan [32] and Liao [42]. Fig. 1 shows the resulting graph of Thomas–Fermi for N = 7 and L = 0.0958885 which tends to zero as x increases by the boundary condition limx→∞ y(x) = 0.
84
A. Kılıçman et al. / Journal of Computational and Applied Mathematics 257 (2014) 79–85 Table 3 Approximation of y(x) for present method, [32,42]. x
Khan [32]
Liao [42]
Present method
0.25 0.50 0.75 1.00 1.25 1.50 1.75 2.00 2.25 2.50 2.75 3.00 3.25 3.50 3.75 4.00 4.25 4.50 4.75 5.00 6.00 7.00 8.00 9.00 10.0 15.0 20.0 25.0 50.0 75.0 100.0 1000.0 5000.0
0.776191000 0.615917000 0.505380000 0.423772000 0.362935000 0.314490000 0.275154000 0.242718000 0.215630000 0.192795000 0.173364000 0.156719000 0.142371000 0.129937000 0.119108000 0.109632000 0.101303000 0.093950400 0.087432000 0.081629600 0.063816200 0.051800500 0.043285900 0.037002300 0.032208100 0.019184300 0.013493700 0.010357000 0.004730890 0.003052460 0.002251000 0.000214641 –
0.755202000 0.606987000 0.502347000 0.424008000 0.363202000 0.314778000 0.275451000 0.243009000 0.215895000 0.192984000 0.173441000 0.156633000 0.142070000 0.129370000 0.118229000 0.108404000 0.099697900 0.091948200 0.085021800 0.078807800 0.059423000 0.046097800 0.036587300 0.029590900 0.024314300 0.010805400 0.005784940 0.003473750 0.000632255 0.000218210 0.000100243 0.000000135 –
0.755903399 0.605270502 0.497823942 0.420343948 0.362814756 0.318737461 0.284014958 0.256010764 0.232974191 0.213705386 0.197357837 0.183318729 0.171134229 0.160461449 0.151036695 0.142653971 0.135150082 0.128394110 0.122279834 0.116720187 0.098752518 0.085573204 0.075494762 0.067538722 0.061098873 0.041370727 0.031271686 0.025135426 0.012687078 0.008484835 0.006373709 0.000640105 0.000128070
4. Conclusion The fundamental goal of this paper has been to construct an approximation to the solution of the nonlinear Thomas–Fermi equation in a semi-infinite interval which has singularity at x = 0 and its boundary condition occurs in infinity. In the above discussion, the combination of tau and pseudospectral methods with RSC functions, which have the property of orthogonality, is employed to achieve this goal. Advantages of this method are that we do not reform the problem to a finite domain and with a small N very accurate results are obtained. There is a good agreement between obtained results for y′ (0) and exact value that demonstrates the validity of the present method for this type of problems and gives the method a wider applicability. Comparing the computed results by this method with the others shows that this method provides more accurate and numerically stable solutions than those obtained by other methods. Several possible extensions of this work can be done for further research. Some of these areas include numerical solutions of integral equations defined on infinite domains, infinite horizon optimal control problems and partial boundary value problems over infinite intervals. References [1] B.Y. Guo, J. Shen, Laguerre–Galerkin method for nonlinear partial differential equations on a semi-infinite interval, Numer. Math. 86 (4) (2000) 635–654. [2] J. Shen, Stable and efficient spectral methods in unbounded domains using Laguerre functions, SIAM J. Numer. Anal. 38 (4) (2000) 1113–1133. [3] Y. Maday, B. Pernaud-Thomas, H. Vandeven, Shock-fitting techniques for solving hyperbolic problems with spectral methods, Rech. Aerosp. 6 (1985) 1–9. [4] H.I. Siyyam, Laguerre Tau methods for solving higher-order ordinary differential equations, J. Comput. Anal. Appl. 3 (2) (2001) 173–182. [5] K. Parand, M. Dehghan, A.R. Rezaei, S.M. Ghaderi, An approximation algorithm for the solution of the nonlinear Lane–Emden type equations arising in astrophysics using Hermite functions collocation method, Comput. Phys. Commun. 181 (6) (2010) 1096–1108. [6] J.P. Boyd, Chebyshev and Fourier Spectral Methods, second ed., Springer, Berlin, 2000. [7] M. Maleki, M. Tavassoli Kajani, I. Hashim, A. Kılıçman, K.A.M. Atan, A nonclassical radau collocation method for nonlinear initial-value problems with applications to Lane–Emden type equations, J. Appl. Math. 2012 (2012), Article ID 103205, 13 pages. http://dx.doi.org/10.1155/2012/103205. [8] B.Y. Guo, Jacobi spectral approximations to differential equations on the half line, J. Comput. Math. 18 (1) (2000) 95–112. [9] M. Maleki, I. Hashim, S. Abbasbandy, Analysis of IVPs and BVPs on semi-infinite domains via collocation methods, J. Appl. Math. 2012 (2012), Article ID 696574, 21 pages. http://dx.doi.org/10.1155/2012/696574. [10] P.M. Lima, M.P. Carpentier, Iterative methods for a singular boundary-value problem, J. Comput. Appl. Math. 111 (12) (1999) 173–186. [11] K.Q. Lan, Multiple eigenvalues for singular Hammerstein integral equations with applications to boundary value problems, J. Comput. Appl. Math. 189 (12) (2006) 109–119. [12] H. Allouche, N. Marhraoui, Numerical solution of singular regular boundary value problems by pole detection with qd-algorithm, J. Comput. Appl. Math. 233 (2) (2009) 420–436.
A. Kılıçman et al. / Journal of Computational and Applied Mathematics 257 (2014) 79–85
85
[13] J.P. Boyd, Orthogonal rational functions on a semi-infinite interval, J. Comput. Phys. 70 (1) (1987) 63–88. [14] John P. Boyd, Rational Chebyshev series for the ThomasFermi function: endpoint singularities and spectral methods, J. Comput. Appl. Math. 244 (2013) 90–101. [15] B.Y. Guo, J. Shen, Z.Q. Wang, A rational approximation and its applications to differential equations on the half line, J. Sci. Comput. 15 (2) (2000) 117–147. [16] J.P. Boyd, C. Rangan, P.H. Bucksbaum, Pseudospectral methods on a semi-infinite interval with application to the hydrogen atom: a comparison of the mapped Fourier-sine method with Laguerre series and rational Chebyshev expansions, J. Comput. Phys. 188 (1) (2003) 56–74. [17] K. Parand, M. Razzaghi, Rational Legendre approximation for solving some physical problems on semi-infinite intervals, Phys. Scr. 69 (5) (2004) 353–357. [18] K. Parand, M. Razzaghi, Rational Chebyshev tau method for solving Volterra’s population model, Appl. Math. Comput. 149 (3) (2004) 893–900. [19] K. Parand, M. Razzaghi, Rational Chebyshev tau method for solving higher-order ordinary differential equations, Int. J. Comput. Math. 81 (1) (2004) 73. [20] M. Tavassoli Kajani, F. Ghasemi Tabatabaei, Rational Chebyshev approximations for solving Lane–Emde equation of index m, in: Proceeding of the International Conference on Computational and Applied Mathematics, Bangkok, Thailand, 29–31 March, 2011, pp. 840–844. [21] M. Dadkhah Tirani, F. Ghasemi Tabatabaei, M. Tavassoli Kajani, Rational second (third) kind Chebyshev approximations for solving Volterras population model, in: Proceeding of the International Conference on Computational and Applied Mathematics, Bangkok, Thailand, 29–31 March, 2011, pp. 835–839. [22] M. Tavassoli Kajani, S. Vahdati, Zulkifly Abbas, M. Maleki, Application of rational second kind Chebyshev functions for system of integrodifferential equations on semi-infinite intervals, J. Appl. Math. 2012 (2012), Article ID 803503, 11 pages. http://dx.doi.org/10.1155/2012/803503. [23] H.T. Davis, Introduction to Nonlinear Differential and Integral Equations, Dover, New York, 1962. [24] S. Chandrasekhar, Introduction to the Study of Stellar Structure, Dover, New York, 1967. [25] E. Di Grezia, S. Esposito, Fermi, Majorana and the statistical model of atoms, Found. Phys. 34 (2004) 1431–1450. [26] C.M. Bender, K.A. Milton, S.S. Pinsky, L.M. Simmons Jr., A new perturbative approach to nonlinear problems, J. Math. Phys. 30 (7) (1989) 1447–1455. [27] B.J. Laurenzi, An analytic solution to the Thomas–Fermi equation, J. Math. Phys. 31 (10) (1990) 2535–2537. [28] A. Cedillo, A perturbative approach to the Thomas–Fermi equation in terms of the density, J. Math. Phys. 34 (7) (1993) 2713–2717. [29] G. Adomian, Solution of the Thomas–Fermi equation, Appl. Math. Lett. 11 (3) (1998) 131–133. [30] A. Wazwaz, The modified decomposition method and Padé approximants for solving the Thomas–Fermi equation, Appl. Math. Comput. 105 (1) (1999) 11–19. [31] S. Liao, An explicit analytic solution to the Thomas–Fermi equation, Appl. Math. Comput. 144 (23) (2003) 495–506. [32] H. Khan, H. Xu, Series solution to the Thomas–Fermi equation, Phys. Lett. A 365 (12) (2007) 111–115. [33] V.B. Mandelzweig, F. Tabakin, Quasilinearization approach to nonlinear problems in physics with application to nonlinear ODEs, Comput. Phys. Commun. 141 (2) (2001) 268–281. [34] J.I. Ramos, Piecewise quasilinearization techniques for singular boundary-value problems, Comput. Phys. Commun. 158 (1) (2004) 12–25. [35] S. Zhu, H. Zhu, Q. Wu, Y. Khan, An adaptive algorithm for the Thomas–Fermi equation, Numer. Algorithms 59 (2012) 359–372. [36] K. Parand, M. Dehghan, A. Pirkhedri, The Sinc-collocation method for solving the Thomas–Fermi equation, J. Comput. Appl. Math. 237 (2013) 244–252. [37] M. Abramowitz, I.A. Stegun, Handbook of Mathematical Functions, Dover, New York, 1972. 10th printing with corrections. [38] J.P. Boyd, The optimization of convergence for Chebyshev polynomial methods in an unbounded domain, J. Comput. Phys. 45 (1) (1982) 43–79. [39] J.P. Boyd, Orthogonal rational functions on a semi-infinite interval, J. Comput. Phys. 70 (1987) 63–88. [40] S. Kobayashi, T. Matsukuma, S. Nagi, K. Umeda, Accurate value of the initial slope of the ordinary TF function, J. Phys. Soc. Japan 10 (1955) 759–762. [41] B. Yao, A series solution to the Thomas–Fermi equation, Appl. Math. Comput. 203 (2008) 396–401. [42] S. Liao, Beyond Perturbation-Introduction to the Homotopy Analysis Method, Chapman & Hall/CRC, Boca Raton, 2003.
Contents lists available at ScienceDirect
Journal of Computational and Applied Mathematics journal homepage: www.elsevier.com/locate/cam
On the rational second kind Chebyshev pseudospectral method for the solution of the Thomas–Fermi equation over an infinite interval A. Kılıçman a,∗ , Ishak Hashim b , M. Tavassoli Kajani c , Mohammad Maleki b a
Department of Mathematics and Institute of Mathematical Research, Universiti Putra Malaysia (UPM), 43400 UPM, Serdang, Selangor, Malaysia b
School of Mathematical Sciences, Universiti Kebangsaan Malaysia, 43600 Bangi Selangor, Malaysia
c
Department of Mathematics, Khorasgan Branch, Islamic Azad University, Isfahan, Iran
article
info
Article history: Received 17 February 2012 Received in revised form 24 June 2013 Keywords: Thomas–Fermi equation Rational second kind Chebyshev functions Pseudospectral method Semi-infinite interval Singular ODE
abstract In this paper, we propose a pseudospectral method for solving the Thomas–Fermi equation which is a nonlinear singular ordinary differential equation on a semi-infinite interval. This approach is based on the rational second kind Chebyshev pseudospectral method that is indeed a combination of tau and collocation methods. This method reduces the solution of this problem to the solution of a system of algebraic equations. The slope at origin is provided with high accuracy. Comparison with some numerical solutions shows that the present solution is effective and highly accurate. © 2013 Elsevier B.V. All rights reserved.
1. Introduction A number of problems arising in science and engineering are set in semi-infinite domains. We can apply different spectral methods that are used to solve problems in semi-infinite domains. The first approach is using Laguerre or Hermite polynomials [1–5]. The second approach is replacing the semi-infinite domain with [0, L] interval by choosing L, sufficiently large. This method is named domain truncation [6,7]. The third approach is reformulating original problem in the semi-infinite domain to singular problem in the bounded domain by variable transformation and then using the Jacobi polynomials to approximate the resulting singular problem [8–12]. The fourth approach of the spectral method is based on rational orthogonal functions. Boyd [13,14] defined a new spectral basis, named rational Chebyshev functions on the semi-infinite interval, by mapping to the Chebyshev polynomials. Guo et al. [15] introduced a new set of rational Legendre functions which are mutually orthogonal in L2 (0, +∞). They applied a spectral scheme using the rational Legendre functions for solving the Korteweg–de Vries equation on the half line. Boyd et al. [16] applied pseudospectral methods on a semi-infinite interval and compared rational Chebyshev, Laguerre and mapped Fourier sine. The authors of [17–19] applied the spectral method to solve nonlinear ordinary differential equations on semi-infinite intervals. Their approach was based on a rational Tau method. They obtained the operational matrices of derivative and product of rational Chebyshev and Legendre functions and then they applied these matrices together with the Tau method to reduce the solution of these problems to the solution of the system of algebraic equations. Furthermore, the authors of [20,21] introduced the rational second and third kind Chebyshev tau method for solving the Lane–Emden equation and Volterra’s population model as nonlinear differential equations over the infinite interval. In [22], rational second kind Chebyshev bases
∗
Corresponding author. E-mail addresses: [email protected], [email protected] (A. Kılıçman).
0377-0427/$ – see front matter © 2013 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.cam.2013.07.050
80
A. Kılıçman et al. / Journal of Computational and Applied Mathematics 257 (2014) 79–85
and Galerkin method are used to obtain the approximate solution of a system of high-order integro-differential equations on the semi-infinite interval. One of the most important nonlinear singular ordinary differential equations that occurs in a semi-infinite interval is the Thomas–Fermi equation, as follows [23,24]: d2 y dx2
1
3
= √ y 2 (x), x
x > 0,
(1)
which appears in the problem of determining the effective nuclear charge in heavy atoms. Boundary conditions for this equation are as follows: y(0) = 1,
lim y(x) = 0.
(2)
x→∞
The Thomas–Fermi equation is useful for calculating form factors and for obtaining effective potentials which can be used as initial trial potentials in self-consistent field calculations. For an understanding of the early developments of the Thomas–Fermi equation, we refer the reader to see the review [25], which especially pointed out the considerable contributions to this field made by Fermi and Majorana. The Thomas–Fermi problem has been solved by different techniques. The authors of [26–28] used a perturbative approach to determine analytic solutions for the Thomas–Fermi equation. Bender et al. [26] replaced the right-hand side of this equation by one which contains the parameter δ , i.e., y′′ (x) = y(x) series in δ y = y0 + δ y1 + δ 2 y2 + δ 3 y3 + · · · .
y(x) x
δ
; the potential is then expanded in a power
(3)
This procedure reduced Eq. (1) into a set of linear equations with associated boundary conditions. Laurenzi [27] applied a perturbative method by combining it with an alternate choice of the nonlinear term of Eq. (1) to produce a rapidly converging analytic solution. Cedillo [28] wrote Eq. (1) in terms of density and then the δ -expansion was employed to obtain an absolute converging series of equations. Adomian [29] applied the decomposition method for solving the Thomas–Fermi equation and then Wazwaz [30] proposed a non-perturbative approximate solution to this equation by using the modified decomposition method in a direct manner without any need to a perturbative expansion or restrictive assumptions. He combined the series obtained with the Padé approximation which provided a promising tool to handle problems on an unbounded domain. Liao [31] solved the Thomas–Fermi equation by the homotopy analysis method. This method provided a convenient way to control the convergence of approximation series and adjusted convergence regions when necessary, which was a fundamental qualitative difference in analysis between the homotopy analysis method and all other reported analytic techniques. Khan [32] used the homotopy analysis method with a new and better transformation which improved the results in comparison with Liao’s work. In [33], the quasilinearization approach was applied for solving Eq. (1). This method approximated the solution of a nonlinear differential equation by treating the nonlinear terms as a perturbation about the linear ones, and unlike perturbation theories it is not based on the existence of some kind of a small parameter. Ramos [34] presented two piecewise quasilinearization methods for the numerical solution of Eq. (1). Both methods were based on the piecewise linearization of ordinary differential equations. The first method (C1-linearization) provided global smooth solutions, whereas the second one (C0-linearization) provided continuous solutions. Yao [34] solved the Thomas–Fermi equation by an analytic technique named the homotopy analysis method with a more generalized set of basis function, and consequential auxiliary linear operator was introduced to provide a series solution. Zhu et al. [35] approximated the original Thomas–Fermi equation by a nonlinear free boundary value problem. They then used an iterative method to solve the free boundary value problem. Parand et al. [36] investigated the Sinc-collocation method on the half-line for solving Eq. (1) by using Sinc functions. Similar, for a singular boundary-value problem we refer to. In this paper, we introduce a combination of tau and pseudospectral methods based on rational second kind Chebyshev (RSC) functions. The proposed method requires the definition of RSC functions, operational matrix of derivative and rational second kind Chebyshev–Gauss collocation points and weights. The application of the method to the Thomas–Fermi equation leads to a nonlinear algebraic system. We employed this method to the Thomas–Fermi equation because first, it is easy to apply and numerically achieve spectral convergence; second, because of singularity in this equation this method can handle this problem; and third, the limit of the RSC functions at infinity is computable and thus the boundary conditions at infinity can easily be handled. This paper is arranged as follows: in Section 2, we describe the formulation and some properties of rational second kind Chebyshev functions required for our subsequent development. Section 3 summarizes the application of this method for solving the Thomas–Fermi equation and a comparison is made with existing methods in the literature. The conclusions are described in the final section. 2. Properties of RSC functions In this section, we present some properties of rational second kind Chebyshev functions and introduce the rational second kind Chebyshev pseudospectral approach.
A. Kılıçman et al. / Journal of Computational and Applied Mathematics 257 (2014) 79–85
81
2.1. RSC functions
√ The second kind Chebyshev polynomials are orthogonal in the interval [−1, 1] with respect to the weight function 1 − x2 and we find that Un (x) satisfies the recurrence relation U0 (x) = 1,
U1 (x) = 2x,
Un (x) = 2xUn−1 (x) − Un−2 (x),
(4)
n ≥ 2.
The RSC functions are defined by
R n ( x) = U n
x−L
x+L
;
thus RSC functions satisfy R 0 ( x) = 1 , R n ( x) = 2
R1 (x) = 2
x−L
x+L
x−L
x+L
,
Rn−1 (x) − Rn−2 (x),
(5) n ≥ 2.
2.2. Function approximation √
Let w(x) = (4x+LLx)3 denote a non-negative, integrable, real-valued function over the interval I = [0, +∞). We define L2w (I ) = {y : I → R | y is measurable and ∥y∥w < ∞} ,
(6)
where ∞
|y(x)| w(x)dx 2
∥y∥w =
12
,
(7)
0
is the norm induced by the scalar product
⟨y, z ⟩w =
∞
y(x)z (x)w(x)dx.
(8)
0
Thus {Rn (x)}n≥0 denote a system which is mutually orthogonal under Eq. (8), i.e., ∞
Rn (x)Rm (x)w(x) dx = 0
π 2
δnm ,
(9)
where δnm is the Kronecker delta function. This system is complete in L2w (I ); as a result, any function y ∈ L2w (I ) can be expanded as follows: y(x) =
∞
ak Rk (x),
(10)
k=0
with ak =
2
π
⟨y, Rk ⟩w .
(11)
The ak ’s are the expansion coefficients associated with the family {Rk (x)}. If the infinite series in Eq. (10) is truncated, then it can be written as y(x) ≃
N
ak Rk (x) = AT R(x),
(12)
k=0
where A = [a0 , a1 , . . . , aN ]T , R(x) = [R0 (x), R1 (x), . . . , RN (x)]T . Moreover, from the recurrence relation in Eq. (5) we have R(0) = [1, −2, 3, −4, . . . , (−1)N (N + 1)]T = e1 , lim R(x) = [1, 2, 3, 4, . . . , N + 1]T = e2 .
x→∞
(13) (14)
82
A. Kılıçman et al. / Journal of Computational and Applied Mathematics 257 (2014) 79–85
2.3. Operational matrix of derivative The derivative of the vector R(x) defined in Eq. (12) can be approximated by R′ (x) ≃ DR(x),
(15)
where D is the n × n operational matrix for derivative. Differentiating Eq. (5) we get R′0 (x) = 0,
5
R′1 (x) =
4
R0 (x) − R1 (x) +
R′n (x) = (R1 (x)Rn−1 (x))′ − R′n−2 (x),
1 4
R2 (x)
n ≥ 2.
(16) (17)
By using Eqs. (16) and (17) the matrix D can be calculated. The matrix D is a lower Hessenberg matrix and can be expressed as D = D1 + D2 , where D1 is a tridiagonal matrix which is obtained from
−2 + 7i
D1 = diag
4
, −i ,
i
4
,
i = 0, . . . , n − 1
and the dij elements of matrix D2 are obtained from
dij =
0,
i≤j+1 i > j + 1.
(−1)i+j+1 (2j),
2.4. RSC collocation points and weights Abramowitz and Stegun [37] introduced rational second kind Chebyshev–Gauss points. Let
RN = span{R0 , R1 , . . . , RN },
(18)
and
tj = − cos
(j + 1)π N +2
,
j = 0, 1, . . . , N ,
(19)
are the N + 1 second kind Chebyshev–Gauss points; thus we define xj = L
1 + tj 1 − tj
,
j = 0, 1, . . . , N ,
(20)
which are named rational second kind Chebyshev–Gauss nodes. The constant L is a user-defined mapping parameter. Boyd [38,39] introduced some sophisticated ways to estimate the best choice of L; however, in practice the criterion for ‘‘optimum’’ is the rate of convergence. The relations between rational second kind Chebyshev orthogonal systems and rational second kind Gauss integration are as follows: ∞
y(x)w(x)dx =
1
y L
−1
0
=
N
1+t
1−t
y(xj )wj ,
ρ(t )dt
∀y ∈ R2N +1 ,
(21)
j =0
where ρ(t ) =
√
1 − t 2 and wj = Nπ+2 sin2
(j+1)π N +2
, 0 ≤ j ≤ N.
3. Numerical solution of the Thomas–Fermi equation In this phase at first we rewrite the Thomas–Fermi equation introduced in Eqs. (1) and (2) as
√
3
xy′′ (x) − y 2 (x) = 0,
y(0) = 1,
(22)
lim y(x) = 0.
x→∞
By applying Eqs. (12) and (15) on Eq. (22) we define Res(x) =
√
3
xAT D2 R(x) − (AT R(x)) 2 .
(23)
A. Kılıçman et al. / Journal of Computational and Applied Mathematics 257 (2014) 79–85
83
Table 1 Approximations of y′ (0) for the present method. N
L
y′ (0)
Relative error
5 7
0.1485860 0.0958885
−1.588071834 −1.588071347
8.34 × 10−7 3.47 × 10−7 y′ (0) = −1.588071
Kobayashi result [40]
Table 2 Comparison between methods in [31,32,41,35] and the present method for y′ (0). Padé
Liao [31]
Khan [32]
Yao [41]
N
Zhu [35]
[20, 20] −1.58281 −1.582901807 −1.585148733 1000 −1.58794357 [30, 30] −1.58606 −1.586494973 −1.588004950 2000 −1.58807411
N
Present method
5 7
−1.588071834 −1.588071347
Fig. 1. Thomas–Fermi graph obtained by the present method.
As in a typical tau method and using Eq. (8) we can write
⟨Res(x), Rk (x)⟩w =
∞
Res(x)Rk (x)w(x)dx = 0,
k = 0, 1, . . . , N − 2.
(24)
0
Now, a pseudospectral method is defined by applying Eq. (21) on (24) to generate (N − 1) algebraic equations as follows: N
Res(xj )Rk (xj )wj = 0,
k = 0, 1, . . . , N − 2.
(25)
j=0
In addition, using Eqs. (12)–(14) the boundary conditions in Eq. (22) can be approximated as AT e1 = 1,
AT e2 = 0.
(26)
Solving the system of N + 1 nonlinear equations in Eqs. (25) and (26) using Newton’s method for the unknown coefficients aj and substituting the obtained results in Eq. (12) the values of y(x) can be obtained. The initial slope y′ (0) is difficult to compute by any means and plays an important role in determining many physical properties of the Thomas–Fermi atom. It determines the energy of a neutral atom in the Thomas–Fermi approximation:
E=
6 7
4π 3
32
7
Z 3 y′ (0),
(27)
where Z is the nuclear charge. Kobayashi [40] provided a highly accurate numerical value of the initial slope as y′ (0) = −1.588071. The approximations of y′ (0) computed by the present method and their relative errors are shown in Table 1. Obviously, this method is convergent by increasing the number of collocation points and obtaining suitable L. The comparison of the initial slope y′ (0) calculated by the present work with values obtained by Liao [31], Khan [32], Yao [41] and Zhu et al. [35] is given in Table 2, which shows that the present solution is highly accurate. Table 3 shows the approximations of y(x) obtained by the method proposed in this paper for N = 7 and L = 0.0958885, and those obtained by Khan [32] and Liao [42]. Fig. 1 shows the resulting graph of Thomas–Fermi for N = 7 and L = 0.0958885 which tends to zero as x increases by the boundary condition limx→∞ y(x) = 0.
84
A. Kılıçman et al. / Journal of Computational and Applied Mathematics 257 (2014) 79–85 Table 3 Approximation of y(x) for present method, [32,42]. x
Khan [32]
Liao [42]
Present method
0.25 0.50 0.75 1.00 1.25 1.50 1.75 2.00 2.25 2.50 2.75 3.00 3.25 3.50 3.75 4.00 4.25 4.50 4.75 5.00 6.00 7.00 8.00 9.00 10.0 15.0 20.0 25.0 50.0 75.0 100.0 1000.0 5000.0
0.776191000 0.615917000 0.505380000 0.423772000 0.362935000 0.314490000 0.275154000 0.242718000 0.215630000 0.192795000 0.173364000 0.156719000 0.142371000 0.129937000 0.119108000 0.109632000 0.101303000 0.093950400 0.087432000 0.081629600 0.063816200 0.051800500 0.043285900 0.037002300 0.032208100 0.019184300 0.013493700 0.010357000 0.004730890 0.003052460 0.002251000 0.000214641 –
0.755202000 0.606987000 0.502347000 0.424008000 0.363202000 0.314778000 0.275451000 0.243009000 0.215895000 0.192984000 0.173441000 0.156633000 0.142070000 0.129370000 0.118229000 0.108404000 0.099697900 0.091948200 0.085021800 0.078807800 0.059423000 0.046097800 0.036587300 0.029590900 0.024314300 0.010805400 0.005784940 0.003473750 0.000632255 0.000218210 0.000100243 0.000000135 –
0.755903399 0.605270502 0.497823942 0.420343948 0.362814756 0.318737461 0.284014958 0.256010764 0.232974191 0.213705386 0.197357837 0.183318729 0.171134229 0.160461449 0.151036695 0.142653971 0.135150082 0.128394110 0.122279834 0.116720187 0.098752518 0.085573204 0.075494762 0.067538722 0.061098873 0.041370727 0.031271686 0.025135426 0.012687078 0.008484835 0.006373709 0.000640105 0.000128070
4. Conclusion The fundamental goal of this paper has been to construct an approximation to the solution of the nonlinear Thomas–Fermi equation in a semi-infinite interval which has singularity at x = 0 and its boundary condition occurs in infinity. In the above discussion, the combination of tau and pseudospectral methods with RSC functions, which have the property of orthogonality, is employed to achieve this goal. Advantages of this method are that we do not reform the problem to a finite domain and with a small N very accurate results are obtained. There is a good agreement between obtained results for y′ (0) and exact value that demonstrates the validity of the present method for this type of problems and gives the method a wider applicability. Comparing the computed results by this method with the others shows that this method provides more accurate and numerically stable solutions than those obtained by other methods. Several possible extensions of this work can be done for further research. Some of these areas include numerical solutions of integral equations defined on infinite domains, infinite horizon optimal control problems and partial boundary value problems over infinite intervals. References [1] B.Y. Guo, J. Shen, Laguerre–Galerkin method for nonlinear partial differential equations on a semi-infinite interval, Numer. Math. 86 (4) (2000) 635–654. [2] J. Shen, Stable and efficient spectral methods in unbounded domains using Laguerre functions, SIAM J. Numer. Anal. 38 (4) (2000) 1113–1133. [3] Y. Maday, B. Pernaud-Thomas, H. Vandeven, Shock-fitting techniques for solving hyperbolic problems with spectral methods, Rech. Aerosp. 6 (1985) 1–9. [4] H.I. Siyyam, Laguerre Tau methods for solving higher-order ordinary differential equations, J. Comput. Anal. Appl. 3 (2) (2001) 173–182. [5] K. Parand, M. Dehghan, A.R. Rezaei, S.M. Ghaderi, An approximation algorithm for the solution of the nonlinear Lane–Emden type equations arising in astrophysics using Hermite functions collocation method, Comput. Phys. Commun. 181 (6) (2010) 1096–1108. [6] J.P. Boyd, Chebyshev and Fourier Spectral Methods, second ed., Springer, Berlin, 2000. [7] M. Maleki, M. Tavassoli Kajani, I. Hashim, A. Kılıçman, K.A.M. Atan, A nonclassical radau collocation method for nonlinear initial-value problems with applications to Lane–Emden type equations, J. Appl. Math. 2012 (2012), Article ID 103205, 13 pages. http://dx.doi.org/10.1155/2012/103205. [8] B.Y. Guo, Jacobi spectral approximations to differential equations on the half line, J. Comput. Math. 18 (1) (2000) 95–112. [9] M. Maleki, I. Hashim, S. Abbasbandy, Analysis of IVPs and BVPs on semi-infinite domains via collocation methods, J. Appl. Math. 2012 (2012), Article ID 696574, 21 pages. http://dx.doi.org/10.1155/2012/696574. [10] P.M. Lima, M.P. Carpentier, Iterative methods for a singular boundary-value problem, J. Comput. Appl. Math. 111 (12) (1999) 173–186. [11] K.Q. Lan, Multiple eigenvalues for singular Hammerstein integral equations with applications to boundary value problems, J. Comput. Appl. Math. 189 (12) (2006) 109–119. [12] H. Allouche, N. Marhraoui, Numerical solution of singular regular boundary value problems by pole detection with qd-algorithm, J. Comput. Appl. Math. 233 (2) (2009) 420–436.
A. Kılıçman et al. / Journal of Computational and Applied Mathematics 257 (2014) 79–85
85
[13] J.P. Boyd, Orthogonal rational functions on a semi-infinite interval, J. Comput. Phys. 70 (1) (1987) 63–88. [14] John P. Boyd, Rational Chebyshev series for the ThomasFermi function: endpoint singularities and spectral methods, J. Comput. Appl. Math. 244 (2013) 90–101. [15] B.Y. Guo, J. Shen, Z.Q. Wang, A rational approximation and its applications to differential equations on the half line, J. Sci. Comput. 15 (2) (2000) 117–147. [16] J.P. Boyd, C. Rangan, P.H. Bucksbaum, Pseudospectral methods on a semi-infinite interval with application to the hydrogen atom: a comparison of the mapped Fourier-sine method with Laguerre series and rational Chebyshev expansions, J. Comput. Phys. 188 (1) (2003) 56–74. [17] K. Parand, M. Razzaghi, Rational Legendre approximation for solving some physical problems on semi-infinite intervals, Phys. Scr. 69 (5) (2004) 353–357. [18] K. Parand, M. Razzaghi, Rational Chebyshev tau method for solving Volterra’s population model, Appl. Math. Comput. 149 (3) (2004) 893–900. [19] K. Parand, M. Razzaghi, Rational Chebyshev tau method for solving higher-order ordinary differential equations, Int. J. Comput. Math. 81 (1) (2004) 73. [20] M. Tavassoli Kajani, F. Ghasemi Tabatabaei, Rational Chebyshev approximations for solving Lane–Emde equation of index m, in: Proceeding of the International Conference on Computational and Applied Mathematics, Bangkok, Thailand, 29–31 March, 2011, pp. 840–844. [21] M. Dadkhah Tirani, F. Ghasemi Tabatabaei, M. Tavassoli Kajani, Rational second (third) kind Chebyshev approximations for solving Volterras population model, in: Proceeding of the International Conference on Computational and Applied Mathematics, Bangkok, Thailand, 29–31 March, 2011, pp. 835–839. [22] M. Tavassoli Kajani, S. Vahdati, Zulkifly Abbas, M. Maleki, Application of rational second kind Chebyshev functions for system of integrodifferential equations on semi-infinite intervals, J. Appl. Math. 2012 (2012), Article ID 803503, 11 pages. http://dx.doi.org/10.1155/2012/803503. [23] H.T. Davis, Introduction to Nonlinear Differential and Integral Equations, Dover, New York, 1962. [24] S. Chandrasekhar, Introduction to the Study of Stellar Structure, Dover, New York, 1967. [25] E. Di Grezia, S. Esposito, Fermi, Majorana and the statistical model of atoms, Found. Phys. 34 (2004) 1431–1450. [26] C.M. Bender, K.A. Milton, S.S. Pinsky, L.M. Simmons Jr., A new perturbative approach to nonlinear problems, J. Math. Phys. 30 (7) (1989) 1447–1455. [27] B.J. Laurenzi, An analytic solution to the Thomas–Fermi equation, J. Math. Phys. 31 (10) (1990) 2535–2537. [28] A. Cedillo, A perturbative approach to the Thomas–Fermi equation in terms of the density, J. Math. Phys. 34 (7) (1993) 2713–2717. [29] G. Adomian, Solution of the Thomas–Fermi equation, Appl. Math. Lett. 11 (3) (1998) 131–133. [30] A. Wazwaz, The modified decomposition method and Padé approximants for solving the Thomas–Fermi equation, Appl. Math. Comput. 105 (1) (1999) 11–19. [31] S. Liao, An explicit analytic solution to the Thomas–Fermi equation, Appl. Math. Comput. 144 (23) (2003) 495–506. [32] H. Khan, H. Xu, Series solution to the Thomas–Fermi equation, Phys. Lett. A 365 (12) (2007) 111–115. [33] V.B. Mandelzweig, F. Tabakin, Quasilinearization approach to nonlinear problems in physics with application to nonlinear ODEs, Comput. Phys. Commun. 141 (2) (2001) 268–281. [34] J.I. Ramos, Piecewise quasilinearization techniques for singular boundary-value problems, Comput. Phys. Commun. 158 (1) (2004) 12–25. [35] S. Zhu, H. Zhu, Q. Wu, Y. Khan, An adaptive algorithm for the Thomas–Fermi equation, Numer. Algorithms 59 (2012) 359–372. [36] K. Parand, M. Dehghan, A. Pirkhedri, The Sinc-collocation method for solving the Thomas–Fermi equation, J. Comput. Appl. Math. 237 (2013) 244–252. [37] M. Abramowitz, I.A. Stegun, Handbook of Mathematical Functions, Dover, New York, 1972. 10th printing with corrections. [38] J.P. Boyd, The optimization of convergence for Chebyshev polynomial methods in an unbounded domain, J. Comput. Phys. 45 (1) (1982) 43–79. [39] J.P. Boyd, Orthogonal rational functions on a semi-infinite interval, J. Comput. Phys. 70 (1987) 63–88. [40] S. Kobayashi, T. Matsukuma, S. Nagi, K. Umeda, Accurate value of the initial slope of the ordinary TF function, J. Phys. Soc. Japan 10 (1955) 759–762. [41] B. Yao, A series solution to the Thomas–Fermi equation, Appl. Math. Comput. 203 (2008) 396–401. [42] S. Liao, Beyond Perturbation-Introduction to the Homotopy Analysis Method, Chapman & Hall/CRC, Boca Raton, 2003.