Long-term orbit dynamics viewed through the yellow main component in the parameter space of a family of optimal fourth-order multiple-root finders

. An analysis based on an elementary theory of plane curves is pre-sented to locate bifurcation points from a main component in the parameter space of a family of optimal fourth-order multiple-root ﬁnders. We explore the basic dynamics of the iterative multiple-root ﬁnders under the M¨obius conjugacy map on the Riemann sphere. A linear stability theory on local bifurcations is developed from the viewpoint of an arbitrarily small perturbation about the ﬁxed point of the iterative map with a control parameter. Invariant conjugacy properties are established for the ﬁxed point and its multiplier. The parameter spaces and dynamical planes are investigated to analyze the underlying dynamics behind the iterative map. Numerical experiments support the theory of locating bifurcation points of satellite and primitive components in the parameter space.

1. Introduction. We can formulate a dynamical system by means of any fixed rule characterizing the time-dependence of an evolving point along with its position in the relevant state(phase)-space. It is then best described by a function whose domain and codomain respectively consist of time as an independent variable and state(phase)-space as a dependent variable. The independent variable time can be measured in terms of integers, real or complex numbers. Some examples of continuous dynamical systems can be seen in differential equations, whereas other examples of discrete dynamical systems can be found in difference equations. This analysis will be limited to a discrete dynamical system whose governing equation represents a difference equation formed by an iterative method: where Φ f is an iteration function or fixed point operator [10]. Such an iteration function can be found in finding the numerical roots for the equation f (x) = 0 encountered in many fields of applied sciences and engineering. The root-finding iterative numerical methods have been developed by many researchers [6,7,17,27].
The main task of this iterative method is to ensure its convergence to the desired root by tracing the sequence {x n+1 = Φ f (x n )} ∞ n=0 within a prescribed error bound. Regarding the iteration index n as the discrete time variable, we keep tracking of the such sequence to get useful information on the long-term behavior of discrete dynamical system (1). The root-finding process for the equation f (x) = 0 leads us to a sequence of images of initial guess x 0 under the action of Φ f : f (x 0 ), · · · , Φ n f (x 0 ), · · · , }, which is regarded as a discrete dynamical system. Iteration function Φ f is usually given by a meromorphic function whose dynamics requires a rather sophisticated analysis. According to Lemma 2.1 (to be shown later), it suffices to treat the dynamics of rational functions.
An iteration function Φ f will be represented by a family of optimal [24] fourthorder modified Newton-like existing multiple-root finders [17,18,32] and their dynamical behavior [3,4,8,9,10,13,14,16,19,25,30] behind the parameter spaces as well as dynamical planes will be pursued on the Riemann sphere C [5] via Möbius map z−b z−a , (a = b) [9] which is conjugated to (z − a) m (z − b) m : where L f : C → C is analytic [1] in a neighborhood of 0 and s is taken on the principal analytic branch of the m-th root. Evidently, Φ f (x n ) will be the right side of the last equation of (2). We are interested in the specific form of L f with a control parameter λ ∈ C: L f (s) = 1+(λ+1)s+(λ+2)s 2 1+λs . ( Subsequent analyses and developments will be treated in additional five sections. Described in Section 2 are preliminary studies on long-term behavior of a dynamical system defined on C. Section 3 fully discusses a linear stability theory based on the analysis of a small perturbation about the fixed point of the iterative map (2) under the Möbius conjugacy and classifies the local bifurcations into three types according to the location of the spectral radius of the conjugated iterative map along the unit circle. In Section 4, we describe some properties of the fixed and critical points related to the dynamics under the Möbius conjugacy map. Section 5 investigates a long-term dynamical behavior of the conjugated iterative map. Parameter spaces and dynamical planes are defined and extensively explored as well as visually being illustrated with a number of self-explanatory figures. Besides, we develop a theory on locating bifurcation points budding from another component (maximal connected set) [1] in the parameter space with both theoretical and numerical methods from a viewpoint of plane geometry involving simple 2-dimensional curves. A numerical algorithm is also established to find approximate bifurcation points under consideration. Finally in Section 6, we draw a overall conclusion and state the future work on the analysis of long-term dynamical behavior of other parameter-controlled iterative maps.
2. Preliminary studies. Lemma 2.1. Meromorphic functions on the Riemann sphere C can be represented by rational functions.
Proof. Note that a meromorphic function on C is analytic (holomorphic) except at a finite number of poles. Let f be meromorphic on C and z 1 , z 2 , · · · , z n ∈ C be its poles with their corresponding integer multiplicities d 1 , d 2 , · · · , d n . Then we construct which has no poles and perhaps has at most a pole at z = ∞. Hence p must reduce to a polynomial, from which we have f = p n i=1 (z−zi) d i as a rational function.
To discuss the orbit dynamics between two functions, we introduce the definition of conjugacy and a relevant theorem.
Then the isomorphism h is called a conjugacy [5].
. The converse is also true. Note that h(ξ) is also a fixed point of F under the topological conjugacy h. (b) Differentiating both sides of with respect to any z ∈ Y , we find its evaluation at ξ: from which F (h(ξ)) = G (ξ) holds true due to the fact that h (ξ) = 0 in view of diffeomorphism h. Remark 1. In view of Theorem 2.3, the conjugacy h preserves the dynamical behavior between F and G, possessing the invariance properties of the fixed point and its multiplier under the conjugacy. Furthermore, we find Hence the topological conjugacy h gives an isomorphism (one-to-one correspondence) between the orbits of F and G, stating that the two orbits behave similarly under h.
The following lemma is useful to the development of linear stability theory to be discussed in the next section.
Lemma 2.4. Let a n = (a (1) n , a (2) n , · · · , a (r) n ) ∈ C r with r ∈ N and n ∈ N ∪ {0}. Further let {a n } ∞ 0 be a sequence satisfying the following recursive relations among its components with a given ω ∈ C: where µ = 0 if r = 1, 1, if r = 1 and a 0 = 0 is assumed. Then the following hold: Proof. (i) For r = 1, (4) reduces to a single scalar equation of the form: a n . After omitting the superscript (1), we find the relation |a n | = |ω| n |a 0 |, which easily proves the corresponding assertion due to a 0 = 0.
(ii) For r ≥ 2, (4) can be written as a vector equation of the form: where B is an r × r matrix given by: Hence we find the resulting relation: a n = B n a 0 , where B n is expressed as: whose derivation can be done without difficulty by mathematical induction on n ∈ N. We now proceed with the proof of (ii) as follows. Close inspection of a general upper off-diagonal entry in (6) leads us to the equation: Remark 2. The result of (i) implies that {a n } ∞ 0 is bounded if and only if |ω| = 1. According to the Bolzano-Weierstrass Theorem [1], the exists a convergent subsequence of {a n }, which will exhibit a complicated limit behavior including a possible bounded oscillation.
Let Ω ⊂ C and f : Ω → Ω be analytic. Further, let f have a fixed point ξ ∈ Ω with |f (ξ)| < 1. Then f has a unique fixed point ξ and the sequence {z n+1 = f (z n )} ∞ 0 converges to ξ, provided that any z 0 ∈ Ω is given. Proof. The property of being analytic at ξ allows us to write where η → 0 as z n → ξ. We now choose a sufficiently large positive integer N and a positive constant β such that |f (ξ)| + |η| = β < 1 whenever n ≥ N . Then for any k ∈ N ∪ {0}, from which we obtain lim k→∞ |z N +k+1 − ξ| = 0, i.e., lim n→∞ z n = ξ. Suppose that µ = ξ is another fixed point satisfying lim n→∞ z n = µ. Since a convergent sequence must have a unique limit, µ = ξ, contradicting the hypothesis. Hence the fixed point ξ is unique.
In Lemma 2.5, z 0 can be chosen as a critical point z * of f , due to the fact that f (z * ) = 0 implies z * ∈ Ω. As a result, Lemma 2.5 with z 0 = z * yields the following corollary.
Corollary 1. Let ξ be the fixed point described in Lemma 2.5. Then every critical orbit (which is an orbit of a critical point) of f tends to ξ.
The notions stated in the following definition [22,23,29] are useful to describe the qualitative behavior of a dynamical system. Definition 2.6. Let z 0 be in the domain of f . The orbit of z 0 is defined to be the sequence {f k (z 0 )} ∞ k=0 with f k as the k-fold composite map of f . Then we say the following: If z 0 is a period-k point (or k-periodic point), then the orbit of z 0 is called a k-periodic orbit (or k-cycle). (c) If z 0 is a point some iterate of which is periodic, i.e., if there exists an integer 1 ≤ ≤ k − 2 satisfying f k (z 0 ) = f k− (z 0 ), then z 0 is called an eventually periodic (or a pre-periodic) point. (d) If the the orbit of z 0 contains a subsequence converging to to a stable periodic point, then z 0 is called an asymptotically periodic point. (e) If z 0 is not of types (a),(c),(d), then z 0 is called an aperiodic (or a non-periodic) point. The orbit of a non-periodic point is said to be "non-periodic, stochastic or chaotic".
3. Linear stability theory and local bifurcations. We first let a discrete dynamical system F : C m × C → C m with m ∈ N be characterized by an iterative scheme: where λ ∈ C is a control parameter and z 0 ∈ C m is given. Under the assumption that ξ ∈ C m is a fixed point of F, we excite an arbitrarily small perturbation about ξ as described by where the initial perturbation δ 0 = 0 is arbitrary and small. For a fixed λ, Taylor expansion of (9) about ξ up to the first-order term in δ n gives us the recurrence relation: with S = S(ξ, λ) as an m × m Jacobian matrix evaluated at ξ. Consequently, we are allowed to discuss the linear stability about the fixed point ξ = F(ξ, λ) for a fixed λ by considering the iterative scheme: Transforming (12) by means of an m × m nonsingular matrix P , we obtain the where b n = P −1 δ n , and J = P −1 SP is the Jordan canonical form [28] of S with k Jordan blocks J 1 , J 2 , · · · , J k . A typical J i is given by an r i × r i upper triangular matrix with ω i 's as all diagonal entries and 1's as all super-diagonal entries as shown below: where ω i is an eigenvalue of S. Without loss of generality, let J i be chosen such that |ω i | = ρ(S) which is the spectral radius of S. For simplicity, we denote J i bỹ J , ω i by ω and r i by r. The limit behavior of the perturbation δ n = P b n will be best described by analyzing a subsystem of the form: whereb n ∈ C r consists of the corresponding r components of b n related withJ . After writing out component relations for (14) and identifyingb n andJ respectively with a n and B, we apply Lemma 2.4 in case a n = δ n to obtain the corollary: According to Corollary 2, the fixed point behavior or long-term orbit behavior of the iterative map F is stable if and only if ρ(S) < 1 (modulus of all eigenvalues of S < 1), and unstable if and only if ρ(S) ≥ 1 (modulus of some eigenvalues of S ≥ 1, with the equality holding when some eigenvalues are repeated). When all eigenvalues of S are simple, according to the Bolzano-Weierstrass Theorem, the result of (i) indicates that there exists a convergent subsequence of {δ n }, which will induce k-periodic orbits of F as well as its non-periodic bounded orbits, if and only if ρ(S) = 1. Since ρ(S) around the unit circle plays a significant role in stability analysis, we had better designate the unit circle as the stability unit circle shown in Figure 1 for further analysis.
The geometrical properties in the parameter space when ρ(S) = 1 would strongly influence the long-term orbit behavior of F, which often causes an abrupt qualitative change when an eigenvalue ω (or Floquet multiplier [21]) of S with the maximum modulus crosses a certain location ω * of the stability unit circle. This kind of qualitative change in the behavior of a dynamical system is called a bifurcation. We classify such bifurcations [26] in the space of control parameters into three types according to the location of ω * along the stability unit circle as follows: (1) The (cyclic) fold(saddle-node) bifurcation occurs when ω * = 1.
The location of the control parameter in the parameter space where qualitative change in a dynamical behavior occurs is called a bifurcation point. After solving the relation |ω| = ρ(S(ξ, λ) = 1 for λ in terms of ω, we can trace the control parameter λ as ω varies along the stability unit circle, from which the bifurcation point λ * in the parameter space can be given by λ(ω * ) based on the types of bifurcation mentioned above.   (3), we will discuss the relevant properties behind fixed and critical points, as the control parameter λ varies in the finite complex plane. [9]. For convenience of analysis, after applying M (z) to f (z) = ((z − a)(z − b)) m , we find the resulting J:

Fixed points and stability. Let
with Θ N and Θ D as polynomials whose coefficients are generally dependent upon parameters m, a, b, λ. The conjugacy M (z) enables us to find that all the coefficients of both Θ N and Θ D are free from parameters m, a, b. Accordingly, we can simplify (15) as In C, we allow special treatment of arithmetic operations 1 0 = ∞ and 1 ∞ = 0 as well as ∞ = ∞ being its conjugate. We state the following lemma to efficiently treat the fixed and critical points of J(z; λ).
holds for any λ ∈ C and any z ∈ C.
Proof. We find that J is conjugate to itself via conjugacy for any z ∈ C and λ ∈ C, completing the proof.
The underlying dynamics behind iterative map (16) will be initiated by investigating the fixed points of J and their stability. In view of Theorem 2.3-(a) and (b), we express: where T (z; λ) = 1 + z(5 + λ) + z 2 (10 + 3λ) + z 3 (5 + λ) + z 4 and q(z) = 1 + z(4 + λ) + z 2 (5 + 2λ). The numerator of the last equation of (17) immediately gives us the fact that z = 0 and z = ∞ are two λ-free fixed points, which give little impact on the dynamics since their orbits simply tend to themselves. Fixed points other than {0, ∞} are called strange fixed points [12] which are not related to the roots of p(z). Clearly z = 1 is a strange fixed point which may influence the relevant dynamics. To locate other strange fixed points dependent on λ, we seek the roots of T (z; λ) = 0 in (17) for given values of λ. The following lemma can be easily proved by means of simple computations.
holds for any any λ ∈ C and any z ∈ C.
Corollary 3. Let J be given by (16). If ξ ∈ C is any fixed point of J, then so is 1 ξ . Proof. Let X = Y = C, F = G = J in Theorem 2.3 and h(z) = 1 z . Then by Theorem 2.3-(a), the proof is complete.
We pay special attention to additional strange fixed points including z = 1 (corresponding to the original convergence to infinity in view of the fact that M −1 (1) = ∞ or M (∞) = 1). Direct computation enables us to locate the strange fixed points from the roots of J(z; λ) = z in terms of λ, whereas their stability is described in the following two Subsections 4.1.1 and 4.1.2. 4.1.1. Fixed points. It would be worthwhile to resolve the inquiries: (i) For what values of λ can T (z; λ) and q(z) possess common divisors? (ii) Is it possible for either T or q to contain a divisor (z − 1) ? (iii) What are the strange fixed points? Such inquiries would be resolved by the following theorem whose proof is elementary and hence omitted. (b) Due to the fact that T and q have possible factors of (z − 1) for some special λ-values, we find , if λ = −10/3.
which are the roots of T (z; λ) = 0 for any given λ, then so are 1 z , in view of Corollary 3. Indeed, these four λ-dependent fixed points will recover the special fixed points in (18).

4.1.2.
Stability of the fixed points. The derivative of J from (16) characterizes the stability of the fixed points given by: where Q(z; λ) = 10 + 4λ + z(20 + 14λ + 3λ 2 ) + 2z 2 (5 + 2λ). Observe that fixed points 0 and ∞, being related to the roots a and b of polynomial p(z) = (z − a) m (z − b) m , are super-attractive and become the critical points of J(z; λ) in view of (19). The multiplier J (z; λ) = 0 at z = ∞ is found from the analyticity in a neighborhood of ∞ on the Riemann sphere in view of Corollary 4 (to be shown later). It is not difficult to prove the following lemma by means of a simple algebraic substitution.
holds for any λ ∈ C and any z ∈ C.
Similar treatment done in Theorem 4.3 successfully leads us to the following theorem. (b) Due to the fact that Q(z; λ) and q(z) have possible factors of z 3 (z + 1) 2 for some special λ-values, we find (c) For any λ ∈ C, the desired stability of a λ-dependent fixed point can be achieved graphically or analytically by varying parameter λ with the aid of Corollary 4. Corollary 4. Let ξ ∈ C be any fixed point of J defined in (16). Then the following holds:  Proof. (a) With the aid of (19), we get J (1; λ) = 8(4+λ) 10+3λ without difficulty. We seek λ satisfying |J (1; λ)| = 1 to give the required relation for a parabolic point. It is then not complicated to prove the remaining part.  Figure 2 illustrates that S is the circular boundary between M and Γ with a radius 16 55 and a center at (− 226 55 , 0). We appropriately call S as the stability circle, by taking into account the fact that the fixed point z = 1 with λ inside or outside of S becomes repulsive or attractive, respectively. Yellow-shaded region M is a bounded component inside of S, whereas non-shaded region Γ is an unbounded component outside of S.
In view of Corollary 4, T (z; λ) can be factored out with four linear factors as described in the following lemma, whose proof is easily made with the help of Mathematica [31] symbolic capability.
To find the super-attractors of conjugated maps J(z; λ) for some values of λ, we need to simultaneously solve T (z; λ) = Q(z; λ) = 0 for z and λ. Eliminating λ in T (z; λ) = Q(z; λ) = 0 easily gives us an 8-degree polynomial in z. The result of Corollary 4 requires only 4 super-attractors for each λ. Indeed, by a numerical technique, we list only 4 pairs of super-attractors in the following lemma:

4.2.
Critical points and free critical points. The roots of J (z; λ) = 0 are said to be the critical points of iterative map (16). At first glance of (19), we find that z = −1 is a critical point for any λ ∈ C. Observe that z = 0 and z = ∞ are critical points respectively associated with the roots a and b of the polynomial (z −a) m (z −b) m . If the critical points are not related to any roots of the polynomial (z − a) m (z − b) m , then they are called free critical points [12]. Moreover, we are interested in exploring the dynamics of strange fixed points under iterative map J(z; λ) related to free critical points.

5.1.
Parameter spaces and dynamical planes. Our further investigation on some properties of the relevant complex dynamics of conjugated map J(z; λ) given by (16) requires notions of parameter spaces and dynamical planes. It is essential to effectively construct parameter spaces and dynamical planes for the analysis of long-term behavior of J(z; λ). A useful task of constructing parameter spaces is preferably to use a critical point and generate its orbit under the action of J(z; λ) which will induce attracting periodic orbits for each given λ according to Corollary 1. The orbit behavior of two critical points z and 1/z of J is best described in Corollary 7 and Remark 4.

Corollary 7.
Let z ∈ C and q ∈ N. If z is a q-periodic point of J, then so is 1 z . Proof. In view of Lemma 4.1 and Remark 1, we find J q (z;

Remark 4. Corollary 7 and Eqn. (22) claim relation
for any integer q ≥ 1 and λ. Let the orbit of ζ 1 approach a q-periodic point of J. Then the orbit of ζ 2 approaches a q-periodic point ζ 2 = 1 ζ1 due to Corollary 7. Therefore, the orbit of ζ 2 behaves quite similarly to that of ζ 1 as follows: (1) If the orbit of ζ 1 converges to a q-cycle, then so does the orbit of ζ 2 .
(2) If the orbit of ζ 1 is divergent but bounded, then so is the orbit of ζ 2 .
(3) If the orbit of ζ 1 converges to ∞, then the orbit of ζ 2 converges to 0; and vise versa.
According to Remark 4, we explore only one branch of the critical points ζ 1 for its typical orbit behavior. We introduce two terms called the parameter space (plane) P and dynamical plane D to describe the underlying dynamics of the iterative map J(z; λ): P = {λ ∈ C : an orbit of a free critical point z(λ) under J(z; λ) tends to σ p ∈ C}, D = {z ∈ C : an orbit of z(λ) under J(z; λ) for a given λ ∈ P tends to σ d ∈ C}.
There exist a finite periodic orbit when σ p or σ d reduces to a finite constant. If no such σ p or σ d exists, then the orbit is non-periodic but bounded or approaches infinity. By definition, D will represent a union of attractor basins. Note that z = ∞ is treated as a fixed point on C. For any λ selected in some component of P, one can find members of proposed family (2) exhibiting similar dynamical orbit behavior in D. The values of λ contained in an appropriate component of P play a role selecting best members of family (2) to locate the desired roots z = 0 or z = ∞ with acceptable numerical stability. Figure 6(a) illustrates a parameter space P for critical point ζ 1 (λ). A point λ ∈ P is painted by the coloring scheme in Table 1. Any point λ ∈ P which is neither cyan (root z = a) nor magenta (root z = b) is not a good choice of λ for convergence behavior. According to a simple computation with λ = 0, we find that the orbit of J(ζ 1 ; 0) approaches ∞ and the orbit of J(ζ 2 ; 0) approaches 0, since ∞ and 0 are the roots of f (z) = (z − a) m (z − b) m .
A closer examination confirms that space P for critical point ζ 1 consists of a union of simply connected components including components A m = {λ : critical orbit of z cr under the map J(z cr ; λ) tends to a fixed point m} with m ∈ {0, 1, ∞}. In view of black * : q = 0 implies a non-periodic but bounded orbit. These 84 colors are explicitly illustrated in Figure 5.
Remark 4, λ-parameter spaces for both ζ 1 and ζ 2 have same components and boundaries in a neighborhood of which different anomalies occur. It is clear that the component associated with fixed point ∞ shares its boundary with the component associated with fixed point 0. Best members of the family can be found with λ selected in proper components of P for critical point ζ 1 (λ). The dynamical plane D is a set of starting points whose orbits under J(z; λ) tend to a value in C for a selected λ ∈ P. While tracing an orbit of a λ-dependent critical point, we may encounter with components of P containing a λ-parameter for which the corresponding dynamical plane D exhibits attractors 0 or ∞.
Lemma 5.1. Let F : C 2 → C be defined by F (λ, z) = λ · P 1 (z) + P 2 (z), where P 1 (z) and P 2 (z) are complex polynomials with real coefficients. Suppose F (λ, z) = 0. Let z denote the complex conjugate of z. Then the following hold: See the proof of Lemma 3.5 in [20].
Theorem 5.2. Let z(λ) be a free critical point of J(z; λ) given by a root of Q(z; λ) = 0 in (19). Then parameter space P is symmetric about its horizontal axis.
If J q (z; λ) tends to σ p ∈ C for a given q ∈ N ∪ {0}, then by induction on q we find: which states that J q (z; λ) tends to σ p . Hence two points λ and λ are painted with same color C q , proving that P related to J(z; λ) is symmetric about its horizontal axis.
Theorem 5.3. For a given λ ∈ R, let z(λ) be a starting point of iterative map J(z; λ). Then dynamical plane D is symmetric about its horizontal axis.
Judging from Figures 6-8 of P and D, we find members of the λ-dependent family exhibiting delicate dynamical orbit behavior. Figure 6 displays P and its components containing λ-parameters for which q-periodic orbits are generated with q ∈ N. In Figures 6(b)-6(c), we can clearly see period-1 components (one in red, the other one in yellow) and period-2 components (in orange color). Typical qperiodic components are identified by arrow lines. Period-q components in Figure  6(b) are born along the boundary of the period-1 red component in the manner of Faray sequence [2] based on Schleicher's Algorithm [15]. Figure 7 illustrates some magnified q-periodic components of P with bifurcation points (to be described in Section 5.2). Figure 6(c) shows period-q components emerging from the boundary of the stability circle S. -3. -2. -1.
(a) Parameter space P In Figure 8, we observe D for members of the family (2) with regions of starting values z 0 whose orbits under J(z 0 ; λ) converge to some of the strange fixed points as well as other k-periodic cycles for k ∈ {2, 3, 4, 5, 6, 7, 8, 9, 14}   In view of the above-mentioned D, there exist regions of starting values z 0 whose orbit under J(z 0 ; λ) converges to period-q points for q ∈ {2, 3, · · · , }. Besides, each complement of such D contains of a union of basins of the fixed point in {0, 1, ∞}.

5.2.
Bifurcation points in parameter space P. It is said that a component H in P is hyperbolic if the orbit of z under J(z; λ) has a finite attracting cycle for each given λ ∈ H ; the word 'hyperbolic' comes from the meaning of 'hyperbolicity' [11] of a rational function, not from that of a fixed point of a dynamical system. We find that H is an open set with infinitely many connected components each of which is characterized by the period of the corresponding cycle. For ease of  [29] of H . If we view this kind of budding phenomenon from the side of W , then the root point ∈ W can be regarded as a bifurcation point of W where it splits into two components meaning the word bifurcation. To locate such bifurcation points λ ∈ P, we define a k-periodic component H k = {λ ∈ C : J(z; λ) has an attracting k-cycle} as a subset of H . For notational convenience in the subsequent discussion, we write J(z; λ) as J λ (z). Indeed, H k is characterized Figure 9. Typical geometries for primitive and satellite components as: It is worth to note that H 1 plays the role of a main component from which finiteperiodic components are born. A choice of λ-dependent critical points ζ(λ) found from Q(ζ; λ) would produce the parameter space P whose λ-value leads us to a longterm dynamical behavior with possible periodic, non-periodic or chaotic orbits. We further characterize H 1 with the fixed point ξ = 1 or λ-dependent strange fixed point ξ = ξ(λ) found from T (ξ; λ) = 0 in (17). The current analysis considers H 1 only related to strange fixed point ξ = 1. We name such H 1 as the yellow main component, which turns out to be exactly M described in Remark 3. For ease of subsequent analysis, let us call M the main component.

Exact bifurcation points.
Observing parameter space P in Figure 6, we should pay attention to various components of finite periods budding from the boundary S of M . We denote the finite boundary of M by ∂M and select a parameter λ ∈ ∂M . Let q ∈ N be given. If an attracting period-q component H q arises at λ, then the point λ is called the period-q bifurcation point. If M is itself a primitive component, then the bifurcation point λ turns to be a cusp point. We now consider two components M and H q (budding from M ) osculating at the common boundary point λ as shown in Figure 9. Then we have the following relations for λ ∈ ∂M ∂H q and a fixed point ξ ∈ C: with µ = J λ (ξ). Let ξ = |ξ|e iθ be a parametric representation for θ ∈ [0, 2π). We will show that µ q = 1 holds at λ where M and H q share the common tangent line.
We solve for λ from relation J λ (ξ(θ)) = ξ(θ) to obtain: where F (z) is a rational function in z in view of (23) and differentiable everywhere except at a finite number of poles. This λ will trace a curve in the complex plane as parameter θ varies. As a result, we find the derivative at the fixed point ξ: from which we have µ q = 1 using the relation J q λ (ξ) = ξ. Hence, Eqn. (25) reduces to: Since dλ dθ represents a tangent line in the direction of θ, Eqn. (27) describes the common tangent line of M and H q . The preceding discussion thus far enables us to introduce the notion of /q-bifurcation point or /q-root point.
Typical values of λ for 1 ≤ q ≤ 10 are given in Table 2. The type of 0/1bifurcation point is fold, while the type of 1/2-bifurcation point is flip. The type of all other /q-bifurcation points is of Neimark-Sacker. Note that 1/q-bifurcation points are positioned clockwise monotonically along the lower half-circle of S as q ≥ 2 increases. On the other hand, /q-bifurcation points are positioned counterclockwise monotonically along the full-circle of S as increases for a fixed q ≥ 2.

Numerical bifurcation points.
When the boundary of M is not exactly known, we resort to a numerical method to find local bifurcation points.
Let q, k ∈ N be given. Consider a period-qk satellite component H qk budding at λ from a period-q component W q . Applying the same argument as done in the preceding subsection with λ ∈ ∂W q ∂H qk and a q-periodic fixed point ξ(θ) ∈ C, we find: J q λ (ξ) = ξ, d dz J qk λ (z)| z=ξ = 1.
Since the first relation J q λ (ξ) = ξ in (30) defines a bivariate rational function in λ and ξ, we can solve for λ as an analytic function of ξ(θ), except at a finite number of poles. Hence, we express λ = F (J q λ (ξ(θ))) implicitly with F as a bivariate rational function of λ(θ) and ξ(θ). Since J q λ (ξ) = nq j=0 λ j P j (ξ)/ mq i=0 λ i Q i (ξ), with P j (ξ) and Q i (ξ) as polynomials in ξ for q-dependent nonnegative integers n q and m q , the relation J q λ (ξ) = ξ defines a polynomial in λ with its coefficients as polynomials of ξ, and hence λ can be solved as an analytic function of ξ except at the corresponding Table 2. Typical /q-bifurcation points λ for 1 ≤ q ≤ 10 and 0 ≤ ≤ 9   (2) with a parameter-controlled family of first-order rational weight functions. A study on the dynamical orbit behavior has been treated with the introduction of the parameter spaces and the dynamical planes. The boundary of the yellow main component has been found to be a circle along which infinitely many periodic satellite components are born. We have employed the linear stability theory based on a small perturbation about the fixed point and the elementary geometrical viewpoint through a parametric representation of the fixed point on the common tangent line of the main and a budding component in order to find the desired bifurcation points. The parameter spaces in Figure 6 display components containing λ-parameter for which the orbit under iterative scheme Φ f approaches a periodic orbit or chaotic [22] orbit. With λ selected in the components colored in magenta or cyan, Φ f stably converges to the desired root. However, for λ selected in other components, Φ f would exhibit unfavorable numerical behavior.
One should note that Theorem 5.5 and Algorithm 5.6 play important roles in locating exact or numerical bifurcation points. For higher-values of qk, one would encounter difficulties in finding such bifurcation points since the composition task of J qk λ (z) increases the degree of complexity. As part of our future study, we will pursue the exact boundary equation of the red main component shown in Figure  6(b) and locate the bifurcation points where period-q components are born along the boundary.