# On Constructing Two-Point Optimal Fourth-Order Multiple-Root Finders with a Generic Error Corrector and Illustrating Their Dynamics

Authors: Young Ik Kim; Young Hee Geum
Journal: Discrete Dynamics in Nature and Society
 (2015)

---

## Abstract

With an error corrector via principal branch of themth root of a function-to-function ratio, we propose optimal quartic-order multiple-root finders for nonlinear equations. The relevant optimal order satisfies Kung-Traub conjecture made in 1974. Numerical experiments performed for various test equations demonstrate convergence behavior agreeing with theory and the basins of attractions for several examples are presented.


---

## Body

## 1. Introduction

Iterative root-finding methods have been constantly developed by numerous researchers [1–3] to solve more accurately the root α of a nonlinear equation f(x)=0 that arises frequently in a scientific world. Classical Newtons method below, (1)xn+1=xn-fxnf′xn,n=0,1,2,…,is widely used for a simple root under normal circumstances, provided that an initial guess x0 is chosen close enough to α. Likewise, modified Newtons method [4] of the form(2)xn+1=xn-mfxnf′xn,n=0,1,2,…,is most popular for a multiple root α with its integer multiplicity m≥1.Definition 1.
Let{x0,x1,x2,…,xn,…} be a sequence converging to α and en=xn-α be the nth iterate error. If there exist real numbers p∈R and b∈R-{0} such that the following error equation holds:(3)en+1=benp+Oenp+1,then b or |b| is called the asymptotic error constant and p is called the order of convergence [5].Note that both methods (1) and (2) are quadratically convergent one-point optimal [3] methods. Two-point higher-order methods for multiple roots can be found in papers [6–11]. Among these, we introduce, here, respectively, in (4), (5), and (6), some interesting works of Soleymani and Babajee [12], Kanwar et al. [13], and Zhou et al. [14] who have recently developed fourth-order multiple-root finders:(4)yn=xn-2mm+2fxnf′xn,xn+1=xn+4mddm2+2m-4-m2v1-m3m-216d2f′ynf′xn-m+2md2·fxnf′xn,d=mm+2m,(5)yn=xn-2mm+2fxnf′xn-pfxn,xn+1=xn-mfxnm23m-2f′yn2-C1f′xnf′yn+C2f′yn22f′xnm-1dmf′xn-mf′xnmf′ynmf′yn-dmm+8f′xn,d=mm+2m,wherep=0,C1=mdm6m2+17m-14,C2=d2m3m3+19m2+16m+16,(6)yn=xn-mfxnf′xn,xn+1=yn-mGw·fxnf′xn=xn-m1+Gwfxnf′xn,w=fynfxnm-1,whereGissufficientlydifferentiableat0withG0=0,G′0=1,G′′0=4.By a close inspection of iterative methods (1), (2), and (4)–(6), we find that an iterative method can be constructed in the following form: (7)en+1=en-Qf,where Qf is implicitly dependent upon xn, for example, with Qf=f(xn)/f′(xn) in (1) and Qf=m{1+G(w)}f(xn)/f′(xn) in (6). We may regard Qf as an error-correcting function. Consequently, it would be natural to call the function Qf as the error corrector. Usually Qf takes the form of Qf=Wf·h, with Wf as a weighting function that is widely known among many researchers. A more generic form of the error corrector Qf will be investigated during the course of developing new quartic-order multiple-root finders.The main aim of this paper is to design new two-step two-point optimal quartic-order multiple-root finders with multiplicity ofm≥1. The first step is to compute yn(xn) using f(xn),f′(xn) usually with a Newton-like method or other. The second step is to update xn in the first step by introducing an error corrector Qf formed by f(xn)/f′(xn) and a principal branch of [f(yn)/f(xn)]1/m. We will check the optimality based on the Kung-Traub conjecture [3] that a multipoint method [15] without memory can achieve its convergence order at most of 2k-1 for k functional evaluations.This paper is comprised of four sections as follows. Following this introductory section, Section2 describes main results with convergence analysis for newly proposed two-point optimal fourth-order multiple-root finders. Principal branch of a logarithmic function plays a crucial role in developing new methods in view of the relation [f(yn)/f(xn)]1/m=exp[1/mlogf(yn)/f(xn)]. The convergence analysis includes the derivation of the error equation for the proposed methods. In Section 3, special cases of error-correcting functions are treated with tabulated results and labeled case numbers. Two types of error-correcting functions are constructed based on bivariate polynomials and rational functions. In the first part of Section 4, with error-correcting functions properly chosen from Section 3, a variety of numerical examples are presented for a wide selection of test functions. A comparison for the convergence behavior is made among the proposed methods and the listed existing methods (4)–(6). The second part of Section 4 discusses related dynamics of maps (8) behind the basins of attraction. Dynamical properties of the proposed methods along with their illustrative basins of attraction are displayed with detailed analyses and comments. Section 5 describes overall conclusion as well as possible future work.


## 2. Main Results

We first assume that a functionf:C→C has a multiple root α with integer multiplicity m≥1 and is analytic [16] in a small neighborhood of α. Then with the concept of error corrector Qf introduced in (7) we propose new two-step iterative multiple-root finders below, given an initial guess x0 sufficiently close to α: for n=0,1,2,…,(8)yn=xn-γ·h,h=fxnf′xn,γ∈R,xn+1=xn-Qfs,h,s=fynfxn1/m,where γ∈N is a parameter and Qf:C2→C is holomorphic [17] in a neighborhood of (λ,0), where λ∈R is to be determined later for optimal quartic-order convergence. Since s is a one-to-m multiple-valued function, we choose s as a principal analytic branch given by s=exp[1/mlog(f(yn)/f(xn))], with Log(f(yn)/f(xn))=log|f(yn)/f(xn)|+iarg(f(yn)/f(xn)) for -π<arg(f′(yn)/f′(xn))≤π; this convention of arg(z) for z∈C agrees with that of log[z] command of Mathematica to be adopted in numerical experiments of Section 4. By means of further inspection of s, we find that λ∈R is characterized in such a way that s=|f(yn)/f(xn)|1/m·exp[i/marg(f(yn)/f(xn))]=λ+O(en), as shown by (15).These suggested methods require one derivative and two functions in order to achieve optimal order of four. In this section, we establish a main theorem describing the convergence analysis regarding proposed methods (8) and find out how to select the parameter γ and the error-correcting function Qf for optimal fourth-order convergence.Theorem 2.
Letf:C→C have a zero α with integer multiplicity m≥1 and be analytic in a small neighborhood of α. Let γ=m and θj=f(m+j)(α)/f(m)(α) for j∈N. Let x0 be an initial guess chosen in a sufficiently small neighborhood of α. Let Qf:C2→C be holomorphic in a neighborhood of (0,0). Let Qij=1/i!j!∂i+j/∂si∂hjQf(s,h)|(s=0,h=0) for 1≤i, j≤4. Suppose that relations Q00=Q02=Q03=Q10=Q12=Q20=Q30=0, Q01=Q11=m, and Q21=2m hold. Then iterative methods (8) are optimal and of order four and possess the following error equation: (9)en+1=-Q04m4-Q13m4m+1θ1-1m2m+12m+2θ1θ2-Q22m4m+12θ12+mm+9-2Q312m4m+13θ13-Q40m4m+14θ14en4+Oen5.Proof.
Three functional evaluations evidently are eligible for optimal convergence order in the sense of Kung-Traub. Hence, it suffices to determine the constant parameterγ and some properties of the error-correcting function Qf for fourth-order convergence. Applying the Taylors series expansion of f about α, we get the following relations:(10)fxn=fmαm!enm1+A1en+A2en2+A3en3+A4en4+Oen5,(11)f′xn=fmαm-1!enm-11+B1en+B2en2+B3en3+B4en4+Oen5,where Ak=m!/(m+k)!θk, Bk=(m-1)!/(m+k-1)!θk, and θk=f(m+k)(α)/f(m)(α) for k=1,2,3,4.
Dividing (10) by (11), we have(12)fxnf′xn=1men-K1en2-K2en3+K3en4+Oen5,where K1=-A1+B1, K2=-A2+A1B1-B12+B2, and K3=-A3+A2B1-A1B12+B13+A1B2-2B1B2+B3.
Lettingt=1-γ/m for convenience with the above relation (12), we obtain(13)yn=α+ten+K11-ten2+K21-ten3+K31-ten4+Oen5.Expanding f(yn) about the multiple root α leads us to relation(14)fyn=fmαenmm!tm-1t+A1t2+K1m1-ten+12tm-2K12m-1mt-12-2A1K1m+1t-1t2+2tA2t3+K2m1-ten2+16tm-3-K13m-2m-1mt-13+3A1m+1t-1t2K12mt-1-2tK2+6t2A3t4+K3m1-t-6K1t-1tA2m+2t3+K2mm-11-ten3+Oen4.Hence, we get(15)s=fynfxn1/m=t+t-1K1m+A1tken-1m2t-1A1K1m2t-1+A12t1+m+m-1t+2mK2m-A2tt+aen2+16m3t-13A12K1m-1-m+4t+3m-1t2+A13tt-12-3mt2-1+2m21+t+t2+6m2-K2m+A3t1+t+t2+A2K11-3t2-6A1mK2m2t-1+K12m1-t+A2t1+m+mt+m-1t2en3+Oen4.Hence, λ mentioned in the description of (8) is selected to be t as desired. It must be emphasized that s denotes a principal analytic branch as mentioned earlier from (8). Noting that O(s-t)=O(h)=O(en), Taylor expansion of Qf(s,h) about (t,0) up to fourth-order terms in both variables yields after retaining up to fourth-order terms in en with Q14=Q23=Q24=Q32=Q33=Q34=Q41=Q42=Q43=Q44=0:(16)Qfs,h=Q00+hQ01+Q11s-t+Q21s-t2+Q31s-t3+h2Q02+Q12s-t+Q22s-t2+h3Q03+Q13s-t+h3Q04+Q10s-t+Q20s-t2+Q30s-t3+Q40s-t4+Oen4.
Hence by substituting (10)–(16) into the proposed method (8), we obtain the error equation as (17)en+1=en-Qfs,h=ψ1en+ψ2en2+ψ3en3+ψ4en4+Oen5,where ψ0=-Q00 and the coefficients ψi, i=1,2,3,4, generally depend on the parameters t, Qij0≤i,j≤4, andθi(i=1,2,…) and the function f(x).
Solvingψ1=0 for Q00, we obtain (18)Q00=0.Substituting (18) into ψ1=0 and simplifying, we have (19)ψ1=m-Q01m-Q10t-12θ1mm+1=0from which m-Q01=0 and Q10(t-1)2=0 must be satisfied independently of θ1. Hence, (20)Q01=m,Q10=0.Substituting (18) and (20) into θ2=0, we get (21)ψ2=-Q02m2-m-Q11t-12θ1m2m+1-Q20t-14θ12m2m+12=0.Solving the above equation independently of θ1 for Q02, Q11 and Q20, we obtain (22)Q11=mt-12,Q02=Q20=0.
Substituting (18), (20), and (22) into ψ3=0 and simplifying, we have (23)ψ3=-Q03m3-Q12t-12m3m+1θ1-Q30t-16m3m+13θ13+-mt-4-2Q21t-14+m2t2m3m+12θ12-tmm+1m+2θ2=0.
To makeψ3=0 independently of θ1 and θ2, we obtain (24)t=Q00=Q10=Q02=Q20=Q12=Q03=Q30=0,Q01=Q11=m,Q21=2m.With the aid of symbolic computation of Mathematica [18], we substitute (18)–(24) into ψ4 to arrive at the following relation: (25)en+1=-Q04m4-Q13m4m+1θ1-1m2m+12m+2θ1θ2-Q22m4m+12θ12+mm+9-2Q312m4m+13θ13-Q40m4m+14θ14en4+Oen5.This completes the proof.


## 3. Special Cases of Error-Correcting Functions

Using relations (18)–(24), the bivariate Taylor polynomial of Qf(s,h) is easily given by (26)Qfs,h=m1+s+2s2+Q31s3h+h2Q22s2+h3Q13s+h4Q04+Q40s4.Here notaions s=[f(yn)/f(xn)]1/m and h=f(xn)/f′(xn) are introduced for simplicity. Special cases of Qf(s,h) are considered here. In each case, relevant coefficients are determined based on relations (18)–(24). If Q22=Q13=Q04=Q40=0, then (26) yields (27)Qfs,h=m1+s+2s2+Q31s3h.Remark 3.
If we take the Taylor polynomial ofG(s) being in the form of G(s)=(1+s+2s2)+Q31/ms3, then G(0)=1,G′(0)=1 and G′′(0)=4 hold. Hence, method (6) is a subcase of this case in which Qf(s,h)=mG(s)·h is satisfied under the restriction of f to R.Although a variety of forms of error-correcting functionsQf(s,h) are available in view of (26), we will limit ourselves to considering two cases of error correctors comprising low-order bivariate polynomials or simple rational functions.Case 1 (Qf with a bivariate polynomial).
In this case,Q31,Q22,Q13,Q04, and Q40 can be regarded as parameters to be chosen to satisfy (28)Qfs,h=m1+s+2s2+Q31s3h+h2Q22s2+h3Q13s+h4Q04+Q40s4.We list typical subcases with selected parameters and Qf(s,h) in Table 1, where SN stands for the corresponding subcase identification number.Table 1
Typical subcases of Case1 with selected parameters.
SN
Q
31
Q
22
Q
13
Q
0
4
Q
4
0
Q
f
s
,
h
1A
0
0
0
0
0
m
1
+
s
+
s
2
⋅
h
1B
1
0
0
0
0
m
1
+
s
+
s
2
+
s
3
⋅
h
1C
2
m
0
0
0
0
m
1
+
s
1
+
2
s
2
⋅
h
1D
0
1
0
0
0
h
3
+
m
1
+
s
+
2
s
2
⋅
h
1E
2
m
0
1
0
0
h
2
s
+
m
1
+
s
1
+
2
s
2
⋅
h
1F
0
0
0
1
0
m
1
+
s
+
h
+
2
m
s
2
⋅
h
1G
2
m
0
0
0
1
s
4
+
m
1
+
s
1
+
s
2
⋅
hCase 2 (Qf with a bivariate polynomial and rational functions).
Consider(29)Qfs,h=a1+a2s+a3s21+b2s+b3s2·h+ds4+qh4+h2sq1h+q2s1+rsh,where a1,a2, and a3 and b2,b3 are determined using (18)–(24) with two of them as free parameters. We find that a1=m,a2=a3-m(1+b3), and b2=-2-b3+a3/m and the remaining parameters d,q,r,q1, and q2 are free to be chosen.
In Table2, we list typical subcases with interesting choices of parameters and Qf(s,h).Table 2
Typical subcases of Case2 with selected parameters.
SN
a
1
a
2
a
3
b
2
b
3
d
q
r
q
1
q
2
Q
f
s
2A
m
-
m
0
−2
0
0
0
0
0
0
m
1
-
s
1
-
2
s
⋅
h
2B
m
0
0
−1
−1
0
0
0
0
0
m
1
-
s
-
s
2
⋅
h
2C
m
0
m
−1
0
0
0
0
0
0
m
1
+
s
2
1
-
s
⋅
h
2D
m
-
m
0
−2
0
1
0
0
0
0
s
4
+
m
1
-
s
1
-
2
s
⋅
h
2E
m
0
0
−1
−1
0
1
0
0
0
h
3
+
m
1
-
s
-
s
2
⋅
h
2F
m
0
m
−1
0
0
0
1
0
1
s
2
h
1
+
s
h
+
m
1
+
s
2
1
-
s
⋅
h
2G
m
-
m
0
−2
0
0
0
1
1
0
s
h
2
1
+
s
h
+
m
1
-
s
1
-
2
s
⋅
h
2H
m
0
0
−1
−1
0
0
0
0
1
s
2
h
+
m
1
-
s
-
s
2
⋅
h


## 4. Numerical Experiments and Dynamics

In this section, we will first describe the computational experiments of proposed methods (8) and then illustrate the complex dynamics [19–26] related to the basins of attraction [27–29] of iterative maps (26) along with comparisons among existing methods.Throughout the experiments, we have moderately assigned 112 significant digits as the minimum number of precision digits, via Mathematica [18] command MinPrecision=112, to achieve the desired accuracy ensuring convergence of the proposed methods. It is necessary to compute en=xn-α with high accuracy for desired numerical results. In case that α is not exact, it is replaced by a more accurate value which has a larger number of significant digits than the assigned $MinPrecision=112.Definition 4 (computational convergence order).
Assume that theoretical asymptotic error constantη=limn→∞(en/|en-1|p) and convergence order p≥1 are known. Define pn=log|en/η|/log|en-1| as the computational convergence order. Note that limn→∞pn=p.
Ifα and xn both mutually have the same accuracy of $MinPrecision=112, then en=xn-α would be nearly zero as n becomes large and thus computing |en+1/en4| would unfavorably cause a numeric overflow. Computed values of xn are accurate up to 112 significant digits. To observe the reliable convergence behavior, we desire α with enough accuracy of 16 digits higher than MinPrecision, which has 128 significant digits. To supply such α, a set of following Mathematica commands are used:(30)sol=FindRootfx,x,x0,PrecisionGoal⟶16+$MinPrecision,WorkingPrecision⟶2∗$MinPrecision;α=sol1,2.Although the number of significant digits of xn and α is 112 and 128, respectively, we list the two values at most up to 15 significant digits due to the limited paper space. We set the error bound ε to 1/2×10-80 satisfying xn-α<ε.Iterative methods (26) with all subcases of both Cases 1 and 2 were, respectively, identified byY1A,Y1B,Y1C,Y1D,Y1E,Y1F, andY1G andY2A,Y2B,Y2C,Y2D,Y2E,Y2F,Y2G, andY2H, beingY-prefixed. Among them, typical three methods have been successfully applied to three test functions shown below:(31)MethodY1D:F1x=4sin2πx+x2-124,m=4,α=-3.MethodY2B:F2x=2x+4i-2π+sinx+2ilogx2+1·sin2x+2i,m=3,α=-2i+π,i=-1.MethodY2E:F3x=sin-1x-12+ex-62,m=2,α≈1.69008688293711,wherelogzz∈Crepresentsaprincipalanalyticbranchsuchthat-π<Imlogz≤π.MethodsY1D,Y2B, andY2E in Table 3 have clearly confirmed quartic-order convergence. Table 3 lists iteration indexes n, approximate zeros xn, residual errors |F(xn)|, errors |en|, and computational asymptotic error constants ηn as well as the theoretical asymptotic error constant η and computational asymptotic convergence order pn. The values of ηn agree up to 10 significant digits with η. Undoubtedly, the computed asymptotic order of convergence well approaches 4. The computational asymptotic error constant reveals a good agreement with the theory developed in Section 3.Table 3
Convergence for test functionsF1(x)-F3(x) with methods  M1F, M2E, and M2F.
M
T
f
i
n
x
n
|
F
(
x
n
)
|
e
n
=
x
n
-
α
|
e
n
/
e
n
-
1
4
|
η
p
n
Y
1
C
F
1
0
−3.3
5.68515
0.3
1
−3.00038943156575
1.212
×
1
0
-
11
3.894
×
1
0
-
4
0.04807797108
0.06814101693
4.28967
2
−3.00000000000000
3.173
×
1
0
-
57
1.566
×
1
0
-
15
0.06811440433
4.00005
3
−3.00000000000000
1.495
×
1
0
-
239
4.104
×
1
0
-
61
0.06814101693
4.00000
4
−3.00000000000000
0.0
×
1
0
-
443
0.0
×
1
0
-
111
Y
1
D
F
2
0
3.1
-
1.95
∗
0.000338602
0.0650381
1
3.14216523070554
2.00080488739465
1.213
×
1
0
-
9
9.877
×
1
0
-
4
55.20587710
0.02458677729
1.17628
2
3.14159265358979
2.00000000000000
1.614
×
1
0
-
41
2.340
×
1
0
-
14
0.02459092337
3.99998
3
3.14159265358979
2.00000000000000
5.066
×
1
0
-
169
4.794
×
1
0
-
57
0.02458677729
4.00000
4
3.14159265358979
2.00000000000000
0.0
×
1
0
-
557
0.0
×
1
0
-
111
Y
1
E
F
3
0
1.6
0.400529
0.0900869
1
1.70503564060025
0.0129610
0.0149488
226.9649843
0.3269289253
1.28174
2
1.69008689870921
1.408
×
1
0
-
14
1.577
×
1
0
-
8
0.3158415015
4.00821
3
1.69008688293711
2.317
×
1
0
-
62
2.023
×
1
0
-
32
0.3269289253
4.00000
4
1.69008688293711
0.0
×
1
0
-
222
0.0
×
1
0
-
111
MT
means method, pn=log⁡|en/η|/log⁡|en-1|, and ∗3.1-1.95=3.1-1.95i.Additional functions below are tested to further confirm the convergence of methods (8): (32)f1x=x2x+log1+x2,α=0,m=3,x0=0.1,f2x=5-6x2+cos2x-sin3x25-6x2+cos2x1-2sinx-sinx2,α≈0.84777738787,m=4,x0=1.0,f3x=1+cosx2-1-xlogx2-π4x2-1-π2,α=1+π,m=6,x0=2.0,f4x=x2+e-x2+sin1x-37,α=1.51115415815562,m=7,x0=1.0,f5x=x-3x3cosπx6+11+x2-115+432x-23,α=2,m=5,x0=1.7,f6x=x-22+316+1x2logx-22+1916x-22+316,α=2-i34,m=2,x0=1.93-0.41i,f7x=x21-logx-1x8,α=1,m=8,x0=0.9;herelogzz∈Crepresentsaprincipalanalyticbranchsuchthat-π<Imlogz≤π.In Table8, we directly compare numerical errors |xn-α| of proposed methodsY1B,Y1C,Y1E, andY1G in Case 1 with those of existing optimal fourth-order multiple-root finders. AbbreviationsKan,Zhou, andSol denote existing optimal fourth-order multiple root finders obtained by Kanwar et al. (5), Zhou et al. (6) with G(w)=(1+w+2w2), and Soleymani and Babajee (4), respectively. The least errors within the prescribed error bound are highlighted in boldface. MethodY1E shows best convergence for f1,f5, and f7 and methodY1G for f2,f4 and methodSol for f3,f6.Likewise, in Table9, we compare numerical errors |xn-α| of proposed methodsY2A,Y2D,Y2F, andY2G in Case 2 with those of existing optimal fourth-order multiple-root finders. MethodY2G shows best convergence for f1,f3,f5, and f7 and methodY2D for f4,f6, while method shows best convergenceY2F for f2.Recall thatlimn→∞xn+1-α/(xn-α)p=gf′(α)/p! for an iterative method xn+1=gf(xn) converging to α with order p; note that the iteration function gf(x) is dependent upon a given nonlinear equation f(x). Consequently, with a given order of convergence p, the local convergence behavior is dependent on the function f(x), an initial value x0, and a multiple zero α itself. Hence, for given test functions, zeros, and initial values, no method is expected to always show better performance than the others.The efficiency index [5] is defined by EI=p1/d, where p is the order of convergence and d is the number of distinct functional or derivative evaluations per iteration. The proposed methods (8) as well as all other listed methods have the same EI=41/3≈ 1.587 agreeing with Kung-Traub optimality conjecture.At this point, we are now ready to discuss the dynamical behavior of iterative maps (8). It is important to properly select initial values close to zero α for guaranteed convergence of iterative methods. It is not, however, an easy task to determine how close the initial values are to zero α. A convenient way of employing stable initial values is to directly use visual basins of attraction. In view of inspection of the area of convergence on the basins of attraction, the larger area of convergence would imply a better method. Clearly measuring the size of area of convergence requires a quantitative analysis. To this end, we provide Tables 5–8 featuring statistical data for the average number of iterations per point and the number of divergent points includingCPU time. In the following 4 examples, we select a 6 by 6 square region centered at the origin and containing all the zeros of the given test polynomial functions. We then take a 600 × 600 uniform grid in the square to display initial points for the iterative methods via basins of attraction. We color the grid point, according to the iteration number for convergence and the root it converged to. This way we can find out if the method converged within the maximum number of iteration allowed and if it converged to the root closer to the initial grid point. To continue our discussion, let us first identify four members of iterative map (8) associated with Cases1D,1G,2D, and2G byY1D,Y1G,Y2D, andY2G, respectively.For illustrating the complex dynamics of (8) with the desired basins of attraction, we take applied to various polynomials having multiple roots with multiplicity m=2,3,4,7. Statistical data for the basins of attraction are tabulated in Tables 4–7. In this tables, abbreviationsCPU,TCON,AVG, andTDIV denote the value ofCPU time for convergence, the value of total convergent points, the value of average iteration number for convergence, and the value of divergent points.Table 4
Typical Example 1 withP1(z)=(9z2+z-7)2, m=2.
Methods
CPU
TCON
AVG
TDIV
Kan
104
360,000
5.3
0
Sol
74
352,462
5.6
7538
Zhou
54
360,000
5.9
0
Y1D
83
360,000
5.5
5758
Y1G
77
359,998
5.7
2
Y2D
100
360,000
5.3
0
Y2G
90
360,000
5.0
0Table 5
Typical Example 2 withP2(z)=z3(z+1), m=3.
Method
CPU
TCON
AVG
TDIV
Kan
108
360,000
5.2
0
Sol
96
359,194
5.1
806
Zhou
69
360,000
5.3
0
Y1D
333
360,000
4.8
0
Y1G
166
360,000
4.8
0
Y2D
526
360,000
4.6
0
Y2G
251
359,996
4.7
4Table 6
Typical Example 3 withP3(z)=(z3-z)4, m=4.
Method
CPU
TCON
AVG
TDIV
Kan
320
360,000
6.4
0
Sol
201
359,860
6.2
140
Zhou
80
359,988
7.4
2
Y1D
112
360,000
6.7
0
Y1G
138
360,000
6.4
0
Y2D
172
360,000
5.7
0
Y2G
187
360,000
5.5
0Table 7
Typical Example 4 withP4(z)=(z2-1)7, m=7.
Method
CPU
TCON
AVG
TDIV
Kan
157
360,000
5.5
0
Sol
146
359,960
5.4
40
Zhou
146
360,000
5.9
0
Y1D
58
360,000
5.4
0
Y1G
61
359,994
5.7
0
Y2D
62
359,998
5.3
2
Y2G
74
360,000
4.7
0Table 8
Comparison of|xn-α| among existing and proposed (Case 1) multiple-root finders.
f,x0; m
|
x
n
-
α
|
Kan
Sol
Zhou
Y1B
Y1C
Y1E
Y1G
f
1, 0.1; 3
|
x
1
-
α
|
1.23e − 5
1.12e − 5
1.54e − 5
1.47e − 5
1.11e − 5
1.02e − 5
1.05e − 5
|
x
2
-
α
|
3.44e − 21
2.15e − 21
1.28e − 20
1.00e − 20
2.31e − 21
1.48e − 21
1.68e − 21
|
x
3
-
α
|
2.08e − 83
2.89e − 84
6.00e − 81
2.10e − 81
4.23e − 84
6.58e − 85
1.09e − 84
f
2, 1.0; 4
|
x
1
-
α
|
3.47e − 5
4.61e − 5
7.32e − 5
6.75e − 5
6.49e − 6
2.53e − 5
1.92e − 5
|
x
2
-
α
|
2.33e − 20
1.32e − 19
2.92e − 18
1.92e − 18
2.23e − 21
9.36e − 21
2.16e − 21
|
x
3
-
α
|
4.75e − 81
9.15e − 78
7.36e − 72
1.26e − 72
3.14e − 83
1.74e − 82
3.46e − 85
f
3, 2.0; 6
|
x
1
-
α
|
8.68e − 6
7.39e − 6
1.64e − 5
1.58e − 5
9.40e − 6
9.39e − 6
1.44e − 5
|
x
2
-
α
|
4.03e − 20
1.87e − 20
9.76e − 19
8.12e − 19
5.59e − 20
5.56e − 20
5.01e − 19
|
x
3
-
α
|
1.87e − 77
7.68e − 79
1.20e − 71
5.53e − 72
6.99e − 77
6.86e − 77
7.29e − 73
|
x
4
-
α
|
0.0e − 111
0.0e − 111
0.0e − 111
0.0e − 111
0.0e − 111
0.0e − 111
0.0e − 111
f
4, 1.0; 7
|
x
1
-
α
|
4.93e − 5
4.64e − 5
5.65e − 5
5.52e − 5
3.81e − 5
3.80e − 5
3.29e − 5
|
x
2
-
α
|
6.98e − 18
5.14e − 18
1.44e − 17
1.28e − 17
1.85e − 18
1.83e − 18
8.25e − 19
|
x
3
-
α
|
2.79e − 69
7.79e − 70
6.19e − 68
3.70e − 68
1.03e − 71
9.87e − 72
3.24e − 73
|
x
4
-
α
|
0.0e − 111
0.0e − 111
0.0e − 111
0.0e − 111
0.0e − 111
0.0e − 111
0.0e − 111
f
5, 1.7; 5
|
x
1
-
α
|
2.89e − 3
2.87e − 3
1.26e − 3
1.25e − 3
1.14e − 3
1.13e − 3
1.16e − 3
|
x
2
-
α
|
2.13e − 10
2.02e − 10
4.83e − 12
4.53e − 12
2.39e − 12
2.30e − 12
2.80e − 12
|
x
3
-
α
|
6.46e − 39
5.09e − 39
1.05e − 45
7.91e − 46
4.59e − 47
3.90e − 47
9.67e − 47
|
x
4
-
α
|
0.0e − 111
0.0e − 111
0.0e − 111
0.0e − 111
0.0e − 111
0.0e − 111
0.0e − 111
f
6, 1.93, −0.41i; 2
|
x
1
-
α
|
1.22e − 4
9.90e − 5
2.51e − 4
2.26e − 4
1.53e − 4
1.57e − 4
1.98e − 4
|
x
2
-
α
|
8.28e − 16
2.91e − 16
3.03e − 14
1.79e − 14
2.48e − 15
2.87e − 15
7.78e − 15
|
x
3
-
α
|
1.74e − 60
2.20e − 62
6.38e − 54
6.98e − 55
1.41e − 29
1.90e − 29
1.39e − 28
|
x
4
-
α
|
0.0e − 111
0.0e − 111
0.0e − 111
0.0e − 111
8.31e − 58
8.31e − 58
0.0e − 111
f
7, 0.9; 8
|
x
1
-
α
|
5.76e − 5
5.42e − 5
6.65e − 5
6.51e − 5
4.33e − 5
4.32e − 5
4.94e − 5
|
x
2
-
α
|
9.45e − 18
7.03e − 18
1.96e − 17
1.75e − 17
2.13e − 18
2.11e − 18
4.29e − 18
|
x
3
-
α
|
6.82e − 69
1.99e − 69
1.48e − 67
9.20e − 68
1.25e − 71
1.20e − 71
2.45e − 70
|
x
4
-
α
|
0.0e − 111
0.0e − 112
0.0e − 112
0.0e − 112
0.0e − 112
0.0e − 112
0.0e − 112
Here, 1.23e − 5 denotes 1.23×10-5.Table 9
Comparison of|xn-α| among existing and proposed (Case 2) multiple-root finders.
f,x0; m
|
x
n
-
α
|
Kan
Sol
Zhou
Y2A
Y2D
Y2F
Y2G
f
1, 0.1; 3
|
x
1
-
α
|
1.23e − 5
1.12e − 5
1.54e − 5
6.34e − 6
5.71e − 6
1.02e − 5
5.38e − 6
|
x
2
-
α
|
3.44e − 21
2.15e − 21
1.28e − 20
1.19e − 22
6.59e − 23
1.48e − 21
5.20e − 23
|
x
3
-
α
|
2.08e − 83
2.89e − 84
6.00e − 81
1.52e − 89
1.16e − 90
6.58e − 85
4.53e − 91
f
2, 1.0; 4
|
x
1
-
α
|
3.47e − 5
4.61e − 5
7.32e − 5
2.86e − 5
3.71e − 5
2.14e − 5
3.11e − 5
|
x
2
-
α
|
2.33e − 20
1.32e − 19
2.92e − 18
5.43e − 10
9.18e − 10
4.47e − 21
6.42e − 10
|
x
3
-
α
|
4.75e − 81
9.15e − 78
7.36e − 72
3.98e − 39
4.10e − 38
8.37e − 84
8.66e − 39
f
3, 2.0; 6
|
x
1
-
α
|
8.68e − 6
7.39e − 6
1.64e − 5
8.53e − 7
5.87e − 6
9.12e − 6
8.45e − 7
|
x
2
-
α
|
4.03e − 20
1.87e − 20
9.76e − 19
5.55e − 25
6.51e − 21
4.99e − 20
5.31e − 25
|
x
3
-
α
|
1.87e − 77
7.68e − 79
1.20e − 71
9.94e − 98
9.82e − 81
4.47e − 77
8.25e − 98
f
4, 1.0; 7
|
x
1
-
α
|
4.93e − 5
4.64e − 5
5.65e − 5
1.59e − 5
1.07e − 5
3.69e − 5
1.58e − 5
|
x
2
-
α
|
6.98e − 18
5.14e − 18
1.44e − 17
2.13e − 20
2.06e − 21
1.60e − 18
2.07e − 20
|
x
3
-
α
|
2.79e − 69
7.79e − 70
6.19e − 68
6.85e − 80
2.86e − 84
5.77e − 72
6.11e − 80
f
5, 1.7; 5
|
x
1
-
α
|
2.89e − 3
2.87e − 3
1.26e − 3
1.00e − 3
1.01e − 3
1.15e − 3
9.93e − 4
|
x
2
-
α
|
2.13e − 10
2.02e − 10
4.83e − 12
8.63e − 13
1.08e − 12
2.45e − 12
8.23e − 13
|
x
3
-
α
|
6.46e − 39
5.09e − 39
1.05e − 45
4.77e − 49
1.43e − 48
5.13e − 47
3.92e − 49
|
x
4
-
α
|
0.0e − 111
0.0e − 111
0.0e − 111
0.0e − 111
0.0e − 111
0.0e − 111
0.0e − 111
f
6, 1.93, −0.41i; 2
|
x
1
-
α
|
1.22e − 4
9.90e − 5
2.51e − 4
4.88e − 5
1.04e − 4
1.41e − 4
5.34e − 5
|
x
2
-
α
|
8.28e − 16
2.91e − 16
3.03e − 14
8.58e − 18
2.95e − 16
1.83e − 15
1.34e − 17
|
x
3
-
α
|
1.74e − 60
2.20e − 62
6.38e − 54
1.69e − 34
1.89e − 62
7.75e − 30
4.14e − 34
|
x
4
-
α
|
0.0e − 111
0.0e − 111
0.0e − 111
0.0e − 111
0.0e − 111
1.38e − 58
0.0e − 111
f
7, 0.9; 8
|
x
1
-
α
|
5.76e − 5
5.42e − 5
6.65e − 5
1.50e − 5
2.11e − 5
4.24e − 5
1.49e − 5
|
x
2
-
α
|
9.45e − 18
7.03e − 18
1.96e − 17
1.06e − 20
6.43e − 20
1.98e − 18
1.03e − 20
|
x
3
-
α
|
6.82e − 69
1.99e − 69
1.48e − 67
2.64e − 81
5.52e − 78
9.47e − 72
2.35e − 81
Here, 1.23e − 5 denotes 1.23×10-5.In the first example, we have taken the following polynomial:(33)P1z=9z2+z-72,whose roots z=-0.939221,0.82811 are both real with multiplicity m=2. Based on Table 4 and Figure 1, we can find out thatY2G is best in view of lowerAVG andTDIV, followed byKan andY2D. As can be seen in Figure 1,Sol andY1D have shown considerable amount of black points. These points causing divergence behavior were expected from the last column of Table 4. The best result forCPU is byZhou and the worst one is byKan.Figure 1
(a) forKan, (b) forSol, and (c) forZhou and (d) for1D, (e) for1G, (f) for2D, and (g) for2G, for the roots of the polynomial P1(z)=(9z2+z-7)2.
(a)
(b)
(c)
(d)
(e)
(f)
(g)Our next sample has triple roots. The polynomial has three identical roots at the origin:(34)P2z=z3z+1.The results are listed in Table 5 and Figure 2. The methodY2D performs best in view of lowerAVG andTDIV. These proposed methodsY1D,Y1G,Y2D, andY2G appear to perform better thanKan,Sol, andZhou. As can be seen in Figure 2,Sol has shown considerable amount of black points, whileY2G has shown a few black ones. The best result forCPU is byZhou and the worst one is byY2D.Figure 2
(a) forKan, (b) forSol, and (c) forZhou and (d) for1D, (e) for1G, (f) for2D, and (g) for2G, for the roots of the polynomial P2(z)=z3(z+1).
(a)
(b)
(c)
(d)
(e)
(f)
(g)As a third example, we take the following polynomial whose roots are all of multiplicity four:(35)P3z=z3-z4.As can be seen in Figure 3, Sol has shown considerable amount of black points, whileZhou has shown several black ones. The best result forCPU is byZhou and the worst one is byKan.Figure 3
(a) forKan, (b) forSol, and (c) forZhou and (d) for1D, (e) for1G, (f) for2D, and (g) for2G, for the roots of the polynomial P3(z)=(z3-z)4.
(a)
(b)
(c)
(d)
(e)
(f)
(g)The results are presented in Table6 and Figure 3. The method2G is best in view of lowerAVG andTDIV. In the last example, we use the polynomial having two roots of unity (36)P4z=z2-17. The results are given in Table 7 and Figure 4. The method2G is best in view of lowerAVG andTDIV. As can be seen in Figure 4,Sol has shown many black points, while has shownY2G a few black ones. The best result forCPU is byY1D and the worst one is byKan.Figure 4
(a) forKan, (b) forSol, and (c) forZhou and (d) for1D, (e) for1G, (f) for2D, and (g) for2G, for the roots of the polynomial P4(z)=(z2-1)7.
(a)
(b)
(c)
(d)
(e)
(f)
(g)


## 5. Conclusion

We have shown that optimal quartic-order multiple-root finders are constructed with an error-correcting function generated byf(xn)/f′(xn) and a principal branch of [f(yn)/f(xn)]1/m along with the derivation of their relevant error equation. To select initial values near zero α for ensured convergence, we should investigate the dynamics behind the basins of attraction of the corresponding iterative maps applied to various polynomials. Our future work will strengthen the current approach to pursue higher-order methods by controlling free parameters of the error-correcting functions that enhance relevant basins of attraction for a wide class of polynomials.



---

*Source: 378517-2015-12-16.xml*