Convergence and Stability of a Parametric Class of Iterative Schemes for Solving Nonlinear Systems

: A new parametric class of iterative schemes for solving nonlinear systems is designed. The third- or fourth-order convergence, depending on the values of the parameter being proven. The analysis of the dynamical behavior of this class in the context of scalar nonlinear equations is presented. This study gives us important information about the stability and reliability of the members of the family. The numerical results obtained by applying different elements of the family for solving the Hammerstein integral equation and the Fisher’s equation conﬁrm the theoretical results.


Design of a Parametric Family of Iterative Methods
The need to find a solutionx of equations or systems of nonlinear equations of the form F(x) = 0, where F : D ⊆ R n → R n , n ≥ 1, is present in many problems of applied mathematics as a basis for solving other more complex ones. In general, it is not possible to find the exact solution to this type of equations, so iterative methods are required in order to approximate the desired solution.
The essence of these methods is to find, through an iterative process and, from an initial approximation x (0) close to a solutionx, a sequence {x (k) } of approximations such that, under different requirements, lim k→∞ x (k) =x.
It is well known that one of the most used iterative methods, due to its simplicity and efficiency, is Newton's scheme, whose iterative expression is where F (x (k) ) denotes the derivative or the Jacobian matrix of function F evaluated in the kth iteration x (k) . In addition, this method has great importance in the study of iterative methods because it presents quadratic convergence under certain conditions and has great accessibility, that is, the region of initial estimates x (0) for which the method converges is wide, at least for polynomials or polynomial systems. Based on Newton-type methods and by using different procedures, many iterative schemes for solving F(x) = 0 have been presented in the last years. Refs. [1,2] compile many of the methods recently designed to solve this type of problem. These books give us good overviews about this area of research.
In this paper, we use a convex combination of the methods presented by Chun et al. in [3] and Maheswari in [4]. As the mentioned schemes are designed for nonlinear equations and they have as the first step Newton's method, we use the following algebraic manipulation in order to extend the mentioned schemes to nonlinear systems: f (y (k) ) f (x (k) ) = f (y (k) ) (x (k) − y (k) ) f (x (k) ) = f (y (k) ) − f (x (k) ) + f (x (k) ) (x (k) − y (k) ) f (x (k) ) = − [x (k) , y (k) ; f ] f (x (k) ) + 1.
The rest of the paper is organized as follows: Section 2 is devoted to analyze the convergence of family (2) in terms of the values of parameter γ. In Section 3, we study the dynamical behavior of the class on quadratic polynomials in the context of scalar equations. This study allows for selecting the members that are more stable in the family. In the numerical section, (Section 4), we apply the proposed class on different examples such as the Hammerstein integral equation and the Fisher's equation in order to confirm the theoretical results obtained in Sections 2 and 3. We finish the work with some conclusions and the references used in it.

Convergence Analysis
Let us consider function F : D ⊆ R n → R n , differentiable in the convex set D ⊂ R n which contains a solutionx of the nonlinear equation F(x) = 0. From the Genochi-Hermite formula (see [5]) of the divided difference operator and by performing the Taylor's expansion of F (x + th) on the point x and integrating, we obtain the following development: which we will use in the proof of the following result, when the order of convergence of family is established. Theorem 1. Let F : D ⊆ R n −→ R n be a sufficiently Fréchet differentiable function in a convex neighborhood D ofx, being F(x) = 0. We suppose the Jacobian matrix F (x) is continuous and non-singular inx. Then, taking an initial estimate x (0) close enough tox, the sequence of iterates {x (k) } generated with family (2) converges tox with the following error equation: (5) where C j = 1 j! F (α) −1 F (j) (α) ∈ L j (R n , R n ), L j (R n , R n ) being the set of j-linear functions of bounded functions, j = 2, 3, . . . and e k = x (k) −x. In the particular case in which γ = 0, the error equation is e k+1 = (4C 3 2 − C 3 C 2 )e 4 k + O(e 5 k ), (6) and so the method has an order of convergence four.
Proof. We consider the Taylor's expansion of F(x (k) ) aroundx: where In a similar way, the derivatives of F(x (k) ) aroundx take the form: From the development of F (x (k) ) aroundx, we calculate the inverse with X 2 , X 3 , X 4 and X 5 satisfying F (x (k) ) −1 F (x (k) ) = I. Therefore, (7) and (9), we obtain Then, we obtain the error equation of the first step of the parametric family (2): Substituting this expression in the Taylor expansion of F(y (k) ) aroundx, we get: Furthermore, Multiplying expressions (9) and (13), we obtain: To obtain the development of the divided difference operator of (2), we use the Taylor series expansion of (4), considering in this case x + h = y and, so, h = y − x = −F (x (k) ) −1 F(x (k) ). Therefore, substituting (8) and (10) in (4), we obtain To calculate the inverse of this operator, we search (9) and (15), we obtain B k , and using (8) and (16), we calculate B −1 K , Substituting the expressions (10), (14), (17), and (18) in the scheme (2), we get the error equation of the parametric family Finally, from the error equation, we conclude that the parametric family (2) has order 3 for all γ = 0 and order 4 for γ = 0, being in this last case the error equation In the next section, we analyze the dynamical behavior of the parametric family (2) on quadratic scalar polynomials.

Complex Dynamics
The dynamical analysis of (2) is performed throughout this section in terms of complex analysis. The order of convergence is not the only important criterion to study when evaluating an iterative scheme. The validity of a method also depends on other aspects such as knowing how it behaves based on the initial estimates that are taken, that is, how wide the set of initial estimations is for which the method is convergent. For this reason, it is necessary to introduce several tools that allow for a more exhaustive study.
The analysis of the dynamics of a method is becoming one of the most investigated parts within the study of iterative methods since it allows for classifying the different iterative schemes, not only from the point of view of their speed of convergence, but also analyzing its behavior based on the initial estimate taken (see, for example, [6][7][8][9][10][11][12][13]). This study allows for visualizing graphically the set of initial approximations that converge to a given root or to points that are not roots of the equation. In addition, it provides important information about the stability and reliability of the iterative method.
In this paper, we focus on studying the complex dynamic of the parametric family (2) on quadratic polynomials of the form p(z) = (z − a)(z − b), where a, b ∈ C. For this study, we need to present the result called the Scaling Theorem, since it allows us to conjugate the dynamical behavior of one operator with the behavior associated with another, conjugated through an affine application, that is, our operator has the same stability on all quadratic polynomials. This result will be of great use to us since we can apply the Möbius transformation on the operator R p,γ associated with our parametric family acting on p(z), assuming that the conclusions obtained will be of general application for any quadratic polynomial used.
Theorem 2 (Scaling Theorem for family (2)). Let f (z) be an analytic function in the Riemann sphereĈ and let T(z) = αz + β be an affine transformation with α = 0. We consider g(z) = λ( f • T)(z), λ = 0. Let R f ,γ and R g,γ be the fixed point operators of the family (2) associated with the functions f and g, respectively, that is to say, where y = z − f (z) f (z) and z ∈ C. Then, R f ,γ is analytically conjugated to R g,γ through T, that is to say, Therefore, substituting these equalities and simplifying, we have , that is to say, R f ,γ and R g,γ are analytically conjugated by T(z). Now, we can apply the Möbius transformation on the operator associated with the parametric family (2) in order to obtain an operator that does not depend on the constants a and b and, thus, be able to study the dynamical behavior of this family for any quadratic polynomial. The Möbius transformation, in this case, is h(z) = z−a z−b and has the following properties: The fixed-point rational operator of family (2) on p(z) has the expression We can also deduce from (23) that the order of the methods for quadratic polynomials is 3 when γ = 0 and 4 when γ = 0.

Fixed Points
The orbit of a point z ∈ C is defined (see, for example, [14,15]) as the set of the successive applications of the rational operator, i.e., The performance of the orbit of z is deduced attending to its asymptotic behavior.
Therefore, a fixed point is one that is kept invariant by the operator O γ , that is, it is one that satisfies the equation O γ (z) = z. All the roots of the quadratic polynomial are, of course, fixed points of the O γ operator. However, it may happen that fixed points appear that do not correspond to any root; we call these points strange fixed points. These points are not desirable from a numerical point of view because when an initial estimate is taken that is in the neighborhood of a strange fixed point, there is a possibility that the numerical method will converge to it, that is, to a point that is not a solution of the equation. Strange fixed points often appear when iterative methods are analyzed and their presence can show the instability of the method.
Fixed points can be classified according to the behavior of the derivative operator on them; thus, a fixed point z * can be: Moreover, the basin of attraction A(z * ) of an attracting fixed point z * is the set of initial guesses whose orbits tend to z * . Therefore, the set of points whose orbit tends to an attracting fixed point defines the Fatou set F (O γ ), while its complement is the Julia set In what follows, we study what are the fixed points of operator O γ and their character depending on the value of parameter γ. The proof of the following result is straightforward, as it only needs to solve the equation O γ (z) = z. (iii) the roots of polynomial which we denote by Ex i (γ), where i = 1, 2, . . . , 6, are also strange fixed points for each value of γ.
We need the expression of the differentiated operator to analyze the stability of the fixed points and to obtain the critical points: O γ (z) = z 2 (z + 1) 4 γ 6z 6 + 8z 5 + 7z 4 + 7z 2 + 8z + 6 + z 16z 4 + 41z 3 + 60z 2 + 41z + 16 It is clear that 0 and ∞ are always superattracting fixed points because they come from the roots of the polynomial, and the order of the iterative methods is higher than 2, but the stability of the other fixed points can change depending on the values of parameter γ. Proof. We obtain that |O γ (1)| = 96 7γ + 29 .
It is not difficult to check that |O γ (1)| cannot be 0, so z = 1 cannot be a superattractor, and, when γ = − 29 7 , z = 1 is not a fixed point. Now, we are going to study when z = 1 is an attracting point. It is easy to check that |O γ (1)| < 1 is equivalent to 96 2 < |29 + 7γ| 2 . Rewriting the last expression, we obtain the following inequality: Let us see when this inequality is verified. , we need Im(γ) to satisfy 8375 < 406Re(γ) + 49Re(γ) 2 + 49Im(γ) 2 , for z = 1 being a superattractor. We are going to study when z = 1 is a parabolic point. z = 1 will be a parabolic point when 8375 − 406Re(γ) − 49Re(γ) 2 = 49Im(γ) 2 , that is, z = 1 is a parabolic point when Re(γ) ∈ − 125 7 , 67 7 and 49Im(γ) 2 = −Re(γ) 2 − 406Re(γ) + 8375. Now, we establish the stability of the strange fixed points that are roots of the polynomial (24). To do this, we calculate these roots noting that this polynomial is a sixth degree symmetric polynomial, that is, it is a polynomial that can be reduced to a third degree one, and that satisfies the following properties: • t = 0 is not the root; • if t = α is the root, t = 1 α is also the root. Performing the reduction of (24), we obtain: To calculate the roots of polynomial (24) from the z i (γ), i = 1, 2, 3, we undo the variable change since t = z i (γ) ± z i (γ) 2 − 4 2 . Therefore, we obtain the roots of the sixth degree polynomial, which are conjugated two by two Now, we study when the roots of the polynomial (24) are superattractors. For them, we solve |O γ (Ex i (γ))| = 0 for all i = 1, . . . , 6, and we get the following relevant values of γ: Next, we are going to study the character of the fixed points by analyzing those values of γ close to the values of the parameter for which some Ex i (γ) is a supertractor. To do this, we study how |O γ (Ex i (γ))| behaves near the four previous values, and we obtain regions where some of the roots will be attractors. These regions are represented in Figure 1.
Let us remark that z = −1 is a preimage of the fixed point z = 1. We can see that q(t) is a symmetric polynomial, so we can obtain the roots of q(t) obtaining roots of a polynomial of degree 3. The polynomial reduced of q(t) is the following one that we obtain analogously to the polynomial (24): In order to calculate the roots z of q(t), we need to obtain the roots ofq(t) and apply the following expression to them z ± √ z 2 − 4 2 . Thus, the roots of q(t) are conjugated. Now, we are going to study the asymptotic behavior of the critical points to establish if there are different convergence basins than those generated by the roots. For the free critical point −1, we have O γ (−1) = 1, who is a strange fixed point, so the parameter plane associated with this critical point is not significative, since we know the stability of z = 1.
The other free critical points are roots of a polynomial that depends on γ; for that, we draw the parameter planes. As we have that the roots are conjugated, we will only draw three planes. We use as an initial estimate a free critical point that depends on γ. We establish a mesh in the complex plane of 500 × 500 points. Each point of the mesh corresponds to a parameter value. In each of them, the rational function is iterated to obtain the orbit of the critical point as a function of γ. If that orbit converges to z = 0 or to z = ∞ in less than 40 iterations, that point of the mesh is painted red; otherwise, the point appears in black.
As we can see, there are many values of the parameter γ that would result in a method in which the free critical points converge to one of the two roots. As it is observed in Figure 2, they are located in the red area on the right side of the plane. Moreover, some black areas can be identified as the regions of stability of those fixed points that can be attracting, such as Figure 1b, whose stability region appears in black on the right side of Figure 2c. Now, we select some stable (in red in parameter planes) and unstable values of γ (in black) in order to show their performance.
In the case of dynamical planes, the value of the parameter γ is fixed. Each point in the complex plane is considered as a starting point of the iterative scheme, and it is painted in different colors depending on the point that it has converged to. In this case, we paint in blue points what converged to ∞, and in orange points what converged to 0. These dynamical planes have been generated with a mesh of 500 × 500 points and a maximum of 40 iterations per point. We mark strange fixed points with white circles, the fixed point z = 0 with a white star, and free critical points with white squares (again, the routines used appear in [6]).
One value of the parameter that would be an interesting value is γ = 0 because it is the only one that obtains order 4. In that case, we obtain the dynamical plane that we can see in Figure 3a. In this case, two free critical points are in each basin of attraction, and the strange fixed points are in the boundary of both basins of attraction, so they are repulsive. In that case, the method is stable, and, as we can see, almost every point converges to 0 or ∞ (Let us notice that, in practice, any initial estimation taken in the Julia set will converge to 0 or to∞, due to the rounding error).
Other value for the parameter that we study is γ = 1, Figure 3b. As we can see, this dynamical plane is similar to that of γ = 0, but, in this case, we obtain less free critical points and less strange fixed points, due to the simplification of the rational function for this value of γ.   Carrying out numerous experiments, we have realized that the simplest dynamics is that of the methods with parameter γ = 0 and γ = 1. Next, we will see other dynamical planes associated with other values of the parameter γ. Some of these planes do not have a bad dynamics, although it is not as simple as the previous ones. This is the case of γ = 2, Figure 4b, or the case of γ = 2i, Figure 4a.
However, values such as γ = −10 + i, γ = −5 or γ = − 29 7 present a dynamical plane with the same number of basins of attraction but with more complicated performance. We can see some of these dynamical planes in Figures 5a,b and 6a. There are also parameter values for which the number of basins of attraction increases, for example, γ = 5 (Figure 6b). These cases should be avoided since our method may not converge to the roots and may end up converging to other points.

Numerical Experiments
In this section, we compare different iterative methods of the parametric family (2), solving two classical problems of applied mathematics: the Hammerstein integral equation and the Fisher partial derivative equation. We are going to use elements for the proposed class for which we have studied the dynamical plane because we want to verify that, although some of them have complicated dynamics, they can be methods that give good numerical results.
For the computational calculations, Matlab R2020b with variable precision arithmetics with 1000 digits of mantissa is used. From an initial estimation x (0) , the different algorithms calculate iterations until the stoping criterium x (k+1) − x (k) < tol is satisfied.
For the different examples and algorithms, we compare the approximation obtained, the norm of the function in the last iterate, the norm of the distance between the last two approximations, the number of iterations needed to satisfy the required tolerance, the computational time and the approximate computational convergence order (ACOC), defined by Cordero and Torregrosa in [16], which has the following expression:

Hammerstein Equation
In this example, we consider the well-known Hammerstein integral equation (see [5]), which is given as follows: where x ∈ C[0, 1], s, t ∈ [0, 1] and the kernel F is We transform the above equation into a finite-dimensional nonlinear problem by using the Gauss-Legendre quadrature formula given as , where the nodes t j and the weights ω j are determined for n = 7 by the Gauss-Legendre quadrature formula. In this case, the nodes and the weights are in Table 1. By denoting the approximations of x(t i ) by x i (i = 1, . . . , 7), one gets the system of nonlinear equations: where i = 1, . . . , 7 and Starting from an initial approximation x (0) = (−1, . . . , −1) T and with a tolerance of tol = 10 −15 , we run the parametric family for different values of the parameter γ. The numerical results are shown in Table 2. In the case of the Hammerstein integral equation, we see that the numerical results of the parametric family (2) for different values of γ are quite similar. The main difference observed between the methods is that the ACOC for γ = 0 is 4, and, for the rest of the methods, it is approximately 3. On the other hand, we note that the method with γ = −10 + i needs to perform a larger number of iterations than the rest of the methods to satisfy the required tolerance, so the time it takes to approximate the solution is also longer. Finally, taking into account the columns that measure the error of the approximation, that is, F(x (k+1) ) 2 and x (k+1) − x (k) 2 , we see that iterative methods that get lower errors are those associated with the parameters γ = 0 and γ = 2. These results confirm the information obtained in the dynamical section.

Fisher Equation
In this second example, we are going to study the equation proposed in [17] by Fisher to model the diffusion process in population dynamics. The analytical expression of this partial derivative equation is as follows: where D ≤ 0 is the diffusion constant, r is the level of growth of the species, and p is the carrying capacity.
In this case, we will study the Fisher equation for the values p = 1, r = 1, and D = 1 in the spatial interval [0, 1] and with the initial condition u(x, 0) = sech 2 (πx) and null boundary conditions.
We transform the problem we just described in a set of nonlinear systems by applying an implicit method of finite differences, providing the estimated solution in the instant t k from the estimated one in t k−1 . We denote the spatial step by h = 1 n x and the temporal step by k = T max n t , where T max is the final instant and n x and n t are the number of subintervals in x and t, respectively. Therefore, we define a mesh of the domain [0, 1] × [0, T max ], composed of points (x i , t j ), as follows: x i = 0 + ih, i = 0, . . . , n x , t j = 0 + jk, j = 0, . . . , n t .
Our objective is to approximate the solution of problem (26) in these points of the mesh, solving as many nonlinear systems as there are temporary nodes t j in the mesh. For this, we use the following finite differences: We observe that, for the time step, we use first order backward divided differences and for the spatial step they are second order centered divided differences.
By denoting u i,j as the approximation of the solution at (x i , t j ), and, by replacing it in the Cauchy problem, we get the system for i = 1, 2, . . . , n x − 1 and j = 1, 2, . . . , n t . The unknowns of this system are u 1,j , u 2,j , . . . , u n x −1,j , that is, the approximations of the solution in each spatial node for the fixed instant t j .
In this example, we are going to work with the parameters T max = 10, n x = 10 and n t = 50. As we have said, it is necessary to solve as many systems as there are temporary nodes t j ; for each of these systems, we use the parametric family (2) to approximate its solution. Thus, starting from an initial approximation u i,0 = sech 2 (πx i ), i = 0, . . . , n x , with a tolerance of 10 −6 , we execute the parametric family for different values of γ so that we get Table 3. In all cases, we obtain as an approximation of the solution of problem (26) the following vector x (k+1) = (0,4.32639,0.708718,0.853425,0.918847,0.93729,0.918847,0.853425,0.708718, 0.432639,0) T .
In this case, it can seen that the results are very similar, although there are subtle differences. For example, the method when γ = 0 uses a smaller number of iterations than the rest to satisfy the required tolerance, although this does not make it much faster than the rest of the methods since the difference in time is seconds. On the other hand, if we look at the time column, we can see that there is a method that stands out for its slowness; this is the case of γ = −10 + i. Again, we note that the ACOC of the methods roughly match the theoretical predictions made throughout the article. Observing the columns of the errors, we find similar results as well and that, in this case, having a higher tolerance than in the first example, no great differences are observed in these results.

Conclusions
A parametric family of iterative methods for solving nonlinear systems is presented. The dynamical analysis of the class on quadratic polynomials is done in order to select the members of the family with better stability properties. We prove that there exist a wide set of real and complex values of the parameter for which the corresponding methods are stable. That is, the set of initial estimations converging to the roots is very wide.In particular, we have stated that those procedures with γ = 0, γ = 1, and γ = 2 are especially stable, although some other ones can also show similar dynamical properties. Two numerical examples related to Hammerstein's equation and Fisher's equation allow us to confirm the theoretical results corresponding to the convergence and the stability of the proposed class.