Closed-Form Solutions to Differential Equations via Differential Evolution

We focus on solving ordinary differential equations using the evolutionary algorithm known as differential evolution (DE). The main purpose is to obtain a closed-form solution to differential equations. To solve the problem at hand, three steps are proposed. First, the problem is stated as an optimization problem where the independent variables are elementary functions. Second, as the domain of DE is real numbers, we propose a grammar that assigns numbers to functions.Third, to avoid truncation and subtractive cancellation errors, to increase the efficiency of the calculation of derivatives, the dual numbers are used to obtain derivatives of functions. Some examples validating the effectiveness and efficiency of our method are presented.


Introduction
Most of the problems in engineering and physics can be modeled as ordinary differential equations (ODEs).For this reason there are many studies addressing their solution.Regarding the deterministic arena, the most used methods are the Runge-Kutta methods [1][2][3][4], predictor-corrector methods [5][6][7], and radial basis functions methods [8][9][10].Recently, some studies dedicated to solve differential equations using nondeterministic methods have been published.In [11], genetic algorithms are used to solve some differential equations appearing in economic sciences.In [12] a variational approach has been used in order to solve elliptic partial differential equations, and a genetic algorithm is used as the optimization method.In all the previously referenced articles-deterministic or not-the solution is given in a numerical approximated form.There are very few studies reporting closed-form solutions to differential equations.For example, using the evolutionary method known as grammatical evolution, [13] reports a method that produces closedform solutions.Another approach that produces closed-form solutions to differential equations is [14], and the method used there is a hybrid method combining grammatical evolution and neural networks.
In this paper we propose a method based on the DE algorithm that obtains solutions to second-order ODEs as closedform expressions.When the exact solution is not reached, the algorithm we propose converges to a solution that minimizes the objective functional of an optimization problem considering both, the differential equation and the boundary conditions.As DE was proposed to minimize real valued functions (i.e., not functionals), we propose a one-to-one grammar that assigns integer numbers to elementary functions; in this way, we are able to use DE to minimize the objective functional, which measures if a candidate function (a candidate function is the result of the evolution process obeying the DE algorithm) satisfies or not the ODE we want to solve.
In order to compute the first and second derivatives of the candidate function, we use the dual number approach.In this way, the derivatives are directly obtained without the use of a limit process, thus avoiding truncation and subtractive cancellation errors.All the programming functions required to solve an ODE are coded in Fortran language and they are included in a folder, which is provided as additional material to this paper, whose download link is as follows: http://www .meca.cinvestav.mx/personal/cacruz/archivos-ccv/.
The rest of the paper is organized as follows.Section 2 states the optimization problem and describes the classical DE algorithm.Section 3 presents the proposal of the oneto-one grammar which allows using DE for minimization of functionals.Section 4 presents the dual number approach to obtain the first and second derivatives of the candidate functions.Section 5 works out several examples and applications of our method.Conclusions are presented in Section 6.Finally, two appendixes close the paper.Appendix A presents the way in which the candidate function is generated and evaluated.Appendix B presents graphs of the behavior of the DE algorithm for each worked-out example.

Statement of the Problem
Let us consider a second-order ordinary differential equation (1) defined on the real interval [, ], with boundary conditions ( 2) and ( 3), where  1 <  2 and  ∈  2 .Note that it is not required that the functions  : R 4 → R, h 1 : R 3 → R 2 , and h 2 : R 3 → R 2 be differentiable: The problem that we address in this paper is to find a closedform expression for () satisfying ( 1), (2), and (3).Therefore, we rewrite the problem as that of minimizing the functional (4) under   (), where  1 > 0 and  2 > 0 are weighting factors (chosen by the user): If there exists a function   () for which  = 0, then   () will be the solution () to (1), satisfying (2) and (3).
As the approach we follow to minimize (4) is evolutive, we use the differential evolution algorithm which is a simple yet powerful evolutionary algorithm for global optimization introduced by Storn and Price [15], which is presented below.

Differential Evolution.
The DE algorithm has gradually become more popular and has been used in many practical cases.It only requires information about the objective function itself, which does not need to be a differentiable function, and the state space of possible solutions can be disjoint and can encompass infeasible regions [16].Below, the original version of the method-known as DE/rand/1/bin-is outlined [17].
(2) Initialization of population is as follows: Vectors b  and b  are the parameter limits and rand  (0, 1) is a random number in [0, 1) generated for each parameter.
(3) Mutation is as follows: x  0 ; is called the base vector which is perturbed by the difference of two other vectors.
is a scale factor greater than zero.
A dual recombination of vectors is used to generate the trial vector: The crossover probability, Cr ∈ [0, 1], is a userdefined value,  rand ∈ [1, ].
The selection is made according to x ; otherwise.
In our study we use the DE/rand/1/bin method with the dither variant, which means that the parameter  is taken to be a random number-in our case  ∈ (0, 1).

Construction of Functions
As we can see, the DE method was designed to seek the optimum individual x on a real continuum domain.Since we are interested in solving differential equations (i.e., minimizing a functional), our individuals are functions.Therefore, we associate a function with a real vector; once this is done, the DE method can be applied as usual.
In order to assign a function to a real vector we propose the grammar shown in Table 1.The relation between a set of numbers and a function can be done by using a parse tree.This parse tree should be read from top to down and from left to right.Table 2 shows an example of a function construction and Figure 1 shows the corresponding parse tree.We can get an easy understanding of the function construction by conceptualizing the operators +, * , /, as functions of two arguments.For example, the expression 1 + 2 can be seen as plus (1,2) where the plus function is defined as plus(, ) =  + .
The set of integers related to a mathematical function can be manipulated by the DE method but the mutation operator will produce a set of real numbers that will not necessarily be a set of integers.This is addressed by taking the integer part of the numbers or by using the floor (ceiling) function.

Evaluation of the Constructed Function and Its Derivatives
For the evaluation of the constructed function we have written a Fortran parse function that receives the integer vector generated by the proposed grammar and produces a candidate function and its derivatives (first and second) evaluated at some specified point.The construction of this programming function is explained in Appendix A. Traditional methods for calculating numerical derivatives (finite-difference) are subject to both truncation and subtractive cancellation errors.These problems are avoided by using dual functions.The approach to obtain first order derivatives by using dual functions is well known [18][19][20].However, in order to make the paper self-contained this section presents the essential ideas as follows [20].
A dual number is a number of the form x =  +  where ,  ∈ R, the field of the real numbers, and  2 = 0. From the Taylor theorem if a function () is analytic, then a substitution of  + ℎ →  +  in the above formula will give The function ( + ) is called a dual function ŷ of the dual variable x = +.So if we substitute all of our real numbers by dual numbers and make the coefficient  of the dual variable equal to one, we end up with a dual function whose real term is the original function and the dual term is its derivative.Another convenient notation for the function ŷ is where  0 = () and  1 =   ().
Applying the chain rule we can dualize the composition of () with the function (): From this, a generalization to second derivatives is straightforward.Using a tilde to denote such a generalization we have where  0 = (),  1 =   (), and  2 =   ().Similarly, for the composition (()), we will have

Experimental Results
This section presents the results that are obtained when we apply our proposal to obtain closed-form solutions to the case studies considered in [13].All the cases with closed-form solutions were reproduced.Below, we present our results for some interesting cases.The experiments were performed on an Intel Core i5-3230M @ 2.60 GHz processor running Debian GNU/Linux and using the Intel Fortran Compiler.In all the cases we used 100 equally spaced points to evaluate (4),  1 =  2 = 10, and a crossover probability of 0.2.Case 1.In this case, we want to find a closed-form solution to the differential equation ( 16), subject to boundary conditions (0) = 0 and   (0) = 1, with  ∈ [0, 1]: Since for the present case lim it is not difficult to prove that ( 16) is equivalent to Even when both equations are equivalent, ( 18) is more adequate to the minimization of (4).The exact solution to ( 16), ( 18) is given by ( 19) and the approximated closed-form solution we have found is given by (20): Since the exact solution is known we can quantify the error in the interval  ∈ [0, 1] as (21), which for this case resulted as  2 = 2.83 × 10 −7 : In Figure 2, we show the exact and approximated solutions in the interval [0, 2.5], since the error in the interval [0, 1] is not enough to recognize any difference between both solutions.For this case we used a population of 100 individuals of dimension 20 and 5 000 generations.Case 2. Let us consider the ODE with  ∈ [0, 1] and initial conditions (0) = 0 and   (0) = 1.For this case a closed-form solution is not known.The approximate solution we found is giving a value of 1.76 × 10 −4 for (4).In order to obtain (23) we used a population of 80 individuals of dimension 20 and 15 000 generations.In Figure 3 we show the approximate solution we found.
Case 3. Let us consider the ODE studied in [21]: with  ∈ [0, 1] and initial condition (0) = 0.The proposed method is able to find the exact solution () =  + tanh , when the population size is 100 individuals of dimension 30 and using 100 generations.Clearly the final value for objective functional (4) is zero.
In Figure 4 we show the approximate solution we found.Case 5. Let us consider the temperature distribution equation on a uniformly thick rectangular fin radiation to free space studied in [23,24]: with  ∈ [0, 1] and boundary conditions (1) = 1 and   (0) = 0.
For this case, the approximate solution that the proposed methodology has obtained is   () = sin () tan −1 (0.25 tan ()) + 0.693147, when we use 300 individuals of dimension 30 and 10 000 generations.A final value of 9 × 10 −4 was obtained for the objective functional (4).In Figure 5 we show the approximate solution we found.
Case 6.Consider the differential equation modeling the cooling process of a lumped system by combined convection and radiation [25]: with  ∈ [0, 1] and initial condition (0) = 1.
Applying the proposed methodology we obtained the approximate solution For this case the value of (4) was 2 × 10 −4 , when using 500 individuals of dimension 30 and 5 000 generations.
In Figure 6 we show the approximate solution we found.
Case 7. Consider the Duffing equation studied in [26]: but with boundary conditions (0.5) = 0, (1) = 1.The approximate solution we found is giving a value of 5.2 × 10 −4 for (4).In order to obtain (32) we used a population of 500 individuals of dimension 30 and 15 000 generations.In Figure 7 we show the approximate solution we found.
Case 8. We close this section considering three examples of the Van der Pol equation: studied in [27].Let us consider  = 1,  = −0.05, = 0.05,  = −1, and  = 0, so we have By applying the proposed methodology to this equation and considering boundary conditions (0) = 0 and   (0) = 0.5, we obtain For this case the final value of the objective functional (4) is 3 × 10 −4 , when using 300 individuals of dimension 30 and 5 000 generations.In Figure 8 we show the approximate solution we have found.Considering boundary conditions (0) = 0.5 and   (0) = 1 and applying the proposed methodology, we obtain the approximate solution + arctan (cos (cosh cosh 0.205087 () For this case, the final value for the objective functional (4) is 7 × 10 −3 , when using a population of 150 individuals of dimension 15 and 70 000 generations.In Figure 9 we show the approximate solution we have found.
Case 10.Taking  = 1000, for (33) and considering the same boundary conditions and values for the other parameters as in the previous case, we have By applying the proposed methodology to solve this equation, we obtain the approximate solution For this case the value of ( 4) is 2 × 10 −2 , when using 300 individuals of dimension 30 and 20 000 generations. Figure 10 shows the found approximate solution.
As a final remark, we want to point out that all the approximate solutions shown in this section were practically the same to those obtained by applying the Runge-Kutta (RK) algorithm.It is well known that the RK algorithm is by far the most used and in our opinion the best choice, for numerically solving a differential equation.Nevertheless, the proposed evolutionary approach is useful for cases when a closed-form solution is needed.Moreover, the introduced algorithms with dual numbers open a door for many applications.

Conclusions
We have demonstrated the use of the differential evolution algorithm in order to obtain closed-form solutions to secondorder differential equations.The main drawback for the application of the differential evolution algorithm to find exact closed-form solutions to differential equations was handled by constructing a function associated with a real vector.This function was constructed using a simple grammar and the concept of parse tree.This evolutionary algorithm along with the use of dual numbers in order to obtain the derivatives of a function produces an approximate solution expressed as a closed-form expression, even if the exact solution cannot be found (or it is not known as a closed-form expression).Thus, the proposed method could be efficient and useful for practical applications.All the programming functions needed for the solution of the differential equations were coded in Fortran language and provided as additional material to this paper.

A. Parse Function
The function that maps the integer vector x associated with a function and evaluates it in the dual number xval is the function parsv(x,xval) and is coded in the module parsef mod.f90 of the additional material.Below we describe the construction of this function.
In what follows, NaN stands for not a number, ñan stands for a dual number with, at least, its real component a NaN.The dual version of a number  is written as ã; for instance, the dual version of the number 2 is written as 2 and in our work is treated as the vector [2, 0, 0].
The parsv function requires the following functions: (i) mynan() (ii) narg(x) Function that calculates the number of arguments for the elements of the grammar.For example, narg The Fortran code for the parsv function is shown in Algorithm 1.
For the sake of concreteness, we will exemplify taking into account only the real component.The dual components will be omitted; notice however that the parsv function coded above takes into account such dual components.
The instruction auxv = mynan() sets auxv to auxv = [N,N,N,N,N,N,N,N,N,N] where N stands for NaN.In fact auxv is a 3 × 10 matrix but as we said before, we are ignoring the dual components.