Developing an Optimal Class of Generic Sixteenth-Order Simple-Root Finders and Investigating Their Dynamics

Developed here are sixteenth-order simple-root-finding optimal methods with generic weight functions. Their numerical and dynamical aspects are investigated with the establishment of a main theorem describing the desired optimal convergence. Special cases with polynomial and rational weight functions have been extensively studied for applications to real-world problems. A number of computational experiments clearly support the underlying theory on the local convergence of the proposed methods. In addition, to investigate the relevant global convergence, we focus on the dynamics of the developed methods, as well as other known methods through the visual description of attraction basins. Finally, we summarized the results, discussion, conclusion, and future work.


Introduction
The governing equations of real-world natural phenomena are often described by nonlinear equations whose exact solutions are infeasible due to their inherent complexities. The attainment of precise numerical approximations to the roots of such complicated nonlinear functions is important for many scientific fields. The classical second-order Newton's method is best known as the numerical root-finder for the governing equations. For several decades, many authors [1][2][3][4][5][6][7][8][9][10][11] have developed higher-order multipoint methods. If an iterative root-finding method satisfies Kung-Traub's conjecture [12], then it is said to be optimal. A few authors [12][13][14] have recently established optimal sixteenth-order methods, despite the lack of applicability to real-life nonlinear governing equations due to their algebraic complexities, not only to emphasize the theoretical importance of developing extremely high-order methods, but also to apply them to root-finding of real-world nonlinear problems, we strongly desire to establish a new optimal family of sixteenth-order simple-root finders that are comparable to or competitive against the existing methods.
For the sake of comparison with the new optimal family of methods to be proposed in this paper, we introduce existing three optimal sixteenth-order Equations [12][13][14] respectively given by Equations (1), (2), and (4) below. where Maroju-Behl-Motsa method (MBM): where G : C 2 → C is an analytic function in a neighborhood of (0, 0) satisfying G 00 = 1, G 01 = 2, G 10 = 1, G 02 = 10 − 4β, G 11 = 4, G 03 = 12(β 2 − 6β + 6), with G ij = ∂ i+j ∂u i ∂s j G(u, s)| (u=0,s=0) for i, j = 0, 1, 2, 3, and θ 5 is given by the following: with u 1 = f (w n )[b 2 n f (x n ) + b n f (x n ) − c n f (z n )] + a n [ f (x n ) − a n f (x n )] f (z n ), u 2 = a n b n c n f (x n )[ f (y n ) − f (x n )] + c n f (y n ) f (x n )(a n − b n ), v 1 = f (y n )[b n f (w n ){b 2 n f (x n ) + b n f (x n ) − c n f (z n )} + {a 3 n f (x n ) + c n a n f (w n ) − a 2 n f (x n )} f (z n )], v 2 = a 2 n b 2 n c n f (x n ) 2 {2 f (y n ) − f (x n )} + a n b n c n (2a n − c n ) f (x n ) f (y n ) f (x n ) + c n {a n b n − a n c n − b 2 n } f (y n ) f (x n ) 2 , a n = x n − z n , b n = w n − x n , c n = w n − z n . • Sharma-Argyros-Kumar method (SAK): n f (x n ) 1 f (y n ) y n f (y n ) y 2 n f (y n ) 1 f (z n ) z n f (z n ) z 2 n f (z n ) 1 f (w n ) w n f (w n ) w 2 n f (w n ) , , with | • | denoting the determinant of •.
In order to develop the desired competitive optimal sixteenth-order simple-root finders, we seek a class of iterative methods with generic weight functions: where , v = f (w n ) f (z n ) ; Q f : C → C is analytic [15] in a neighborhood of 0, K f : C 2 → C holomorphic [16,17] in a neighborhood of (0, 0), and J f : C 3 → C holomorphic in a neighborhood of (0, 0, 0).
Our main aim is to devise an optimal class of sixteenth-order methods by characterizing the algebraic structure of weight functions Q f (s), K f (s, u), and J f (s, u, v), as well as to explore their dynamics through basins of attractions [19] behind the extraneous fixed points [20] with applications to f (z) = (z − a) m (z − b) m . The right side of final substep of (5) conveniently locates extraneous fixed points from the roots of the combined weight function 1 + sQ f (s) + suK f (s, u) + suvJ f (s, u, v).
It is important that we seek appropriate parameters for which the attractor basin contains larger regions of convergence. A motivation undertaking this research was to extensively study the dynamics behind the extraneous fixed points, which would impact on the relevant dynamics of the iterative methods by producing attractive, indifferent, repulsive, and other chaotic orbits. The entire complex plane is composed of two symmetrical half-planes whose boundary is the imaginary axis. We display the convergence behavior in the dynamical planes through the attractor basins within a square region centered at the origin. We also want to make the dynamics behind the extraneous fixed points on the imaginary axis less influenced by the possible periodic or chaotic attractors. Therefore, in addition to general cases with free parameters leading us to simple weight functions, we preferably include some interesting cases with free parameters, chosen for the purely imaginary extraneous fixed points.
In Section 2, the main theorem is presented with three weight functions, Q f , K f , and J f , containing free parameters. Section 3 considers special cases of weight functions. Section 4 discusses the extraneous fixed points, including purely imaginary ones and investigates their stabilities. Section 5 presents numerical experiments as well as illustrates the relevant dynamics and summarizes the overall work together with future work.

Methods and Materials
The main theorem is established by describing the error equation and the asymptotic error constant with relationships among generic weight functions Q f (s), K f (s, u), and J f (s, u, v): Theorem 1. Suppose that f : C → C has simple root α and is analytic in a neighborhood of α. Let c j = f (j) (α) j! f (α) for j = 2, 3, · · · . Let x 0 be an initial guess selected in a sufficiently small region containing α. Assume L f : C → C is analytic in a neighborhood of 0.
. Let K f : C 2 → C be holomorphic in a neighborhood of (0, 0). Let J f : C 3 → C be holomorphic in a neighborhood of (0, 0, 0).
are fulfilled, then scheme (5) leads to an optimal class of sixteenth-order root-finders possessing the following error equation: with e n = x n − α for n = 0, 1, 2, · · · , Proof. Since Scheme (5) employs five functional evaluations, namely, f (x n ), f (x n ), f (y n ), f (z n ), and f (w n ), optimality can be achieved if the corresponding convergence order is 16. In order to induce the desired order of convergence, we begin by the 16th-order Taylor series expansion of f (x n ) about α: It follows that For brevity of notation, we abbreviate e n with e. Using Mathematica [21], we find: (15) where Since f (y n ) = f (x n )| e n →(y n −α) , we are led to an expression: where D i = D i (c 2 , c 3 , · · · , c 16 ) for 5 ≤ i ≤ 16. Hence, we have: where In the third substep of Scheme (5), w n = O(e 8 ) can be achieved based on Kung-Traub's conjecture. To reflect the effect on w n from z n in the second substep, we need to expand z n up to eighth-order terms; hence, we carry out a sixth-order Taylor expansion of Q f (s) about 0 by noting that s = O(e) and f (y n ) where Q j = 1 j! d j ds j Q f (s) for 0 ≤ j ≤ 6. As a result, we come up with: where W i = W i (c 2 , c 3 , · · · , c 16 , Q 0 , · · · , Q 6 ) for 4 ≤ i ≤ 16. Selecting Q 0 = 1 and Q 1 = 2 leads us to an expression: On account of the fact that f (z n ) = f (x n )| e n →(z n −α) , we deduce: where F i = F i (c 2 , c 3 , · · · , c 16 , Q 2 , · · · , Q 6 ) for 5 ≤ i ≤ 16. Consequently, we find: In the last substep of Scheme (5), x n+1 = O(e 16 ) can be achieved based on Kung-Traub's conjecture. To reflect the effect on x n+1 from w n in the third substep, we need to expand w n up to sixteenth-order terms; hence, we carry out a 12th-order Taylor expansion of K f (s, u) about (0, 0) by noting that: K f (s, u) = K 00 + K 10 s + K 20 s 2 + K 30 s 3 + K 40 s 4 + K 50 s 5 + K 60 s 6 + K 70 s 7 + K 80 s 8 + K 90 s 9 + K 100 s 10 + K 110 s 11 + K 120 s 12 + (K 01 + K 11 s + K 21 s 2 + K 31 s 3 + K 41 s 4 + K 51 s 5 + K 61 s 6 + K 71 s 7 + K 81 s 8 + K 91 s 9 + K 101 s 10 )u+ Substituting , and K f (s, u) into the third substep of (5) leads us to: where Γ i = Γ i (c 2 , c 3 , · · · , c 16 , Q 2 , · · · , Q 6 , K j ), for 5 ≤ i ≤ 16, 0 ≤ j ≤ 12 and 0 ≤ ≤ 6. Thus K 00 = 1 immediately annihilates the fourth-order term. Substituting K 00 = 1 into Γ 5 = 0 and solving for K 10 , we find: Continuing the algebraic operations in this manner at the i-th (6 ≤ i ≤ 7) stage with known values of K j , we solve Γ i = 0 for remaining K j to find: Using K 00 = 1, K 10 = 2, K 20 = 1 + Q 2 , K 01 = 1 yields where τ, µ, ν are described in (12). Consequently, we obtain: where η 0 and η 1 are described in (12) and T i = T i (c 2 , c 3 , · · · , c 16 , Q 2 , · · · , Q 6 ) for 5 ≤ i ≤ 16.
To compute the last substep of Scheme (5), it is necessary to have an eighth-order Taylor expansion of J f (s, u, v) about (0, 0, 0) due to the fact that f (w n ) f (x n ) = O(e 8 ). It suffices to expand J f up to eighth-, Substituting (5), we arrive at: Substituting J 000 = 1 into Ω 9 = 0 and solving for J 100 , we find: Continuing the algebraic operations in the same manner at the i-th (10 ≤ i ≤ 15) stage with known values of J jk , we solve Ω i = 0 for remaining J jk to find: Upon substituting Relation (31) into Ω 16 , we finally obtain: where τ, ν, µ, and Ψ as described in (12). This completes the proof.
Among various possible weight functions Q f (s), K f (s, u), and J f (s, u, v), we restrict the current study to simple ones employing polynomials, as well as low-order rational functions. We consider the first case with simple polynomial weight functions by setting all available free parameters to the desired values, as follows:

Case 1: Polynomial weight functions
After setting all free parameters to zero, we obtain: After setting , and all other free parameters to zero, we obtain: After setting Q 2 = −1 and all other free parameters to zero, we obtain: After setting Q 2 = −1, Q 3 = 6 and all other free parameters to zero, we obtain: As a second case, we restrict ourselves to considering all three rational-type weight functions with real coefficients: where a i , b i , r i , q i ∈ R are to be determined for optimal sixteenth-order convergence; the coefficients of Q f and K f are already selected to satisfy the constraints stated in Theorem 1, while the coefficients of J f should satisfy the constraints J 100 = 2 and affine relations described by Relation (31). These 25 constraints determine 25 relations among the 56 coefficients r i , q i , (1 ≤ i ≤ 28), from which 25 out of 56 coefficients may be solved as an appropriate affine combination of the remaining 31 coefficients. For ease of analysis, we employ some simple forms of K f by appropriate choices of the free parameters as follows: To consider simple forms of J f connected with K f via Relation (31), we first conveniently set all 31 free parameters q 8 , q 12 -q 15 , q 17 -q 21 , q 24 -q 28 , r 6 -r 8 , r 11 , r 13 -r 15 , r 17 , r 19 -r 21 , r 23 -r 26 , r 28 to zero. Then, we get seven forms of J f matching with seven forms, (A)-(G), of (39) in order as follows: where A 0 = −725 − 760s + 1725s 2 + 4140s 3 − 5625s 4 − 1914s 5 − 9743s 6 − 22,386s 7 and B 0 = −725 + 690s + 3970s 2 + 1450s 3 − 15,775s 4 + 4986s 5 . We denote these seven subcases described by weight functions Q f (s) = 1 1−2s , K f (s, u) in (39) and J f (s, u, v) in (40) by Cases 2A-2G in order. For example, Case 2A takes the form of: As the final case, with four second-order rational weight functions K f in (39), i.e., with K f given by (39)-(A), (39)-(C), (39)-(D), (39)-(G), we will pursue possible forms of J f (s, u, v) whose all extraneous fixed points are purely imaginary, when prototype polynomial f (z) = z 2 − 1 is applied. According to our previous studies [22,23] on the dynamics of root finders for nonlinear equations behind the purely imaginary extraneous fixed points, the relevant convergence behavior is improved compared to the usual root finders. This convergence advantage inspires us to investigate Case 3 below underlying the presence of purely imaginary extraneous fixed points.
i=15 r i s i−15 and the specific choice of parameter β for K f and determination of the 64 coefficients q i , r i of J f are described below. Relationships were sought among all free parameters of J f (s, u, v), giving us a simple governing equation for extraneous fixed points of the proposed Family of Methods (5).
To this end, we first express s and u for f (z) = z 2 − 1 as follows: In order to obtain a simple form of J f (s, u, v), we needed to closely inspect how it is connected with K f (s, u). In view of Relation (31), it is appropriate to select a form of K f (s, u) that reduces to a lower-order rational function in t. Such a lower-order rational weight function K f (s, u) would eventually lead us to obtaining a simplified J f (s, u, v). When applying to f (z) = z 2 − 1, we find K f (s, u) with t = z 2 as shown below: where β should be selected in such a way that the order of rational function K f (s, u) is minimized. For such minimization, we first conveniently let Since K 2 (1) = 256 guarantees that K 2 does not have a factor (t − 1) 4 for any value of β, we need to check if K 1 and K 2 have common factors for some values of β reducing the order of rational function K f . By eliminating β between K 1 and K 2 , we find (1 + 3t) 4 = 0, from which (1 + 3t) j = 0 with some j ∈ {1, 2, 3, 4} is found to be a common factor yielding a unique β = 2 in view of the fact that . Hence, with this β = 2 employed, (44) reduces to a desired simple rational weight function: Using the two selected weight functions Q f , K f ( with β = 2), we continue to determine coefficients q i , r i of J f yielding a simple governing equation for extraneous fixed points of the proposed methods when f (z) = z 2 − 1 is applied. As a result of tedious algebraic operations reflecting the 25 constraints (with possible rank deficiency) given by (30) and (31), we find only 21 effective relations, as follows: The two relations J 300 = −4 + 2Q 2 + Q 3 and J 400 = K 40 equivalently represent r 3 = q 3 , while three relations, J 500 = K 50 , J 600 = K 60 , and J 700 = K 70 are identically satisfied and give us no valuable information at all.
In what follows, we classify various subcases based on interesting selections of the 43 remaining free parameters. For each subcase, the concrete forms of J f are shown without displaying weight functions Q f and K f since they remain the same as in (42) with β = 2.

Remark 1.
The above Case 3E represents a method obtained with a different approach by Sharma et. al [14].
For ease of analysis for interesting subcases of Case 3F, we first impose further constraints: and then we seek parameter relationships yielding purely imaginary extraneous fixed points of the proposed family of methods when f (z) = z 2 − 1 is applied.
It still remains for us to select values of λ for the desired purely imaginary extraneous fixed points. To this end, we need to determine the conditions for negative roots of the cubic equation Φ(t) = 0 whose discriminant should be nonnegative and coefficients should have the same sign. (See Lemma 4.1 in Reference [23]). As a result, the values of λ should satisfy the following constraints: As a final case, we consider subcase Case 3G, leading us to purely imaginary extraneous fixed points. To easily proceed with our investigation, we initially impose three more constraints than Case 3F as follows: We simply let G(t) have a factor (1 + 3t) by taking the effect of 1 + sQ f (s) = 1+3t 2(1+t) into account, when f (z) = z 2 − 1 is applied. This case will give us the governing equation H(z) on the extraneous fixed points as follows: where A is a constant factor; G(t) and W (t) are 18-degree and 16-degree polynomials respectively with their coefficients containing free parameters. The two polynomials G(t) and W (t) shall reduce to 9-degree polynomial G 9 (t) and 7-degree polynomial W 7 (t), respectively, after lowering their degrees by annihilating their relevant coefficients gradually. As a result of this process of factorization and annihilation, we obtain a set of relations among the desired coefficients with 9 additional free parameters set to zero as follows: , r 1 = q 1 , r 2 = q 2 , r 3 = q 3 , r 4 = q 4 , r 5 = q 5 , r 6 = q 6 − 1, r 7 = q 7 − q 1 − 2, r 8 = q 8 + q 3 2 , r 9 = q 9 , r 11  , r 21 = q 21 2 , r 22 = −1, r 23 = −q 1 , r 24 = −q 2 , where λ = r 28 . In addition, the corresponding J f is given by: where J (s, u) is given by (42). With the coefficients given by (62), we also find H(z), G 9 (t), W 7 (t) as follows: In Proposition 1 of Section 4, it is shown that G 9 (t) and W 7 (t) have such factorizations as well as a common factor γ(t; λ). Consequently, after cancelling out the common factor, H(z) reduces to: Holding true regardless of λ. Hence, it is regretfully infeasible to directly use γ(t; λ) for obtaining the desired purely imaginary extraneous fixed points with a possible λ. Nevertheless, it is interesting to observe that the roots of the right-hand side of (65)  These subcases further simplify J f with integer coefficients q i , r i of J f with r 28 = 0, (q 1 = r 1 = r 23 = 0), (q 2 = r 2 = r 24 = 0), (q 3 = r 3 = r 25 = 0), (q 4 = r 4 = 0, r 26 = −1), (q 5 = r 5 = 0), (q 9 = r 9 = 0), (q 15 = r 15 = 0), (q 21 = r 21 = 0) in order.
In the next section, we investigate how to select appropriate free parameters giving purely imaginary extraneous fixed points.

Extraneous Fixed Points and Their Dynamics
The dynamics behind the extraneous fixed points [20] of iterative map (5) have been investigated by many authors with the aid of the relevant basins of attraction. Such dynamics were studied by Stewart [ [38][39][40], and Scott et al. [41].
To find a root α of f (x) under consideration, we usually locate a fixed point ξ of the iterative map: where R f is the iteration function associated with f . Usually, R f is expressed with weight function . Thus the zeros of H f are other fixed points ξ = α called extraneous fixed points of R f . The presence of extraneous fixed points may induce attractive, indifferent, or repulsive, and other periodic or chaotic orbits influencing the underlying dynamics of R f . The dynamics behind the extraneous fixed points motivates the current analysis. To make our analysis more feasible, we rewrite Iterative Map (66) in a more specific form: where H f (x n ) = 1 + sQ f (s) + suK f (s, u) + suvJ f (s, u, v) plays a role of a weight function in the classical Newton's method. It is clear that α is a fixed point of R f . Points ξ = α for which H f (ξ) = 0 are extraneous fixed points of R f . The influence of extraneous fixed points on the convergence behavior of the iterative dynamical system was extensively demonstrated for simple zeros via König functions and Schröder functions [20] with applications to a family of For ease of dynamics behind the extraneous fixed points of Iterative Maps (67), we select a simple member f (z) = (z 2 − 1). By a similar approach made by Chun et al. [42,43] and Neta et al. [38,40,44], we are able to construct (67). Applying f (z) = (z 2 − 1) to H f , we construct a rational function H(z) with t = z 2 in the form: where both D(t) and N (t) are coprime polynomial functions of t. The underlying dynamics of Iterative Map (67) can be favorably investigated on the Riemann sphere [45] with possible fixed points "0(zero)" and "∞". As can be seen in Section 5, the relevant dynamics will be illustrated in a 6 × 6 square region centered at the origin. Indeed, the roots t of N (t) provide the extraneous fixed points ξ of R f in Map (67) by the relation: (69)
Proposition 2. Let f (z) = z 2 − 1. Then the extraneous fixed points ξ for Cases 1-3 discussed in Section 3 are all found to be indifferent.
Proof. All subcases of Cases 1 and 2 have the same procedure for their stability proofs. Hence, it suffices to show the relevant proof for typical Subcases 1A and 2A.

(a)
Case 1A: H f (z) = By direct computation of R f (z) with f (z) = z 2 − 1, we write it as with t = z 2 : (1 + t)N 1 288,230,376,151,711,744t 43 , from which we find R f (ξ) = 1 due to the fact that N 1 = 0 at an extraneous fixed point ξ.

Remark 2.
Although not described here in detail due to limited space, by means of a similar proof as shown in Proposition 2, extraneous fixed points ξ for methods KT16, MBM (with G(u, s) = ) and SAK (Case 3E) were also all found to be indifferent.
If f (z) = p(z) is a generic polynomial other than z 2 − 1, then required dynamical analysis may not be feasible due to the increased algebraic complexity. Despite that, we explore such dynamics by means of Iterative Map (67) applied to f (z) = p(z), which is denoted by R p as follows: Basins of attraction for various polynomials are illustrated in Section 5 to observe the complicated dynamics behind the fixed points or the extraneous fixed points. In order to better describe the numerical and dynamical aspects of Iterative Maps (67), the letter W was conveniently prefixed to each case number to designate the relevant map in Table 1. The highlighted methods are investigated for their dynamics in Section 5.

Results and Discussion on Numerical and Dynamical Aspects
For various test functions, we begin by numerical aspects of (5), as well as three existing methods, , β = 0), and SAK; then we explore the dynamics underlying extraneous fixed points based on Iterative Maps (71) from the illustrated relevant attractor basins. Numerical experiments have been implemented for all listed methods in Table 1. Computational experiments on dynamical aspects have carried out with only 17 members of (5) and three methods KT16, MBM, and SAK. For each map, numerical experiments have strongly confirmed the desired convergence properties.
Throughout the computational experiments with the aid of Mathematica, $MinPrecision = 400 has been assigned to maintain 400 digits of minimum number of precision. If α is not exact, then it should be given by an approximate value with more precision digits than $MinPrecision.
The value of α is approximately given with 450 precision digits unless exact. Limited paper space allows us to list x n and α with up to 15 significant digits. We set error bound = 1 2 × 10 −360 to meet |x n − α| < .
Methods W3G1, W2D, and W1C successfully located desired zeros of test functions F 1 − F 3 : Ensured in Table 2 is sixteenth-order convergence. The computational asymptotic error constant |e n |/|e n−1 | 16 is in agreement with η = lim n→∞ |e n |/|e n−1 | 16 up to 10 significant digits. The computational convergence order p n = log |e n /η|/log |e n−1 | well approaches 16. Additional test functions in Table 3 confirm the convergence of scheme (5). The nth iterate errors |x n − α| are listed in Table 4 for comparison among methods W1A-W3G9 and three methods KT16, MBM, SAK. Bold-face numbers indicate the least errors for the selected test functions. In the current experiments, MBM has slightly better convergence for f 4 , f 6 , W3E or SAK for f 2 , f 3 , W1D for f 5 as well as W1C for f 1 and f 7 . According to the definition of the asymptotic error constant 16 , the convergence is dependent on iterative map R f (x n ), f (x), x 0 , α and the weight functions Q f , K f and J f . Consequently, it is hard to believe that one method always achieves better convergence than the others for any test functions.
The proposed Family of Methods (5) has efficiency index EI [11], which is 16 1/5 ≈ 1.741101 and larger than that of Newton's method. In general, the local convergence of iterative methods (71) is guaranteed with good initial values x 0 that are close to α. Selection of good initial values is a difficult task, depending on precision digits, error bound, and the given function f (x).
The influence of initial values x 0 on the global convergence is effectively described by means of a basin of attraction that is the set of initial values leading to long-time behavior approaching the attractors under the iterative action of R f . Basins of attraction contain information about the region of convergence. A method occupying a larger region of convergence is likely to be a more robust method. A quantitative analysis undoubtedly measures the size of region of convergence.
The basins of attraction, as well as the relevant statistical data, are constructed in a similar manner shown in the work of Geum-Kim-Neta [22].
Owing to the limited space, in Table 1, we explore the dynamics of 17 maps, W1A, W1C, W2A, W2D, W3A, W3C, W3F2, W3F3, W3G1, W3G2, W3G3, W3G4, W3G5, W3G6, W3G7, W3G8, W3G9, and three existing methods KT16, MBM, SAK, with applications to p k (z), (1 ≤ k ≤ 6) and one nonpolynomial equation through the following seven examples. In each example, we have shown dynamical planes for the convergence behavior of iterative maps x n+1 = R f (x n ) (67) with f (z) = p k (z) by illustrating the relevant basins of attraction through Figures 1-6. Example 1. We have taken p 1 (z) as a quadratic polynomial with all real roots: The roots are obviously ±1. Note that the extraneous fixed points are computed based on this example. Clearly the best methods have basins separated by the imaginary axis. Basins of attraction for W1A-W3G9, KT16, MBM with β = 0 and SAK are given in Figure 1. Methods W3A, W3G1-W3G8 show basins separated by the imaginary axis. Consulting Tables 5-7, we find that the methods W3G7-W3G9 use the least number of iterations per point on average (ANIP), followed closely by W3C, W3G1, W3G3-W3G6. The fastest method is W3A. The methods KT16, MBM, and SAK have the least number of black points. Methods W1A and W1C have the highest number of black points. Notice that these methods use polynomial weight functions.
(1) W1A (2) W1C (3) W2A (4) W2D  Example 2. We have taken p 2 (z) as a cubic polynomial: Basins of attraction are given in Figure 2. Notice that the basins for W1A and W1C have many black points and therefore will not be considered in the rest of the examples. Besides that, in view of a close inspection that the basins for W3G3-W3G6 have similarities to the other remaining members of the listed W3G-family, we will omit them in the rest of the examples. Consulting Tables 5-7, we find that the method with the fewest ANIP is MBM with 2.22 iteration. All the others require between 2.71 and 12.12. In terms of CPU timein seconds, the fastest is W3A (472.715 s) and the slowest is W1C (2399.093 s). The methods W1C and W1A have the most black points (59,904 and 58,910, respectively) and SAK has the least number (6 points). We will not consider W1A and W1C any further.

Example 3.
We have taken p 3 (z) as another cubic polynomial: All roots were easily found to be real. The basins for this example are plotted in Figure 3. The basins for W2A, W2D, KT16, and SAK are too chaotic. Based on Table 5 we see that W3C has the lowest ANIP followed closely by MBM. The fastest method is again W3A (406.102 s) and the slowest is W3G8 (2041.663 seconds). The methods KT16, MBM, and SAK have no black points, and the rest have between 58 and 204 black points.
The basins are given in Figure 4. We now see that W2A is the worst, followed by W2D. The best are those with smaller lobes along the diagonals. In terms of ANIP, W3C is the best (2.58), followed by MBM (2.64), and the worst is W2A (25.22). The fastest is again W3A (497.284 s), followed by W3C (542.509 s), and the slowest is W2A (5117.503 s). Most of the methods have between 1201 and 1889 black points with the worst being W2A with 218,849 points and W2D with 10,785 black points. We remove W2A from further consideration.  .

Example 5.
We have taken p 5 (z) as a quintic polynomial: The basins for the best methods left are plotted in Figure 5. The worst are W3A, W2D, and SAK. In terms of ANIP, the best is W3C (2.64), followed closely by MBM (2.76), and the worst are W3A (7.91) and SAK (5.40). The fastest is W3C using 615.892 s, followed by KT16 using 888.784 s, and the slowest was W3G1 (2409.685 s). SAK has 18 black points, but the basins are chaotic. The highest number of black points is for W3A (49,176), preceded by W2D with 13,828 black points. We remove W2D from further consideration. We now average all these results across the seven examples to try and pick the best method. MBM had the lowest ANIP (2.34), followed by W3C with 2.88. The fastest method was W3C (695.931 s), followed by W3A (847.596 s). MBM has the lowest number of black points on average (705), followed by KT16 (1101 black points).
Based on these seven examples we see that MBM and W3C have three examples with the lowest ANIP, W3G7, W3G8, and W3G9 each with one example. W3A is the fastest in five examples and W3C in two examples. Thus, we recommend W3C and W3G7, since W3C is in the top four and W3G7 in the top five in all three categories. MBM was at the top in only two categories, KT16 is in top three in two categories, and W3F3 and W3G9 were in the top six in two categories.

Conclusions
Both numerical and dynamical aspects of Iterative Maps (5) support the main theorem well through a number of test equations and examples. The W3F and W3G methods were observed to occupy relatively slower CPU time due to the intensive use of rational coefficients for weight function J f . If less number of the rational coefficients were employed, it would take less CPU time to build the relevant basins of attraction.
The proposed Family of Methods (5) favorably cover most of optimal sixteenth-order simple-root finders with certain weight functions developed (or to be developed), since they employ fairly generic weight functions. The dynamics behind the purely imaginary extraneous fixed points will choose best members of the family with improved convergence behavior. However, due to the order of convergence is rather high, the required algebra encounters difficulty resolving its increased complexity. The current work is limited to univariate nonlinear equations; its extension to multivariate ones becomes another task. In future work, as a follow-up study, we will not only extend Case 3 with other combinations of simple coefficients q i , r i , but also investigate different types of weight functions possessing less number of rational coefficients to obtain purely imaginary extraneous fixed points, and strengthen the desired computational as well as dynamical behavior.