Next Article in Journal / Special Issue
On the Kung-Traub Conjecture for Iterative Methods for Solving Quadratic Equations
Previous Article in Journal
A New Smoothing Conjugate Gradient Method for Solving Nonlinear Nonsmooth Complementarity Problems
Previous Article in Special Issue
On the Local Convergence of a Third Order Family of Iterative Processes
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

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

Department of Mathematics and Computer Sciences, University of La Rioja, Logroño 26004, Spain
Algorithms 2015, 8(4), 1210-1218; https://0-doi-org.brum.beds.ac.uk/10.3390/a8041210
Submission received: 28 October 2015 / Revised: 10 December 2015 / Accepted: 14 December 2015 / Published: 17 December 2015
(This article belongs to the Special Issue Numerical Algorithms for Solving Nonlinear Equations and Systems)

Abstract

:
This paper is dedicated to the study of continuous Newton’s method, which is a generic differential equation whose associated flow tends to the zeros of a given polynomial. Firstly, we analyze some numerical features related to the root-finding 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.

1. 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
z n + 1 = z n h p ( z n ) p ( z n ) , n 0 , h C
It is well known, since the works of Cayley and Schröder at the end of the 19th century, and Fatou and Julia in the early 20th century, that the basins of attraction of Newton’s method (obtained for h = 1 in 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.
Figure 1. Basins of attraction of the relaxed Newton’s method given in Equation (1) applied to the polynomial p ( z ) = z 3 1 for h = 1 , h = 2 / 3 and h = 1 / 3 , respectively.
Figure 1. Basins of attraction of the relaxed Newton’s method given in Equation (1) applied to the polynomial p ( z ) = z 3 1 for h = 1 , h = 2 / 3 and h = 1 / 3 , respectively.
Algorithms 08 01210 g001
We can identify iterative method Equation (1) as an Euler approximation of the differential equation
z ( 0 ) = z 0 , z ( t ) = p ( z ( t ) ) p ( z ( t ) )
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 z0 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. [3], if p ( z ) = z 3 1 in Equation (2), the explicit solution is
z ( t ) = ( z 0 3 1 ) e t + 1 3
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.

2. 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
d z d t = f ( t , z ( t ) ) z ( t 0 ) = z 0
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:
  • Euler’s method:
    y n + 1 = y n + h f ( t n , y n )
  • Refined Euler’s method:
    y * = y n + h / 2 f ( t n , y n ) y n + 1 = y n + h f ( t n + h / 2 , y * )
  • Heun’s method:
    y * = y n + h f ( t n , y n ) y n + 1 = y n + h / 2 f ( t n , y n ) + f ( t n + h , y * )
  • Runge-Kutta method of order 2:
    y * = y n + 2 h / 3 f ( t n , y n ) y n + 1 = y n + h / 4 f ( t n , y n ) + 3 f ( t n + 2 h / 3 , y * )
  • Runge-Kutta method of order 4:
    k 1 = h f ( t n , y n ) k 2 = h f ( t n + h / 2 , y n + k 1 / 2 ) k 3 = h f ( t n + h / 2 , y n + k 2 / 2 ) k 4 = h f ( t n + h , y n + k 3 ) y n + 1 = y n + 1 / 6 k 1 + 2 k 2 + 2 k 3 + k 4
  • Adams-Bashforth method of order 2:
    y * = y 0 + 2 h / 3 f ( t n , y 0 ) y 1 = y 0 + h / 4 f ( t n , y 0 ) + 3 f ( t n + 2 h / 3 , y * ) y n + 1 = y n + y n + h 3 / 2 f ( t n , y n ) 1 / 2 f ( t n 1 , y n 1 )
Note that, in the case of continuous Newton’s method, the function f ( t , z ( t ) ) that appears in differential equation Equation (4) is
f ( t , z ( t ) ) = p ( z ( t ) ) p ( z ( t ) )
So, the application of Euler’s method Equation (5) gives rise to the root-finding algorithm
z n + 1 = F 1 ( z n ) = z n h p ( z n ) p ( z n ) , n 0
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
F 1 ( z * ) = z * , F 1 ( z * ) = 1 h
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
z n + 1 = F 2 ( z n ) = z n h p ( ω 2 ( z n ) ) p ( ω 2 ( z n ) ) , n 0 , ω 2 ( z ) = z h 2 p ( z ) p ( z )
If z * is a simple root of p ( z ) = 0 , we have ω ( z ) = 1 h ( 1 L p ( z ) ) / 2 , where
L p ( z ) = p ( z ) p ( z ) / p ( z ) 2
Therefore, ω ( z * ) = z * , ω ( z * ) = 1 h / 2 and then
F 2 ( z * ) = z * , F 2 ( z * ) = 1 + ( 1 h ) 2 2
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
F 3 ( z ) = z h 2 p ( ω 3 ( z ) ) p ( ω 3 ( z ) ) + p ( z ) p ( z ) n 0 , ω 3 ( z ) = z h p ( z ) p ( z )
F 4 ( z ) = z h 4 p ( ω 4 ( z ) ) p ( ω 4 ( z ) ) + p ( z ) p ( z ) , n 0 , ω 4 ( z ) = z 2 h 3 p ( z ) p ( z )
In both cases, we have
F j ( z * ) = z * , F j ( z * ) = 1 + ( 1 h ) 2 2 , j = 3 , 4
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
F 5 ( z ) = z h 6 p ( z ) p ( z ) + 2 p ( z + k 1 ( z ) / 2 ) p ( z + k 1 ( z ) / 2 ) + 2 p ( z + k 2 ( z ) / 2 ) p ( z + k 2 ( z ) / 2 ) + p ( z + k 3 ( z ) ) p ( z + k 3 ( z ) )
and
k 1 ( z ) = h p ( z ) p ( z ) , k 2 ( z ) = h p ( z + k 1 ( z ) / 2 ) p ( z + k 1 ( z ) / 2 ) , k 3 ( z ) = h p ( z + k 2 ( z ) / 2 ) p ( z + k 2 ( z ) / 2 )
Note that, for a simple root z * , we have k 1 ( z * ) = 0 , k 1 ( z * ) = h , k 2 ( z * ) = 0 , k 2 ( z * ) = h ( 1 h / 2 ) , k 3 ( z * ) = 0 , k 3 ( z * ) = h ( 1 h / 2 + h 2 / 4 ) and then
F 5 ( z * ) = z * , F 5 ( z * ) = h 4 4 h 3 + 12 h 2 24 h + 24 24
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 + 12 h 4 h 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:
z n + 1 = z n h 2 3 p ( z n ) p ( z n ) p ( z n 1 ) p ( z n 1 ) , n 1
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,
e n + 1 = e n h 2 ( 3 u ( z * + e n ) u ( z * + e n 1 ) ) 1 3 h 2 e n + h 2 e n 1
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
λ 2 1 3 h 2 λ h 2
Consequently, the method Equation (14) is consistent if the two roots of the previous equation,
λ = 1 4 3 h + 2 9 h 2 4 h + 4 , λ + = 1 4 3 h + 2 + 9 h 2 4 h + 4
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).
Table 1. Some numerical properties of the root-finding methods obtained after applying methods Equations (5)–(9) to continuous Newton’s method Equation (2).
MethodA.E.C.I.C. h *
Euler Equation (5) 1 h ( 0 , 2 ) 1
Refined Euler Equation (6) ( 1 + ( 1 h ) 2 ) / 2 ( 0 , 2 ) 1
Heun Equation (7) ( 1 + ( 1 h ) 2 ) / 2 ( 0 , 2 ) 1
Runge-Kutta 2 Equation (8) ( 1 + ( 1 h ) 2 ) / 2 ( 0 , 2 ) 1
Runge-Kutta 4 Equation (9) ( h 4 4 h 3 + 12 h 2 24 h + 24 ) / 24 ( 0 , 2 . 7853 ) 1 . 5961

3. 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
z n + 1 = z n + h f ( t n , z n ) + h 2 2 f t ( t n , z n ) + f z ( t n , z n ) f ( t n , z n )
we obtain the following iterative method
z n + 1 = F 6 ( z n ) = z n h p ( z n ) p ( z n ) + h 2 2 ( 1 L p ( z n ) ) p ( z n ) p ( z n )
where L p ( z ) has been defined in Equation (13). But, once again,
F 6 ( z * ) = z * , F 6 ( z * ) = 1 + ( 1 h ) 2 2
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 ) :
z n + 1 = F ( z n ) = z n H ( z n ) p ( z n ) p ( z n ) , n 1
As in the previous section, if we consider a simple root z * of p ( z ) = 0 , we have F ( z * ) = z * . In addition,
F ( z * ) = 0 if H ( z * ) = 1
F ( z * ) = 0 if H ( z * ) = L p ( z * ) / 2
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
H ( z ) = 1 + 1 2 L p ( z )
that leads us to the well-known Chebyshev’s method ([4,5]):
z n + 1 = F 7 ( z n ) = z n 1 + 1 2 L p ( z n ) p ( z n ) p ( z n )
Even more, if we choose
H ( z ) = 1 + L p ( z ) 2 ( 1 α L p ( z ) )
we deduce the well-known Chebyshev-Halley family of methods ([4,5]):
z n + 1 = F α ( z n ) = z n 1 + L p ( z n ) 2 ( 1 α L p ( z n ) ) p ( z n ) p ( z n ) , α C
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 .
Figure 2. Basins of attraction of the Chebyshev, Halley and super-Halley methods ( α = 0 , α = 1 / 2 and α = 1 , respectively, in Equation (19)) applied to the polynomial p ( z ) = z 3 1 .
Figure 2. Basins of attraction of the Chebyshev, Halley and super-Halley methods ( α = 0 , α = 1 / 2 and α = 1 , respectively, in Equation (19)) applied to the polynomial p ( z ) = z 3 1 .
Algorithms 08 01210 g002
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
z n + 1 = z n H ( L p ( z n ) ) p ( z n ) p ( z n )
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].

4. 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.

Acknowledgments

This work has been supported by the research project Ref. MTM2014-52016-C2-1-P by the Spanish Ministry of Economy and Competitiveness. The graphics showed in Figure 1 and Figure 2 have been made in Mathematica following the patterns given by Varona [9]. The author thanks the referees for their helpful suggestions.

Conflicts of Interest

The author declares no conflict of interest.

References

  1. Neuberger, J.W. Continuous Newton’s method for polynomials. Math. Intell. 1999, 21, 18–23. [Google Scholar] [CrossRef]
  2. Epureanu, B.I.; Greenside, H.S. Fractal basins of attraction associated with a damped Newton’s method. SIAM Rev. 1998, 40, 102–109. [Google Scholar] [CrossRef]
  3. Jacobsen, J.; Lewis, O.; Tennis, B. Approximations of continuous Newton’s method: An extension of Cayley’s problem. Electron. J. Diff. Equ. 2007, 15, 163–173. [Google Scholar]
  4. Gutiérrez, J.M.; Hernández, M.Á. A family of Chebyshev-Halley type methods in Banach spaces. Bull. Aust. Math. Soc. 1997, 55, 113–130. [Google Scholar] [CrossRef]
  5. Werner, W. Some improvements of classical iterative methods for the solution of nonlinear equations. Lect. Notes Math. 1981, 878, 426–440. [Google Scholar]
  6. Cordero, A.; Torregrosa, J.R.; Vindel, P. Dynamics of a family of Chebyshev-Halley type methods. Appl. Math. Comput. 2013, 219, 8568–8583. [Google Scholar] [CrossRef]
  7. Gander, W. On Halley’s iteration method. Am. Math. Mon. 1985, 92, 131–134. [Google Scholar] [CrossRef]
  8. Ezquerro, J.A.; Hernández, M.Á.; Romero, N. An extension of Gander’s result for quadratic equations. J. Comput. Appl. Math. 2010, 234, 960–971. [Google Scholar] [CrossRef]
  9. Varona, J.L. Graphic and numerical comparison between iterative methods. Math. Intell. 2002, 24, 37–46. [Google Scholar] [CrossRef]

Share and Cite

MDPI and ACS Style

Gutiérrez, J.M. Numerical Properties of Different Root-Finding Algorithms Obtained for Approximating Continuous Newton’s Method. Algorithms 2015, 8, 1210-1218. https://0-doi-org.brum.beds.ac.uk/10.3390/a8041210

AMA Style

Gutiérrez JM. Numerical Properties of Different Root-Finding Algorithms Obtained for Approximating Continuous Newton’s Method. Algorithms. 2015; 8(4):1210-1218. https://0-doi-org.brum.beds.ac.uk/10.3390/a8041210

Chicago/Turabian Style

Gutiérrez, José M. 2015. "Numerical Properties of Different Root-Finding Algorithms Obtained for Approximating Continuous Newton’s Method" Algorithms 8, no. 4: 1210-1218. https://0-doi-org.brum.beds.ac.uk/10.3390/a8041210

Article Metrics

Back to TopTop