A new numerical approach to solve Thomas–Fermi model of an atom using bio-inspired heuristics integrated with sequential quadratic programming

In this study, a novel bio-inspired computing approach is developed to analyze the dynamics of nonlinear singular Thomas–Fermi equation (TFE) arising in potential and charge density models of an atom by exploiting the strength of finite difference scheme (FDS) for discretization and optimization through genetic algorithms (GAs) hybrid with sequential quadratic programming. The FDS procedures are used to transform the TFE differential equations into a system of nonlinear equations. A fitness function is constructed based on the residual error of constituent equations in the mean square sense and is formulated as the minimization problem. Optimization of parameters for the system is carried out with GAs, used as a tool for viable global search integrated with SQP algorithm for rapid refinement of the results. The design scheme is applied to solve TFE for five different scenarios by taking various step sizes and different input intervals. Comparison of the proposed results with the state of the art numerical and analytical solutions reveals that the worth of our scheme in terms of accuracy and convergence. The reliability and effectiveness of the proposed scheme are validated through consistently getting optimal values of statistical performance indices calculated for a sufficiently large number of independent runs to establish its significance.

For practical analysis of the dynamics of TFE, most of the studies are conducted for restricted inputs instead of unbounded domain.
After a profound study of the literature regarding Thomas-Fermi equation (TFE) it is observed that only deterministic solvers are applied to analyze its dynamics and no one yet applied the stochastic solvers. Recently, stochastic numerical techniques based on artificial intelligence techniques are effectively used to calculate the accurate solutions for initial and boundary value problems (BVPs) of differential equations involving both integer and fractional derivatives (Parand et al. 2013b;Arqub and Abo-Hammour 2014;Abo-Hammour et al. 2013;Abu Arqub et al. 2012. For instance, few potential application of stochastic solvers are solution of nonlinear Van-der-Pol oscillatory systems (Khan et al. 2015b), inverse Kinematics problem (Momani et al. 2016), problem arising in Electromagnetic theory (Khan et al. 2015a), nonlinear singular systems (Abo-Hammour et al. 2014a), plasma physics problems (Raja 2014a), nonlinear Navier-Stokes problems (Abo-Hammour et al. 2014b), nanotechnology problems involving carbon nanotubes (Raja et al. 2016a), magnetohydrodynamic problems ), fuel ignition model of combustion theory (Raja 2014b), fluid dynamics problem of thin film flow (Raja et al. 2016d), mathematical models of electrically conducting solids (Raja et al. 2016b), Schrödinger equation for the hydrogen atom (Caetano et al. 2011), nonlinear Jeffery-Hamel flow in the presence of high magnetic field , and strong nonlinear systems based on Painlevé, Bratu, Emden-Fowler, Riccati, Bagley-Torvik, Troesch's, Lane-Emden, Flierl-Petviashivili and pantograph models (see Mall and Chakraverty 2013, 2015Raja et al. , 2015bRaja et al. , c, 2016c and references therein). Authors motivated from these studies to carry out exploration and exploitation in the field on stochastic numerical solvers to solve governing model of nonlinear TFE based on discritization with finite difference scheme and trained through bio-inspired computing integrated with sequential quadratic programming procedures.
The most incorporated stochastic solvers are normally based on bio-inspired computational heuristics through Genetic Algorithms (GAs), a kind of effective global search methodology. GAs are used extensively to solve variety of the problems arising in various applications in physical sciences (Homayouni et al. 2014;Chiroma et al. 2014Chiroma et al. , 2015Toledo et al. 2014;Chiroma et al. 2016) which motivates the authors to exploit the strength of these techniques to study the dynamics of nonlinear singular TFE. The advantages of these methodologies are reflected through simplicity of the concept, ease in implementation processes, wider applicability, avoid divergence, stability, robustness and reliability, which make them impressive to be exploited for challenging models of mathematical physics like TFE. With the advent of modern computer architectures based on signal processing platform, an immense increase in computational power of the machines is achieved which gives a rebirth to population based meta-heuristic methodologies to be used for stiff problems of mathematical physics. The significance of the present research is a step forward in designing the machines learning algorithms for providing the solution of highly nonlinear and singular system for Thomas-Fermi model of an atom given in the form of boundary value problem of TFE for unbounded domain.
The rest of the paper is organized as follows: in "Methods" section, the proposed design methodology based on discretization of differential equation into system of difference equations by using a finite difference scheme is provided along with the optimization procedure for solving system of nonlinear equations; in "Numerical experimentation and results" section, the results of numerical simulations for different cases of TFE are presented in a number of graphs and numerical illustrations; conclusions are listed in the last section along with few suggested research directions.

Methods
Proposed methodology for Thomas-Fermi Eq. (1) is presented here that consists of two parts. First part is the formulation of optimization problem with the construction of overall individual residual error with the help of finite different schemes satisfying the constrained boundary conditions, while is the second part, a hybrid computing framework based on Genetic Algorithms (GAs) supported with Sequential Quadratic Programming (SQP) is exploited for minimization of the overall residual error function.

Discretization through finite difference scheme
The simplest and effective technique for solving differential equation is based on finite difference schemes which are used for approximation of derivative terms in the system by using difference quotients.
To obtain the approximate solution of the Thomas-Fermi equation on equally distributed mesh points in the finite interval t ∈ [0, T ], one can proceed by taking t i = ih, i = 0, 1, . . . , N for h = 1/N, hence the approximated equation is given for interior mesh points, t i , i = 1, 2, . . . , N − 1 as: while the boundary conditions are given as: Here F (y(t i )) = y(t i ) 3 t i , β 1 = 1, and β 2 = 0 for Thomas-Fermi equation.
The difference quotients approximation formulates based on 5 interior mesh points is used to closely approximate y ′′ (t i ), i = 1, 2, . . . , N − 1, by taking small step size h with the error on the order of 0 (h 3 ) and are written mathematically for forward Δ, central Π and backward ∇ differences operators, respectively, as: where forward Δ, central Π and backward ∇ differences are defined as: The finite difference schemes are used for solving the Thomas-Fermi differential equations by the procedure of discretization.

Fitness function formulation
Discretization procedure of differential equations converts the equation to the system of algebraic equations which are then solved by construction of fitness function based on individual residual errors of each equation. The necessary details for the construction of fitness function is given here.
The finite difference approximation formulae for y ′′ (t i ), i = 1, 2, . . . , N − 1, given in (5), (6) and (7) are used in (3) to transform the Thomas-Fermi equation as: � y(t i−2 ), y(t i+2 ) = − 1 12 Equations (11), (12) and (13) are the system of algebraic equations with N dependent variables, i.e., y(t 0 ), y(t 1 ), …,y(t N ). In order to formulate a fitness function the residual errors R err are defined as: The overall individual residual function O R is defined, similar to l 2 the norm of the residuals of all nodes and it is given mathematically as: Now the requirement is to minimize the fitness function O R for which the value of individual residual errors for each equation decreases. Consequently, the desired results or optimized solutions of the problem are achieved when O R approaches zero.

Learning methodology
Residual error function (17) is minimized through hybrid computing approach consisting of GAs integrated with SQP method, i.e., GA-SQP algorithms.
First real application of GAs has been given by Holland (Holland 1975) in early 70's of the last century and afterwards GAs is used as one of the premier derivative free solver for both constrained and non-constrained optimization problems. GAs belongs to the class of global search methods formulated through mathematical modeling of natural genetic mechanism. The standard operation of GAs is based on its reproduction operators, which are selection, crossover, and mutation. GAs are applied effectively as a good optimization mechanism in diverse fields such as electronics, optics, electromagnetism, controls, digital communication, robotics, astrophysics, chemical industry, materials, signal processing, nuclear power systems, bioinformatics, economics, and financial mathematics etc. (see Haupt and Haupt 2004;Kumar et al. 2010;Dasgupta and Michalewicz 2013;Grefenstette 2013, and references therein). Few recently reported potential applications of GAs are optimization in orbital maneuvers (dos Santos and da Silva Formiga 2015), overlapping community detection in complex networks (Yuxin et al. 2015), formulation of a public bicycle-sharing system (Askari and Bashiri 2015), preemptive identical parallel machines scheduling problem (Aalaei et al. 2015) and spacecraft guidance and control system (Shirazi and Mazinan 2015).
The GAs is implemented through built-in functions available in the MATLAB optimization toolbox and for effective optimization. GAs are hybridized with SQP for rapid local search. In this manner, a hybrid meta-heuristic optimization mechanism is designed for training of weights of ANNs based on GAs integrated with SQP to solve Thomas-Fermi equation. Detailed workflow of the proposed design scheme, in terms of the problem, modeling, optimization procedures and comparison, are shown in Fig. 1. Optimization of the fitness function given in Eq. (17) is carried out with the help of GA-SQP, which are implemented through built-in routines of GAs and FMINCON with algorithm SQP available in the Matlab optimization toolbox. The parameter settings applied for GAs and SQP are given in Table 1. These settings are made with care, after a lot of experimentation. A slight variation in these settings may result in pre-mature convergence of the algorithms.
The detailed description of the procedural steps of proposed GA-SQP algorithm is given as follows: Firstly, the initial chromosome for GAs are created with bounded real values randomly having genes or elements equal to the number of unknown variables in the residual error function. These chromosomes formulate an initial population P for the algorithm. Mathematically the population P based on chromosomes is given as:

Thomas Fermi Model of an Atom
Nonlinear, singular, second order differential equation

Initialize of population, Random Assignment, Bounds and Declarations
Setting of "gaoptimset" parameters

Best of GA as Initial weight, Bounds and Declarations
Best weights of GA-SQP

Yes No
No

Fitness Evalution
Setting of "optimset" parameters Yes Fig. 1 Graphical abstract of proposed methodology for solving Thomas-Fermi equation where M is the total number of chromosomes in the population, while each chromosome has M elements or genes which represents the discritization points of finite difference scheme. The rest of parameters are set for GA as listed in Table 1. In second step, the value of fitness O R is determined for each chromosome c of the population P using Eq. (17) and its constituent parts given in Eqs. (14)(15)(16). Each chromosome c of the population P with a minimum value of fitness O R is ranked high and vice versa. Ranking each chromosome of the population accordingly. Thirdly, algorithm stops its execution in case of fulfillments of the fitness limit; total number of generations/cycles executed, tolerance limits are attained such as function tolerance (TolFun) and nonlinear constraint tolerance (TolCon). If termination criteria satisfied just go SQP algorithm otherwise reproduced population using crossover, mutation and selection operations by invoking the built-in functions for these genetic operators, as listed in Table 1. Repeat the procedure accordingly. The rapid refinement of the results is carried out using SQP algorithm by using 'fmincon' routine with initial weights which are the global best solution of GAs. Initial declarations, setting and boundes for the algorithm are listed in Table 1. Determined the value of fitness O R , as given in the Eq. (17) for each updated weights vector by SQP procedure and terminate the cyclic updating of weights if execution of total number of iterations, fitness, tolerances limits are achieved and refined design parameter are obtained. Finally, store the final optimized weights of both GA and GA-SQP algorithms along with their fitness, time consumed, generations executed, and functions evaluated.

Numerical experimentation and results
The results of numerical experiments are presented in this section for solving Thomas-Fermi equation by taking five scenarios based on the size of input interval, while in each scenario three cases are taken with different values of the step size parameter. The scenarios with different cases are given as follows: (1) for this scenario along with related boundary condition is given as: Scenario 2 In this scenario, TFE for input span t ∈ [0, 5] with cases 1, 2 and 3 based on step size parameter h = 0.5, 0.25, and 0.1, respectively is taken and mathematically is given as: Scenario 3 TFE (1) for input span t ∈ [0, 25] with cases 1, 2 and 3 based on step size parameter h = 2.5, 1.25, and 0.5, respectively is taken in this scenario and it is written as: Scenario 4 Study the dynamics of Thomas-Fermi model for relatively larger input span t ∈ [0, 50] with cases 1, 2 and 3 based on step size parameter h = 5.0, 2.5, and 1.0, respectively. Mathematically the model Eq.
(1) for this scenario is given as: Scenario 5 Solution of TFE is analyzed for larger input span t ∈ [0, 100] with cases 1, 2 and 3 based on step size parameter h = 10.0, 5.0, and 2.0, respectively. Mathematically the model Eq.
(1) for this scenario is written as: Design methodology is applied to obtain the solution of Thomas-Fermi equation for all three cases of each scenario as per procedure given in the last section. Fitness function as given in Eq. (17) for inputs t ∈ [0, 1] with step size h = 0.1, 0.05, 0.02 i.e., N = 10, 20, 50, for cases 1, 2 and 3, are formulated, respectively as: (R err (i)) 2 + (R err (9)) 2 , Accordingly, fitness functions O R are developed for all three cases of scenarios 2-5. Optimization of fitness function O R for each case of all five scenarios are carried out with the help of a hybrid computing approach based on GA-SQP in case of 100 independent runs using parameter settings as listed in Tables 1. Optimization output plots for the fitness function O R by GAs in terms of fitness values, current best chromosome, fitness scaling, selection, average distance between the individuals, best, mean and worst scores are shown in Fig. 2 for the case 1 of scenario 1. Accordingly, for the same case optimization plots of SQP algorithms in terms of current best point, function counts, learning curves, constraints violation, step size parameter, and first-order optimality, are given in Fig. 3. Similarly the optimization outputs for other cases are determined and results obtained with one of the runs of GA-SQP algorithms for each case of scenarios 1 and 2, are presented in Fig. 4, along with the values of fitness plotted against 100 independent runs of the GA-SQP algorithm. The fitness values are plotted on a semi-log scale in order to observe the small variation in the results. Accordingly, results of proposed solutions for all three cases of 3, 4 and 5 scenarios are given in Fig. 5, while result of statistical analysis are plotted in Fig. 6 for each scenario. In case of Fig. 4a the plots of all three cases based on values of step sizes h = 0.1, 0.05 and 0.02 are consistently overlapping. In case of Fig. 4b, it is seen that with the decrease in step size, the value of fitness also decreases, which is due to the fact that with the decreased step size, the discretization of the system using the finite difference scheme based on increased number of nonlinear equations. Consequently, the system becomes stiff with a decrease in step size and hence the solution is determined with less accuracy, generally by all methods. The same inferences and trend have been observed for scenarios 2-5 but the level of matching the results degraded because with more number of mesh points, i.e., for smaller values of h, the smooth results are obtained which is not possible for few mesh points. Additionally, it seems that small variations in the results are observed in each case of all five scenarios in the study, but closely seen reveals that for larger input span the variation in the results is rather more frequent.
The statistical indicator based on the minimum (MIN), mean and standard deviation (STD) values are calculated for 100 independent runs of the proposed scheme for each case of all five scenarios of the problem and results are given in Table 2 for five inputs. While in Tables 3 and 4 the statistical indices are given for more intermediate inputs to analyze dynamics of the problem for few selected cases and scenarios of TFE. From the values presented in these tables, no noticeable difference is apparently observed between the MIN and mean values because of very small values of STD for each case of all five scenarios. Additionally, it is seen that from smaller input span t ∈ [0, 1] to larger inputs t ∈ [0, 100], the values of STD degraded but still remains of the order of 10 −09 , which established the consistency of the proposed methodology for solving TFE.
Comparative study for the proposed solution of TFE is made with the results of existing techniques based on Rational Chebyshev PseudoSpectral Method (RCPSM) (Parand and Shahini 2009) Table 5 for inputs t ∈ [0, 5]. It can be seen that trend of the proposed results is aligned with the similar patterns of state of the art numerical and analytical solutions. Moreover, comparison of the results is presented on larger inputs span t ∈ [0, 100] in Table 6 for the proposed and reference solver based on HAMTA (Khan and Xu 2007), HAM (Liao 2003b) and Chebyshev PseudoSpectral Method (CPSM) (Kılıçman et al.

Comparative analysis on global performance operators
Comparative analysis of the results is made on the basis of the mean value of fitness to analyze the accuracy and convergence and in case of examining the complexity operator based on mean time, mean generations and mean function counts are incorporated. The mean and STD values of fitness, time, generations and function counts are given in Table 7 for each case of all five scenarios of TFE. It is seen quite apparently that with the increase in length of input span, i.e., moving from scenario 1 to scenario 5, the accuracy of the algorithms decreases due to the fact that for large input span, handling of nonlinearity and singularity associated with TFE is rather more difficult.
It is seen from the results presented in Table 7 that the values of mean time for optimization of fitness function by GA-SQP algorithm in case of 11, 21, 51 and 101 input grid points are 6 ± 1, 11 ± 1, 22 ± 1, and 50 ± 1 s. While the number of generations is around 440 ± 10, 480 ± 10, 600 ± 20, and 720 ± 25 for 11, 21, 51 and 101 input grid points and respective values for function counts are around 60,950 ± 60, 63,400 ± 200, Fig. 3 Optimization output plots of GA-SQP to solve scenario 1: case 1 of Thomas-Fermi equation 76,500 ± 700 and 125,000 ± 14,000. With an increase in the number of grid points, the complexity of the fitness function increases due to which larger values of complexity operator are obtained.

Conclusions
Although there has been achieved a number of solutions of Thomas-Fermi Equation but the highly impressive potent outcome of our study are summarized by the following concluding remarks: • A novel design is presented for solving nonlinear singular Thomas-Fermi equation with the help of finite difference method for discritization of problem and resultant system of nonlinear equations are solved by exploiting the strength of bio-inspired   • Comparison of the results on the basis of mean values of residual error or fitness function shows that for wider input span the performance of the proposed scheme is slightly degraded but still achieving the values of fitness of the order 10 −09 . • Computational complexity of the proposed scheme is examined through an average values of time, generations and function counts being consumed by GA-SQP algorithm for optimization of fitness functions and it is found that with increase in the number of grid points, the values of complexity operators are on the higher side because in these cases more computational drive is required for solving the decretization model of the equation. • Beside the consistent accuracy and convergence, other perks of the proposed scheme are simplicity of the concept, easy implementation and a good alternate avenue to be exploited for solving the nonlinear and singular problems for which the conventional methodologies fail.
In future, present research may prove to be a beacon for researchers working in the domain of application of artificial intelligence techniques to stiff problems arising in physical models of practical importance.