Numerical Properties of Different Root-Finding Algorithms Obtained for Approximating Continuous Newton’s Method

: This paper is dedicated to the study of continuous Newton’s method, which is a generic differential equation whose associated ﬂow tends to the zeros of a given polynomial. Firstly, we analyze some numerical features related to the root-ﬁnding methods obtained after applying different numerical methods for solving initial value problems. The relationship between the step size and the order of convergence is particularly considered. We have analyzed both the cases of a constant and non-constant step size in the procedure of integration. We show that working with a non-constant step, the well-known Chebyshev-Halley family of iterative methods for solving nonlinear scalar equations is obtained.


Introduction
We can find the origins of continuous Newton's method in the seminal paper of Neuberger [1].Actually, it appears in the study of the basins of attraction related with the relaxed Newton's method for solving a complex equation p(z) = 0 Equation ( 1)) have an intricate fractal structure.Neuberger, as well other authors ( [2,3]), has realized that this fractal structure shrinks away when h → 0, as we can appreciate in Figure 1.
We can identify iterative method Equation (1) as an Euler approximation of the differential equation with step size h.Initial value problem Equation (2), or its improved version of finding a function z : [0, ∞) → C such that z(0) = z 0 , p(z) (t) = −p(z(t)) for a given z 0 ∈ C, where p is a non-constant complex polynomial is known as continuous Newton's method.We refer to [1] for the theoretical basis of continuous Newton's method.In particular, it is shown that solutions z(t) of Equation (2) (or Equation ( 3)) flow to a zero of p while keeping the argument of p(z(t)) constant at arg(p(z 0 )).
For instance, as it was pointed out by Jacobsen et al.
The choice of the appropriate branch of the cube root, defined by the rays θ = π/3, θ = π and θ = −π/3, determines the ternary division of the complex plane that can be seen in the three graphics of Figure 1.According to these authors, the fractal boundaries in the basins of attraction of the roots can be originated by numerical errors inherent to the discretization of initial value problem Equation (2).
This work is two fold.In Section 2, we consider other different strategies (not only Euler's method) for solving initial value problem Equation (2).The efficiency of the iterative processes obtained in this way is analyzed.

Numerical Algorithms Applied to Continuous Newton's Method
In [3] the dynamical properties of six numerical methods for solving differential equations are considered when they are applied to continuous Newton's method Equation (2).In particular, they compare the basins of attraction of these methods applied to Equation (2) for p(z) = z 3 − 1.They estimate the corresponding fractal dimension using a box-counting algorithm.The main conclusion is that higher-order algorithms are not necessarily related with basin boundaries with smaller fractal dimension.
In this section we are concerned with some numerical properties of root-finding methods derived from the application of numerical methods for solving differential equations as continuous Newton's method Equation (2).So, as it has been done in [3], we use y n for the approximations of z(t n ), where t n+1 = t n + h and h for the step size.The six methods considered by Jacobsen et al. in [3] are: 1. Euler's method: 4. Runge-Kutta method of order 2: 5. Runge-Kutta method of order 4: 6. Adams-Bashforth method of order 2: Note that, in the case of continuous Newton's method, the function f (t, z(t)) that appears in differential equation Equation ( 4) is So, the application of Euler's method Equation ( 5) gives rise to the root-finding algorithm that is the wel-known relaxed Newton's method.Note that the role of the step size h in the methods for solving differential equations is moved to a relaxing parameter in Equation ( 12).As it was stated in [3] and we have previously mentioned, h has a clear influence in the dynamical properties of an iterative method.We analyze now the influence of h in the numerical properties of an iterative method.Let F 1 be the iteration map related to the relaxed Newton's method Equation ( 12) and let z * be a simple root of p(z) = 0.It is a straightforward calculation to show that So, we have that method Equation (12) is consistent (the roots of p are attracting fixed points of the iteration map F 1 ) only for h ∈ (0, 2).In addition, we obtain iterative methods with just linear convergence, except for the case h = 1, that is the classical Newton's method with quadratic convergence.Now, we see what happens with the refined Euler's method defined in Equation ( 6).The corresponding root-finding algorithm can be written as As a consequence, iterative method Equation ( 12) is consistent only for h ∈ (0, 2).All the methods deduced in this case have linear convergence (the lower value for the asymptotic constant error is obtained for h = 1).This fact together with a number of function evaluations equals to four, makes this method uninteresting from an efficiency point of view.
Something similar happens with Heun's method given in Equation ( 7) and Runge-Kutta method of second order given in Equation (8).They have different iteration maps, respectively z n+1 = F 3 (z n ) and z n+1 = F 4 (z n ), where p(z) p (z) In both cases, we have for a simple root z * of p(z) = 0. So, as in the case of the refined Euler's method, only linear convergence is achieved for h ∈ (0, 2) and then these two methods are inefficient from a numerical perspective.
The study of the Runge-Kutta method of order fourth defined in Equation ( 9) applied to Equation (4) leads us to the iteration scheme z n+1 = F 5 (z n ), where Note that, for a simple root z * , we have There are no real values of h such that F 5 (z * ) = 0, then only linear convergence can be attained.In addition, we obtain consistent methods, that is z * is an attracting fixed point of F 5 if |F 5 (z * )| < 1.These inequalities happens for h ∈ (0, 2.7853) where 2.7853 is, approximately, the only real root of the polynomial −24 + 12h − 4h 2 + h 3 .The optimum value for h, for which the asymptotic error constant has the lowest value, is h * = 1.5961.In this case F 5 (1.5961) = 0.2704.
The analysis of the Adams-Bahforth second order method given in Equation (10) leads us to the following two-step multipoint method: To study the local order of convergence related to Equation ( 14), we consider a simple root z * of p(z) = 0 and the error in the n-th step e n = z n − z * .Then, where the terms of order higher than two in the error have been neglected.This approached formula for the errors generates a linear recurrence of second order, whose characteristic equation is Consequently, the method Equation ( 14) is consistent if the two roots of the previous equation, have module less than one.Note that λ + < 1 for all h > 0 but λ − ≤ −1 for h ≥ 1.So, the consistency condition is satisfied only if h ∈ (0, 1).In this case, λ + ∈ (1/3, 1) and λ − ∈ (−1, 0) .Taking into account that λ + is a decreasing function for h > 0 and −λ − is an increasing function for h > 0, we can then obtain the best convergence rate when λ + = −λ − , that is, for h = 2/3.
In Table 1 we show, in brief, the numerical information that we have deduced for methods Equations ( 5)-( 9) considered by Jacobsen et al. in [3].In fact, for each method, we show the asymptotic error constant (A.E.C.), the interval of values of h for which consistent methods are obtained (I.C.) and the optimal value for h regarding the order of convergence, h * .In addition, the Adams-Basforth method given in Equation ( 10) is consistent for h ∈ (0, 2) with h * = 2/3.As a conclusion, we can say that, despite the order of the numerical methods for solving differential equations and the size of the step h, the root-finding algorithms derived from their application to continuous Newton's method Equation ( 2) have only a linear order of convergence, with the exception of Euler's method with h = 1, where the convergence is quadratic.In this sense, we conclude that the numerical efficiency of these methods is poor.However, the study of the aforementioned methods can be interesting from other points of view, that reveals other topological or dynamical aspects as, for instance, the impact of the stepsize on the fractal dimension of the boundaries of the basins of attraction for the associated roots.Table 1.Some numerical properties of the root-finding methods obtained after applying methods Equations ( 5)-( 9) to continuous Newton's method Equation (2).Refined Euler Equation ( 6) Heun Equation ( 7) Runge-Kutta 4 Equation ( 9) (h 4 − 4h 3 + 12h 2 − 24h + 24)/24 (0, 2.7853) 1.5961

Numerical Algorithms with Non-Constant Step Size
As we can see, the application of the numerical algorithms considered in [3] to the continuous Newton's method Equation (2) generates iterative root finding methods that does not present a high numerical efficiency.The situation does not improve if we consider higher order numerical methods for the initial value problem Equation (2).In fact, if we consider the second order Taylor's method we obtain the following iterative method where L p (z) has been defined in Equation (13).But, once again, so that the convergence of the method is just linear although the number of functional evaluations have been increased.
To avoid these difficulties, we can consider iterative methods with a non-constant step size.For instance, instead of considering method Equation (12) derived from Euler's method, we consider methods where the step h is substituted by a non-constant function H(z n ): As in the previous section, if we consider a simple root z * of p(z) = 0, we have F (z * ) = z * .In addition, So, if we consider a function H that satisfies conditions Equations ( 16) and ( 17), we can obtain iterative methods with at least cubic convergence.One choice (not the only one) for function H that fulfills Equations ( 16) and ( 17) is that leads us to the well-known Chebyshev's method ( [4,5]): Even more, if we choose we deduce the well-known Chebyshev-Halley family of methods ( [4,5]): Note that Chebyshev's iterative method defined in Equation ( 18) is one of the methods defined in Equation (19), actually for α = 0. Other methods belonging to this family are Halley's method (α = 1/2) or super-Halley's method (α = 1).The dynamical behavior of methods in the family Equation (19) have been studied in detail in the work of Cordero et al. [6].The structure and dynamical properties of the basins of attraction of the roots of a given function changes depending on the choice of the iterative method considered for solving the equation.As a visual sample we show in Figure 2 the basins of attraction of three methods in family Equation (19) when they applied to the polynomial p(z) = z 3 − 1.We can generalize family Equation ( 19) to obtain the general form of methods with cubic convergence given by Gander in [7].In fact, the iteration given by where H is a function satisfying H(0) = 1, H (0) = 1/2 and H (0) < ∞, has cubic convergence to simple roots of p.The task of extending in this way Gander's result to a higher order of convergence seems to be difficult.An approach for the quadratic case was given by Romero et al. in [8].

Conclusions
We have made a numerical incursion in some properties of the root-finding methods obtained when different numerical procedures have been applied to the initial value problem known as continuous Newton's method.In particular, we have shown that, for a constant step size, the deduced iterative root-finding methods have a poor efficiency behavior.However, if we consider numerical methods with non-constant step size, a plethora of root-finding methods can be constructed.In this work we are only concerned with the famous Chebyshev-Halley family of methods, but many other iterative methods can be built up in this way.We must take into account that just Euler's method with non-constant step size has been considered in this work.