Nature Inspired Computational Technique for the Numerical Solution of Nonlinear Singular Boundary Value Problems Arising in Physiology

We present a hybrid heuristic computing method for the numerical solution of nonlinear singular boundary value problems arising in physiology. The approximate solution is deduced as a linear combination of some log sigmoid basis functions. A fitness function representing the sum of the mean square error of the given nonlinear ordinary differential equation (ODE) and its boundary conditions is formulated. The optimization of the unknown adjustable parameters contained in the fitness function is performed by the hybrid heuristic computation algorithm based on genetic algorithm (GA), interior point algorithm (IPA), and active set algorithm (ASA). The efficiency and the viability of the proposed method are confirmed by solving three examples from physiology. The obtained approximate solutions are found in excellent agreement with the exact solutions as well as some conventional numerical solutions.


Introduction
The numerical treatment of nonlinear singular boundary value problems (BVPs) has been considered by several authors in the last few decades due to their varied applications in the fields of engineering and science such as gas dynamics, atomic structures, atomic calculations, and chemical reactions [1]. Many methods including finite difference method, Chebyshev polynomial, B-spline method, and nonpolynomial cubic spline have been employed to handle singular boundary value problems. The reader may find a comprehensive survey of computational techniques utilized for the numerical solution of singular boundary value problems in [1].
The key objective of this paper is to present a stochastic computing technique for the numerical solution of nonlinear singular boundary value problems of the following form [2]: ( ) + ( + ) ( ) = ( , ) , 0 ≤ ≤ 1, (1) The assumptions normally applied on ( , ) are that it is continuous, / exists and is continuous, and / ≥ 0, for all 0 ≤ ≤ 1. The singular boundary value problem (1)-(3) occurs in numerous applications, especially with = 0, 1, 2 and = 0 in the study of many tumor growth problems [3,4] and with = 2 and = 0 in the study of steady state oxygen diffusion in a spherical cell with Michaelis-Menten uptake kinetics [5][6][7][8]. A similar problem for = 2 and = 0 also arises in the study of the distribution of heat sources in the human head [9][10][11].
An incredible amount of research work has been invested for the study of the singular boundary value problems (BVPs) of the form (1)- (3), and different analytical and numerical methods have been utilized [2,[12][13][14][15][16]. Among many authors, Rashidinia et al. employed nonpolynomial cubic spline method [2], Pandey and Singh used finite difference method [12], Khuri and Sayfy recently proposed a method based 2 The Scientific World Journal on the combination of modified decomposition method and cubic B-spline collocation [13], Ç aglar et al. applied B-spline method [14], Mittal and Nigam used Adomian decomposition method (ADM) [15], and very recently Motsa and Sibanda proposed a numerical scheme based on successive linearization method (SLM) [16] for the approximate numerical solutions of singular BVPs of the form (1)- (3).
In last few decades many researchers have employed evolutionary computing based methods for handling nonlinear problems arising in engineering and science and particularly for the numerical solution of nonlinear systems of ordinary differential equations (ODEs). The efficiency and the viability of evolutionary computing (EC) based methods have been demonstrated by several authors [17][18][19][20][21][22]. Although a large number of nonlinear ODEs have been solved using the evolutionary computation based methods, only a few are narrated here. Khan [21] for the numerical solution of nonlinear Riccati differential equation of fractional order. Behrang et al. employed PSO based NN for the solution of nonlinear differential equation arising in fluid flow and heat transfer of vertical cone embedded in porous media [22].
The main focus of this work is to present an approximate numerical technique for the solution of nonlinear singular boundary value problems (1)- (3). The technique is stochastic in nature which is based on the hybrid approach of GA with two local search algorithms such as IPA and ASA. GA, IPA, ASA, and two hybrid schemes combining GA with IPA and ASA here called as GA-IPA and GA-ASA have been employed for the optimization of a fitness function which is the main thrust of the presented method. The efficiency and the viability of the presented method are demonstrated by solving three examples arising in physiology. To prove the accuracy and the validity of our results, comparisons have been carried out with the exact solutions and some conventional numerical solutions given in [2,[12][13][14].
The rest of the paper is arranged as follows. In Section 2 we give a brief description of the proposed methodology and heuristic search algorithms such as GA, IPA, and ASA. In Section 3 we provide numerical results that are followed by the discussion. Finally in Section 4 we give some concluding remarks and future work.

Materials and Method
2.1. Proposed Methodology. We may assume that the solution ( ) and its first and second derivatives ( ) and ( ) can be approximated by a linear combination of some basis functions which can be represented as follows: where ( ) is taken as the log sigmoid function which is given by where , , and are real valued unknown adjustable parameters, and is the number of basis functions.
The approximate numerical solution ( ) represented by (4) is conveniently obtained once the optimum values of the adjustable parameters ( , , and ) are acquired. To attain the optimum values of the adjustable parameters, a problem exclusive fitness function is formulated. This fitness function consists of the sum of mean square error associated with the given ODE ( 1 )and the mean square error related to the initial conditions ( 2 ), which is given as follows: where is the cycle index and 1 and 2 are given by (9) and (10), respectively, where ( ), ( ), and ( ), are given by (4)- (6). The fitness function given by (8) contains unknown adjustable parameters ( , , and ). The minimization of (8) is performed by utilizing the heuristic algorithms GA, IPA, and ASA, and two hybrid schemes such as GA-IPA, and GA-ASA to attain the optimal values of , , and . Consequently the approximate numerical solution of the given problem is straightforward determined from (4).

Brief Description of Heuristic Search Algorithms.
Genetic algorithm (GA) is one of the most popular stochastic global search methods in evolutionary algorithms. GA finds the optimal solution of a problem using the evolutionary ideas of natural selection and genetics. GA starts from a randomly generated population of individuals called chromosome.
The Scientific World Journal 3 Each individual within a population is regarded as a possible solution to the problem. The individuals within a population are evaluated using an exclusive fitness function of the problem at hand. The algorithm evolves population iteratively by mean three primary operations: selection, crossover, and mutation to reach the optimal solution [23].
Interior point algorithm (IPA) is a popular local search method which is widely used in varied optimization problems including linear and nonlinear, convex and nonconvex. The algorithm navigates through the interior feasible region following a central path until it reaches an optimal solution. At each iteration the algorithm applies a direct step also called Newton step or a conjugate gradient step to solve a system of Karush-Kuhn-Tucker (KKT) equations [24,25].
Active set methods are iterative algorithms that solve a sequence of subproblems. The algorithm usually predicts and preserves a set of active and inactive constraints at an optimal solution. Generally active set methods work in two separate phases such as feasibility phase and optimality phase. In the feasibility phase the method attempts to find a feasible point for the constraints, while the objective function is ignored. In the optimality phase the method preserves the feasibility, while it attempts to find an optimal point [26].
The hybridization of GA with local search methods can provide improved performance in many optimization problems [27]. In this work we have used the hybridization of GA with two local search methods such as interior point algorithm (IPA) and active set algorithm (ASA). The GA has been used as global optimization which provides the global optimal solution which is subsequently fed into IPA and ASA for local search fine tuning to achieve the improved results.
The procedural steps of the hybrid schemes GA-IPA and GA-ASA are provided in Pseudocode 1, while the parameter settings for the implementation of these algorithms are given in Table 1.

Results and Discussion
In this section the proposed methodology is implemented on three problems arising in physiology. For the accuracy, efficacy, and viability of the proposed method, comparisons of the results are made with the exact solutions and some conventional methods including modified decomposition method (MDM) combined with B-spline collocation technique [13], B-spline functions [14], finite difference method [12], and nonpolynomial cubic spline method [2].
Example 1. We consider the following special case of (1) which arises in thermal explosions [13,14]: subject to the boundary conditions The exact solution of (11) is given by The approximate numerical solution of (11) with the given boundary conditions (12) using the proposed method is achieved by formulating its fitness function described in Section 2. Assuming the number of basis functions = 10, the fitness function is developed as follows: The fitness function given by (13) is minimized by employing the algorithms GA, IPA, and ASA and two hybrid schemes GA-IPA and GA-ASA for the determination of the optimal values of unknown adjustable parameters ( , , and ). MATLAB 7.6 has been used for the implementation of the algorithms in this study. The parameter settings used for the implementation of the algorithms GA, IPA, GA-IPA, and GA-ASA are given in Table 1. The length of the chromosome, that is, the number of unknown adjustable parameters ( , , and ), are chosen equal to 30. The values of these unknown adjustable parameters are bounded between −15 and +15. This was observed by performing many simulation experiments that, by restricting the values of unknown adjustable parameters to the interval [−15, +15], we get better results.
The algorithms are executed according to the prescribed settings in Table 1. The optimal values of the unknown adjustable parameters corresponding to the minimum fitness are acquired. The optimal values achieved by the hybrid schemes GA-IPA and GA-ASA are given in Table 2. The optimal values of the adjustable parameters from Table 2 are used in (3) and consequently the approximate numerical solution ( ) of Example 1 is obtained. The approximate numerical results obtained by the proposed method are given in Table 3. For the accuracy of numerical results and the potency of our method, we also present the absolute errors and maximum absolute errors by our method in Tables 4  and 5, respectively. The comparisons are made with the exact solutions and the approximate numerical results obtained by other conventional methods including modified decomposition method (MDM) combined with B-spline collocation technique [13] and B-spline functions [14]. It is observed from the comparison of the absolute errors in Table 4 and the maximum absolute errors in Table 5 that the proposed method based on heuristic computing yields the approximate numerical solution of the problem given by (11)-(12) with greater accuracy. Furthermore it is evident from comparison of Table 4 that the errors relative to the exact solutions by the proposed heuristic hybrid schemes are much smaller than the errors by the approach I method used in [13], whereas they are relatively smaller than approach II method errors given in [13]. Moreover it is observed from the comparison of Table 5 that the maximum errors by the proposed method using hybrid schemes are comparable to B-spline method given in [14]. 4 The Scientific World Journal Step 1. Population Initialization A population of individuals or chromosome is generated using random number generator. The length of the chromosome represents the number of unknown adjustable parameters to be optimized.

Step 2. Fitness Evaluation
A problem relevant fitness function is used to compute the fitness of each individual in the current population.

Step 3. Stoppage Criteria
The algorithm stops if the maximum number of generations (cycles) has exceeded or a certain level of fitness value has reached. If the stopping criterion is fulfilled then go to step 6 for local search fine tuning, else continue and repeat steps 2 to 5.

Step 4. Selection and Reproduction
The chromosomes from the current population are selected on the basis of their fitness which acts as parents for new generation. These parents produce children (offsprings) with a probability to their fitness through crossover operation.

Step 5. Mutation
Mutation operation introduces random alterations in the genes to retain the genetic diversity to find a good solution.

Step 6. Local Search Fine Tuning
The optimal chromosome achieved by the GA is fed to IPA for fine tuning and improvement.
Pseudocode 1: Hybridization of GA with IPA and ASA.
Example 2. We consider the following nonlinear singular boundary value problem [2]: with the boundary conditions To obtain the approximate numerical solution of (16) subject the boundary conditions (17) using the proposed method, its fitness function with = 10 is developed as follows: The Scientific World Journal 5   The Scientific World Journal Table 5: Comparison of maximum absolute error for Example 1 between the proposed method and the methods given in [13,14].
Proposed method Method in [13] Approach I Method in [13] Approach II B-spline method [14] 3 The heuristic algorithms GA, IPA, GA-IPA, and GA-ASA are executed with the same parameter settings given in Table 1 for the minimization of (20). To prove the effectiveness and the viability of the proposed method we have obtained the approximate numerical solutions of (16)-(17) for various 8 The Scientific World Journal  values of the parameter (0.25, 1, 2, and 8). The optimal values attained by the hybrid schemes corresponding to the minimum fitness are provided in Table 6 for = 0.25, while for the rest of the values of , these optimal values of adjustable parameters have been omitted here. For comparison maximum absolute errors have been computed corresponding to all the specified values of (0.25, 1, 2, and 8). The comparison of the maximum absolute errors between the proposed heuristic method and the standard numerical methods such as finite difference method [12] and nonpolynomial cubic spline method [2] are presented in Table 7 for = 0.25 and = 1 and in Table 8 for = 2 and = 8, respectively. The comparison noticeably reveals the potency and the accuracy of the proposed method. The comparison also reveals the improved performance of hybrid schemes GA-IPA and GA-ASA.
Example 3. We consider Example 2 again with a change in boundary condition as follows [12]: The fitness function of this example is given as follows: The fitness functions given by (24) is subject to minimization using the heuristic algorithms GA, IPA, and ASA and two hybrid schemes GA-IPA and GA-ASA with the same parameter settings prescribed in Table 1 for the determination of the unknown adjustable parameters. The optimal values of  In Tables 11 and 12 our results are compared with the exact solutions when = 0.25 and = 0.75. Moreover in Table 13 we also present a comparison of the maximum absolute errors between our method and the finite difference method given in [12]. It is observed that the proposed method yields the approximate solutions fairly comparable to the finite difference method given in [12].

Conclusions and Future Work
In this study a hybrid heuristic computational approach has been successfully implemented for the approximate numerical solution of nonlinear singular boundary value problems (BVPs) arising in physiology. It can be concluded on the basis of the comparisons of the results made with the exact solutions and some of the standard approximate numerical solutions that the proposed method possesses a great potential and viability for solving nonlinear singular boundary value problems (BVPs) arising in diverse fields of engineering and science. The strength of proposed method has been illustrated by solving three nonlinear problems appearing in physiology. Moreover the proposed methodology can provide the approximate numerical solution straightforward and on a continuous grid of time once the optimal values of the unknown adjustable parameters are attained.
In future we intend to employ the proposed methodology to other such nonlinear singular boundary value problems, nonlinear ordinary differential equations (ODEs), and nonlinear coupled ODEs arising in various fields of engineering and applied science. We also seek to use other evolutionary algorithms and different basis functions for the approximate solutions of such problems.