Isonormal surfaces: A new tool for the multidimensional dynamical analysis of iterative methods for solving nonlinear systems

The dynamical behavior of the rational vectorial operator associated with a multidimensional iterative method on polynomial systems gives us interesting information about the stability of the iterative scheme. The stability of fixed points, dynamic planes, bifurcation diagrams, etc. are known tools that provide us this information. In this manuscript, we introduce a new tool, which we call isonormal surface, to complement the information about the stability of the iterative method provided by the dynamical elements mentioned above. These dynamical instruments are used for analyzing the stability of a parametric family of multidimensional iterative schemes in terms of the value of the parameter. Some numerical tests confirm the obtained dynamical results.


INTRODUCTION
The design of iterative processes for solving nonlinear systems, F(x) = 0, F ∶ D ⊑ R n → R n , with n unknowns and n equations, is an interesting challenge of numerical analysis (see, for example, 1,2 ). Many problems in Science and Engineering need the solution of a nonlinear system in any step of the process.
In the absence of analytical methods to solve this type of problem, the solutions of these systems must be approached through iterative fixed-point schemes, , k = 0, 1, … , being x (0) the initial estimation.
The best known method for approximating a solution ∈ D is Newton's procedure, being F ′ (x (k) ) the Jacobian of function F evaluated in the kth iteration. This method reaches quadratic convergence, under some conditions.
In recent literature in this area of research, we can find other methods that reach higher orders (see, for example, previous studies [3][4][5][6][7], and they are designed by using different techniques like: Adomian decomposition, [8][9][10] multidimensional Steffensen-type schemes, [11][12][13][14][15][16] and weight function techniques. [17][18][19] Once a class of methods has been designed, it is interesting to carry out a multidimensional real dynamical study (see 20,21 ) in order to obtain the values most suitable of the parameters for setting stable schemes.
In 22 presented a family of iterative schemes for solving a nonlinear system F(x) = 0. Being F a real Fréchet differentiable function and H ∶ R n×n → R n×n a matrix weight function with variable t (k) = I−[F ′ (x (k) )] −1 [x (k) , (k) ; F], where I denotes the identity matrix of size n × n, the divided difference operator of F on R n is a mapping [·, ·; F] ∶ Ω × Ω ⊂ R n × R n → (R n ) (see 23  Then, the following three-step class of iterative methods is designed: x (k+1) = z (k) − H(t (k) )[F ′ (x (k) )] −1 F(z (k) ), k ≥ 0.
Let X = R n×n be the Banach space of real square matrices of size n × n, and function H : X → X can be defined in the way that its Fréchet derivatives satisfies where (X) denotes the set of linear operators defined in X. When k tends to infinity, variable t (k) tends to the zero matrix 0. So, there exist real H 1 , H 2 such that H can be expanded around 0 as The following result was proven in, 22 regarding the convergence analysis of (1). [ (H 2 2 − 22H 2 + 120)C 5 2 + (−24 + 2H 2 )C 2 2 C 3 C 2 + (−20 + 2H 2 )C 3 C 3 2 + 4C 2 3 C 2 ] e (k) 6 + O(e (k) 7 ), We select a particular weight function H(t) = I + 2t + 1∕2 t 2 where ∈ R is a free parameter. This function satisfies the hypothesis of Theorem 1, and the resulting one-parameter family, called PSH 6 , is In Capdevila et al, 22 we saw that in PSH 6 family, the method obtained with = 0 was the most computationally efficient, although in practice, other parameter values provided even better results in terms of precision or numerical estimation of the order of convergence, with the same number of iterations.
This fact leads us to deepen the analysis of the stability of the schemes belonging to that family, with the aim of establishing which are the values that can provide better results and more stable methods. In the process, we will define a new tool that has proven to be very useful in achieving this goal: the isonormal surface. With this tool that we will introduce in Section 2, we will establish the regions of the plane that have the same qualitative behavior (set of initial estimates that converge to periodic orbits of the same period, regions of chaotic behavior, dense orbits, etc.) which are identified with different colors for any value selected from the bifurcation diagrams.
Thus, the set formed by the Feigenbaum diagram (or the parameter line) and the isonormal surface provides us accurate and reliable information about the behavior of the members of any kind of iterative methods for solving systems of nonlinear equations. Now, let us denote as G a vectorial rational operator. We define the orbit of x (0) as a set of the successive images by G, that is, {x (0) , G(x (0) ), … , G m (x (0) )}. Indeed, the dynamical behavior of a point x ∈ R n can be classified examining its asymptotic performance, thus the point x * such that G(x * ) = x * is denominated a fixed point of G. In the same way, a periodic point Some results about stability of fixed points are summarized in the next result ( 24 p. 558 ).
Those fixed point whose eigenvalues satisfy | j | ≠ 1 are called hyperbolic, and they are saddle hyperbolic points if some of their eigenvalues hold | j | > 1 and the rest | i | < 1, i, j ∈ {1, 2, … , n}. If a fixed point is not a zero of F ∶ D ⊆ R n → R n , then it is called a strange fixed point, and its stability is classified according to Theorem 2.
Additionally, if x * is an attracting fixed point of the rational function G, we can define its basin of attraction (x * ) as the set of pre-images of any order such that And finally, x c ∈ R n is called critical point if all eigenvalues of G ′ (x c ) are equal to zero. These points are important because a classic result states that in every basin of attraction of a fixed or periodic point, there is at least one critical point. In fact, it yields in the immediate basin of attraction, that is, in the same connected component as the periodic or fixed point. Those critical points different of the roots of F(x) = 0 are called free critical points.
The manuscript is structured as follows: In Section 2, the dynamical behavior of the rational vectorial operator obtained when family (1) is applied on a polynomial system is analyzed. In Section 3, some numerical tests are presented to confirm the results obtained in the dynamical section. Finally, some conclusions and the references used in this paper are shown. 6 We analyze the performance of the vectorial rational function resulting from applying the iterative family PSH 6 on the function

DYNAMICAL ANALYSIS OF THE CLASS PSH
In order to select the most stable members of this class of iterative algorithms, we will study the existence of fixed points different from the roots we are looking for with attracting character or another attracting elements that can be considered pathological.
By applying the iterative expression of PSH 6 on the polynomial system p(x) = 0, we obtain its rational multidimensional operator associated U(x, ) = {u 1 (x, ), u 2 (x, )}, whose jth coordinate is From the last expression, it is possible to formulate the following result about the stability of fixed points related to the iterative family (2). Theorem 3. The rational function U(x, ) associated to the family of iterative methods PSH 6 has, as superattracting fixed points, (1,1), (1 and also  is composed of (± 1, l j ) and (l i , ± 1). Therefore, the number of fixed points included in  depends on : i) There is no real strange fixed point for ∈ (n * , 0); however, if ∈ (− ∞, n * ) ∪ (0, m * ), the total amount of strange fixed points is 12, being n * ≈ −93.210875 and m * ≈ 327.44373142 the only real roots of polynomial n(t) = t 3 − 112t 2 + 3190784 + 15104t and m(t) = 100000t 9 − 33232500t 8 + 162157876t 7 − 769658421t 6 + 2887352848t 5 −7815029120t 4 +14060488832t 3 −17146950656t 2 +12837074944t −5333716992, respectively. Indeed, four of them are repelling, while the eight elements left are saddle points. ii) If ∈ (m * , + ∞),  is composed of sixty strange fixed points, whose character depends on the value of . Therefore, two different situations can be found regarding their stability: a) When ∈ (m * , * ), being * ≈ 369.97117, then  is composed by twelve attracting, sixteen repelling and thirty-two saddle strange fixed points. b) Finally, if ∈ ( * , + ∞) then  is composed by twenty repelling and forty saddle fixed points.
Proof. Let us remark that, due to the fact that the polynomial system has separated variables, the coordinates of U(x, ) have the same expression with exception of the sub-indices. So, being a fixed point the solution of equation u j (x, ) = x j , = 1, 2, we obtain an equivalent form Then, we state that the values x = ±1 satisfy this expression, and (1, 1), (1, − 1), (− 1, 1), (− 1, − 1) are fixed points of the rational operator U(x, ) and the roots of p(x), simultaneously. To analyze their stability, we calculate the associate Jacobian matrix U ′ (x, ) of the rational multidimensional operator, with diagonal shape It is straightforward that J (±1, ) = 0 for any and = 1, 2. So, the roots of p(x) are superattracting fixed points as the eigenvalues of U ′ ((± 1, ± 1), ) are zero.
The rest of fixed points can be found through l(t, ) by means of s = t 2 . A ninth degree polynomial is then obtained The fixed points of U(x, ) must have real components, which are defined as ± √ L i , being L i any real and positive root of L(s). Then, the number of elements of  depends on the number of roots of L(s) that can be simultaneously real and positive, as well as their combinations with ±1. It can be checked that there is no real positive root of L(s) for ∈ (n * , 0). Then, In this case, the roots , and the members of  are {(l 1 , l 1 ), (l 1 , l 4 ), (l 4 , l 1 ), (l 4 , l 4 ), (l 1 , ± 1), (± 1, l 1 ), (l 4 , ± 1), (± 1, l 4 )}. The information about stability of these fixed points of U(x, ) can be inferred from the analysis of the absolute value of the eigenvalues of matrix U ′ ((l i , l k )), ), i, k ∈ {1, 4}; these functions of are called stability functions of the respective fixed points. Due to the nature of the polynomial system, the eigenvalues ; however, if any of the components of the fixed point is ±1, the corresponding eigenvalue is always null. We can see in Figure 1 the values of | | associated to the Jacobian matrix and evaluated at pairs composed of l 1 and l 4 . Let us remark that they are greater than one; therefore, the behavior of the strange fixed points ( are saddle points as one of the eigenvalues is zero and the another one is greater than one. ii) There exist three real positive roots L 1 , L 2 , and L 3 for ∈ (m * , + ∞), being n * ≈ −93.210875 the only real root of polynomial n(t) = t 3 − 112t 2 + 3190784 + 15104t and m * ≈ 327.44373142 the only real root of m(t) = are denoted by l i for i ∈ {1, 2, 3, 4, 5, 6}. So, the set of all strange fixed points is obtained combining in pairs l i , i = 1, 2, … , 6 with themselves and with 1 and −1. In this case, there exist values of that allow some of the fixed points to be attracting. To calculate them, we solve the equation | ((l i , l k ), ( ))| = 1, for i, k ∈ {1, 2, 3, 4, 5, 6} and = 1, 2. We found that only strange fixed points (l i , l k ) with i, k ∈ {3, 6} satisfy this equation at = m * and = * ≈ 369.97117. a) Therefore, they are not hyperbolic for these values of and attracting for ∈ (m * , * ), as can be seen in Figure 2A. Moreover, the stability of fixed points (l i , ± 1) or (± 1, l k ), for i, k ∈ {3, 6} and ∈ (m * , * ), corresponds to a Jacobian matrix whose eigenvalues take values | 1 ((l i , ± 1), ( )) |<1 and | 2 ((l i , ±1), ( ))| = 0 (respectively, | 1 (±1, (l k ), ( ))| = 0 and | 2 (± 1, (l k ), ( )) | < 1) if ∈ (m * , * ), so they are attracting, or if ∈ ( * , + ∞), then | 1 ((l i , ± 1), ( )) | > 1 and | 2 ((l i , ±1), ( ))| = 0 (respectively, | 1 (±1, (l k ), ( ))| = 0 and | 2 (± 1, (l k ), ( )) | > 1), and these components of  are saddle points. b) For ∈ ( * , ∞), all the 36 pairs (l i , l j ) are repelling fixed points, as absolute values of | j ((l i , l k ), ( )) |>1, for i, k ∈ {1, 2, 3, 4, 5, 6} and = 1, 2. The stability of one of these fixed points can be inferred from Figure 2B. Now, the stability of strange fixed points is determined, and we know the range of values of that can be used in order to avoid attracting strange points. However, there are other attracting elements that should be avoided in order to assure an stable performance of the iterative methods. As we know that critical points appear in each attracting basin of attraction, we analyze the existence of critical points different from the roots of p(x) = 0 (called free critical points), depending on the value of the parameter, in the following statement.

Theorem 4.
Let  be the collection of all free real critical points of the rational operator U(x, ) associated with the iterative family (2). Then  is composed of the pairs (c i , c j ), (c i , ± 1), and (± 1, c j ) for i, j≤10 whose entries different from ±1 are the real roots of polynomial r(t) described in (6). The amount of free critical points depends on the value of parameter > 0, as the roots must be real. Therefore, the composition of set  is as follows: being r(t) = (−1003520 + 15104 − 112 2 + 3 )t 10 + (319488 + 170752 − 640 2 + 15 3 )t 8 By definition, the critical points are found by solving the equation (x, ) = 0, for j ∈ {1, 2}. It is clear that change t 2 = s reduces to the half the degree of polynomial r(t). Therefore, its positive real roots, denoted by C j , will derive the components of the free critical points, obtained as c k = ± √ C . Again, the amount of positive C j depends on the value of parameter : i) Forcing the roots of r(s) to be real, no result is found for negative values of . Moreover, there exist three real roots for ∈ (0, t * ), being t * ≈ 1.578466 the only real root of polynomial U(t), but only one of them, C 1 , is positive. Then,  is composed of (c 1 , c 1 ), (c 1 , c 4 ), (c 4 , c 1 ), and (c 4 , c 4 ), where c 1,4 = ± √ C 1 , and also of (c j , ± 1) and (± 1, c j ), j ∈ {1, 4}. ii) If ∈ (t * , 80), there exist only three real positive roots of polynomial r(t), C 1 , C 2 , and C 3 . So,  is composed of elements of kind (c i , c j ), where i, j ∈ {1, 2, … , 6}, being c 1,4 = ± √ C 1 , c 2,5 = ± √ C 2 , and c 3,6 = ± √ C 3 . Moreover, mixed points (± 1, c j ) and (c i , ± 1) where i, j ∈ {1, 2, … , 6} also belong to . iii) Finally, for ∈ (80, + ∞), C 1 and C 2 are the only two positive real roots of r(t).
, holding 32 different free critical points.
From Theorems 3 and 4, we state that for ∈ (n * , 0), there exist neither free critical points nor strange fixed points. Hence, only stable dynamical behaviors can be found for values of the parameter in this interval, that is, for iterative methods of the class PSH 6 with ∈ [n * , 0]. In Figure 3, some dynamical planes for values of inside this interval can be observed. These figures have been obtained following the routines described in. 25 We have built a mesh with step equal to 0.01, and every initial estimation is iterated 100 times with an estimation of the error lower than 10 −3 as a stopping criterium.
Each point of the mesh, used as initial estimation, is colored depending on the root (represented with circles) it converges to. The color is brighter when lower is the number of iterations needed. If all the iterations are completed and no convergence to any roots is reached, then the point is painted in black. We can see in Figure 3 the basins of attraction of all roots and how the black areas of no convergence (divergence or maximum number of iterations reached) are narrower as the value of increases. All the dynamical planes in this paper have been made with the same features.

Finding chaos: New and known tools
Once the most stable area has been detected, our aim is to find the worst areas in terms of convergence to attracting elements different from the roots of the nonlinear system or even chaos. To get this aim, we take into account the results of Theorems 3 and 4: for ∈ (t * , 80) and ∈ (80, + ∞), a big amount of free critical and fixed points have been found, some of the later ones with attracting behavior. The analysis of the performance of the rational function in these areas will yield a rich dynamical performance.
Due to the existence of critical points in the immediate basin of attraction of any attracting fixed or periodic point (see 26 ), the analysis of the asymptotic behavior of free critical points gives us relevant information about the performance of the rational operator and the related class of iterative methods. To study the orbits of critical points, we use the parameter line, firstly introduced in. 27 To plot it, we define a mesh of 500 × 500 points in a certain domain of parameter . A point of the mesh is painted in red color if the critical point, evaluated at this value of and used as initial estimation, converges to any of the roots of the polynomial system before a maximum of 200 iterations; in another case, it is painted in black color. Each one of these red or black points is thickened by means of a multiplication by the unit interval [0, 1]. The tolerance for the error estimation is equal to 10 −3 . Figure 4A,B shows the limit behavior of the free critical points corresponding to ∈ (t * , 80). In Figure 4A, desirable behavior is observed, and free critical points (c i , c j ) with i, j ∈ {1, 2, 4, 5} converge to any of the roots of polynomial system; meanwhile, in Figure 4B, unstable behavior is detected for (c i , c j ) with i, j ∈ {3, 6} around the values ≈ 74, ≈ 77 and ≈ 79 in the parameter line. We make a first approach to the analysis of the unstable behavior in these three dark zones through the Feigenbaum diagrams. These bifurcation diagrams are obtained using each one of the free critical points of the rational operator as starting point in a mesh of 3000 subintervals for ∈ (t * , 80) and observing their performance in the last 100 of a total amount of 1000 iterations. The bifurcation diagram of critical point (c 3 , c 6 ) is shown in Figure 5.
Due to the nature of polynomial system p(x), that has separate variables, the coordinate functions of the rational operator have the same shape and are symmetrical. Taking advantage of this feature, we have plotted in red color the x 1 dimension and in blue color x 2 , and for the details, we have chosen the positive part in all the cases. In Figure 5A, we observe the convergence to the roots of p(x) far away of ≈ 74.05 where the unstable behavior lies. Analogous performances appear in Figure 5C,E, for different scales. On the other hand, in the details given by Figure 5B,D,F, the same pattern of bifurcations as period-doubling cascades, periodical orbits, and chaotical behavior can be observed around ≈ 74.05, ≈ 77.01, and ≈ 79.36, respectively. On the other hand, in Figure 4C,D, the parameter lines corresponding to the limit performance of the free critical points with components different from ±1 are shown, for ∈ (80, + ∞). Although they have been plotted for pairs (c 1 , c 3 ) ( Figure 4C) and (c 2 , c 4 ) ( Figure 4D), the rest of pairs have the same behavior. Let us remark that for a wide interval around = 150, only stable behavior is found, and for > m * , there is no possibility of stable performance, as all critical free points show black color in the parameter lines. To deepen in these unstable performances, we propose a new tool, based on the norm of the nonlinear function at the orbit of any initial estimation, for a fixed value of .

The new tool: Isonormal surfaces
In Figure 6, the dynamical planes related with the shaded areas found in the parameter lines appear (see Figure 4B), where the critical points show the unstable behavior for ∈ (t * , 80).
In Figure 6A, white squares represent the free critical points. Some of them are located in the black zone where unstable behavior takes place and the rest lie in the basin attraction of the roots of p(x). In Figure 6B, an infinite number of connected components of the basins of attraction of the roots can be observed and also the proximity between four critical and four repelling (red circles) that lie in the Julia set, with coordinates { (− 0.618, − 0.618), (− 0.618, 0.618), (0.618, − 0.618), (0.618, 0.618)} and To gain more information about the dynamical performance into the black areas, in the designed tool that we call isonormal surface, every initial estimation of a dynamical plane for a fixed value of is iterated 1000 times, and the Euclidean norm of the last iteration is calculated. Therefore, every point is painted in a different color depending on this calculated value of the norm. As a result, certain patterns emerge; see Figure 7. In this picture, magenta area corresponds to norm equal to √ 2, that is, the basins of attraction of the roots of system p(x); higher values of the norm yield to orange and yellow areas that correspond to the basins of attraction of periodical or dense orbits; the yellow area corresponds with the basin of a diagonal orbit with period 16 that has been plotted with blue points and lines in Figure 7A Figure 5B for the specific value of = 74.056. A symmetrical performance can be found in the other diagonal corners. In Figure 7B, a detail of these yellow areas can be observed, where also small regions of convergence to the roots (in magenta color) or orange areas are found. These orange areas correspond to the basin of attraction of vertical and horizontal periodical orbits of the same period. In Figure 7C,D, the isonormal surface for = 74.06 and a detail, respectively, are represented. With this slight modification of the value of , the yellow area of periodic behavior has become into a darker area including an strange attractor: a dense orbit. In Figure 7D, the trajectory of the point (− 0.57638250, 0.57679642) already contained in the last periodic orbit analyzed has been plotted, and it can be observed that the iterates are filling the corner areas as far as the maximum number of iterations is reached. For this value of , this dense orbit not only fills the diagonal areas shown in Figure 7D but also the small areas around the points (− 0.5, 0.5) and (0.5, − 0.5) that were also basin of attraction of the periodic orbit for = 74.056.  Figure 4C,D) that the free critical points converge to attracting elements different from the roots. In Figure 8, the dynamical planes for the values = 350 and = 370 are plotted. In the first of them, the basins of attraction of the attracting fixed points can be observed in different colors. In these basins, the attracting strange fixed points appear as red circles. When = 370, these points are no longer attracting but repelling. So, they lay in Figure 8B in the black area of the dynamical plane. However, the strange fixed points that are repelling should appear in the Julia set of the dynamical plane, that is, in the boundary between different basins of attraction that do not appear as do not correspond to the roots of the system. These attracting elements must be found in order to fully understand the performance of the methods. These two cases are representative of the unstable performance of the methods corresponding to values of > 80.
This valuable information can be found in the isonormal surfaces shown in Figure 9. Let us remark that 12 of the strange fixed points that were attracting for ∈ (m * , * ), after the bifurcation value = * , have not only become repelling but have created a periodic orbit of period 2, whose elements are very close to the fixed point. Due to the symmetry of the system, the value of the norm coincides when a change of sign happens in the components of the limit points. So, we The reason is the high number of iterations used in the isonormal surface, 1000. As the Julia set is the unstable manifold, an orbit of an arbitrary initial point close to it will be "repelled" towards the closer attracting element, in this case any of the periodical orbits. So, even the repelling strange fixed points appear in the isonormal surface as initial guesses that converge to the periodical orbits.

NUMERICAL RESULTS
In this section, the numerical performance depending on of the studied class of iterative methods PSH 6 is tested, using the dynamical information deduced in the previous section. The information regarding its behavior in comparison with other methods can be found in 22 but using values of arbitrarily chosen.  Figure 6). 3. Also unstable performance is expected for = 330, = 350, and = 369 (see Figure 8A), belonging to the interval (327.44373, 369.97117). 4. Finally, for = 370, = 380, and = 390 (see Figure 8B), chosen inside the interval (369.97117, ∞), unstable behavior is supposed to be found, too.
For values of in groups 2-4, a high number of critical points and strange fixed points of different character can be found according to the stability analysis previously realized. To perform the numerical experiments, we have employed MATLAB computer algebra system with 2000 digits of mantissa in variable precision arithmetics. The stopping criterion used is ||x (k + 1) − x (k) || < 10 −200 or ||F(x (k + 1) ) || < 10 −200 ; the initial values employed and the searched solutions are symbolized as x (0) and , respectively. The implementation of our method involves the evaluation of a divided difference operator that is calculated by using the first-order estimation of the Jacobian matrix whose elements are (see 23 ) [ , where f i , i = 1, 2, … , n, are the coordinate functions of the nonlinear function F describing the system to be solved.
We have tested the different iterative methods for solving two nonlinear systems. The information displayed in the corresponding tables is organized as follows: "Group" and show the number associated to the interval that includes the particular value of the parameter used; k is the number of iterations needed, whose maximum value is 25 ("NaN" appears in the table if undefined numeric results appear, such as 0/0); the value of the stopping residuals ||x (k + 1) − x (k) || and ||F(x (k + 1) )|| at the last iteration and also the approximated computational order of convergence (defined in), (if the value of for the last iterations is not stable, then "-" appears in the table). In this way, once the process has stopped without reaching the maximum number of iterations, it can be checked if the iterative process has reached the root (||F(x (k + 1) ) || < 10 −200 is achieved) or it is only a very slow convergence with a no significant difference between the two last iterates (||x (k + 1) − x (k) || < 10 −200 but ||F(x (k + 1) ) || > 10 −200 ). Example 1. The first nonlinear system, whose solution is = (0.0, 2.0) T , is described as The initial estimation is x (0) = (0.5, 3.0) T , and the results appear in Table 1.
In Table 1, the results of PSH 6 on F 1 (x 1 , x 2 ) are shown. A good performance is obtained for the iterative schemes in group 1; all of them converge to the solution in five iterations with very small residual ||F(x (k + 1) )||. As it is expected, after the maximum number of iterations, no convergence is observed for iterative methods belonging to groups 3 and 4, as can be observed by the difference between last two iterates (||x (k + 1) − x (k) ||) and the residual (||F(x (k + 1) )||) at the 25th iteration. In this case, the performance of the members of second group is unexpectedly good, and the computational approximated order of convergence is unstable for all groups.
As we expected, a good performance is obtained for the members of the family given by group 1; all of them converge to the solution in few iterations, being the residual null in all cases. The approximated computational order of convergence is stable enough to be near the theoretical one in two cases. Regarding the second group, our expectations have been fulfilled, and no convergence to the solution is obtained after 25 iterations. Similar results are obtained for groups 3 and 4.

CONCLUSIONS
When a family of iterative methods is applied on a real vectorial nonlinear problem, a dynamical analysis of the class is important to avoid chaotical or unstable performance. In this work, we have made a deep dynamical analysis to select the most stable elements of a parametric family of efficient sixth-order schemes. To get this aim, a new tool has been designed, the isonormal surface, that allows us to detect the areas of a dynamical plane with the same kind of unstable behavior and to find areas of chaotic performance, for a fixed value of the parameter. Once the characterization of the members of the class of iterative methods has been made in terms of their qualitative properties, these results have been checked with other nonlinear systems, with good results. In the numerical section, the effect of a bad election in the value of the parameter involved can be observed, due to the unstability or no convergence of the related iterative method in these cases.