Stability Anomalies of Some Jacobian-Free Iterative Methods of High Order of Convergence

: In this manuscript, we design two classes of parametric iterative schemes to solve nonlinear problems that do not need to evaluate Jacobian matrices and need to solve three linear systems per iteration with the same divided difference operator as the coefﬁcient matrix. The stability performance of the classes is analyzed on a quadratic polynomial system, and it is shown that for many values of the parameter, only convergence to the roots of the problem exists. Finally, we check the performance of these methods on some test problems to conﬁrm the theoretical results.


Introduction
Systems of nonlinear equations must usually be solved when nonlinear models, appearing in Science and Engineering, are discretized. There are no analytical techniques for solving these systems, so we approach their solutions by using iterative schemes. Although the most known iterative procedure is Newton's scheme, in recent years, the focus of this area of research has been in constructing new iterative methods, trying to improve Newton's one, in terms of convergence, efficiency, and stability (see, for example, some third-order schemes in References [1][2][3][4][5][6], or higher-order ones in References [7][8][9][10][11][12]). The key fact to get the most efficient methods is to evaluate as few Jacobian matrices as possible, per iteration (see Reference [13] and the references therein).
In 2016, the authors of [14] proposed a parametric class of iterative schemes with fourth-order of convergence, including a very efficient fifth-order procedure. The efficiency index used is defined as p 1/d+op , where p is the order of convergence of the method, d the number of functional evaluations per iteration, and op the number of products/quotients used for obtaining a new iterate. This class combined three evaluations of the nonlinear function and only one of the Jacobian matrix, per iteration. In addition to the study of the local convergence and the computational cost of an iterative method, we must analyze its dependence on the set of initial estimations used. This analysis was made in Reference [15].
In this manuscript, we modified this class of iterative methods to avoid the calculation of Jacobian matrices, getting two different families, whose difference is based on the order of the estimation of the Jacobian matrix. Then, using the real multidimensional dynamics, we studied the behavior of the elements of the family (usually with the same order of convergence) and decided which ones are the best, from this point of view. Specifically, we studied the multipliers of the strange fixed points, analyzed the asymptotical behavior of the free critical points, and found different kinds of strange attractors using bifurcation diagrams. At the end, the dynamical planes are shown to check the theoretical results.

Introductory Concepts
Now, let us recall some concepts that we will use along this manuscript. Let F(x) = 0 be a system of n equations with n variables, where F : D ⊆ R n → R n , with coordinate functions f i , i = 1, 2, . . . , n. Our aim is to find a solutionx ∈ D as a fixed point of the iteration method: x (k+1) =Ḡ(x (k) ), k = 0, 1, . . . , (1) x (0) being the initial guess. Let us denote with G(x) the vectorial fixed-point rational function associated to the iterative method applied to n-variable polynomial system p(x) = 0, p : R n → R n , x ∈ R n p(x). Most of the following concepts are direct extensions of those considered in complex dynamics (see, for example, [16,17]). The orbit of x (0) ∈ R n is defined as and its stability is characterized in a result from Robinson ([16], p. 558), stating that the character of a period-k point x * depends on the eigenvalues of G (x * ), λ 1 , λ 1 , . . . , λ n . It is attracting if all |λ j | < 1, repelling if all |λ j | > 1 (j = 1, 2, . . . , n)), and unstable or saddle if at least one λ j 0 exists such that |λ j 0 | > 1. In addition, a fixed point is called hyperbolic if all the eigenvalues λ j of G (x * ) have |λ j | = 1.
On the other hand, if an iterative scheme has at least order two, then the zeroes of the function are superattracting fixed points of G. We call a strange fixed point to that not being a zero of p(x). The basin of attraction A(x * ) of an attracting fixed point . . , n}, then x is a critical point. If it is not a root of polynomial p(x), it is called a free critical point.
In this manuscript, we apply these techniques (as it has been previously made, for example, in Reference [17]), to a couple of Jacobian-free families of iterative schemes. Specifically, we study the strange fixed points and their stability, prove the existence of free critical points, and find different chaotic behaviors (strange attractors, period-doubling cascades, etc.). In the last section, some basins of attraction are shown to visualize the performance of the elements of the family. Some numerical results are presented.

Jacobian-Free Parametric Families
We propose the following Jacobian-free version of a fourth-order iterative class designed by the authors of [14], denoted by M41, whose iterative expression is: Because of the technique used in the proof, these methods have an order of convergence four in a neighborhood ofx, where F (x) is nonsingular. The following result shows its local convergence, using the notation about the Taylor expansion of F(x (k) ) and F (x (k) ) aroundx, introduced in Reference [18]. To get the Taylor expansion of the divided difference operator [·, ·; F] : D × D ⊂ R n × R n → L(R n ) (see ), its integral expression is used (known as the Genochi-Hermite formula): Theorem 1. Letx ∈ D be a zero of F : D ⊂ R n → R n , a sufficiently differentiable Frechét function in an open convex set D. Let us assume that F (x) is nonsingular and let x (0) be an initial estimation near tox. Class (2) has fourth-order of convergence for any β ∈ R, β = 0. The error equation of the family is: ., e k = x (k) −x, k ≥ 0, and I denotes the n × n identity matrix.
Proof. The Taylor development of F(x (k) ) aroundx is: being also the expansion of its successive derivatives: Then, denoting by w (k) = x (k) + F(x (k) ) and e w = w (k) −x: Then, the expansion of error at the first step of the iterative process e y = y (k) −x is: and, therefore: F(y (k) ) = F (x) e y + C 2 e 2 y + O(e 3 y ).
Finally, the error equation of the method can be expressed as: that, expressed in terms of powers of e k and after some algebraic simplifications, results in: proving the fourth-order of convergence of all the elements of the family for any real β = 0. Now, if a second-order estimation of the Jacobian matrix is made in the original M4 family of Reference [14] in a similar way to that in Reference [20], then it can be expressed as: ; F , which we denote as M42. The following result shows its local order of convergence.
Theorem 2. Letx ∈ D be a zero of F : D ⊂ R n → R n , a sufficiently differentiable Frechét function in an open convex set D. Let us assume that F (x) is nonsingular, and let x (0) be an initial estimation near tox. If β = 0, the family defined by (3) has, at least, fourth-order of convergence for any β ∈ R, being fifth-order convergent when β = 5. The error equation of the class is: and in a similar way as in the proof of Theorem 1, the Taylor expression of the second-order divided difference and its inverse can be expressed as: Then, the error at the first step is: and again: Thus, the expansion of the error at the second step is: Therefore: Finally, the error equation of the method can be expressed as: that, expressed in terms of e k and simplifying: proving that all the members of the class have, at least, fourth-order of convergence, except the element corresponding to β = 5, which is fifth-order convergent.
In the following sections, we use the properties of the discrete dynamical systems associated to these families in order to select their most stable members.

Multidimensional Dynamical Analysis
In the following, we denote by . . , M i n (x, β)), i = 1, 2 the fixed point function associated to M41 and M42 classes, respectively, applied on n-dimensional quadratic polynomial p(x) = 0, where:

Class M41
As the polynomial system has separated variables, all the components of M i (x, β) (for a fixed i ∈ {1, 2}) have the same expression, with the only difference of the sub-index corresponding to the component of x. In the case of method M41, the coordinate functions of the multidimensional rational operator can be expressed as: for j = 1, 2, . . . , n.
All the information about the stability analysis of these fixed points appears in the next result.

Family M42
When the rational operator associated to method M42 on p(x) is analyzed, its coordinate functions can be expressed as: To calculate the fixed points of the multidimensional rational function, we get the roots of M 2 (x, β) = x. We obtain again x j = ±1 and also the real roots of the polynomial r 2 (t) = (β + 88)t 6 + (−3β − 32)t 4 + (3β + 8)t 2 − β. It can be checked that r 2 (t) only has, at most, two real roots, which are denoted by r 1 2 (β) and r 2 2 (β). However, the rational function M 2 (x, β) is the same as that analyzed in Reference [15], corresponding to the original fourth-order iterative class designed by the authors of [14], who use Jacobian matrices in its iterative expression applied on p(x).
Although the classes of iterative methods are different (that of Reference [14] uses Jacobian matrices and the proposed one is Jacobian-free), the second order estimation of the Jacobian matrix has the same effect on quadratic polynomial systems. So, the following result coincides with that of [14]. Theorem 4. The 2 n roots of p(x) are superattracting fixed points of the rational function associated to M42. This operator also has a different number of real strange fixed points whose values depend on β: If β > 0 or β < −88, r 1 2 (β) and r 2 2 (β) are real, being their respective eigenvalues of the associate Jacobian For −88 ≤ β ≤ 0, the only real fixed points are the roots z = ±1 of p(z); thus, there are no strange fixed points.
In the next section, we delve into the bifurcation analysis of these free critical points in order to locate other undesirable behaviors, such as attracting periodic orbits, strange attractors, etc.

Critical Points and Bifurcation Diagrams
Firstly, we analyze the case of class M41, that is, the rational function M 1 (x, β) and its critical points. We use Feigenbaum diagrams to analyze the bifurcations of the map related to each family on system p(x) using each one of the free critical points of the map as a starting point and observing their behavior for different ranges of β. When the rational function is iterated on these critical points, different behavior can be found after 1000 iterations of the method corresponding to each value of β in a mesh of 3000 subintervals. The resulting behavior is from convergence to the roots to periodic orbits or even other attractors. In the case of class M41, we observe in Figure 2 the orbits of these critical points, in the range of β where each one of them is real. A clear convergent behavior to the roots in the intervals is observed where the critical points are real, except inside the interval [60, 70], where a chaotic region can be seen where a strange fixed point has bifurcated (for β = 62.0613, there is change in the number of strange fixed points, two of them parabolic, see Theorem 3) into period-doubling cascades. Further, some blue regions appear associated with chaotic behavior. In them, strange attractors can be found (in the bidimensional case, n = 2). To do it, we plotted in the (x 1 , x 2 )-space the iteration of M (x 1 , x 2 ), for values of β in the blue area. For each fixed value of parameter β, 1000 different initial estimations have been used and, for each of them, the iterates have been plotted following this code color: The first 100 iterations are not plotted, while the following 400 are plotted in blue color, the next 400 in green color, and the last hundred in magenta color. This can be observed in Figure 3a,b, as the strange fixed point bifurcates in periodic orbits of increasing period, until the chaotic behavior appears (Figure 3c,d, where the orbits are dense in a small rectangular region of the x 1 , x 2 space).  Figure 4a,b, a clear symmetry in the convergence to the roots is observed, and we can also see nonconvergent behavior to the roots for elements corresponding to β ∈ [35, 40]. We see in Figure 4c four period-doubling bifurcation areas where the repulsive fixed points have bifurcated into period-doubling orbits.
In order to better visualize this behavior, we plotted in the (x 1 , x 2 )-plane the iteration of M 2 (x, β), for β in the blue area of Figure 5a near to β = 36. Then, symmetric strange attractors appear (see Figure 5b,c).
As in this case, the four strange attractors must be separately identified, the way these pictures have been obtained is slightly different than that used in case M41 family: Fixing the value of parameter β, 1000 different initial estimations have been taken in small rectangles close to the origin. For each specific value of β, the method M 2 (x, β) has been used on each of them, plotting one point per iteration. The code color used is as follows: The first 800 iterations do not appear, while the following 100 appear in blue color and the last hundred are shown in magenta color. The resulting plots show how the four attracting fixed points change into attracting areas, being disjoint or not depending on β. Nevertheless, the set of starting guesses belonging to their respective basins of attraction is very small, as well as the interval of real values of β inducing this performance.

Dynamical Planes for n = 2
In this section, we compare the sensitivity of the proposed classes to the initial estimation, depending on some values of the parameter β that have resulted to provide stable or unstable behavior, for each one of the families.
In Figures 6 and 7, we plotted the dynamical planes corresponding to M41 for several values of β. They were obtained using the routines appearing in Reference [21]. A mesh of 400 × 400 points was used, the maximum number of iterations employed was 40, and the stopping criterium has a tolerance of 10 −3 . We plotted a point of this mesh of different colors depending on the root they converge to. Moreover, the color is brighter when the number of iterations used is lower. It is plotted in black color if it reaches the maximum number of iterations without converging to any of the roots.
We observe in Figure 6 that the four roots of the vectorial polynomial have their respective wide basins of attraction (colored in orange, cyan, blue, and purple), with several connected components for each root, separated by the Julia set. There are also four repulsive fixed points and 16 saddle points. In case β = 70, there exist four superattracting fixed points, sixteen repulsive, and other sixteen saddle fixed points. It can be observed in Figure 6b that the number of connected components in the plotted area is much higher, and the immediate basin of attraction (the component holding the fixed point) is smaller in general. On the other hand, unstable behavior appears in Figure 7, corresponding to the values β = −130 and β = −0.5, where nine strange attracting fixed points had been found, in both cases. They are located in the colored areas of no convergence to the roots, presented in yellow, brown, gray, pink, and green, among others. There are also 25 repulsive and 30 saddle fixed points for β = −130 and, in the case of β = −0.5, there are 9 repulsive and 18 saddle fixed points. Regarding the symbols, fixed points are marked with a white circle, while those attracting are shown with a white star. Let us remark that repulsive and nonhyperbolic points are always in the Julia set, and attracting ones lay in their respective basin of attraction.
To sum up, the stability of the members of Jacobian-free class M41 of iterative schemes has been analyzed for a quadratic multivariate polynomial. Unstable behavior, in terms of attracting strange fixed points, has been located for β ≤ −0.4991, the rest of elements being stable, with wide sets of convergent initial estimations.

Class M42
The dynamical planes related to the M42 class on p(x) for different values of parameter β are shown in Figures 8 and 9.  We see in Figure 8a, for β = −10, that the basins of attraction of the four zeroes of the polynomial form a balanced division of the real plane into four parts (with only a connected component for each root). The fifth-order case β = 5 (Figure 8b) is also stable, but the number of basins of attraction increases slightly with the value of the parameter, and they are divided in several connected components.
Finally, the unstable behavior is shown in Figure 9a,b, for β = 36.15 and β = 36.175, respectively, where strange attractors appear. They are located in the black regions, showing no convergence to the zeroes.
To summarize, the stability of the elements of classes M41 and M42 of iterative methods were studied for quadratic multivariate polynomial. The only unstable behavior was located in small intervals where attracting periodic orbits or dense attracting regions were found, and the rest of the elements were stable, the most stable ones being clearly stated. This is numerically checked in the following section.

Numerical Performance
The numerical results shown in this section were obtained using Matlab R2014b, with variable precision arithmetics with 200 digits of mantissa. In each table of results, the residual F(x (k) and the difference between consecutive iterations x (k) − x (k−1) are presented for the first three iterations. The approximated computational order of convergence (ACOC) is also calculated, by means of (see [5]): The nonlinear systems used in these tests, F(x 1 , x 2 , . . . , x n ) = 0, are defined by the following coordinate functions, joint with the corresponding size and the estimated solution to be found: For each Jacobian-free class, M41 and M42, two stable elements were used (β = 5 and β = 10 for M41 and β = −10 and β = 5 for M42), and also two respective unstable elements (β = 68 and β = −130 for M41 and β = 36 and β = 36.15 for M42), expecting to see their performance on nonpolynomial functions.
In both examples of nonpolynomial systems, it is observed that, being good results in general, those obtained using values of parameter β that have shown to be stable in the dynamical analysis give more precise results than those that are considered unstable, as can be checked in Tables 1 and 2. Class   Class

Conclusions
Two highly efficient Jacobian-free classes of iterative methods were designed, and it has been proven that their order of convergence is four, except for a fifth-order particular case at β = 5 for the second class. The real multidimensional stability of the members of both classes of iterative schemes was analyzed for quadratic n × n polynomial systems. For the first family, only two real components of strange fixed points that can yield attracting points different from the roots appear. In the second one, there are no attracting strange fixed points. On the other hand, free real critical points fall in chaotic behavior only in small intervals for both cases. The rest of the members are stable, being the most similar performance to Newton's scheme located around β = 5. Finally, numerical experiments confirm the theoretical results, showing good performance in all cases but having lower residual errors in those cases that correspond to the stable values of the parameter.