Application of Differential Evolution in the Solution of Stiff System of Ordinary Differential Equations

In recent times, the adaptation of evolutionary optimization algorithms for obtaining optimal solutions of many classical problems is gaining popularity. To accurately solve stiff systems in Ordinary Differential Equation (ODEs), the use of more accurate numerical techniques is of great importance. In this paper, optimal approximate solutions of initial--valued stiff system of ODEs are obtained by converting the ODE into constrained optimization problem. The latter is then solved via differential evolution algorithm. To illustrate the efficiency of the proposed approach, two numerical examples were considered. This approach showed significant improvement on the accuracy of the results produced compared with existing methods discussed in literature.


Introduction
Initial value problems involving systems of first-order ordinary differential equation are often encountered in many fields. Most of these systems exhibit a phenomenon known as stiffness. Definition 1.1 The linear system ̇= + ( ) is said to be stiff if R < 0 , = 1,2, ⋯ , , and max =1,2,⋯, |R | >> min =1,2,⋯, |R |, where , = 1,2, ⋯ , are the eigenvalues of . The ratio is called the stiffness ratio. [15] � max =1,2,⋯, Such systems are always very difficult to solve when using classical techniques [9]. The difficulty of solving stiff initial value problems was clearly identified in early 1950s when the authors in [9] published one of the first papers in which the problem of stiffness was stated. Subsequently, a whole lot of methods and algorithms have been proposed for solving problems that exhibit stiffness. The author in [6] introduced a class of extended backward differentiation formulae suitable for the integration of stiff systems of autonomous initial value problems. In a later work, [8], the author proposed classes of predictor-corrector method involving the second derivative using the extended backward differentiation formulae. Obtaining the solution of stiff differential equations using exponentially fitted formula was first proposed by the authors in [16]. Following the work of the authors in [16], several method based on the concept of exponential fitting have been proposed by various authors [1,3,6,7,8,10,12,19,20].
In recent years however, interest in the application of evolutionary algorithm for solving many problems is on the increase. Evolutionary algorithm has also been adopted in solving ODE and its application to real life problems is on the increase. Problems in ordinary differential equations are converted to optimization problems and then solved using optimization techniques. In this direction, the classical genetic algorithm was used to obtain approximate solutions of second-order initial value problems by the authors in [11]. The author in [18] computed approximate solutions of first-order initial value problem via the combination of collocation method (finite elements) and genetic algorithms. The combination of genetic algorithm with the Nelder-Mead method for solving the second-order initial value problem of the form ′′ = ( , ) was proposed in [17]. The use of neural network for the approximate solution of secondorder initial value problems was proposed by the authors in [14]. In the work of the authors in [21], the use of continuous genetic algorithm in obtaining the solution of the two-point second-order ordinary differential equation was presented. The use of differential evolution algorithm in obtaining approximate solution of second-order initial value problem of the form ′′ + ( ) ′ + ( ) = ( ) was proposed by the authors in [5]. In [23], approximate solutions of the second-order two-point boundary value problem ′′ = ( , ); ( ) = 1 ; ( ) = 2 with oscillatory/periodic behaviour was obtained via the application of differential evolution algorithm. In this work, we show that the problem of finding approximate solution of a stiff system of first-order ordinary differential equation can be converted to an optimization problem. The resulting optimization problem is then solved using differential evolution. It is assumed that the solution can be approximated by a linear combination of exponential terms. The differential evolution algorithm is used to optimize the coefficients of the terms of the solutions.

Differential Evolution Algorithm
Differential evolution is a simple stochastic function minimizer which details can be found in many standard texts. However, we give a brief overview of the algorithm as described in [22]. Formally, let : ℝ → ℝ be the function to be optimized. The function takes a candidate solution as argument in the form of a vector of real numbers and produces a real number as output which indicates the fitness of the given candidate solution. The gradient of f is not known. The goal is to find a solution m for which ( ) ≤ ( ) for all p in the search-space, which would mean m is the global minimum. Maximization can be performed by considering the function ℎ: = − instead.
Let ∈ ℝ designate a candidate solution (agent) in the population. The basic differential evolution algorithm can then be described as follows: Initialize all agents with random positions in the search-space.
Until a termination criterion is met (e.g. number of iterations performed, or adequate fitness reached), repeat the following: For each agent in the population do: Pick three agents , , and from the population at random, they must be distinct from each other as well as from agent If ( ) < ( ) then replace the agent in the population with the improved candidate solution, that is, replace with in the population.
Pick the agent from the population that has the highest fitness or lowest cost and return it as the best found candidate solution.
Note that ∈ [0,2] is called the differential weight and CR ∈ [0,1] is called the crossover probability, both these parameters are selectable by the practitioner along with the population size NP ≥ 4.

Procedure for Proposed Technique
In this work, we considered the system of r equations = ( , 1 , 2 , ⋯ , ), . It is assumed that the solution ( ) of (2) can be written as where , , , , = 1,2 ⋯ , are real constants which values are to be determined by the proposed procedure.
Substituting the derivative of (3) into (2) gives Using the initial conditions given in (2) we have, From (4), the error expression ℰ , ( ) at each node point is obtained as Now we need to find values of , , , , = 1,2 ⋯ , , that minimizes the expression where ℎ is the steplength.
Using (7) and (5), the problem of finding approximate solution of (2) is now formulated as an optimization problem in the following manner: In this work, we solved (8) using the differential evolution algorithm. Suitable values of , , , , = 1,2 ⋯ , which satisfies ((8)) were obtained using this technique. We shall refer to this scheme as DE4stiffODE.

Numerical Examples
In this section, the efficiency of the proposed method is investigated on two test problems. The results obtained by our approach are compared with those obtained with the methods proposed by authors in the literature.

Example 1
The two dimensional nonlinear Kaps problem with exact solution is considered. Using our proposed technique, the problem is solved over the interval [0,1] with a steplength ℎ = 0.1. As expected, our method approximated the solution almost exactly as A comparison of our results with those obtained by the methods of authors in [4,13,24] is presented in Table 1. We shall respectively denote by "SDM10" and "SDM14" the methods of order 10 and order 14 proposed in [24]. Similarly, the method proposed in [4] shall be denoted by "NMTD". The 8th-order method proposed in [13] shall be denoted by "SDAM8". Table 1. The absolute error of our "DE4stiffODE" method compared with "SDM10", "SDM14" and "NMTD" methods at some values of t on Example 1. It is clear from Table 1 and Table 2 that our proposed approach gave better solution compared with existing methods in literature.

Example 2
The second problem considered is the following system of two linear equations Problem (13) has been studied by authors in [2,4,6,7,10,12]. Using a steplength of ℎ = 0.7, we solve the problem in the interval [0,1] to obtain the following solution: Again, as expected, our approach approximated the solution exactly. For comparison purpose, we shall denote respectively by "ABOT", "NMTD", "Cash4", "Cash5", "SDEBDF", and "J-K" the methods proposed by the authors in [2], [4], [6], [7], [10], and [12]. The efficiency of our proposed scheme can be seen from the choice of the steplength. Even with a larger steplength, our scheme gave a better result compared with the methods discussed in the literature.
7 Table 3. Comparison of absolute error of Example 2 at t = 1 with methods discussed in literature

Conclusion
The basic steps involved in converting a stiff system of differential equation problem to an optimization problem and obtaining approximate solution via differential evolution have been shown. Results of numerical implementation show that our approach can be used to obtain almost-exact solution. The comparison of our method with those discussed in literature clearly shows the degree of accuracy of our method.