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

,


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 subjected to the following boundary conditions 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 semiinfinite domain.Significant research has been conducted on applying SM on infinite and semi-infinite domains [7][8][9][10][11][12][13][14][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 semiinfinite domains [13][14][15].The use of a similar basis set was initially presented in the 1970s by Raffenetti, Bardo, and Ruedenberg [25][26][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 (  (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.

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/ 3  [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  [28][29][30][31][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   ∈ R that best satisfy the following equation: In ( 3), the values of   ∈ 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   are strictly positive.When using such a basis set, it is possible to directly specify the values for   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   presented in [13], which is given in In (4)  is used to define the level of exponential decay of the different values of   .In the same equation   is used to give a distribution for the values of each   .This distribution is defined in (5), where  is the number of basis functions and   and   represent the smallest and largest values for   .All the values for   are uniformly distributed in the interval [  ,   ].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   ,   ∈ C. The method itself gives a very efficient algorithm for finding the values of   ,   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][36][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.

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: After substituting the new variables into (1) we have a new nonlinear ODE given in 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 functions   is used as a test function.As our goal is to find the coefficients   that best satisfy (7), it is also necessary to present the   and   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 The new algebraic system of equations can be written in form given in (10): 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 M  to simplify the notation.Each element of vectors M  corresponds to one of the test functions.More precisely, each of the vector elements for  , are calculated in the following way: As previously mentioned one of the main advantages of the proposed basis set is that, after substituting ,   ,   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: In ( 12)  is a nonnegative integer.To be exact, after incorporating the expansions of ,   ,   into (11) and using the analytic solutions of integrals given in ( 12) the values of  , can be easily calculated from the summations given in the following equations: In addition to the  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 In this work, the Tau approach is used where the first  − 1 rows of the system are given in ( 10) and ( 14) is inserted in the th row.

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-in fsolve function.This method uses a minimization algorithm based on trust regions.In practice fsolve uses an iterative and self-correcting procedure to solve the system given in (10) with  variables   corresponding to the number of bases, for which it needs an initial solution.
To validate the method we have first used the initial values   = 1 to confirm that convergence is being achieved.It is important to point out that although for a wide range of initial values of   the same behavior occurs, in some specific cases like when   ≫ 1 for most of , 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 slope   (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]: where  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 of   in (1) can directly be calculated using the values of coefficients   calculated for the transformed equation.This can done using some simple derivation.Starting from after incorporating the variable substitutions given in (6) and using the proposed expansion, we have  2.
Finally we have the formula for calculating the values for   () as follows: 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 Table 1, 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:   (0) = 1.5880710226113753127186845 [19].In the case of our method we have used the following parameter values,  = 5,   = −4,   = 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  in Table 2 for a wide range of values of the variable .These values have also been collected from the same published works.
When we observe the results in Table 1, it is noticeable that the methods RC-P, RC-SK, and RC-TK have a much higher speed of convergence for   (0) than the other methods.We believe that this high precision of   (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  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 () 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  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.The approximations for the function  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  = 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 calculating   (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 2 − 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   is done by summing 10 5 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. Due to the fact that the proposed method starts exhibiting numerical instability for a relatively low number of basis functions , 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 Figure 4 we present the speed of convergence for several different basis sets.The basis sets are generated using various sets of   , which are acquired for different values of   ,   .It is notable that, for  = 5, values of   ranging from −5 to −1 and for   from 1 to 4 good precision are achieved.From this figure it is notable that the method manages to get accuracy of below 10 − 9 in many cases, which reassures us in assumption for the source of error.

Conclusion
In this paper, a new approach for solving the Thomas-Fermi equation using spectral methods has been presented.We have Journal of Applied Mathematics 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 for   (0).This has been accomplished while maintaining a high precision for the function , 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.

Figure 1 :
Figure 1: Graphical representation of results given in Table2.

Figure 2 :
Figure 2: Comparison of convergence speeds for different methods, corresponding to the number of basis functions () used, for calculating the initial slope   (0).The exponential basis set is defined with the following parameters:  = 5,   = −4, and   = 2.

Figure 3 :Figure 4 :
Figure 3: Detailed presentation of the convergence rate of the initial slope   (0) for an exponential basis set defined with the following parameters:  = 5,   = −4, and   = 2.

Table 1 :
Comparison of the effectiveness of the application of spectral methods based on different basis sets for calculating the initial slope   (0) of the TF equation.

Table 2 :
Comparison of the approximate solutions  of the TF equation acquired by spectral methods based on different basis sets.