A Family of Multiple-Root Finding Iterative Methods Based on Weight Functions

: A straightforward family of one-point multiple-root iterative methods is introduced. The family is generated using the technique of weight functions. The order of convergence of the family is determined in its convergence analysis, which shows the constraints that the weight function must satisfy to achieve order three. In this sense, a family of iterative methods can be obtained with a suitable design of the weight function. That is, an iterative algorithm that depends on one or more parameters is designed. This family of iterative methods, starting with proper initial estimations, generates a sequence of approximations to the solution of a problem. A dynamical analysis is also included in the manuscript to study the long-term behavior of the family depending on the parameter value and the initial guess considered. This analysis reveals the good properties of the family for a wide range of values of the parameter. In addition, a numerical test on academic and engineering multiple-root functions is performed.


Introduction
Solving nonlinear equation f (x) = 0 has been a difficult problem to handle for a long time in several branches of Science and Engineering. Although Newton's method dates from the 17th century, it is at the end of the 20th century that computing drives research in the discipline of iterative methods. Hundreds of papers can be found covering this topic in the most recent literature.
An initial classification splits the problem in methods without or with memory for single or multiple roots. The variety of publications is vast. Furthermore, the stability of the iterative procedures is also analyzed in several publications. Focusing on methods without memory, we can find complete analysis for single-root finding on simple two-step methods [1][2][3][4][5], as well as those that reach orders of convergence of value 16 [6][7][8][9] most of them collected in [10]. In a similar manner, multiple-root finding procedures-recently very aptly named as root-ratio methods [11]-are also present in the literature with a large variety of orders of convergence (see [12,13] and the references therein). In this sense, the higher order of convergence, the higher computational cost.
Iterative methods for finding single roots of a nonlinear equation often lose their order of convergence when used for multiple roots. A classic example is Newton's method, with quadratic order of convergence for single roots, but it has linear convergence when it is applied for finding multiple roots. However, several modifications have been proposed in order to keep its quadratic order. Among them, we find Newton's accelerated method (or Rall's method [14]) and the modified Newton's scheme (or Schröder's method [15]), based on applying Newton's iterative scheme to the functions u 1 (x) = f (x) 1/m and u 2 (x) = f (x) f (x) , respectively. In both cases, the problem becomes finding a simple root of a nonlinear equation, so the order of convergence is again quadratic.
Linear combinations of these two previous methods can also generate schemes that hold the order of convergence two. Mora [16] proposes a one-parametric family. In particular, the method of the family obtained when α = 1 is the well-known SuperHalley's method that has order of convergence three.
In this paper we focus on a simple idea for generating a class of one-step iterative methods for finding multiple roots of nonlinear equations. A specific case of weight function is developed, including a parameter, and therefore a family of iterative methods is generated. In order to select the most stable methods of the family, a dynamical study is carried out. To this end, we perform a complex dynamics analysis, in which fixed and critical points-and their asymptotic behavior-allows the knowledge of the stability of the family. As a consequence, we will be able to select the members of the family whose basins of attraction are wider, ensuring the convergence of the method for a wide set of initial estimates. This analysis is more common for single root methods; however, there are papers that include it for multiple root methods [17][18][19][20][21].
This paper is organized as follows. In Section 2 we present the iterative one-step family based on a weight function, as well as its convergence analysis. The stability analysis of the family is performed in Section 3, where we find that compliance with the Scaling Theorem simplifies the study. Section 4 covers the numerical analysis for specific members of the family on different nonlinear functions with multiple roots. Finally, Section 5 collects the main conclusions.

Convergence Analysis of the Parametric Family
Mora [16] introduced the one-parameter family Family (1) can also be written as Following this guideline and in order to design iterative schemes of order three for multiple roots, we propose the following triparametric family where the weight function and its variable are given by Theorem 1 shows the conditions that the parameters b 1 and b 2 must satisfy for the family to reach the order of convergence three. The proof of the theorem is based on Taylor series developments around the solution since the initial estimates are taken close to the solution. The only requirement for the function is to be sufficiently differentiable. Theorem 1. Let us consider r ∈ I a multiple root with multiplicity m > 1 of a sufficiently differentiable function f : I ⊂ R −→ R defined in an open interval I. If the initial estimation x 0 is close enough to r, then iterative family (3) has order of convergence 3 when parameters satisfy In addition, its error equation is given by , i = 1, 2, . . ., and e k = x k − r denotes the error in each iteration.
Proof. Considering the Taylor series expansions of f (x k ), f (x k ) and f (x k ) around r and taking into account that f (j) (r) = 0, j = 0, 1, . . . , m − 1 and f (m) (r) = 0, we have being

Using Equation (5) and developing the weight function variable
where the coefficients K i , i = 1, 2, 3, are given by .
In order to achieve order of convergence 3, we must solve equations K 1 = 0 and K 2 = 0. Then, we obtain the values for parameters b 1 and b 2 : Finally, by using parameters Equation (7) the error Equation (6) turns into so the order of convergence is three.
If we set b 3 = α and parameters b 1 and b 2 from (4), the weight function of family (3) turns into After some algebraic manipulations, (8) can be written as a Chebyshev-Halley type expression: Let us note that m = 1 gives Chebyshev-Halley's family of iterative schemes for simple roots .
From now on, we consider the third-order iterative family resulting after applying the previous conditions. This parametric family can be expressed as follows: Family Equation (10) collects all Chebyshev-Halley's methods for single or multiple roots and only these. Some of the well-known methods of this family can be deduced from different values of α: 2. Halley (α = 1 2 ): 3. Super-Halley (α = 1): 4. Osada (α = ∞):

Dynamical Analysis of the Family
In this section, we conduct a dynamical analysis of the iterative family under study trough complex dynamics. After a review of the basics on complex dynamics, we apply the scaling theorem. The rational function after Möbius conjugation gives us an expression from which the fixed and critical points are obtained. Finally, some well-known representations-such as the stability, the parameters or the dynamical plane-provides us with important information about the stability of the family.
First, the fundamental tools needed to begin this complex dynamics study are described in Section 3.1. This dynamical tools are defined for a rational function, which is obtained later after applying the corresponding iterative method on a polynomial.

Preliminaries on Complex Dynamics
In this part, the basics of the complex dynamics are revisited. For a more extensive explanation, see [22][23][24].
Let Q :Ĉ →Ĉ be a rational function, whereĈ is the Riemann sphere. The orbit of a point z 0 ∈Ĉ is the set consisting of the successive applications of Q on the point z 0 , i.e., A point that remains unaltered with the application of Q is a fixed point, that satisfies Q(z f ) = z f . The asymptotic behavior of the fixed points sorts them into attracting, repelling, neutral or superattracting, depending on the value of the multiplier |Q (z f )|. Table 1 gathers the conditions for classifying a fixed point inside one of the above-mentioned groups. Table 1. Classification of fixed points.
Asymptotic behavior Attracting Repelling Neutral Superattracting The evaluation of the infinity as a fixed point and its asymptotic behavior can be obtained with the technique presented in [25] and widely applied ( [26] or [27], among others).
When the rational function is the result of the application of a function f (z) on an iterative method, the fixed points match with the roots of f . In addition, more fixed points may appear, which are called strange fixed points. The presence of these additional points may alter the stability of the iterative method, as long as its asymptotic behavior is not repelling.
A critical point z c is a point that satisfies Q (z c ) = 0. Again, the roots of f match with the critical points, and the additional points that do not match are named free critical points.
Finally, the basin of attraction of an attracting fixed point is the set of initial estimations whose orbit converge to this attracting fixed point, i.e., When the iterative expression includes a parameter, it is worth mentioning the associated graphic tools [28]. On the one hand, the unified stability plane gathers in one representation the regions of the parameter where the strange fixed points behave as attracting points, altering the stability of the method. On the other hand, the unified parameter plane represents the regions of the parameter where the free critical points converge either to a root of the function f or to another point.

The Rational Function
In the following, we are going to analyze the stability of family (10) on cubic polynomials with double roots. We will use for this study the generic polynomial Previously, we will apply the Scaling Theorem in order to simplify this study.
For this purpose, let us recall the fixed point operator obtained when family (10) is applied to an analytical function f (z) inĈ: where Theorem 2 (Scaling Theorem for family (10)). Let f (z) be an analytic function inĈ and Q f (z) the fixed point operator obtained when family (10) is applied to f (z). Let us consider an affine map A(z) = βz + γ with , λ = 0, and Q g (z) the corresponding fixed point operator. Then: , so Q f and Q g are affine conjugated by A.
Proof. First, let us consider the following equalities: Then we have (16) and From (16) and (17), we obtain The Scaling Theorem (Theorem 2) proves that the fixed point operators associated to family (10) applied to different analytic functions are affine conjugated by an affine map. Thus, the dynamical study of the roots of a cubic polynomial with a double root can be reduced by an affine change of coordinates to the dynamics associated with the generic cubic polynomial p(z) = (z − a) 2 (z − b).
The rational operator obtained when family (10) is applied on p(z) is which depends not only on α but also on parameters a and b.
In order to work with a rational operator with dependence only on a parameter, we apply Möbius transformation and we obtain a rational operator that only depends on parameter α Fixed point operator Op(z, α) is affine conjugated to Q p (z, α, a, b) by M, so the study of the dynamics is equivalent.
Since M(a) = 0, M(b) = ∞ and M(∞) = 1, Möbius transformation has sent the roots a and b to 0 and ∞, respectively, while the divergence is set in 1. Now the fixed points of Op(z, α) are 0 and ∞, which correspond to the roots of p(z) with multiplicity 2 and 1, respectively.
In order to study the asymptotic behavior of the fixed points and to obtain the critical points of the rational operator, we consider its derivative, whose expression is However, the study of the fixed point ∞ requires a different analysis [25]. The point ∞ is a fixed point of Op(z, α) if, and only if, z = 0 is a fixed point of the function F(z, α) = 1 Op(1/z,α) . In addition, ∞ is classified as attracting when |F (0, α)| < 1, repelling if |F (0, α)| > 1 and neutral when |F (0, α)| = 1. In this analysis, rational functions F(z, α) and F (z, α) are given by For all value of α it is satisfied Op(0, α) = 0 and F(0, α) = 0, so the roots z = 0 and z = ∞ are fixed points of the fixed point operator. Now, we analyze the set of strange fixed points and free critical points of Op(z, α) and also the asymptotic behavior of the fixed points. Lemma 1. The fixed point operator Op(z, α) has the following three strange fixed points when α ∈ C − {−2, 2, 2/3, 14/9}.

Proof.
Solving the equation Op(z, α) = z, we obtain the strange fixed points z e 0 , z e 1 and z e 2 . On the one hand, the values α = ±2 are roots of −6 + 5α − √ 28 − 60α + 27α 2 , so they cancel the numerator of z e 1 and z e 2 . On the other hand, α = 2/3 and α = 14/9 are roots of 28 − 60α + 27α 2 . In addition, α = 14/9 is solution of the equations z e 1 = 1 and z e 2 = 1. Then, we consider this set of values of α and the associated fixed point operator is analyzed for each case.
Op(z, 2) = z 2 2 − z and z = 1 is the only strange fixed point.
• If α = 14/9, the fixed point operator is The only strange fixed point for this case is z = 1.
From (20), solving equation Op (z, α) = 0 we obtain the critical points as shown the following result.

Lemma 2.
The critical points of the fixed point operator Op(z, α) are z = 0, corresponding to the double root of the polynomial, and four free critical points, corresponding to the roots of polynomial In particular, when α = 2 there is only a free critical point, for α ∈ {∞, 2/3, 3/2, 4/3} there exist two different free critical points, and if α ∈ {0, 10/9, −5.66328, 1.59444, 1.58399}, the fixed point operator has three free critical points.

Asymptotic Behavior of Fixed and Critical Points
In order to design iterative schemes with stability depending on the initial estimations, it is desirable that the strange fixed points are repelling, so the iterative process does not converge to points different to the roots of the function.
According to Lemma 1, rational operator Op(z, α) has strange fixed points, so we must study its asymptotic behavior to determine the values of parameter α that correspond to the methods of the family with the best dynamic properties. Proof. On the one hand, from (20) we have Op (0, α), so z = 0 is a superattracting fixed point. On the other hand, from (22), F (0, α) = − α 2 . Then, z = ∞ is: • Attracting if − α 2 < 1, that is, when |α| < 2 (and superattracting when α = 0).
• Superattracting if 6(−3 + 2α) −4 + 3α = 0, that is when α = 3 2 . The asymptotic behavior of the strange fixed points z e 1 and z e 2 is set through the derivatives However, the resulting expressions make difficult the analytical study. Instead of it, we represent the stability regions in the complex plane for each strange fixed point depending on the value of α.  In addition, Figure 2 shows the value of α in the complex plane where the strange fixed points and z = ∞ behave as attracting points. Let us remark that Figure 3b is a detail of the red region represented in Figure 3a where the root z = ∞ is attracting. According to Figure 3a, the only attracting point outside the green and red regions is the double root z = 0. In order to choose the methods of the family with the best properties in terms of stability, we must avoid the areas represented with any color other than white in Figure 3b.
In a family of iterative schemes with dependence on a parameter, the parameters plane provides an overview of the stability of the family. In this paper, parameters planes are generated using a mesh of 800 × 800 points in the complex plane, where the real and imaginary part of parameter α are represented in the abscissas and ordinate axis, respectively. Each point in the plane corresponds to a value of α, that is, with a method of the iterative family. Taking a free critical point as initial estimation to iterate each method associated to a value of α, the point is represented in red when there is convergence to any root of the polynomial and otherwise, the point is plotted in black. Therefore, there exist points in the black regions, different to the roots 0 and ∞, that behave as attracting points. These points can be strange fixed points or attracting periodic points. Then, the most stable methods of the family are located in the plane in the red regions.
According to Lemma 2, rational operator Op(z, α) has four free critical points. We can see in Figure 4 the parameters plane associated to each free critical point for ( (α), (α)) ∈ [−3, 3] × [−3, 3]. They have been generated using Matlab and the implementation presented in [29]. The iterative process finishes when the difference between the point and some of the roots is lower than 10 −3 with a maximum of 150 iterations. Figure 4 shows wide red areas for the four free critical points. This areas allow us to select values of parameter α where the associated scheme converges to a root taking as initial iterate a free critical point. In addition, we have choose from Figure 4 different values of α to generate the corresponding dynamical planes and then compare the stability of methods that belong to the same family. The dynamical planes show the basins of attraction of an attracting point. Each point in the complex plane is taken as initial estimation to iterate the method and, according to the root to which its orbit converges, it is represented with a color. In this paper, the dynamical planes have been generated in Matlab according to [29]. A mesh of 800 × 800 points in the complex plane are taken as initial estimations, being the stopping criteria a maximum of 50 iterations or a difference between the point and an attracting point lower than 10 −3 .
We have denoted the basins of attraction of z = 0 and z = ∞ in orange or green, respectively. Moreover, colors blue and red denote the convergence to a strange fixed point. Finally, the point in the plane is depicted in black when its orbit is divergent or converges to any other attracting point different to the roots or the strange fixed points. Moreover, colors are lighter when less iterations are required until the convergence. Figures 5 and 6 show the dynamical planes associated to values of α where the number of strange fixed points is lower than three (Lemma 1), that is α ∈ {−2, 2/3, 14/9, 2}. In all the cases the strange fixed points are repelling, so we do not observe red or blue regions. In particular, for α = 2 there is only convergence to the double root z = 0 or to the periodic orbit {−20001, 19999}.  In Figure 7 we can observe the dynamical planes corresponding to Chebyshev (α = 0), Halley (α = 1 2 ), SuperHalley (α = 1) and Osada's method (α = ∞) for multiple roots taking m = 2. Although each of these methods has three strange fixed points, they are not attracting. As for the first three methods α < 2, they can converge to the single root. Moreover, Osada's method has a wide black area of convergence to the attracting periodic orbit {−6.5102, 2.9057}. Finally, we show in Figure 8 the dynamical planes associated with values of the parameter with attracting strange fixed points. For α = 3/2, the strange fixed point z e 0 = 1 is superattracting, while if α = 1.58, z e 1 = 0.6277 and z e 2 = 1.4952 are attracting, so we can observe their basins of attraction represented in red and blue. Both values of α correspond to points represented in black in the parameters plane, so they are unstable methods.

Numerical Test
In order to analyze the behavior of the family for solving multiple-root nonlinear equations, this section performs a numerical test on different problems.
The computations have been performed using MATLAB R2017b. We use of variable precision arithmetic with 50 digits of mantissa. We consider that the solution has been obtained when the absolute difference of two iterates |x k+1 − x k | is lower than 10 −30 . If the iterative method needs more than 50 iterations to reach the solution, we consider that the procedure does not converge (N.C.).
Tables 2-4 collect the results of the numerical test on different nonlinear functions with multiple roots. They display the method, the initial guess (x 0 ), the amount of iterates to reach the multiple root (# iter), the difference between the two last iterates (|x k+1 − x k |) and the approximated computational order of convergence (ACOC) [30].
The first nonlinear function under study is Let us note that this function has a double real root for x = 2 and a single root for x = −3. Table 2 shows the results of the iterative methods for two different initial estimations for reaching the double real root.
The second nonlinear function under study is that has two double real roots: x = 1 and x = 12.3402. Table 3 collects the values resulting from the application of the iterative methods to solve f 2 (x) = 0. The third nonlinear function is , that has a root x = 1.228 with multiplicity 9. Table 4 displays the resulting results.

Conclusions
A new family of iterative methods for solving multiple-root nonlinear equations based on weight functions with parameters has been introduced.
The introduced family satisfies the scaling theorem. Therefore, the stability analysis is simplified. A complete dynamical analysis of the family -including the analysis of the fixed and critical points, or the parameters and stability planes -has been performed. This study allows the selection of the best members of the family in terms of stability, with the help of the stability and the parameters planes. For a wide set of values of α, we can guarantee the convergence to the desired root, as the dynamical planes show.
Selected members of the family have been tested in a numerical comparison with other well-known methods. The use of nonlinear functions with multiple roots evidences in the numerical results that the proposed family is competitive with the rest of the methods.