# Combining Legendre's Polynomials and Genetic Algorithm in the Solution of Nonlinear Initial-Value Problems

Authors: Osvaldo Guimarães; José Roberto C. Piqueira; Marcio Lobo Netto
Journal: Mathematical Problems in Engineering
 (2011)

---

## Abstract

This work develops a method for solving ordinary differential equations, that is, initial-value problems, with solutions approximated by using Legendre's polynomials. An iterative procedure for the adjustment of the polynomial coefficients is developed, based on the genetic algorithm. This procedure is applied to several examples providing comparisons between its results and the best polynomial fitting when numerical solutions by the traditional Runge-Kutta or Adams methods are available. The resulting algorithm provides reliable solutions even if the numerical solutions are not available, that is, when the mass matrix is singular or the equation produces unstable running processes.


---

## Body

## 1. Introduction

Mathematical modeling is a constant in engineering daily life, and expressing phenomena by using differential equations is almost mandatory. Successful modeling implies ability to make simplifying assumptions but, even under these assumptions, the solution of differential equations implies a lot of analytical and numerical efforts [1].Considering these difficulties and sometimes the impossibility of obtaining analytical solutions, several numerical methods were proposed and, nowadays, the great majority of these methods are part of symbolic and numeric softwares as MATLAB and Mathematica, for instance [2]. These methods are accurate and fast, but each one has its limitations, depending on the type of equation to be solved. The main difficulties appear when the mass-matrix is singular or the equation generates instability in the running process [3].Here, trying to deal with this kind of problem for equations with continuous solutions, a method for initial-value problems in ordinary differential equations is developed, assuming that the space of solutions is a Hilbert space, equipped with a Legendre's polynomial basis [4].The central idea is to look for the best polynomial coefficient combination, in order to satisfy the original differential equation in the whole interval, instead of searching a polynomial that better fits a set of points in the given interval.Consequently, a new error criterion must be defined and a program by using differential evolution methods with elitism [5–7] is developed, based on this new criterion.The determination of the search space developed here is performed by using sequential contractions, allowing fast evolutions with noticeable reductions in the error bands when compared with random search in large populations, as proposed by [5].Several examples are solved by the method presented and by using the classical numerical methods, with the accuracy of the results being compared.In Section2, the error criterion is discussed and a combination of a global and local criteria is proposed, providing optimal fitting of the polynomial, regarding the function and its derivatives. Section 3 presents the analytical basis of the method, by using Legendre's polynomials to optimize the solutions of an initial-value problem.The use of the differential genetic algorithm with elitism, guaranteeing stability, applied to the solution of ordinary differential equations, emphasizing the choice of the coefficients under the optimizing criteria developed in Section2, is described in Section 4. Finally, in Section 5, some examples are presented and the numerical quality of the solutions is discussed.


## 2. Error Criterion

Usually, the adopted optimization criterion to approximate a function by a polynomial in a given interval is to minimize the square error integral [8]. This kind of approach minimizes the global error along the whole interval, but produces a large value for the local error at the starting and ending points of the interval.Consequently, when this criterion is applied to an initial-value problem, the approximation does not produce the best result, as the error for the starting point of the interval remains large. To solve this problem, an error criterion considering the local and global results is proposed, minimizing the product of two factors: the maximum error at the starting point of the interval and the global integral error.To give an idea of this procedure, the differential equation given by (2.1) with initial condition y(-1)=e-1, that admits the solution given by (2.2) in the [-1,1] interval, is studied.(2.1)dydx+2xe-x2=0,(2.2)y=e-x2.If one tries to obtain, by using MATLAB, the 8th-order Legendre's polynomial that better fits to (2.2), the result, with double-precision format, is given by (2.3). Considering the substitution of y, given by (2.3), in the original differential equation (2.1), the integral error is 1.9.10-8 for the whole interval(2.3)y=.0263x8-.1551x6+.4963x4-.9996x2+1,0000.If one continues the process, with the result given by (2.3) as an initial solution candidate to the original equation (2.1), by using the method developed in the next Sections, a solution given by (2.4), with an integral error 8.7.10-9 for the whole interval, is obtained(2.4)y=.0242x8-.1511x6+0.4940x4-.9992x2+1,0000.This result is related to the fact that the polynomial given by (2.3) is adequate to fit the function given by (2.2), without considering its derivatives. In order to adequate the solution to the differential equation, given by (2.1), the derivatives must be adjusted in the same way.


## 3. Analytic Foundations

In this section, it is considered the Hilbert space of piecewise continuous functions, defined in a closed interval[a,b], with a scalar product between f and g, belonging to the space, defined as [4]
(3.1)〈f,g〉=∫abfgdx.Considering the real interval[-1,1], the Legendre's polynomials set is an orthogonal basis for the space of piecewise continuous function S, defined in this interval, that is, for any function f∈S(3.2)lim⁡n→∞∫-11[f(x)-∑k=0nakPk(x)]2dx=0,
with Pk being the k-degree normalized Legendre's polynomial and ak=〈f,Pk〉 [4, 9].It is important to notice that either the integral or the derivative of a Legendre's polynomial can be expressed as a linear combination of the Legendre's polynomial basis [10], as the following theorems, with proofs presented in Appendix A, show.Theorem 3.1.
CallingPn(x) the  n-degree Legendre's polynomial,
(3.3)∫-1xPn(ξ)dξ=Pn+1(x)-Pn(x)2n+1,n∈N*.Theorem 3.2.
CallingPn(x) the  n-degree Legendre's polynomial,
(3.4)Pn′=∑k=0n-1(2k+1)2[1-(-1)k+n]Pk.Considering the theorems above, about integration and differentiation of Legendre's polynomials, it is interesting to express these results in a matrix form.


### 3.1. Operational Matrices

As shown in [11], Theorem 3.3 allows the construction of matrices expressing any linear operation executed over the Hilbert space of the continuous functions. Here, numerical methods for solving differential equations by using Legendre's polynomials are developed. Therefore, constructing these matrices for integration and differentiation operations expressed in a Legendre's polynomials basis is useful.Letf(u) a linear combination of a generic basis, {Bk}, in the Hilbert space of continuous functions(3.5)f(u)=∑k=0nckBk(u).Theorem 3.3.
LetZ a  (n+1)×(n+1) matrix describing the resulting coefficients from a linear operation,  α, over the generic basis,  {Bk}, as functions of the same basis. Calling  V=[v0v1…vn] the resulting coefficients vector from the linear operation,  α, over  f(u):  VT=ZTCT, with  C being the vector of the coefficients described in (3.5).The procedure described by Theorem3.3, proved in Appendix A, can be used for any basis in a Hilbert space, with the calculation of matrix Z representing the main operation to be done.Combining Theorems3.1 and 3.3, taking Legendre's polynomials as the basis in the Hilbert space of continuous functions, if ZLi is (n+2)×(n+2) and considering α the integral operation(3.6)ZLi(1,1)=1;ZLi(1,2)=1,ZLi(k,k-1)=12k-1;k=2,…n+1,ZLi(k,k+1)=-12k-1;k=2,…n+1,ZLi(n+2,n+1)=-12n+3.All the elements ofZLi, not defined by (3.6), are equal to zero. For instance, if n=2, ZLi can be written as(3.7)ZLi=[1100-1301300-1501500-170].It can be observed that the elementZLi(n+2,n+3)=ZLi(4,5) does not appear, because the third degree Legendre's polynomial was not integrated.Combining Theorems3.2 and 3.3, taking Legendre's polynomials as the basis in the Hilbert space of continuous functions, if Zdi is (n+1)×(n+1) and considering α the derivative operation(3.8)ZLd(k+1,j+1)=[1-(-1)k+j]2j+12;k=1,…n-1;j=0,…k.All the elements ofZLd, not defined by (3.8), are equal to zero. For instance, for n=3, ZLd can be written as(3.9)ZLd=[0000100003001050].The integration matrixZLi is nonsingular, and, therefore, its inverse always exists. On the other hand, ZLd needs a correction if c0 in (3.5) is such that c0≠-∑n=1k+1cn(-1)n. The correction to be done is to modify c0, imposing the condition described, without changing the derivative of f(u).In [10], the operational matrix was derived for the particular case where α is an integration operation. Theorem 3.3, developed here, permits general matrix expressions for all linear α-operations and considering any kind of basis in the Hilbert space of continuous functions. [12] shows, additionally, this development for Fourier trigonometric and canonical basis.


### 3.2. Conservation of Initial Conditions in a Legendre's Basis

When using evolution methods to treat initial-value problems, expressing solutions in a Legendre's basis, the initial conditions must be conserved during all algorithm running process. This conservation depends on the correct combination of the coefficients,ci, to express the function in any step of the running process. This combination generates an algebraic linear system of equations that admits infinite solutions. It is solved choosing a subset of q coefficients, with q representing the order of the original differential equation, and expressing then as a function of the others.In order to obtain the coefficients of a function derivative, the differentiation matrix, built as explained in the former subsection, has to be used.As an example, if a third-order equation with solutions in the[-1,1] interval is studied with seven coefficients for the expression in Legendre's basis, one of the possible relationship among the coefficients is
(3.10)c0=5c3-21c4+56c5-120c6+23y0(2)+y0(1)+y0,c1=9c3-35c4+90c5-189c6+y0(2)+y0(1),c2=5c3-15c4+35c5-70c6+13y0(2).
where (y0,y0(1),y0(2)) is the vector of the initial conditions. If one needs to obtain the coefficients of the integral of the function, the integration matrix, built as explained in the former subsection, must be used, but the initial condition must be considered, as explained below.Considering thatF(u)=∫-1uf(ξ)dξ and its Legendres expansion is given by: F(u)=∑k=0n+1bkPk(u), with Pk being the k-degree Legendre's polynomial and that the initial condition for the original differential equation is F(-1)=β, b0 obtained by the integration matrix must be replaced by b¯0=b0+β. For q differential equations, this process must be repeated q times.


### 3.3. Domain Transposition

The purpose of this work is to develop numerical methods to solve differential equations with solutions defined in a generic closed real interval [a,b], but the exposed analytical foundations and other developed methods consider only particular intervals [13, 14]. Here, a generic transformation providing tools to transpose any general closed domain to a particular one is studied for up to fourth-order initial-value problems. Besides, the impact of this transformation in original differential equation and initial conditions is analyzed.Letu=T(x);T:[a,b]→[c,d] be a bijective transformation and its inverse T-1, and y=f(x);f:[a,b]→R a solution of an initial value problem, with a vector Y0 representing the initial condition.The interval transposed functionyt(u):[c,d]→R is defined by yt(u)=y(x), with x=T-1(u) and, for the first derivative, calling D0=y(x)(3.11)D1=y(1)(x)=yt(1)(u)⋅dudx.Callingλ=du/dx=1/(dx/du) and D1t=yt(1),(3.12)D1=D1tλ.Consequently, the second and third derivatives are(3.13)D2=D2t⋅λ2+D1t⋅dλdu⋅λ,D3=D3t⋅λ3+3D2tλ2dλdu+D1tλ[(dλdu)2+λd2λd2u].For the sake of clarity, the forth derivative is written by grouping the terms expressed in function ofD1t, D2t, D3t, and D4t as follows:(i)
terms inD4t: D4tλ4,(ii)
terms inD3t: 6D3tλ3(dλ/du),(iii)
terms inD2t: D2tλ2(7(dλ/du)2+4λ(d2λ/du2)),(iv)
terms inD1t: D1tλ[(dλ/du)3+4λ(dλ/du)(d2λ/d2u)+λ2(d3λ/du3)].Combining the expressions above,(3.14)D4=D4tλ4+6D3tλ3dλdu+D2tλ2(7(dλdu)2+4λd2λdu2)+D1tλ[(dλdu)3+4λdλdud2λd2u+λ2d3λdu3].As the solution of a differential equation is tied to the initial conditions, it is necessary to establish domain transformations as described above, in an inverse way. Therefore, callingρ=dx/du=1/(du/dx) and considering equations (3.12), (3.13), and (3.14) in an inverse way, with the replacement of λ by ρ, the initial conditions, y(0), y(1)(0), y(2)(0), and y(3)(0) are given by(3.15)y0t=y0,y1t=y1ρ,y2t=y2⋅ρ2+y1⋅dρdx⋅ρ,y3t=y3⋅ρ3+3y2ρ2dρdx+y1ρ[(dρdx)2+ρd2ρd2x].IfT is a linear transformation, expressions (3.12), (3.13), (3.14), and (3.15), are considerably simplified, as λ and ρ are constants.In spite of linear transformation being simple, it does not allow to work with equations defined in an infinite domain. In order to take care of this problem, the transformationT, given by(3.16)u=T(x)=(d-c)(arctan(x)-v1)v2-v1+c,x=T-1(u)=tan((v2-v1)(u-c)d-c+v1),
where v1=arctana and v2=arctanb can be used.


## 3.1. Operational Matrices

As shown in [11], Theorem 3.3 allows the construction of matrices expressing any linear operation executed over the Hilbert space of the continuous functions. Here, numerical methods for solving differential equations by using Legendre's polynomials are developed. Therefore, constructing these matrices for integration and differentiation operations expressed in a Legendre's polynomials basis is useful.Letf(u) a linear combination of a generic basis, {Bk}, in the Hilbert space of continuous functions(3.5)f(u)=∑k=0nckBk(u).Theorem 3.3.
LetZ a  (n+1)×(n+1) matrix describing the resulting coefficients from a linear operation,  α, over the generic basis,  {Bk}, as functions of the same basis. Calling  V=[v0v1…vn] the resulting coefficients vector from the linear operation,  α, over  f(u):  VT=ZTCT, with  C being the vector of the coefficients described in (3.5).The procedure described by Theorem3.3, proved in Appendix A, can be used for any basis in a Hilbert space, with the calculation of matrix Z representing the main operation to be done.Combining Theorems3.1 and 3.3, taking Legendre's polynomials as the basis in the Hilbert space of continuous functions, if ZLi is (n+2)×(n+2) and considering α the integral operation(3.6)ZLi(1,1)=1;ZLi(1,2)=1,ZLi(k,k-1)=12k-1;k=2,…n+1,ZLi(k,k+1)=-12k-1;k=2,…n+1,ZLi(n+2,n+1)=-12n+3.All the elements ofZLi, not defined by (3.6), are equal to zero. For instance, if n=2, ZLi can be written as(3.7)ZLi=[1100-1301300-1501500-170].It can be observed that the elementZLi(n+2,n+3)=ZLi(4,5) does not appear, because the third degree Legendre's polynomial was not integrated.Combining Theorems3.2 and 3.3, taking Legendre's polynomials as the basis in the Hilbert space of continuous functions, if Zdi is (n+1)×(n+1) and considering α the derivative operation(3.8)ZLd(k+1,j+1)=[1-(-1)k+j]2j+12;k=1,…n-1;j=0,…k.All the elements ofZLd, not defined by (3.8), are equal to zero. For instance, for n=3, ZLd can be written as(3.9)ZLd=[0000100003001050].The integration matrixZLi is nonsingular, and, therefore, its inverse always exists. On the other hand, ZLd needs a correction if c0 in (3.5) is such that c0≠-∑n=1k+1cn(-1)n. The correction to be done is to modify c0, imposing the condition described, without changing the derivative of f(u).In [10], the operational matrix was derived for the particular case where α is an integration operation. Theorem 3.3, developed here, permits general matrix expressions for all linear α-operations and considering any kind of basis in the Hilbert space of continuous functions. [12] shows, additionally, this development for Fourier trigonometric and canonical basis.


## 3.2. Conservation of Initial Conditions in a Legendre's Basis

When using evolution methods to treat initial-value problems, expressing solutions in a Legendre's basis, the initial conditions must be conserved during all algorithm running process. This conservation depends on the correct combination of the coefficients,ci, to express the function in any step of the running process. This combination generates an algebraic linear system of equations that admits infinite solutions. It is solved choosing a subset of q coefficients, with q representing the order of the original differential equation, and expressing then as a function of the others.In order to obtain the coefficients of a function derivative, the differentiation matrix, built as explained in the former subsection, has to be used.As an example, if a third-order equation with solutions in the[-1,1] interval is studied with seven coefficients for the expression in Legendre's basis, one of the possible relationship among the coefficients is
(3.10)c0=5c3-21c4+56c5-120c6+23y0(2)+y0(1)+y0,c1=9c3-35c4+90c5-189c6+y0(2)+y0(1),c2=5c3-15c4+35c5-70c6+13y0(2).
where (y0,y0(1),y0(2)) is the vector of the initial conditions. If one needs to obtain the coefficients of the integral of the function, the integration matrix, built as explained in the former subsection, must be used, but the initial condition must be considered, as explained below.Considering thatF(u)=∫-1uf(ξ)dξ and its Legendres expansion is given by: F(u)=∑k=0n+1bkPk(u), with Pk being the k-degree Legendre's polynomial and that the initial condition for the original differential equation is F(-1)=β, b0 obtained by the integration matrix must be replaced by b¯0=b0+β. For q differential equations, this process must be repeated q times.


## 3.3. Domain Transposition

The purpose of this work is to develop numerical methods to solve differential equations with solutions defined in a generic closed real interval [a,b], but the exposed analytical foundations and other developed methods consider only particular intervals [13, 14]. Here, a generic transformation providing tools to transpose any general closed domain to a particular one is studied for up to fourth-order initial-value problems. Besides, the impact of this transformation in original differential equation and initial conditions is analyzed.Letu=T(x);T:[a,b]→[c,d] be a bijective transformation and its inverse T-1, and y=f(x);f:[a,b]→R a solution of an initial value problem, with a vector Y0 representing the initial condition.The interval transposed functionyt(u):[c,d]→R is defined by yt(u)=y(x), with x=T-1(u) and, for the first derivative, calling D0=y(x)(3.11)D1=y(1)(x)=yt(1)(u)⋅dudx.Callingλ=du/dx=1/(dx/du) and D1t=yt(1),(3.12)D1=D1tλ.Consequently, the second and third derivatives are(3.13)D2=D2t⋅λ2+D1t⋅dλdu⋅λ,D3=D3t⋅λ3+3D2tλ2dλdu+D1tλ[(dλdu)2+λd2λd2u].For the sake of clarity, the forth derivative is written by grouping the terms expressed in function ofD1t, D2t, D3t, and D4t as follows:(i)
terms inD4t: D4tλ4,(ii)
terms inD3t: 6D3tλ3(dλ/du),(iii)
terms inD2t: D2tλ2(7(dλ/du)2+4λ(d2λ/du2)),(iv)
terms inD1t: D1tλ[(dλ/du)3+4λ(dλ/du)(d2λ/d2u)+λ2(d3λ/du3)].Combining the expressions above,(3.14)D4=D4tλ4+6D3tλ3dλdu+D2tλ2(7(dλdu)2+4λd2λdu2)+D1tλ[(dλdu)3+4λdλdud2λd2u+λ2d3λdu3].As the solution of a differential equation is tied to the initial conditions, it is necessary to establish domain transformations as described above, in an inverse way. Therefore, callingρ=dx/du=1/(du/dx) and considering equations (3.12), (3.13), and (3.14) in an inverse way, with the replacement of λ by ρ, the initial conditions, y(0), y(1)(0), y(2)(0), and y(3)(0) are given by(3.15)y0t=y0,y1t=y1ρ,y2t=y2⋅ρ2+y1⋅dρdx⋅ρ,y3t=y3⋅ρ3+3y2ρ2dρdx+y1ρ[(dρdx)2+ρd2ρd2x].IfT is a linear transformation, expressions (3.12), (3.13), (3.14), and (3.15), are considerably simplified, as λ and ρ are constants.In spite of linear transformation being simple, it does not allow to work with equations defined in an infinite domain. In order to take care of this problem, the transformationT, given by(3.16)u=T(x)=(d-c)(arctan(x)-v1)v2-v1+c,x=T-1(u)=tan((v2-v1)(u-c)d-c+v1),
where v1=arctana and v2=arctanb can be used.


## 4. Evolution-Adaptive Process for Solving Initial-Value Problems

In this section, the main question to be answered is if it is possible, with the computation power available, to generate random coefficients for ann-degree polynomial and, under a sequential algorithm, like differential evolution algorithm, over the coefficients, to obtain the exact solution of an initial value problem.Due to the randomness of the original genetic algorithms, the solution space has its regions visited in a random way, and there is no convergence. Nevertheless, when differential versions of genetic algorithms with elitism are applied, some stable processes and satisfactory results are obtained for particular cases [15–18], including complex biological problems [19, 20] and quantum mechanics [21, 22].Here, a differential genetic algorithm with an implicit learning process, given by a good choice of the region of the space of solutions based on a restriction scheme that continuously contracts the volume of the space of research, guaranteeing stability, is developed.Considering that an-degree polynomial is the best candidate to the solution of an initial value problem, in a particular step of a sequential algorithm, under a certain error criterion, it is natural to ask how far of the exact solution it is. In order to treat this problem, the original differential equation is written with the higher-order derivative term expressed as a function of the lower-order derivatives, the original function, and the independent variable as:(4.1)dqyduq=G(yt,yt(1),…,yt(q-1),u).Replacing the best solution candidate asy, into G given by (4.1), it is possible to evaluate G in the defined interval and fitting the obtained results by using a Legendre's polynomial up to (n-q)-degree. The fitting process follows Gauss-Legendre quadrature method [23, 24]. After that, the obtained polynomial is q-times integrated, by using the integration matrix given by [ZLi]T, defined in Section 3.1, and considering the given initial conditions.Comparing this result with the original candidate, the dispersion is obtained. This dispersion defines the band variation for the coefficients of the polynomial candidate, in the next step of the algorithm. For the first step, the candidate is the null polynomial, providing the first band variation for each coefficient of the Legendre's polynomial.To perform the differential evolution algorithm, a population is randomly mounted inside the first band variation. Then, by using a combination of crossing-over, random and mutation, new off-springs are generated feeding the sequence of the algorithm.The program starts a tournament selection. If the best output of the tournament overcomes the former winner, a new band variation is constructed. Otherwise, by using elitism with random cross-over and mutation [25], the algorithm creates a new generation with each gene (coefficient) varying inside the new band, avoiding stagnation in local minima.The described procedure is possible, even when the mass-matrix is singular [3] because the Gauss-Legendre's abscissae do not include the interval extremes [23, 24]. If the numerical integration solution is available, it can be the starting point of the algorithm, instead of the null-polynomial, shortening the running time of the procedure.


## 5. Results and Analysis

Based on the analytical development presented in Section3, by using a differential evolution algorithm, a method for solving initial value problems, from now on abbreviated by SIV and detailed in Appendix B, was developed, considering that the solutions are expressed as linear combinations of Legendre's polynomials, in the transposed domain.Starting with an interval transposition, theSIV algorithm runs, adjusting the coefficients of the polynomial solution, minimizing the residual error when replacing the polynomial approximation and derivatives and transposed initial conditions in the transposed differential equation. After a chosen number of generations, applying the inverse transposition, the solution is obtained, in the original domain.Here, the application of theSIV algorithm is shown considering four examples of nonlinear up to fourth-order equations. The obtained solutions are compared with the analytic ones, when they exist.The algorithm is implemented with 66 chromosomes and5% of mutation. The CPU is a PC Dual-core with 2 GB RAM and running MATLAB-7.In the cases where it is not possible to express an analytic solution by a finite combination of elementary functions, the solution obtained bySIV is compared with the obtained by Adams' method, when the maximum MATLAB precision is used [2, 3].


### 5.1. Linear Transposition: Zero-Order Bessel Equation

The first chosen example is the zero-order Bessel equation, corresponding to(5.1)xD2+D1+D0x=0,
with Di denoting the ith derivative of the searched function, related to the independent variable x, x∈[0,1], and with the initial condition given by the vector [1,0].Equation (5.1) has analytic solution given by(5.2)J0(x)=∑k=0∞(-1)k(k!)2(x2)2k.In order to transpose the domain from the[0,1] to [-1,1], a linear transformation is used. For this kind of transformation, considering that a function is described in an interval [a,b] and will be transposed to an interval [c,d], the relation (u-c)/(d-c)=(x-a)/(b-a), with u=T(x)=(d-c)/(b-a)·x+(bc-ad)/(b-a), is used.TheSIV algorithm was applied to (5.1), with a 26-coefficient polynomial approximation, starting with the null-polynomial.The final result was compared with the analytic solution (5.2). As, in this case, the mass matrix is singular [3, 23, 24], the Runge-Kutta and Adams algorithms cannot be initialized and, consequently, they are not used for comparisons.After two hundred generations, an evolved Legendre's series approximation was obtained with very low error margins, as shown in Figure1 for the obtained solution of (5.1) and the corresponding derivatives.Zero-order Bessel equation: error in the function and derivatives after 200 generations, comparingSIV and analytic solution.
(a)(b)(c)To give an idea about the precision of the method, Figure2 shows the residue by substituting the approximation of the function given by (5.2) by using Legendre's polynomial fitting process [2], and the residue for the same degree evolved solution. Observing Figure 2, it can be noticed that the SIV algorithm improves the accuracy in, at least, two order of magnitude. Besides, the problem of lack of accuracy at the interval extremes is considerably reduced.Figure 2
Zero-order Bessel equation: residue due to the best fitting polynomial compared with the residue bySIV polynomial.


### 5.2. Linear Interval Transposition: Fourth-Order Differential Equation

Now, it will be considered an example of a fourth-order differential equation. DenotingDi the i-derivative of the searched function, the equation to be solved is(5.3)D4+D3⋅D1-D02+3⋅sin2x+(4D0+3)cos⁡x+(D1-sinx)2+D22=0,x∈[0,π], and with the initial condition given by the vector: [0,0,2,0].Equation (5.3) has analytic solution given by(5.4)f(x)=x⋅sin⁡(x).
In order to transpose the domain from the [0,π] to [-1,1], a linear transformation, as shown in Section 5.1, is used.Starting with the solution obtained by using MATLAB (ode113) with maximum precision [2], the SIV method with 22 coefficients gives a solution with the error figures shown in Figure 3, after two hundred generations.A fourth-order nonlinear equation with linear domain transposition: error in the function and derivatives after 200 generations, comparingSIV and analytic solution.
(a)(b)(c)(d)(e)To give an idea about the precision of the method, Figure4 shows the residue by substituting the approximation of the function given by (5.4) by using Legendre's polynomial fitting process [2], and the residue for the same degree evolved solution. Observing Figure 4, it can be noticed that the SIV algorithm improves the accuracy in, at least, two order of magnitude. Besides, the problem of lack of accuracy at the interval extremes is considerably reduced.Figure 4
A fourth-order nonlinear equation with linear domain transposition: residue due to the best fitting polynomial compared with the residue due toSIV polynomial.


### 5.3. Trigonometric Interval Transposition

Another example where an interval transposition is necessary is described here. DenotingDi the i-derivative of the searched function, the equation to be solved is(5.5)D4+4⋅D1⋅[3D1⋅x⋅(x4-1)-9⋅x2+2](1+x2)3⋅tan-1(x)=0,x∈[0,π], and with the initial condition given by the vector: [0,0,2,0].Equation (5.3) has analytic solution given by(5.6)f(x)=[tan-1(x)]2.In order to transpose the domain from[0,π] to [-1,1], a trigonometric transformation is used. For this kind of transformation, considering that a function of x is defined in an interval [a,b], firstly it is transposed with a new function u=tan(v), defined for an interval [v1,v2] of v, such that v1=tan-1(a) and v2=tan-1(b).Then, this interval is transposed to[c,d] for the variable u, by using the transformation:u=(d-c)·(tan-1(x)-v1)/(v2-v1)+c.Starting with the null-polynomial, with seven coefficients and running theSIV process for six hundred generations, taking about 1 minute, a solution with the error figures shown in Figure 5 was obtained.A fourth-order nonlinear equation with trigonometric domain transposition: error in the function and derivatives after 600 generations, comparingSIV and analytic solution.
(a)(b)(c)(d)(e)As a 6th-order approximation for the function given by (5.6) has a very small error, Figure 6 shows the residue of the evolved series, because it is practically the exact solution.Figure 6
A fourth-order nonlinear equation with trigonometric domain transposition: residue due toSIV polynomial.


### 5.4. Equation without Analytic Solution

There are equations such that the analytic solution cannot be expressed by a finite combination of elementary functions. Here, an example of this kind of case is presented. DenotingDi the i-derivative of the searched function, the equation to be solved is(5.7)D3+D2⋅D0=0,x∈[0,5], and with the initial condition given by the vector [0,0,.4699].After a linear interval transposition, theSIV algorithm started with a 21-degree null-polynomial. After four hundred generations, a solution with residual errors shown in Figure 7 was obtained.A third-order nonlinear equation without analytic solution: error in the function and derivatives after 600 generations, comparingSIV and ode 113.
(a)(b)(c)As the analytic solution cannot be expressed as a finite combination of elementary functions, in order to establish comparisons, the Adams' method, implemented with the maximum MATLAB precision [2], was used to perform the Legendre's polynomial fitting for the function. Figure 8 shows that the SIV algorithm improves the precision in, at least, one order of magnitude, attenuating the lack of accuracy at the domain extremes.Figure 8
A third-order nonlinear equation without analytic solution: comparing residues in ode 113 andSIV.


## 5.1. Linear Transposition: Zero-Order Bessel Equation

The first chosen example is the zero-order Bessel equation, corresponding to(5.1)xD2+D1+D0x=0,
with Di denoting the ith derivative of the searched function, related to the independent variable x, x∈[0,1], and with the initial condition given by the vector [1,0].Equation (5.1) has analytic solution given by(5.2)J0(x)=∑k=0∞(-1)k(k!)2(x2)2k.In order to transpose the domain from the[0,1] to [-1,1], a linear transformation is used. For this kind of transformation, considering that a function is described in an interval [a,b] and will be transposed to an interval [c,d], the relation (u-c)/(d-c)=(x-a)/(b-a), with u=T(x)=(d-c)/(b-a)·x+(bc-ad)/(b-a), is used.TheSIV algorithm was applied to (5.1), with a 26-coefficient polynomial approximation, starting with the null-polynomial.The final result was compared with the analytic solution (5.2). As, in this case, the mass matrix is singular [3, 23, 24], the Runge-Kutta and Adams algorithms cannot be initialized and, consequently, they are not used for comparisons.After two hundred generations, an evolved Legendre's series approximation was obtained with very low error margins, as shown in Figure1 for the obtained solution of (5.1) and the corresponding derivatives.Zero-order Bessel equation: error in the function and derivatives after 200 generations, comparingSIV and analytic solution.
(a)(b)(c)To give an idea about the precision of the method, Figure2 shows the residue by substituting the approximation of the function given by (5.2) by using Legendre's polynomial fitting process [2], and the residue for the same degree evolved solution. Observing Figure 2, it can be noticed that the SIV algorithm improves the accuracy in, at least, two order of magnitude. Besides, the problem of lack of accuracy at the interval extremes is considerably reduced.Figure 2
Zero-order Bessel equation: residue due to the best fitting polynomial compared with the residue bySIV polynomial.


## 5.2. Linear Interval Transposition: Fourth-Order Differential Equation

Now, it will be considered an example of a fourth-order differential equation. DenotingDi the i-derivative of the searched function, the equation to be solved is(5.3)D4+D3⋅D1-D02+3⋅sin2x+(4D0+3)cos⁡x+(D1-sinx)2+D22=0,x∈[0,π], and with the initial condition given by the vector: [0,0,2,0].Equation (5.3) has analytic solution given by(5.4)f(x)=x⋅sin⁡(x).
In order to transpose the domain from the [0,π] to [-1,1], a linear transformation, as shown in Section 5.1, is used.Starting with the solution obtained by using MATLAB (ode113) with maximum precision [2], the SIV method with 22 coefficients gives a solution with the error figures shown in Figure 3, after two hundred generations.A fourth-order nonlinear equation with linear domain transposition: error in the function and derivatives after 200 generations, comparingSIV and analytic solution.
(a)(b)(c)(d)(e)To give an idea about the precision of the method, Figure4 shows the residue by substituting the approximation of the function given by (5.4) by using Legendre's polynomial fitting process [2], and the residue for the same degree evolved solution. Observing Figure 4, it can be noticed that the SIV algorithm improves the accuracy in, at least, two order of magnitude. Besides, the problem of lack of accuracy at the interval extremes is considerably reduced.Figure 4
A fourth-order nonlinear equation with linear domain transposition: residue due to the best fitting polynomial compared with the residue due toSIV polynomial.


## 5.3. Trigonometric Interval Transposition

Another example where an interval transposition is necessary is described here. DenotingDi the i-derivative of the searched function, the equation to be solved is(5.5)D4+4⋅D1⋅[3D1⋅x⋅(x4-1)-9⋅x2+2](1+x2)3⋅tan-1(x)=0,x∈[0,π], and with the initial condition given by the vector: [0,0,2,0].Equation (5.3) has analytic solution given by(5.6)f(x)=[tan-1(x)]2.In order to transpose the domain from[0,π] to [-1,1], a trigonometric transformation is used. For this kind of transformation, considering that a function of x is defined in an interval [a,b], firstly it is transposed with a new function u=tan(v), defined for an interval [v1,v2] of v, such that v1=tan-1(a) and v2=tan-1(b).Then, this interval is transposed to[c,d] for the variable u, by using the transformation:u=(d-c)·(tan-1(x)-v1)/(v2-v1)+c.Starting with the null-polynomial, with seven coefficients and running theSIV process for six hundred generations, taking about 1 minute, a solution with the error figures shown in Figure 5 was obtained.A fourth-order nonlinear equation with trigonometric domain transposition: error in the function and derivatives after 600 generations, comparingSIV and analytic solution.
(a)(b)(c)(d)(e)As a 6th-order approximation for the function given by (5.6) has a very small error, Figure 6 shows the residue of the evolved series, because it is practically the exact solution.Figure 6
A fourth-order nonlinear equation with trigonometric domain transposition: residue due toSIV polynomial.


## 5.4. Equation without Analytic Solution

There are equations such that the analytic solution cannot be expressed by a finite combination of elementary functions. Here, an example of this kind of case is presented. DenotingDi the i-derivative of the searched function, the equation to be solved is(5.7)D3+D2⋅D0=0,x∈[0,5], and with the initial condition given by the vector [0,0,.4699].After a linear interval transposition, theSIV algorithm started with a 21-degree null-polynomial. After four hundred generations, a solution with residual errors shown in Figure 7 was obtained.A third-order nonlinear equation without analytic solution: error in the function and derivatives after 600 generations, comparingSIV and ode 113.
(a)(b)(c)As the analytic solution cannot be expressed as a finite combination of elementary functions, in order to establish comparisons, the Adams' method, implemented with the maximum MATLAB precision [2], was used to perform the Legendre's polynomial fitting for the function. Figure 8 shows that the SIV algorithm improves the precision in, at least, one order of magnitude, attenuating the lack of accuracy at the domain extremes.Figure 8
A third-order nonlinear equation without analytic solution: comparing residues in ode 113 andSIV.


## 6. Conclusions

A method that combines genetic algorithms with Legendre's polynomial approximation was proposed, in order to solve nonlinear initial-value problems, up to 4th-order differential equations.The method presents several advantages.(i)
It provides analytical expressions for the solutions and their derivatives for the whole domain, instead of a low-degree piecewise polynomial;(ii)
It runs even if the mass matrix is singular, differently from the traditional Runge-Kutta and Adams' methods;(iii)
It is able to start the process of finding solution from a null-polynomial;(iv)
It considerably improves the precision of the solutions, solving the problem related to the lack of accuracy in the interval extremes;(v)
The spectral evolution runs fast once the evolution of each coefficient does not affect the others, due to orthogonality.



---

*Source: 521342-2011-07-26.xml*