# Closed-Form Solutions to Differential Equations via Differential Evolution

Authors: L. Mex; Carlos A. Cruz-Villar; F. Peñuñuri
Journal: Discrete Dynamics in Nature and Society
 (2015)

---

## 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 areelementary 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.


---

## Body

## 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. Section2 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 [
a
,
b
], with boundary conditions (2) and (3), where x
1
<
x
2 and y
∈
C
2. Note that it is not required that the functions f
:
R
4
→
R, h
1
:
R
3
→
R
2, and h
2
:
R
3
→
R
2 be differentiable:(1)
f
(
x
,
y
(
x
)
,
y
′
(
x
)
,
y
′′
(
x
)
)
=
0
,
(2)
h
1
x
1
,
y
x
1
,
y
′
x
1
=
0
,
(3)
h
2
x
2
,
y
x
2
,
y
′
x
2
=
0
.The problem that we address in this paper is to find a closed-form expression for y
(
x
) satisfying (1), (2), and (3). Therefore, we rewrite the problem as that of minimizing the functional (4) under y
s
(
x
), where λ
1
>
0 and λ
2
>
0 are weighting factors (chosen by the user):(4)
R
=
λ
1
h
1
x
1
,
y
s
x
1
,
y
s
′
x
1
2
+
λ
2
h
2
x
2
,
y
s
x
2
,
y
s
′
x
2
2
+
∫
x
1
x
2
f
2
x
,
y
s
x
,
y
s
′
x
,
y
s
′′
x
d
x
.If there exists a function y
s
(
x
) for which R
=
0, then y
s
(
x
) will be the solution y
(
x
) 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(5)
P
x
,
g
=
x
i
,
g
,
i
=
1
,
…
,
m
;
g
=
0
,
…
,
g
max
⁡
,
x
i
;
g
=
x
i
;
g
j
,
j
=
1
,
…
,
D
,
whereD, m, and g
max
⁡ represent the dimensionality of x, the number of individuals, and the number of generations, respectively.
(2)
Initialization of population is as follows:(6)
x
i
;
0
j
=
rand
j
(
0,1
)
·
(
b
U
j
-
b
L
j
)
+
b
L
j
.
Vectorsb
U and b
L are the parameter limits and rand
j
(
0,1
) is a random number in [
0,1
) generated for each parameter.
(3)
Mutation is as follows:(7)
v
i
;
g
=
x
r
0
;
g
+
F
·
(
x
r
1
;
g
-
x
r
2
;
g
)
.
x
r
0
;
g is called the base vector which is perturbed by the difference of two other vectors.
r
0
,
r
1
,
r
2
∈
{
1,2
,
…
,
m
}, r
1
≠
r
2
≠
r
3
≠
i. F is a scale factor greater than zero.
(4)
Crossover is as follows.
A dual recombination of vectors is used to generate the trial vector:(8)
u
i
;
g
=
u
i
;
g
j
=
v
i
;
g
j
if
rand
j
0,1
⩽
Cr
or
j
=
j
rand
x
i
;
g
j
otherwise
.
The crossover probability,Cr
∈
[
0,1
], is a user-defined value, j
rand
∈
[
1
,
D
].
(5)
Selection is as follows.
The selection is made according to(9)
x
i
;
g
+
1
=
u
i
;
g
if
f
(
u
i
;
g
)
⩽
f
(
x
i
;
g
)
x
i
;
g
otherwise
.In our study we use the DE/rand/1/bin method with thedither variant, which means that the parameter F is taken to be a random number—in our case F
∈
(
0,1
).


## 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(5)
P
x
,
g
=
x
i
,
g
,
i
=
1
,
…
,
m
;
g
=
0
,
…
,
g
max
⁡
,
x
i
;
g
=
x
i
;
g
j
,
j
=
1
,
…
,
D
,
whereD, m, and g
max
⁡ represent the dimensionality of x, the number of individuals, and the number of generations, respectively.
(2)
Initialization of population is as follows:(6)
x
i
;
0
j
=
rand
j
(
0,1
)
·
(
b
U
j
-
b
L
j
)
+
b
L
j
.
Vectorsb
U and b
L are the parameter limits and rand
j
(
0,1
) is a random number in [
0,1
) generated for each parameter.
(3)
Mutation is as follows:(7)
v
i
;
g
=
x
r
0
;
g
+
F
·
(
x
r
1
;
g
-
x
r
2
;
g
)
.
x
r
0
;
g is called the base vector which is perturbed by the difference of two other vectors.
r
0
,
r
1
,
r
2
∈
{
1,2
,
…
,
m
}, r
1
≠
r
2
≠
r
3
≠
i. F is a scale factor greater than zero.
(4)
Crossover is as follows.
A dual recombination of vectors is used to generate the trial vector:(8)
u
i
;
g
=
u
i
;
g
j
=
v
i
;
g
j
if
rand
j
0,1
⩽
Cr
or
j
=
j
rand
x
i
;
g
j
otherwise
.
The crossover probability,Cr
∈
[
0,1
], is a user-defined value, j
rand
∈
[
1
,
D
].
(5)
Selection is as follows.
The selection is made according to(9)
x
i
;
g
+
1
=
u
i
;
g
if
f
(
u
i
;
g
)
⩽
f
(
x
i
;
g
)
x
i
;
g
otherwise
.In our study we use the DE/rand/1/bin method with thedither variant, which means that the parameter F is taken to be a random number—in our case F
∈
(
0,1
).


## 3. Construction of Functions

As we can see, the DE method was designed to seek the optimum individualx 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 Table1. 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
(
x
,
y
)
=
x
+
y.Table 1
Proposed grammar.
String
Associated number
Integer from one to ten
1 : 10
x
11
+
12
*
13
/
14
pow()
15
-
(
)
16
sin
⁡
⁡
(
)
17
cos
⁡
(
)
18
exp
⁡
(
)
19
log
⁡
(
)
20
(
)
21
tan
⁡
(
)
22
arcsin
⁡
(
)
23
arccos
⁡
(
)
24
arctan
⁡
(
)
25
erf
⁡
(
)
26
sinh
⁡
(
)
27
cosh
⁡
(
)
28
tanh
⁡
(
)
29Table 2
Construction of the functionf
(
x
)
=
sin
⁡
⁡
x
+
cos
⁡
x.
Chromosome
String
12:
+
12, 17:
sin
⁡
⁡
(
)
+
12, 17, 11:
sin
⁡
⁡
(
x
)
+
12, 17, 11, 18:
sin
⁡
⁡
(
x
)
+
cos
⁡
(
)
12, 17, 11, 18, 11:
sin
⁡
⁡
(
x
)
+
cos
⁡
(
x
)Figure 1
Parse tree for the functionf
(
x
)
=
sin
x
+
cos
⁡
x.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 AppendixA.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 formx
^
=
α
+
ϵ
β where α
,
β
∈
R, the field of the real numbers, and ϵ
2
=
0. From the Taylor theorem if a function y
(
x
) is analytic, then (10)
y
x
+
h
=
y
x
+
y
′
x
h
+
y
′′
x
2
h
2
+
⋯
;a substitution of x
+
h
→
x
+
ϵ in the above formula will give (11)
y
(
x
+
ϵ
)
=
y
(
x
)
+
y
′
(
x
)
ϵ
.The function y
(
x
+
ϵ
) is called a dual function y
^ of the dual variable x
^
=
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 y
^ is (12)
y
^
(
x
^
)
=
{
y
0
,
y
1
}
,where y
0
=
y
(
x
) and y
1
=
y
′
(
x
).Applying the chain rule we can dualize the composition ofy
(
x
) with the function u
(
x
): (13)
y
^
(
u
^
)
=
{
y
0
(
u
0
)
,
y
1
(
u
0
)
u
1
}
.From this, a generalization to second derivatives is straightforward. Using a tilde to denote such a generalization we have (14)
x
~
=
x
,
1,0
,
y
~
(
x
~
)
=
{
y
0
,
y
1
,
y
2
}
,where y
0
=
y
(
x
), y
1
=
y
′
(
x
), and y
2
=
y
′′
(
x
). Similarly, for the composition y
(
u
(
x
)
), we will have (15)
y
~
(
u
~
)
=
{
y
0
(
u
0
)
,
y
1
(
u
0
)
u
1
,
y
2
(
u
0
)
u
1
2
+
y
1
(
u
0
)
u
2
}
.


## 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), λ
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 y
(
0
)
=
0 and y
′
(
0
)
=
1, with x
∈
[
0,1
]:(16)
y
′′
+
y
′
x
-
1
x
cos
⁡
x
=
0
.Since for the present case (17)
lim
x
→
0
y
′
x
-
1
x
cos
⁡
x
=
0
,it is not difficult to prove that (16) is equivalent to (18)
x
y
′′
+
y
′
-
cos
⁡
x
=
0
.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): (19)
y
(
x
)
=
∫
0
x
sin
t
t
d
t
,
(20)
y
s
(
x
)
=
x
cos
⁡
x
5
cos
⁡
(
cos
⁡
9
)
.Since the exact solution is known we can quantify the error in the interval x
∈
[
0,1
] as (21), which for this case resulted as E
2
=
2.83
×
1
0
-
7: (21)
E
2
=
∫
a
b
y
(
x
)
-
y
s
(
x
)
2
d
x
.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.Figure 2
Exact and approximate solution fory
′′
+
y
′
/
x
-
(
1
/
x
)
cos
⁡
x
=
0 with y
(
0
)
=
0 and y
′
(
0
)
=
1.Case 2.
Let us consider the ODE(22)
x
2
+
1
y
′′
-
2
x
y
-
x
2
-
1
=
0with x
∈
[
0,1
] and initial conditions y
(
0
)
=
0 and y
′
(
0
)
=
1. For this case a closed-form solution is not known. The approximate solution we found is (23)
y
s
(
x
)
=
e
x
-
e
sin
4
[
x
/
(
x
+
sin
8
)
]
,giving a value of 1.76
×
1
0
-
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.Figure 3
Approximate solution forx
2
+
1
y
′′
-
2
x
y
-
x
2
-
1
=
0 with y
(
0
)
=
0 and y
′
(
0
)
=
1.Case 3.
Let us consider the ODE studied in [21]: (24)
y
′
+
y
2
=
2
+
x
2
+
2
x
tanh
(
x
)with x
∈
[
0,1
] and initial condition y
(
0
)
=
0. The proposed method is able to find the exact solution y
(
x
)
=
x
+
tanh
x, 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]: (25)
y
′′
-
y
y
′′
-
0.5
y
′
2
-
0.5
=
0with x
∈
[
0,1
] and boundary conditions y
(
0
)
=
0 and y
(
1
)
=
-
0.5. The approximate solution we obtained with a value of 3
×
1
0
-
3 for (4) using 100 individuals of dimension 30 and 10 000 generations was (26)
y
s
x
=
0.125
x
e
-
sinh
⁡
⁡
erf
⁡
⁡
cosh
⁡
⁡
x
-
erf
⁡
0.999978
-
0.47942
6
x
.In Figure 4 we show the approximate solution we found.Figure 4
Approximate solution fory
′′
-
y
y
′′
-
0.5
y
′
2
-
0.5
=
0 with y
(
0
)
=
0 and y
(
1
)
=
-
0.5.Case 5.
Let us consider the temperature distribution equation on a uniformly thick rectangular fin radiation to free space studied in [23, 24]: (27)
y
′′
-
2
y
4
=
0with x
∈
[
0,1
] and boundary conditions y
(
1
)
=
1 and y
′
(
0
)
=
0.For this case, the approximate solution that the proposed methodology has obtained is(28)
y
s
(
x
)
=
sin
(
x
)
tan
-
1
(
0.25
tan
(
x
)
)
+
0.693147
,when we use 300 individuals of dimension 30 and 10 000 generations. A final value of 9
×
1
0
-
4 was obtained for the objective functional (4). In Figure 5 we show the approximate solution we found.Figure 5
Approximate solution fory
′′
-
2
y
4
=
0 with y
(
1
)
=
1 and y
′
(
0
)
=
0.Case 6.
Consider the differential equation modeling the cooling process of a lumped system by combined convection and radiation [25]:(29)
y
′
+
y
4
+
y
=
0with x
∈
[
0,1
] and initial condition y
(
0
)
=
1.Applying the proposed methodology we obtained the approximate solution(30)
y
s
(
x
)
=
1
erf
⁡
(
x
)
e
x
+
e
x
.For this case the value of (4) was 2
×
1
0
-
4, when using 500 individuals of dimension 30 and 5 000 generations.In Figure6 we show the approximate solution we found.Figure 6
Approximate solution fory
′
+
y
4
+
y
=
0 with y
(
0
)
=
1.Case 7.
Consider the Duffing equation studied in [26]: (31)
y
′′
(
x
)
+
0.4
y
′
(
x
)
+
1.1
y
(
x
)
+
y
3
(
x
)
=
2.1
cos
⁡
(
1.8
x
)
,but with boundary conditions y
(
0.5
)
=
0, y
(
1
)
=
1. The approximate solution we found is (32)
y
s
x
=
arctan
ln
⁡
x
tanh
cos
⁡
5
tan
9
+
3
+
x
+
ln
⁡
sin
arctan
3
x
+
4
+
x
+
e
-
4
,giving a value of 5.2
×
1
0
-
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.Figure 7
Approximate solution fory
′′
(
x
)
+
0.4
y
′
(
x
)
+
1.1
y
(
x
)
+
y
3
(
x
)
=
2.1
cos
⁡
(
1.8
x
), with y
(
0.5
)
=
0 and y
(
1
)
=
1.Case 8.
We close this section considering three examples of the Van der Pol equation:(33)
ε
y
′′
+
δ
+
β
y
2
y
′
-
μ
y
=
F
sin
(
ω
t
)studied in [27].Let us considerε
=
1, δ
=
-
0.05, β
=
0.05, μ
=
-
1, and F
=
0, so we have (34)
y
′′
-
0.05
1
-
y
2
y
′
+
y
=
0
.By applying the proposed methodology to this equation and considering boundary conditions y
(
0
)
=
0 and y
′
(
0
)
=
0.5, we obtain (35)
y
s
(
x
)
=
0.00182376
+
0.5
sin
x
.For this case the final value of the objective functional (4) is 3
×
1
0
-
4, when using 300 individuals of dimension 30 and 5 000 generations. In Figure 8 we show the approximate solution we have found.Figure 8
Approximate solution fory
′′
-
0.05
1
-
y
2
y
′
+
y
=
0 with y
(
0
)
=
0 and y
′
(
0
)
=
0.5.Case 9.
If we now considerε
=
0.1, δ
=
1, β
=
-
1, μ
=
-
1, and F
=
0 for (33), we obtain (36)
0.1
y
′′
+
1
-
y
2
y
′
+
y
=
0
.Considering boundary conditions y
(
0
)
=
0.5 and y
′
(
0
)
=
1 and applying the proposed methodology, we obtain the approximate solution(37)
y
s
x
=
tanh
arctan
⁡
x
+
arctan
cos
⁡
cosh
⁡
cosh
⁡
0.205087
(
x
)
fggggggggfgfg
cosh
⁡
cosh
⁡
0.205087
(
x
)
·
arctan
3.49384
x
.For this case, the final value for the objective functional (4) is 7
×
1
0
-
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.Figure 9
Approximate solution for0.1
y
′′
+
1
-
y
2
y
′
+
y
=
0 with y
(
0
)
=
0.5 and y
′
(
0
)
=
1.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 (38)
1000
y
′′
+
1
-
y
2
y
′
+
y
=
0
.By applying the proposed methodology to solve this equation, we obtain the approximate solution (39)
y
s
(
x
)
=
x
+
0.498241
0.99918
.For this case the value of (4) is 2
×
1
0
-
2, when using 300 individuals of dimension 30 and 20 000 generations. Figure 10 shows the found approximate solution.Figure 10
Approximate solution for1000
y
′′
+
1
-
y
2
y
′
+
y
=
0 with y
(
0
)
=
0.5 and y
′
(
0
)
=
1.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.



---

*Source: 910316-2015-03-11.xml*