- Email: [email protected]
The Sinc-collocation method for solving the Thomas–Fermi equation
Journal of Computational and Applied Mathematics 237 (2013) 244–252
Contents lists available at SciVerse ScienceDirect
Journal of Computational and Applied Mathematics journal homepage: www.elsevier.com/locate/cam
The Sinc-collocation method for solving the Thomas–Fermi equation K. Parand a , Mehdi Dehghan b,∗ , A. Pirkhedri c a
Department of Computer Sciences, Shahid Beheshti University, G.C. Tehran, Iran
b
Department of Applied Mathematics, Faculty of Mathematics and Computer Science, Amirkabir University of Technology, No.424, Hafez Avenue, Tehran 15914, Iran c Department of Computer Sciences, Islamic Azad University Branch of Marivan, Marivan, Iran
article
info
Article history: Received 8 November 2009 Received in revised form 22 December 2011 MSC: 33C45 65L60 34A34
abstract A numerical technique for solving nonlinear ordinary differential equations on a semiinfinite interval is presented. We solve the Thomas–Fermi equation by the Sinc-Collocation method that converges to the solution at an exponential rate. This method is utilized to reduce the nonlinear ordinary differential equation to some algebraic equations. This method is easy to implement and yields very accurate results. © 2012 Elsevier B.V. All rights reserved.
Keywords: Thomas–Fermi equation Sinc functions Collocation method Nonlinear ODE Semi-infinite
1. Introduction Recently, spectral methods [1,2] have been successfully applied in the approximation of boundary value problems defined in the semi-infinite domains. We can apply different approaches using spectral methods to solve problems in semi-infinite domains. Various spectral methods for treating semi-infinite domains have been proposed. One approach is using Laguerre polynomials [3–6]. A second approach is reformulating the original problem in the semi-infinite domain to a singular problem in a bounded domain by variable transformation and then using the Jacobi polynomials to approximate the resulting singular problem [7–9]. A third approach of the spectral method is based on rational orthogonal functions, for example, Christov [10] and Boyd [11,12] developed some spectral methods on unbounded intervals by using mutually orthogonal systems of rational functions. Boyd [12] defined a new spectral basis, named rational Chebyshev functions on the semiinfinite interval, by mapping it to the Chebyshev polynomials [13]. Guo et al. [14] proposed and analyzed a set of Legendre rational functions which are mutually orthogonal in Lχ 2 (0, ∞) with a non-uniform weight function χ (x) = (x + 1)−2 . Parand et al. [15–19] applied the rational Chebyshev, rational Legendre and rational scaled generalized Laguerre functions with tau [20,21] and collocation methods to solve nonlinear ordinary differential equations on semi-infinite intervals. A fourth approach is replacing the semi-infinite domain with [0, L] interval by choosing L, sufficiently large, this method is named as the domain truncation [22]. The Sinc-collocation and Sinc–Galerkin methods for solving differential equations are based on Sinc approximation [23]. In this paper, we investigate the Sinc-collocation [24,25] method on the half line by using Sinc functions. The Sinc-collocation
∗
Corresponding author. Fax: +98 21 66497930. E-mail addresses: [email protected] (K. Parand), [email protected], [email protected] (M. Dehghan), [email protected] (A. Pirkhedri). 0377-0427/$ – see front matter © 2012 Elsevier B.V. All rights reserved. doi:10.1016/j.cam.2012.08.001
K. Parand et al. / Journal of Computational and Applied Mathematics 237 (2013) 244–252
245
method for the numerical solution of the initial value problems is developed in [26] and it is approved that it converges to the solution at an exponential rate. In [27,28] the Sinc-collocation method is discussed for solving the Blasius and Planar Coulomb Schrödinger equations on the half line. The Sinc–Galerkin method was applied widely by [29–31]. The author of [29] applied this method to solve certain class of singular two-point boundary value problems and expressed the exact solution of the differential equations via the use of Green’s functions as an integral type. Saadatmandi et al. [30] applied the Sinc–Galerkin method for solving nonlinear two-point boundary value problems [32] for the second order differential equations. 2. The Thomas–Fermi equation The Thomas–Fermi theory establishes a functional relation between the energy of an electronic system, E, and the electronic density, ρ , namely, E [ρ] =
9 10B
ρ(r )dτ +
1
2
ρ(r )ρ(r ′ ) ′ dτ dτ + |r − r ′ |
ρ(r )ν(r )dτ ,
(2.1)
2
where ν(r ) is the external potential and B = 3(3π 2 )− 3 . The density can be obtained by minimizing the energy functional with respect to ρ , subject to the normalization restriction ρ(r )dτ = N where N is a number of electrons [33]. Then the density must satisfy the following integral equation: 3 2B
2
ρ(r ) 3 +
ρ(r ′ ) dτ ′ + ν(r ) = µ, |r − r ′ |
(2.2)
where µ is the Lagrange multiplier related to the normalization restriction. Poisson’s equation can be used to remove the density, and a change of variables leads to the Thomas–Fermi equation d2 y dx2
1
3
= √ y 2 (x),
(2.3)
x
with boundary conditions as follows: y(0) = 1,
lim y(x) = 0.
x→∞
(2.4)
This equation describes the charge density in atoms of high atomic number and appears in the problem of determining the effective of nuclear charge in heavy atoms [34,35]. It is useful for calculating form-factors and for obtaining effective potentials which can be used as initial trial potentials in self-consistent field calculations. It is also applicable to the study of nucleons in the atom and electrons in the metal. If we use the Runge–Kutta method [36] to solve Eq. (2.3) it can be found numerically with great difficulty, in that, to integrate from x = 0 we must assume a value for y′ (0), if y′ (0) is chosen too small, the solution will cross below the x axis at some finite values of x and becomes complex and if y′ (0) is chosen too large the solution will eventually become singular at some finite values of x [37]. It is long known that the solution of this equation is very sensitive to a value of the first derivative at zero which ensures smooth and monotonic decay from y(0) = 1 to y(∞) = 0 as demanded by boundary conditions [38]. Cedillo [33] wrote the Thomas–Fermi equation in terms of density and then the δ -expansion was employed to obtain an absolutely convergent series of equations. Mandelzweig and Tabakin [39] employed a quasilinearization method to solve the Thomas–Fermi equation written as y′′ = x−1/2 y3/2 for 0 ≤ x ≤ 40 by means of a quasilinearization method starting with the initial guess y0 (x) = 1, and showed that the convergence starts at the boundaries and expands with each iteration to a wider range of values of x. In this paper we approximate y(x) by the Sinc-collocation method because first, it is easy to apply and numerically achieve exponential convergence, second, because of singularity in this equation Sinc functions can handle this problem, third, the limit of the Sinc function at infinity is zero and thus one of the boundary conditions in Thomas–Fermi equation implicitly becomes true. The organization of this paper is as follows: in Section 3, we explain the formulation of Sinc functions required for our subsequent development. In Section 4, we summarize the application of the method of Sinc functions for solving the Thomas–Fermi equation and comparing it with the existing methods in the literature. Section 5 is devoted to summary and conclusions. 3. Sinc function properties The Sinc function is defined for all x ∈ R by
sinc(x) =
sin(π x)
1,
πx
,
x ̸= 0, x = 0.
(3.1)
246
K. Parand et al. / Journal of Computational and Applied Mathematics 237 (2013) 244–252
For each integer k and the mesh size h, the Sinc basis functions are defined on R by Sk (h, x) ≡ sinc
x − kh
π sin( h (x − kh)) , π = (x − kh) h 1,
h
x ̸= kh,
(3.2)
x = kh.
Approximations can be constructed for infinite, semi-infinite and finite intervals. Define the function w = φ(x) = ln (sinh(x)), which is a conformal mapping from DE , the eye-shaped spatial domain in the z-plane, onto the infinite strip, DS in the w -plane, where
DE = z ∈ R : |arg(sinh(z ))| < d ≤
π 2
DS = w ∈ R : w = x + iy, |y| < d ≤
, π 2
.
The basis functions on (0, +∞) are taken to be the composite translated Sinc functions Sk (x) ≡ S (k, h) ◦ φ(x) = sinc
φ(x) − kh h
,
(3.3)
where S (k, h) ◦ φ(x) is defined by S (k, h)(φ(x)). The inverse map of ω = φ(x) is x = φ −1 (w) = ln(ew +
1 + e2w ).
(3.4)
Thus we may define the inverse images of the real line and of the evenly spaced nodes {kh}
k=+∞ k=−∞
as
Γ = {φ −1 (w) : w ∈ DS } = (0, +∞),
(3.5)
and xk = φ −1 (kh) = ln(ekh +
1 + e2kh ),
k = 0, ±1, ±2, . . . .
(3.6)
Let w(x) denotes a non-negative, integrable, real-valued function over the interval Γ , we define L2w (Γ ) = {v : Γ → R | v is measurable and ∥v∥w < ∞},
(3.7)
where ∞
|v(x)| w(x)dx 2
∥v∥w =
21
,
(3.8)
0
is the norm induced by the inner product of the space L2w (Γ ),
⟨u, v⟩w =
∞
u(x)v(x)w(x)dx.
(3.9)
0
Thus {Sk (x)}k∈Z with constant h denote a system which is mutually orthogonal under Eq. (3.9), i.e.,
⟨Skn (x), Skm (x)⟩w(x) = hδnm ,
(3.10)
where w(x) = coth(x) and δnm is the Kronecker delta function. This system is complete in Lw (Γ ), and for any function f ∈ L2w (Γ ) the following expansion holds 2
f (x) ∼ =
N
fk Sk (x),
(3.11)
k=−N
with fk =
⟨f (x), Sk (x)⟩w(x) . ∥Sk (x)∥2w(x)
(3.12)
The fk = f (xk ) are the discrete expansion coefficients associated with the family {Sk (x)} and the mesh size is given by
h=
2π d
αN
,
where N is suitably chosen, α and d depend [40] on the asymptotic behavior of f (x).
K. Parand et al. / Journal of Computational and Applied Mathematics 237 (2013) 244–252
247
Definition. Let H 2 (DE ) be the class of functions f which are analytic in DE , satisfy
|f (z )|dz → 0,
φ −1 (x+P )
x → ±∞,
where P = {iy : |y| < d ≤ π2 }, and on the boundary of DE (denoted ∂ D ) satisfy N (f ) =
∂D
|f (z )dz | < ∞.
Interpolation for the function in H 2 (DE ) is defined in the following theorem. Theorem ([41]). Assume that f φ ′ ∈ H 2 (DE ) then for all x ∈ Γ
∞ N (f φ ′ ) −π d/h N (f φ ′ ) E (f , h)(x) = f (x) − ≤2 e . f (kh)S (k, h) ◦ φ(x) ≤ 2π d sinh(π d/h) πd k=−∞ Moreover, if |f (x)| ≤ Ce−α|φ(x)| , x ∈ Γ , for some positive constants C and α , and if the selection h =
πd αN
≤
2π d , ln(2)
then
N √ √ f (kh)S (k, h) ◦ φ(x) ≤ C2 Ne(− π dα N ) , f (x) − k=−N where C2 depends only on f , d and α . 3.1. The matrix representation of the derivatives of Sinc basis functions at nodal points The nth derivative of function f at the nodal point at xk is denoted by dr 1 (r ) δk,j = [S (k, h) ◦ φ(x)]|x=xj . r h dφ r
(3.13) (r )
The expressions in Eq. (3.13) for each k and j can be stored in a matrix I (r ) = [δk,j ] where (0)
I (0) = [δk,j ],
δk(0,j) = [S (k, h) ◦ φ(x)]|x=xj =
(1)
δk(1,j)
(2)
δk(2,j)
I (1) = [δk,j ],
I (2) = [δk,j ],
1, 0,
k = j, k ̸= j,
k = j, 0, = h [S (k, h) ◦ φ(x)]|x=xj = (−1)j−k , k ̸= j, dφ j−k −π 2 2 , k = j, d = h2 2 [S (k, h) ◦ φ(x)]|x=xj = −32(−1)j−k dφ , k ̸= j. (j − k)2 d
(3.14)
(3.15)
(3.16)
The size of the matrix is determined by the limits on row and column indices k, j. 4. Solving the Thomas–Fermi equation For the boundary conditions in Eq. (2.4), the Sinc basis functions in Eq. (3.3) do not have a derivative and when the value of x tends to zero, thus we modify the Sinc basis functions as x x2 + 1
Sk (x).
(4.1)
In order to discretize Eq. (2.3) by using Sinc-collocation, first of all, we approximate y(x) as yN (x) = uN (x) + p(x),
(4.2)
where uN (x) =
N k=−N
ck
xSk (x) x2 + 1
.
(4.3)
248
K. Parand et al. / Journal of Computational and Applied Mathematics 237 (2013) 244–252
In order to approximate the solution of Eq. (2.3) with boundary conditions, we construct a polynomial p(x) that satisfies Eq. (2.4), this polynomial is given by p(x) =
λ x+λ
.
(4.4)
In Eq. (4.4), λ is constant to be determined. It is noted that the approximate solution yN (x), satisfies the boundary conditions in Eq. (2.4), since lim p(x) = 1,
lim uN (x) = 0,
x→0
x→0
lim p(x) = 0,
(4.5)
lim uN (x) = 0.
x→∞
x→∞
We refer the interested reader to [42,43] which address these issues. Now we display equations for yN (x), y′N (x) and y′′N (x). Using Eqs. (3.14)–(3.16) and (4.3) with substitution the result at the Sinc points:
xj = ln ejh +
1 + ln(e2jh ) ,
j = −N − 1 . . . N ,
(4.6)
we get yN (xj ) =
λ xj + λ
+
cj xj x2j
+1
,
(4.7)
N λ ck + yN (xj ) = − (xj + λ)2 k=−N
1
′
1 + x2j
−
2x2j
(1 + x2j )2
δk(0,j)
+
xj φ ′ (xj )
1 + x2j
δk(1,j)
,
(4.8)
and
8x3j −6xj ck + yN (xj ) = + δ (0) (xj + λ)3 k=−N (1 + x2j )2 (1 + x2j )3 k,j 2 4x2j φ ′ (xj ) 2φ ′ (xj ) xj φ ′ (xj ) xj φ ′′ (xj ) (1) (2) + − + δk,j + δk,j . 1 + x2j (1 + x2j )2 1 + x2j 1 + x2j 2λ
′′
N
(4.9)
4.1. The matrix representation of the Sinc-collocation method for solving the equation Suppose that
8x3 −6x + , (1 + x2 )2 (1 + x2 )3 2 xφ ′ (x) 1 Y(x) = , Z(x) = √ , 1 + x2 x f (x )
W(x) =
−N −1
D (f) =
..
.
X(x) =
T(x) = 0
2φ ′ (x) 1 + x2 x
x2 + 1
−
..
.
0
(1 + x2 )2
+
xφ ′′ (x) 1 + x2
,
f (xN )
,
,
f (x0 )
4x2 φ ′ (x)
0
c−N
⃗ = . C . .
,
(x−N −1 + λ) .. . λ (xN + λ)
⃗ = Λ
cN
λ
.
(4.10)
Substituting Eqs. (4.7) and (4.9) in Eq. (2.3) and subject to Eq. (4.2) we obtain:
⃗+ D (W)I (0) + D (X)I (1) + D (Y)I (2) C
∂2 3 ⃗ + D (Z)(Λ ⃗ + D (T)I (0) C ⃗ ) 2 = 0. Λ 2 ∂x
(4.11)
N The 2N + 1 coefficients {ck }kk= =−N and the unknown λ determined by solving Eq. (4.11) and evaluating the result at the Sinc points:
xj = ln ejh +
1 + ln(e2jh ) ,
j = −N − 1 . . . N .
(4.12)
K. Parand et al. / Journal of Computational and Applied Mathematics 237 (2013) 244–252
249
When −N − 1 ≤ j ≤ N , −N ≤ k ≤ N we remove the last column, the first and last rows of I (0) , I (1) and I (2) to arrive at the non-square matrices:
... .. . ...
−1
0
I (0) = 1 0
I (1)
=
0
0 , 1 1 2
..
.
..
.
1
..
.
..
.
..
.
..
.
...
..
.
1 2
(−1)2N +1 2N
I
(2)
...
0
−
(4.13)
−
2
π2 − 3 2 = 2 − 2 2 .. . −2(−1)2N +1 (2N + 1)
2
2 22
(−1)2N +1 2N + 1 .. . 1 , 2 −1
(4.14)
0
...
..
.
..
.
..
.
..
.
..
.
..
.
..
.
..
.
..
.
..
.
−2(−1)2N +1 (2N + 2)1 .. . 2 − 2 2 . 2 2 2 π −
(4.15)
3
Eq. (4.11) gives 2N + 2 nonlinear algebraic equations which can be solved for the unknown coefficients ck and λ by using the well known Newton’s method by Maple programming and we use
ck = 0, λ = 0,
k = −N . . . N ,
as starting points to obtain convergence of the method, consequently, y(x) given in Eq. (2.3) can be calculated. One measure of the rapidity of convergence of the procedure is provided by calculation of the value of the initial slope y′ (0) of the Thomas–Fermi potential [44]. 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),
(4.16)
where Z is the nuclear charge. Kobayashi et al. [45] provided a highly accurate numerical solution of the Thomas–Fermi as y′ (0) = −1.588071. Table 1 shows the comparison of y′ (0), between Padé approximation used by [46,47] and the y′N (0) = − λ1 obtained by Sinc-collocation method and choosing α = 1, d = (h =
2π 3
, N = 32 and the mesh size is given by El-Gamel et al. [40]
2π d ). αN
The approximations of y(x) obtained by Sinc-collocation method and the methods introduced in [46–48] are listed in Table 2 by choosing α = 1, d = 23π , N = 32 and h = √2π . 3N
250
K. Parand et al. / Journal of Computational and Applied Mathematics 237 (2013) 244–252 Table 1 Results for the value of y′ (0) obtained by Sinc-collocation method for α = 1, d = 23π , N = 32 and h = √2π . 3N
N
Sinc-collocation method
h = √2π
Padé
Liao [46]
Khan [47]
4 8 16 32
−1.552409174 −1.581426241 −1.587396189 −1.588070339
1.813799365 1.282549831 0.906899682 0.641274915
[10, 10] [20, 20] [30, 30] −
−1.51508 −1.58606 −1.58281 −
−1.573824678 −1.586494973 −1.582901907 −
3N
Table 2 Approximations of y(x) for Sinc-collocation method, by choosing α = 1, d = 23π , N = 32 and h = √2π . 3N
x
Sinc-collocation method
Liao [46]
Khan [47]
Bender [48]
1.00 2.00 4.00 6.00 8.00 10.00 50.00
0.4240642728 0.2430282344 0.1084800242 0.0594010318 0.0365353166 0.0243069553 0.0006213281
0.424008000 0.243009000 0.108404000 0.059423000 0.036587300 0.024314300 0.000632255
0.423772000 0.242718000 0.109632000 0.063816200 0.043285900 0.032208100 0.004730890
0.4240 0.2430 0.1084 0.0594 0.0365 0.0243 0.0006
Fig. 1. Thomas–Fermi graph obtained by Sinc-collocation method (solid line) in comparison with numerical solution (solid circles) [48] for N = 32.
Fig. 1 shows the result graph of Thomas–Fermi for N = 32 and it compares the new results with numerical results in [48]. Fig. 2 shows the logarithmic scale for the relative errors versus N. We refer the interested reader to [49,50] for more applications of spectral techniques. 5. Summary and conclusions The Sinc-collocation method is used to solve the Thomas–Fermi equation which is defined in the semi-infinite interval and has singularity at x = 0 and its boundary condition occurs at infinity. The numerical results demonstrate the reliability and efficiency of using the new method to solve such problems. The results of the present method for this type of problem clearly indicate that the new method is accurate even when singularity occurs at the boundary. Acknowledgments The authors would like to thank the referees for their comments and suggestions that have improved the paper. Also they would like to thank professor Taketomo Mitsui, the principal editor and professor Masaaki Sugihara, the advisory editor of
K. Parand et al. / Journal of Computational and Applied Mathematics 237 (2013) 244–252
251
Fig. 2. Logarithmic scale of the relative errors versus N for the value of y′ (0).
the journal for their useful comments and for managing the review process for this paper during 8 November 2009–1 August 2012. The first author (KP) would like to thank Shahid Beheshti University for the awarded grant. References [1] M. Dehghan, F. Fakhar-Izadi, The spectral collocation method with three different bases for solving a nonlinear partial differential equation arising in modeling of nonlinear waves, Math. Comput. Modelling 53 (2011) 1865–1877. [2] F. Fakhar-Izadi, M. Dehghan, The spectral methods for parabolic Volterra integro–differential equations, J. Comput. Appl. Math. 235 (2011) 4032–4046. [3] Y. Maday, B. Pernaud-Thomas, H. Vandeven, Reappraisal of Laguerre type spectral methods, La Recherche Aerospatiale 6 (1985) 13–35. [4] D. Funaro, Computational aspects of pseudospectral Laguerre approximations, Appl. Numer. Math. 6 (1990) 447–457. [5] B.Y. Guo, J. Shen, Laguerre-Galerkin method for nonlinear partial differential equations on a semi-infinite interval, Numer. Math. 86 (2000) 635–654. [6] J. Shen, Stable and efficient spectral methods in unbounded domains using Laguerre functions, SIAM J. Numer. Anal. 38 (2000) 1113–1133. [7] B.Y. Guo, Gegenbauer approximation and its applications to differential equations on the whole line, J. Math. Anal. Appl. 226 (1998) 180–206. [8] B.Y. Guo, Jacobi spectral approximation and its applications to differential equations on the half line, J. Comput. Math. 18 (2000) 95–112. [9] B.Y. Guo, Jacobi approximations in certain Hilbert spaces and their applications to singular differential equations, J. Math. Anal. Appl. 243 (2000) 373–408. [10] C.I. Christov, A complete orthogonal system of functions in L2 (−∞, ∞) space, SIAM J. Appl. Math. 42 (1982) 1337–1344. [11] J.P. Boyd, Spectral methods using rational basis functions on an infinite interval, J. Comput. Phys. 69 (1987) 112–142. [12] J.P. Boyd, Orthogonal rational functions on a semi–infinite interval, J. Comput. Phys. 70 (1987) 63–68. [13] M.R. Eslahchi, M. Dehghan, S. Amani, The third and fourth kinds of Chebyshev polynomials and best uniform approximation, Math. Comput. Modelling 55 (2012) 1746–1762. [14] 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 (2000) 117–147. [15] K. Parand, M. Razzaghi, Rational Chebyshev tau method for solving Volterra’s population model, Appl. Math. Comput. 149 (2004) 893–900. [16] K. Parand, M. Razzaghi, Rational Chebyshev tau method for solving higher-order ordinary differential equations, Int. J. Comput. Math. 81 (2004) 73–80. [17] K. Parand, M. Razzaghi, Rational Legendre approximation for solving some physical problems on semi-infinite intervals, Phys. Scr. 69 (2004) 353–357. [18] K. Parand, M. Shahini, Rational Chebyshev pseudospectral approach for solving Thomas–Fermi equation, Phys. Lett. A 373 (2009) 210–213. [19] K. Parand, A. Taghavi, Rational scaled generalized Laguerre function collocation method for solving the Blasius equation, J. Comput. Applied. Math. 233 (2009) 980–989. [20] A. Saadatmandi, M. Dehghan, A tau approach for solution of the space fractional diffusion equation, Comput. Math. Applic. 62 (2011) 1135–1142. [21] F. Shakeri, M. Dehghan, A hybrid Legendre tau method for the solution of a class of nonlinear wave equations with nonlinear dissipative terms, Numer. Methods Partial Differential Equations 27 (2011) 1055–1071. [22] J.P. Boyd, Chebyshev and Fourier Spectral Methods, second ed., Dover, New York, 2000. [23] T. Okayama, T. Matsuo, M. Sugihara, Sinc–collocation methods for weakly singular Fredholm integral equation of the second kind, J. Comput. Appl. Math. 234 (2010) 1211–1227. [24] A. Saadatmandi, M. Dehghan, M.R. Azizi, The Sinc-Legendre collocation method for a class of fractional convection-diffusion equations with variable coefficients, Commun. Nonlinear Sci. Numer. Simul. 17 (2012) 4125–4136. [25] A. Saadatmandi, M. Dehghan, The use of Sinc–collocation method for solving multi–point boundary value problems, Commun. Nonlinear Sci. Numer. Simul. 17 (2012) 593–601. [26] S. Timothy Carlson, Jack Dockery, John Lund, A Sinc–collocation method for initial value problems, Math. Comp. 217 (1997) 215–235. [27] K. Parand, M. Dehghan, A. Pirkhedri, Sinc-collocation methods for solving the Blasius equation, Phys. Letl. A 373 (2009) 4060–4065. [28] V.G. Koures, F.E. Harris, Sinc collocation in quantum chemistry: Solving the planar coulomb Schrödinger equation, International Journal of Quantum Chemistry: Quantum Chemistry Symposium 30 (1996) 1311–1318. [29] Kamel Al-Khaled, Theory and computation in singular boundary value problems, Chaos Solitons Fractals 33 (2007) 678–684. [30] A. Saadatmandi, M. Razzaghi, M. Dehghan, Sinc–Galerkin solution for nonlinear two-point boundary value problems with application to chemical reactor theory, Math. Comput. Modelling 42 (2005) 1237–1244. [31] F. Stenger, Sinc–Galerkin method of solution of boundary value problems, Math. Comp. 33 (1979) 85–109.
252
K. Parand et al. / Journal of Computational and Applied Mathematics 237 (2013) 244–252
[32] M. Tatari, M. Dehghan, An efficient method for solving multi–point boundary value problems and applications in physics, Journal of Vibration and Control 18 (2012) 1116–1124. [33] A. Cedillo, A perturbative approach to the Thomas–Fermi equation in terms of the density, J. Math. Phys. 34 (1993) 2713–2717. [34] H.T. Davis, Introduction to Nonlinear Differential and Integral Equations, Dover, New York, 1962. [35] S. Chandrasekhar, Introduction to the Study of Stellar Structure, Dover, New York, 1967. [36] M. Dehghan, R. Jazlanian, A fourth–order central Runge–Kutta scheme for hyperbolic conservation laws, Numer. Methods Partial Differential Equations 26 (2010) 1675–1692. [37] C.M. Bender, K.A. Milton, S.S. Pinsky, L.M. Simmons Jr, A new perturbative approach to nonlinear problems, J. Math. Phys. 30 (1989) 1447–1455. [38] H.A. Bethe, R.W. Jackiw, Intermediate Quantum Mechanics, W.A. Benjamin Inc, New York, 1968. [39] V.B. Mandelzweig, F. Tabakin, Quasilinearization approach to nonlinear problems in physics with application to nonlinear ODEs, Comput. Phys. Commun. 141 (2001) 268–281. [40] M. El-Gamel, S.H. Behiry, H. Hashish, Numerical method for the solution of special nonlinear fourth-order boundary value problems, Appl. Math. Comput 145 (2003) 717–734. [41] F. Stenger, Integration formulas via the trapezoidal formula, Journal of Institute Mathematics and its Applications 12 (1973) 103–114. [42] J. Lund, K. Bowers, Sinc Methods for Quadrature and Differential Equations, SIAM, Philadelphia, 1992. [43] F. Stenger, Numerical Methods Based on Sinc and Analytic Functions, Springer-Verlag, New York, 1993. [44] B.J. Laurenzi, An analytic solution to the Thomas–Fermi equation, J. Math. Phys. 10 (1990) 2535–2537. [45] S. Kobayashi, T. Matsukuma, S. Nagi, K. Umeda, Accurate value of the initial slope of the ordinary T-F function, J. Phys. Soc. Jpn. 10 (1955) 759–762. [46] S.J. Liao, An explicit analytic solution to the Thomas–Fermi equation, Appl. Math. Comput. 144 (2003) 495–506. [47] H. Khan, H. Xu, Series solution to the Thomas–Fermi equation, Phys. Lett. A 365 (2007) 111–115. [48] C.M. Bender, S.A. Orszag, Advanced Mathematical Methods for Scientists and Engineers, McGraw-Hill, New York, 1978. [49] J.P. Boyd, The optimization of convergence for Chebyshev polynomial methods in an unbounded domain, J. Comput. Phys. 45 (1982) 43–79. [50] M. Shamsi, M. Dehghan, Determination of a control function in three-dimensional parabolic equations by Legendre pseudospectral method, Numer. Methods Partial Differential Equations 28 (2012) 74–93.
Contents lists available at SciVerse ScienceDirect
Journal of Computational and Applied Mathematics journal homepage: www.elsevier.com/locate/cam
The Sinc-collocation method for solving the Thomas–Fermi equation K. Parand a , Mehdi Dehghan b,∗ , A. Pirkhedri c a
Department of Computer Sciences, Shahid Beheshti University, G.C. Tehran, Iran
b
Department of Applied Mathematics, Faculty of Mathematics and Computer Science, Amirkabir University of Technology, No.424, Hafez Avenue, Tehran 15914, Iran c Department of Computer Sciences, Islamic Azad University Branch of Marivan, Marivan, Iran
article
info
Article history: Received 8 November 2009 Received in revised form 22 December 2011 MSC: 33C45 65L60 34A34
abstract A numerical technique for solving nonlinear ordinary differential equations on a semiinfinite interval is presented. We solve the Thomas–Fermi equation by the Sinc-Collocation method that converges to the solution at an exponential rate. This method is utilized to reduce the nonlinear ordinary differential equation to some algebraic equations. This method is easy to implement and yields very accurate results. © 2012 Elsevier B.V. All rights reserved.
Keywords: Thomas–Fermi equation Sinc functions Collocation method Nonlinear ODE Semi-infinite
1. Introduction Recently, spectral methods [1,2] have been successfully applied in the approximation of boundary value problems defined in the semi-infinite domains. We can apply different approaches using spectral methods to solve problems in semi-infinite domains. Various spectral methods for treating semi-infinite domains have been proposed. One approach is using Laguerre polynomials [3–6]. A second approach is reformulating the original problem in the semi-infinite domain to a singular problem in a bounded domain by variable transformation and then using the Jacobi polynomials to approximate the resulting singular problem [7–9]. A third approach of the spectral method is based on rational orthogonal functions, for example, Christov [10] and Boyd [11,12] developed some spectral methods on unbounded intervals by using mutually orthogonal systems of rational functions. Boyd [12] defined a new spectral basis, named rational Chebyshev functions on the semiinfinite interval, by mapping it to the Chebyshev polynomials [13]. Guo et al. [14] proposed and analyzed a set of Legendre rational functions which are mutually orthogonal in Lχ 2 (0, ∞) with a non-uniform weight function χ (x) = (x + 1)−2 . Parand et al. [15–19] applied the rational Chebyshev, rational Legendre and rational scaled generalized Laguerre functions with tau [20,21] and collocation methods to solve nonlinear ordinary differential equations on semi-infinite intervals. A fourth approach is replacing the semi-infinite domain with [0, L] interval by choosing L, sufficiently large, this method is named as the domain truncation [22]. The Sinc-collocation and Sinc–Galerkin methods for solving differential equations are based on Sinc approximation [23]. In this paper, we investigate the Sinc-collocation [24,25] method on the half line by using Sinc functions. The Sinc-collocation
∗
Corresponding author. Fax: +98 21 66497930. E-mail addresses: [email protected] (K. Parand), [email protected], [email protected] (M. Dehghan), [email protected] (A. Pirkhedri). 0377-0427/$ – see front matter © 2012 Elsevier B.V. All rights reserved. doi:10.1016/j.cam.2012.08.001
K. Parand et al. / Journal of Computational and Applied Mathematics 237 (2013) 244–252
245
method for the numerical solution of the initial value problems is developed in [26] and it is approved that it converges to the solution at an exponential rate. In [27,28] the Sinc-collocation method is discussed for solving the Blasius and Planar Coulomb Schrödinger equations on the half line. The Sinc–Galerkin method was applied widely by [29–31]. The author of [29] applied this method to solve certain class of singular two-point boundary value problems and expressed the exact solution of the differential equations via the use of Green’s functions as an integral type. Saadatmandi et al. [30] applied the Sinc–Galerkin method for solving nonlinear two-point boundary value problems [32] for the second order differential equations. 2. The Thomas–Fermi equation The Thomas–Fermi theory establishes a functional relation between the energy of an electronic system, E, and the electronic density, ρ , namely, E [ρ] =
9 10B
ρ(r )dτ +
1
2
ρ(r )ρ(r ′ ) ′ dτ dτ + |r − r ′ |
ρ(r )ν(r )dτ ,
(2.1)
2
where ν(r ) is the external potential and B = 3(3π 2 )− 3 . The density can be obtained by minimizing the energy functional with respect to ρ , subject to the normalization restriction ρ(r )dτ = N where N is a number of electrons [33]. Then the density must satisfy the following integral equation: 3 2B
2
ρ(r ) 3 +
ρ(r ′ ) dτ ′ + ν(r ) = µ, |r − r ′ |
(2.2)
where µ is the Lagrange multiplier related to the normalization restriction. Poisson’s equation can be used to remove the density, and a change of variables leads to the Thomas–Fermi equation d2 y dx2
1
3
= √ y 2 (x),
(2.3)
x
with boundary conditions as follows: y(0) = 1,
lim y(x) = 0.
x→∞
(2.4)
This equation describes the charge density in atoms of high atomic number and appears in the problem of determining the effective of nuclear charge in heavy atoms [34,35]. It is useful for calculating form-factors and for obtaining effective potentials which can be used as initial trial potentials in self-consistent field calculations. It is also applicable to the study of nucleons in the atom and electrons in the metal. If we use the Runge–Kutta method [36] to solve Eq. (2.3) it can be found numerically with great difficulty, in that, to integrate from x = 0 we must assume a value for y′ (0), if y′ (0) is chosen too small, the solution will cross below the x axis at some finite values of x and becomes complex and if y′ (0) is chosen too large the solution will eventually become singular at some finite values of x [37]. It is long known that the solution of this equation is very sensitive to a value of the first derivative at zero which ensures smooth and monotonic decay from y(0) = 1 to y(∞) = 0 as demanded by boundary conditions [38]. Cedillo [33] wrote the Thomas–Fermi equation in terms of density and then the δ -expansion was employed to obtain an absolutely convergent series of equations. Mandelzweig and Tabakin [39] employed a quasilinearization method to solve the Thomas–Fermi equation written as y′′ = x−1/2 y3/2 for 0 ≤ x ≤ 40 by means of a quasilinearization method starting with the initial guess y0 (x) = 1, and showed that the convergence starts at the boundaries and expands with each iteration to a wider range of values of x. In this paper we approximate y(x) by the Sinc-collocation method because first, it is easy to apply and numerically achieve exponential convergence, second, because of singularity in this equation Sinc functions can handle this problem, third, the limit of the Sinc function at infinity is zero and thus one of the boundary conditions in Thomas–Fermi equation implicitly becomes true. The organization of this paper is as follows: in Section 3, we explain the formulation of Sinc functions required for our subsequent development. In Section 4, we summarize the application of the method of Sinc functions for solving the Thomas–Fermi equation and comparing it with the existing methods in the literature. Section 5 is devoted to summary and conclusions. 3. Sinc function properties The Sinc function is defined for all x ∈ R by
sinc(x) =
sin(π x)
1,
πx
,
x ̸= 0, x = 0.
(3.1)
246
K. Parand et al. / Journal of Computational and Applied Mathematics 237 (2013) 244–252
For each integer k and the mesh size h, the Sinc basis functions are defined on R by Sk (h, x) ≡ sinc
x − kh
π sin( h (x − kh)) , π = (x − kh) h 1,
h
x ̸= kh,
(3.2)
x = kh.
Approximations can be constructed for infinite, semi-infinite and finite intervals. Define the function w = φ(x) = ln (sinh(x)), which is a conformal mapping from DE , the eye-shaped spatial domain in the z-plane, onto the infinite strip, DS in the w -plane, where
DE = z ∈ R : |arg(sinh(z ))| < d ≤
π 2
DS = w ∈ R : w = x + iy, |y| < d ≤
, π 2
.
The basis functions on (0, +∞) are taken to be the composite translated Sinc functions Sk (x) ≡ S (k, h) ◦ φ(x) = sinc
φ(x) − kh h
,
(3.3)
where S (k, h) ◦ φ(x) is defined by S (k, h)(φ(x)). The inverse map of ω = φ(x) is x = φ −1 (w) = ln(ew +
1 + e2w ).
(3.4)
Thus we may define the inverse images of the real line and of the evenly spaced nodes {kh}
k=+∞ k=−∞
as
Γ = {φ −1 (w) : w ∈ DS } = (0, +∞),
(3.5)
and xk = φ −1 (kh) = ln(ekh +
1 + e2kh ),
k = 0, ±1, ±2, . . . .
(3.6)
Let w(x) denotes a non-negative, integrable, real-valued function over the interval Γ , we define L2w (Γ ) = {v : Γ → R | v is measurable and ∥v∥w < ∞},
(3.7)
where ∞
|v(x)| w(x)dx 2
∥v∥w =
21
,
(3.8)
0
is the norm induced by the inner product of the space L2w (Γ ),
⟨u, v⟩w =
∞
u(x)v(x)w(x)dx.
(3.9)
0
Thus {Sk (x)}k∈Z with constant h denote a system which is mutually orthogonal under Eq. (3.9), i.e.,
⟨Skn (x), Skm (x)⟩w(x) = hδnm ,
(3.10)
where w(x) = coth(x) and δnm is the Kronecker delta function. This system is complete in Lw (Γ ), and for any function f ∈ L2w (Γ ) the following expansion holds 2
f (x) ∼ =
N
fk Sk (x),
(3.11)
k=−N
with fk =
⟨f (x), Sk (x)⟩w(x) . ∥Sk (x)∥2w(x)
(3.12)
The fk = f (xk ) are the discrete expansion coefficients associated with the family {Sk (x)} and the mesh size is given by
h=
2π d
αN
,
where N is suitably chosen, α and d depend [40] on the asymptotic behavior of f (x).
K. Parand et al. / Journal of Computational and Applied Mathematics 237 (2013) 244–252
247
Definition. Let H 2 (DE ) be the class of functions f which are analytic in DE , satisfy
|f (z )|dz → 0,
φ −1 (x+P )
x → ±∞,
where P = {iy : |y| < d ≤ π2 }, and on the boundary of DE (denoted ∂ D ) satisfy N (f ) =
∂D
|f (z )dz | < ∞.
Interpolation for the function in H 2 (DE ) is defined in the following theorem. Theorem ([41]). Assume that f φ ′ ∈ H 2 (DE ) then for all x ∈ Γ
∞ N (f φ ′ ) −π d/h N (f φ ′ ) E (f , h)(x) = f (x) − ≤2 e . f (kh)S (k, h) ◦ φ(x) ≤ 2π d sinh(π d/h) πd k=−∞ Moreover, if |f (x)| ≤ Ce−α|φ(x)| , x ∈ Γ , for some positive constants C and α , and if the selection h =
πd αN
≤
2π d , ln(2)
then
N √ √ f (kh)S (k, h) ◦ φ(x) ≤ C2 Ne(− π dα N ) , f (x) − k=−N where C2 depends only on f , d and α . 3.1. The matrix representation of the derivatives of Sinc basis functions at nodal points The nth derivative of function f at the nodal point at xk is denoted by dr 1 (r ) δk,j = [S (k, h) ◦ φ(x)]|x=xj . r h dφ r
(3.13) (r )
The expressions in Eq. (3.13) for each k and j can be stored in a matrix I (r ) = [δk,j ] where (0)
I (0) = [δk,j ],
δk(0,j) = [S (k, h) ◦ φ(x)]|x=xj =
(1)
δk(1,j)
(2)
δk(2,j)
I (1) = [δk,j ],
I (2) = [δk,j ],
1, 0,
k = j, k ̸= j,
k = j, 0, = h [S (k, h) ◦ φ(x)]|x=xj = (−1)j−k , k ̸= j, dφ j−k −π 2 2 , k = j, d = h2 2 [S (k, h) ◦ φ(x)]|x=xj = −32(−1)j−k dφ , k ̸= j. (j − k)2 d
(3.14)
(3.15)
(3.16)
The size of the matrix is determined by the limits on row and column indices k, j. 4. Solving the Thomas–Fermi equation For the boundary conditions in Eq. (2.4), the Sinc basis functions in Eq. (3.3) do not have a derivative and when the value of x tends to zero, thus we modify the Sinc basis functions as x x2 + 1
Sk (x).
(4.1)
In order to discretize Eq. (2.3) by using Sinc-collocation, first of all, we approximate y(x) as yN (x) = uN (x) + p(x),
(4.2)
where uN (x) =
N k=−N
ck
xSk (x) x2 + 1
.
(4.3)
248
K. Parand et al. / Journal of Computational and Applied Mathematics 237 (2013) 244–252
In order to approximate the solution of Eq. (2.3) with boundary conditions, we construct a polynomial p(x) that satisfies Eq. (2.4), this polynomial is given by p(x) =
λ x+λ
.
(4.4)
In Eq. (4.4), λ is constant to be determined. It is noted that the approximate solution yN (x), satisfies the boundary conditions in Eq. (2.4), since lim p(x) = 1,
lim uN (x) = 0,
x→0
x→0
lim p(x) = 0,
(4.5)
lim uN (x) = 0.
x→∞
x→∞
We refer the interested reader to [42,43] which address these issues. Now we display equations for yN (x), y′N (x) and y′′N (x). Using Eqs. (3.14)–(3.16) and (4.3) with substitution the result at the Sinc points:
xj = ln ejh +
1 + ln(e2jh ) ,
j = −N − 1 . . . N ,
(4.6)
we get yN (xj ) =
λ xj + λ
+
cj xj x2j
+1
,
(4.7)
N λ ck + yN (xj ) = − (xj + λ)2 k=−N
1
′
1 + x2j
−
2x2j
(1 + x2j )2
δk(0,j)
+
xj φ ′ (xj )
1 + x2j
δk(1,j)
,
(4.8)
and
8x3j −6xj ck + yN (xj ) = + δ (0) (xj + λ)3 k=−N (1 + x2j )2 (1 + x2j )3 k,j 2 4x2j φ ′ (xj ) 2φ ′ (xj ) xj φ ′ (xj ) xj φ ′′ (xj ) (1) (2) + − + δk,j + δk,j . 1 + x2j (1 + x2j )2 1 + x2j 1 + x2j 2λ
′′
N
(4.9)
4.1. The matrix representation of the Sinc-collocation method for solving the equation Suppose that
8x3 −6x + , (1 + x2 )2 (1 + x2 )3 2 xφ ′ (x) 1 Y(x) = , Z(x) = √ , 1 + x2 x f (x )
W(x) =
−N −1
D (f) =
..
.
X(x) =
T(x) = 0
2φ ′ (x) 1 + x2 x
x2 + 1
−
..
.
0
(1 + x2 )2
+
xφ ′′ (x) 1 + x2
,
f (xN )
,
,
f (x0 )
4x2 φ ′ (x)
0
c−N
⃗ = . C . .
,
(x−N −1 + λ) .. . λ (xN + λ)
⃗ = Λ
cN
λ
.
(4.10)
Substituting Eqs. (4.7) and (4.9) in Eq. (2.3) and subject to Eq. (4.2) we obtain:
⃗+ D (W)I (0) + D (X)I (1) + D (Y)I (2) C
∂2 3 ⃗ + D (Z)(Λ ⃗ + D (T)I (0) C ⃗ ) 2 = 0. Λ 2 ∂x
(4.11)
N The 2N + 1 coefficients {ck }kk= =−N and the unknown λ determined by solving Eq. (4.11) and evaluating the result at the Sinc points:
xj = ln ejh +
1 + ln(e2jh ) ,
j = −N − 1 . . . N .
(4.12)
K. Parand et al. / Journal of Computational and Applied Mathematics 237 (2013) 244–252
249
When −N − 1 ≤ j ≤ N , −N ≤ k ≤ N we remove the last column, the first and last rows of I (0) , I (1) and I (2) to arrive at the non-square matrices:
... .. . ...
−1
0
I (0) = 1 0
I (1)
=
0
0 , 1 1 2
..
.
..
.
1
..
.
..
.
..
.
..
.
...
..
.
1 2
(−1)2N +1 2N
I
(2)
...
0
−
(4.13)
−
2
π2 − 3 2 = 2 − 2 2 .. . −2(−1)2N +1 (2N + 1)
2
2 22
(−1)2N +1 2N + 1 .. . 1 , 2 −1
(4.14)
0
...
..
.
..
.
..
.
..
.
..
.
..
.
..
.
..
.
..
.
..
.
−2(−1)2N +1 (2N + 2)1 .. . 2 − 2 2 . 2 2 2 π −
(4.15)
3
Eq. (4.11) gives 2N + 2 nonlinear algebraic equations which can be solved for the unknown coefficients ck and λ by using the well known Newton’s method by Maple programming and we use
ck = 0, λ = 0,
k = −N . . . N ,
as starting points to obtain convergence of the method, consequently, y(x) given in Eq. (2.3) can be calculated. One measure of the rapidity of convergence of the procedure is provided by calculation of the value of the initial slope y′ (0) of the Thomas–Fermi potential [44]. 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),
(4.16)
where Z is the nuclear charge. Kobayashi et al. [45] provided a highly accurate numerical solution of the Thomas–Fermi as y′ (0) = −1.588071. Table 1 shows the comparison of y′ (0), between Padé approximation used by [46,47] and the y′N (0) = − λ1 obtained by Sinc-collocation method and choosing α = 1, d = (h =
2π 3
, N = 32 and the mesh size is given by El-Gamel et al. [40]
2π d ). αN
The approximations of y(x) obtained by Sinc-collocation method and the methods introduced in [46–48] are listed in Table 2 by choosing α = 1, d = 23π , N = 32 and h = √2π . 3N
250
K. Parand et al. / Journal of Computational and Applied Mathematics 237 (2013) 244–252 Table 1 Results for the value of y′ (0) obtained by Sinc-collocation method for α = 1, d = 23π , N = 32 and h = √2π . 3N
N
Sinc-collocation method
h = √2π
Padé
Liao [46]
Khan [47]
4 8 16 32
−1.552409174 −1.581426241 −1.587396189 −1.588070339
1.813799365 1.282549831 0.906899682 0.641274915
[10, 10] [20, 20] [30, 30] −
−1.51508 −1.58606 −1.58281 −
−1.573824678 −1.586494973 −1.582901907 −
3N
Table 2 Approximations of y(x) for Sinc-collocation method, by choosing α = 1, d = 23π , N = 32 and h = √2π . 3N
x
Sinc-collocation method
Liao [46]
Khan [47]
Bender [48]
1.00 2.00 4.00 6.00 8.00 10.00 50.00
0.4240642728 0.2430282344 0.1084800242 0.0594010318 0.0365353166 0.0243069553 0.0006213281
0.424008000 0.243009000 0.108404000 0.059423000 0.036587300 0.024314300 0.000632255
0.423772000 0.242718000 0.109632000 0.063816200 0.043285900 0.032208100 0.004730890
0.4240 0.2430 0.1084 0.0594 0.0365 0.0243 0.0006
Fig. 1. Thomas–Fermi graph obtained by Sinc-collocation method (solid line) in comparison with numerical solution (solid circles) [48] for N = 32.
Fig. 1 shows the result graph of Thomas–Fermi for N = 32 and it compares the new results with numerical results in [48]. Fig. 2 shows the logarithmic scale for the relative errors versus N. We refer the interested reader to [49,50] for more applications of spectral techniques. 5. Summary and conclusions The Sinc-collocation method is used to solve the Thomas–Fermi equation which is defined in the semi-infinite interval and has singularity at x = 0 and its boundary condition occurs at infinity. The numerical results demonstrate the reliability and efficiency of using the new method to solve such problems. The results of the present method for this type of problem clearly indicate that the new method is accurate even when singularity occurs at the boundary. Acknowledgments The authors would like to thank the referees for their comments and suggestions that have improved the paper. Also they would like to thank professor Taketomo Mitsui, the principal editor and professor Masaaki Sugihara, the advisory editor of
K. Parand et al. / Journal of Computational and Applied Mathematics 237 (2013) 244–252
251
Fig. 2. Logarithmic scale of the relative errors versus N for the value of y′ (0).
the journal for their useful comments and for managing the review process for this paper during 8 November 2009–1 August 2012. The first author (KP) would like to thank Shahid Beheshti University for the awarded grant. References [1] M. Dehghan, F. Fakhar-Izadi, The spectral collocation method with three different bases for solving a nonlinear partial differential equation arising in modeling of nonlinear waves, Math. Comput. Modelling 53 (2011) 1865–1877. [2] F. Fakhar-Izadi, M. Dehghan, The spectral methods for parabolic Volterra integro–differential equations, J. Comput. Appl. Math. 235 (2011) 4032–4046. [3] Y. Maday, B. Pernaud-Thomas, H. Vandeven, Reappraisal of Laguerre type spectral methods, La Recherche Aerospatiale 6 (1985) 13–35. [4] D. Funaro, Computational aspects of pseudospectral Laguerre approximations, Appl. Numer. Math. 6 (1990) 447–457. [5] B.Y. Guo, J. Shen, Laguerre-Galerkin method for nonlinear partial differential equations on a semi-infinite interval, Numer. Math. 86 (2000) 635–654. [6] J. Shen, Stable and efficient spectral methods in unbounded domains using Laguerre functions, SIAM J. Numer. Anal. 38 (2000) 1113–1133. [7] B.Y. Guo, Gegenbauer approximation and its applications to differential equations on the whole line, J. Math. Anal. Appl. 226 (1998) 180–206. [8] B.Y. Guo, Jacobi spectral approximation and its applications to differential equations on the half line, J. Comput. Math. 18 (2000) 95–112. [9] B.Y. Guo, Jacobi approximations in certain Hilbert spaces and their applications to singular differential equations, J. Math. Anal. Appl. 243 (2000) 373–408. [10] C.I. Christov, A complete orthogonal system of functions in L2 (−∞, ∞) space, SIAM J. Appl. Math. 42 (1982) 1337–1344. [11] J.P. Boyd, Spectral methods using rational basis functions on an infinite interval, J. Comput. Phys. 69 (1987) 112–142. [12] J.P. Boyd, Orthogonal rational functions on a semi–infinite interval, J. Comput. Phys. 70 (1987) 63–68. [13] M.R. Eslahchi, M. Dehghan, S. Amani, The third and fourth kinds of Chebyshev polynomials and best uniform approximation, Math. Comput. Modelling 55 (2012) 1746–1762. [14] 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 (2000) 117–147. [15] K. Parand, M. Razzaghi, Rational Chebyshev tau method for solving Volterra’s population model, Appl. Math. Comput. 149 (2004) 893–900. [16] K. Parand, M. Razzaghi, Rational Chebyshev tau method for solving higher-order ordinary differential equations, Int. J. Comput. Math. 81 (2004) 73–80. [17] K. Parand, M. Razzaghi, Rational Legendre approximation for solving some physical problems on semi-infinite intervals, Phys. Scr. 69 (2004) 353–357. [18] K. Parand, M. Shahini, Rational Chebyshev pseudospectral approach for solving Thomas–Fermi equation, Phys. Lett. A 373 (2009) 210–213. [19] K. Parand, A. Taghavi, Rational scaled generalized Laguerre function collocation method for solving the Blasius equation, J. Comput. Applied. Math. 233 (2009) 980–989. [20] A. Saadatmandi, M. Dehghan, A tau approach for solution of the space fractional diffusion equation, Comput. Math. Applic. 62 (2011) 1135–1142. [21] F. Shakeri, M. Dehghan, A hybrid Legendre tau method for the solution of a class of nonlinear wave equations with nonlinear dissipative terms, Numer. Methods Partial Differential Equations 27 (2011) 1055–1071. [22] J.P. Boyd, Chebyshev and Fourier Spectral Methods, second ed., Dover, New York, 2000. [23] T. Okayama, T. Matsuo, M. Sugihara, Sinc–collocation methods for weakly singular Fredholm integral equation of the second kind, J. Comput. Appl. Math. 234 (2010) 1211–1227. [24] A. Saadatmandi, M. Dehghan, M.R. Azizi, The Sinc-Legendre collocation method for a class of fractional convection-diffusion equations with variable coefficients, Commun. Nonlinear Sci. Numer. Simul. 17 (2012) 4125–4136. [25] A. Saadatmandi, M. Dehghan, The use of Sinc–collocation method for solving multi–point boundary value problems, Commun. Nonlinear Sci. Numer. Simul. 17 (2012) 593–601. [26] S. Timothy Carlson, Jack Dockery, John Lund, A Sinc–collocation method for initial value problems, Math. Comp. 217 (1997) 215–235. [27] K. Parand, M. Dehghan, A. Pirkhedri, Sinc-collocation methods for solving the Blasius equation, Phys. Letl. A 373 (2009) 4060–4065. [28] V.G. Koures, F.E. Harris, Sinc collocation in quantum chemistry: Solving the planar coulomb Schrödinger equation, International Journal of Quantum Chemistry: Quantum Chemistry Symposium 30 (1996) 1311–1318. [29] Kamel Al-Khaled, Theory and computation in singular boundary value problems, Chaos Solitons Fractals 33 (2007) 678–684. [30] A. Saadatmandi, M. Razzaghi, M. Dehghan, Sinc–Galerkin solution for nonlinear two-point boundary value problems with application to chemical reactor theory, Math. Comput. Modelling 42 (2005) 1237–1244. [31] F. Stenger, Sinc–Galerkin method of solution of boundary value problems, Math. Comp. 33 (1979) 85–109.
252
K. Parand et al. / Journal of Computational and Applied Mathematics 237 (2013) 244–252
[32] M. Tatari, M. Dehghan, An efficient method for solving multi–point boundary value problems and applications in physics, Journal of Vibration and Control 18 (2012) 1116–1124. [33] A. Cedillo, A perturbative approach to the Thomas–Fermi equation in terms of the density, J. Math. Phys. 34 (1993) 2713–2717. [34] H.T. Davis, Introduction to Nonlinear Differential and Integral Equations, Dover, New York, 1962. [35] S. Chandrasekhar, Introduction to the Study of Stellar Structure, Dover, New York, 1967. [36] M. Dehghan, R. Jazlanian, A fourth–order central Runge–Kutta scheme for hyperbolic conservation laws, Numer. Methods Partial Differential Equations 26 (2010) 1675–1692. [37] C.M. Bender, K.A. Milton, S.S. Pinsky, L.M. Simmons Jr, A new perturbative approach to nonlinear problems, J. Math. Phys. 30 (1989) 1447–1455. [38] H.A. Bethe, R.W. Jackiw, Intermediate Quantum Mechanics, W.A. Benjamin Inc, New York, 1968. [39] V.B. Mandelzweig, F. Tabakin, Quasilinearization approach to nonlinear problems in physics with application to nonlinear ODEs, Comput. Phys. Commun. 141 (2001) 268–281. [40] M. El-Gamel, S.H. Behiry, H. Hashish, Numerical method for the solution of special nonlinear fourth-order boundary value problems, Appl. Math. Comput 145 (2003) 717–734. [41] F. Stenger, Integration formulas via the trapezoidal formula, Journal of Institute Mathematics and its Applications 12 (1973) 103–114. [42] J. Lund, K. Bowers, Sinc Methods for Quadrature and Differential Equations, SIAM, Philadelphia, 1992. [43] F. Stenger, Numerical Methods Based on Sinc and Analytic Functions, Springer-Verlag, New York, 1993. [44] B.J. Laurenzi, An analytic solution to the Thomas–Fermi equation, J. Math. Phys. 10 (1990) 2535–2537. [45] S. Kobayashi, T. Matsukuma, S. Nagi, K. Umeda, Accurate value of the initial slope of the ordinary T-F function, J. Phys. Soc. Jpn. 10 (1955) 759–762. [46] S.J. Liao, An explicit analytic solution to the Thomas–Fermi equation, Appl. Math. Comput. 144 (2003) 495–506. [47] H. Khan, H. Xu, Series solution to the Thomas–Fermi equation, Phys. Lett. A 365 (2007) 111–115. [48] C.M. Bender, S.A. Orszag, Advanced Mathematical Methods for Scientists and Engineers, McGraw-Hill, New York, 1978. [49] J.P. Boyd, The optimization of convergence for Chebyshev polynomial methods in an unbounded domain, J. Comput. Phys. 45 (1982) 43–79. [50] M. Shamsi, M. Dehghan, Determination of a control function in three-dimensional parabolic equations by Legendre pseudospectral method, Numer. Methods Partial Differential Equations 28 (2012) 74–93.