New Predictor-Corrector Methods with High Efﬁciency for Solving Nonlinear Systems

A new set of predictor-corrector iterative methods with increasing order of convergence is proposed in order to estimate the solution of nonlinear systems. Our aim is to achieve high order of convergence with few Jacobian and/or functional evaluations. Moreover, we pay special attention to the number of linear systems to be solved in the process, with di ﬀ erent matrices of coe ﬃ cients. On the other hand, by applying the pseudocomposition technique on each proposed scheme we get to increase their order of convergence, obtaining new e ﬃ cient high-order methods. We use the classical e ﬃ ciency index to compare the obtained procedures and make some numerical test, that allow us to conﬁrm the theoretical results.


Introduction
Many relationships in nature are inherently nonlinear, which according to these effects are not in direct proportion to their cause.Approximating a solution ξ of a nonlinear system, F x 0, is a classical problem that appears in different branches of science and engineering see, e.g. 1 .In particular, the numerical solution of nonlinear equations and systems is needed in the study of dynamical models of chemical reactors 2 or in radioactive transfer 3 .Moreover, many of numerical applications use high precision in their computations; in 4 , high-precision calculations are used to solve interpolation problems in astronomy; in 5 the authors describe the use of arbitrary precision computations to improve the results obtained in climate simulations; the results of these numerical experiments show that the high-order methods associated with a multiprecision arithmetic floating point are very useful, because it yields a clear reduction in iterations.A motivation for an arbitrary precision in interval methods can be found in 6 , in particular for the calculation of zeros of nonlinear functions.
Recently, many robust and efficient methods with high convergence order have been proposed to solve nonlinear equations, but in most of cases the schemes cannot be extended to multivariate problems.Few papers for the multidimensional case introduce methods with high order of convergence.The authors design in 7 a modified Newton-Jarrat scheme of sixth order; in 8 a third-order method is presented for computing real and complex roots of nonlinear systems; Shin et.al. compare in 9 Newton-Krylov methods and Newton-like schemes for solving big-sized nonlinear systems; the authors in 10 and A. Iliev and I. Iliev in 11 show general procedures to design high-order methods by using frozen Jacobian and Taylor expansion, respectively.Special case of sparse Jacobian matrices is studied in 12 .
Dayton et al. in 13 formulate the multiplicity for the general nonlinear system at an isolated zero.They present an algorithm for computing the multiplicity structure, propose a depth-deflation method for accurate computation of multiple zeros, and introduce the basic algebraic theory of the multiplicity.
In this paper, we present three new Newton-like schemes, of order of convergence four, six, and eight, respectively.After the analysis of convergence of the new methods, we apply the pseudocomposition technique in order to get higher-order procedures.This technique see 14 consists of the following: we consider a method of order of convergence p as a predictor, whose penultimate step is of order q, and then we use a corrector step based on the Gaussian quadrature.So, we obtain a family of iterative schemes whose order of convergence is min{q p, 3q}.This is a general procedure to improve the order of convergence of known methods.
To analyze and compare the efficiency of the proposed methods we use the classic efficiency index I p 1/d due to Ostrowski 15 , where p is the order of convergence and d is the number of functional evaluations at each iteration.
The convergence theorem in Section 2 is demonstrated by means of the n-dimensional Taylor expansion of the functions involved.Let F : D ⊆ R n → R n be sufficiently Frechet differentiable in D. By using the notation introduced in 7 , the qth derivative of F at u ∈ R n , q ≥ 1, is the q-linear function F q u v 1 , . . ., v q , for all permutation σ of {1, 2, . . ., q}.
So, in the following we will denote: It is well known that, for ξ h ∈ R n lying in a neighborhood of a solution ξ of the nonlinear system F x 0, Taylor's expansion can be applied assuming that the jacobian matrix F ξ is nonsingular , and where C q 1/q! F ξ −1 F q ξ , q ≥ 2. We observe that C q h q ∈ R n since

Journal of Applied Mathematics 3
In addition, we can express the Jacobian matrix of F, F as where I is the identity matrix.Therefore, qC q h q−1 ∈ L R n .From 1.2 , we obtain where , is called the error equation and p is the order of convergence.
The rest of the paper is organized as follows: in the next section, we present the new methods of order four, six, and eight, respectively.Moreover, the convergence order is increased when the pseudocomposition technique is applied.Section 3 is devoted to the comparison of the different methods by means of several numerical tests.

Design and Convergence Analysis of the New Methods
Let us introduce now a new Jarratt-type scheme of five steps which we will denote as M8.We will prove that its first three steps define a fourth-order scheme, denoted by M4, and its four first steps become a sixth-order method that will be denoted by M6.The coefficients involved have been obtained optimizing the order of the convergence, and the whole scheme requires three functional evaluations of F and two of F to attain eighth order of convergence.Let us also note that the linear systems to be solved in first, second, and last step have the same matrix and also have the third and fourth steps, so the number of operations involved is not as high as it can seem.Theorem 2.1.Let F : Ω ⊆ R n → R n be a sufficiently differentiable in a neighborhood of ξ ∈ Ω which is a solution of the nonlinear system F x 0, and let x 0 be an initial estimation close enough to the solution ξ.One also supposes that F x is continuous and nonsingular at ξ.Then, the sequence {x k } k≥0 obtained by converges to ξ with order of convergence eight.The error equation is

2.2
Proof.From 1.1 and 1.2 we obtain

2.3
As where X 1 I and X s − s j 2 jX s−j 1 C j , for s 2, 3, . . .So, where where

2.7
We also obtain the Taylor expansion of F x k − 3F y k : where where where , and the error equation of the method at this step is

2.12
So, where N s L s − 1/2 M s and s 6, 7, . . .This shows that the first four steps of the method define a sixth-order scheme.Indeed, and the error equation is proving that the order of convergence of the analyzed method is eight.
In 14 the authors presented a new procedure to design higher-order schemes.This technique, called pseudocomposition, uses the two last steps of the predictor method to obtain a corrected scheme with higher order of convergence.Theorem 2.2 see 14 .Let F : Ω ⊆ R n → R n be differentiable enough Ω, let ξ ∈ Ω be a solution of the nonlinear system F x 0, and let x 0 be an initial estimation close enough to the solution ξ.We suppose that F x is continuous and nonsingular at ξ. Let y k and z k be the penultimate and final steps of orders q and p, respectively, of a certain iterative method.Taking this scheme as a predictor we get a new approximation x k 1 of ξ given by where 1 − τ i y k and τ i , ω i i 1, 2, . . ., m are the nodes and weights of the orthogonal polynomial corresponding to the Gaussian quadrature used.Then, 1 the obtained set of families will have an order of convergence at least q; 2 if σ 2 is satisfied, then the order of convergence will be at least 2q; 3 if, also, σ 1 0, the order of convergence will be min{p q, 3q}, where n i 1 ω i σ and n i 1 ω i τ j i /σ σ j with j 1, 2.
Depending on the orthogonal polynomial corresponding to the Gaussian quadrature used in the corrector step, this procedure will determine a family of schemes.Furthermore, it is possible to obtain different methods in these families by using distinct number of nodes corresponding to the orthogonal polynomial used see Table 1 .However, according to the proof of Theorem 2.2 the order of convergence of the obtained methods does not depend on the number of nodes used.
Let us note that these methods, obtained by means of Gaussian quadratures, seem to be known interpolation quadrature schemes such as midpoint, trapezoidal, or Simpson's method see 16 .It is only a similitude, as they are not applied on the last iteration x k , and the last step of the predictor z k , but on the two last steps of the predictor.In the following, we will use a midpoint-like as a corrector step, which corresponds to a Gauss-Legendre quadrature with one node; for this scheme the order of convergence will be at least min{q p, 3q}, by applying Theorem 2.2.
The pseudocomposition can be applied to the proposed scheme M8 with iterative expression 2.1 , but also to M6.By pseudocomposing on M6 and M8 there can be obtained two procedures of order of convergence 10 and 14 denoted by PsM10 and PsM14 ,  respectively.Let us note that it is also possible to pseudocompose on M4, but the resulting scheme would be of third order of convergence, which is worse than the original M4, so it will not be considered.Following the notation used in 2.1 , the last step of PsM10 is and the last three steps of psM14 can be expressed as

2.24
In Figure 1, we analyze the efficiency indices of the proposed methods, compared with Newton and Jarrat's schemes and between themselves.There can be deduced the following conclusions: the new methods M4, M6, and M8 and also the pseudocomposed PsM10 and PsM14 improve Newton and Jarratt's schemes in fact, the indices of M4 and Jarratt's are equal .Indeed, for n ≥ 3 the best index is that of M8.Nevertheless, none of the pseudocomposed methods improve the efficiency index of their original partners.Nevertheless, as we will see in the following section, the pseudocomposed schemes will show a very stable behavior that makes them worth.

Numerical Results
In order to illustrate the effectiveness of the proposed methods, we will compare them with other known schemes.Numerical computations have been performed in MATLAB R2011a by using variable-precision arithmetic, which uses floating-point representation of 2000 decimal digits of mantissa.The computer specifications are Intel R Core TM i5-2500 CPU @ 3.30 GHz with 16.00 GB of RAM.Each iteration is obtained from the former by means of an iterative expression x k 1 x k − A −1 b, where x k ∈ R n , A is a real matrix n × n and b ∈ R n .The matrix A and vector b are different according to the method used, but in any case, we calculate A −1 b as the solution of the linear system Ay b, with Gaussian elimination with partial pivoting.The stopping criterion used is Firstly, let us consider the following nonlinear systems of different sizes:

3.1
When n is odd, the exact zeros of Table 2 presents results showing the following information: the different iterative methods employed Newton NC , Jarratt JT , the new methods M4, M6, and M8, and the pseudocomposed PsM10 and PsM14 , the number of iterations Iter needed to converge to the solution Sol, the value of the stopping factors at the last step, and the computational order of convergence ρ see 17 approximated by the formula: The value of ρ which appears in Table 2 is the last coordinate of the vector ρ when the variation between their coordinates is small.Also the elapsed time, in seconds, appears in Table 2, being the mean execution time for 100 performances of the method the command cputime of MATLAB has been used .We observe from Table 2 that not only the order of convergence and the number of new functional evaluations and operations are important in order to obtain new efficient iterative methods to solve nonlinear systems of equations.A key factor is the range of applicability of the methods.Although they are slower than the original methods when the initial estimation is quite good, when we are far from the solution or inside a region of instability, the original schemes do not converge or do it more slowly, the corresponding pseudocomposed procedures usually still converge or do it faster.

Figure 1 :
Figure 1: Efficiency index of the different methods for different sizes of the system.

Figure 2 :Figure 3 :
Figure 2: Real dynamical planes for system F 2 x 0 and methods M6 and PsM10.

Figure 4 :Figure 5 :
Figure 4: Real dynamical planes for system F 3 x 0 and methods M6 and PsM10.

Table 1 :
Values of the parameters for the different quadratures used.

Table 2 :
Numerical results for functions F 1 to F 4 .