Iterative Methods with Memory for Solving Systems of Nonlinear Equations Using a Second Order Approximation

: Iterative methods for solving nonlinear equations are said to have memory when the calculation of the next iterate requires the use of more than one previous iteration. Methods with memory usually have a very stable behavior in the sense of the wideness of the set of convergent initial estimations. With the right choice of parameters, iterative methods without memory can increase their order of convergence signiﬁcantly, becoming schemes with memory. In this work, starting from a simple method without memory, we increase its order of convergence without adding new functional evaluations by approximating the accelerating parameter with Newton interpolation polynomials of degree one and two. Using this technique in the multidimensional case, we extend the proposed method to systems of nonlinear equations. Numerical tests are presented to verify the theoretical results and a study of the dynamics of the method is applied to different problems to show its stability.


Introduction
This paper deals with iterative methods for approximating the solutions of a nonlinear system of n equations and n unknowns, F(x) = 0, where F : D ⊆ R n −→ R n is a nonlinear vectorial function defined in a convex set D. The main aim of this work is to design iterative methods with memory for approximating the solutions ξ of F(x) = 0. An iterative method is said to have memory when the fixed point function G depends on more than one previous iteration, that is, the iterative expression is x k+1 = G(x k , x k−1 , . . . ).
A classical iterative scheme with memory for solving scalar equations f (x) = 0, n = 1, is the known secant method, whose iterative expression is with x 0 and x 1 as initial estimations. This method can be obtained from Newton's scheme by replacing the derivative by the first divided difference f [x k , x k−1 ]. For n > 1 the secant method has the expression (k) , x (k−1) ; F] −1 F(x (k) ), k = 1, 2, . . . , (2) where [x (k) , x (k−1) ; F] is the divided difference operator, defined in [1]. It is possible to design methods with memory starting from methods without memory by introducing accelerating parameters. The idea is to introduce parameters in the original scheme and study the error equation of the method to see if any particular values of these parameters allow to increase the order of convergence, see for example [2][3][4] and the references therein. Several authors have constructed iterative methods with memory for solving nonlinear systems, by approximating the accelerating parameters by Newton's polynomial interpolation of first degree, see for example the results presented by Petković and Sharma in [5] and by Narang et al. in [6]. As far as we know, in the multivariate case, higher-degree Newton polynomials have never been used. In this paper, we will approximate the accelerating parameters with multidimensional Newton's polynomial interpolation of second degree.
In [1], Ortega and Rheinboldt introduced a general definition of the order of convergence, called R-order defined as follows where e k = x k − ξ, k = 0, 1, . . ., and R p is called R-factor. The R-order of convergence of an iterative method (IM) at the point ξ is Let (IM) be an iterative method with memory that generates a sequence {x k } of approximations to the root ξ, and let us also assume that this sequence converges to ξ. If there exists a nonzero constant η and nonnegative numbers t i , 0 ≤ i ≤ m such that the inequality holds, then the R-order of convergence of (IM) satisfies the inequality where s * is the unique positive root of the equation We can find the proof of this result in [1]. Now, we introduce some notations that we use in this manuscript. Let r be the order of an iterative method, then e k+1 ∼ D k,r e r k , where D k,r tends to the asymptotic error constant of the iterative method when k → ∞. To avoid higher order terms in the Taylor series that do not influence the convergence order, we will use the notation used by Traub in [2]. If { f k } and {g k } are null sequences (that is, sequences convergent to zero) and where C is a nonzero constant, we can write Let {x (k) } k≥0 be a sequence of vectors generated by an iterative method converging to a zero of F, with R-order greater than or equal to r, then according to [4] we can write where {D (k,r) } is a sequence which tends to the asymptotic error constant D r of the iterative method when k → ∞, so e k+1 ∼ e k .

Accelerating Parameters
We illustrate the technique of the accelerating parameters, introduced by Traub in [2], starting from a very simple method with a real parameter α This method has first order of convergence, with error equation where e k = x k − ξ, k = 0, 1, . . .
As it is easy to observe, the order of convergence can increase up to 2 if α = 1 f (ξ) . Since ξ is unknown, we estimate f (ξ) by approximating the nonlinear function with a Newton polynomial interpolation of degree 1, at points (x k , f (x k )) and (x k−1 , f (x k−1 )) To construct the polynomial N 1 (t), it is necessary to evaluate the nonlinear function f in two points; so, two initial estimates are required, x 1 and x 0 . The derivative of the nonlinear function is approximated by the derivative of the interpolating polynomial, that is, , and the resulting scheme is substituting α k in the iterative expression we get the secant method [7], with order of convergence p = 1 + √ 5 2 ≈ 1.6180.

Modified Secant Method
A similar process can be done now using a Newton polynomial interpolation of second degree for approximating parameter α. As it is a better approximation of the nonlinear function, it is expected that the order of convergence of the resulting scheme is greater than the obtained in the case of a polynomial of degree 1. For the modified secant method, the approximation of the parameter is α k ≈ 1 being and replacing this expression in an iterative method with memory is obtained Now, three initial points are necessary, x k , x k−1 and x k−2 , to calculate the value of parameter α k . In the following result, we present the order of convergence of the iterative scheme with memory (19). Proof. By substituting the different terms of α k by their Taylor expansions, we have By Writing e k and e k−1 in terms of e k−2 , we obtain where p denotes the order of convergence of method (19). But, if p is the order of (19), then e k+1 ∼ e p 3 k−2 . Therefore, p must satisfy p 3 = p 2 + p + 1. The unique positive root of this cubic polynomial is p ≈ 1.8393 and, by applying the result of Ortega and Rheinbolt, this is the order of convergence of scheme (19).
The efficiency index, defined by Ostrowski in [8], depends on the number of functional evaluations per iteration, d, and on the order of convergence p, in the way so, as in the modified secant method there is only one new functional evaluation per iteration, its efficiency index is 1.8393.

Nonlinear Systems
Method (12) can be extended for approximating the roots of a nonlinear system F(x) = 0, where F : D ⊆ R n −→ R n is a vectorial function defined in a convex set D, It is easy to prove that this method has linear convergence, with error equation being In the multidimensional case, we can approximate the parameter α = [F (ξ)] −1 with a multivariate Newton polynomial interpolation. If we use a polynomial of first degree that corresponds to the secant method for the multidimensional case.
To extend the modified secant method (19) to this context, we need a multidimensional Newton polynomial interpolation of second degree where where B denotes the set of bi-linear mappings. Moreover, evaluating the derivative in x (k) , So, the iterative expression of the resulting method with memory is The divided difference operator can be expressed in its integral form by using the Genocchi-Hermite formula for divided difference of first order [9] [ as and, in general, the divided difference of k-order can be calculated by [10] [x 0 , x 1 , . . . , where

Proof.
For the second order divided difference we write (35) for k = 2 in the following way where h = x k − x k−1 and k = x k−2 − x k−1 . We can write h and k in terms of the errors as h = e k − e k−1 and k = e k−2 − e k−1 . Substituting the approximation of α, α k , in the error equation (24), we obtain e k+1 ∼ c 3 e k−1 e k−2 e k .
Following the same steps as in the unidimensional case, the order of convergence p of method (31) must satisfy the equation p 3 = p 2 + p + 1, which has the unique positive root p = 1.8393.

Dynamical and Numerical Study
When a new method is designed, it is interesting to analyze the dynamical behavior of the rational function that appears when the scheme is applied on polynomials of low degree. For the secant method, this analysis has been made by different researchers, so we are going to study the behavior of the modified secant method. As this scheme would be accurate for quadratic polynomials, due to the approximation we make, we use 3rd degree polynomials in this analysis. In order to make a general study, we use four polynomials that combined can generate any third degree polynomial, that is, p 1 (x) = x 3 − x, p 2 (x) = x 3 + x, p 3 (x) = x 3 and p 4 (x) = x 3 + γx + 1, where γ is a free real parameter.
The rational operators obtained when our method (31) is applied on the mentioned polynomials are: In order to analyze the fixed points, we introduce the auxiliary operators G i , i = 1, 2, 3, 4, defined as follows . We can prove that there are no fixed points different from the roots of the polynomials. So, the method is very stable. On the other hand, critical points are also interesting because a classical result of Fatou and Julia states that each basin of attraction of an attracting fixed point contains at least a critical point ( [11,12]), so it is important to determine to which basin of attraction belongs each critical point. The critical points are determined by the calculation of the determinant of the Jacobian matrix G i , i = 1, 2, 3, 4 ( [13]). For the polynomials p 1 to p 4 the critical points are the roots of the polynomials and the points ( The bifurcation diagram ( [13]) is a dynamical tool that shows the behavior of the sequence of iterates depending on a parameter. In this case, we use parameter γ from p 4 (x). Starting from an initial estimation close to zero, we use a mesh of 500 subintervals for γ in the x axis and we draw the last 100 iterates. In Figure 1a, we plot the real roots of p 4 and in Figure 1b we draw the bifurcation diagram. We can see how the bifurcation diagram always matches the solutions.  Now, we are going to compare Newton's method and the secant and modified secant schemes on the following test functions, under different numerical characteristics, as well as the dynamical planes associated to them. ( Variable precision arithmetics has been used with 100 digits of mantissa, with the stopping criterion |x k+1 − x k | < 10 −25 or | f (x k+1 )| < 10 −25 with a maximum of 100 iterations. That means the iterative method will stop when any of these conditions are met. These tests have been executed on a computer with 16GB of RAM using MATLAB 2014a version. As a numerical approximation of the order of convergence of the method, we use the approximated computational order of convergence (ACOC), defined as [14]: Let us remark that we use the same stopping criterium and calculation of ACOC in the multidimensional case, only replacing the absolute value by a norm. As we only use one initial point to be able to compare methods with memory with Newton's method, to calculate the initial points needed for the methods with memory we use an initial estimation of α. We take α 1 = 0.01 in the unidimensional case. In the multivariate case, we use the initial value α −1 1 = 5I for the secant method, and two different values for the two first iterations, α −1 1 = 5I and α −1 2 = 3I for the modified secant method, where I denotes the identity matrix of size n × n. We have observed that, taking two different approximations of α leads to a more stable behavior of the modified secant method.
On the other hand, for each test function we determine the associated dynamical plane. This tool of the dynamical analysis for methods with memory was introduced in [15]. The dynamical plane is a visual representation of the basins of attraction of an specific problem. It is constructed by defining a mesh of points, each of which is taken as a initial estimation of the iterative method. The complex plane is created showing the real part of the initial estimate in the x axis and the imaginary part on the y axis ( [16]). In a similar way as before, we use an initial estimation of α to calculate the necessary initial points. This approach makes possible to draw the performance of iterative schemes with memory in the complex plane, thus allowing to compare the performance of methods with and without memory. To draw the dynamical planes we have used a mesh of 400 × 400 initial estimations, a maximum of 40 iterations and a tolerance of 10 −3 . We used α = 0.01 to calculate the initial estimations. Each point used as an initial estimate is painted in a certain color, depending on the root to which the method converges. If it does not converge to any root, it is painted black.
In Table 1, we show the results obtained by Newton's method, secant and modified secant schemes for the scalar functions f 1 (x), f 2 (x) and f 3 (x). These tests confirm the theoretical results with a very stable ACOC when the method is convergent.    On the other hand, in Figure 2 we see that the three methods behave well on function f 1 (x), with no significant differences between their basins of attraction. In Figure 3, we see that in case of function f 2 (x) the dynamical plane of the modified secant is better than the secant and very similar to Newton's one. The secant method presents black regions that, in this case, mean slow convergence. In Figure 4, we see the basins of attraction of the methods on f 3 (x). We observe that the wider basins of attraction correspond to the modified secant method. In this case, the black regions mean points where the methods are divergent.    In Table 2 the numerical results corresponding to the vectorial functions F 4 (x 1 , x 2 ), F 5 (x 1 , x 2 ) and F 6 (x 1 , x 2 , x 3 ) are presented. There is nothing significant in these results, they are those expected taking into account the theoretical results. In this case, the definition of ACOC is the same as before, replacing the absolute value by the norm. For systems F 4 (x) = 0 and F 5 (x) = 0, we show in Figures 5 and 6 the basins of attraction of each root. In the first case, the three methods behave similarly, but in the second one, Newton's scheme presents some black regions that do not appear in the modified secant method.

Conclusions
The technique of introducing accelerating parameters allows us to generate new methods with higher convergence order than the original one, with the same number of functional evaluations. The increase is more significant if we start from a method with high order of convergence. As far as we know, this is the first time that an iterative method with memory for solving nonlinear systems is designed by using Newton polynomial interpolation with several variables of second degree. In addition to obtaining the expression of what we called modified secant method, a dynamical study was carried out adapting tools from the dynamics of methods without memory to methods with memory. Although numerically computing the divided differences is hard, methods with memory show a very stable dynamical behavior, even more than other known methods without memory.