Error Upper Bounds for a Computational Method for Nonlinear Boundary and Initial-Value Problems

Differential equations are ubiquitous in engineering daily life but their solutions are sometimes very difficult, mainly if they are nonlinear. Several of them do not have an analytical solution describable by a finite combination of elementary functions, or even by an unlimited series with a determinable recurrence relation. In previous works, analytical and numerical results were obtained for nonlinear differential equations 1–3 , and an algorithm called SIV Solving Initial Value was developed. Here, the problem of determining the error limits for SIV is emphasized, providing quality parameters for the method. In Section 2, some general theoretical considerations about a numerical method for differential equations, by using expansions in Hilbert space, are presented. Error upper bounds estimation is discussed in Section 3, including some examples. Section 4 concludes the reasoning.


Introduction
Differential equations are ubiquitous in engineering daily life but their solutions are sometimes very difficult, mainly if they are nonlinear.Several of them do not have an analytical solution describable by a finite combination of elementary functions, or even by an unlimited series with a determinable recurrence relation.
In previous works, analytical and numerical results were obtained for nonlinear differential equations 1-3 , and an algorithm called SIV Solving Initial Value was developed.Here, the problem of determining the error limits for SIV is emphasized, providing quality parameters for the method.
In Section 2, some general theoretical considerations about a numerical method for differential equations, by using expansions in Hilbert space, are presented.Error upper bounds estimation is discussed in Section 3, including some examples.Section 4 concludes the reasoning.

Methodology
Problems described by differential equations can be submitted to initial conditions and, in those cases, they are called initial-value problems IVP and there are a number of softwares dedicated to their solutions, based on Runge-Kutta and Adams methods 4 constantly improved 5, 6 , mainly considering applications to physics problems 7-10 .
Lipschitz theorem gives a sufficient condition for a differential equation to have a unique solution and, under the theorem assumptions, mass matrix is guaranteed to be nonsingular, providing that the classical methods can start running and obtaining, for sequential intervals, continuous expansions for limited functions 11 .The precision is only limited by the computational apparatus.
Besides IVP problems, there are the boundary value problems BVP that, mainly, have been solved by using numerical methods derived from Galerkin and variational methods 12, 13 .Previous works 1-3 developed analytical and numerical methods, based on a spectral approach, that can be applied either to IVP or to BVP.
Details of the method, called SIV, are described mainly in 3 .Here, some hints are presented for the sake of clarity, in order to discuss precision issues.
The implemented evolutive algorithm is genetic and differential 14 , conceived for an evolution on a Hilbert space, with which gene is represented by the coefficients of a possible expansion in terms of a basis of this space and corresponding to a sixty-individual vector population chromosomes .Each generation is submitted to ranking contests.
The error obtained substituting the vector in the equation is the criterium to determine the individual ranking position.The lesser the error is, the better is the individual ranking position.Selecting the reproducers for the next generation is performed by attributing greater choice probability to the best ranked individuals 3 .
Heritage is simulated by taking two vectors and combining them to generate a descendent that is a new vector with the same dimension of the original vectors.The criteria to choose if the gene in a locus is given by one or another original vector are the cross-over with the loci randomly chosen 15 .
The introduction of learning rules in the algorithm 1, 3 reduced about a hundred times the processing time because the search space is continually contracted with error margins decreasing.
Considering the solution f u to be continuous and expressed as a linear combination of the elements from an orthogonal basis B k u in the Hilbert functional space is the main assumption of the method presented in 3 .
Under these conditions, f u 1, the Legendre polynomial basis is obtained with g k 2/ 2k 1 .For the same interval and by using generic w u , Jacobi polynomials are obtained.Legendre and Chebyshev polynomials are particular cases of Jacobi polynomials.For half-limited intervals, the Laguerre polynomials are used; for unlimited intervals, the Hermite polynomials are adequate.
Solving differential equations with orthogonal series approximations have been thoroughly studied in the last three decades, providing efficient algorithms 16-19 .In this approach, it is possible to transform a differential equation into an algebraic one, giving the coefficients of the series expansion.This transformation is performed by using operational matrices, extensively discussed in 3 .
Completeness and orthogonality in Hilbert spaces are related to certain domains, depending on the basis chosen.Consequently, it is important to develop basis transformations to consider this fact.In 1, 3 a method for transposing domains changing a variable that transforms the differential equation was developed.Besides, initial and boundary conditions need to be transposed in the same way and the program described in 1, 3 takes care of this transposition, too.

Error Upper Bounds
One of the main problems in the numerical solutions of differential equations is how to estimate the error margins.In the majority of the cases, it is possible to stipulate the maximum tolerable error for each step but, as the integration for the whole interval comprehends thousands of steps, the cumulative error is difficult to be obtained 20 .
Here, a method to determine the upper bounds for the absolute error in an integration process is described, even if the algorithm and the analytical solution are unknown.These ideas can be used, in order to estimate the confidence for numerical approximations of the solution 21 .
Starting with the differential equation written in the following form: , u , evaluating G for the whole interval.This reasoning allows the calculation of the relative dispersion of the coefficients of the series that approximates the solution.In order to illustrate the process, Example 1 is developed.

Example 1
As an example of dispersion calculation, the boundary value problem given by is considered, with D i representing the i-derivative of the function to be found.The domain is x ∈ 0, π ; the boundary conditions are y 0 y π y 0 0 and y π −π, with analytical solution being y x x sin x.Transposing the domain to the interval −1, 1 , by using the methodology described in 3, Sections 3.3 and 5 , 3.2 becomes with boundary conditions being y t −1 y t 1 y t −1 0, and y t 1 −π 2 /2.Then, isolating the major order derivative, choosing an order 12 Legendre series given in Table 1 from a first interaction of the SIV algorithm, described in 3, appendix B , as a solution candidate, and applying the same SIV procedure 3 to integrate it four times, the function shown in Figure 1 is obtained.After the integrations, the coefficients of the result are compared with the coefficients of the candidate and the dispersion for each coefficient is determined.In spite of not having a statistical meaning, the term dispersion is used, for the sake of simplicity.
It is worth noticing that, for the Legendre basis 1 −1 P k u 2 du 2/ 2k 1 , and, consequently, increasing the order, this term decreases meaning that higher-order coefficients with the same absolute errors present greater relative errors.
Consequently, the dispersion to be exhibited follows the following relation: Figure 2 shows the relative dispersion for the coefficients of the approximated solution presented in Example 1.This dispersion calculation is used to estimate the error upper bounds, as shown in the following.

Solution Error Upper Bound
Considering f u given by a Legendre series approximation f * u , up to order n, one can write f u ∞ k 0 c k P k and f * u n k 0 c k δ k P k , with δk being the error affecting c k .Consequently, the total error up to order n is given by

3.5
Considering the worst scenario, all δ k will contribute to increase the error exactly where P k is maximum.As a matter of fact, some increase the error and some decrease it.In any case, the error upper bound is given by 3.5 and the maximum value of P k , for any k, is 1.
Consequently, the maximum error in the polynomial approximation is Mathematical Problems in Engineering

Derivative Error Upper Bounds
For each derivative, the error bound calculation is analogous to the procedure for the solution error.For the j-order derivatives, the error can be written as δf j u f * j u −f j u , for j from 1 to the differential equation order.
This calculation can be made by using the differentiation matrix M LD 2 .Denoting the exponent of the differentiation matrix by the number of times that was applied to the coefficients vector, for j from 1 to the differential equation order.Considering P k,max 1,

3.8
The error upper bounds in the derivatives are slightly greater than those in the solutions, because the derivatives are expressed by lower order polynomials.In spite of this, the error margins are very low.As a matter of illustration, if the solution of order 7 is used, the error upper bounds for the first, second, and third derivatives are,

Transposing Error Upper Bounds
The expressions for the error upper bounds developed in the former sections are for the work domain, but it is important to transpose them for the domain of the original equation.The domain transposition method developed in 1, 3 can be applied to the errors.
For instance, in order to transpose the first-order derivative error, it can be written as D 1 D 1t λ, with du/dx λ and, therefore δ 1 δ 1t λ.

3.10
For the second-order derivative,

3.11
For superior-order derivatives, the reasoning is analogous.

Finishing Example 1
For Example 1, the upper bounds of the errors for the Legendre series in the interval −1, 1 can now be calculated by using the transposing procedures described previously and are shown in Table 2. Additionally, Figure 3 shows the local errors along the interval.Again, the SIV procedure described in 3 is used, in order to search the solutions.

Example 2
In the following example, the tools described here and in 1-3 will be applied to a boundary value problem, running in 4 GB RAM processor, and taking 20 s for the whole processing.The differential equation to be solved is 1 x 2 3 arctan x 0, 3.12 with x ∈ 0, 1 and y 0 y 0; y 1 π 2 /16 and y 1 π/4, with D i representing the i-derivative of the function to be found.
Under these conditions, the analytical solution is y arctan 2 x , and, after 100 generations of the SIV algorithm described in 3, appendix B , in the Jacobi polynomial domain, the dispersion of the coefficients is shown in Figure 4.The residual values when the approximated solution is replaced in 3.12 are shown in Figure 5.
The error upper bounds are given in Table 3, the obtained solution is f * x arctan 2 x 3.5x10 −17 arctan x 5.7x10 −16 arctan 3 x 2.8x10 −18 , 3.13 and considering that the analytical solution is given, the local error for the interval can be found and is shown in Figure 6.

Conclusions
Since the seminal work of Chen and Hsiao 22 , the research about operational matrices and the computational capacity advances permitted a strong development in numerical solutions for differential equations.The SIV algorithm developed in 1, 3 gives an interesting combination of operational matrices and computational techniques, by using polynomial approximations for the solutions.

Mathematical Problems in Engineering
Error in D 0     Several advantages are inherent to the method: i the basis of the Hilbert spaces is complete and the solutions are square integrable; consequently, increasing the order of the series, the approximation is improved; ii as the elements of the basis are orthogonal, changing the coefficients of a component does not interfere with the precision of the other components; iii the learning process of the algorithm reduces the search space, reducing the computational costs and decreasing the error margins; iv the error bounds can be determined in a very simple way.
j , and the symbol • | • represents the internal product in the Hilbert space 3 .Defining f | B k b a f • B k w u du with w u being the nonnegative weight function, considering the interval −1, 1 and w u

Figure 2 :
Figure 2: Relative dispersion of the coefficients Example 1 .
representing the coefficient errors.

Figure 4 : 14 Error
Figure 4: Relative dispersion of the coefficients Example 2 .

Table 2 :
Error upper bounds for Example 1.

Table 3 :
Error upper bounds for Example 2.