Modiﬁed Potra-Pták Multi-step Schemes with Accelerated Order of Convergence for Solving Systems of Nonlinear Equations

: In this study, an iterative scheme of sixth order of convergence for solving systems of nonlinear equations is presented. The scheme is composed of three steps, of which the ﬁrst two steps are that of third order Potra-Pták method and last is weighted-Newton step. Furthermore, we generalize our work to derive a family of multi-step iterative methods with order of convergence 3 r + 6, r = 0, 1, 2, . . . . The sixth order method is the special case of this multi-step scheme for r = 0. The family gives a four-step ninth order method for r = 1. As much higher order methods are not used in practice, so we study sixth and ninth order methods in detail. Numerical examples are included to conﬁrm theoretical results and to compare the methods with some existing ones. Different numerical tests, containing academical functions and systems resulting from the discretization of boundary problems, are introduced to show the efﬁciency and reliability of the proposed methods.


Introduction
Many applied problems in Science and Engineering [1][2][3] are reduced to solve nonlinear systems F(x) = 0 numerically, that is, for a given nonlinear function F(x) : D ⊂ R m −→ R m , where F(x) = ( f 1 (x), f 2 (x), ..., f m (x)) T and x = (x 1 , x 2 , ..., x m ) T , to find a vector α = (α 1 , α 2 , ..., α m ) T such that F(α) = 0.The most widely used method for this purpose is the classical Newton's method [3,4], which converges quadratically under the conditions that the function F is continuously differentiable and a good initial approximation x (0) is given.It is defined by where F (x (k) ) −1 is the inverse of Fréchet derivative F (x (k) ) of the function F(x) at x (k) .In order to improve the order of convergence of Newton's method, several methods have been proposed in literature, see, for example [5][6][7][8][9][10][11] and references therein.In particular, the third order method by Potra-Pták [11] for systems of nonlinear equations is given by It is quite clear that this scheme requires the evaluation of two functions, one derivative and one matrix inversion per iteration, that is usually avoided by solving a linear system.This algorithm is illustrious not only for its simplicity but also for its efficient character.
In this paper, based on Potra-Pták method (2), we develop a three-step scheme with increased order of convergence and still maintaining the efficient character.With these considerations, we propose a three-step iterative method with accelerated sixth order of convergence; of the three steps, the first two are those of Potra-Pták method whereas the third is a weighted Newton-step.Then, based on this scheme, a multi-step family with increasing order of convergence 3r + 6, r = 0, 1, 2, . . ., is developed.The sixth order method is the special case of this multi-step scheme for r = 0.The family gives a four-step ninth order scheme for r = 1.As much higher order methods are not used in practice, so we study sixth and ninth order methods in particular.
The rest of the paper is organized as follows: In Section 2, we present the new three-step scheme for solving nonlinear systems and we analyze its order of convergence.In Section 2.1, the order of this scheme is improved in three units, by adding another step that needs a new functional evaluation and to solve a linear system with the same matrix of coefficients as before.This idea can be generalized for designing an iterative method with arbitrary order of convergence.The computational efficiency of the proposed schemes is studied in Section 3, doing a comparative analysis with the efficiency of other known methods.Section 4 is devoted to the numerical experiments with academical multivariate functions and nonlinear systems resulted of the discretization of boundary problems.For some cases, the basins of attraction of the methods are also showed.The paper finishes with some conclusions and the references used in it.

The Method and Analysis of Convergence
We start with the following iterative scheme: where the first two-steps are those of Potra-Pták scheme [11] for nonlinear systems and [z (k) , y (k) ; F] is the first order divided difference of F. In order to discuss the behavior of scheme (3), we consider the following expression of divided difference operator [•, By expanding F (x + th) in Taylor series at the point x and integrating, we have where where Inversion of F (x (k) ) yields, We are in a position to analyze the behavior of scheme (3).Thus, the following result is proven.
Employing Equations ( 10) and (11) in the second step of (3), we get Using Equations ( 7)-( 9) in ( 5) With the help of Equations ( 10) and ( 13), we can write Expanding F(z (k) ) about α, we obtain Using (10), ( 14) and (15) in the last step of (3), we get In order to achieve sixth order of convergence it is clear that terms 1 − a − b − c, a + 2b + 3c and 4a + 13(b + 2c) must vanish for some values of a, b, and c.This happens when a = 13  4 , b = − 7 2 and c = 5 4 .Thus, the error equation ( 16) becomes This proves the sixth order of convergence.
Clearly this formula uses three functional evaluations, one evaluation of the Jacobian matrix, one of a divided difference, and one matrix inversion per iteration.We denote this scheme by H 6,1 .
2.1.Multi-step Method with Order 3r + 6 In this section, we improve the H 6,1 by adding a functional evaluation per each new step to get the multi-step version called H 3r+6,1 method.The method is defined as Let us note that case r = 0 is H 6,1 method given by ( 18).This multi-step version has the order of convergence 3r + 6, r 1, which we shall prove through the following result.
open convex set D containing the zero α of F(x) and F (x) is continuous and nonsingular at α.Then, sequence {x (k) } k 0 , x (0) ∈ D obtained by using method (19) converges to α with convergence order 3r + 6.
As a special case, this family when r = 1 gives a ninth order method: It is clear that this scheme requires four functional evaluations, one Jacobian matrix, one divided difference, and one matrix inversion per full iteration.We denote this scheme by H 9,1 .
Remark 1. Multi-step version H 3r+6 (r 0), utilizes r + 3 functional evaluations of F, one evaluation of F , and one divided difference.Also, the scheme (19) requires only one matrix inversion per iteration.

Computational Efficiency
To obtain an estimation of the efficiency of the proposed methods we shall make use of efficiency index, according to which the efficiency of an iterative method is given by E = ρ 1 C , where ρ is the order of convergence and C is the computational cost per iteration.For a system of m nonlinear equations and m unknowns, the computational cost per iteration is given by (see [9,10]) where P 0 (m) represents the number of evaluations of scalar functions ( f 1 , f 2 , ....., f m ) used in the evaluations of F and [y, x ; F].The divided difference [y, x ; F] of F is an m × m matrix with elements (see [11,12]) The number of evaluations of scalar functions of F , i.e.
, 1 i, j m, is P 1 (m).P(m) represents the number of products or quotients needed per iteration, and µ 0 and µ 1 are ratios between products and evaluations required to express the value of C(µ 0 , µ 1 , m) in terms of products.
To compute F in any iterative method we calculate m scalar functions and if we compute the divided difference [y, x ; F] then we evaluate 2m(m − 1) scalar functions, where F(x) and F(y) are computed separately.We must add m 2 quotients from any divided difference.The number of scalar evaluations is m 2 for any new derivative F .In order to compute an inverse linear operator we solve a linear system, where we have m(m − 1)(2m − 1)/6 products and m(m − 1)/2 quotients in the LU decomposition and m(m − 1) products and m quotients in the resolution of two triangular linear systems.We suppose that a quotient is equivalent to l products.Moreover, we add m 2 products for the multiplication of a matrix with a vector or of a matrix by a scalar and m products for the multiplication of a vector by a scalar.
Denoting the efficiency indices of H ρ,i (ρ = 6, 9 and i = 1, 2, 3, 4) by E ρ,i and computational cost by C ρ,i , then taking into account the above considerations, we obtain ) and E 6,1 = 6 1/C 6,1 . (28) To check the performance of the new sixth order method, we compare it with some existing sixth order which belongs to the same class.So, we choose the sixth order methods presented in [10,13].The sixth order methods presented in [10] are given by and The sixth order method proposed in [13] is given by Per iteration these methods utilize the same number of function evaluations as that of H 6,1 .The computational cost and efficiency for the methods H 6,2 , H 6,3 and H 6,4 is given as below: ) and E 6,2 = 6 1/C 6,2 , (33)

Comparison among the Efficiencies
To compare the iterative methods H ρ,i , we consider the ratio It is clear that if R p,i;q,j > 1, the iterative method H p,i is more efficient than H q,j .Moreover, if we need to compare the methods having the same order then using (36) we can say that the iterative method H p,i is more efficient than H q,j , if C q,j > C p,i .This means that the comparison of the methods possessing same order can be done by just comparing their computational costs.
In addition to the above comparisons we compare the proposed methods H 6,1 and H 9,1 with each other.H 9,1 versus H 6,1 case: The particular boundary R 9,1;6,1 = 1 expressed by µ 0 written as a function of µ 1 and m is where r = ln(3) and s = ln(2).This function has a vertical asymptote for m = s/(r − s) = 1.7095.Note that for m 44, the numerator of (37) is positive and the denominator is negative, which shows that µ 0 is always negative for m 44.That is, the boundary is out of admissible region for m 44 and we have E 9,1 > E 6,1 ∀ (µ 1 , µ 0 ) ∈ (0, +∞) × (0, +∞) and l 1.
We summarize the above results in following theorem: Theorem 3.For all µ 0 > 0, µ 1 > 0 and l 1 we have: Otherwise, the efficiency comparison depends on m, µ 0 , µ 1 and l.

Numerical Results
This section is devoted to check the effectiveness and efficiency of some of our proposed methods on different types of applications.In all cases, we apply our proposed schemes H 6,1 and H 9,1 and compare the results with those obtained by known methods H 6,2 , H 6,3 , and H 6,4 .We consider academic examples, a special case of a nonlinear conservative problem which is transformed in a nonlinear system by approximating the derivatives by divided differences and, also the approximation of the solution of an elliptic partial differential equation that model a chemical problem.
All the experiments have been carried out in Matlab 2017 with variable precision arithmetics with 1000 digits of mantissa.These calculations have been made with an Intel Core processor i7-4700HQ with a CPU of 2.40GHz and 16.0 GB al RAM memory.In the tables we include the number of iterations (iter), the residual error of the corresponding function ( F(x (k+1) ) ) in the last iterated and the difference between the two last iterates ( x (k+1) − x (k) ).We also present the approximated computational order of convergence (ACOC) defined in [14] with the expression When the components of this vector are stable, it is an approximation of the theoretical order of convergence.In other case, it does not give us any information and we will denote it by −.Also the execution time (in seconds) has been calculated (by means of "cputime" Matlab command) with the mean value of 100 consecutive executions, for big-sized systems corresponding to Examples 1 to 3.

Example 1
We consider the case of a nonlinear conservative problem described by the differential equation with the boundary conditions y(0) = y(1) = 0.
We transform this boundary problem into a system of nonlinear equations by approximating the second derivative by a divided difference of second order.We introduce points t i = 0 + ih, i = 0, 1, . . ., m + 1, where h = 1 m+1 and m is an appropriate positive integer.A scheme is then designed for the determination of numbers y i as approximations of the solution y(t) at point t i .By using divided differences of second order we transform the boundary problem into the nonlinear system Introducing vectors y = (y 1 , y 2 , . . ., y m ) T , Φ y = (Φ(y 1 ), Φ(y 2 ), . . ., Φ(y m )) T and the matrix of size m × m 38) can be written in the form In this case, we choose the law Φ(y(t)) = 1 + y(t) 3 for the heat generation in the boundary problem and we solve system (38) by using iterative methods H 6,1 , H 6,2 , H 6,3 , H 6,4 , and H 9,1 .In all cases, we use as initial estimation x (0) = (0.5, 0.5, . . ., 0.5) T and the stoping criterium x (k+1) − x (k) < 10 −100 or F(x (k+1) ) < 10 −100 .We can see in Table 1 the results obtained for m = 20 and m = 50.There are no significant differences among the results obtained for different values of step h, i.e. for different sizes of the nonlinear system resulting from discretization.Let us remark that, for the lowest number of iterations, the best results in terms of lowest residuals have been obtained by methods H 9,1 and H 6,1 .

Example 2
Let us consider the system of nonlinear equations m ∑ j=1,j =i for an arbitrary positive integer m.We solve this system whit the same schemes as before, using as initial guess x (0) = (1, 1, . . ., 1) T and two values for the size, m = 20 and m = 50, being the solution α = (0.05, 0.05, . . ., 0.05) T and α = (0.02, 0.02, . . ., 0.02) T , respectively.The stopping criterium used is again x (k+1) − x (k) < 10 −100 or F(x (k+1) ) < 10 −100 , that is, the process finishes when one of them is satisfied.The obtained results are shown in Table 2.In this case, all the methods use the same number of iterations to achieve the solution, but the lowest residual are those of proposed schemes H 9,1 and H 6,1 .

Example 3
Gas dynamics can be modeled by a boundary problem described by the following elliptic partial differential equation and the boundary conditions: By using central divided differences and step h = 1/5 in both variables, this problem is discretized in the nonlinear system To solve problem (39) we have used as initial guess x (0) = (1, 1, ..., 1) T .Also the stopping criteria F(x (k+1) ) < 10 −100 or x (k+1) − x (k) < 10 −100 have been used and the process finishes when one of them is satisfied or the number of iterations reachs to 50.
The results are shown in Table 3.The first column shows the numerical aspects (approximated computational order of convergence ACOC, last difference between consecutive iterates x (k+1) − x (k) and residual F(x (k+1) ) ) analyzed for the schemes used to solve the problem.In the rest of the columns we show the numerical results obtained by the methods H 6,1 , H 6,2 , H 6,3 , H 6,4 , and H 9,1 .The low number of iterations needed justify that the ACOC does not approximate properly the theoretical order of convergence.The methods giving lowest exact error at the last iteration are again H 9,1 and H 6,1 .

Example 4
Let us consider the nonlinear two-dimensional system F(x 1 , x 2 ) = 0 with coordinate functions This system has two real solutions at points, approximately, (0.954811, 0.301815) T and (−0.954811, −0.301815) T .By using x (0) = (1, 0.5) T as initial estimation, we have obtained the results appearing in Table 4, in which the residuals x (k+1) − x (k) and F(x (k+1) ) appear, for the three first iterations.It can be observed that, although the difference between two consecutive iterations is not very small, the precision in the estimation of the root is very high, being the best for proposed schemes H 6,1 and H 9,1 .
Moreover, dynamical planes help us to get global information about the convergence process.In Figure 1 we can see the dynamical planes of the proposed methods H 6,1 and H 9,1 and known schemes H 6,2 , H 6,3 , and H 6,4 when they are applied on system F(x 1 , x 2 ) = 0.This figures are obtained by using the routines described in [15].To draw these images, a mesh of 400 × 400 initial points has been used, 80 was the maximum number of iterations involved and 10 −3 the tolerance used as the stopping criterium.In this paper, we have used a white star to show the roots of the nonlinear system.A color is assigned to each initial estimation (each point of the mesh) depending on where they converge Blue and orange correspond to the basins of the roots of the system (40) (brighter as lower is the number of iterations needed to converge) and black color is adopted when the maximum number of iterations is reached or the process diverges.In Figure 1, the shape and wideness of the basins of attraction show that H 6,1 , H 6,2 , H 6,3 , and H 6,4 can found any of both roots by using a great amount of initial estimations, some of them far from the roots.However, H 6,4 is hardly able to find the roots and their basins are very small.

Example 5
Let us consider the nonlinear two-dimensional system G(x 1 , x 2 ) = 0 with coordinate functions This system has four real solutions at points 1 2 , T , and T .Let us remark that they are symmetric as the system shows the intersection between two conical curves, an ellipse and a circumference.By using x (0) = (1, 1) T as initial estimation, we have obtained the results appearing in Table 5, in which the residuals x (k+1) − x (k) and F(x (k+1) ) also appear, for the three first iterations.The numerical results appearing in Table 5 show as the estimation of the roots has similar errors in all sixth order methods, although H 6,2 and H 6,3 are slightly better than H 6,1 .Indeed, the ninth order scheme has the best precision in the approximation of the roots of G(x 1 , x 2 ), as corresponds to its high order of convergence.In terms of computational time, the mean of one hundred consecutive iterations have been used in order to get good estimations of the real time and, in case of proposed methods, they use similar times as the existing schemes getting better accuracy.
In Figure 2 the colors blue, orange, green and purple correspond to the basins of the roots of the system (41).We can see the dynamical planes of the proposed method methods H 6,1 and H 9,1 and known schemes H 6,2 , H 6,3 , and H 6,4 when they are applied on system G(x 1 x 2 ) = 0. Regarding the stability of proposed and known iterative processes on G(x 1 , x 2 ), it can be observed in Figure 2 that, in all cases except H 6,4 , the connected component of the basins of attraction that holds the root (usually known as immediate basin of attraction) is very wide.For all iterative methods black areas of low convergence or divergence appear, being more diffuse in case of H 6,2 and wider in H 6,4 .Methods H 6,1 , H 6,3 and H 9,1 present black areas that are mainly regions of slower convergence, that would be colored if a higher number of iterations is fixed as a maximum; meanwhile, the biggest black area of the dynamical plane associated to method H 6,4 corresponds to divergence.

Concluding Remarks
From Potra-Pták third order scheme, an efficient sixth order method for solving nonlinear system is proposed.Moreover, it has been extended by adding subsequent steps with the same structure that adds a new functional evaluation per new step.Then the order of the resulting procedure increases in three units per step.The proposed methods have been shown to be more efficient that several known methods of the same order.Some numerical tests with academic and real-life problems have been made to check the efficiency and applicability of the designed procedures.The numerical performance have been opposed to the stability shown by the dynamical planes of new and known methods on two-dimensional systems.

Theorem 1 .
Let F : D ⊂ R m → D be sufficiently differentiable function in an open neighborhood D of its zero α.Let us suppose that the Jacobian matrix F (x) is continuous and nonsingular at α.If an initial approximation x (0) is sufficiently close to α, then the local order of convergence of method (3) is at least 6, provided a = 13/4, b = −7/2 and c = 5/4.Proof.Let e (k) y = y (k) − α is the local error of Newton's method given by e (k)

Table 1 .
Numerical results for conservative boundary problem.