- Email: [email protected]
Chebyshev series solution of the Thomas-Fermi equation
Computer Physics Communications
Computer Physics Communications 67 (1992) 389—391 North-Holland
Chebyshev series solution of the Thomas—Fermi equation Allan J. MacLeod Department of Mathematics and Statistics, Paisley College of Technology, High Street, Paisley, Scotland, UK Received 15 March 1991
The numerical solution of the Thomas—Fermi equation is considered. Chebyshev series for small and large values of x are derived. Values of the coefficients to 1OD are given.
1. Introduction
while for x large, the following is the asymptotic series,
The Thomas—Fermi equation is defined by d2y
—
p3/2
with the conditions, y
y=144x3(1+diz+d2z2+ (1)
=
1 at x
=
0 and y
=
(3)
•..),
where z= ~
0 at with
x00.
Numerous tables of values of y exist, but there do not seem to be many approximations that can be used for any value of x. Mason [1] gives coefficients for a rational approximation with a claimed accuracy of 4 decimal places. Krutter [2] gives coefficients for a power series and for an asymptotic series, but the coefficients decline in size very slowly. Anderson and Arthurs [3], and Burrows and Core [4], derive approximations from a variational method but these do not seem to be very accurate. In this note, Chebyshev series coefficients for y are derived which can be used to calculate y to an accuracy of up to 10 decimal places.
fl
=
(v~i 7)/2 —
=
0.772001 8726587656.
This unusual combination of behaviour suggests that two differing approximations are used, one for small x, one for large x. Accordingly we have ~‘g(l), y (x)
)
0
X
=
2(x/L)~2
t
=
2( L/x)
—
1.0,
L,
144x 3/i ( t), L < x < 00, -
t
—
1.0, (4)
where g( t) and h (t) are Chebyshev expansions, 2. Approximations
g(t)=a 0T0(t)+a1T1(e)+ h(t)=b0T0(t)+b1T1(t)±
For x small, y is known to have the following power series y = 1 + Bx + c 3’~’2+ C 2 + C 5”2 + ..., (2) 3x 4X 5X OO10-4655/92/$05.O0 © 1992
—
+aNTN(t), +bNTN(t)
()
1t) being the ith Chebyshev with T,(t) = cos(i polynomial. Usingcos simple analysis it is easy to
Elsevier Science Publishers B.V. All rights reserved
390
A.J. MacLeod
/ Chebyshev series solution of the Thomas
show that g and h satisfy the following non-linear equations:
—
Fermi equation
Table 2 Chebyshev coefficients
_____________________________________________ a-
(1 ~
—
n2(1 =
+
t)2~
dt
dg +
—
L3/2(
2g3/2=0, 1 +t)
(n2 + 7n)(1 +
dh
t)~
+
12h
(6)
12h3”2.
The boundary conditions on y give the following conditions: g= 1 when t —1, and h 1 when —1. Consideration of y’ at x 0 also shows that ~ 0 when t —1. Continuity of y and at x L is also preserved giving the following two conditions: =
=
=
=
=
=
~‘
=
g(1)
=
144 _L3_h(l)~ (7)
-144 L3 (2nh(1)+3h(1)).
0
0.5090350232
E0
0.3239038967
E0
I
—0.4822560319
E0
—0.4048513212
E0
2 3
0.65697 6538 0.31187 5770
E—1 E—1
0.17761 12639 —0.64425 3632
E0 E—l
4
—0.200088032
E—1
E—2 E—3 E—3 E—4 E—4 E—5 E—6 E—6
0.207349604
—0.61454934 0.17147744 —0.45690 65 0.11739 00 —0.292836 0.71290 —0.17003 0.3985
E—1
12
0.54185 416 —0.6020586 —0.1305680 0.783231 —0.17920 1 0.14460 0.4860 —0.2244
E—2 E—2 E—3 E—3 E—4 E—5 E—5 E—6
13 14
0.451 —0.30
E—7 E—8
—0.920 0.210
E—7 E—7
E—8 E—9
—0.47 0.11 —0.2
E—7 E—8 E—9
5 6 7 8 9
10 11
15 —0.12 16 0.5 17 ____________________
where K 3. Solution of differential equations
The two non-linear equations are solved by using an iterative Newton method described by Norton [5]. This method leads to the following iteration equations, 2)(1 + t)2gr+ (1 + t)~r±i ~r+i (3Kg~” 1 —
2g~2, =
n2(1 +
—K(1
+ t) + t)2)~r+i+
(12
—
(n2 + 7n)(1 + t)/;r±i 18hy2 ) h,.±~ 6h3~2 =
(8)
—
r
‘
Table I Values of m1 and m2 for various L L mi 1 12 2 15 3 16
~L3~2. These equations are now linear
for g,.~ At
—
=
m2
1and hr± ~ each iteration
the coefficients of the
Chebyshev expansions are determined by (1) satisfaction of the five conditions from section 2, (2) collocation at selected points. If we define t, cos(i’rr/N), i 0 N, then we use t1 t~2 as collocation points for ~ and 1~,...,tN_i as collocation points for h,~1. This gives a total of 2N + 2 linear equations which are solved by Gaussian elimination with partial pivoting. Iteration is continued until successive iterates change by less than some specified tolerance. As initial estimates g 0 and h0, the rational =
=
approximation of Mason is used to generate val1N• This gives (N + linear systems for the ues oftwo g and h 1) at xthe(NN ++1)1 points 10 initial Chebyshev coefficients.
28 23 20
4
17
18
5
18
17
6
19
16
7 19 15 8_____________________________________________ 20 15
4. Numerical results All the calculations were performed on a
PRIME 6350 using the F77 compiler which has the facility for quadruple precision arithmetic. This
A.J. MacLeod
/ Chebyshev series solution ofthe Thomas
gives a machine precision of about 10 30• The first stage was to determine a suitable break-point L. This was done using N = 30, and iteration continued until successive iterates remained the same to within the rounding error. If we define m1 and m2 as follows, a,
10_b,
I
m~,
lb1 I
10~°,I
m~,
then table 1 contains values of m1 and m2 for various values of L. The table shows that a reasonable choice for L would be in the 4—5 range. The value L = 4 was chosen since this gives a simple t transform on [0, L], namely t = 1. With this value of L the calculations were repeated using N = 50 and the coefficients to 1OD are given in table 2. Finally, as a by-product of the calculation, the values of B in eq. (2) and d1 in eq. (3) can be determined. These are important constants in the theory of the Thomas—Fermi equation. The
—
Fermi equation
391
quadruple-precision results for L = 4 give the values B = —1.5880710226 and d1
=
—
13.2709738480
to 1OD. The value of B agrees with that found by Krutter while d1 agrees with the 6D value found by Kobayashi et al. [6].
References
—
[1] J.C.
Mason, Proc. Phys. Soc. 84 (1964) 357. [2] H. Krutter, J. Comput. Phys. 47 (1982) 308. [3] N. Anderson and A.M. Arthurs, Q. Appl. Math.
39 (1981)
127. [4] B.L. Burrow and P.W. Core, Q. AppI. Math. 42 (1984) 73. [5] H.J. Norton, Comput. J. 7 (1964) 76. [6] S. Kobayashi, T. Matsukuma, S. Nagai and K. Umeda, J. Phys. Soc. Japan 10 (1955) 759.
Computer Physics Communications 67 (1992) 389—391 North-Holland
Chebyshev series solution of the Thomas—Fermi equation Allan J. MacLeod Department of Mathematics and Statistics, Paisley College of Technology, High Street, Paisley, Scotland, UK Received 15 March 1991
The numerical solution of the Thomas—Fermi equation is considered. Chebyshev series for small and large values of x are derived. Values of the coefficients to 1OD are given.
1. Introduction
while for x large, the following is the asymptotic series,
The Thomas—Fermi equation is defined by d2y
—
p3/2
with the conditions, y
y=144x3(1+diz+d2z2+ (1)
=
1 at x
=
0 and y
=
(3)
•..),
where z= ~
0 at with
x00.
Numerous tables of values of y exist, but there do not seem to be many approximations that can be used for any value of x. Mason [1] gives coefficients for a rational approximation with a claimed accuracy of 4 decimal places. Krutter [2] gives coefficients for a power series and for an asymptotic series, but the coefficients decline in size very slowly. Anderson and Arthurs [3], and Burrows and Core [4], derive approximations from a variational method but these do not seem to be very accurate. In this note, Chebyshev series coefficients for y are derived which can be used to calculate y to an accuracy of up to 10 decimal places.
fl
=
(v~i 7)/2 —
=
0.772001 8726587656.
This unusual combination of behaviour suggests that two differing approximations are used, one for small x, one for large x. Accordingly we have ~‘g(l), y (x)
)
0
X
=
2(x/L)~2
t
=
2( L/x)
—
1.0,
L,
144x 3/i ( t), L < x < 00, -
t
—
1.0, (4)
where g( t) and h (t) are Chebyshev expansions, 2. Approximations
g(t)=a 0T0(t)+a1T1(e)+ h(t)=b0T0(t)+b1T1(t)±
For x small, y is known to have the following power series y = 1 + Bx + c 3’~’2+ C 2 + C 5”2 + ..., (2) 3x 4X 5X OO10-4655/92/$05.O0 © 1992
—
+aNTN(t), +bNTN(t)
()
1t) being the ith Chebyshev with T,(t) = cos(i polynomial. Usingcos simple analysis it is easy to
Elsevier Science Publishers B.V. All rights reserved
390
A.J. MacLeod
/ Chebyshev series solution of the Thomas
show that g and h satisfy the following non-linear equations:
—
Fermi equation
Table 2 Chebyshev coefficients
_____________________________________________ a-
(1 ~
—
n2(1 =
+
t)2~
dt
dg +
—
L3/2(
2g3/2=0, 1 +t)
(n2 + 7n)(1 +
dh
t)~
+
12h
(6)
12h3”2.
The boundary conditions on y give the following conditions: g= 1 when t —1, and h 1 when —1. Consideration of y’ at x 0 also shows that ~ 0 when t —1. Continuity of y and at x L is also preserved giving the following two conditions: =
=
=
=
=
=
~‘
=
g(1)
=
144 _L3_h(l)~ (7)
-144 L3 (2nh(1)+3h(1)).
0
0.5090350232
E0
0.3239038967
E0
I
—0.4822560319
E0
—0.4048513212
E0
2 3
0.65697 6538 0.31187 5770
E—1 E—1
0.17761 12639 —0.64425 3632
E0 E—l
4
—0.200088032
E—1
E—2 E—3 E—3 E—4 E—4 E—5 E—6 E—6
0.207349604
—0.61454934 0.17147744 —0.45690 65 0.11739 00 —0.292836 0.71290 —0.17003 0.3985
E—1
12
0.54185 416 —0.6020586 —0.1305680 0.783231 —0.17920 1 0.14460 0.4860 —0.2244
E—2 E—2 E—3 E—3 E—4 E—5 E—5 E—6
13 14
0.451 —0.30
E—7 E—8
—0.920 0.210
E—7 E—7
E—8 E—9
—0.47 0.11 —0.2
E—7 E—8 E—9
5 6 7 8 9
10 11
15 —0.12 16 0.5 17 ____________________
where K 3. Solution of differential equations
The two non-linear equations are solved by using an iterative Newton method described by Norton [5]. This method leads to the following iteration equations, 2)(1 + t)2gr+ (1 + t)~r±i ~r+i (3Kg~” 1 —
2g~2, =
n2(1 +
—K(1
+ t) + t)2)~r+i+
(12
—
(n2 + 7n)(1 + t)/;r±i 18hy2 ) h,.±~ 6h3~2 =
(8)
—
r
‘
Table I Values of m1 and m2 for various L L mi 1 12 2 15 3 16
~L3~2. These equations are now linear
for g,.~ At
—
=
m2
1and hr± ~ each iteration
the coefficients of the
Chebyshev expansions are determined by (1) satisfaction of the five conditions from section 2, (2) collocation at selected points. If we define t, cos(i’rr/N), i 0 N, then we use t1 t~2 as collocation points for ~ and 1~,...,tN_i as collocation points for h,~1. This gives a total of 2N + 2 linear equations which are solved by Gaussian elimination with partial pivoting. Iteration is continued until successive iterates change by less than some specified tolerance. As initial estimates g 0 and h0, the rational =
=
approximation of Mason is used to generate val1N• This gives (N + linear systems for the ues oftwo g and h 1) at xthe(NN ++1)1 points 10 initial Chebyshev coefficients.
28 23 20
4
17
18
5
18
17
6
19
16
7 19 15 8_____________________________________________ 20 15
4. Numerical results All the calculations were performed on a
PRIME 6350 using the F77 compiler which has the facility for quadruple precision arithmetic. This
A.J. MacLeod
/ Chebyshev series solution ofthe Thomas
gives a machine precision of about 10 30• The first stage was to determine a suitable break-point L. This was done using N = 30, and iteration continued until successive iterates remained the same to within the rounding error. If we define m1 and m2 as follows, a,
10_b,
I
m~,
lb1 I
10~°,I
m~,
then table 1 contains values of m1 and m2 for various values of L. The table shows that a reasonable choice for L would be in the 4—5 range. The value L = 4 was chosen since this gives a simple t transform on [0, L], namely t = 1. With this value of L the calculations were repeated using N = 50 and the coefficients to 1OD are given in table 2. Finally, as a by-product of the calculation, the values of B in eq. (2) and d1 in eq. (3) can be determined. These are important constants in the theory of the Thomas—Fermi equation. The
—
Fermi equation
391
quadruple-precision results for L = 4 give the values B = —1.5880710226 and d1
=
—
13.2709738480
to 1OD. The value of B agrees with that found by Krutter while d1 agrees with the 6D value found by Kobayashi et al. [6].
References
—
[1] J.C.
Mason, Proc. Phys. Soc. 84 (1964) 357. [2] H. Krutter, J. Comput. Phys. 47 (1982) 308. [3] N. Anderson and A.M. Arthurs, Q. Appl. Math.
39 (1981)
127. [4] B.L. Burrow and P.W. Core, Q. AppI. Math. 42 (1984) 73. [5] H.J. Norton, Comput. J. 7 (1964) 76. [6] S. Kobayashi, T. Matsukuma, S. Nagai and K. Umeda, J. Phys. Soc. Japan 10 (1955) 759.