Design of High-Order Iterative Methods for Nonlinear Systems by Using Weight Function Procedure

We present two classes of iterative methods whose orders of convergence are four and five, respectively, for solving systems of nonlinear equations, by using the technique of weight functions in each step. Moreover, we show an extension to higher order, adding only one functional evaluation of the vectorial nonlinear function. We perform numerical tests to compare the proposed methods with other schemes in the literature and test their effectiveness on specific nonlinear problems. Moreover, some real basins of attraction are analyzed in order to check the relation between the order of convergence and the set of convergent starting points.


Introduction
To find the solution of systems of nonlinear equations ( ) = 0, : R → R , is a common, ancient, and important problem in mathematics and engineering (see, e.g., [1]). Previous experimentation on some of these applied problems shows that high-order methods, associated with floating point arithmetics of multiprecision, are very useful because they carry a clear reduction in the number of required iterations (see, e.g., [2,3]).
Over the last years there have been numerous contributions of different authors that have designed iterative methods trying to solve nonlinear systems, making several modifications to the classical schemes to accelerate the convergence and reduce the number of operations and functional evaluations at each step of the iterative process. The extension of the variants of Newton's method for scalar problems described by Weerakoon and Fernando in [4], byÖzban in [5], and by Gerlach in [6], to functions of several variables has been developed in [7][8][9][10]. All these processes have yield generating multipoint methods with Newton's method as a predictor. However, a general procedure called pseudocomposition is designed in [11,12] by using a generic predictor and the Gaussian quadrature as the corrector.
On the other hand, one of the most used techniques to generate high-order methods for solving nonlinear equations or systems was introduced by Traub in [13], that is, the composition of two iterative methods of orders 1 and 2 , respectively, to obtain a method of order 1 2 . With the direct application of this technique, it is usually necessary to evaluate the nonlinear function and its associated Jacobian matrix to increase the order of convergence and the new method usually ends up being less efficient than the original ones.
So, for designing a more efficient method, it is usual to estimate the new evaluation of the Jacobian matrix, in terms of Jacobian matrices previously used. With this in mind, the new method generally has a lower order of convergence than It is known that this method has second order of convergence and it uses one functional evaluation of the vectorial function ( ) and also of the associated Jacobian matrix ( ) at the previous iterate, to generate a new one.
In the numerical section, we will use some multipoint iterative schemes to compare with ours. These have been developed by different authors: the first one is the fourth-order method designed by Jarratt in [14], whose iterative expression is In the following, we will denote it by JM4. Let us also consider the fourth-order method designed by Sharma et al. in [15] whose iterative formula is denoted by SHM4. We will also use some fifth-order methods to compare with ours: the scheme obtained by Vassileva [12], generalizing Jarratt's method, improves this convergence by adding functional evaluations (one functional evaluation of in the second step of the process), being its iterative expression: where = 1, 1 = 1, 2 = 5, 1 = −3, and 2 = −1 and denoted by VM5, and the method designed by Sharma and Gupta in [16] whose iterative formula is where = 2 and = −1. We denote this method as SHM5.
In this work we have designed some iterative methods with orders of convergence four and five and presented an analysis in Section 2. We compare them in Section 3 with other known processes on academic examples, with the aim that the results of this comparison will give us an idea of how robust our methods are. The stability is studied under the point of view of the real dynamics on an specific 2-dimensional problem. An applied problem, the partial differential equation of molecular interaction, is used to check the conclusions reached in previous sections.

Design and Analysis of Convergence
Let : Ω ⊆ R → R , > 1, be a sufficiently differentiable function on a convex set Ω ∈ R and let be a solution of the nonlinear system of equations ( ) = 0. We propose the following iterative scheme: where with 1 , 2 , and 3 being real parameters. This scheme has been designed as a modified double-Newton method with a frozen evaluation of and and a matrix weight function ( ) at the second step. We denote this class of methods (whose order of convergence is four under certain conditions) by MS4. As it has been stated, is a matrix function with matrix variable. Specifically, if = × denotes the Banach space of real square × matrices, then we can define : → such that its Frechet derivatives satisfy Let us note that L( ) denotes the space of linear mappings of on itself.
To prove the local convergence of the proposed iterative method to the solution of ( ) = 0, we will use a Taylor series expansion of the involved functions around the solution, by assuming that Jacobian matrix ( ) is nonsingular in a neighborhood of , Abstract and Applied Analysis 3 So, can be expressed (in a neighborhood of ) as where is the identity matrix. This technique, including the notation used, has been detailed in [17].
Finally, the error equation is So, the proof is completed.
Starting from the general family of iterative methods (6) and under the hypothesis of Theorem 1, we develop two specific fourth-order iterative schemes.
(i) By requiring 1 = 3 = 0 and 2 = 1, it is directly obtained from Theorem 1 that ( ) = (1/3) and . So, we get the following iterative expression denoted by MS4A: (ii) The second iterative scheme is denoted by MS4B and is defined by imposing the following conditions: 1 = 1/4, 2 = 0, and 3 = 3/4 so that Theorem 1 implies that ( ) = (1/3) , 1 = −1, and 2 = 4. The resultant iterative expression that we call MS4B is At this point we wonder if we can obtain a method of order of convergence higher than four slightly modifying the iterative scheme (6), just adding one functional evaluation of in its second step: Abstract and Applied Analysis We denote this class of methods by MS5. In the following, we prove that its order of convergence is five under certain conditions. One supposes that ( ) is continuous and nonsingular at . Let where is the identity matrix. Let one consider (0) as an initial guess, sufficiently close to . Then, the sequence { ( ) } ≥0 obtained using the expression (30) converges to with order of convergence five if = 1 and 1 + 2 + 3 = 1. Moreover, the error equation is Proof. Let us remark that the only difference between methods (30) and (6) is that, in its second step, ( ( ) ) is replaced by ( ( ) ). So, we need the Taylor series expansion of Then, we obtain the following expressions: where 1 = 1 , = + ∑ =1 +1 − , 1 = ℎ 1 1 , and = ∑ =1 ℎ − +1 , , = 2, 3, . . .. So, the expression of the error equation in the last step is where 1 = − 1 and = − − , = 2, 3, . . .. Forcing 1 = 0, 2 = 0, 3 = 0, and 4 = 0, we obtain order of convergence five. Solving this system of equations we obtain the conditions that guarantee order of convergence five as they appear in the hypothesis of this theorem. Then, the error equation takes the form and the proof is completed.
From the family of iterative methods (30) under the conditions imposed in Theorem 2, we get two particular fifthorder iterative schemes.
(i) By requiring that 1 = 3 = 1/2 and 2 = 0, conditions from Theorem 2 yield ( ) = , 1 = −2, and 2 = 10. The expression of the resulting iterative method that we call MS5A is (ii) The second scheme will be denoted by MS5B and results from setting the parameters 1 = 1/4, 2 = 0, and 3 = 3/4 in Theorem 2. Then, ( ) = , 1 = 1, and 2 = 4. Then the iterative expression is It is interesting that, based on Theorem 2, we can generate methods whose order of convergence is higher than 5. This process consists of adding a new step to scheme (30), keeping the weight function ( ( ) ) unaltered. In this way, the iterative expression is Assuming that hypotheses of Theorem 2 are satisfied, the error in the second step of (30) can be expressed as = ( ) − = 5 + 6 + O( 7 ). Then, the error equation of (39) takes the form and the iterative method has order of convergence five. If = 1, the order of convergence of the iterative scheme (39) is seven.
Based on this result, we state the following theorem which generalizes this procedure.

Theorem 3.
Under the same hypothesis of Theorem 2 and adding = 1, sequence { ( ) } ≥0 obtained using the following expression converges to with order of convergence + 2: , and the order of convergence of the penultimate step is with error equation = ( )− = + +1 +O( +2 ).
Proof. In general, the error equation of the last step takes the form This shows that the order of convergence of the method will be + 2 if = 1.
In order to calculate the different efficiency indices, it is necessary to take into account not only the order of convergence, but also the number of functional evaluations ( 2 : each Jacobian matrix and : each vectorial function ) per iteration and the amount of product-quotients per step. This is obtained by observing the iterative expression of the method and using the fact that the number of productquotients needed for solving a set of linear systems with the same coefficient matrix is (1/3) 3 + 2 − (1/3) and the product between a matrix and a vector implies 2 productsquotients.
With respect to the classical efficiency index, it is clear that one of all (both proposed and known) fourth-order methods is 4th = 4 1/(2 2 + ) and, for all fifth-order schemes, 5th = 5 1/(2 2 +2 ) . It is also easy to check that 4th < 5th , for any > 2.
In Table 1, the computational efficiency indices are showed, and also the needed information is showed: the number of functional evaluations of matrices (FEM, 2 evaluations) or nonlinear functions (FEF, evaluations), amount of linear systems with different Jacobian matrices (NS1, (1/3) 3 + 2 − (1/3) products-quotients), the number of linear systems with the same Jacobian matrix (NS2, 2 products-quotients), and also how many times a matrix-vector product appears (MxV, 2 products).  Figure 1 shows the computational efficiency index, for different sizes of the system, of fourth-(MS4A and MS4B) and fifth-order (MS5A and MS5B) proposed methods, joint with the known ones (JM4, SHM4, VM5, and SHM5), for comparison purposes. It can be observed that, in both cases, known methods have better computational efficiency indices for small sizes of the system, but in larger cases the new methods show better behavior.

Numerical Results
In this section we test the developed methods to check their effectiveness compared with other known ones. Here we can see the list of systems of nonlinear equations used in these numerical tests, conducted in Matlab R2010a by using variable precision arithmetics with 2000 digits of mantissa and ‖ ( +1) − ( ) ‖ < 10 −500 or ‖ ( ( +1) )‖ < 10 −500 as a stopping criterion. Consider We have two comparative tables in which iterative methods are grouped according to the order of convergence of the involved schemes: in Table 2, we can observe the behavior of the fourth-order methods and in Table 3, the schemes of order 5 can be found. We perform numerical tests for each of the selected systems of nonlinear equations either for the proposed methods MS4A, MS4B, MS5A, and MS5B or for consolidated methods JM4, SHM4, VM5, and SHM5.
The columns of the tables correspond from left to right to nonlinear system to be solved with the initial approximation, iterative method to be used, solutions or roots found, absolute value of the difference between the two last iterations (component by component), | ( +1) − ( ) |, absolute value of each coordinate function evaluated at the last iteration, | ( ( +1) ) |, the number of iterations used in the process, and approximated computational order of convergence, ACOC, according to (see [8]) Let us remark that the value of ACOC that is presented in Tables 2 and 3 is the last coordinate of vector ACOC when the variation between its values is small. Observing Tables 2 and 3, we note that the approximated computational order of convergence confirms the theoretical results. The results in terms of accuracy in the approximation of the roots and number of iterations are kept in the range expected for the orders of convergence of the methods, as is apparent from a comparison with established methods. We Table 2: Numerical results obtained with the fourth-order schemes.
Methods believe that, in general, the proposed methods are competitive on each of the systems of nonlinear equations used. With respect to the extension to high-order methods, let us see in the following how the increasing order does not reduce the real region of good starting points, as it is usual in iterative methods. We have selected a rectangular region of R 2 that contains some of the solutions of system 1 ( ) = 0 and have used the points of this region as starting ones for the fifth-order iterative methods and their partners of seventh and ninth order.
So, in Figure 2, the dynamical planes associated with known and proposed methods on 1 ( ) are showed. These planes have been generated by slightly modifying the routines described in [18]. In them, a mesh of 400 × 400 points has been used, 80 has been the maximum number of iterations involved, and 10 −3 has been the tolerance used as a stopping criterium. Then, if a starting point of this mesh converges to one of the solutions of the system (marked with white stars), it is painted in the color assigned to the root which it has converged to. The color used is brighter when the number of iterations is lower. If it reaches the maximum number of iterations without converging to any of the roots, it is painted in black. At the sight of Figure 2, it can be concluded that the areas of converging starting points remain with slight variations, in spite of the increasing order of convergence; this makes this process especially interesting, as it can increase the speed of convergence to a solution with no need of being closer to it.

Molecular Interaction Problem.
To solve the equation of molecular interaction (see [19]) we need to deal with a boundary value problem with a nonlinear partial differential equation of second order. To estimate its solution numerically, we have used central divided differences in order to transform the problem in a nonlinear system of equations, which is solved by using the proposed methods (of orders four and five) and the extensions of family MS5 up to order nine. The discretization process yields to the nonlinear system of equations, where , denotes the estimation of the unknown ( , ), = ℎ and = (with = 0, 1, . . . , and = 0, 1, . . . , ) are the nodes in both variables, where ℎ = 1/ and = 1/ .
In this case, we fix = = 4, so a mesh of 5 × 5 is generated. As the boundary conditions give us the value Table 4: Numerical tests for (0) = (1, . . . , 1) .  Method So, the system can be expressed as and is the 3 × 3 identity matrix and = (7/4, 1, 27/8, 1, 0, 2, 27/8, 2, 4) . Now, we will check the performance of the methods by means of some numerical tests, by using variable precision arithmetics of 2000 digits of mantissa. In Tables 4 and 5, we show the numerical results obtained for the problem of molecular interaction (45), with different initial estimations. We show, for the first three iterations, the residual of the function at the last iteration, ‖ ( ( +1) )‖, and the difference between the last iteration and the preceding one ‖ ( +1) − ( ) ‖. We can observe in Tables 4 and 5 that all the new methods converge to the solution of the problem that appears in Table 6. It can be noticed that the error of the test is lower when the order of the iterative method is higher, even at the first iterations. Indeed, if the initial estimation is far from the solution, the proposed methods converge with reasonable results. This is especially important as in real problems, where good initial estimations are not always known.

Conclusions
As far as we know, the weight functions procedure has been used only for designing iterative schemes to solve nonlinear equations. In this paper, by using matrix functions, the mentioned procedure is applied to obtain iterative methods for solving nonlinear systems. Fourth-and fifth-order methods are obtained and a technique for designing iterative schemes of any order is presented.
By using different academic test problems and the discretization of a partial differential equation modeling the molecular interaction problem, we compare our methods with several known ones such as Jarratt's method and Sharma's method, some of them optimal in the context of nonlinear equations. The numerical tests confirm the theoretical results and show that our methods are more competitive than those used for comparing. In addition, the calculus of the computational efficiency index of the different schemes allows us to Finally, with respect to the extension to high-order methods, we show that the increasing order does not reduce the real region of good starting points, as it is usual in iterative schemes.