A novel neuroevolutionary paradigm for solving strongly nonlinear singular boundary value problems in physiology

In the present research work, we designed a hybrid stochastic numerical solver to investigate nonlinear singular two-point boundary value problems with Neumann and Robin boundary conditions arising in various physical models. In this method, we hybridized Harris Hawks Optimizer with Interior Point Algorithm named HHO-IPA. We construct artificial neural networks (ANNs) model for the problem, and this model is tuned with the proposed scheme. This scheme overcomes the singular behavior of problems. The accuracy and applicability of the method are illustrated by finding absolute errors in the solution. The outcomes are compared with the results present in the literature to demonstrate the effectiveness and robustness of the scheme by considering four different nonlinear singular boundary value problems. Further, the convergence of the scheme is proved by performing computational complexity analysis. Moreover, the graphical overview of statistical analysis is added to our investigation to elaborate further on the scheme’s stability, accuracy, and consistency.


I. INTRODUCTION
The class of nonlinear singular differential equations performs a vital role in many fields of science and technology. The mathematical modeling of real-life problems arises in several physics branches, specifically astrophysics, thermodynamics, physical chemistry, nuclear technology, atomic energy, and engineering. For example, the thermal explosion [1]- [3] and heat conduction [4], are modeled in the form of singular boundary value problems. Many nonlinear problems are modeled in the form of boundary value problems based on given boundary conditions. The boundary conditions used to model real-life problems are ordinarily Neumann boundary conditions, mixed boundary conditions, and Dirichlet boundary conditions. The most physically and reasonable choice for the modeling of real-world problems is Neumann boundary conditions [5]. We consider the following general class of nonlinear singular, two-point boundary value prob-lem (BVP): (p(r)f (r)) = p(r)g(r, f (r)), 0 < r ≤ 1 (1) subject to the following Robin and Neumann boundary conditions: where µ ≥ 0, η ≥ 0 and ρ ≥ 0 are finite real constants and g is a nonlinear term. If p(r) = r and the non-linear term g of the following type: g(r, f (r)) = −ke f (r) (3) where k is real constant, then equation (1) represents isothermal explosion [1], [6], [7]. If p(r) = r and the nonlinear term g of the following form: g(r, f (r)) = e f (r) (4) then equation (1) represents equilibrium of thermal explosion in cylindrical vessel [1], [3]. If p(r) = r and the nonlinear term g of the following form: g(r, f (r)) = −f 5 (r) (5) then equation (1) arising in the study of equilibrium of isothermal gas sphere [1], [2], [8].
If p(r) = r 2 and the nonlinear term g of the following form: where c 1 and c 2 are real constant. Then equation (1) represent steady-state oxygen diffusion in a spherical cell [2], [9]- [11]. The numerical solution of nonlinear singular boundary value problems (SBVPs) of ordinary differential equations is a difficult assignment due to the singular behavior of the problem at the origin. As almost every real word problem is modeled in the form of nonlinear SBVPs [12] and the analytical solutions of most of the problems are not available, to handle this type of problem researchers make contact with approximation theory. The recent literature proved that various researchers proposed numerous approximation methods to tackled SBVPs numerically and analytically like the finite difference method (FDM) [1], the monotonic iterative method of Bessel functions [12], [13], an improved iterative technique [14], finite-element method (FEM), monotonic iterative scheme containing the expansion of eigenfunction [15], differential quadrature method [16], modified adomian decomposition method (MADM) [17], collocation method [18], Borel-Laplace transformation technique [19], Adomian decomposition method (ADM) [20], variational iteration method (VIM) [21], approximate power series solution method [22] and homotopy perturbation method (HPM) [23] etc.
The algorithm like VIM [21], ADM [20], MADM [17] and HPM [24] provide solution of an under consideration problem in the form of series. The convergence of the series solution obtained by approximation techniques is not guaranteed. In 2017 Biswal and Raoul proposed an algorithm [24] based on the combination of an integral equation formalism and optimal homotopy analysis method for the solution of problem (1) by replacing g(r) subjected to boundary conditions in (2). In addition, the analysis of different approximation solution methodologies, like VIM [21], ADM [20], or MADM [17] shows that these algorithms required the undetermined coefficients. At every iteration of approximation, it needs to solve a transcendental or an algebraic equation. The solution of transcendental or an algebraic equation at each iteration is a time-consuming process, and as a result, the computational complexity of the algorithm increase.
The advanced adomian decomposition method (AADM) is a numerical method that is developed based on a numerical technique named Adomian decomposition method proposed by Adomian [25]. The main reason for designing AADM is to deal with nonlinear problems without considering unrealistic norms like discretization and linearization and improve the convergence as the author claimed in [26]. But to achieve this goal, the computational complexity of the algorithm increase. In addition, the Advance adomian decomposition method needs the solution of adomian polynomial at each iteration which is also a time-consuming task.
Some other modifications are made to improve the rate of convergence and accuracy of adomian decomposition method (ADM) Adomian [25] and named it as modified adomian decomposition method (MADM) [27]. The adomian decomposition method is modified to overcome the singular behavior of initial and boundary value problems at x = 0. For this, the author introduced the Laplace operator into ADM, which caused the algorithm more complicated. In addition, adomian decomposition methods provide a small region of convergence [28].
In the list of numerical solvers, the compact finite difference method (CFDM) is a more favorable numerical technique for solving problems like hyperbolic problems, linear and nonlinear boundary, and initial value problems. The compact finite difference method is also utilized for the approximation of singular boundary value problem [29]. The compact finite difference method is used for the approximations of finite-difference. These approximations are more ideal when their dissipative error and dispersive error are compared with explicit algorithms.
But a disadvantage of the compact finite difference algorithm is that it needs the solution of the system of a diagonal matrix for the derivatives at all grid points or the calculation of interpolations.
Most of the real-life problems are nonlinear and contain singularity [28]. To remove the singularity of the problem, researchers decompose the domain of the problem into two subdomains. The singular behavior of the problem in the first subdomain is removed with the application of the decomposition method based on the integral operator. The resulting system is dealt with cubic B-spline collocation method [30], which increases the computational complexity of the algorithm. The Cubic B-spline collocation method also needs the solution of a matrix.
From the above discussion, we conclude that the MDM-CBC method has a small region for convergence due to the involvement of the decomposition method, which may cause the divergence of the algorithm. Moreover, to remove the singularity of the problem, MDM-CBC decomposes the given interval into two subintervals that modify the original problem.
The disadvantages related to the convergence of the series solution and computational complexity of the above listed numerical approximation algorithms motivated us to develop a stochastic algorithm based on artificial neural networks (ANNs) for the numerical treatment of the problem given in equation (1) subjected to two different kinds of boundary conditions given in (2). Besides, the perturbation algorithms

Input
Hidden Laxer Output r r (0,1) depend on small parameters, but the proposed scheme does not require any small parameter present in the given problem. Also, it does not require any prior information about the derivatives of the differential equations. The proposed stochastic algorithm provides a set of optimized weights for a given problem with easily computable components and provides a solution in the form of convergent series. The advantage of our designed scheme over the numerical method listed above is that it approximates the probability of getting optimal solutions. Furthermore, our designed scheme has less computational complexity. The proposed algorithm is developed based on a minimizing function of mean square error in the differential equation together with mean square error in boundary conditions called fitness function to get the optimal solution of a given problem. The elementary structure of our proposed scheme is quite simple, which comprises of three parts, in the varied first part of the method, we designed ANNs model for under consideration SBVP, in the second part, a global search algorithm is applied to get the global optimal solution and in last part of the proposed scheme a local search technique is used to tune the optimal results obtained from global search algorithm finely.
The main goal is to develop a soft computing algorithm to analyze the singular boundary value problem. We hybridized a global search optimizer, namely Harris Hawks Optimizer (HHO), through a local search paradigm called Interior Point Algorithm (IPA) and labeled this hybrid scheme as HHO-IPA. Our proposed hybrid optimizer globally tunes the unknown parameter of the ANNs model by global search algorithm HHO. The best weights of HHO are finely trained through local search optimizer IPA.
The accuracy and efficiency of our designed technique are demonstrated by considering four singular boundary value problems examined by the different researchers using numerical techniques. The comparison of numerical outcomes of our designed technique is made with some existing results in the literature. This comparison of results proves that our designed technique provides a much better solution than the semi-numerical techniques.
The outcomes of HHO-IPA are also compared with an advanced artificial intelligence-based algorithm named Slime Mould Algorithm (SMA). Slime Mould Algorithm is a bioinspired algorithm and was first introduced by Shimin Li and Huiling Chen in 2020. Slime mold algorithm mimic morphological changes and foraging behavior of slime mold Physarum polycephalum [31].
The rest of the paper is structured as, in Section 2, the designed hybrid algorithm HHO-IPA is discussed in detail. The modeling of ANNs for a given problem and the for-mulation of the fitness function is reported in Section 3. The accuracy and efficiency of the proposed algorithm are narrated in Section 4 by considering four different problems. Section 5 consists of our conclusion about this research.

II. HARRIS HAWKS OPTIMIZER
Harris Hawks Optimizer (HHO) is a nature-inspired, gradient-free optimization, population-based algorithm. HHO was first familiarized by Ali Asghar Heidari in 2019 in his unique work in the class of intelligent computing [32]. The exploration and exploitation phases of Harris' Hawks Optimizer is modeled on the basis of the cooperative behaviors, exploring a prey, surprise pounce and chasing style of hawks. HHO can be implemented to optimization problem with proper formulation. Different phases of HHO are described briefly in later subsections and shown in figure 2.

A. EXPLORATION PHASE
In exploration phase of HHO, hawks track and spot the prey. The position of hawks is the candidate solution.In HHO, hawks randomly perch on the basis of two main strategies. In the first strategy, hawks perch on the basis of location of prey and other group members following the condition q < 0.5 while in later strategy the hawks perch on random position (tall trees) following the condition q ≥ 0.5. The mathematical modeling of both strategies are given in equation (7).
2r 2 X(t)| q ≥ 0.5 (X rabbit (t) − X m (t))− r 3 (LB + r 4 (U B − LB)) q < 0.5 (7) where X(t + 1) in the above equation represents the location of hawks in iteration t, X m , X rand (t), X(t) and X rabbit (t) represent average position of hawks, randomly selected hawk, position vector of hawks and location of prey (rabbit) respectively. U B and LB represent the upper and lower bounds while q, r 1 , r 2 , r 3 and r 4 are random numbers between (0,1). The average position X m of hawks in above equation (7) is calculated through equation (8) given below: where t is current iteration while N and X i (t) in equation (8) represent the total number of Harris' hawks and the position of hawk respectively.

B. TRANSITION FROM EXPLORATION PHASE TO EXPLOITATION PHASE
The HHO has ability to transfer to exploitation from exploration. The change of HHO among various exploitative behaviors depend upon the escaping energy of the rabbit (prey). Due to the escaping behavior the escaping energy of rabbit significantly decrease which are model as following: where maximum number of iterations are represented by T . E 0 and E in above equation indicate initial state of prey's energy and its escaping energy respectively.

C. EXPLOITATION PHASE
In exploitation phase, the hawks try to catch the detected prey while the prey attempt to escape. Four possible strategies of hawks occurs due to the chasing style of hawks and escaping style of prey. Which are modeled in the following subsections.

1) Soft besiege
This strategy occur when the prey has enough energy |E| ≥ 0.5 to escape and r ≥ 0.5. The mathematical modelling of this behavior is given as: where J = 2(1 − r 5 ) indicate the strength of random jump of rabbit, r 5 is random number inside the interval (0,1), the value of J randomly vary in iteration t and ∆X(t) is difference of current position in iteration t and position of rabbit.

2) Hard besiege
This strategy occur when the prey so tired and has not enough energy |E| <0.5 to escape and r ≥0.5. The mathematical modelling of this behavior is given as: 3) Soft besiege with progressive rapid dives This strategy occur when the prey has still enough energy |E| ≥ 0.5 to escape but r < 0. 5. This indicate that the prey can successfully escape and hence the Levy flight concept introduced to model this behavior: where D indicate the problem's dimension, S is a random victor of size 1 × D and LF is Levy flight function given as following: where β is a default constant of value 1.5. u and v are random numbers in (0 1). In soft besiege the locations of hawks can be update using equation (16).
where Y and Z victors can be evaluate using equations. (13) and (14).

4) Hard besiege with progressive rapid dives
This strategy occur when the prey has not enough energy |E| < 0.5 to escape but r < 0.5. In this case, the hawks attempt to decrease the distance of their average position with the escaping prey. This situation is modeled as: where Y and Z victors can be evaluate using the following new rules in equations (18) and (19).
where X m (t) is calculated by using equation (8).

III. ARTIFICIAL NEURAL NETWORK
The solution of complex multivariate non-linear relationships motivate the researchers to design an intelligent system. They utilized the architecture of the human brain to designed these intelligent system. In the way of development of intelligent systems, Warren McCullough and Walter Pitts publish his mathematical model of neural networks (NN) [33]. Neural network is the set of neurons. A network which consist of artificial neurons is known as artificial neural network (ANN) [34] or biological neurons called biological neural network [35]. ANNs use the computational or mathematical model which is used for processing available information is an interconnected set/group/collection of artificial neurons, in other words artificial neural networks are utilize to model the complex relationship among the input and output data. The strength of artificial neural networks is utilized by various scientists to get the solution of different real-life problems [36] based on artificial intelligence (AI) models.
The construction of ANN model is based on the human brain's layout which are modeled through artificial neurons. These networks emulate the human brain layout [37]but they use a reduced group of biological neural network.
The structure of ANNs consist of three or more layers, where these layers are interconnected through connections. The first layer is the input layer which is some time called receptor that receive information in the form of text, numbers etc. The neurons of input layer calculate the weighted sum of all input data and send this weighted sum to the hidden layer(neural zero layer). The hidden layer is core of neural system, which receives data from input layer in the form of weighted sum and produce the output through an activation function. Activation function is sometime called transfer function [38], due to the property of transferring data from hidden layer to output layer. The selection of activation function has a great impact on the capability and performance of neural networks [39]. Generally an activation function is a nonlinear function [40], [41] or orthogonal polynomials [42]. Hyperbolic tangent function, logistic activation function, binary step function [38], [43], [44], Hermite polynomials [45], Legendre polynomials [46], Chebyshev polynomials [47] etc are some examples of activation functions. These activation functions decide whether send data to output layer or not. The output layer (results layer) is the final layer of ANNs which generate the final result. The structure of artificial neural network is depicted in figure (3).

IV. PROPOSED SCHEME
The class of evolutionary algorithms (EAs) [48] play a vital role in the solution of real-world problems [49] in different research areas like chemistry [50], machine learning [51], finance [52], mechanics, economics, optics, biology, geophysics, biochemistry, etc. Stochastic optimizers are currently widely applied as global search optimization methodologies. These algorithms mimic biological or social behavior or nature, that's why these are known as nature-inspired algorithms, Grey Wolf optimizer (GWO) [53], Whale Optimisation Algorithm (WOA) [54], Harris Hawks Optimizer (HHO) [32], etc, are some common examples of natureinspired algorithms. Nature-inspired algorithms are almost population-based approaches that learn from the set or group of data.
The hybridization of algorithms is a trend of 21 st century. Hybridization improves the ability and capability of algorithms such as accuracy, convergence, computational speed, etc. In a hybrid scheme, the advantages of each algorithm are combined to overcome the disadvantages of each algorithm. Hybridization is fallen into different kinds: one of these is qualified as sequential hybridization [55]. In this kind of hybridization, a set of algorithms implement one after another. In other words, in sequential hybridization, the outputs of one algorithm are used as input of another algorithm. These algorithms are also categorized as sequential hybrid metaheuristics [56]. The popularity of hybrid algorithms increases because of their better performance when dealing with imprecision, uncertainty, vagueness, nonlinearity, and singularity, etc. The advantage/importance of hybrid evolutionary approaches increase when finding the global optimal value in complex and, vast search domains and especially when the gradient information of the problem is not available. From an application point of view, these algorithms are relatively simple and the concepts are easy to understand. The parameters of these algorithms are reasonably flexible, in other words, parameters can be variate for better performance of the algorithm. Besides, the implementation of hybrid meta-heuristic algorithms is easy and straightforward.

A. PROBLEM FORMULATION BASED ON ARTIFICIAL NEURAL NETWORKS
The concept of Artificial Neural Networks (ANNs) was initiated by a United States logician in the area of computational neuroscience Walter Pitts and an American cybernetician and neurophysiologist Warren McCulloch through his computational model [57] in 1943. This model provides a way to divide research into two branches: the first approach concentrated on the biological process while the second one focused on the applications of neural networks (NNs) to artificial intelligence (AI). The artificial neural network is the component of a soft computing system that simulates the human brain layout. In other words, the artificial neural network is the foundation of artificial intelligence. The subject of artificial neural networks has great attention of researchers in the community of machine learning due to the capability and ability to solve challenging problems, particularly the practical applications of ANNs in the area of industry, medical diagnosis, personal communication, speech recognition, encompassing finance, object recognition, education and image processing, etc. Artificial neural networks are more interesting in computational mathematics and approximation theory due to their ability to successfully approximate the solution of almost every kind of differential equations (ordinary and partial differential equations), arbitrary functions, and system of differential equations such as appearing in chemistry, engineering, mechanics and biotechnology. Recently many articles from several kinds of research are published on this topic. Most of the researchers considered meta-heuristic and heuristic schemes which are based on artificial intelligence. These artificial intelligence-based algorithms are applied to governing equations and their boundary conditions in a unified search zone.
In this work, we employ ANNs formulation in the form of continuous mapping for solving 2nd order singular boundary value problems. The n th order derivatives of these networks

1) Fitness Function
A fitness function is operated to summarise designed solution for achieving the set desired aims. The fitness functions for singular boundary value problem in Equations (1)-(2) guide simulations towards optimal design solutions. It is formulated in the form of square of residual errors (mean square errors) as minimization objective function as following: In above Equation (22), ε 1 represent mean square error in singular boundary value problem (1)-(2) which can be written as: wheref and h is the step size. While mean square error in boundary conditions of given problem is represented by ε 2 and formulated as: The unknown parameters w = [α, β, γ] of artificial neural networks model for getting the solution of singular boundary value problem are trained by using proposed hybrid scheme HHO-IPA. The parameters are trained until the fitness value of SBVPs approximately equal to zero. The solution of SBVPs, in this case near the exact solution. i.e., ε → 0 then f (r) → f (r).   'Fitness' less than or equal to 10 −13 'total iterations' less than or equal to 3000 'TolFun' less than or equal to 10 −18 'TolX' less than or equal to 10 −20 'TolCon' less than or equal to 10 −18 'MaxFunEvals' less than or equal to 200,000 while (Terminate required achieved)

B. HYBRID SCHEME HHO-IPA
We provide a detailed description of our proposed hybrid soft computing algorithm HHO-IPA for the investigation of nonlinear singular boundary value problems arising in physiology. Generally, our soft computing scheme is narrated in two stages, in the first stage of HHO-IPA we model the SB-VPs (1) in the form of fitness function based on ANNs, while the unknown parameters of ANNs model, detail steps of learning methodology and setting of the essential parameter are presented for proposed hybrid scheme of "Harris Hawks Optimizer" and "Interior Point Algorithm" is presented in a later stage.
Harris Hawks Optimizer (HHO) is a nature-inspired population-based algorithm. HHO was first familiarized by Ali Asghar Heidari in 2019 in his unique work in the class of intelligent computing [32]. The main inspiration of Harris' Hawks Optimizer is the cooperative behaviors and chasing style of hawks. It is considered one of the most intelligent birds.
Harris Hawks Optimizer can solve constraint, uncon-straint, linear, nonlinear, singular boundary value problems (BVP), and initial value problems (IVPs) arising in different fields of science and engineering. The systematic flowchart of HHO is given in Figure (4). Harris Hawks Optimizer finds near best candidate solution of under consideration problem in a given domain of the problem because of its better global search capability. This global search capability of HHO is boosted through Interior Point Algorithm (IPA). IPA is an individual-based algorithm that follows a single path to approach the results. The global best results of HHO are considered as starting point of IPA to get the solution rather quickly. The graphical representation of the Interior point method is in Figure (5). Different kinds of problems are effectively solved in optimization theory with the help of IPA, like hyperbolicity cone problems [58] the approximation of parameters of discrete-time infective disease model [59] and nonlinear non-convex programming [60]. The accuracy and convergence of the outcomes toward the exact solutions of the problems inspired us to implement IPA for the solution of nonlinear singular boundary value problems of the second order with Robin and Neumann boundary conditions.
We proposed a new soft computing scheme by combining the global search capability of HHO and the local search strength of IPA namely HHO-IPA. The working procedure of the hybrid scheme is demonstrated in Table (3) in the form of pseudo-code.
HHO-IPA is implemented to train the unknown parameters of the artificial neural networks model for the solution of SB-VPs. The global search of HHO-IPA is performed by using MATLAB script file on the other hand we use Matlab's builtin function "fmincon" routine algorithm "interior-point" for IPA.

C. ERRORS ESTIMATION
The accuracy of our proposed algorithm is tested by calculating absolute errors in our solution. These absolute errors are compared with some existing in the literature. The absolute error is represented by e n (r) and given as: here n represent is the number of solution points andf (r) is the solution approximated through the designed method while f (r) in the above equation is an exact solution of the problem at the corresponding point.

V. NUMERICAL EXPERIMENTATION
This section of the paper composes four different kinds of real-life nonlinear, singular boundary value problems. The comparison of experimental outcomes of the designed scheme is made with some existing solutions in the literature to demonstrate the accuracy, applicability, and efficiency of the designed technique. The solutions are given in the form of tables and also illustrated graphically. Example-1: Consider the singular boundary value problem emerges in the isothermal explosions: subjected to boundary conditions The problem (26) - (27) corresponds to (1) -(2) with p(r) = r, g(r, f ) = −ke f (r) , µ = 1, η = 0 and ρ = 0. Exact solution for this problem found in the literature is given as: The problem in Equation (26) (1) and (2). The formulation of fitness function (FF) (22) for given problem is below: The formulated fitness function (28) is tuned by using HHO-IPA and a set of optimized weights is obtained with the fitness value of 2.1165E − 14. This set of optimized weights is used for the derivation of solution. These weights are plotted in Figure ( The solutions obtained for singular boundary value problem given in Equation (26) subjected to boundary conditions in (27) are cross examined by getting the absolute errors in our solution. The absolute errors are calculated by using Equation (25). To prove accuracy of HHO-IPA these absolute errors are compared with some existing in the literature for the same problem. The absolute errors are presented in Table (4) along with some existing absolute errors obtained by the modified decomposition method and cubic B-spline collocation (MDM-CBC) [30], modified Adomian decomposition method (MADM) [27] and advanced Adomian decomposition method (AADM) [26]. Table (  confirms that the magnitude of errors in HHO-IPA solution is two orders lesser than the magnitude of errors of AADM. The accuracy of 10 −08 to 10 −11 achieved by our proposed scheme based on hundred independent runs. The comparison of errors proved that our proposed algorithm is efficient and accurate as compare to advanced Adomian decomposition method described in [26]. For instance, the logarithmic plot of the absolute errors obtained from the designed scheme are depicted in Figure ( Example-2: Consider the problem appears in the investi-gation of thermal explosion of gas in cylindrical container. This is a nonlinear boundary value with a singularity [1]- [3]: with boundary conditions The problem (30) -(31) corresponds to (1) -(2) with p(r) = r, g(r, f ) = e f (r) , µ = 1, η = 0 and ρ = 0. Exact solution for this problem found in the literature is given as: The problem given in Equation (30) is singular ordinary differential equation subjected to boundary conditions in (31) form a nonlinear, singular boundary value problem (SBVP). We solve problem (30) -(31) using the proposed HHO-IPA by using ANNs model with in the interval of [0 1] by taking step size of 0.1. The setting of essential parameter are tabulated in Tables (1) and (2). The given fitness function (FF) in Equation (22) is formulated for this problem as below: The formulated fitness function (32) is tuned by using proposed hybrid algorithm HHO-IPA and the set of finely optimized weights is obtained with the fitness value of 1.9727E − 12. This set of optimized weights are used for the derivation of solution. These weights are plotted in Figure  (8b) in the form of bar graph and given in Equation (33). The solutions obtained for singular boundary value problem given in Equation (30) subjected to boundary conditions in (31) are cross examined by getting the absolute errors in our solution. The absolute errors are calculated by using Equation (25). To prove accuracy of HHO-IPA these absolute errors are compared with some existing in the literature for the same problem. The absolute errors of the problem are listed in Table (5) together with the existing absolute errors attained by finite difference method of order four (4 th -FDM) [1] and compact finite difference method (CFDM) [29]. The comparison between our designed scheme and compact finite difference method (CFDM) proved that HHO-IPA is highly accurate. Table (5) confirms that the errors in HHO-IPA solution is smaller in magnitude than the CFDM. [29]. Furthermore, the absolute error for m = 16, in HHO-IPA solution is 3.8619E − 08, however, for m = 16 CFDM reduced absolute error about 5.3227E − 07. The accuracy of 10 −08 to 10 −10 achieved by our proposed scheme based on hundred independent runs. For instance, the logarithmic plot of the absolute errors accomplished by our designed scheme are depicted in Figure ( Example-3: Consider the problem arising in the study of equilibrium of isothermal gas sphere form a nonlinear boundary value problem: (r 2 f (r)) = −r 2 f 5 (r), (34) subjected to boundary conditions The problem (34) -(35) corresponds to (1) -(2) with p(r) = r 2 , g(r, f ) = −f 5 (r), µ = 1, η = 0 and ρ = 3 4 . Exact solution for this problem is given as: The system in Equation (34) (1) and (2). The fitness function (FF) given in Equation (22) for this problem is constructed as: The formulated fitness function (36) is tuned by using HHO-IPA and a set of optimized weights is obtained with the fitness value of 2.4045E − 13. This set of optimized weights are used for the derivation of solution. These weights are plotted in Figure (11a) in the form of bar graph and given in Equation (37).     (6) along with the existing absolute errors obtained by the advanced Adomian decomposition method (AADM) [26] and He's variational iteration method (HVIM) [61]. Table (6) confirms that the errors in HHO-IPA solution is smaller in magnitude than the AADM [26]. The accuracy of 10 −06 to 10 −08 achieved by our proposed scheme based on hundred independent runs. The comparison proved that HHO-IPA is highly accurate than advanced Adomian decomposition method (AADM) [26] and He's variational iteration method (HVIM) [61]. For instance, the logarithmic plot of the absolute errors obtained by the designed scheme are depicted in Figure ( subjected to boundary conditions: with c 1 > 0, c 2 > 0 and ν > 0. Here, ν, c 1 and c 2 are the permeability of cell membrane, Michaeli's constant and the rate of maximum reaction respectively. Michaeli's constant shows the half-saturation concentration. In above system, the radial distance of spherical cell is represented by r while f shows the oxygen concentration. This problem have no exact solution in the literature The problem in equations (38) and (39) corresponds to (1) and (2) with p(r) = r 2 , g(r, f ) = c1f (r) f (r)+c2 , µ = ν, η = 1 and ρ = ν. We solve problem (20) -(20) using the proposed method for the values parameters c 1 = 0.76129, c 2 = 0.03119 and ν = 5.
The problem in Equation (38) is also a singular ordinary differential equation subjected to the boundary conditions in (39) form a nonlinear, singular boundary value problem (SBVP). We solve problem (38) -(39) using the proposed HHO-IPA by using ANNs model with in the solution domain [0 1] by taking the step size of 0.1. The setting of essential parameter are tabulated in Tables (1) and (2). The fitness function (FF) (22) is formulated for this nonlinear singular boundary value problem as following: where c 1 = 0.76129, c 2 = 0.03119 and ν = 5.
The formulated fitness function (40) is tuned by using proposed paradigm namely HHO-IPA, which achieved the fitness of 1.1080E − 12 and provide us a set of finely optimized weights. This set of optimized weights are used for the derivation of solution. These weights are plotted in Figure (11b) in the form of bar graph and given in Equation (41).
The solutions obtained for singular boundary value problem given in Equation (38) subjected to boundary conditions in (39) are tabulated in Table (7) together with some existing solutions approximated by the compact finite difference method (CFDM) [29], advanced Adomian decomposition method (AADM) [26], differential transform method (DTM) [16], quartic B-spline collocation method (QB S CM) [62] and improved decomposition method (IDM) [20]. The designed scheme is highly accurate when compared to the advanced Adomian decomposition method (AADM) [26]. The solutions of HHO-IPA are presented in Figure ( Figure (10a). The solution graph of the problem shows that the oxygen concentration f is monotonically increasing.
The solutions obtained for singular boundary value problem given in Equation (38) subjected to boundary conditions in (39) are cross examined by getting the absolute errors in our solution. The absolute errors are calculated by using Equation (25). For instance, the absolute errors achieved by the designed scheme are plotted on logarithmic scale in Figure (10b) in form of Minimum absolute error, Mean absolute error and Maximum absolute error. The graph of absolute errors shows that the it is little increases in the direction of the point of singularity but at that point (singular point) i.e. r = 0, the absolute error is dramatically falls to zero and the error decrease away from the point of singularity and minimum at i.e. r = 1.

VI. STATISTICAL ANALYSIS
This section of the paper consist the comparative study of our results obtained through proposed algorithm HHO-IPA. The performance of HHO-IPA is proved through performance indicators to proved the accuracy, consistency, efficiency and reliability of our soft computing scheme. The results obtained from HHO-IPA are analyzed through indices like Fitness (FIT), Mean Absolute Deviation (MAD), Theil's inequality      (22), (42), (43), (44), (45) respectively as following: where  (12), (14), (16), (18) and (20) respectively while the graphs for example-3 and example-4 are shown respectively in Figures (13), (15), (17), (19) and (21). All the outcomes tabulated in the term of Mean and standard deviation (STD) in table (8). In addition to prove the consistency and stability of proposed scheme the results of performance indicators are plotted on semi-log scale in the form of bar graph in term of Best, Mean and worst result for example-1, example-2, example-3 and example-4 respectively in Figures (22a), (22b), (23a) and (23b).
This performance measurement is further extend by finding the values of global performance indices like global Fitness (GbFIT), global mean absolute deviation (GbMAD), global root mean square error (GbRMSE), global error in nash-sutcliffe efficiency (GbENSE) and global Theil's inequality coefficient (GbTIC). The mathematical formulation of GbFIT, GbMAD, GbENSE, GbRMSE and GbTIC for singular boundary value problems arising in physiology are listed in equations (46), (47), (48), (49) and (50), respectively as following: where in above listed equations, R n represents total numbers of executed runs, n is current number of run, m is the current solution point, ε n is objective value of n th experiment and N is shows the total solution points of the problem. The termf (r m ) in the above systems represent the approximate solution obtained by using HHO-IPA while the term f (r m ) is exact/standard solution of the given problem. The results of global performance indicators obtained by considering input between r [0, 1] by taking the step size of 0.1. In present case the values of parameters R n and N are taken 100 and 11 respectively. The results of global performance indicators, i.e., GbTIC, GbRMSE, GbMAD, GbENSE and GbFIT are tabulated in table (9) in the form of standard deviation (STD) and Mean. Additionally it is clear from the table that the global performance indicators has less values, which proved the consistency and accuracy of proposed algorithm for the solution of nonlinear singular boundary value problem considered in this research article.
The consistency of proposed algorithm is further analyzed by performing the computational complexity analysis (CCA) test of the results. CCA is implemented for the results obtained from hundred independent runs of the designed technique. The experimental results of CCA are based on population creation, function evaluation and average time consumed by our algorithm for the training of unknown factors of designed model of artificial neural networks. The outcomes are tabulated in table (10) in the terms of mean and standard deviation (STD) for all problem. Table (10) shows that the values of population creation for problem-1,       In this research, we have developed an intelligent soft computing technique, namely, the HHO-IPA algorithm. To validate the efficiency and application of our algorithm, we have analyzed a problem involving nonlinear and singular differential equations that have applications in physiology. Four different problems are considered in this paper. Our experimental outcome shows that artificial neural networks based HHO-IPA procedure is capable of solving nonlinear and singular differential equations modeling real-life issues. Below, we present key advances/ contributions in this paper, • We have constructed a generalized series solution using feed-forward artificial neural networks architecture. • An error function is suggested for the problem considered in this paper. We have developed a hybrid optimization technique by combining Harris Hawks Optimizer and Interior Point Algorithm to minimize the error function. • The detailed steps of our proposed scheme are presented in the table (3). • From tables (4) and (6), it is obvious that our approach improves the quality of the solution. • From table (5), it is clear that by increasing the total number of solution points, the quality of the solution is improved by our designed scheme. • Table (7) evident that the solution obtained by our approach is much better than other methods in the literature. • Performance indicators are used to show the quality of solutions, efficiency, and reliability of our algorithm. We have calculated values of FIT, MAD, RMSE, ENSE, TIC, GbFIT, GbMAD, GbRMSE, GbENSE and GbTIC, see tables (8) and (9). • Computational complexity is measured through statistical values, success rates, population sizes, function evaluations, and time taken (see tables (10).
It is established that the neural networks-based HHO-IPA algorithm is a strong candidate as a solver for problems involving singularity and nonlinearity.