EUROPEAN JOURNAL OF PURE AND APPLIED MATHEMATICS Vol. 16, No. 4, 2023, 2775-2785 ISSN 1307-5543 – ejpam.com Published by New York Business Global Study on the dynamical analysis of a family of third-order multiple zero finders Young Hee Geum Department of Mathematics, Dankook University, Cheonan, Korea 330-714 Abstract. We study the complex dynamics by analyzing the dynamical planes associated with the family of third order multiple root finders and the parameter spaces related to the free critical points. The conjugacy maps with the theoretical results of dynamical analysis for the iterative schemes are investigated. In addition, the various experiments are implemented to draw the dy- namics of the cubic-order schemes as well as existing methods. 2020 Mathematics Subject Classifications: 65H05, 65H99 Key Words and Phrases: Conjugacy map, basins of attraction, multiple root, nonlinear equa- tion, parameter space, Blood Rheology model 1. Introduction Many real-world problems involve finding solutions [1, 4, 17, 19] to equations. These equations can represent physical, biological, economic, and other systems. The root-finding problem[3, 10, 14] is a fundamental and crucial concept in mathematics, engineering, computer science, and various scientific disciplines. Root-finding algorithms[3] help us find the points where these equations cross zero, which often correspond to critical points or solutions of the underlying problem. The modified Newton scheme is the effective method to compute multiple roots given by xn+1 = xn −m f(x) f ′(x) . The scholars [2, 5, 7, 16, 20] are studying the dynamical problems of finding multiple zeros of nonlinear equations. The dynamical behavior of nonlinear equations[15] are found in many fields, such as artificial intelligence, engineering, medicine, and satellites. The following third-order method[11] given by{ xn+1 = xn − mf(yn) tmf ′(xn) , yn = xn −m(1− t) f(xn)f ′(xn) , (1) DOI: https://doi.org/10.29020/nybg.ejpam.v16i4.4986 Email address: conpana@empal.com (Y.H. Geum) https://www.ejpam.com 2775 © 2023 EJPAM All rights reserved. Y. H. Geum / Eur. J. Pure Appl. Math, 16 (4) (2023), 2775-2785 2776 is pursued its dynamics by constructing the conjugacy map. Osada’s third-order method[16] is written as xn+1 = xn − 1 2 m(m+ 1) f(xn) f ′(xn) + 1 2 (m− 1)2 f ′(xn) f ′′(xn) . (2) Dong’s third-order method [7] is given by{ xn+1 = zn − (1− 1√ m )1−m f(zn) f ′(zn) , zn = xn − √ m f(xn) f ′(xn) . (3) Let g : X → X and h : Y → Y be two analytic functions. The functions g and h are said to be topologically conjugate if there exists a homeomorphism k : X → Y such that k ◦ g = h◦k, where ◦ denotes function composition. Then the map k is called a conjugacy [4]. Then, h = k ◦ g ◦ k−1 and hn = k ◦ gn ◦ k−1 . If g is topologically conjugate to h via k and ν is a fixed point of h, then k−1(ν) is a fixed point of g. If g and h are invertible, then the topological conjugacy k maps an orbit of g onto an orbit of h and the order of points is preserved. In this paper, we study the dynamics of the parameter space of the following cubic order iterative method developed in [9] for finding multiple roots using the conjugacy map.{ yn = xn −m(1− t) f(xn)f ′(xn) , xn+1 = xn − f(yn)+(m−tm)f(xn) f ′(xn) . (4) The rest of this work is organized as follows. In section 2, we construct the conjugacy map and the stability surfaces are shown. Section 3 describes the complex dynamical analysis including the parameter spaces and the basins of attraction. In the last section, we describe the future work by investigating the dynamical analysis from the diverse viewpoints. 2. Dynamical Analysis A nonlinear equation (4) is reformed as a discrete dynamical system: xn+1 = If (xn), (5) where If is the iterative method. We have a complex discrete dynamical system zn+1 = If (zn) = zn − f(yn) + γf(zn) f ′(zn) , (6) where yn = zn − µh(zn), µ = m(1− t), h(zn) = f(zn) f ′(zn) and γ = m− tm. For M(z) = z−A z−B and f(z) = (z −A)m(z −B)m with m ∈ N [6], we have J(z, t) =M ◦ If ◦M−1(z) = z(r1r2 − r3µ1) zr1r2 − r3µ2 , (7) Y. H. Geum / Eur. J. Pure Appl. Math, 16 (4) (2023), 2775-2785 2777 where r1 = (t+ z)m, r2 = (1 + tz)m, r3 = (1 + z)2m, µ1 = tm +mz and µ2 = m+ tmz and A, B ∈ C ∪ {∞}, A ̸= B. Then If (z, t) is conjugate to J(z, t). Using the inverse of M(z), that is M−1(z) = Bz−A z−1 , we have the following form: J(z, t) = { (−2+t−z)z2(t+z) (−1+(−2+t)z)(1+tz) , m=1, z(t2(1+z)4+2z(1+z)4−(t+z)2(1+tz)2) 2(1+z)4+t2z(1+z)4−z(t+z)2(1+tz)2 , m=2. (8) We find out that z = 0 and z = ∞ are the fixed points of the conjugate map J(z, t), regardless of t-values. By a lengthy and accurate computation [18], we have that z = 1 is a strange fixed point of J which is not a zero of f(z) = (z − A)m(z − B)m, regardless of t-values. We will find the fixed points of the iterative scheme J(z, t). Let ϕ(z, t) = J(z, t) − z, whose roots are the fixed points of J . ϕ(z, t) = −z(z − 1)(r1r2 − r3µ3) zr1r2 − r3µ3 , (9) with µ3 = tm −m. We find that z = 0 and z = 1 are the roots of ϕ. To study the dynamical analysis behind iterative map (1) applied to a quadratic poly- nomial raised to the power ofm, f(z) = (z−A)m(z−B)m, we will consider the fixed points of J and their stability. For m ∈ {1, 2}, we have the explicit form of ϕ(z, t) satisfying ϕ(z, t) = { z(z−1)ψ1(z) q1(z) , if m=1, z(z−1)ψ2(z) q2(z) , if m=2, (10) where ψ1(z) = −(1 + (3− 2t+ t2)z + z2), q1(z) = −1− 2z + (−2t+ t2)z2, ψ2(z) = 2 + (8 + 2t− 4t2 + 2t3)z + (13− 2t2 + t4)z2 + (8 + 2t− 4t2 + 2t3)z3 + z4, q2(z) = 2 + 8z + (12− 2t+ 4t2 − 2t3)z2 + (7 + 2t2 − t4)z3 + (2− 2t+ 4t2 − 2t3)z4. Suppose ψ1(z) = 0 and q1(z) = 0, for some value of z. By eliminating t from two polynomials, we get W (z) = (z + 1)3. Substituting the roots of W (z) into ψ1(z) = 0 and q1(z) = 0, we get the relations for t. Solving them for t, we have t = 1. For t ̸= 1, ψ1( 1 z ) = z−2ψ1(z). If z ̸= 0 is a root of ψ1(z), then so is 1 z . In case of m = 2, we have W (z) = (z + 1)5 and ψ2( 1 z ) = z−4ψ2(z) through a similar process. We have that J ′(z, t) can be reformed a fraction as follows: J ′(z, t) = mr3(mz(−(−1 + t)2(−1 + z)2r −1+m m 1 r −1+m m 2 + 2r3) + (1 + z2)(−r1r2 + tmr3)) (zr1r2 − (m+ tmz)r3)2 . (11) We will draw the stability surfaces of the fixed points. The stability surfaces of the fixed points are shown by conical surfaces in Figures 1 - 2. Y. H. Geum / Eur. J. Pure Appl. Math, 16 (4) (2023), 2775-2785 2778 Computing the derivative of J , we obtain the following relations J ′(z, t) = { 2z(z+1)2Q1(z) w1(z)2 , if m=1, 2z(z+1)4Q2(z) w2(z)2 , if m=2, (12) where Q1(z) = −((−2 + t)t) + (3− 2t+ t2)z − (−2 + t)tz2, Q2(z) = −4(−1 + t− 2t2 + t3) + (13 + 8t− 10t2 + 8t3 − 3t4)z + 4(7− 4t+ 6t2 − 4t3 + t4)z2 + (13 + 8t− 10t2 + 8t3 − 3t4)z3 − 4(−1 + t− 2t2 + t3)z4, w1(z) = 1 + (4− t)z − (−2 + t)(2 + t)z2 + (−2 + t)2tz3, w2(z) = −2− 8z + 2(−6 + t− 2t2 + t3)z2 + (−7− 2t2 + t4)z3 + 2(−1 + t− 2t2 + t3)z4. The critical points of the numerical method refers to a point where the derivative of a function is zero, that is J ′(z, t) = 0. The points z = 0 and z = ∞ are critical points associated with (z−A)(z−B). The critical points that are not any roots of the polynomial (z −A)(z −B) are said to be free critical points. (a) |ℜ(t)| ≤ 50000, |ℑ(t)| ≤ 50000 (b) −1 ≤ ℜ(t) ≤ 1, |ℑ(t)| ≤ 1 (c) 0 ≤ |ℜ(t)| ≤ 2, − 1 ≤ |ℑ(t)| ≤ 1 (d) 0 ≤ |ℜ(t)| ≤ 2, − 1 ≤ |ℑ(t)| ≤ 1 Figure 1: Stability surfaces for m = 1. We describe the complex dynamical behavior from the viewpoint of parameter spaces and dynamical plane. The following theorems will be useful to confirm the properties of symmetry on dynamical planes and parameter spaces. LetD = {z ∈ C : an orbit of z under J(z, t) tends to a number ηd in C} as the dynam- ical plane and let P = {t ∈ C : a critical orbit of z under J(z, t) tends to a number ηp in C} as the parameter space [12, 13]. If the number ηd or ηp is a finite constant, there exist finite periods in the orbit. Otherwise, the orbits go to infinity or are periodic but bounded. We consider z = ∞ as a fixed point in Riemann sphere. Y. H. Geum / Eur. J. Pure Appl. Math, 16 (4) (2023), 2775-2785 2779 (a) |ℜ(t)| ≤ 5000, |ℑ(t)| ≤ 5000 (b) −2 ≤ ℜ(t) ≤ 3, |ℑ(t)| ≤ 3 (c) |ℜ(t)| ≤ 1, |ℑ(t)| ≤ 1 (d) |ℜ(t)| ≤ 5, |ℑ(t)| ≤ 5 (e) |ℜ(t)| ≤ 3, |ℑ(t)| ≤ 3 (f) |ℜ(t)| ≤ 3, |ℑ(t)| ≤ 3 Figure 2: Stability surfaces for m = 2. Theorem 1. Let z(t) be a point of iterative map J(z, t) where t ∈ R. Then the dynamical plane is symmetric with respect to its horizontal axis. Proof. Considering t = t, We have J(z, t): |J(z, t)| = |J(z, t)| = |J(z, t| = |J(z, t)|, which means that the magnitude of the orbit of z(t) is the same as that of the orbit of z(t). Then the dynamical plane related to map J(z, t) is symmetric with respect to its horizontal axis if t is real. Theorem 2. Let z(t) be a free critical point of J(z, t) depending on parameter t by a root of Q(z) = 0. Then the parameter space is symmetric with respect to its horizontal axis. Proof. For m = 2, we have z(t) is a zero of Q2(z) for a given t. Then z(t) is a root of Q(z) at t. For a free critical point z(t), we have the conjugated map J(z, t) from (8). Then |J(z, t)| = |J(z(t), t)| = |J(z(t), t)| = |J(z(t), t| = |J(z((t), t)|, which implies that the magnitude of the orbit of free critical point z is the same as that of the orbit of free critical point z at t. It states that the corresponding parameter space is symmetric with respect to its horizontal axis. We develop a systematic color palette in Table 1 to paint a point t ∈ P according to the orbital period of z of J(z, t) for t ∈ P. Then the point t is determined by the corresponding color Ck if t induces a k-periodic orbit with k ∈ N ∪ {0} under J(z, t). We use a tolerance of 10−6 after up to 1000 iterations to allow for desired k periodic convergence of an orbit[1, 18] associated with P. In this experiments, we draw the color Cq using the color palette in Table 1 Y. H. Geum / Eur. J. Pure Appl. Math, 16 (4) (2023), 2775-2785 2780 Table 1: Color palette for a q-periodic orbit with q ∈ N ∪ {0} q Cq q = 1 C1 =  magenta, for fixed point ∞ cyan, for fixed point 0 yellow, for fixed point 1 red, for other strange fixed point , 2 ≤ q ≤ 68 C2 = orange, C3 = light green, C4 = dark red, C5 = dark blue, C6 = dark green, C7 = dark yellow, C8 = floral white, C9 = light pink, C10 = khaki, C11 = dark orange, C12 = turquoise, C13 = lavender, C14 = thistle, C15 = plum, C16 = orchid, C17 = medium orchid, C18 = blue violet, C19 = dark orchid, C20 = purple, C21 = power blue, C22 = sky blue, C23 = deep sky blue, C24 = dodger blue, C25 = royal blue, C26 = medium spring green, C27 = spring green, C28 = medium sea green, C29 = sea green, C30 = forest green, C31 = olive drab, C32 = bisque, C33 = moccasin, C34 = light salmon, C35 = salmon, C36 = light coral, C37 = Indian red, C38 = brown, C39 = fire brick, C40 = peach puff, C41 = wheat, C42 = sandy brown, C43 = tomato, C44 = orange red, C45 = chocolate, C46 = pink, C47 = pale violet red, C48 = deep pink, C49 = violet red, C50 = gainsboro, C51 = light gray, C52 = dark gray, C53 = gray, C54 = charteruse, C55 = electric indigo, C56 = electric lime, C57 = lime, C58 = silver, C59 = teal, C60 = pale turquoise, C61 = sandy brown, C62 = honeydew, C63 = misty rose, C64 = lemon chiffon, C65 = lavender blush, C66 = gold, C67 = crimson, C68 = tan. q = 0∗ or q > 69 Cq = black. ∗: q = 0 : the orbit is non-periodic but bounded. 2. 2.1 2.2 2.3 2.4 2.5 2.6 2.7 2.8 2.9 3. -1. -0.9 -0.8 -0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1 0. 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 2.5 2.55 2.6 2.65 2.7 2.75 2.8 2.85 2.9 2.95 3. 3.05 3.1 3.15 3.2 3.25 -0.25 -0.2 -0.15 -0.1 -0.05 0. 0.05 0.1 0.15 0.2 (a) 2 ≤ ℜ(t) ≤ 3, |ℑ(t)| ≤ 1 (b) 2.5 ≤ ℜ(t) ≤ 3.25, |ℑ(t)| ≤ 0.8 3.01 3.02 3.03 3.04 3.05 -0.02 -0.01 0. 0.01 1.1 1.15 1.2 1.25 1.3 1.35 1.4 -0.15 -0.1 -0.05 2.77556×10 -17 0.05 (c) 3.01 ≤ ℜ(t) ≤ 3.05, |ℑ(t)| ≤ 0.2 (d) 1.1 ≤ ℜ(t) ≤ 1.4, |ℑ(t)| ≤ 0.15 Figure 3: Parameter space for m = 1 . In Figures 3–4, the related parameter spaces P for m = 1, 2 are displayed. A point t ∈ P is determined by the color palette in Table 1. In terms of numerical experiments, every Y. H. Geum / Eur. J. Pure Appl. Math, 16 (4) (2023), 2775-2785 2781 point of the parameter space P whose color is none of cyan(root z = A), magenta(root z = B), yellow or red is not a better choice of t. Let Pi denote the parameter space associated with branch cpi for 1 ≤ i ≤ 4. We figure out that n ∈ {2, 3, · · · } periodic orbit is budding at the main component(period-1 component) and 4-periodic component is budding at period-2 component. 2. 2.1 2.2 2.3 2.4 2.5 2.6 2.7 2.8 2.9 3. -0.5 -0.4 -0.3 -0.2 -0.1 0. 0.1 0.2 0.3 0.4 0.5 2.28 2.281 2.282 2.283 2.284 2.285 2.286 2.287 2.288 2.289 2.29 2.291 2.292 2.293 2.294 2.295 2.296 2.297 2.298 2.299 2.3 -0.01 -0.009 -0.008 -0.007 -0.006 -0.005 -0.004 -0.003 -0.002 -0.001 0. 0.001 0.002 0.003 0.004 0.005 0.006 0.007 0.008 0.009 (a) 2 ≤ ℜ(t)| ≤ 3, |ℑ(t)| ≤ 0.5 (b) 2.28 ≤ ℜ(t) ≤ 2.3, |ℑ(t)| ≤ 0.01 -1.5 -1.4 -1.3 -1.2 -1.1 -1. -0.9 -0.3 -0.2 -0.1 5.55112×10 -17 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1. 1.1 -0.4 -0.3 -0.2 -0.1 0. 0.1 0.2 0.3 (c) −1.5 ≤ ℜ(t) ≤ 0.9, |ℑ(t)| ≤ 0.3 (d) 0.3 ≤ ℜ(t) ≤ 1.1, |ℑ(t)| ≤ 0.4 Figure 4: Parameter space for m = 2 . -6 -4 -2 0 2 4 6 -6 -4 -2 0 2 -6 -4 -2 0 2 4 6 -6 -4 -2 0 2 -6 -4 -2 0 2 4 6 -6 -4 -2 0 2 -6 -4 -2 0 2 4 6 -6 -4 -2 0 2 4 (a) y1 (b) y2 (c) d1 (d) o1 Figure 5: Basins of attraction of Blood Rheology model We compare the proposed scheme(y1) in (4) with the third order method(y2) in (1), the Dong’s third-order method(d1) in (2) and the Osada’s scheme(o1) in (3) based on the dynamical analysis. We have taken the test functions having multiple roots with multiplicity m = 2, 3, 4, 5, 6 to draw the complex dynamics of the proposed method[9] y1, Y. H. Geum / Eur. J. Pure Appl. Math, 16 (4) (2023), 2775-2785 2782 Table 2: Test functions pδ method time con no div y1 78.937 359,555 7.63739 445 y2 75.516 359,554 7.4628 446 p1 d1 150.359 352,420 10.7719 7,580 o1 74.046 339,770 12.244 20,230 y1 126.313 359,998 5.97738 2 y2 43.64 359,910 6.00031 90 p2 d1 117.625 359,963 11.2723 37 o1 48.094 360,000 6.20214 0 y1 101.485 359,894 6.36682 106 y2 115.391 359,757 8.37467 243 p3 d1 165.921 355,194 11.243 4,806 o1 52.954 356,058 10.3829 3,942 y1 83.797 360,000 6.05845 0 y2 54.844 358.860 5.9645 1,140 p4 d1 544.218 5,263 20.4971 354,737 o1 46.219 360,000 6.17781 0 y1 104.625 359,937 6.84216 63 y2 62.781 358,867 5.97515 1133 p5 d1 607.906 1,057 8.1211 358,943 o1 676.015 188 12.6383 359,812 p1(z) = (z7 − 9)2, p2(z) = (z2 − 3z + 5)3 , p3(z) = (z5 − 6)4, p4(z) = (z2 + 3z + 7)5, p5(z) = (z2 + 3z + 5)6. Table 3: Blood Rheology model bδ method time con no div y1 486.844 359,920 9.54763 80 b1 y2 631.36 351,124 8.70651 8,876 d1 1014.52 355,602 17.8899 4,398 o1 233.89 359,910 9.85888 90 y2, Dong’s scheme and Osada’s method with the basins of attraction. In this statistical Table 2 for the basins of attraction, abbreviations time, con, no and div denote the value of CPU time for convergence, the value of total convergent points, the value of average iteration number for convergence and the value of divergent points. As the first example, we select the polynomial p1(z) = (z7 − 9)2 with roots z = −1.23319 ± 0.593873i,−0.304573 ± 1.33442i, 0.853394 ± 1.07012i, 1.36874 of multiplicity m = 2. The method y1 is better in view of con. As the next instance, the polynomial p2(z) = (z2 − 3z + 5)3 has the roots z = 1/2(3± i √ 11). The method y1 is better in view of no. We select p3(z) = (z5 − 6)4 with the roots z = −1.15768 ± 0.841103i, 0.442194 ± 1.36093i, 1.43097 and p4(z) = (z2 +3z +7)5, z = 1 2(−3± i √ 19). As the last test function, Y. H. Geum / Eur. J. Pure Appl. Math, 16 (4) (2023), 2775-2785 2783 -3 -2 -1 0 1 2 3 -3 -2 -1 0 1 2 3 -3 -2 -1 0 1 2 3 -3 -2 -1 0 1 2 3 -3 -2 -1 0 1 2 3 -3 -2 -1 0 1 2 -3 -2 -1 0 1 2 3 -3 -2 -1 0 1 2 (a) y1,m = 2 (b) y2,m = 2 (c) d1,m = 2 (d) o1,m = 2 -3 -2 -1 0 1 2 3 -3 -2 -1 0 1 2 -3 -2 -1 0 1 2 3 -3 -2 -1 0 1 2 -3 -2 -1 0 1 2 3 -3 -2 -1 0 1 0 100 200 300 400 500 600 0 100 200 300 400 500 600 (e) y1,m = 3 (f) y2,m = 3 (g) d1,m = 3 (h) o1,m = 3 -3 -2 -1 0 1 2 3 -3 -2 -1 0 1 0 100 200 300 400 500 600 0 100 200 300 400 500 600 -3 -2 -1 0 1 2 3 -3 -2 -1 0 1 2 (i) y1,m = 4 (j) y2,m = 4 (k) d1,m = 4 (l) o1,m = 4 -3 -2 -1 0 1 2 3 -3 -2 -1 0 1 2 -3 -2 -1 0 1 2 3 -3 -2 -1 0 1 -3 -2 -1 0 1 2 3 -3 -2 -1 0 1 2 -3 -2 -1 0 1 2 3 -3 -2 -1 0 1 (m) y1,m = 5 (n) y2,m = 5 (o) d1,m = 5 (p) o1,m = 5 -3 -2 -1 0 1 2 3 -3 -2 -1 0 1 -3 -2 -1 0 1 2 3 -3 -2 -1 0 1 2 -3 -2 -1 0 1 2 3 -3 -2 -1 0 1 2 3 -3 -2 -1 0 1 2 3 -3 -2 -1 0 1 (q) y1,m = 6 (r) y2,m = 6 (s) d1,m = 6 (t) o1,m = 6 Figure 6: Basins of attraction for test functions. we use p5(z) = (z2 + 3z + 5)6 with the zeros z = 1 2(−3± i √ 11). In Figure 6, the picture (d) has shown some black point and we find the considerable black point in (o), (s) and (t). The illustrative results are listed in Table 2 and Figure 6 . REFERENCES 2784 We choose the equation z = ( x 8 441 + 8x5 63 − 2857144357x4 50000000000 + 16x2 9 − 906122449x 250000000 + 3 10) 4 in Blood Rheology model[8] to carry out the experiments. The method y1 is better in view of con. 3. Discussion Using an Möbius conjugacy map applied to a polynomial of the form f(z) = (z − A)m(z − B)m with the multiplicity m, the complex dynamical analysis is investigated including the stability surfaces, the dynamical planes and the parameter spaces. We compare the proposed method with the existing iterative schemes. Based on the theoretical result, we experiment the test functions and draw the basins of attraction. The basins of attraction of this model are shown in Figure 5. We will improve the current study to accurately address the visualization of diverse iterative schemes. In addition, the parameter space and basins of attraction of the de- veloped multiple zero solver will be investigated in more detail. We will figure out the beautiful but complicated fractal generated by numerical iterative schemes from various perspectives. Acknowledgements The author is supported by Basic Science Research Program through the National Research Foundation of Korea funded by the Ministry of Education under the research grant(Project Number:NRF-2021R1A2C1012922 ) Conflicts of Interest : The author declares no conflict of interest. References [1] L Ahlfors. Complex Analysis. McGraw-Hill Book Inc., 1979. [2] S Amat, S Busquier, and S Plaza. Iterative root-finding methods. unpublished manuscript., 2004. [3] I Argyros and A Magre nan. On the convergence of an optimal fourth-order family of methods and its dynamics. Appl. Math. Comput., 252:336–346, 2015. [4] A Beardon. Iteration of Rational Functions. Springer-Verlag, New York, 1991. [5] R Behl, A Cordero, S Motsa, and J Torregrosa. On developing fourth-order optimal families of methods for multiple roots and their dynamics. Appl. Math. Comput., 265:520–532, 2015. [6] P Blanchard. The dynamics of newton’s method. Proceedings of Symposia in Appl. Math., 49:139–154, 1994. REFERENCES 2785 [7] C.C.Dong. A basic theorem of constructing an iterative formula of the higher order for computing multiple roots of an euqation. Mathematica Numerica Sinica, 11:445–450, 1982. [8] R.L. Fournier. Basic transport phenomena in biomedical engineering. CRC press., 2017. [9] Y.H Geum and Y.I Kim. A cubic-order variant of newton’s method for finding mul- tiple roots of nonlinear equations. Comput. Math. Appl., 62:1634–1640, 2011. [10] V Kanwar, S Bhatia, and M Kansal. New optimal class of higher-order methods for multiple roots, permitting f ′(xn) = 0. Appl. Math. Comput., 222(1):564–574, 2013. [11] H.J. Kim, A.K. Rathie, and Y.H. Geum. A family of optimal cubic-order multiple- root solvers and their dynamics. European Journal of Pure and Applied Mathematics, 16(3):1902–1912, 2023. [12] Y.I Kim and Y.H Geum. A triparametirc family of optimal fourth-order multiple-root finders and their dynamics. J. Appl. Math., 8436759:1–23, 2016. [13] Y.I Kim and Y.H Geum. Dynamical analysis via möbius conjugacy map on a uni- parametric family of optimal fourth-order multiple-zero solvers with rational weight functions. Discrete Dyn. Nature Soc., 7486125:1–19, 2018. [14] S Kumar, V Kanwar, and S Singh. On some modified families of multipoint iterative methods for multiple roots of nonlinear equations. Appl. Math. Comput., 218:7382– 7394, 2012. [15] S Li, X Liao, and L Cheng. A new fourth-order iterative method for finding multiple roots of nonlinear equations. Appl. Math. Comput., 215:1288–1292, 2009. [16] N Osada. An optimal multiple root-finding mehtod of order three. J. Comput. Appl. Math., 51:131–133, 1994. [17] A. M. Ostrowski. Solutions of Equations and System of Equations. Academic Press, New York, 1960. [18] A. M. Ostrowski. The Mathematica Book. Wolfram Media, 5 edition, 2003. [19] J. Traub. Iterative Methods for the Solution of Equations. Chelsea Publishing Com- pany, 1997. [20] J Wright, J Deane, M Bartucceli, and G Gentile. Basins of attraction in forced systems with time-varying dissipation. Commun. Nonlinear Sci. Numer. Simul., 29(1-3):72– 87, 2015.