A note on a new cubically convergent one-parameter root solver
aa r X i v : . [ m a t h . NA ] M a y A note on a new cubically convergent one-parameter rootsolver
L. D. Petkovi´c , ∗ , M. S. Petkovi´c , Faculty of Mechanical Engineering, University of Niˇs,A. Medvedeva 14, 18000 Niˇs, Serbia Faculty of Electronic Engineering, University of Niˇs,A. Medvedeva 14, 18000 Niˇs, Serbia
Abstract
A new one-parameter family of iterative method for solving nonlinear equations is con-structed and studied. Two variants, both with cubic convergence, are developed, one forfinding simple zeros and other for multiple zeros of known multiplicities. This family gen-erates a variety of different third order methods, including Halley-like method as a specialcase. Four numerical examples are given to demonstrate convergence properties of the pro-posed methods for multiple zeros and various values of the parameter.
AMS Mathematical Subject Classification (2010):
Key words and phrases:
Solving nonlinear equations; Parametric iterative methods; Con-vergence analysis; Multiple zeros.
Approximating zeros of a given scalar function f belongs to the most important problems thatoccur not only in applied mathematics but also in many disciplines of engineering branches,computer science, physics, finance, and so on. Since there is a vast number of papers and booksdevoted to iterative methods for finding simple and multiple roots of nonlinear equations, see,e.g., [1]–[12], we will not discuss in details characteristics of existing methods.The main goal of this paper is to present two new one-parameter families of iterativemethods for finding simple or multiple zeros of a given function. The main advantages of thisfamily are: 1) the ability to generate a variety of different cubically convergent methods; theproposed family can serve for the construction of very fast iterative methods for approximatingall zeros of a polynomial, see [13].The paper is organized as follows. In Section 2 we construct a one-parameter family ofiterative methods for finding simple roots of nonlinear equations and show that its order ofconvergence is three. In Section 3 the iterative formula for simple zeros is directly used forthe construction of a one-parameter family for finding multiple zeros of the known multiplicity.Results of numerical experiments for several values of the parameter p through three iterationsteps are displayed in Section 4 using four test functions. ∗ Corresponding author
E-mail addresses: [email protected] (L. D. Petkovi´c), [email protected] (M. S. Petkovi´c) One-parameter family for simple zeros
We begin this section with Traub’s result given in [1, Theorem 2.5].
Theorem 1.
Let ψ ( x ) be an iteration function which defines an iterative method forfinding a zero α of multiplicity m of a given function f. Then for these values of m there existsa function ω ( x ) such that ψ ( x ) = x − u ( x ) ω ( x ) , u ( x ) = f ( x ) f ′ ( x ) , ω ( α ) = 0 . (1)In this paper we will restrict our attention to iterative methods with cubic convergence( r = 3). We will often use the abbreviations A λ ( x ) = f ( λ ) ( α ) λ ! f ′ ( α ) ( λ = 1 , , . . . ) . For brevity, we will write sometimes only u instead of u ( x ) for short. The abbreviation AEC (IM) will denote asymptotic error constant of the iterative method (IM). First we presenttwo well known cubically convergent methods free of squares: C ( x ) = x − u ( x ) (cid:16) A ( x ) u ( x ) (cid:17) (Chebyshev’s method) , (2) H ( x ) = x − u ( x )1 − A ( x ) u ( x ) (Halley’s method) . (3)Regarding (1) we note that ω ( u ) = 1 + A u for Chebyshev’s method (2) and ω ( u ) =1 / (1 − A u ) for Halley’s method (3). Therefore, ω ( u ) is a polynomial approximation in (2), while ω ( u ) is a rational approximation in (3). In this paper we will consider a rational approximationto construct a new cubically convergent iterative method in the form G ( u ) = x − u ( x ) · a + p · u ( x )1 + c · u ( x ) . (4)We allow that the coefficients a and c in (4) can be constants as well as some functions of theargument x, while p is a real or complex parameter.First, we start from Schr¨oder-Traub basic sequence { E k } defined recursively by E ( x ) = x − u ( x ) ,E k +1 ( x ) = E k ( x ) − u ( x ) k E ′ k ( x ) ( k ≥ , (5)which defines the generalized iterative method of order k + 1 in the form of a power series, seeTraub [1, Sec. 5.1]. For example (suppressing the argument x ) E ( x ) = x − u − A u (Chebyshev’s method (2)) ,E ( x ) = x − u − A u − (2 A − A ) u ,E ( x ) = x − u − A u − (2 A − A ) u − (5 A − A A + A ) u , etc.We will employ the following assertion. 2 heorem 2. (Schr¨oder [14]) Any root-finding algorithm F n of the order n can be presentedin the form F n ( x ) = E n ( x ) + f ( x ) n η n ( x ) , (6) where η n is a function bounded in α and depending on f and its derivatives. Let G ( x ; p ) be the root-solver to be constructed. According to Theorem 2 we seek for thecoefficients a, p, c in (4) so that G ( x ; p ) = E ( x ) + f ( x ) η ( x ) (7)holds.For two real or complex numbers z and w we will write z = O M ( w ) if | z | = O ( | w | ) (thesame order of their moduli), where O represents the Landau symbol. After the developmentin geometric series we have G ( u ) = x − u · a + pu cu = x − u ( a + pu )(1 − cu + c u + · · · ) = x − au + ( ac − p ) u + · · · . Using this relation and (7), and applying the method of undetermined coefficients, we obtain a = 1 and c = p − A . In this way we have constructed the following one-parameter cubicallyconvergent iterative methodˆ x = G ( x ; p ) = x − u ( x ) (cid:16) p u ( x ) (cid:17) (cid:16) p − A ( x ) (cid:17) u ( x ) , (8)where ˆ x is a new approximation to the zero α of f. Let ε = x − α be the approximation error. To find asymptotic error constant ( AEC forshort), we use the developments in Taylor’ series: f ( x ) = f ′ ( α ) (cid:16) ε + A ε + A ε + O M ( ε ) (cid:17) ,f ′ ( x ) = f ′ ( α ) (cid:16) A ε + 3 A ε + O M ( ε ) (cid:17) ,f ′′ ( x ) = f ′ ( α ) (cid:16) A + 6 A ε + O M ( ε ) (cid:17) . Hence u ( x ) = f ( x ) f ′ ( x ) = ε − A ε + (2 A − A ) ε + O M ( ε ) ,A ( x ) = f ′′ ( x )2 f ′ ( x ) = A + (3 A − A ) ε + (4 A − A A ) ε + O M ( ε ) . Substituting the expressions for u ( x ) and A ( x ) in (8), we obtainˆ ε = ˆ x − α = ( A − A + pA ) ε + O M ( ε ) . (9)From (9) we immediately state the following assertion. Theorem 3.
Assume that x is sufficiently close initial approximation to the zero α of atleast three-time differentiable function f. Then the one-parameter family of iterative methods x k +1 = x k − u ( x k ) (cid:16) p u ( x k ) (cid:17) (cid:16) p − A ( x k ) (cid:17) u ( x k ) ( k = 0 , , . . . ) (10)3 as the order of convergence three for any real or complex parameter p, bounded in magnitude,and AEC (10) = lim k →∞ (cid:12)(cid:12)(cid:12)(cid:12) x k +1 − α ( x k − α ) (cid:12)(cid:12)(cid:12)(cid:12) = (cid:12)(cid:12) A ( α ) − A ( α ) + pA ( α ) (cid:12)(cid:12) is valid. Remark 1.
In a special case when p = 0, from (9) we obtain Halley’s method (3) with AEC (3) = | A ( α ) − A ( α ) | , which is well-known result. Furthermore, when p → ±∞ , then themethod (10) reduce to quadratically convergent Newton’s method x k +1 = x k − u ( x k ) . For thisreason, one should avoid the choice of the parameter p large in magnitude . Remark 2.
We restrict ourselves that the parameter p is a real or complex constant. It isinteresting to consider another special case p = A ( x k ) which could happen accidentally in the k -iteration. Then the iterative process (10) switches to Chebyshev’s method x k +1 = x k − u ( x k ) (cid:16) A ( x k ) u ( x k ) (cid:17) ( k = 0 , , . . . ) , see (2). Since the probability of this case is 0, we will not discuss Chebyshev’s method in whatfollows. However, the described case can be helpful in finding suitable range od the parameter p. Remark 3.
If we choose p = (cid:0) A ( x k ) − A ( x k ) (cid:1) /A ( x k ) in (10), the iterative method (11)becomes x k +1 = x k − u ( x k ) (cid:18) A ( x k ) u ( x k ) A ( x k ) + ( A ( x k ) − A ( x k )) u ( x k ) (cid:19) ( k = 0 , , . . . ) . (11)From (9) we find that the iterative method (11) has the order of convergence equal to 4.Developing the denominator of the expression in the above parenthesis around the point u = 0 , the iterative formula (11) reduced to the fourth order method E ( x k ) given above. Let us now consider the case when α is the zero of f of the known order of multiplicity m ≥ . Note that α is a simple zero for the function F ( x ) = f ( x ) /m . We find the first two derivatives of F : F ′ ( x ) = F ( x ) f ′ ( x ) mf ( x ) , F ′′ ( x ) = F ( x ) · (1 − m ) f ′ ( x ) + mf ( x ) f ′′ ( x ) m f ( x ) . (12)Taking into account the expressions for the derivatives F ′ and F ′′ given by (12), let usreplace u ( x ) and A ( x ) , appearing in (8), by new functions v ( x ) and d ( x ) defined by v ( x ) := F ( x ) F ′ ( x ) = mf ( x ) f ′ ( x ) , d ( x ) := F ′′ ( x )2 F ′ ( x ) = (1 − m ) f ′ ( x ) + mf ( x ) f ′′ ( x )2 mf ( x ) f ′ ( x ) . (13)Then the iteration function (8) becomesˆ x = G m ( x ; p ) = x − v ( x ) (cid:0) pv ( x ) (cid:1) (cid:0) p − d ( x ) (cid:1) v ( x ) . (14)4eplacing the expressions (13) for v ( x ) and d ( x ) we can modify the iteration function (8)for simple zeros to the following iteration function for finding multiple zeros x k +1 = G m ( x k ; p ) = x k − mu ( x k ) (cid:0) mpu ( x k ) (cid:1) m + 2 m (cid:0) p − A ( x k ) (cid:1) u ( x k ) ( k = 0 , , . . . ) . (15)Note that the choice of p = 0 in (15) gives Halley-like method for finding multiple zeros [15] x k +1 = x k − u ( x k ) m + 12 m − A ( x k ) u ( x k ) ( k = 0 , , . . . ) . (16) Theorem 4.
Let x be sufficiently close initial approximation to the zero α of the knownmultiplicity m ≥ of a given function f. Then the iterative method (15) is cubically convergentand
AEC (15) = lim k →∞ (cid:12)(cid:12)(cid:12)(cid:12) x k +1 − α ( x k − α ) (cid:12)(cid:12)(cid:12)(cid:12) = (cid:12)(cid:12)(cid:12)(cid:12) pB m +1 mB m − B m +2 mB m + ( m + 1) B m +1 m B m (cid:12)(cid:12)(cid:12)(cid:12) (17) is valid, where B r = f ( r ) ( α ) /r ! . Proof.
To prove the theorem we will use the iteration function (14). Introduce the errorsof approximations ε = x − α, ˆ ε = ˆ x − α and coefficients C r = m !( m + r )! f ( m + r ) ( α ) f ( m ) ( α ) ( r = 1 , , . . . ) . Then the following developments in Taylor serious are valid: f ( x ) = B m ε m (cid:16) C ε + C ε + C ε + O M ( ε ) (cid:17) ,f ′ ( x ) = B m ε m − (cid:16) m + ( m + 1) C ε + ( m + 2) C ε + ( m + 3) C ε + O M ( ε ) (cid:17) ,f ′′ ( x ) = B m ε m − (cid:16) m ( m −
1) + m ( m + 1) C ε + ( m + 1)( m + 2) C ε +( m + 2)( m + 3) C ε + O M ( ε ) (cid:17) . Using (13) and these expressions, we find v ( x ) = ε − C ε m + (( m + 1) C − mC ) ε m + O M ( ε ) ,d ( x ) = C m − (( m + 1) C − mC ) ε m + O M ( ε ) . (18)According to (18) and Taylor’s series it follows11 + ( p − d ( x )) v ( x ) = 1+ (cid:18) C m − p (cid:19) ε + (cid:0) − C (3 m + 1) − mpC + 6 mC + 2 m p (cid:1) ε m + O (cid:0) ε (cid:1) . Using the last expression and (14), we find after short arrangementˆ ε = ˆ x − α = (cid:16) ( m + 1) C + 2 mpC − mC (cid:17) ε m + O M ( ε ) . (19)With regard to (19) it follows that the order of the iterative method (15) is three. Since C r = B m + r /B m , from (19) we obtain the asymptotic error constant AEC (15) , given by (17). (cid:3) Numerical results
The theoretical order of convergence of the iterative method (15) is three, see Theorem 4. How-ever, it is always convenient to check the convergence behavior in practice. For this reason, inour numerical experiments we have calculated the so-called computational order of convergence r c (COC, for brevity) using the approximate formula r c = log | f ( x k +1 ) /f ( x k ) | log | f ( x k ) /f ( x k − ) | . (20)Note that the formula (20) is a special case of a general formula given in [16]. The testedfunctions are given in Table 1. f ( x ) m x αf ( x ) = (cid:0) x sin x − ( x/ √ (cid:1)(cid:0) x + x + 100 (cid:1) − . f ( x ) = (cid:0) xe x − sin x + 3 cos x + 5) − − . . . .f ( x ) = (cid:0) e x +4 x +5 − (cid:1) sin ( t + 2 − i ) 5 − . . i − if ( x ) = (cid:0) x − sin x (cid:1)
12 0 . Table 1: Tested functions for f − f In Table 2 we have presented the errors of approximations ε k = | z k − α | ( k = 1 , , p. The denotation A ( − h ) means A × − h . The most accurate approximations, obtained after the third iterative step, are boxedin Table 2. We observe that the best results are obtained taking p = − f , p = 0 for f , and p = 1 for f and f . Except the functions listed in Table 1, we have also tested a number of functions of variousstructure. However, we have not found the value of p which defines approximately the bestmethod from the family (15). The influence of the parameter p to the accuracy of approxima-tions to the zeros of a given function is very complex and it is hard to find its optimal valueeven within a particular class of functions. From the discussion given in Remark 1 we canconclude that large values of p are not convenient since the order of convergence decreases andtends to 2. Furthermore, for p = 0 the method (10) reduces to Halley’s method which belongsto the group of cubically convergent methods with very good convergence behavior.Our numerical experiments have shown that optimal parameter p for some classes of func-tions takes negative values belonging to the interval [ − b,
0] ( b > . According to all factsmentioned above and a number of tested functions (see Table 2 for demonstration), we haveconcluded that p should be taken from the interval [ − a, a ] ( a >
0) for relatively small a, say, a ≤ . Following Remark 2, there follows that the choice of p very close to A ( x k ) could be alsogood (taking p ≈ A ( x k ) before running the k -iteration). However, if an initial approximation x is not sufficiently close to the sought zero, the values A ( x ) can be rather crude, producingslow convergence at the beginning of iterative process. Remark 4.
The values of COC r c in Table 2, taken with 3 decimal digits of mantissa(thus, 3.000 does not mean 3) mainly match well the theoretical order 3. However, in somecases unexpected values of r c appear. The explanation is simple: formula (20) works well whenthe approximations x k − , x k , x k +1 are sufficiently close to the zero. One additional iterationmore would give more realistic value of r c . ( x ) = (cid:0) x sin x − ( x/ √ (cid:1)(cid:0) x + x + 100 (cid:1) p | x − α | | x − α | | x − α | r c (20) − . −
2) 1 . −
7) 2 . −
23) 3 . − . −
4) 7 . −
12) 3 . −
36) 3.0000 7 . −
2) 3 . −
6) 3 . −
19) 3 . .
111 1 . −
2) 3 . −
8) 3 . .
172 1 . −
5) 1 . −
17) 2 . f ( x ) = ( xe x − sin x + 3 cos x + 5) − . −
2) 4 . −
4) 2 . −
10) 3 . − . −
2) 1 . −
5) 2 . −
15) 3.0130 7 . −
4) 1 . −
10) 5 . −
31) 3 . . −
2) 1 . −
6) 5 . −
18) 2 . . −
2) 2 . −
5) 2 . −
14) 2 . f ( x ) = ( e x +4 x +5 − sin ( t + 2 − i ) − . −
2) 1 . −
4) 3 . −
12) 3 . − . −
2) 1 . −
5) 1 . −
15) 3.0070 1 . −
2) 2 . −
7) 5 . −
20) 3 . . −
2) 1 . −
7) 9 . −
22) 2 . . −
2) 7 . −
7) 2 . −
19) 2 . f ( x ) = (cid:0) x − sin x (cid:1) − . −
2) 4 . −
8) 1 . −
24) 3 . − . −
3) 5 . −
10) 2 . −
30) 3.0010 1 . −
3) 2 . −
11) 1 . −
34) 3 . . −
4) 6 . −
14) 4 . −
42) 3 . . −
4) 7 . −
13) 6 . −
39) 3 . f − f References [1] J.F. Traub, Iterative Methods for the Solution of Equations, Prentice Hall, New York, 1964.[2] A.M. Ostrowski, Solution of Equations in Euclidean and Banach space, Academic Press, New York,1973.[3] E. Hansen, M. Patrick, A family of root finding methods, Numer. Math. 27 (1977), 257–269.[4] B. Sendov, A. Andreev, N. Kyurkchiev, Numerical Solution of Polynomial Equations, in: Handbookof Numerical Analysis, Vol. 3 (eds P. Ciarlet, J. Lions), Elsevier, Amsterdam, 1993.[5] N.V. Kyurkchiev, Initial Approximations and Root Finding Methods, Wiley-VCH, Berlin, 1998.[6] J.M. McNamee, Numerical Methods for Roots of Polynomials, Part I, Elsevier, Amsterdam, 2007.[7] M.S. Petkovi´c, Point Estimation of Root Finding Methods, Springer, Berlin-Heidelberg, 2008.[8] B. Neta, A.N. Johnson, High-order nonlinear solver for multiple roots, Comp. Math. Appls 55(2008), 2012–2017.
9] C. Chun, B. Neta, A third-order modification of Newton’s method for multiple roots, Appl. Math.Comput. 211 (2009), 474–479.[10] X. Zhou, X. Chen, Y. Song, Construction of higher order methods for multiple roots of nonlinearequations, J. Comput. Appl. Math. 235 (2011), 4199–4206.[11] M. S. Petkovi´c, B. Neta, L. D. Petkovi´c, J. Dˇzuni´c, Multipoint Methods for Solving NonlinearEquations, Elsevier/Academic Press, Amsterdam, 2013.[12] B. Neta, C. Chun, On a family of Laguerre methods to find multiple roots of nonlinear equations,Appl. Math. Comput. 219 (2013), 10987–11004.[13] L. D. Petkovi´c, M. S. Petkovi´c, On a high-order one-parameter family for the simultaneous deter-mination of polynomial roots, Appl. Math. Letters (to appear), DOI:10.1016/j.aml.2017.05.013.[14] E. Schr¨oder, ¨Uber unendlich viele Algorithmen zur Aufl¨osung der Gleichungen, Math. Ann. 2(1870), 317–365.[15] E. Bodewig, Sur la m´ethode Laguerre pour l’approximation des racines de certaines ´equationsalg´ebriques et sur la critique d’Hermite, Indag. Math. 8 (1946), 570–580.[16] L.O. Jay, A note on Q-order of convergence, BIT 41 (2001), 422–429.9] C. Chun, B. Neta, A third-order modification of Newton’s method for multiple roots, Appl. Math.Comput. 211 (2009), 474–479.[10] X. Zhou, X. Chen, Y. Song, Construction of higher order methods for multiple roots of nonlinearequations, J. Comput. Appl. Math. 235 (2011), 4199–4206.[11] M. S. Petkovi´c, B. Neta, L. D. Petkovi´c, J. Dˇzuni´c, Multipoint Methods for Solving NonlinearEquations, Elsevier/Academic Press, Amsterdam, 2013.[12] B. Neta, C. Chun, On a family of Laguerre methods to find multiple roots of nonlinear equations,Appl. Math. Comput. 219 (2013), 10987–11004.[13] L. D. Petkovi´c, M. S. Petkovi´c, On a high-order one-parameter family for the simultaneous deter-mination of polynomial roots, Appl. Math. Letters (to appear), DOI:10.1016/j.aml.2017.05.013.[14] E. Schr¨oder, ¨Uber unendlich viele Algorithmen zur Aufl¨osung der Gleichungen, Math. Ann. 2(1870), 317–365.[15] E. Bodewig, Sur la m´ethode Laguerre pour l’approximation des racines de certaines ´equationsalg´ebriques et sur la critique d’Hermite, Indag. Math. 8 (1946), 570–580.[16] L.O. Jay, A note on Q-order of convergence, BIT 41 (2001), 422–429.