Chaos and Stability in a New Iterative Family for Solving Nonlinear Equations

: In this paper, we present a new parametric family of three-step iterative for solving nonlinear equations. First, we design a fourth-order triparametric family that, by holding only one of its parameters, we get to accelerate its convergence and ﬁnally obtain a sixth-order uniparametric family. With this last family, we study its convergence, its complex dynamics (stability), and its numerical behavior. The parameter spaces and dynamical planes are presented showing the complexity of the family. From the parameter spaces, we have been able to determine different members of the family that have bad convergence properties, as attracting periodic orbits and attracting strange ﬁxed points appear in their dynamical planes. Moreover, this same study has allowed us to detect family members with especially stable behavior and suitable for solving practical problems. Several numerical tests are performed to illustrate the efﬁciency and stability of the presented family.


Introduction
Many problems in Computational Sciences and other disciplines can be stated in the form of a nonlinear equation or nonlinear systems using mathematical modeling. In particular, a large number of problems in Applied Mathematics and Engineering are solved by finding the solutions of these equations.
In the literature, there are many methods and families of iterative schemes that have been designed by using different procedures to approximate the simple roots of a nonlinear equation f (x) = 0, where f : I ⊆ R → R is a real function defined in an open interval I. We can find in [1][2][3] several surveys and overviews of the iterative schemes published in the last years. Each method has a different behavior. This behavior is characterized with the efficiency criteria and the complex dynamics tools.
In this paper, we introduce a new family of multistep iterative schemes to solve nonlinear equations, which contains as an element of this family, a particular method presented in [4]. This family is built from the Ostrowski's scheme, adding a Newton step with a "frozen" derivative and using a divided difference operator. Therefore, the family has a three-step iterative expression. Furthermore, it has three arbitrary parameters named α, β, and γ, which can take real or complex values, and an order of convergence of at least four. The order of convergence will be discussed in Section 2.
From the error equation we observe by fixing two parameters in function of the third one, an uniparametric family of sixth-order iterative methods is obtained. We analyze the dynamical behavior of this family in terms of values of the parameter, in order to detect its elements with good stability properties and others with chaotic behavior. The concept of chaos has been widely discussed (see, for example, in [5]) and it is commonly understood as the presence of complex orbit structure and extreme sensitivity of orbits to small perturbations. Moreover, the presence of unstable periodic orbits of all periods is also included in the concept of chaotic system. For this study, we use tools of discrete complex dynamics that we introduce in Section 3.
In Section 4, we present the performance of the presented schemes on several test functions. These numerical tests allow us to confirm the results obtained in the dynamical section and to compare our schemes with other known ones. The manuscript finishes with some conclusions and the references used in it.

Convergence of the New Family
In this section, we perform the convergence analysis of the new triparametric iterative family. Furthermore, we propose a strategy to reduce the triparametric scheme to an uniparametric scheme in order to accelerate the convergence.
Theorem 1. Let f : I ⊆ R → R be a sufficiently differentiable function on an open interval I and ξ ∈ I a simple root of the nonlinear equation f (x) = 0. Suppose that f (x) is continuous and sufficiently differentiable in an environment of the simple root ξ, and x 0 is an initial estimate close enough to ξ. Then, the sequence {x k } k≥0 obtained by using the expression (1) converges to ξ with an order of convergence of four, being its error equation and q = 2, 3, ...
Using Taylor expansion of f (x k ) and f (x k ) around ξ, we have and where (4), we get
From Theorem 1, it follows that the new triparametric family of iterative methods has an order of convergence of four for any real or complex values of the parameters α, β, and γ. However, convergence can be speed-up if only one parameter is held and the family is reduced to an uniparametric iterative scheme. The latter can be seen in Theorem 2.
Theorem 2. Let f : I ⊆ R → R be a sufficiently differentiable function on an open interval I and ξ ∈ I a simple root of the nonlinear equation f (x) = 0. Suppose that f (x) is continuous and sufficiently differentiable in an environment of the simple root ξ, and x 0 is an initial estimate close enough to ξ. Then, the sequence {x k } k≥0 obtained by using the expression (1) converges to ξ with an order of convergence of six, provided that β = 1 + α and γ = 1 − α, being its error equation and q = 2, 3, ...
Using Taylor expansion of f (x k ) and f (x k ) around ξ, we have and where (16), we get

Dividing (15) by
Replacing (17) in the first step of family (1), we have Using Taylor expansion again, similar to (15), to expand f (y k ) around ξ, we get With (15), (18), and (19), we calculate the divided difference operator defined in (2), obtaining Then, substituting (15), (16), (18), and (20) in the second step of family (1), we have Using Taylor series once again, similar to (15), to expand f (z k ) around ξ, we get Replacing (16) and (20) in u k and v k of family (1), we have Finally, substituting (16) and (21)-(24) in the third step of family (1), we get being the error equation To cancel the factors accompanying e 4 k and e 5 k in (26), it must be satisfied that α + γ = 1, 6α − β + 5γ = 4 and 10α − β + 9γ = 8. It is easy to show that this system of equations has infinite solutions for where α is a free parameter. Therefore, replacing (27) in (26), we obtain and the proof is finished.
From Theorem 2, it follows that, if we only hold parameter α in (1), the new triparametric family of iterative methods is reduced to an uniparametric family with an order of convergence of six for any real or complex values of the parameters α, β and γ, as long as (27) is satisfied. Therefore, the iterative expression of the new uniparametric family, dependent only on parameter α and which we will call CMT(α) family, is defined as where , and k = 0, 1, 2, ...
Because of the results obtained with the convergence analysis carried out, from now on we will only work with CMT(α) family of iterative methods and, to select the best members of this family, we will use the complex dynamics tools discussed in Section 3.

Complex Dynamics Behavior
This topic refers to the study of the behavior of a rational function associated with an iterative family or method. From the numerical point of view, the dynamical properties of the referred rational function give us important information about its stability and reliability. The parameter spaces of a family of methods, built from the critical points, allow us to understand the performance of the different members of the family, helping us in the election of a particular one. The dynamical planes show the behavior of these particular methods in terms of the basins of attraction of their fixed points, periodic points, etc. A basin of attraction provides us to visually interpret how a method works based on several initial estimates.
In this section, we present the study of the complex dynamics of CMT (α) family given in (29). To do this, we construct a rational operator associated with the family, on a generic low-degree nonlinear polynomial, and we analyze the stability and convergence of the corresponding fixed and critical points. Then, we construct the parameter spaces of the free critical points and generate dynamical planes of some methods of the family for good and bad values of α, in terms of stability.

Rational Operator
The rational operator can be built on any nonlinear function; however, we construct this operator on quadratic polynomials, as the criterion of stability or instability of a method applied to these polynomials can be generalized for other nonlinear functions.
be a generic quadratic polynomial with roots a, b ∈ R. Therefore, the rational operator R α (x) associated with CMT(α) family given in (29) and applied on p(x), is with α ∈ C an arbitrary parameter.
be a generic quadratic polynomial with roots a, b ∈ R. We apply the iterative scheme given in (29) on p(x) and obtain a rational function A p,α (x) which depends on the roots a, b ∈ R and a parameter α ∈ C. Then, if we use Möbius transformation (see in [7][8][9] which only depends on an arbitrary parameter α ∈ C. Furthermore, if we factor numerator and denominator of (35), it is easy to show that for α ∈ {−77, −1, 1, 5} some roots coincide and simplify R α (x), as it is observed in Equations (31)-(34), and the proof is finished.
From Proposition 1, for four values of α the rational operator R α (x) is simpler, so there will be fewer fixed and critical points that can improve the stability of the associated methods. This will be seen in Sections 3.2 and 3.3.

Analysis and Stability of Fixed Points
We calculate the fixed points of the rational operator R α (x) given in (30) and analyze their stability.

Proposition 2.
The fixed points of R α (x) are the roots of the equation R α (x) = x. That is, x = 0, x = ∞ and the following strange fixed points: • ex i (α) that correspond to the 10 roots of polynomial x 10 + 6x 9 The total number of different fixed points varies with the value of α: • If α ∈ C and α / ∈ {−77, −1, 1, 5, 307}, then R α (x) has 13 fixed points.
The pairs of conjugated strange fixed points, satisfying ex i = 1 ex j for i = j, are ex 2 and ex 3 , ex 4 and ex 5 , ex 6 and ex 9 , ex 7 and ex 8 , and ex 10 and ex 11 .
From Proposition 2, we establish there is a minimum of 11 and a maximum of 13 fixed points. Of these, 0 and ∞ correspond to the roots of the original quadratic polynomial p(x), and the strange fixed point ex 1 = 1 (if α = −77) corresponds to the divergence of the original method, before Möbius transformation. (ii) If 384 77 + α > 1, then ex 1 is a repulsor.
The repulsive fixed points, which always satisfy |R α (x)| > 1, are the strange fixed points ex 2 and ex 3 .
It is clear that 0 and ∞ are always superattracting fixed points, but the stability of the rest of fixed points depends on the values of parameter α. From Proposition 3, there are 6 strange fixed points that can become superattractors for certain values of α. This means that there would be a basin of attraction of the strange fixed point and it could cause the method not to converge to the solution. Figure 1 shows the stability surface of the strange fixed point ex 1 . In this figure, the zones of attraction (yellow surface) and repulsion (gray surface) are observed, being the first one much greater than the second one. Note that for values of α inside disk, ex 1 is a repulsor; and, for off-disk values of α, ex 1 is an attractor. Therefore, it is in our interest to always work inside the disk because the strange fixed point ex 1 = 1 comes from the divergence of the original method and, therefore, it is better for the performance of the iterative method that the divergence is repulsive.
From Proposition 2, the study of the stability of strange fixed points is reduced by a half. This is because each pair of conjugated strange fixed points exhibits the same stability characteristics. Furthermore, due to Proposition 3, ex 2 and ex 3 are always repulsors regardless of the value of α. Thus, Figure 2 shows the stability surfaces of the remaining 8 strange fixed points, which can be attracting or repulsive depending on the value of α, for analysis.

Analysis of Critical Points
We calculate the critical points of the rational operator R α (x) given in (30).
The pairs of conjugated free critical points, satisfying cr i = 1 cr j for i = j, are cr 2 and cr 3 , cr 4 and cr 5 , cr 6 and cr 7 , and cr 8 and cr 9 .
From Proposition 4, we establish there is a minimum of 7 and a maximum of 11 critical points. Of these, 0 and ∞ correspond to the roots of the original quadratic polynomial p(x). The free critical points cr 1 = −1, cr 2 = −i, and cr 3 = i are pre-images of the strange fixed point ex 1 = 1. Therefore, the stability of cr 1 , cr 2 , and cr 3 will correspond to the stability of ex 1 (see Section 3.2). Moreover, the dynamical study of the free critical points is reduced by a half because each pair of conjugated free critical points presents the same stability characteristics. This will be seen in Section 3.4.

Parameter Spaces
The dynamical behavior of operator R α (x) depends on the values of parameter α. The parameter space is defined as a mesh in the complex plane where each point of this mesh corresponds to a different value of α. Its graphical representation shows the convergence analysis of a method of CMT (α) family associated with this α using one of the free critical points cr(α) given in Proposition 4 as initial estimate. The resulting graphic is made in Matlab R2020a programming package with a resolution of 1000 × 1000 pixels. If a method converges to any of the roots starting from cr(α) in a maximum of 80 iterations with a tolerance of 10 −3 , the pixel is colored red; in other cases, the pixel is colored black.
Each value of α that belongs to the same connected component of the parameter space results in subsets of schemas with similar dynamical behavior. Therefore, it is interesting to find regions of the parameter space as stable as possible (red regions), because these values of α will give us the best members of the family in terms of numerical stability.
CMT(α) family has a maximum of 9 free critical points. Of these, cr 1 , cr 2 , and cr 3 have the same parameter space which corresponds to the stability surface of ex 1 (see Figure 1), because they are pre-images of this point. The remaining free critical points, cr 4 to cr 9 , are conjugated in pairs (see Proposition 4), which gives rise to 3 different parameter spaces. These parameter spaces, named P 1 (for x = cr 4 , cr 5 ), P 2 (for x = cr 6 , cr 7 ), and P 3 (for x = cr 8 , cr 9 ), are shown in Figure 3.
From Figure 3b,c, we observe that the parameter spaces P 2 and P 3 have similar characteristics; then, we can select any of them for analysis.
On the one hand, if we choose values of α inside the stability regions (red regions) of the parameter spaces, for example, α = −1, 0, 1, the methods associated with these parameters will show good dynamical behavior in terms of numerical stability. Furthermore, note that these particular values of α simplify the iterative scheme of CMT(α) family given in (29) by canceling a term in its third step. This is especially useful to improve the computational efficiency of the associated method because the processing times required to reach the solution are reduced (see Section 4).
On the other hand, if we choose values of α outside the stability regions (black regions) of the parameter spaces, for example α = −300, 200, 400, the methods associated with these parameters will show poor dynamical behavior in terms of numerical stability.
The methods associated with the values of α treated above are discussed in Section 3.5.

Dynamical Planes
We begin this section by presenting how we generate a dynamical plane that will allow us to see the stability of a method for a specific value of α. This is defined as a mesh in the complex plane where each point of this mesh corresponds to a different value of the initial estimate x 0 . Its graphical representation shows the convergence of the method to any of the roots starting from x 0 with a maximum of 50 iterations and a tolerance of 10 −3 . Fixed points are illustrated with a white circle " ", critical points with a white square " " and attractors with a white asterisk " * ". Moreover, the basins of attraction are depicted in different colors. The resulting graphic is made in Matlab R2020a with a resolution of 1000 × 1000 pixels.
Here, we study the stability of some CMT(α) family methods through the use of dynamical planes. We will consider the methods proposed in Section 3.4 for values of α inside and outside the stability regions of the parameter spaces.
On the one hand, examples of methods inside the stability region are given for α = −1, 0, 1. Their dynamical planes with some convergence orbits in yellow are shown in Figure 4. Note that all three methods present only two basins of attraction associated with the roots: the basin of 0 colored in orange and the basin of ∞ colored in blue. Furthermore, there are no black areas of non-convergence to the solution. Consequently, these methods show good dynamical behavior: they are very stable. Of these methods, the best member of CMT(α) family is for α = 1, as it has fewer strange fixed points and free critical points.   On the other hand, examples of methods outside the stability region are given for α = −300, 200, 400. Their dynamical planes with some convergence orbits in yellow are shown in Figure 5. Note that all three methods present more than two basins of attraction, that is, there are other basins of attraction that do not correspond to the roots. The basins of 0 and ∞ are colored in orange and blue, respectively, and the other basins are colored in black, red, and green. Figure 5a shows the convergence to an attracting periodic orbit of period 2.

Numerical Results
Here, we perform several numerical tests in order to check the theoretical convergence and stability results of CMT(α) family obtained in previous sections. To do this, we use some stable and unstable methods of (29). These methods are applied on five nonlinear test functions, whose expressions and corresponding roots are Thus, we performed two experiments. In a first experiment, we carried out a efficiency analysis of CMT(α) family through a comparative study between one of its stable methods and five different methods given in the literature: Newton of order 2, Ostrowski of order 4, and three other methods of order 6 proposed by Alzahrani et al. in [10] (ABA), Chun and Ham in [11] (CH), and Amat et al. in [12] (AHR). In a second experiment, we carried out a stability analysis of CMT (α) family using six of its methods obtained with three good and three bad values of parameter α, in terms of stability.
In the development of the numerical tests we start the iterations with different initial estimates: close (x 0 ≈ ξ), far (x 0 ≈ 10ξ), and very far (x 0 ≈ 100ξ) to the root ξ, respectively. This allows us to measure, to some extent, how demanding the methods are relative to the initial estimation for finding a solution.
The calculations are developed in Matlab R2020a programming package using variable precision arithmetics with 200 digits of mantissa. For each method, we analyze the number of iterations (iter) required to converge to the solution, so that the stopping criteria |x k+1 − x k | < 10 −100 or | f (x k+1 )| < 10 −100 are satisfied. Note that |x k+1 − x k | represents the error estimation between two consecutive iterations and | f (x k+1 )| is the residual error of the nonlinear test function. This stopping criterium does not need the exact solution, on the contrary of absolute error, and differs from recent ones as CESTAC (see in [13]) in the absence of additional calculations or functional evaluations, as f (x k+1 ) is needed for the following iteration and its absolute values is an efficient control element of the proximity to the exact root, where f is zero. Indeed, although a precision of one hundred exact digits is not usually necessary in the applications, we employ this value in the stopping criterium as it is useful to check the robustness and effectiveness of the numerical methods.
To check the theoretical order of convergence (p), we calculate the approximate computational order of convergence (ACOC) given by Cordero and Torregrosa in [14]. In the numerical results presented below, if the ACOC vector inputs do not stabilize their values throughout the iterative process, it is marked as "-"; and, if any of the methods used does not reach convergence in a maximum of 50 iterations, it is marked as "nc".
To illustrate the computational efficiency of each used method, the processing time (tcpu) in seconds required by the iterative scheme to converge to the solution is measured. This value is determined as the arithmetic mean of 10 runs of the method.

First Experiment: Efficiency Analysis of CMT (α) Family
In this experiment, we carried out a comparative study between a stable method of CMT (α) family and the methods of Newton, Ostrowski, ABA, CH and AHR, in order to contrast their numerical performances in nonlinear equations. We consider as a stable member of CMT (α) family the method associated with α = 1, that is, CMT(1). Tables 1-3 we show the numerical results of the six known methods, considering close, far, and very far initial estimates. Furthermore, in Figure 6 we show graphics that summarize these results for the number of iterations (iter) and the processing time (tcpu).

Thereby, in
Therefore, from the results of the first experiment we conclude that CMT(α) family has an excellent numerical performance considering a stable member (α = 1) as a representative. This conclusion has been made based on the following aspects from Tables 1-3: CMT (1) method has the lowest error and lowest number of iterations (iter). However, the mean of the execution time (tcpu) varies according to the nonlinear test function used and the inherent complexity that the iterative scheme of the method presents on the nonlinear function. In several cases, the tcpu of the CMT (1) method is significantly lower than the 6th order ABA, CH and AHR methods. The theoretical convergence order is also verified by the ACOC, which is close to 6.

Second Experiment: Stability Analysis of CMT(α) Family
In this experiment, we carried out a stability analysis of CMT (α) family considering some values of α inside the stability regions of the parameter spaces (α = −1, 0, 1) and outside of them (α = −300, 200, 400).
Thus, in Tables 4-9 we show the numerical performance of iterative methods associated with these values of α for close, far, and very far initial estimations. The results for α = 1 were already presented in the first experiment; however, these are presented again due to the different conditions in which each experiment was performed.
On the one hand, from Tables 4-6 we observe that the methods associated with α = −1, 0, 1 always converge to the solution, although the number of iterations (iter) needed differs for any initial estimate and nonlinear test function. Thus, in estimations close to the root, the methods converge to ξ with a minimum iter of 3 and a maximum of 7. When the initial guess is far from the root, they converge to ξ with a minimum iter of 4 and a maximum of 22. When the starting estimations are very far from the root, the iterative schemes converge to ξ with a minimum iter of 6 and a maximum of 37.  On the other hand, from the results shown in Tables 7-9, we see that the methods associated with α = −300, 200, 400 do not always converge to the solution, confirming the conclusions obtained in the dynamical analysis. The convergence highly depends on the initial estimation and the nonlinear test function used. Thus, for estimations close to the root, these methods do not converge to the solution in up to 2 test functions. Moreover, for estimations far and very far from the root, they do not converge to the solution even for any function.
Consequently, we conclude that the methods for α = −1, 0, 1 are stable, have the lowest processing times (tcpu), and always converge to the solution for any initial estimate and nonlinear test function used. The methods for α = −300, 200, 400 are unstable, chaotic, have the highest tcpu, and tend not to converge to the solution according to the initial estimate and the nonlinear test function used. With this, the theoretical results obtained in previous sections about the dynamical behavior of CMT(α) family are verified. In the dynamical study, we constructed parameters spaces of the free critical points of the rational operator associated with the uniparametric family. These parameter spaces allowed us to understand the performance of the different members of the family, helping us to choose stable (for α = −1, 0, 1, ...) and unstable (for α = −300, 200, 400, ...) methods. Furthermore, we generated dynamical planes to show the behavior of these particular methods.
From numerical results, the order of convergence is verified by the ACOC, which is close to 6. The CMT (α) family proved to have an excellent numerical performance considering stable members as representatives. In general, this family has low errors and number of iterations to converge to the solution. However, the processing time (tcpu) varies depending on the nonlinear test functions used and the inherent complexity that the iterative schemes of the methods present when they are applied to said functions. In several cases, the tcpu of stable methods is significantly lower than other sixth-order methods developed so far. Furthermore, the methods for α = −1, 0, 1 proved to be stable, have the lowest tcpu, and always converge to the solution for any initial estimate and nonlinear test function used. The methods for α = −300, 200, 400 proved to be unstable, chaotic, have the highest tcpu, and tend not to converge to the solution according to the initial estimate and the nonlinear test function used. This verifies the theoretical results obtained in convergence analysis and dynamical study of CMT(α) family.