# Spectral Method for Solving the Nonlinear Thomas-Fermi Equation Based on Exponential Functions

Authors: Raka Jovanovic; Sabre Kais; Fahhad H. Alharbi
Journal: Journal of Applied Mathematics
 (2014)

---

## Abstract

We present an efficient spectral methods solver for the Thomas-Fermi equation for neutral atoms in a semi-infinite domain. The ordinary differential equation has been solved by applying a spectral method using an exponential basis set. One of the main advantages of this approach, when compared to other relevant applications of spectral methods, is that the underlying integrals can be solved analytically and numerical integration can be avoided. The nonlinear algebraic system of equations that is derived using this method is solved using a minimization approach. The presented method has shown robustness in the sense that it can find high precision solution for a wide range of parameters that define the basis set. In our test, we show that the new approach can achieve a very high rate of convergence using a small number of bases elements. We also present a comparison of recently published results for this problem using spectral methods based on several different basis sets. The comparison shows that our method is highly competitive and in many aspects outperforms the previous work.


---

## Body

## 1. Introduction

In this paper we focus on solving the Thomas-Fermi (TF) equation, which is of great importance for a wide range of physical problems. Some of its applications are the determination of effective nuclear charge in heavy atoms and in finding effective potentials for self-consistent calculations. The equation is a nonlinear ordinary differential equation that is solved on a semi-infinite domain. The scaled TF equation [1] is given in
(1)d2ydx2=1xy3/2
subjected to the following boundary conditions
(2)y0=1,lim⁡x→∞yx=0.Due to its importance this equation has been solved by many different methods like the perturbative approach [2], homotopy analysis method [3, 4], quasilinearization approaches [5], and Padé approximations [6]. The complexity of solving this relatively simple looking equation is that it is singular at both endpoints.One of the complexities in applying spectral methods (SM) to the TF equation is the fact that it is defined on a semi-infinite domain. Significant research has been conducted on applying SM on infinite and semi-infinite domains [7–15]. This has been achieved by implementing a wide range of approaches varying from using suitable basis sets and truncating the numerical window to forcing size scaling. Very good results for such problems have been achieved by using nonclassical orthogonal basis sets for systems [9], mapped orthogonal systems [16, 17], Laguerre functions [10], mapped Legendre functions [11], and mapped Fourier sine series [12].Recently several approaches have been developed for solving the TF equation using pseudo spectral methods. The greatest effort in applying this type of approach has been done by using different versions of rational Chebyshev functions. In the work of Parand and Shahini [18] the original equation has been solved by using this basis set. In his article, Boyd [19] presents a spectral method based on the same basis set but applying it to the transformed version of the equation as given in [20] to get solutions of a much higher accuracy. His method has been able to use a very high number of basis functions while avoiding numerical instability and achieving precision of 10-25. The TF equation has also been solved by using the second [21] and third [22] kind of rational Chebyshev functions by Kilicman. These two approaches manage to get a significantly faster rate of convergence than the mentioned work of Boyd but have a significant level of parameter tuning that diminishes the robustness of the approach. The equation of interest has also been solved using Hermite [23] and Sinc functions [24]. In all of the mentioned applications of spectral methods, the collocation approach has been used.In our approach to solve the TF equation, we apply the spectral method based on the exponential basis set. This basis set and its polynomial version have been recently used for solving several differential equations on semi-infinite domains [13–15]. The use of a similar basis set was initially presented in the 1970s by Raffenetti, Bardo, and Ruedenberg [25–27] for self-consistent field wave functions. It is important to mention that the use of such a basis set is closely connected to Prony analysis. The inspiration for applying such a basis set for acquiring a numerical approximation for the solution of the TF equation comes from the fact that a combination of exponential functions was used for finding approximate analytic solutions for this problem. More precisely, such solutions were calculated by the use of the variational principle in [28, 29]. An improvement with similar but not strictly exponential functions has recently been presented by Oulne [30, 31]. Research has also been conducted on finding the analytic approximation to the solution by using the homotopy analysis method combined with a polynomial exponential basis [32].Similar to the mentioned application of SM to the TF equation, when using the exponential basis set the nonlinear ODE is converted to a set of algebraic equations. One of the main advantages of such an approach is that there is no need for numerical integration since all the underlining integrals can be analytically solved. In our implementation we apply the exponential basis set to the transformed version of the equation as presented in [20]. Another positive aspect of using such a basis set is the direct analytical calculation of the first derivative of the solution of TF equation. The importance of the value of the derivative at point 0 (y′(0)) is due to its major role in determining the energy of a neutral atom in the TF equation. Although the transformed equation has been solved, the derivative of the original equation can be acquired analytically from it and in this way the creation of additional error has been avoided. In our numerical experiments we show that the new approach is highly competitive with existing methods and in many aspects outperforms the previous work.The paper is organized as follows. In the second section we present the exponential basis set. In the next section we show details of applying this basis to the TF equation. In the fourth section we give a comparison of the new method to previously published research on the application of spectral methods for the problem of interest.


## 2. The Exponential Basis Set

It has been shown that the exact particular solution of (1), without the boundary condition at 0, is the function 144/x3 [33]. This function also gives us the asymptotic behavior of the solution of TF equation at infinity. The solution of the TF equation has been effectively approximated using variational approaches based on a combination of exponential functions. Although an approximate solution having exponential decay will tend to zero faster than the exact solution it has been shown that such approximations can give a high level of accuracy even for relatively large values of x [28–32]. This gives us incentive to represent its approximation using an exponential basis set. The goal of the SM approach is to find the values of coefficients ai∈R that best satisfy the following equation:
(3)ft≈∑i=1N‍aiRi,Ri=e-βit.In (3), the values of βi∈R are selected in an intuitive way to cover all the possible decay rate. For the proposed basis set to converge at infinity it is also necessary that the values of βi are strictly positive. When using such a basis set, it is possible to directly specify the values for βi depending on the problem specifics, or use some more general method. The later is preferential since it makes the basis more robust and also avoids unnecessary “fine tuning” when a new problem is solved. We use a slight modification of the method for generating values for βi presented in [13], which is given in
(4)βi=spi(5)pi=ds+i-1N-1dl-ds.In (4)  s is used to define the level of exponential decay of the different values of βi. In the same equation pi is used to give a distribution for the values of each βi. This distribution is defined in (5), where N is the number of basis functions and ds and dl represent the smallest and largest values for pi. All the values for pi are uniformly distributed in the interval [ds,dl]. We wish to point out that when applying this basis set to a specific problem some simple rough tuning is necessary in the sense of specifying the values for these three parameters. In the last part of the paper it will be shown that there is a wide range of their values in which the convergence rate is high. It is important to mention that this is a nonorthonormal basis set, but it can be easily converted to an orthonormal one using the Gram-Schmid orthogonalization procedure.The proposed use of the exponential basis set is closely related to Prony’s analysis which is used for fitting a sum of exponential functions to equally spaced data points and extended the model to interpolate at intermediate points. In case of Prony’s method a function is approximated in the same form as (3) but βi,ai∈C. The method itself gives a very efficient algorithm for finding the values of βi,ai if a set of equally spaced data points and corresponding function values are known. The original algorithm and the modified version [34] have been successfully applied in a wide range of problems [35–37]. Although Prony’s method has extensively been used in the field of signal processing only limited research has been conducted for its application to solving ODEs. Osborne has presented a method for applying a modified version of this method to linear homogeneous difference equation [38]. Maergoiz has done research for inhomogeneous linear ordinary differential equations using the modified method in case of one [39] or more dimensions [40].The proposed use of the exponential basis set can be understood as a specialized version that can be applied to ODEs which have solutions without oscillating decay behavior at infinity. This fact makes the use of a real valued exponent basis set significantly different in several aspects. While the use of only real coefficients is less robust and cannot easily be applied to multimodal functions, it is still suitable for the Thomas-Fermi equation. On the other hand due to it simpler form it has some specific properties. More precisely it can be used as a real spectral method and can efficiently be applied to nonlinear problems. We wish to emphasize that the proposed method does not use the main aspect of Prony’s algorithm, the computationally efficient method for calculating the parameters that specify the expansion from values at equidistant points. But due to the similarity of the basis sets used, there exists a strong potential for adapting Prony’s method to make the proposed one more computationally efficient.It is important to point out that, unlike the standard spectral basis sets, no rigorous convergence theory which proves that an arbitrarily large accuracy can be achieved using a high enough number of basis functions has not been developed, to the best of our knowledge. On the other hand the exponential basis set and its polynomial version have been successfully applied to a wide range of problems [13–15, 25–27, 32], which gives incentive to attempt developing such theory.


## 3. Solving the Thomas-Fermi Equation

The direct solution of the equation of interest given by (1) is possible but a preferred form, due to the simplified differential form where the radicals are avoided and hence the integrals are simplified, is the one given in [20]. The new form of the ODE is acquired by using the following variable substitutions:
(6)t=x1/2,ft=yt21/2.
After substituting the new variables into (1) we have a new nonlinear ODE given in
(7)A=tftf′′t+f′t2-ftf′t-2t2ft3=0.
The transformed equation (7) has the same boundary conditions as the original equation.The new ODE is transformed into a set of algebraic equations using the standard spectral method approach where each of the basis functionsRi is used as a test function. As our goal is to find the coefficients ai that best satisfy (7), it is also necessary to present the f′ and f′′ in the proposed basis set. This can easily be done by using the properties of the exponential function and their expansions are given in (8). A new algebraic system can easily be formed when
(8)f′t≈-∑i=1N‍aiβie-βit,f′′t≈∑i=1N‍aiβi2e-βit.
The new algebraic system of equations can be written in form given in (10):
(9)R=R1,R2,…,Rn⊤(10)∫0∞‍ARdt=M1+M2-M3-2M4=0.Equation (10) represents the simplified representation of all the equations that appear in the new system. With the intention of having a more comprehensible form of the equation we have used vector Mi to simplify the notation. Each element of vectors Mi corresponds to one of the test functions. More precisely, each of the vector elements for Mi,k are calculated in the following way:
(11)M1,k=∫0∞‍tftf′′te-βktdtM2,k=∫0∞‍tf′t2e-βktdtM3,k=∫0∞‍ftf′te-βktdtM4,k=∫0∞‍t2ft3e-βktdt.As previously mentioned one of the main advantages of the proposed basis set is that, after substitutingf,f′,f′′ with their expanded forms, given in (3), (8), all the integrals in (11) can be analytically solved. This is due to the fact that each of these integrals can be presented as a summation of integrals having the following form:
(12)∫0∞‍tke-atdt=k!ak+1.In (12)  k is a nonnegative integer. To be exact, after incorporating the expansions of f,f′,f′′ into (11) and using the analytic solutions of integrals given in (12) the values of Mi,k can be easily calculated from the summations given in the following equations:
(13)M1,k=∑i=1N‍∑j=1N‍aiajβj21βi+βj+βk2M2,k=∑i=1N‍∑j=1N‍aiajβiβj1βi+βj+βk2M3,k=-∑i=1N‍∑j=1N‍aiajβj1βi+βj+βkM4,k=∑i=1N‍∑j=1N‍∑m=1N‍aiajam2βi+βj+βm+βk3.In addition to theN equation, for each of the test functions that are acquired by incorporating the substitutions given in (13) into (10), it is necessary to include the boundary condition at zero given in
(14)∑i=1N‍ai=1.
In this work, the Tau approach is used where the first N-1 rows of the system are given in (10) and (14) is inserted in the Nth row.


## 4. Results

To evaluate our method, we have implemented the numerical algorithm using MATLAB R2013b. The nonlinear algebraic system of equations, presented in the previous section, was solved using the built-infsolve function. This method uses a minimization algorithm based on trust regions. In practicefsolve uses an iterative and self-correcting procedure to solve the system given in (10) with n variables ai corresponding to the number of bases, for which it needs an initial solution. To validate the method we have first used the initial values ai=1 to confirm that convergence is being achieved. It is important to point out that although for a wide range of initial values of ai the same behavior occurs, in some specific cases like when ai≫1 for most of i, the method does not manage to converge. The summations given in (13) which represent the algebraic system, if not given in vector/matrix form, can be very computationally expensive if the method is implemented using MATLAB. All the mentioned equations can be efficiently be converted in such a vector/matrix form, which has been done in our implementation of the proposed method.We have evaluated the quality of the approximation function acquired by the use of an exponential basis set in two ways. First, the value of the initial slopey′(0) is computed due to its important role in determining many physical properties. Specifically, it is essential in the calculation of the the energy of a neutral atom in the TF approximation [1]:
(15)E=674π32/3Z7/3y′0,
where Z is the nuclear charge.It is important to point out that the chosen basis set is specially suitable for this problem due to the fact that the derivative ofy′ in (1) can directly be calculated using the values of coefficients ai calculated for the transformed equation. This can done using some simple derivation. Starting from
(16)y′x=∫‍y′′xdx=∫‍1xy3/2dx,
after incorporating the variable substitutions given in (6) and using the proposed expansion, we have
(17)y′x=∫‍1xfx23/2dx=-2∫‍∑i=1N‍aie-βit3dt.
Finally we have the formula for calculating the values for y′(x) as follows:
(18)y′x=-2∑i=1N‍∑j=1N‍∑k=1N‍aiajake-(βi+βj+βk)xβi+βj+βk.
In this way we avoid adding any additional numerical error at this step. The use of an exponential basis set will give a similar advantage when solving similar problems.In Table1, we give a comparison of our method (Exp) to recently published results on the application of spectral methods for solving the TF equations. More precisely we compare it to the use of several different basis sets: rational Chebyshev functions in the implementation of Parand (RC-P) [18] and Boyd (RC-B)[19], rational Chebyshev functions of the second (RC-SK) and third kind (RC-TK), Hermite functions (Her) [23], and Sinc functions (Sinc) [24]. We compare the error to the best known approximate solution given by Boyd: y′(0)=1.5880710226113753127186845 [19]. In the case of our method we have used the following parameter values, s=5,ds=-4,dl=2, to specify the basis set. These values have been chosen due to the fact that they have produced the approximation of the highest quality. We wish to emphasize that there is a wide range of values of these parameters that produce a similar level of precision, so this choice cannot be considered as fine tuning. In the last part we will show that the proposed method is not highly sensitive to the choice of these parameters and can be considered robust. In the goal of having a better analysis of the effect of using different basis sets we have presented the values of the function y in Table 2 for a wide range of values of the variable x. These values have also been collected from the same published works.Table 1
Comparison of the effectiveness of the application of spectral methods based on different basis sets for calculating the initial slopey′(0) of the TF equation.
N
RC-P
RC-SK
RC-TK
RC-B
Her
Sinc
Exp
4
—
—
—
—
3.57
e
-
02
9.89
e
-
02
5
—
8.11
e
-
07
1.59
e
-
09
—
—
—
2.76
e
-
01
7
—
3.24
e
-
07
3.46
e
-
11
—
—
—
1.70
e
-
02
8
2.34
e
-
06
—
1.33
e
-
15
—
—
6.64
e
-
03
3.85
e
-
03
10
9.91
e
-
07
—
—
—
—
—
6.48
e
-
04
12
7.26
e
-
07
—
—
—
—
—
5.13
e
-
05
15
—
—
—
—
<5e-07
—
9.03
e
-
06
16
—
—
—
—
—
6.75
e
-
04
6.03
e
-
07
20
—
—
—
5.26
e
-
6
—
—
6.09
e
-
08
25
—
—
—
—
—
—
2.00
e
-
10
32
—
—
—
—
—
6.84
e
-
07
4.13
e
-
09
40
—
—
—
3.90
e
-
11
—
—
3.88
e
-
09Table 2
Comparison of the approximate solutionsy of the TF equation acquired by spectral methods based on different basis sets.
x
RC-P
RC-SK
RC-TK
RC-B
Her
Sinc
Exp
0.25
0.755880759
0.755903399
0.755455402
—
0.754795330
—
0.755202096
0.50
0.606700008
0.605270502
0.602998554
—
0.606658908
—
0.606986951
0.75
0.502964042
0.497823942
0.494347872
—
0.502110510
—
0.502348140
1.00
0.424333179
0.420343948
0.416399658
—
0.423811203
0.4240642728
0.424010148
1.25
0.363227937
0.362814756
0.358770806
—
0.363027725
—
0.363203991
1.50
0.314660642
0.318737461
0.314761643
—
—
—
0.314780118
1.75
0.275233848
0.284014958
0.280179962
—
—
—
0.275453712
2.00
0.242678587
0.256010764
0.252344355
—
0.242918233
0.2430282344
0.243010373
2.25
0.215439334
0.232974191
0.229482688
—
0.215819818
—
0.215895823
2.50
0.192406328
0.213705386
0.210384924
—
0.192917948
—
0.192984580
2.75
0.172758691
0.197357837
0.194199930
—
0.173379623
—
0.173440997
3.00
0.155871862
0.183318729
0.180313058
—
0.156573773
—
0.156631657
3.25
0.141260504
0.171134229
0.168270054
—
0.142013368
—
0.142067963
3.50
0.128541381
0.160461449
0.157728304
—
0.129316613
—
0.129367328
3.75
0.117408054
0.151036695
0.148424721
—
0.118180209
—
0.118226225
4.00
0.107612958
0.142653971
0.140154047
—
0.108360441
0.1084800242
0.108401057
4.25
0.098954329
0.135150082
0.132753853
—
—
—
0.099694306
4.50
0.091266456
0.128394110
0.12609396
—
—
—
0.091944333
4.75
0.084412289
0.122279834
0.120068868
—
—
—
0.085017755
5.00
0.078277758
0.116720187
0.114592127
—
—
—
0.078803669
6.00
0.059236422
0.098752518
0.096904158
—
—
0.0594010318
0.059418888
7.00
0.046231022
0.085573204
0.083941323
—
—
—
0.046094375
8.00
0.036980654
0.075494762
0.074034822
—
0.036580427
0.0365353166
0.036584707
9.00
0.030180901
0.067538722
0.066218399
—
—
—
0.029589375
10.00
0.025044744
0.061098873
0.059894055
0.024314292
—
0.0243069553
0.024313708
15.00
0.011731220
0.041370727
0.040533524
—
0.010803774
—
0.010808302
20.00
0.006585633
0.031271686
0.030630632
0.005784941
0.005792831
—
0.005789307
25.00
0.004104406
0.025135426
0.024616163
—
—
—
0.003478434
30.00
—
—
—
0.002255836
0.002252634
—
0.002260351
50.00
0.000761317
0.012687078
0.012420906
0.000632254
—
0.0006213281
0.000632255
75.00
0.000193597
0.008484835
0.008305908
—
—
—
0.000219970
100.00
0.000023409
0.006373709
0.006238954
0.000100242
—
—
0.000101341
500.00
—
—
—
0.000001034
—
—
0.000001068When we observe the results in Table1, it is noticeable that the methods RC-P, RC-SK, and RC-TK have a much higher speed of convergence for y′(0) than the other methods. We believe that this high precision of y′(0), for a small number of basis sets, does not correctly represent the actual effectiveness of these methods. If we observe the values for function y given in Table 2, or in a graphical representation in Figure 1, we can see that the precision of the approximations acquired by using these basis sets is much lower. This is most notable for approximate solutions RC-SK, RC-TK, which have been developed by the same authors, for which the values y(x) strongly deviate from the results of the other methods. In case of RC-P a similar deviation exists but to a much lower extent. We believe this is a consequence of fine tuning of the method. More precisely, in all of these papers the parameter L that is utilized for specifying the rational Chebyshev functions which are used as a basis set has been presented with precision of 10-6 and higher.Figure 1
Graphical representation of results given in Table2.The approximations for the functiony acquired with RC-B, Her, Sinc, and our method on the other hand have a very similar behavior. We believe that the results presented by Boyd are of the highest accuracy. It can be seen from the values in Table 2 that our method, with 25 basis functions, gives the result that is the closest to the one acquired by Boyd [41] who used 600 basis functions. Even for x=500 the relative difference between these two approaches is around 3% for values close to 10-6.In the case of observing the the speed of convergence for calculatingy′(0), our method proves to be very efficient. Figure 2 gives us a graphical representation of convergence rates of RC-B, Sinc, and Exp, which shows that the new method has a higher speed of convergence with the increase of the number of basis functions. In this figure the results for the other methods based on rational Chebyshev functions have been excluded due to the reason explained in the previous paragraphs. The maximal precision we have achieved with our approach is 2e-10 with 25 basis functions. When a higher number of basis functions has been used our method has started to exhibit numeric instability and the precision of the method has decreased. We assume that the cause for this is due to our implementation of the method using MATLAB where variables are stored using the 64 bit IEEE Standard 754 floating point numbers, for which the precision corresponds to the use of 16 significant decimal digits. In our implementation the calculation of y′ is done by summing 105 elements as presented in (18). This results in the problem of error accumulation. It is expected that a higher accuracy can be achieved by developing the method using software that allows multiprecision arithemtic. We show a more detailed behavior of the error for this approximation in Figure 3.Figure 2
Comparison of convergence speeds for different methods, corresponding to the number of basis functions (N) used, for calculating the initial slope y′(0). The exponential basis set is defined with the following parameters: s=5, ds=-4, and dl=2.Figure 3
Detailed presentation of the convergence rate of the initial slopey′(0) for an exponential basis set defined with the following parameters: s=5, ds=-4, and dl=2.Due to the fact that the proposed method starts exhibiting numerical instability for a relatively low number of basis functionsN, the true convergence rate cannot be truly recognized form Figure 3. Although that, from observing this figure, we get the impression that the method has an exponential convergence rate, from our previous work with the exponential basis [14] and analysis of convergence rates for series approximations for functions with singularities [42], we expect that the convergence speed is lower and closer to root exponential.To show the robustness of this approach, in Figure4 we present the speed of convergence for several different basis sets. The basis sets are generated using various sets of βi, which are acquired for different values of ds,dl. It is notable that, for s=5, values of ds ranging from -5 to -1 and for dl from 1 to 4 good precision are achieved. From this figure it is notable that the method manages to get accuracy of below 10e-9 in many cases, which reassures us in assumption for the source of error.Figure 4
Convergence rate of the initial slopey′(0) of the Thomas-Fermi equation calculated using different parameters for defining the exponential basis set. In all the tests the exponent had the value s=5.


## 5. Conclusion

In this paper, a new approach for solving the Thomas-Fermi equation using spectral methods has been presented. We have shown that by using an exponential basis set a very high level of accuracy of the approximate solution can be achieved by employing a relatively small number of basis functions. One of the main advantages of this method is that there is no need for numerical integration, due to the positive properties of the selected basis set.In our tests we have shown that the use of this basis set has a higher speed of convergence compared to other methods of this type. The approximate solutions acquired by this approach have given a high precision of the value fory′(0). This has been accomplished while maintaining a high precision for the function y, which has not been achieved by several methods of this type. In our test we have also shown that the use of an exponential basis set for this problem is highly robust. More precisely, it achieves a high level of convergence for a wide range of values of parameters that define the basis set.In our work we did not focus on developing a suitable orthonormal version of the basis set, since the Thomas-Fermi equation is one-dimensional, and it was not a problem to solve the underlying system of algebraic equations. In the future we plan to work on adapting the basis set in this way. Another direction of our future research will be in attempting to incorporate some aspects of Prony’s algorithm to have a better method for selecting the parameters that specify the basis set.



---

*Source: 168568-2014-11-12.xml*