Abstract
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.
1. 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–4], predictor-corrector methods [5–7], and radial basis functions methods [8–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 closed-form 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 closed-form 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 one-to-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.
2. 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 and . Note that it is not required that the functions , , and be differentiable:The problem that we address in this paper is to find a closed-form expression for satisfying (1), (2), and (3). Therefore, we rewrite the problem as that of minimizing the functional (4) under , where and are weighting factors (chosen by the user):If there exists a function for which , 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.
2.1. 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].(1)The population is described by where , , and represent the dimensionality of , the number of individuals, and the number of generations, respectively.(2)Initialization of population is as follows: Vectors and are the parameter limits and is a random number in generated for each parameter.(3)Mutation is as follows: is called the base vector which is perturbed by the difference of two other vectors. , . is a scale factor greater than zero.(4)Crossover is as follows. A dual recombination of vectors is used to generate the trial vector: The crossover probability, , is a user-defined value, .(5)Selection is as follows. The selection is made according to 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 .
3. Construction of Functions
As we can see, the DE method was designed to seek the optimum individual 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 can be seen as where the function is defined as .
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.
4. 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–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 where , the field of the real numbers, and . 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 . 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 and .
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 , , and . Similarly, for the composition , we will have
5. 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), , 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 and , with :Since for the present case 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 as (21), which for this case resulted as : In Figure 2, we show the exact and approximated solutions in the interval , since the error in the interval 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 ODEwith and initial conditions and . For this case a closed-form solution is not known. The approximate solution we found is giving a value of 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 and initial condition . The proposed method is able to find the exact solution , when the population size is 100 individuals of dimension 30 and using 100 generations. Clearly the final value for objective functional (4) is zero.
Case 4. Let us consider the nonlinear differential equation [22]: with and boundary conditions and . The approximate solution we obtained with a value of for (4) using 100 individuals of dimension 30 and 10 000 generations was 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 and boundary conditions and .
For this case, the approximate solution that the proposed methodology has obtained is when we use 300 individuals of dimension 30 and 10 000 generations. A final value of 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 and initial condition .
Applying the proposed methodology we obtained the approximate solution For this case the value of (4) was , 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 , . The approximate solution we found is giving a value of 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 , , , , and , so we have By applying the proposed methodology to this equation and considering boundary conditions and , we obtain For this case the final value of the objective functional (4) is , when using 300 individuals of dimension 30 and 5 000 generations. In Figure 8 we show the approximate solution we have found.
Case 9. If we now consider , , , , and for (33), we obtain Considering boundary conditions and and applying the proposed methodology, we obtain the approximate solutionFor this case, the final value for the objective functional (4) is , 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 , 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 , 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.
6. Conclusions
We have demonstrated the use of the differential evolution algorithm in order to obtain closed-form solutions to second-order 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 the Fortran language and provided as additional material to this paper.
Appendices
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, 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 and in our work is treated as the vector .
The parsv function requires the following functions:(i) mynan A NaN implementation.(ii) narg(x) Function that calculates the number of arguments for the elements of the grammar. For example, narg() = 2.(iii) operf(x,xval) Function that evaluates the functions of one argument in the dual point xval and then changes the real component of xval to NaN. For instance, operf(17,xval) will return sindual(xval) and the xval number is changed to .(iv) operf2(x,vecval) Function that evaluates the functions of two arguments, that is, +, , /, power. vecval is an array with dual numbers, possibly with NaN components. This function operates on the first two dual no- components of the vecval array and then changes its real components to NaN components. For instance, operf2(12,[,sindual(), ,]) will return sindual() + and then vecval = [, sindual(),,] is changed to [, , , ].The Fortran code for the parsv function is shown in Algorithm 1.
| ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Let us analyze the above code for the case when x = [,,,,,,,,,] and xval = . In this particular case, the dimension length of the vector x is length = 10 and x represents the function .
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 matrix but as we said before, we are ignoring the dual components.
The first loop evaluates the functions of zero arguments changing auxv to auxv = [N,N,N,2,1.1,N,N,N,1.1,5] as shown in Algorithm 2.
| ||||||||
The second loop is a little more complicated but essentially evaluates the functions of one argument and two arguments. Let us analyze it in detail as shown in Algorithm 3.
| ||||||||||||||||||||||
For = 9 the if condition is satisfied; thus the auxv vector is not changed.
For = 8 the second elseif condition is satisfied and the auxv vector turns out to be
auxv = [N,N,N,2,1.1,N,N,1.1/5,N,N], and so forth.
Finally for = 1 we have
auxv = [ + ,N,N,N,N,N,N,N,N,N].
For the case of taking into account the dual components, auxv(1,1) = , auxv(2,1) = , and auxv(3,1) = which means that parsv(1.1,xval)= .
Notice that if the vector x does not represent a valid function, the function parsv will look from right to left for a “subset” of trying to construct a valid function; if it succeeds, the parsv function will return the value corresponding to that part of corresponding to a valid function. For example, if x = [,,,,,,,,,,], we cannot form a valid function. However, the parsv function will return a right value corresponding to xs = [,,,,]. Notice however that as long as the DE method is concern, is always treated as [,,,,,,,,,,]. This behavior for the function parsv could appear undesirable to first look. Nevertheless it is a nice property if we are thinking to use such a function in connection with an evolutionary algorithm. The reason is that if we discard those vectors that do not correspond to a valid function, the evolutionary method will have not enough individuals to evolve. Of course at the end we need to extract the subset of x which correspond to a valid function; this is done by the function filterx of the parsef_mod module.
B. Performance Graphs
This appendix shows the performance graphs of the DE method for all the studied cases as shown in Figures 11, 12, 13, 14, 15, 16, 17, 18, 19, and 20.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.