A New High-Order Jacobian-Free Iterative Method with Memory for Solving Nonlinear Systems

: We used a Kurchatov-type accelerator to construct an iterative method with memory for solving nonlinear systems, with sixth-order convergence. It was developed from an initial scheme without memory, with order of convergence four. There exist few multidimensional schemes using more than one previous iterate in the very recent literature, mostly with low orders of convergence. The proposed scheme showed its efﬁciency and robustness in several numerical tests, where it was also compared with the existing procedures with high orders of convergence. These numerical tests included large nonlinear systems. In addition, we show that the proposed scheme has very stable qualitative behavior, by means of the analysis of an associated multidimensional, real rational function and also by means of a comparison of its basin of attraction with those of comparison methods.


Introduction
New and efficient iterative techniques are needed for obtaining the solution ξ of a system of nonlinear equations of the form where F : Ω ⊂ R n → R n , Ω ⊆ R n being an open convex set, which is present in scientific, engineering and various other models (details can be found in the articles [1][2][3][4][5]). The search for the solution ξ ∈ R n of the system F(x) = 0 is a much more complex problem than finding a solution of a scalar equation f (x) = 0. As in the scalar case, we can transform the original system into an equivalent system of the form: for a given vector function G : R n → R n , whose coordinate functions we will denote by g i and i = 1, 2, . . . , n. Starting from an initial approximation x (0) = (x x (j+1) = G(x (j) ), j = 0, 1, . . . , n, where x (j) = (x (j) n ) T ∈ R n . We say that the process is convergent if {x (j) } → ξ when j → ∞; then, ξ = (ξ 1 , ξ 2 , . . . , ξ n ) T ∈ R n will be, under certain conditions of the function G, a solution of the system x = G(x). The vector ξ is called a fixed point of the function G and the algorithm described by the Equation (2) is a fixed point method.
A very important concept of iterative methods is their order of convergence, which provides a measure of the speed of convergence. Let x (j) j≥0 be a succession of vectors of R n such that they tend to the solution of the nonlinear system ξ when j tends to infinity. The convergence of this sequence is said to be (i) linear, if there exist M, 0 < M < 1 and j 0 ∈ N such that x (j+1) − ξ ≤ M x (j) − ξ for all j ≥ j 0 . (ii) of order p, if there exist M > 0 and k 0 ∈ N such that x (j+1) − ξ ≤ M x (j) − ξ p for all j 0 .
This definition is independent of the norm of mathbbR n used. We denote by e (j) = x (j) − ξ the error in the j-th iteration and call an error equation to the expression e (j+1) = Le (j) p + O e (j) p+1 , where L is a a p-linear function L ∈ Lle f t(R n , R n , R n right), and p is the order of convergence of the method. The most known fixed-point iterative method is Newton's scheme: where F denotes the Jacobian matrix of operator F. However, there are many practical situations where the calculations of a Jacobian matrix are computationally expensive, and/or it requires a great deal of time for them to be given or calculated. Therefore, derivative-free methods are quite popular for finding the roots of nonlinear equations and systems of nonlinear equations. The first multidimensional derivative-free method was proposed by Samanskii in [6], by replacing the Jacobian matrix F with the divided difference operator: This scheme keeps the quadratic order of convergence of Newton's procedure. It is the vectorial extension of scalar Steffensen's method.
Later on, Traub defined a class of iterative methods (known as Traub-Steffensen's family) [7], given by where u (j) = x (j) + βF(x (j) ). The class (3) can be easily recovered from Newton's wellknown method [7] by replacing the Jacobian matrix with operator [u (j) , x (j) ; F] ≈ F (x (j) ). Let us remark that, for the particular value of β = 1 in expression (3), the deduced scheme is Samanskii's one.
In recent years, different scalar iterative schemes with memory have been designed (a good overview can be found in [8]), mostly derivative-free ones. These have been constructed with increasing orders of convergence, and therefore, with increasing computational complexity. In terms of stability, some researchers compared the amplitude of the set of initial points converging to the same attractor, using complex discrete dynamics techniques. In [9], the authors observed that iterative schemes with seventh-order memory convergence showed better stability properties than many eighth-order optimal procedures without memory. This graphical comparison was subsequently used by different authors; see, for example, the work of Wang et al. in [10] and Cordero et al. [11] in 2016 or the investigations of Bakhtiari et al. [12] in 2016 and Howk et al. [13] in the following years.
Regarding nonlinear vectorial problems, some methods with memory have been developed which improve the convergence rate of Steffensen's method or Steffensen-type methods at the expense of additional evaluations of vector functions, divided difference or changes in the points of iterations. In past and recent years, a few high-order multi-point extensions of Steffensen's method or Steffensen-type methods have been proposed and analyzed in the available literature [14][15][16][17] for solving nonlinear systems of equations. All these modifications are in the direction of increasing the local order of convergence with the view of increasing their efficiency indices, as they usually do not involve new functional evaluations. Therefore, these constructions occasionally possess a better order of convergence and efficiency index, but there are very few iterative schemes of this kind in the literature, due in part to their recent development and also due to the difficulty of design and convergence analysis.
In 2020, Chicharro et al. [18] proposed such an extension, which is given by and its higher-order version is The versions have third and fifth order of convergence, respectively. An extension of this type was first reported by Chicharro et al. [18] based on Kurchatov's divided difference operator.
The authors developed in [19] a technique that, using multidimensional real discrete dynamics tools, is able to analyze the stability of iterative schemes with memory, not only in graphical terms, but essentially in analytical terms. Using this technique, the stability of the fixed and critical points of secant, Steffensen' and Kurchatov's methods (among others) were studied in [19]. It was also used to analyze other procedures, such as those described in [20], the one defined by Choubey et al. in [21] and those by Chicharro et al. in [22][23][24].
The aim of this work was to produce two new schemes without and with memory of orders four and six, respectively. Our scheme is also based on Kurchatov's divided difference operator. However, our scheme does not have only higher-order convergence, unlike the recent scheme (5). However, we did not use any additional functional evaluation of F or F (Jacobian matrix of F) or another iterative substep. We also provide a deep analysis of the suggested scheme regarding the order of convergence (Section 2) and its stability properties, constructing an associated multidimensional discrete dynamical system. Therefore, the good performance in terms of convergence to the searched roots and wideness of their basins of attraction is proven, as on polynomial functions as on non-polynomial ones, in Section 3. In addition, we compare in Section 4 our methods to the existing recent methods on several numerical problems with similar iterative structures. On the basis of the results, we found that our methods perform better than the existing ones when dealing with residual errors, the difference between two consecutive iterations and stable computational order of convergence. Finally, some conclusions and the references used bring this manuscript to an end.

Construction and Convergence of New Iterative Schemes
Combining the Traub-Steffensen family of methods and a second step with different divided-difference operators, we propose the class of iterative schemes described as where u (j) = x (j) + βF(x (j) ), β ∈ R. In order to analyze the convergence of schemes (6), we need the definition of the divided difference operator as well as its Taylor's expansion (more details can be found in [5]).

Lemma 1. Suppose F
: Ω ⊂ R n → R n is k-times Fréchet differentiable in an open convex set Ω.
In addition, we assume that for any x, h, ∈ R j , the following expression holds: where Now, we can obtain the following Taylor series expansion of the divided difference operator, by adopting the Genocchi-Hermite formula [5]: Then, we have Now, we are in a position to analyze the convergence order of proposed schemes (6), as we can see in the following result. In it, I denotes the identity matrix of size n × n and F (k) (x) can be considered as a k-linear operator: Todeepen one's understanding of the concepts of Taylor expansion using several variables, we suggest references [5,25].

be a sufficiently differentiable function in an open convex
neighborhood Ω of a zero ξ. Suppose that F (x) is a continuous and nonsingular at x = ξ and the initial guess x (0) is close enough to ξ. Then, the iterative schemes defined by (6) have a fourth order of convergence for every β, β = 0. They satisfy the following error equation:

Proof.
Let e (j) = x (j) − ξ be the error of the jth-iteration and ξ ∈ R n be a solution of F(x). Then, developing F(x (j) ) in the neighborhood of ξ, we have and By adopting the expressions (8)-(12), we obtain Inversion of u (j) , x (j) ; F is given by where Using expressions (9) and (14), we get Now, use the expression (15) in the first substep of our proposed scheme (6). We have In a similar fashion as we did in expression (13), we can develop F(y (j) ) and its derivatives in the neighborhood of ξ, given by By using the property u (j) , and similarly, we have By using the expressions (14) and (17)- (19) in the second substep of our scheme (6), the error equation of our scheme (6) becomes: Therefore, (6) is a parametric family of fourth-order iterative methods.

Extension to a Higher-Order Scheme with Memory
In this section, we construct a new scheme with memory based on our scheme (6), without using any new functional evaluation. It is straightforward to say from error Equation (20) that we can obtain a higher order of convergence by choosing β = −F (ξ) −1 . However, we also know that the required solution ξ is unknown. We want to increase the order of convergence without additional values of vector function or Jacobian matrix. Therefore, we use one of the most efficient Kurchatov's divided difference operators, ; F , to approximate the value of F (ξ), which is given by Hence, from our scheme (6) we can deduce the following iterative method with memory: where ). In the following result, we analyze the convergence order of scheme (22) with memory, denoted by PM 6 . Theorem 2. Let F : Ω ⊆ R n → R n be a sufficiently differentiable function in an open convex neighborhood Ω of a zero of F, ξ. Suppose that F (x) is continuous and non-singular in ξ. In addition, the initial guesses x (0) and x (1) are close enough to the required solution ξ. Then, the iterative scheme defined by (22) has sixth order of convergence.
Proof. In a similar way as we did in expressions (13) and (17), we expand where O e (j−1) , e (j) 3 stands for the sum of the exponents of e (j−1) and e (j) is at minimum 3. We obtain the inverse of 2x (j) − x (j−1) , x (j−1) ; F is given by which further yields Therefore, it is clear that we have I + A (j) F (ξ) ∼ e (j) and adopt this value in the expression (20). Then, we have Hence, the proof is complete, which demonstrates that the proposed scheme with memory (22) has sixth-order convergence.

A Qualitative Study of Iterative Methods with Memory: New and Known
In this section, we analyze the stability of the proposed scheme with memory in its scalar form, as this kind of study yields multidimensional operators to be analyzed and it does not support the use of vectorial schemes. However, its performance on systems of nonlinear equations is checked in Section 5 and it is compared with other known schemes with memory.
The expression of a scalar fixed-point iterative method with memory, using two previous iterations to calculate the following estimation, is x 0 and x 1 being the starting estimations. We use the technique presented in [19,20] to describe any method with memory as a discrete real vectorial dynamical system, in order to analyze its qualitative behavior.
In order to calculate the fixed points of an iterative method with iteration function Φ, an auxiliary fixed point multidimensional function Γ : R 2 −→ R 2 can be defined, related to Φ by means of: x 0 and x 1 being, again, the initial estimations. Therefore, a fixed point of this operator will be obtained when not only x j+1 = x j , but also x j−1 = x j .
From function Γ : R 2 → R 2 , a discrete dynamical system in R 2 can be defined by Γ x j−1 , x j = (x j , x j+1 ), where Φ is the operator of the iterative method with memory. The fixed points (z, x) of Γ satisfy z = x and x = Φ(z, x). This notation implies z = x j−1 and x = x j . In the following, we recall some basic dynamical concepts that are direct extensions of those used in complex discrete dynamics analysis (see [26]).
Let us consider the vectorial rational function Γ : R 2 → R 2 , usually obtained by applying an iterative method on a scalar polynomial p(x). Then, if a fixed point (z, x) of operator Γ is different from (r, r), where r is a zero of p(x), it is called a strange fixed point. On the other hand, the orbit of a pointx ∈ R 2 is defined as the set of successive images fromx by the vector function-that is, orbit(x) = {x, Γ(x), . . . , Γ n (x), . . .}. Indeed, a point x * ∈ R 2 is called k-periodic if Γ k (x * ) = x * and Γ p (x * ) = x * , for p = 1, 2, . . . , k − 1.
The qualitative performance of a point of R 2 is classified depending on its asymptotic behavior. Thus, the stability of fixed points for vectorial operators satisfies the statements appearing in the following result (see, for instance, [27]).

Theorem 3.
Let Γ from R m to R m be of class C 2 . Assume x * is a k-periodic point. Let λ 1 , λ 2 , . . . , λ m be the eigenvalues of the Jacobian matrix Γ (x * ). Then, it holds that (a) If all the eigenvalues λ j verify |λ j | < 1, then x * is attracting. (b) If one eigenvalue λ j 0 verifies |λ j 0 | > 1, then x * is unstable-that is, repelling or saddle. (c) If all the eigenvalues λ j verify |λ j | > 1, then x * is repelling.
Moreover, a fixed point is said to be not hyperbolic if all the eigenvalues λ j of Γ (x * ) satisfy |λ j | = 1. Specifically, if there exist an eigenvalue λ i satisfying |λ i | < 1 and another one λ j such that |λ j | > 1, then it is called saddle point.
There is a key difference between the study of the stability of a fixed point x * in scalar and vectorial dynamics: In the first case, if |Γ (x * )| < 1 then x * is attracting; in particular it is superattracting if |Γ (x * )| = 0 and it is repelling when |Γ (x * )| > 1, Γ being the scalar rational function related to the iterative scheme on a low-degee polynomial p(x)). In the vectorial case, the character of the fixed points is calculated by means of the eigenvalues of the Jacobian matrix Γ (see Theorem 3). Nevertheless, sometimes the Jacobian is not well-defined at the fixed points. In these cases, we impose on the rational operator G the condition w = z = x, so that it is reduced to a real-valued function. Therefore, the stability of the fixed point can be inferred from the absolute value of its first derivative at the fixed point.
By considering x * an attracting fixed point of function Γ, we define its basin of attraction A(x * ) as the set of preimages of any order: A(x * ) = x 0 ∈ R 2 : Γ m (x 0 ) = x * , for some m ∈ N .
A key element in the stability analysis of an iterative method is the set of critical points of its associated rational function Γ: if Γ (x) satisfies det (Γ (x)) = 0, x is called a critical point. A critical point (c, c) such that c is not a root of p(x) is called a free critical point. Another way to get critical points is finding those points that make null the eigenvalues of Γ . As an extension of the scalar case, if they are ant composed by the roots of polynomial p(x), they are named free critical points. Indeed, Julia and Fatou [26] proved that there is at least one critical point associated with each basin of attraction. Therefore, by studying the orbit of the free critical points, all the attracting elements can be found. This result is valid for both complex and real iterative functions.
These schemes have several similarities: they are vectorial schemes with memory; they include divided differences operators in their iterative expressions; they are mainly used to define their respective accelerating parameters. However, we are going to see that the use of these elements does not determine the wideness of the sets of initial estimations converging to the roots, when real multidimensional discrete dynamics tools are used.
In order to extend the results to any quadratic polynomial, the first analysis is shown for PM 6 on p(x) = x 2 − c, so that the value of c yields a situation with real, complex or multiple roots depending on c > 0, c < 0 or c = 0, respectively. The multidimensional rational function Γ in this particular case will be denoted in what follows by PM. This analysis can be summarized in the following result.
and it is Proof. Let us remark that operator PM(z, x) can be obtained by directly applying method PM 6 to polynomial p(x). Moreover, we know that fixed points of PM will have equal components. This is the reason why, when we force the three consecutive iterates to be equal (x = z), then the only fixed points are composed by the roots . Regarding the critical points, the Jacobian matrix PM is By definition, the components of critical points are those values making null the eigenvalues of PM . By using the change of variables t = x 2 on the second factor of the numerator of the not null eigenvalue, we get s(t) = c 4 + 32c 3 t + 210c 2 t 2 + 360ct 3 + 165t 4 . This polynomial only has real roots for c < 0, denoted by s i , i = 1, 2, 3, 4. Then, if c < 0, there exist eight different componts of free critical points (± A very useful tool to visualize the analytical results is the dynamical plane of the system, composed by the set of the different basins of attraction. It can be drawn by means of the programs presented in [28], after some changes to adapt them to schemes with memory. The dynamical plane of a method is built by calculating the orbit of a mesh of 400 × 400 starting points (z, x) (although z does not appear in the rational function PM) and painting each of them in a different color (orange and green in this case) depending on the attractor they converge to (marked as a white star), with a tolerance of 10 −3 . Additionally, they appear in black if the orbit has not reached any attracting fixed point in a maximum of 80 iterations. In Figure 1, we show the dynamical planes of this method for selected values of c, in order to show its performance. Let us remark that, as by definition all the fixed points have equal components, they will always appear in the bisector of the first and third quadrants of the dynamical plane. It can be observed that, when there are no real roots (c < 0, Figure 1a), no other attracting element appears; when c = 0, the only root is multiple and the convergence is linear, so there is global convergence to x = 0, as can be seen in Figure 1b. In Figure 1c, the convergence to the roots is also observed to be global, their basins of attraction being two symmetrical half-planes, which is exactly the same behavior as Newton's method on quadratic polynomials.
Moreover, let us remark that when c > 0 (real simple roots case), there are not free critical points, as in this case the only possible performance of the method was the convergence to the roots. The reason for this behavior is that in each basin of attraction there must be a critical point; if the only critical points are the roots of that basin, then there is no other possible convergence.

Comparisons with Other Methods with Memory for Nonlinear Problems
Here we compare the stability of the proposed method with that of known AM 5 and SM 4.45 schemes. Firstly, we show their performances on quadratic polynomials: PM 6 . The dynamical planes are plotted in the complex plane, by starting the iterative methods with memory with an initial value of the accelerating parameter of −0.01 om each case, for any initial guess z ∈ C, defined in a mesh of 400 × 400 points and with a maximum of 80 iterations.
As shown in Figures 2 and 3, all three methods have been used to estimate the complex roots of the unity of second and third degrees. It can be observed that the performances with PM 6 and AM 5 were very similar for quadratic polynomials, showing global convergence, similar to that of Newton's scheme.  However, this global convergence is held on cubic polynomials in the case of PM 6 , but in case of AM 5 , black areas of no convergence to the roots appear. Iterative method SM 4.45 shows these black regions both for quadratic and cubic polynomials, which are bigger in the cubic case.
This good performance of the set of converging initial estimations of proposed method PM 6 is also shown in the case of non-polynomial equations: let us notice the performances of all known and new methods on the complex rational function f (z) = z 2 − 1 with a simple zero at ξ = 0. The basins of attraction of the root appearing in Figure 4, being similar, are wider in case of PM 6 than for AM 5 or SM 4.45 .
Additionally, in Figure 5 the basin of attraction of ξ ≈ 0.25753, the simple root of g(z) = z 2 − exp(z) − 3z + 2, is very wide in the case of PM 6 , with small areas of no convergence to the root, in comparison with those of known methods AM 5 and SM 4.45 .
Thus, it can be concluded that proposed method PM 6 had very stable performance, as in real and in complex spaces, both for polynomial and non-polynomial functions, in the scalar case. In spite of having the same accelerating parameter as the AM 5 method, the new scheme has proven to be better than known ones in terms of order of convergence and also in terms of stability. In the next section, its numerical performance on nonlinear systems of different sizes is demonstrated.

Numerical Experiments
This section presents the validity of the proposed scheme with memory on some numerical problems. The proposed method (22) was applied to numerical problems and compared with other existing techniques with memory [15,18], presented as AM 5 and SM 4.45 , respectively. In all numerical problems, initially the parameter β = B (0) = −0.01I, where I is the identity matrix, and for method SM 4.45 , another parameter c = −0.01 was considered. All the numerical tests were conducted by using Mathematica 10, with 400 multiple precision digits of mantissa. For all the examples, we have included in the respective tables the following information: x (j+1) − x (j) ; j = 0, 1, 2; the residual F(x (j) ) at j = 0, 1, 2; approximated computational order of convergence (ACOC) denoted as ρ [29] , for each j = 2, 3, . . . ; (27) and CPU − time.
Further on, the iterative procedure was stopped after three iterations and problems were tested on three different initial values. Notice that the meaning of b(±a) is b × 10 ±a in all the tables. Example 1. Consider the following system of nonlinear equations in four unknown variables: The approximate solution is (0.57735, 0.57735, 0.57735, −0.28867) T . Table 1 depicts that the proposed scheme with all different initial guesses converged to a solution much faster than the methods AM 5 and SM 4.45 . Clearly, the residual error, functional error and computational time of proposed method PM 6 were superior to those of AM 5 and SM 4.45 . Table 1. Convergence behavior of the schemes for Example 1.
This problem has an approximated solution of (0.0736323 . . . , 0.0736323 . . . , . . . , 0.0736323 . . . ) T . The performance of the proposed method was better than those of the other methods ones, as shown in Table 3. Table 3. Convergence behavior of schemes for Example 3. Example 4. We have tested the proposed method on another well-known nonlinear problem Fisher's equation [30] which has many applications in chemistry, heat and mass transfer, biology and ecology. Basically, it determines the process of interaction between diffusion and reaction. This nonlinear problem with homogeneous Neumann's boundary conditions and diffusion coefficient D can be defined as: u t = Du xx + u(1 − u) = 0, u(x, 0) = 1.5 + 0.5cos(πx), 0 ≤ x ≤ 1, u x (0, t) = 0, ∀t ≥ 0, u x (1, t) = 0, ∀t ≥ 0. (31) Applying the finite difference discreitization on the Equation (31) leads to a system of nonlinear equations. Suppose w i,j = u(x i , t j ) is its approximate solution at the grid points of the mesh. Let M and N be the numbers of steps in x and t directions, and h and k be the respective step sizes. Use central difference to approximate the second order partial derivative u xx (x i , t j ) = (w i+1,j − 2w i,j + w i−1,j )/h 2 ; backward difference approximation for the first-order derivative with respect to "t" as u t (x i , t j ) = (w i,j − w i,j−1 )/k; and forward difference for first order derivative with respect to x as u x (x i , t j ) = (w i+1,j − w i,j )/(h). The solution of the system is obtained by taking steps along x-axis, M = 21, and t-axis, N = 21, which form a nonlinear system of size 400, with the initial vector x (0) = (i/(M − 1) 2 ) T , i = 1, 2, . . . , M − 1. The results have been computed by different methods and are shown in Table 4. It can be noticed that the lowest execution time and residual error at the third iteration correspond to the proposed method PM 6 . Moreover, the approximate solution has been plotted in Figure 6.  Ap pr ox im ate so l

Conclusions
The number of iterative schemes with memory for solving multidimensional problems is low, in part due to the difficulty of the task. Additionally, it is in part due to the lack of efficiency of the resulting method, if the usual techniques employed in the design of scalar methods with memory are employed as high-degree interpolation polynomials. With the procedure used, the iterative expression of the scheme with memory remains simple, and the order of the original method is increased by 50%. Thus, the efficiency is highly improved. Moreover it has been proven, by means of the associated multidimensional discrete dynamical system, that it is a very stable scheme with wide basins of attraction and global convergence on quadratic polynomials. Its performance on other nonlinear functions was also found to be very stable in comparison with other known schemes. The numerical tests confirmed these results, even for large nonlinear systems and applied problems, such as Fisher's partial differential equation.