- Email: [email protected]
Rational Chebyshev series for the Thomas–Fermi function: Endpoint singularities and spectral methods
Journal of Computational and Applied Mathematics 244 (2013) 90–101
Contents lists available at SciVerse ScienceDirect
Journal of Computational and Applied Mathematics journal homepage: www.elsevier.com/locate/cam
Rational Chebyshev series for the Thomas–Fermi function: Endpoint singularities and spectral methods John P. Boyd Department of Atmospheric, Oceanic and Space Science, University of Michigan, 2455 Hayward Avenue, Ann Arbor, MI 48109, United States
article
abstract
info
Article history: Received 21 September 2011 Received in revised form 19 February 2012 Keywords: Pseudospectral Chebyshev polynomials Rational Cebyshev functions Thomas–Fermi Orthogonal rational functions
√
We solve the Thomas–Fermi problem for neutral atoms, uyy − (1/ y)u3/2 = 0 on y ∈ [0, ∞] with u(0) = 1 and u(∞) = 0, using rational Chebyshev functions TLn (y; L) to illustrate some themes in solving differential equations on a semi-infinite interval. L is a user-choosable numerical parameter. The Thomas–Fermi equation is singular at the origin, giving a TL convergence rate of only √fourth order, but this can be removed by the change √ of variables, z = y with v(z ) = u(y(z )). The function v(z ) decays as z → ∞ with a term in z −3 , which is consistent with a geometric rate of convergence. However, the asymptotic series has additional terms with irrational fractional powers beginning with z −4.544 . In spite of the faster spatial decay, the irrational powers degrade the convergence rate to slightly larger than tenth order. This vividly illustrates the subtle connection between the spatial decay of u(x) and the decay-with-degree of its rational Chebyshev series. The TL coefficients an (L) are hostages to a tug-of-war between a singularity on the negative real axis, which gives a geometric rate of convergence that slows with increasing L, and the slow inverse power decay for large z, which gives quasi-tenth order convergence with a proportionality constant that decreases inversely as a power of L. For L = 2, we can approximate uy (0) (= vzz (0)) to 1 part in a million with a truncation N of only 20. L = 64 and N = 600 gives uy (0) = −1.5880710226113753127186845, correct to 25 decimal places. © 2012 Elsevier B.V. All rights reserved.
1. Introduction The Thomas–Fermi problem is 1 uyy − √ u3/2 = 0, y
u(0) = 1,
u(∞) = 0,
y ∈ [0, ∞]
(1)
where a subscripted coordinate denotes differentiation with respect to that variable. Because of its importance to theoretical physics, computing its solutions has attracted the attention of the Nobel laureates John Slater (chemistry) [1] and Richard Feynman (physics) [2] and of course Enrico Fermi [3]. Vannevar Bush, who created an analog computer known as the Differential Analyser, thus solved the Thomas–Fermi equation by computer more than a decade before the first electronic computer was built [4]; later Bush was effectively the first Presidential Science Advisor and the father of the National Science Foundation. Computing pioneer Nicholas Metropolis and physicist Edward Teller, a major figure in the development of the hydrogen bomb, also attacked this problem [2]. Hans Bethe [5] and Julian Schwinger [6,7] are other Nobel laureates who have enlarged Thomas–Fermi theory. This attention is in some ways remarkable because the neutral atom Thomas–Fermi problem is parameter-free: the Thomas–Fermi function u(y) is completely described by a single graph or table. Arnold Sommerfeld, the central figure in the
E-mail addresses: [email protected], [email protected]. 0377-0427/$ – see front matter © 2012 Elsevier B.V. All rights reserved. doi:10.1016/j.cam.2012.11.015
J.P. Boyd / Journal of Computational and Applied Mathematics 244 (2013) 90–101
91
‘‘old quantum’’ theory between Bohr’s 1913 hydrogen atom model and the dawn of quantum mechanics in 1926, gave the first analytical approximation [8,9]: 1 u(y) = 3.886 . 1 + (y/5.2415)0.772
(2)
John Mason [10] gave the improved approximation u(y) ≈
√
1 + 1.81061 y + 0.60112y
√
2
1 + 1.81061 y + 1.39515y + 0.77112y3/2 + 0.21465y2 + 0.04793y5/2
(3)
which has an absolute error for all y ∈ [0, ∞] which is smaller than 0.000012. Allan McLeod [11] has computed the Chebyshev coefficients of two series, one for y < 4 and the other for larger y, that yield the Thomas–Fermi function to about ten decimal place accuracy using only 17 terms of either series. However, precisely because of its simplicity, the Thomas–Fermi equation has been used as a testbed for a wide variety of semi-analytical and numerical schemes as exemplified by the delta-perturbation theory [12–15], Adomian decomposition [16–19], variational methods [20–23], homotopy perturbation analysis [24–29], Padé approximants [30,31] and even a technique from seventy-year old notes of Ettore Majorana [32]. We shall not review these efforts in detail because much of these merit the harsh label of Oscar Zariski as ‘‘spiccia matematica’’, that is, ‘‘small change mathematics’’. A modest portion of this work will be retroactively justified when techniques first tested on the Thomas–Fermi problem are successfully applied at the research frontier. Unfortunately, most semi-analytical methods applied to the Thomas–Fermi problem fail for anything more elaborate than a simple nonlinear ordinary differential equation. Fernandez [31,33,34] and Lipscombe [35] offer detailed criticism. The author of [22] is skewered for inflating the impact factor of his journal and of his personal citation record in [36,37]. Our reason for adding to the already vast literature on the Thomas–Fermi equation is that this problem is an illuminating example for Chebyshev spectral methods because u(y) is singular at both endpoints. For functions analytic on the entire interval, spectral methods typically converge exponentially fast with the number N of numerical degrees of freedom. However, endpoint singularities wreck this ‘‘spectral accuracy’’. McLeod [11] successfully finessed both difficulties, but did not discuss the singularities or the full range of options for coping with them. Parand and Shahini [38] obtained good accuracy using rational Chebyshev functions, as here, but were sketchy in their treatment of the singularity at the origin, completely ignored the difficulty as y → ∞, and only used a maximum of N = 12. In this article, we shall systematically treat both difficulties and obtain a single series of rational Chebyshev functions TLn [39] that converges rapidly for all of y ∈ [0, ∞]. We further shall provide something lacking from previous discussions of the TL basis: how decay as an inverse power of y for large y and a singularity at the origin jointly control the rate of convergence of the spectral series. 2. Rational Chebyshev functions for the semi-infinite interval The merits and mechanics of pseudospectral methods in general are discussed in [40–42]. As the author has previously explained in [39,43], rational Chebyshev functions are a good spectral basis for the semi-infinite interval because their connection with ordinary Fourier series yields simple theory, an exponential rate of convergence for functions analytic on the interval with sufficiently fast decay, and simplicity of programming. These basis functions are defined by TLn (y; L) ≡ Tn (x) = cos(nt )
(4)
where the argument of the rational Chebyshev functions, y, is related to the argument of the Chebyshev polynomials, x, and the argument of the trigonometric functions, t, via y=
L(1 + x) 1−x
y = L cot
2
↔x=
y−L
t
2
(5)
y+L
↔ t = 2 arccot
y L
.
(6)
The first few functions are given explicitly in Table 1. The pseudospectral interpolation points are yi = L cot
2
ti
2
↔ ti ≡
(2i − 1) π 2N + 2
i = 1, . . . , N + 1.
(7)
The Thomas–Fermi solution is represented as u(y) ≈
N −1 n=0
an TLn (y; L)
(8)
92
J.P. Boyd / Journal of Computational and Applied Mathematics 244 (2013) 90–101 Table 1 Rational Chebyshev functions for the semi-infinite interval: TLn (y). n
TLn (y; L)
0 1 2 3 4
(y − L)/(y + L) (y2 − 6yL + L2 )/(y + L)2 (y − L)(y2 − 14yL + L2 )/(y + L)3 (y4 − 28Ly3 + 70L2 y2 − 28L3 y + L4 )/(y + L)4
1
where the pseudospectral procedure for computing the an is described in [39,43] and the Maple code in the Appendix. Derivatives of the basis functions are computed applying the chain rule to derivatives of cos((n − 1)t ) to convert t-derivatives to y-derivatives. The constant L is a user-choosable ‘‘map parameter’’. Although there are sophisticated ways to estimate the best choice of L [39], in practice it is usual to begin with L equal to the expected dominant length scale of u(y), and then experiment. The criterion for ‘‘optimum’’ is rate of convergence, which can be evaluated from a plot of the computed coefficients. Note that both theory and the experiments below show that (i) the optimum L increases significantly with increasing truncation N of the TL series and (ii) one can obtain very accurate approximations for a given N with L anywhere within a rather large range. The following pair of theorems have not been previously published. Collectively, they describe how behavior at the endpoints, y = 0 and y → ∞, controls the rate of TL series convergence. Theorem 1 (Singularities at the Endpoint y = 0). Singularities at the endpoint y = 0 have precisely the same effect on convergence of the series of rational Chebyshev functions TL as they have the convergence of the Chebyshev polynomial series on a finite interval x ∈ [0, b]. A list of various species of singularities and their rates of convergence is given at the end of Chapter 2 of [42] and also [44]. In particular, if u(y) has a branch point at the origin with the singular term proportional to yψ , then its rational Chebyshev coefficients will satisfy an ∼ constant /n2ψ+1
(9)
as follows from [45]. Proof. When y is sufficiently small, the rational Chebyshev to Chebyshev polynomial change-of-coordinate, y = L(1 + x)/(1 − x), simplifies through power series expansion to y ≈ (1 + x)
L 2
+ O([1 + x]2 ).
(10)
Thus, the relationship between the two coordinates is linear near the origin and is therefore no map at all except for an irrelevant rescaling through the map parameter L. In other words, TLn (y; L) ∼ Tn
2 −1 + y , L
y≪1
(11)
for small y. (Fig. 1 illustrates this asymptotic relation.) It then follows that the effects of singularities on rational Chebyshev and Chebyshev polynomials must be the same since for small y, they are the same. Theorem 2 (Convergence Effects of Inverse Power Decay). 1. If a function u(y), not restricted to the Thomas–Fermi function, decays asymptotically as O(1/yk ) where k is an integer and is free of singularities everywhere on y ∈ [0, ∞] except at infinity and also that the asymptotic expansion for large y contains only non-positive integral powers, then the rational Chebyshev series will converge exponentially fast. 2. If u(y) has an asymptotic expansion with y−k as the power of smallest non-integral k, then the rational Chebyshev coefficients will decay as an ∼
constant n2k+1
.
(12)
Proof. The point y = ∞ is mapped into t = 0. In this neighborhood, the mapping simplifies to y = L cot2
t
2
≈L
4 t2
,
t ≪ 1.
(13)
This implies that for any positive k y−k ∼ t 2k .
(14)
J.P. Boyd / Journal of Computational and Applied Mathematics 244 (2013) 90–101
93
TL & TN, n=100 1
0.5
0
-0.5
-1 0.05
0
0.1 y
0.15
0.2
Fig. 1. Comparison of TL100 (y; L) (solid) with T100 (−1 + (2/L)y) (dashed) and the difference between them (dotted) for L = 3. The two functions closely approximate each other only for rather small y, but this is sufficient to give both bases identical responses to singularities at y = 0.
Thus, a function whose leading asymptotic term in a non-positive integer k will spawn a zero of order 2k at t = 0. The function t 2k is symmetric about t = 0 and thus poses no problems for approximating u(y(t )) as a Fourier cosine series whose terms are all individually symmetric, too. This proves the first proposition, which is that the TL series will converge exponentially if k is an integer, and there are no singularities elsewhere on the real axis. If k is not an integer, however, then u(y(t )) will have a branch point at t = 0. Elliott’s analysis of the consequences of a branch point of exponent 2k in t proves the second proposition [45]. One might suppose that half-integral powers of k would also yield exponential convergence because the endpoint behavior of u(y(t )) is an integral power of t and therefore not a branch point. However, odd powers of t are antisymmetric with respect to t = 0, i.e., t 2m+1 = −(−t )2m+1 , and this is incompatible with a cosine basis. As illustrated for the related case on an infinite interval [46], exponential convergence can be retrieved by converting to the trigonometric basis and then applying a sine basis. Fig. 2 confirms this analysis: integral powers of y for large y allow an exponential rate of convergence whereas the convergence rate is algebraic for fractional powers of y. The example u(y) = 1/(1 + y)k is atypical in that it is not singular at infinity so that the rate of convergence for integer k is geometric, that is, proportional to exp(−qn) for some constant q. It is more common in applications for functions to be singular at infinity, or to have an infinite number of singularities with a limit at infinity, and then the convergence is ‘‘subgeometric’’ in the terminology of [42], that is, proportional to exp(−qnr ) for some constant r < 1, typically r = 1/2 or r = 2/3. We shall not elaborate further except to refer to Chapter 17 of [42,46] because the Thomas–Fermi problem has relatively strong singularities at infinity that make subgeometric convergence irrelevant. 3. Branch point at the origin 3.1. Transforming the problem The Thomas–Fermi neutral atom problem is 1 uyy = √ u3/2 = 0, y
u(0) = 1,
u(∞) = 0,
y ∈ [0, ∞].
(15)
Recalling that u(y) cannot vanish at the origin because of the boundary condition u(0) = 1, it is obvious that the second √ derivative of u(y) must be unbounded in the limit y → 0 because the second derivative is proportional to 1/ y. In 1930, a graduate student at the University of Michigan, Edward Baker, gave a complete description of the singularity by deriving the small-y series [47]: u = 1 + sy +
4 3
y3/2 +
2 5
s y5/2 + s2
3 70
y7/2 +
1 3
y3 + · · ·
where s is the value of the first derivative at the origin, shown numerically to be roughly s ≈ −1.58.
(16)
94
J.P. Boyd / Journal of Computational and Applied Mathematics 244 (2013) 90–101
SEMILOGY, integer k 100
100
LOGLOG, k=fraction
k=1 k=2 k=3
k=3/2 k=5/2 k=7/2
10-5
10-5
10-10
10-10
10-15
0
10
10-15 100
20
102 degree n
degree n
Fig. 2. The TL coefficients of the function (1 + y)−k for integer k (left graph) and half-integer k (right). With a logarithmic vertical scale and linear horizontal scale, geometric convergence appears linear; the left plot therefore shows the exponential convergence when k is an integer. On a log–log plot, decay proportional to a power of k is linear. The dotted guidelines show that the coefficients for y3/2 (top curve), y5/2 (middle curve) and y7/2 (bottom curve) asymptote to n−4 , n−6 and n−8 , respectively.
TL coefficients of the Thomas-Fermi function u(y) 100 coefficients, L=2 1.8/n**4
10-5
10-10 100
101 degree n
102
Fig. 3. Coefficients of the TL rational Chebyshev expansion for u(y). The dashed line is a guideline with slope proportional to n−4 .
Theorem 1 shows that if a function u(y) has a branch point of order k at the origin, then the TL coefficients will decay proportionally to n−(2k+1) . The leading singular term in Baker’s series is y3/2 , so the predicted rate of convergence is fourth order, that is, |an | ∼ constant/n4 . The TL series for the Thomas–Fermi function (8) does indeed show the expected fourth order rate of convergence as illustrated in Fig. 3. As explained in [42], a spectral series with coefficients decreasing as n−m will yield a maximum pointwise error decreasing one order slower as N −(m−1) where N is the truncation of the spectral series. Thus, the rational Chebyshev series for u(y) has only third order convergence. The first derivative of the n-th TL function at the origin is O(n2 ), so the rate of convergence to the correct value of s = uy (0) is only first order: numerical experiments showed that with L = 4, the error in uy (0) was about 0.174/N where N is the truncation of the TL series. One strategy adopted by McLeod [11] is to apply a spectral method in the coordinate z≡
√
y.
(17)
We prefer the approach of Fernandez [31], who also adopts the new unknown
v(z ) ≡
u(z 2 ).
(18)
The problem becomes z vvzz + (vz )2 − vvz − 2z 2 v 3 = 0,
This is the form we solved numerically.
v(0) = 1,
v(∞) = 0.
(19)
J.P. Boyd / Journal of Computational and Applied Mathematics 244 (2013) 90–101
95
TL coefficients of the Mason approx. to sqrt(u(z)) 100 coefficients, L=1 8 exp(−n) 10-5
10-10
10-15 0
10
20 degree n
30
40
Fig. 4. Solid: absolute values of the rational Chebyshev series for Mason’s quadratic-polynomial-divided-by-quintic polynomial approximation to v(z ); this approximation was used to initialize the Newton–Kantorovich iteration. The dashed line is a guideline that is the graph of exp(−n).
3.2. Power series The transformed differential equation is still singular at the origin; it is merely the solution that is desingularized. Because (19) is singular, the standard Maple command to solve a differential equation as a series returns only a blank. When a series of the form
v(z ) = 1 + a1 z + (1/2)sz 2 + a3 z 3 + a4 z 4 + · · ·
(20)
is substituted into (19), the residual is of the from z vvzz + (vz )2 − vvz − 2z 2 v 3 = −a1 + {3a3 + (3/2)a1 s − 2} z 2 + 8a4 + 8a1 a3 + s2 − 6a1 z 3 + · · · .
(21)
Because the constant in the residual is −a1 , it is impossible to satisfy the differential equation with a power series unless a1 = 0 ↔ vz (0) = 0.
(22)
The remaining coefficients can be found recursively in terms of s ≡ vzz (0); the specific choice s = −1.58807 . . . is the only choice that allows the sum of the power series to satisfy the boundary condition v(∞) = 0. Because v(z ) is not uniquely specified by the initial conditions v(0) = 1, vz (0) = 0, it is not possible to begin a standard initial value calculation at the origin itself. The best one can do is to march away from a point near the origin where the necessary initial values can be obtained from the power series. Remarkably, however, the TL pseudospectral method has no difficulties with the boundary value problem for v(z ). As is typical when the differential equation is singular at an endpoint, but the solution is not, the spectral method automatically picks out the ‘‘nice’’ solution without special tricks, or the imposition of additional, nonstandard boundary conditions. 4. Asymptotic series for large y and its consequences Arnold Sommerfeld, who supervised more Nobel laureates than any other physicist, derived his analytic approximation by investigating the asymptotic behavior of u(y) for large y. Charles Coulson, professor of theoretical chemistry at Oxford, and a student elaborated this into an asymptotic series of the form [48] u(y) ∼
144 y3
F F2 F3 1 − λ + 0.62569 2λ + 0.3138 3λ + · · · y y y
(23)
where F = 13.27094 [49] and
λ = (1/2)(−7 + −3
√
73) = 0.772001872.
(24)
The leading order y does not prevent even a geometric rate of convergence. This is illustrated in Fig. 4, which shows the geometric convergence of the TL series for Mason’s approximation to v(z ), which reproduces the leading order asymptotic term but omits all fractional powers. (By geometric decay, we mean convergence proportional to exp(−µn) for some n constant µ, like that of the terms in the geometric series n x .) However, Theorem 2 predicts that the next term of −3.772 O(y ) will, ignoring the effects of singularities at y = 0, slow the convergence of the TL coefficients to O(n−8.544 ).
96
J.P. Boyd / Journal of Computational and Applied Mathematics 244 (2013) 90–101
abs(TL coeffs)
TL coefficients for L=1 & L=64
10-10
10-20
L=1
L=64
10-30 100
200
300 degree n
400
500
Fig. 5. The absolute of the TL coefficients an of the Thomas–Fermi function for L = 1 and L = 64.
However, we actually solve Fernandez’ form. The asymptotic series for v(z ) is
v(z ) ∼ ∼ ∼
144
z6 12
1−
z3 12 z3
−
F
1−
z 2λ
F z 2λ
6F z 3+2λ
+ 0.62569
+ 0.62569
+ ···.
F2 z 4λ
F2 z 4λ
+ ···
+ ···
(25)
(26) (27)
The leading non-integral power is 1/z 3+2λ = 1/z 4.544 . This implies that the TL series for v(z ) has coefficients asymptotically proportional to an ∼
constant n10.088
.
(28)
Thus, the TL series for v(z ) has higher order convergence than the corresponding series for u(y). It is important to note that (28) is an asymptotic relation that holds only for sufficiently large n. Fig. 5 shows that the TL coefficients computed with two different values of L. On this plot with a logarithmic vertical scale and a linear horizontal scale, geometric decay appears linear. For both choices of L, the coefficients of small degree do in fact decay geometrically. For large n, however, the decay with degree slows dramatically. Fig. 6 shows the same two sets of coefficients, but with a logarithmic scale on both axes. Power law decay appears linear on a log–log plot. For both values of L, the coefficients asymptote to lines paralleling the theoretical prediction of n−10.088 which is shown as the dashed guideline. Fig. 7 is a different graph showing more precisely how the coefficients asymptote to a simple power law. Why does L affect the coefficients this way? And why is there a dramatic break in slope on the coefficient graphs? Darboux’ Principle, discussed in Section 2.10 of [50] and also [51], asserts that the asymptotic coefficients an of a power series are a sum of contributions, one from each of the singularities in the complex plane of the function represented by the power series. This principle extends to Fourier and all species of Chebyshev series also. Singularities on the expansion interval contribute terms that decay more slowly than geometrically—here, the fractional power decay of v(z ) is a branch point at z = ∞ that yields a contribution to an that decays slightly faster than tenth order. Singularities off the expansion interval contribute terms which decay exponentially at a rate which is controlled entirely by the location of the singularity; the exponential in degree n is multiplied by a slower-than-exponential factor, such as powers or logarithms of n or something more complicated, that depends only on the type of singularity (i.e., pole of some order or logarithm or other branch point). It is easy to compute the power series of v(z ) to high order. The coefficients strictly alternate in sign, suggesting that the convergence-limiting singularity is on the negative real axis. The usual ‘‘ratio test’’, an−1 /an = zs + O(1/n), suggests that the singularity is at z = zs ≈ −0.665. However, Hille [52] wrote (p. 158): ‘‘Formula (4.2) as well as Fowler’s result suggest that y(x) has a pole of order four . . . A direct substitution of a Laurent series . . . leads to unexpected complications and unmanageable computations . . . the singularity is more complicated than a pole’’. The same is true of the singularities of the Blasius function and could probably be analyzed by the methods of [53], but this is beyond the scope of this article.
J.P. Boyd / Journal of Computational and Applied Mathematics 244 (2013) 90–101
97
Log−log plot: TL coeffs. for L=1 & L=64
L=1
L=64
abs(TL coeffs)
10-10
0.0 01 /n 10
L=1
.08
10-20
8
L=64 10-30 100
101
102 degree n
Fig. 6. Same as previous figure, but using a logarithmic scale for both axes. The dashed line is a guideline with a slope of n−10.088 .
a(n)* n**10.088, normalized L=1
scaled coefficients
1.05
1
0.95 0
200
400
600
degree n Fig. 7. The solid curve is n10.088 an for L = 1, normalized by its value for n = 350. The coefficients of low degree do not follow the predicted asymptotic pattern because the asymptotic approximation only applies for sufficiently large n. The coefficients on the extreme right deviate from the predicted scaled magnitude because of aliasing errors, a general phenomenon with interpolatory spectral methods.
Fortunately, the location of the singularity controls exponential terms in the spectral coefficients [42,44]; the order and type of the singularity are weak (non-exponential) influences on the rate of convergence. If the singularity in z is sufficiently close to the origin, then z = L cot2 (t /2) ≈ (L/4)(t − π )2 + O([L/24])(t − π )4 .
(29)
If the TL-convergence-limiting singularity is on the negative real axis, then
ts = π ± i
4
|zs |. L This will generate a contribution to the asymptotic TL coefficients proportional to exp(−nℑ(ts )), that is, an ∼ A(n) exp(−n
4|zs |/L) + terms from other singularities
(30)
(31)
where A(n) is a factor that varies more slowly with n than the exponential, such as some combination of constants, powers of n and logarithms of n. This shows explicitly that increasing the map parameter L slows the convergence of the contributions of this term. It predicts that in the geometrically-converging regime, the L = 64 curve should fall proportional to about exp(−0.2), or in other words diminish by about 10−9 for each additional 100 terms in the series. In Fig. 5, this is close to what is observed for n < 300 where the coefficients transition to dominance by the tenth order contributions from the singularity at z = ∞.
98
J.P. Boyd / Journal of Computational and Applied Mathematics 244 (2013) 90–101 Table 2 Selected calculations of uy (0) (=vzz (0)); inaccurate digits claimed to be correct by their authors are in bold face. Entries up to 1955 from Kobayashi et al. Fermi (1928) [3] Baker (1930) [47] Bush and Caldwell (1931) [4] Miranda (1934) [57] Slater and Krutter (1935) [1] Feynman, Metropolis and Teller (1949) [2] Kobayashi et al. (1955) [49] Fernandez (2008) [31] Boyd (2011) [this article] (25 decimal places) ’’ [N = 20, L = 2] ’’ [N = 40, L = 4] ’’ [N = 100, L = 16] ’’ [N = 200, L = 32]
−1.58 −1.588558 −1.589 (analog computer [56]) −1.5880464 −1.58808 −1.58875 −1.5880710 −1.58807102261137531 −1.5880710226113753127186845 −1.58807 −1.5880710226 −1.588071022611375 −1.588071022611375312719
5. Iteration and elimination Because the differential equation is nonlinear, we employed a Newton–Kantorovich iteration [42], also known as ‘‘quasilinearization’’ [54], to reduce the problem to a sequence of linear differential equations: z vn δn,zz + (2z vn,z − v) δn,z + z vn,zz − vn,z − 6z 2 vn2 δn = z vn vn,zz + vn,z vn,z − vn vn,z − 2z 2 vn3
(32)
vn+1 = vn − δn
(33)
subject to the boundary conditions δn (0) = vn (0) − 1 and δn (∞) = 0. The differential equation was solved by truncating the TL series to N terms and then applying collocation, also known as the pseudospectral method, at (N − 2) points. The collocation conditions plus the two boundary conditions gives an N × N dense matrix problem which can be solved by Gaussian elimination, alias Crout reduction, for the vector of TL coefficients. We omit the details since these are discussed with many examples in [39,42]. The computational cost and memory can be greatly reduced by using a preconditioned iteration as explained in [42]. Since an N = 600 computation using direct matrix solves required only about 2.5 h on a 2009 laptop running Maple 13 in Windows Vista, we did not bother with such refinements. However, an inner iteration is very valuable for more difficult problems. The Newton–Kantorovich iteration requires an initialization or ‘‘first guess’’ v0 . We employed Mason’s rational approximation, and found that the residual (the right hand side of the equation) for δn was always reduced to a maximum pointwise value smaller than 10−31 after five iterations without any underrelaxation. The program is so simple we have included the entire Maple code in the Appendix. 6. Approximating the slope at the origin The first derivative of u(y) at left endpoint has always been of great interest because this one number is sufficient to convert the boundary value problem into a second order initial value problem with the initial conditions u(0) = 1 and uy (0) = s where s is a very accurate approximation to the first derivative at the origin. In principle, the initial value problem can be solved by a simple march from y = 0 to y → ∞ using one’s favorite Runge–Kutta or Adams–Bashforth scheme. In practice, such marching is ill-conditioned as reviewed most recently by Zaitsev, Matyushkin, and Shamonov [55]. Table 2 is only a partial list of the vast number of calculations of the Thomas–Fermi function, but the entries show that accuracy has increased with time but not monotonically; some very distinguished authors in fact obtained fewer correct digits than they thought. Kobayashi et al. and Fernandez deserve credit because all the digits they claimed are correct. The last four rows show the best approximation obtained from a given number N of terms of the TL basis. Note that the optimum L increases with N; this is quite typical in calculations with rational Chebyshev functions on both the infinite and semi-infinite intervals. Although we actually computed v(z ) and not u(y), it is easy to show that s ≡ uy (0) ≡ vzz (0).
(34)
Table 3 shows the errors in vzz (0) for various L and N. It is clear that small L is best for small N while large L is better for large N but awful for small N. Because the first derivative of TLn at z = 0 is O(n2 ) and the second derivative is O(n4 ), it is not surprising that the error in vzz (0) converges at a slower rate than the coefficients an . The table shows a seventh order rate of convergence instead of the expected sixth order. We do not have an explanation. As is usual when an exact solution is unknown, a solution of very high N, optimum L, and large value of the Maple precision parameter ‘‘Digits’’ was taken as the ‘‘exact’’ solution for computing errors. Our benchmark employed N = 650, L = 64, and Digits = 45 and required three and a quarter hours computation on our laptop. It was judged, by consistency with other
J.P. Boyd / Journal of Computational and Applied Mathematics 244 (2013) 90–101
99
Table 3 Computed values of uy (0) = vzz (0) for various L and N; the optimum values in a given row or column are in boldface. N
L=1
L=2
L=4
L=8
L = 16
L = 32
20 40 100 200 400 600
0.101e−3 0.303e−6 0.415e−9 0.293e−11 0.210e−13 0.118e−14
0.142e−5 0.343e−8 0.458e−11 0.319e−13 0.227e−15 0.127e−16
0.526e−5 0.390e−10 0.505e−13 0.349e−15 −0.245e−17 0.137e−18
0.961e−4 −0.237e−8 0.563e−15 −0.378e−17 −0.265e−19 0.147e−20
−0.526e−3 −0.226e−5 0.506e−16 −0.414e−19 −0.287e−21 −0.159e−22
−0.145e−3 −0.886e−11 −0.452e−21 −0.313e−23 −0.181e−24
0.0536
L = 64
−0.216 −0.720e−2 −0.342e−7 −0.627e−16 < 1.e − 25 < 1.e − 25
runs with different N , L and Digits, that at least the first twenty-five decimal places are correct, these being invariant among our highest resolution computations. 7. Summary The singularity of the Thomas–Fermi function at the origin, which would otherwise degrade convergence of the rational Chebyshev series to fourth order, can be eliminated by a simple transformation of the coordinate and unknown. The term z −4.544 in the asymptotic series for v(z ) slows the convergence to tenth order. The TL functions depend up a user-choosable constant, the ‘‘map parameter’’ L. The asymptotic TLn coefficient an as n → ∞ is the sum of contributions from the singularities of v(z ). Increasing L slows the rate of geometric convergence of the term in an due to the pole-and-branch-point on the negative real axis, but it greatly diminishes the proportionality constant in the tenth order contribution from the singularity at z = ∞. The result of this tension between two singularities which are affected oppositely by L is that a graph of the TL coefficients shows two regimes: geometrically fast (exp(−nµ)) decay for small n and slower, tenth order decay for large n. The ‘‘crossover point’’ between regimes grows rapidly with L. The result is that when the TL series is truncated to just twenty terms, the best choice is L = 2. When N = 600, however, a fast exponential decay is not important because N is so huge. It is better to pick L = 64 or larger to diminish the effect of the singularity at z = ∞. It is then possible to compute the second derivative of u(y) at the origin to twenty-five decimal places. The insensitivity to L, in the sense that accurate calculations can be made when L is many times larger or smaller than optimum, is a generic feature of rational Chebyshev series [46,39]. It is also common that the optimum L is not fixed, but rather grows with the truncation N of the series. The reason is the tension between resolving the behavior of v(z ) for large |z | and simultaneously resolving singularities off the expansion interval, z ∈ [0, ∞]. Increasing L places grid points at larger and larger z: good for resolving the singularity at infinity. However, a rational Chebyshev expansion can be always be transformed into an equivalent Fourier cosine series through the identity TLn (z ; L) ≡ cos(nt [z ; L]) where the trigonometric coordinate is a function of z and L only as z /L, which implies that increasing L effectively moves singularities at finite z closer to the real axis in t. Thus increasing L worsens the contributions of these singularities to the rational Chebyshev coefficients of high degree. Because the singularity at infinity usually generates a slower-than-geometric contribution to the asymptotic an – tenth order here, subgeometric proportional to exp(−qn2/3 ) more commonly – increasing N requires, for optimality, increasing L to better cope with the slower-decreasing contributions. In theory, the tenth order convergence can be replaced by asymptotic exponential convergence by using an exponential mapping. Our experiments here produced no improvement and are therefore omitted. As discussed in [58], mappings introduce additional variability that slows the convergence of the low degree coefficients while accelerating coefficients an for n > ncrossover for some ‘‘crossover’’ value ncrossover . Because the TL series for v is already converging at better than tenth order, the ‘‘crossover’’ point here is very large, and an exponential mapping does no good until the degree is in the thousands. Acknowledgments This work was supported by the National Science Foundation through grants OCE 0451951 and ATM 0723440. I thank Paolo Amore for his careful reading and comments. I also thank the three reviewers for their constructive suggestions. Appendix. Maple code to solve the Thomas–Fermi equation restart; L:= 2; # map parameter N:= 600; # number of basis functions & grid points; Digits:= 40; # number of decimal digits of precision itermax:= 6; # number of Newton iterations; ncol:= N−2; # number of interior collocations points; Mason25:= proc(y) local yh,vn,vd,Mas; # evaluates Mason’s approximation; yh:= sqrt(y); vn:= 1+1.81061*yh+0.60112*y;
100
J.P. Boyd / Journal of Computational and Applied Mathematics 244 (2013) 90–101
vd:= 1+1.81061*yh+1.39515*y+0.77211*y*yh+0.21465*y*y+ 0.04793*yh*y*y; vM:= (vn/vd); end proc: with(LinearAlgebra); fa:= Vector(N,orientation=column): v:= Vector(N,orientation=column): za:= Vector(N,orientation=column): a:= Vector(N,orientation=column): a0:= Vector(N,orientation=column); # TL differentiation matrices D0,D1,D2 for 0th, d/dz, d**2/dz**2; D0:= Matrix(ncol,N): D1:= Matrix(ncol,N): D2:= Matrix(ncol,N): # collocation points za; TL computed by trigonometric identity; for ii from 1 to ncol do ta[ii]:= evalf(Pi*(2*ii−1) / (2*ncol)); tt:= ta[ii]; za[ii]:= evalf(L*cot(tt/2)**2); for j from 1 to N do TL:= evalf(cos((j−1)*tt)); # Compute derivatives using trigonometric definition and chain rule; PT:= evalf(− (j−1)*sin((j−1)*tt)); PTT:= evalf(− (j−1)*(j−1)*TL); S := evalf(sin(tt/2)); C := evalf(cos(tt/2)); TLz:= ev alf(− S * S * S * PT/ (L * C)); TLzz:= evalf((S**5)*(2*C*S*PTT+(3-2*S*S)*PT)/(L*L * C*C*C* 2)); D0[ii,j] := TL ; D1[ii,j]:= TLz ; D2[ii,j]:= TLzz; od: od: # Expand initialization v0 as TL series with coeffs a0=TLEXPMat*Mason25(za); for ii from 1 by 1 to ncol do v0[ii]:= evalf(Mason25(za[ii]*za[ii])); od: TLEXPMat:= Matrix(ncol,ncol): for ii from 1 to ncol do for j from 1 to ncol do TLEXPMat[ii,j]:= evalf((2/ncol)*cos((ii−1)*ta[j])); od; od: for j from 1 to ncol do TLEXPMat[1,j]:= evalf(TLEXPMat[1,j]*0.5); od: for ii from 1 by 1 to ncol do a0[ii]:= 0; for j from 1 by 1 to ncol do a0[ii]:= evalf(a0[ii]+TLEXPMat[ii,j]*v0[j]); od: od: a0[N−1]:= 0; a0[N]:= 0; a:= a0: # Initialize Jacobian matrix and compute its last two rows; Jacobian:= Matrix(N,N):for j from 1 by 1 to N do Jacobian[N−1,j]:= evalf(cos((j−1)*Pi));Jacobian[N,j]:= 1; od: # Start of Newton–Kantorovich iteration; counter is ‘‘iter’’; for iter from 1 by 1 to itermax do for ii from 1 by 1 to ncol do v:= 0; vz:= 0; vzz:= 0; for j from 1 by 1 to N do v:= evalf(v + D0[ii,j]*a[j]); vz:= evalf(vz + D1[ii,j]*a[j]); vzz:= evalf(vzz + D2[ii,j]*a[j]); od: # fa is residual of TF ODE; fa[ii]:= evalf((v*vzz + vz*vz)*za[ii] − v*vz − 2*za[ii]*za[ii]*v*v*v); for j from 1 by 1 to N do Jacobian[ii,j]:= evalf(za[ii]*v * D2[ii,j] + (2*za[ii]*vz − v)*D1[ii,j] +(za[ii]*vzz − vz − 6*za[ii]*za[ii]*v*v)*D0[ii,j]); od: od: # Apply boundary conditions: v(0)=1, v(infinity)=0; vatzero:= 0; vatinfinity:= 0; for j from 1 by 1 to N dovatzero:= evalf(vatzero + cos((j−1)*Pi)*a[j]); vatinfinity:= evalf(vatinfinity + a[j]); od; fa[N−1]:= vatzero − 1; fa[N]:= vatinfinity; delta_a:= LinearSolve(⟨Jacobian|fa⟩); # compute TL coeffs. delta_aof δn (z ). for j from 1 by 1 to N do a[j]:= evalf(a[j] − delta_a[j]); od: residualLoo[iter]:= max(seq(abs(fa[k]),k=1..N)); print(iter,residualLoo[iter]); od: # end of Newton iteration loop; # Only plots & output below; anpoints := [seq([j,abs(a[j])],j=1..N]: plot2:= loglogplot(anpoints,style=line):guidepoints:= [seq([j,j**(−10.022)],j=1..N)]: plotguide:= loglogplot(guidepoints): plots[display](plot2,plotguide); vzzat0:= 0; for j from 1 to N do j1:= j−1; vzzat0:= evalf(vzzat0+ a[j]*4*cos(j1*Pi)*((1/3)*j1**4+(2/3)*j1**2)/L**2); od: vxzat0ezact:= −1.5880710226113753127186845; errvzzat0:= evalf(abs(vzzat0 vzzat0exact));lprint(N,vzzat0,errvzzat0); References [1] J.C. Slater, H.M. Krutter, The Thomas–Fermi method for metals, Phys. Rev. 47 (1935) 559–568. [2] R.P. Feynman, N. Metropolis, E. Teller, Equations of state of elements based on the generalized Fermi–Thomas theory, Phys. Rev. 75 (1949) 1561–1573.
J.P. Boyd / Journal of Computational and Applied Mathematics 244 (2013) 90–101
101
[3] E. Fermi, Eine statistische Methode zur Bestimmung einiger Eigenschaften des Atoms und ihre Anwendung auf die Theorie des periodischen Systems der Elemente, Z. Phys. 48 (1928) 73–79. [4] V. Bush, S.H. Caldwell, Thomas–Fermi equation solution by the differential analyzer, Phys. Rev. 38 (1931) 1898–1902. [5] H. Bethe, Thomas–Fermi theory of nuclei, Phys. Rev. 167 (1968) 879–907. [6] J. Schwinger, Thomas–Fermi model: the leading correction, Phys. Rev. A 22 (1980) 1827–1832. [7] J. Schwinger, Thomas–Fermi model: the second correction, Phys. Rev. A 24 (1981) 2353–2361. [8] A. Sommerfeld, Asymptotische Integration der Differentialgleichung des Thomas Fermischen Atoms, Z. Phys. 78 (1932) 283. [9] A. Sommerfeld, Integrazione asintotica dell’equazione differentiale di Thomas–Fermi, Rend. R. Accad. Lincei 15 (1932) 293–308. [10] J.C. Mason, Rational approximations to the ordinary Thomas–Fermi function and its derivatives, Math. Proc. Cambridge Philos. Soc. 84 (1964) 357–360. [11] A.J. McLeod, Chebyshev series solution of the Thomas–Fermi equation, Comput. Phys. Comm. 67 (1992) 389–391. [12] 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. [13] B. Laurenzi, An analytic solution to the Thomas–Fermi equation, J. Math. Phys. 31 (1990) 2535–2537. [14] A. Cedillo, A perturbative approach of the Thomas–Fermi equation in terms of density, J. Math. Phys. 34 (1993) 2713–2717. [15] D.A. Morales, An N-dimensional interpretation of the delta-expansion for the Thomas–Fermi equation, J. Math. Phys. 8 (1994) 3916–3921. [16] G. Adomian, Solution of the Thomas–Fermi equation, Appl. Math. Lett. 11 (1998) 131–133. [17] A.M. Wazwaz, The modified decomposition method and Padé approximants for solving the Thomas–Fermi equation, Appl. Math. Comput. 105 (1999) 11–19. [18] J.I. Ramos, Piecewise-adaptive decomposition methods, Chaos Solitons Fractals 40 (2009) 1623–1636. [19] H. Chu, Y. Liu, The new ADM-Pade technique for the generalized Emden–Fowler equations, Modern Phys. Lett. B 24 (2010) 1237–1254. [20] N. Anderson, A. Arthurs, Variational solution of the Thomas–Fermi equation, Quart. Appl. Math. 39 (1981) 127–129. [21] B. Burrows, P. Core, A variational-iterative approximation solution of the Thomas–Fermi equation, Quart. Appl. Math. 42 (1984) 73–76. [22] J.H. He, Variational approach to the Thomas–Fermi equation, Appl. Math. Comput. 143 (2003) 533–535. [23] M. Desaixj, D. Anderson, M. Lisak, Variational approach to the Thomas–Fermi equation, European J. Phys. 25 (2004) 699–705. [24] S. Liao, An explicit analytic solution to the Thomas–Fermi equation, Appl. Math. Comput. 144 (2003) 495–506. [25] H. Khan, H. Xu, Series solution to the Thomas–Fermi equation, Phys. Lett. A 365 (2007) 111–115. [26] A. El-Nahhas, Analytic approximations for Thomas–Fermi equation, Acta Phys. Pol. A 114 (2008) 913–918. [27] M.A. Noor, S.T. Mohyud-Din, Homotopy perturbation method for solving nonlinear higher-order boundary value problems, Int. J. Nonlinear Sci. Numer. Simul. 9 (2008) 395–408. [28] B. Yao, A series solution to the Thomas–Fermi equation, Appl. Math. Comput. 203 (2008) 396–401. [29] S.T. Mohyud-Din, A. Yildirim, Solving nonlinear boundary value problems using He’s polynomials and Pade approximants, Math. Probl. Eng. (2009) p. Art. No. 690547. [30] L.N. Epele, H. Fanchiotti, C.A.G. Canal, J.A. Ponciano, Padé approximant approach to the Thomas–Fermi problem, Phys. Rev. A 60 (1999) 280–283. [31] F.M. Fernandez, Comment on ‘Series solution to the Thomas–Fermi equation’, Phys. Lett. A 365 (2007) 111; Phys. Lett. A 372 (2008) 5258–5260. [32] S. Esposito, Majorana solution of the Thomas–Fermi equation, Amer. J. Phys. 70 (2002) 852–856. [33] F.M. Fernandez, On some approximate methods for nonlinear models, Appl. Math. Comput. 215 (2009) 168–174. [34] F.M. Fernandez, On two new applications of the homotopy analysis method and the homotopy perturbation method, Phys. Scr. 81 (2010) p. Art. No. 037002. [35] T.C. Lipscombe, Comment on ‘Application of the homotopy method for analytical solution of non-Newtonian channel flows’, Phys. Scr. 81 (2010) p. Art. No. 037001. [36] D.N. Arnold, Integrity under attack: the state of scholarly publishing, SIAM News 42 (2009) 3–4. [37] D.N. Arnold, K.K. Fowler, Nefarious numbers, Notices Amer. Math. Soc. 58 (2011) 434–437. [38] K. Parand, M. Shahini, Rational Chebyshev pseudospectral approach for solving Thomas–Fermi equation, Phys. Lett. A 373 (2009) 210–213. [39] J.P. Boyd, Orthogonal rational functions on a semi-infinite interval, J. Comput. Phys. 70 (1987) 63–88. [40] L.N. Trefethen, Spectral Methods in Matlab, Society for Industrial and Applied Mathematics, Philadelphia, 2000. [41] B. Fornberg, A Practical Guide to Pseudospectral Methods, Cambridge University Press, New York, 1996. [42] J.P. Boyd, Chebyshev and Fourier Spectral Methods, second ed., Dover, Mineola, New York, 2001, p. 665, Heavily revised and updated second edition of Boyd (1989). [43] J.P. Boyd, C. Rangan, P.H. Buchsbaum, 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 expansion, J. Comput. Phys. 188 (2003) 56–74. [44] J.P. Boyd, Large-degree asymptotics and exponential asymptotics for Fourier coefficients and transforms, Chebyshev and other spectral coefficients, J. Engrg. Math. 63 (2009) 355–399. [45] D. Elliott, The evaluation and estimation of the coefficients in the Chebyshev series expansion of a function, Math. Comp. 18 (1964) 274–284. [46] J.P. Boyd, Spectral methods using rational basis functions on an infinite interval, J. Comput. Phys. 69 (1987) 112–142. [47] E.B. Baker, The application of the Fermi–Thomas statistical model to the calculation of potential distribution in positive ions, Quart. Appl. Math. 36 (1930) 630–647. [48] C.A. Coulson, N.H. March, Momenta in atoms using the Thomas–Fermi method, Proc. R. Phys. Soc. A 63 (1950) 367–374. [49] S. Kobayashi, T. Matsukuma, S. Nagi, K. Umeda, Accurate value of the initial slope of the ordinary T–F function, J. Phys. Soc. Japan 10 (1955) 759–762. [50] F.W.J. Olver, D.W. Lozier, R.F. Boisvert, C.W. Clark (Eds.), NIST Handbook of Mathematical Functions, Cambridge University Press, New York, 2010. [51] R. Wong, Y.Q. Zhao, On a uniform treatment of Darboux’s method, Constr. Approx. 21 (2005) 225–255. [52] E. Hille, Some aspects of the Thomas–Fermi equation, J. Anal. Math. 23 (1970) 147–170. [53] J.P. Boyd, The Blasius function in the complex plane, J. Experiment. Math. 8 (1999) 381–394. [54] V.B. Mandelzweig, F. Tabakin, Quasilinearization approach to nonlinear problems in physics with application to nonlinear ODEs, Comput. Phys. Comm. 141 (2001) 268–281. [55] N.A. Zaitsev, I.V. Matyushkin, D.V. Shamonov, Numerical solution of the Thomas–Fermi equation for the centrally symmetric atom, Russian Microelectron. 33 (2004) 303–309. [56] V. Bush, The differential analyzer. A new machine for solving differential equations, J. Franklin Inst. 212 (1931) 447–488. [57] C. Miranda, Teoremi e metodi per l’integrazione numerica della equazione differenziale di Fermi, Memorie della Reale Accademia d’Italia, Classe di scienze fisiche, Mat. Nat. 5 (1934) 285–322. [58] J.P. Boyd, The asymptotic Chebyshev coefficients for functions with logarithmic endpoint singularities, Appl. Math. Comput. 29 (1989) 49–67.
Contents lists available at SciVerse ScienceDirect
Journal of Computational and Applied Mathematics journal homepage: www.elsevier.com/locate/cam
Rational Chebyshev series for the Thomas–Fermi function: Endpoint singularities and spectral methods John P. Boyd Department of Atmospheric, Oceanic and Space Science, University of Michigan, 2455 Hayward Avenue, Ann Arbor, MI 48109, United States
article
abstract
info
Article history: Received 21 September 2011 Received in revised form 19 February 2012 Keywords: Pseudospectral Chebyshev polynomials Rational Cebyshev functions Thomas–Fermi Orthogonal rational functions
√
We solve the Thomas–Fermi problem for neutral atoms, uyy − (1/ y)u3/2 = 0 on y ∈ [0, ∞] with u(0) = 1 and u(∞) = 0, using rational Chebyshev functions TLn (y; L) to illustrate some themes in solving differential equations on a semi-infinite interval. L is a user-choosable numerical parameter. The Thomas–Fermi equation is singular at the origin, giving a TL convergence rate of only √fourth order, but this can be removed by the change √ of variables, z = y with v(z ) = u(y(z )). The function v(z ) decays as z → ∞ with a term in z −3 , which is consistent with a geometric rate of convergence. However, the asymptotic series has additional terms with irrational fractional powers beginning with z −4.544 . In spite of the faster spatial decay, the irrational powers degrade the convergence rate to slightly larger than tenth order. This vividly illustrates the subtle connection between the spatial decay of u(x) and the decay-with-degree of its rational Chebyshev series. The TL coefficients an (L) are hostages to a tug-of-war between a singularity on the negative real axis, which gives a geometric rate of convergence that slows with increasing L, and the slow inverse power decay for large z, which gives quasi-tenth order convergence with a proportionality constant that decreases inversely as a power of L. For L = 2, we can approximate uy (0) (= vzz (0)) to 1 part in a million with a truncation N of only 20. L = 64 and N = 600 gives uy (0) = −1.5880710226113753127186845, correct to 25 decimal places. © 2012 Elsevier B.V. All rights reserved.
1. Introduction The Thomas–Fermi problem is 1 uyy − √ u3/2 = 0, y
u(0) = 1,
u(∞) = 0,
y ∈ [0, ∞]
(1)
where a subscripted coordinate denotes differentiation with respect to that variable. Because of its importance to theoretical physics, computing its solutions has attracted the attention of the Nobel laureates John Slater (chemistry) [1] and Richard Feynman (physics) [2] and of course Enrico Fermi [3]. Vannevar Bush, who created an analog computer known as the Differential Analyser, thus solved the Thomas–Fermi equation by computer more than a decade before the first electronic computer was built [4]; later Bush was effectively the first Presidential Science Advisor and the father of the National Science Foundation. Computing pioneer Nicholas Metropolis and physicist Edward Teller, a major figure in the development of the hydrogen bomb, also attacked this problem [2]. Hans Bethe [5] and Julian Schwinger [6,7] are other Nobel laureates who have enlarged Thomas–Fermi theory. This attention is in some ways remarkable because the neutral atom Thomas–Fermi problem is parameter-free: the Thomas–Fermi function u(y) is completely described by a single graph or table. Arnold Sommerfeld, the central figure in the
E-mail addresses: [email protected], [email protected]. 0377-0427/$ – see front matter © 2012 Elsevier B.V. All rights reserved. doi:10.1016/j.cam.2012.11.015
J.P. Boyd / Journal of Computational and Applied Mathematics 244 (2013) 90–101
91
‘‘old quantum’’ theory between Bohr’s 1913 hydrogen atom model and the dawn of quantum mechanics in 1926, gave the first analytical approximation [8,9]: 1 u(y) = 3.886 . 1 + (y/5.2415)0.772
(2)
John Mason [10] gave the improved approximation u(y) ≈
√
1 + 1.81061 y + 0.60112y
√
2
1 + 1.81061 y + 1.39515y + 0.77112y3/2 + 0.21465y2 + 0.04793y5/2
(3)
which has an absolute error for all y ∈ [0, ∞] which is smaller than 0.000012. Allan McLeod [11] has computed the Chebyshev coefficients of two series, one for y < 4 and the other for larger y, that yield the Thomas–Fermi function to about ten decimal place accuracy using only 17 terms of either series. However, precisely because of its simplicity, the Thomas–Fermi equation has been used as a testbed for a wide variety of semi-analytical and numerical schemes as exemplified by the delta-perturbation theory [12–15], Adomian decomposition [16–19], variational methods [20–23], homotopy perturbation analysis [24–29], Padé approximants [30,31] and even a technique from seventy-year old notes of Ettore Majorana [32]. We shall not review these efforts in detail because much of these merit the harsh label of Oscar Zariski as ‘‘spiccia matematica’’, that is, ‘‘small change mathematics’’. A modest portion of this work will be retroactively justified when techniques first tested on the Thomas–Fermi problem are successfully applied at the research frontier. Unfortunately, most semi-analytical methods applied to the Thomas–Fermi problem fail for anything more elaborate than a simple nonlinear ordinary differential equation. Fernandez [31,33,34] and Lipscombe [35] offer detailed criticism. The author of [22] is skewered for inflating the impact factor of his journal and of his personal citation record in [36,37]. Our reason for adding to the already vast literature on the Thomas–Fermi equation is that this problem is an illuminating example for Chebyshev spectral methods because u(y) is singular at both endpoints. For functions analytic on the entire interval, spectral methods typically converge exponentially fast with the number N of numerical degrees of freedom. However, endpoint singularities wreck this ‘‘spectral accuracy’’. McLeod [11] successfully finessed both difficulties, but did not discuss the singularities or the full range of options for coping with them. Parand and Shahini [38] obtained good accuracy using rational Chebyshev functions, as here, but were sketchy in their treatment of the singularity at the origin, completely ignored the difficulty as y → ∞, and only used a maximum of N = 12. In this article, we shall systematically treat both difficulties and obtain a single series of rational Chebyshev functions TLn [39] that converges rapidly for all of y ∈ [0, ∞]. We further shall provide something lacking from previous discussions of the TL basis: how decay as an inverse power of y for large y and a singularity at the origin jointly control the rate of convergence of the spectral series. 2. Rational Chebyshev functions for the semi-infinite interval The merits and mechanics of pseudospectral methods in general are discussed in [40–42]. As the author has previously explained in [39,43], rational Chebyshev functions are a good spectral basis for the semi-infinite interval because their connection with ordinary Fourier series yields simple theory, an exponential rate of convergence for functions analytic on the interval with sufficiently fast decay, and simplicity of programming. These basis functions are defined by TLn (y; L) ≡ Tn (x) = cos(nt )
(4)
where the argument of the rational Chebyshev functions, y, is related to the argument of the Chebyshev polynomials, x, and the argument of the trigonometric functions, t, via y=
L(1 + x) 1−x
y = L cot
2
↔x=
y−L
t
2
(5)
y+L
↔ t = 2 arccot
y L
.
(6)
The first few functions are given explicitly in Table 1. The pseudospectral interpolation points are yi = L cot
2
ti
2
↔ ti ≡
(2i − 1) π 2N + 2
i = 1, . . . , N + 1.
(7)
The Thomas–Fermi solution is represented as u(y) ≈
N −1 n=0
an TLn (y; L)
(8)
92
J.P. Boyd / Journal of Computational and Applied Mathematics 244 (2013) 90–101 Table 1 Rational Chebyshev functions for the semi-infinite interval: TLn (y). n
TLn (y; L)
0 1 2 3 4
(y − L)/(y + L) (y2 − 6yL + L2 )/(y + L)2 (y − L)(y2 − 14yL + L2 )/(y + L)3 (y4 − 28Ly3 + 70L2 y2 − 28L3 y + L4 )/(y + L)4
1
where the pseudospectral procedure for computing the an is described in [39,43] and the Maple code in the Appendix. Derivatives of the basis functions are computed applying the chain rule to derivatives of cos((n − 1)t ) to convert t-derivatives to y-derivatives. The constant L is a user-choosable ‘‘map parameter’’. Although there are sophisticated ways to estimate the best choice of L [39], in practice it is usual to begin with L equal to the expected dominant length scale of u(y), and then experiment. The criterion for ‘‘optimum’’ is rate of convergence, which can be evaluated from a plot of the computed coefficients. Note that both theory and the experiments below show that (i) the optimum L increases significantly with increasing truncation N of the TL series and (ii) one can obtain very accurate approximations for a given N with L anywhere within a rather large range. The following pair of theorems have not been previously published. Collectively, they describe how behavior at the endpoints, y = 0 and y → ∞, controls the rate of TL series convergence. Theorem 1 (Singularities at the Endpoint y = 0). Singularities at the endpoint y = 0 have precisely the same effect on convergence of the series of rational Chebyshev functions TL as they have the convergence of the Chebyshev polynomial series on a finite interval x ∈ [0, b]. A list of various species of singularities and their rates of convergence is given at the end of Chapter 2 of [42] and also [44]. In particular, if u(y) has a branch point at the origin with the singular term proportional to yψ , then its rational Chebyshev coefficients will satisfy an ∼ constant /n2ψ+1
(9)
as follows from [45]. Proof. When y is sufficiently small, the rational Chebyshev to Chebyshev polynomial change-of-coordinate, y = L(1 + x)/(1 − x), simplifies through power series expansion to y ≈ (1 + x)
L 2
+ O([1 + x]2 ).
(10)
Thus, the relationship between the two coordinates is linear near the origin and is therefore no map at all except for an irrelevant rescaling through the map parameter L. In other words, TLn (y; L) ∼ Tn
2 −1 + y , L
y≪1
(11)
for small y. (Fig. 1 illustrates this asymptotic relation.) It then follows that the effects of singularities on rational Chebyshev and Chebyshev polynomials must be the same since for small y, they are the same. Theorem 2 (Convergence Effects of Inverse Power Decay). 1. If a function u(y), not restricted to the Thomas–Fermi function, decays asymptotically as O(1/yk ) where k is an integer and is free of singularities everywhere on y ∈ [0, ∞] except at infinity and also that the asymptotic expansion for large y contains only non-positive integral powers, then the rational Chebyshev series will converge exponentially fast. 2. If u(y) has an asymptotic expansion with y−k as the power of smallest non-integral k, then the rational Chebyshev coefficients will decay as an ∼
constant n2k+1
.
(12)
Proof. The point y = ∞ is mapped into t = 0. In this neighborhood, the mapping simplifies to y = L cot2
t
2
≈L
4 t2
,
t ≪ 1.
(13)
This implies that for any positive k y−k ∼ t 2k .
(14)
J.P. Boyd / Journal of Computational and Applied Mathematics 244 (2013) 90–101
93
TL & TN, n=100 1
0.5
0
-0.5
-1 0.05
0
0.1 y
0.15
0.2
Fig. 1. Comparison of TL100 (y; L) (solid) with T100 (−1 + (2/L)y) (dashed) and the difference between them (dotted) for L = 3. The two functions closely approximate each other only for rather small y, but this is sufficient to give both bases identical responses to singularities at y = 0.
Thus, a function whose leading asymptotic term in a non-positive integer k will spawn a zero of order 2k at t = 0. The function t 2k is symmetric about t = 0 and thus poses no problems for approximating u(y(t )) as a Fourier cosine series whose terms are all individually symmetric, too. This proves the first proposition, which is that the TL series will converge exponentially if k is an integer, and there are no singularities elsewhere on the real axis. If k is not an integer, however, then u(y(t )) will have a branch point at t = 0. Elliott’s analysis of the consequences of a branch point of exponent 2k in t proves the second proposition [45]. One might suppose that half-integral powers of k would also yield exponential convergence because the endpoint behavior of u(y(t )) is an integral power of t and therefore not a branch point. However, odd powers of t are antisymmetric with respect to t = 0, i.e., t 2m+1 = −(−t )2m+1 , and this is incompatible with a cosine basis. As illustrated for the related case on an infinite interval [46], exponential convergence can be retrieved by converting to the trigonometric basis and then applying a sine basis. Fig. 2 confirms this analysis: integral powers of y for large y allow an exponential rate of convergence whereas the convergence rate is algebraic for fractional powers of y. The example u(y) = 1/(1 + y)k is atypical in that it is not singular at infinity so that the rate of convergence for integer k is geometric, that is, proportional to exp(−qn) for some constant q. It is more common in applications for functions to be singular at infinity, or to have an infinite number of singularities with a limit at infinity, and then the convergence is ‘‘subgeometric’’ in the terminology of [42], that is, proportional to exp(−qnr ) for some constant r < 1, typically r = 1/2 or r = 2/3. We shall not elaborate further except to refer to Chapter 17 of [42,46] because the Thomas–Fermi problem has relatively strong singularities at infinity that make subgeometric convergence irrelevant. 3. Branch point at the origin 3.1. Transforming the problem The Thomas–Fermi neutral atom problem is 1 uyy = √ u3/2 = 0, y
u(0) = 1,
u(∞) = 0,
y ∈ [0, ∞].
(15)
Recalling that u(y) cannot vanish at the origin because of the boundary condition u(0) = 1, it is obvious that the second √ derivative of u(y) must be unbounded in the limit y → 0 because the second derivative is proportional to 1/ y. In 1930, a graduate student at the University of Michigan, Edward Baker, gave a complete description of the singularity by deriving the small-y series [47]: u = 1 + sy +
4 3
y3/2 +
2 5
s y5/2 + s2
3 70
y7/2 +
1 3
y3 + · · ·
where s is the value of the first derivative at the origin, shown numerically to be roughly s ≈ −1.58.
(16)
94
J.P. Boyd / Journal of Computational and Applied Mathematics 244 (2013) 90–101
SEMILOGY, integer k 100
100
LOGLOG, k=fraction
k=1 k=2 k=3
k=3/2 k=5/2 k=7/2
10-5
10-5
10-10
10-10
10-15
0
10
10-15 100
20
102 degree n
degree n
Fig. 2. The TL coefficients of the function (1 + y)−k for integer k (left graph) and half-integer k (right). With a logarithmic vertical scale and linear horizontal scale, geometric convergence appears linear; the left plot therefore shows the exponential convergence when k is an integer. On a log–log plot, decay proportional to a power of k is linear. The dotted guidelines show that the coefficients for y3/2 (top curve), y5/2 (middle curve) and y7/2 (bottom curve) asymptote to n−4 , n−6 and n−8 , respectively.
TL coefficients of the Thomas-Fermi function u(y) 100 coefficients, L=2 1.8/n**4
10-5
10-10 100
101 degree n
102
Fig. 3. Coefficients of the TL rational Chebyshev expansion for u(y). The dashed line is a guideline with slope proportional to n−4 .
Theorem 1 shows that if a function u(y) has a branch point of order k at the origin, then the TL coefficients will decay proportionally to n−(2k+1) . The leading singular term in Baker’s series is y3/2 , so the predicted rate of convergence is fourth order, that is, |an | ∼ constant/n4 . The TL series for the Thomas–Fermi function (8) does indeed show the expected fourth order rate of convergence as illustrated in Fig. 3. As explained in [42], a spectral series with coefficients decreasing as n−m will yield a maximum pointwise error decreasing one order slower as N −(m−1) where N is the truncation of the spectral series. Thus, the rational Chebyshev series for u(y) has only third order convergence. The first derivative of the n-th TL function at the origin is O(n2 ), so the rate of convergence to the correct value of s = uy (0) is only first order: numerical experiments showed that with L = 4, the error in uy (0) was about 0.174/N where N is the truncation of the TL series. One strategy adopted by McLeod [11] is to apply a spectral method in the coordinate z≡
√
y.
(17)
We prefer the approach of Fernandez [31], who also adopts the new unknown
v(z ) ≡
u(z 2 ).
(18)
The problem becomes z vvzz + (vz )2 − vvz − 2z 2 v 3 = 0,
This is the form we solved numerically.
v(0) = 1,
v(∞) = 0.
(19)
J.P. Boyd / Journal of Computational and Applied Mathematics 244 (2013) 90–101
95
TL coefficients of the Mason approx. to sqrt(u(z)) 100 coefficients, L=1 8 exp(−n) 10-5
10-10
10-15 0
10
20 degree n
30
40
Fig. 4. Solid: absolute values of the rational Chebyshev series for Mason’s quadratic-polynomial-divided-by-quintic polynomial approximation to v(z ); this approximation was used to initialize the Newton–Kantorovich iteration. The dashed line is a guideline that is the graph of exp(−n).
3.2. Power series The transformed differential equation is still singular at the origin; it is merely the solution that is desingularized. Because (19) is singular, the standard Maple command to solve a differential equation as a series returns only a blank. When a series of the form
v(z ) = 1 + a1 z + (1/2)sz 2 + a3 z 3 + a4 z 4 + · · ·
(20)
is substituted into (19), the residual is of the from z vvzz + (vz )2 − vvz − 2z 2 v 3 = −a1 + {3a3 + (3/2)a1 s − 2} z 2 + 8a4 + 8a1 a3 + s2 − 6a1 z 3 + · · · .
(21)
Because the constant in the residual is −a1 , it is impossible to satisfy the differential equation with a power series unless a1 = 0 ↔ vz (0) = 0.
(22)
The remaining coefficients can be found recursively in terms of s ≡ vzz (0); the specific choice s = −1.58807 . . . is the only choice that allows the sum of the power series to satisfy the boundary condition v(∞) = 0. Because v(z ) is not uniquely specified by the initial conditions v(0) = 1, vz (0) = 0, it is not possible to begin a standard initial value calculation at the origin itself. The best one can do is to march away from a point near the origin where the necessary initial values can be obtained from the power series. Remarkably, however, the TL pseudospectral method has no difficulties with the boundary value problem for v(z ). As is typical when the differential equation is singular at an endpoint, but the solution is not, the spectral method automatically picks out the ‘‘nice’’ solution without special tricks, or the imposition of additional, nonstandard boundary conditions. 4. Asymptotic series for large y and its consequences Arnold Sommerfeld, who supervised more Nobel laureates than any other physicist, derived his analytic approximation by investigating the asymptotic behavior of u(y) for large y. Charles Coulson, professor of theoretical chemistry at Oxford, and a student elaborated this into an asymptotic series of the form [48] u(y) ∼
144 y3
F F2 F3 1 − λ + 0.62569 2λ + 0.3138 3λ + · · · y y y
(23)
where F = 13.27094 [49] and
λ = (1/2)(−7 + −3
√
73) = 0.772001872.
(24)
The leading order y does not prevent even a geometric rate of convergence. This is illustrated in Fig. 4, which shows the geometric convergence of the TL series for Mason’s approximation to v(z ), which reproduces the leading order asymptotic term but omits all fractional powers. (By geometric decay, we mean convergence proportional to exp(−µn) for some n constant µ, like that of the terms in the geometric series n x .) However, Theorem 2 predicts that the next term of −3.772 O(y ) will, ignoring the effects of singularities at y = 0, slow the convergence of the TL coefficients to O(n−8.544 ).
96
J.P. Boyd / Journal of Computational and Applied Mathematics 244 (2013) 90–101
abs(TL coeffs)
TL coefficients for L=1 & L=64
10-10
10-20
L=1
L=64
10-30 100
200
300 degree n
400
500
Fig. 5. The absolute of the TL coefficients an of the Thomas–Fermi function for L = 1 and L = 64.
However, we actually solve Fernandez’ form. The asymptotic series for v(z ) is
v(z ) ∼ ∼ ∼
144
z6 12
1−
z3 12 z3
−
F
1−
z 2λ
F z 2λ
6F z 3+2λ
+ 0.62569
+ 0.62569
+ ···.
F2 z 4λ
F2 z 4λ
+ ···
+ ···
(25)
(26) (27)
The leading non-integral power is 1/z 3+2λ = 1/z 4.544 . This implies that the TL series for v(z ) has coefficients asymptotically proportional to an ∼
constant n10.088
.
(28)
Thus, the TL series for v(z ) has higher order convergence than the corresponding series for u(y). It is important to note that (28) is an asymptotic relation that holds only for sufficiently large n. Fig. 5 shows that the TL coefficients computed with two different values of L. On this plot with a logarithmic vertical scale and a linear horizontal scale, geometric decay appears linear. For both choices of L, the coefficients of small degree do in fact decay geometrically. For large n, however, the decay with degree slows dramatically. Fig. 6 shows the same two sets of coefficients, but with a logarithmic scale on both axes. Power law decay appears linear on a log–log plot. For both values of L, the coefficients asymptote to lines paralleling the theoretical prediction of n−10.088 which is shown as the dashed guideline. Fig. 7 is a different graph showing more precisely how the coefficients asymptote to a simple power law. Why does L affect the coefficients this way? And why is there a dramatic break in slope on the coefficient graphs? Darboux’ Principle, discussed in Section 2.10 of [50] and also [51], asserts that the asymptotic coefficients an of a power series are a sum of contributions, one from each of the singularities in the complex plane of the function represented by the power series. This principle extends to Fourier and all species of Chebyshev series also. Singularities on the expansion interval contribute terms that decay more slowly than geometrically—here, the fractional power decay of v(z ) is a branch point at z = ∞ that yields a contribution to an that decays slightly faster than tenth order. Singularities off the expansion interval contribute terms which decay exponentially at a rate which is controlled entirely by the location of the singularity; the exponential in degree n is multiplied by a slower-than-exponential factor, such as powers or logarithms of n or something more complicated, that depends only on the type of singularity (i.e., pole of some order or logarithm or other branch point). It is easy to compute the power series of v(z ) to high order. The coefficients strictly alternate in sign, suggesting that the convergence-limiting singularity is on the negative real axis. The usual ‘‘ratio test’’, an−1 /an = zs + O(1/n), suggests that the singularity is at z = zs ≈ −0.665. However, Hille [52] wrote (p. 158): ‘‘Formula (4.2) as well as Fowler’s result suggest that y(x) has a pole of order four . . . A direct substitution of a Laurent series . . . leads to unexpected complications and unmanageable computations . . . the singularity is more complicated than a pole’’. The same is true of the singularities of the Blasius function and could probably be analyzed by the methods of [53], but this is beyond the scope of this article.
J.P. Boyd / Journal of Computational and Applied Mathematics 244 (2013) 90–101
97
Log−log plot: TL coeffs. for L=1 & L=64
L=1
L=64
abs(TL coeffs)
10-10
0.0 01 /n 10
L=1
.08
10-20
8
L=64 10-30 100
101
102 degree n
Fig. 6. Same as previous figure, but using a logarithmic scale for both axes. The dashed line is a guideline with a slope of n−10.088 .
a(n)* n**10.088, normalized L=1
scaled coefficients
1.05
1
0.95 0
200
400
600
degree n Fig. 7. The solid curve is n10.088 an for L = 1, normalized by its value for n = 350. The coefficients of low degree do not follow the predicted asymptotic pattern because the asymptotic approximation only applies for sufficiently large n. The coefficients on the extreme right deviate from the predicted scaled magnitude because of aliasing errors, a general phenomenon with interpolatory spectral methods.
Fortunately, the location of the singularity controls exponential terms in the spectral coefficients [42,44]; the order and type of the singularity are weak (non-exponential) influences on the rate of convergence. If the singularity in z is sufficiently close to the origin, then z = L cot2 (t /2) ≈ (L/4)(t − π )2 + O([L/24])(t − π )4 .
(29)
If the TL-convergence-limiting singularity is on the negative real axis, then
ts = π ± i
4
|zs |. L This will generate a contribution to the asymptotic TL coefficients proportional to exp(−nℑ(ts )), that is, an ∼ A(n) exp(−n
4|zs |/L) + terms from other singularities
(30)
(31)
where A(n) is a factor that varies more slowly with n than the exponential, such as some combination of constants, powers of n and logarithms of n. This shows explicitly that increasing the map parameter L slows the convergence of the contributions of this term. It predicts that in the geometrically-converging regime, the L = 64 curve should fall proportional to about exp(−0.2), or in other words diminish by about 10−9 for each additional 100 terms in the series. In Fig. 5, this is close to what is observed for n < 300 where the coefficients transition to dominance by the tenth order contributions from the singularity at z = ∞.
98
J.P. Boyd / Journal of Computational and Applied Mathematics 244 (2013) 90–101 Table 2 Selected calculations of uy (0) (=vzz (0)); inaccurate digits claimed to be correct by their authors are in bold face. Entries up to 1955 from Kobayashi et al. Fermi (1928) [3] Baker (1930) [47] Bush and Caldwell (1931) [4] Miranda (1934) [57] Slater and Krutter (1935) [1] Feynman, Metropolis and Teller (1949) [2] Kobayashi et al. (1955) [49] Fernandez (2008) [31] Boyd (2011) [this article] (25 decimal places) ’’ [N = 20, L = 2] ’’ [N = 40, L = 4] ’’ [N = 100, L = 16] ’’ [N = 200, L = 32]
−1.58 −1.588558 −1.589 (analog computer [56]) −1.5880464 −1.58808 −1.58875 −1.5880710 −1.58807102261137531 −1.5880710226113753127186845 −1.58807 −1.5880710226 −1.588071022611375 −1.588071022611375312719
5. Iteration and elimination Because the differential equation is nonlinear, we employed a Newton–Kantorovich iteration [42], also known as ‘‘quasilinearization’’ [54], to reduce the problem to a sequence of linear differential equations: z vn δn,zz + (2z vn,z − v) δn,z + z vn,zz − vn,z − 6z 2 vn2 δn = z vn vn,zz + vn,z vn,z − vn vn,z − 2z 2 vn3
(32)
vn+1 = vn − δn
(33)
subject to the boundary conditions δn (0) = vn (0) − 1 and δn (∞) = 0. The differential equation was solved by truncating the TL series to N terms and then applying collocation, also known as the pseudospectral method, at (N − 2) points. The collocation conditions plus the two boundary conditions gives an N × N dense matrix problem which can be solved by Gaussian elimination, alias Crout reduction, for the vector of TL coefficients. We omit the details since these are discussed with many examples in [39,42]. The computational cost and memory can be greatly reduced by using a preconditioned iteration as explained in [42]. Since an N = 600 computation using direct matrix solves required only about 2.5 h on a 2009 laptop running Maple 13 in Windows Vista, we did not bother with such refinements. However, an inner iteration is very valuable for more difficult problems. The Newton–Kantorovich iteration requires an initialization or ‘‘first guess’’ v0 . We employed Mason’s rational approximation, and found that the residual (the right hand side of the equation) for δn was always reduced to a maximum pointwise value smaller than 10−31 after five iterations without any underrelaxation. The program is so simple we have included the entire Maple code in the Appendix. 6. Approximating the slope at the origin The first derivative of u(y) at left endpoint has always been of great interest because this one number is sufficient to convert the boundary value problem into a second order initial value problem with the initial conditions u(0) = 1 and uy (0) = s where s is a very accurate approximation to the first derivative at the origin. In principle, the initial value problem can be solved by a simple march from y = 0 to y → ∞ using one’s favorite Runge–Kutta or Adams–Bashforth scheme. In practice, such marching is ill-conditioned as reviewed most recently by Zaitsev, Matyushkin, and Shamonov [55]. Table 2 is only a partial list of the vast number of calculations of the Thomas–Fermi function, but the entries show that accuracy has increased with time but not monotonically; some very distinguished authors in fact obtained fewer correct digits than they thought. Kobayashi et al. and Fernandez deserve credit because all the digits they claimed are correct. The last four rows show the best approximation obtained from a given number N of terms of the TL basis. Note that the optimum L increases with N; this is quite typical in calculations with rational Chebyshev functions on both the infinite and semi-infinite intervals. Although we actually computed v(z ) and not u(y), it is easy to show that s ≡ uy (0) ≡ vzz (0).
(34)
Table 3 shows the errors in vzz (0) for various L and N. It is clear that small L is best for small N while large L is better for large N but awful for small N. Because the first derivative of TLn at z = 0 is O(n2 ) and the second derivative is O(n4 ), it is not surprising that the error in vzz (0) converges at a slower rate than the coefficients an . The table shows a seventh order rate of convergence instead of the expected sixth order. We do not have an explanation. As is usual when an exact solution is unknown, a solution of very high N, optimum L, and large value of the Maple precision parameter ‘‘Digits’’ was taken as the ‘‘exact’’ solution for computing errors. Our benchmark employed N = 650, L = 64, and Digits = 45 and required three and a quarter hours computation on our laptop. It was judged, by consistency with other
J.P. Boyd / Journal of Computational and Applied Mathematics 244 (2013) 90–101
99
Table 3 Computed values of uy (0) = vzz (0) for various L and N; the optimum values in a given row or column are in boldface. N
L=1
L=2
L=4
L=8
L = 16
L = 32
20 40 100 200 400 600
0.101e−3 0.303e−6 0.415e−9 0.293e−11 0.210e−13 0.118e−14
0.142e−5 0.343e−8 0.458e−11 0.319e−13 0.227e−15 0.127e−16
0.526e−5 0.390e−10 0.505e−13 0.349e−15 −0.245e−17 0.137e−18
0.961e−4 −0.237e−8 0.563e−15 −0.378e−17 −0.265e−19 0.147e−20
−0.526e−3 −0.226e−5 0.506e−16 −0.414e−19 −0.287e−21 −0.159e−22
−0.145e−3 −0.886e−11 −0.452e−21 −0.313e−23 −0.181e−24
0.0536
L = 64
−0.216 −0.720e−2 −0.342e−7 −0.627e−16 < 1.e − 25 < 1.e − 25
runs with different N , L and Digits, that at least the first twenty-five decimal places are correct, these being invariant among our highest resolution computations. 7. Summary The singularity of the Thomas–Fermi function at the origin, which would otherwise degrade convergence of the rational Chebyshev series to fourth order, can be eliminated by a simple transformation of the coordinate and unknown. The term z −4.544 in the asymptotic series for v(z ) slows the convergence to tenth order. The TL functions depend up a user-choosable constant, the ‘‘map parameter’’ L. The asymptotic TLn coefficient an as n → ∞ is the sum of contributions from the singularities of v(z ). Increasing L slows the rate of geometric convergence of the term in an due to the pole-and-branch-point on the negative real axis, but it greatly diminishes the proportionality constant in the tenth order contribution from the singularity at z = ∞. The result of this tension between two singularities which are affected oppositely by L is that a graph of the TL coefficients shows two regimes: geometrically fast (exp(−nµ)) decay for small n and slower, tenth order decay for large n. The ‘‘crossover point’’ between regimes grows rapidly with L. The result is that when the TL series is truncated to just twenty terms, the best choice is L = 2. When N = 600, however, a fast exponential decay is not important because N is so huge. It is better to pick L = 64 or larger to diminish the effect of the singularity at z = ∞. It is then possible to compute the second derivative of u(y) at the origin to twenty-five decimal places. The insensitivity to L, in the sense that accurate calculations can be made when L is many times larger or smaller than optimum, is a generic feature of rational Chebyshev series [46,39]. It is also common that the optimum L is not fixed, but rather grows with the truncation N of the series. The reason is the tension between resolving the behavior of v(z ) for large |z | and simultaneously resolving singularities off the expansion interval, z ∈ [0, ∞]. Increasing L places grid points at larger and larger z: good for resolving the singularity at infinity. However, a rational Chebyshev expansion can be always be transformed into an equivalent Fourier cosine series through the identity TLn (z ; L) ≡ cos(nt [z ; L]) where the trigonometric coordinate is a function of z and L only as z /L, which implies that increasing L effectively moves singularities at finite z closer to the real axis in t. Thus increasing L worsens the contributions of these singularities to the rational Chebyshev coefficients of high degree. Because the singularity at infinity usually generates a slower-than-geometric contribution to the asymptotic an – tenth order here, subgeometric proportional to exp(−qn2/3 ) more commonly – increasing N requires, for optimality, increasing L to better cope with the slower-decreasing contributions. In theory, the tenth order convergence can be replaced by asymptotic exponential convergence by using an exponential mapping. Our experiments here produced no improvement and are therefore omitted. As discussed in [58], mappings introduce additional variability that slows the convergence of the low degree coefficients while accelerating coefficients an for n > ncrossover for some ‘‘crossover’’ value ncrossover . Because the TL series for v is already converging at better than tenth order, the ‘‘crossover’’ point here is very large, and an exponential mapping does no good until the degree is in the thousands. Acknowledgments This work was supported by the National Science Foundation through grants OCE 0451951 and ATM 0723440. I thank Paolo Amore for his careful reading and comments. I also thank the three reviewers for their constructive suggestions. Appendix. Maple code to solve the Thomas–Fermi equation restart; L:= 2; # map parameter N:= 600; # number of basis functions & grid points; Digits:= 40; # number of decimal digits of precision itermax:= 6; # number of Newton iterations; ncol:= N−2; # number of interior collocations points; Mason25:= proc(y) local yh,vn,vd,Mas; # evaluates Mason’s approximation; yh:= sqrt(y); vn:= 1+1.81061*yh+0.60112*y;
100
J.P. Boyd / Journal of Computational and Applied Mathematics 244 (2013) 90–101
vd:= 1+1.81061*yh+1.39515*y+0.77211*y*yh+0.21465*y*y+ 0.04793*yh*y*y; vM:= (vn/vd); end proc: with(LinearAlgebra); fa:= Vector(N,orientation=column): v:= Vector(N,orientation=column): za:= Vector(N,orientation=column): a:= Vector(N,orientation=column): a0:= Vector(N,orientation=column); # TL differentiation matrices D0,D1,D2 for 0th, d/dz, d**2/dz**2; D0:= Matrix(ncol,N): D1:= Matrix(ncol,N): D2:= Matrix(ncol,N): # collocation points za; TL computed by trigonometric identity; for ii from 1 to ncol do ta[ii]:= evalf(Pi*(2*ii−1) / (2*ncol)); tt:= ta[ii]; za[ii]:= evalf(L*cot(tt/2)**2); for j from 1 to N do TL:= evalf(cos((j−1)*tt)); # Compute derivatives using trigonometric definition and chain rule; PT:= evalf(− (j−1)*sin((j−1)*tt)); PTT:= evalf(− (j−1)*(j−1)*TL); S := evalf(sin(tt/2)); C := evalf(cos(tt/2)); TLz:= ev alf(− S * S * S * PT/ (L * C)); TLzz:= evalf((S**5)*(2*C*S*PTT+(3-2*S*S)*PT)/(L*L * C*C*C* 2)); D0[ii,j] := TL ; D1[ii,j]:= TLz ; D2[ii,j]:= TLzz; od: od: # Expand initialization v0 as TL series with coeffs a0=TLEXPMat*Mason25(za); for ii from 1 by 1 to ncol do v0[ii]:= evalf(Mason25(za[ii]*za[ii])); od: TLEXPMat:= Matrix(ncol,ncol): for ii from 1 to ncol do for j from 1 to ncol do TLEXPMat[ii,j]:= evalf((2/ncol)*cos((ii−1)*ta[j])); od; od: for j from 1 to ncol do TLEXPMat[1,j]:= evalf(TLEXPMat[1,j]*0.5); od: for ii from 1 by 1 to ncol do a0[ii]:= 0; for j from 1 by 1 to ncol do a0[ii]:= evalf(a0[ii]+TLEXPMat[ii,j]*v0[j]); od: od: a0[N−1]:= 0; a0[N]:= 0; a:= a0: # Initialize Jacobian matrix and compute its last two rows; Jacobian:= Matrix(N,N):for j from 1 by 1 to N do Jacobian[N−1,j]:= evalf(cos((j−1)*Pi));Jacobian[N,j]:= 1; od: # Start of Newton–Kantorovich iteration; counter is ‘‘iter’’; for iter from 1 by 1 to itermax do for ii from 1 by 1 to ncol do v:= 0; vz:= 0; vzz:= 0; for j from 1 by 1 to N do v:= evalf(v + D0[ii,j]*a[j]); vz:= evalf(vz + D1[ii,j]*a[j]); vzz:= evalf(vzz + D2[ii,j]*a[j]); od: # fa is residual of TF ODE; fa[ii]:= evalf((v*vzz + vz*vz)*za[ii] − v*vz − 2*za[ii]*za[ii]*v*v*v); for j from 1 by 1 to N do Jacobian[ii,j]:= evalf(za[ii]*v * D2[ii,j] + (2*za[ii]*vz − v)*D1[ii,j] +(za[ii]*vzz − vz − 6*za[ii]*za[ii]*v*v)*D0[ii,j]); od: od: # Apply boundary conditions: v(0)=1, v(infinity)=0; vatzero:= 0; vatinfinity:= 0; for j from 1 by 1 to N dovatzero:= evalf(vatzero + cos((j−1)*Pi)*a[j]); vatinfinity:= evalf(vatinfinity + a[j]); od; fa[N−1]:= vatzero − 1; fa[N]:= vatinfinity; delta_a:= LinearSolve(⟨Jacobian|fa⟩); # compute TL coeffs. delta_aof δn (z ). for j from 1 by 1 to N do a[j]:= evalf(a[j] − delta_a[j]); od: residualLoo[iter]:= max(seq(abs(fa[k]),k=1..N)); print(iter,residualLoo[iter]); od: # end of Newton iteration loop; # Only plots & output below; anpoints := [seq([j,abs(a[j])],j=1..N]: plot2:= loglogplot(anpoints,style=line):guidepoints:= [seq([j,j**(−10.022)],j=1..N)]: plotguide:= loglogplot(guidepoints): plots[display](plot2,plotguide); vzzat0:= 0; for j from 1 to N do j1:= j−1; vzzat0:= evalf(vzzat0+ a[j]*4*cos(j1*Pi)*((1/3)*j1**4+(2/3)*j1**2)/L**2); od: vxzat0ezact:= −1.5880710226113753127186845; errvzzat0:= evalf(abs(vzzat0 vzzat0exact));lprint(N,vzzat0,errvzzat0); References [1] J.C. Slater, H.M. Krutter, The Thomas–Fermi method for metals, Phys. Rev. 47 (1935) 559–568. [2] R.P. Feynman, N. Metropolis, E. Teller, Equations of state of elements based on the generalized Fermi–Thomas theory, Phys. Rev. 75 (1949) 1561–1573.
J.P. Boyd / Journal of Computational and Applied Mathematics 244 (2013) 90–101
101
[3] E. Fermi, Eine statistische Methode zur Bestimmung einiger Eigenschaften des Atoms und ihre Anwendung auf die Theorie des periodischen Systems der Elemente, Z. Phys. 48 (1928) 73–79. [4] V. Bush, S.H. Caldwell, Thomas–Fermi equation solution by the differential analyzer, Phys. Rev. 38 (1931) 1898–1902. [5] H. Bethe, Thomas–Fermi theory of nuclei, Phys. Rev. 167 (1968) 879–907. [6] J. Schwinger, Thomas–Fermi model: the leading correction, Phys. Rev. A 22 (1980) 1827–1832. [7] J. Schwinger, Thomas–Fermi model: the second correction, Phys. Rev. A 24 (1981) 2353–2361. [8] A. Sommerfeld, Asymptotische Integration der Differentialgleichung des Thomas Fermischen Atoms, Z. Phys. 78 (1932) 283. [9] A. Sommerfeld, Integrazione asintotica dell’equazione differentiale di Thomas–Fermi, Rend. R. Accad. Lincei 15 (1932) 293–308. [10] J.C. Mason, Rational approximations to the ordinary Thomas–Fermi function and its derivatives, Math. Proc. Cambridge Philos. Soc. 84 (1964) 357–360. [11] A.J. McLeod, Chebyshev series solution of the Thomas–Fermi equation, Comput. Phys. Comm. 67 (1992) 389–391. [12] 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. [13] B. Laurenzi, An analytic solution to the Thomas–Fermi equation, J. Math. Phys. 31 (1990) 2535–2537. [14] A. Cedillo, A perturbative approach of the Thomas–Fermi equation in terms of density, J. Math. Phys. 34 (1993) 2713–2717. [15] D.A. Morales, An N-dimensional interpretation of the delta-expansion for the Thomas–Fermi equation, J. Math. Phys. 8 (1994) 3916–3921. [16] G. Adomian, Solution of the Thomas–Fermi equation, Appl. Math. Lett. 11 (1998) 131–133. [17] A.M. Wazwaz, The modified decomposition method and Padé approximants for solving the Thomas–Fermi equation, Appl. Math. Comput. 105 (1999) 11–19. [18] J.I. Ramos, Piecewise-adaptive decomposition methods, Chaos Solitons Fractals 40 (2009) 1623–1636. [19] H. Chu, Y. Liu, The new ADM-Pade technique for the generalized Emden–Fowler equations, Modern Phys. Lett. B 24 (2010) 1237–1254. [20] N. Anderson, A. Arthurs, Variational solution of the Thomas–Fermi equation, Quart. Appl. Math. 39 (1981) 127–129. [21] B. Burrows, P. Core, A variational-iterative approximation solution of the Thomas–Fermi equation, Quart. Appl. Math. 42 (1984) 73–76. [22] J.H. He, Variational approach to the Thomas–Fermi equation, Appl. Math. Comput. 143 (2003) 533–535. [23] M. Desaixj, D. Anderson, M. Lisak, Variational approach to the Thomas–Fermi equation, European J. Phys. 25 (2004) 699–705. [24] S. Liao, An explicit analytic solution to the Thomas–Fermi equation, Appl. Math. Comput. 144 (2003) 495–506. [25] H. Khan, H. Xu, Series solution to the Thomas–Fermi equation, Phys. Lett. A 365 (2007) 111–115. [26] A. El-Nahhas, Analytic approximations for Thomas–Fermi equation, Acta Phys. Pol. A 114 (2008) 913–918. [27] M.A. Noor, S.T. Mohyud-Din, Homotopy perturbation method for solving nonlinear higher-order boundary value problems, Int. J. Nonlinear Sci. Numer. Simul. 9 (2008) 395–408. [28] B. Yao, A series solution to the Thomas–Fermi equation, Appl. Math. Comput. 203 (2008) 396–401. [29] S.T. Mohyud-Din, A. Yildirim, Solving nonlinear boundary value problems using He’s polynomials and Pade approximants, Math. Probl. Eng. (2009) p. Art. No. 690547. [30] L.N. Epele, H. Fanchiotti, C.A.G. Canal, J.A. Ponciano, Padé approximant approach to the Thomas–Fermi problem, Phys. Rev. A 60 (1999) 280–283. [31] F.M. Fernandez, Comment on ‘Series solution to the Thomas–Fermi equation’, Phys. Lett. A 365 (2007) 111; Phys. Lett. A 372 (2008) 5258–5260. [32] S. Esposito, Majorana solution of the Thomas–Fermi equation, Amer. J. Phys. 70 (2002) 852–856. [33] F.M. Fernandez, On some approximate methods for nonlinear models, Appl. Math. Comput. 215 (2009) 168–174. [34] F.M. Fernandez, On two new applications of the homotopy analysis method and the homotopy perturbation method, Phys. Scr. 81 (2010) p. Art. No. 037002. [35] T.C. Lipscombe, Comment on ‘Application of the homotopy method for analytical solution of non-Newtonian channel flows’, Phys. Scr. 81 (2010) p. Art. No. 037001. [36] D.N. Arnold, Integrity under attack: the state of scholarly publishing, SIAM News 42 (2009) 3–4. [37] D.N. Arnold, K.K. Fowler, Nefarious numbers, Notices Amer. Math. Soc. 58 (2011) 434–437. [38] K. Parand, M. Shahini, Rational Chebyshev pseudospectral approach for solving Thomas–Fermi equation, Phys. Lett. A 373 (2009) 210–213. [39] J.P. Boyd, Orthogonal rational functions on a semi-infinite interval, J. Comput. Phys. 70 (1987) 63–88. [40] L.N. Trefethen, Spectral Methods in Matlab, Society for Industrial and Applied Mathematics, Philadelphia, 2000. [41] B. Fornberg, A Practical Guide to Pseudospectral Methods, Cambridge University Press, New York, 1996. [42] J.P. Boyd, Chebyshev and Fourier Spectral Methods, second ed., Dover, Mineola, New York, 2001, p. 665, Heavily revised and updated second edition of Boyd (1989). [43] J.P. Boyd, C. Rangan, P.H. Buchsbaum, 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 expansion, J. Comput. Phys. 188 (2003) 56–74. [44] J.P. Boyd, Large-degree asymptotics and exponential asymptotics for Fourier coefficients and transforms, Chebyshev and other spectral coefficients, J. Engrg. Math. 63 (2009) 355–399. [45] D. Elliott, The evaluation and estimation of the coefficients in the Chebyshev series expansion of a function, Math. Comp. 18 (1964) 274–284. [46] J.P. Boyd, Spectral methods using rational basis functions on an infinite interval, J. Comput. Phys. 69 (1987) 112–142. [47] E.B. Baker, The application of the Fermi–Thomas statistical model to the calculation of potential distribution in positive ions, Quart. Appl. Math. 36 (1930) 630–647. [48] C.A. Coulson, N.H. March, Momenta in atoms using the Thomas–Fermi method, Proc. R. Phys. Soc. A 63 (1950) 367–374. [49] S. Kobayashi, T. Matsukuma, S. Nagi, K. Umeda, Accurate value of the initial slope of the ordinary T–F function, J. Phys. Soc. Japan 10 (1955) 759–762. [50] F.W.J. Olver, D.W. Lozier, R.F. Boisvert, C.W. Clark (Eds.), NIST Handbook of Mathematical Functions, Cambridge University Press, New York, 2010. [51] R. Wong, Y.Q. Zhao, On a uniform treatment of Darboux’s method, Constr. Approx. 21 (2005) 225–255. [52] E. Hille, Some aspects of the Thomas–Fermi equation, J. Anal. Math. 23 (1970) 147–170. [53] J.P. Boyd, The Blasius function in the complex plane, J. Experiment. Math. 8 (1999) 381–394. [54] V.B. Mandelzweig, F. Tabakin, Quasilinearization approach to nonlinear problems in physics with application to nonlinear ODEs, Comput. Phys. Comm. 141 (2001) 268–281. [55] N.A. Zaitsev, I.V. Matyushkin, D.V. Shamonov, Numerical solution of the Thomas–Fermi equation for the centrally symmetric atom, Russian Microelectron. 33 (2004) 303–309. [56] V. Bush, The differential analyzer. A new machine for solving differential equations, J. Franklin Inst. 212 (1931) 447–488. [57] C. Miranda, Teoremi e metodi per l’integrazione numerica della equazione differenziale di Fermi, Memorie della Reale Accademia d’Italia, Classe di scienze fisiche, Mat. Nat. 5 (1934) 285–322. [58] J.P. Boyd, The asymptotic Chebyshev coefficients for functions with logarithmic endpoint singularities, Appl. Math. Comput. 29 (1989) 49–67.