Design, Convergence and Stability Analysis of a New Parametric Class of Fourth-Order Optimal Iterative Schemes for Solving Nonlinear Equations

1. Introduction and Preliminary Concepts

Numerical methods play a fundamental role in scientific and engineering disciplines due to their ability to solve mathematical equations, enabling the analysis of physical phenomena and complex systems. These algorithms provide approximate solutions when exact solutions are challenging to obtain, while also offering computational efficiency in terms of time and computational resources by effectively performing complex calculations.

In the field of numerical analysis, there is a current trend of designing families of numerical methods that generalize existing approaches. Notable examples include King’s family [23] and that developed by Hueso et al. in [20]. Some of these methods incorporate weight functions in their design process (see, for example, [11,12,14]).

The efficiency index of an iterative method is defined by Ostrowski in [26] as $EI=p^{1/d}$, where p is the order of convergence and d is the number of functional evaluations per iteration. This concept is directly related with Kung-Traub’s Conjecture [24], that states the order of convergence of any iterative scheme cannot be greater than $2^{d-1}$ (called optimal order).

Several authors have generalized iterative schemes optimally according to Kung-Traub’s conjecture [24]. However, since there is no difference in the number of function evaluations and the order of the members within a family of iterative methods, it is necessary to conduct a stability study on some simple nonlinear functions, such as quadratic polynomials. This is because in practice, it has been observed that iterative schemes that are stable for such functions tend to perform better compared to other methods that exhibit some pathology when applied to more complicated functions. To this end, complex discrete dynamics are employed to analyze stability in functions like quadratic polynomials (see, for example, the work of Amat et al. in [1,2], Behl et al. in [3], Chicharro et al. in [9], Cordero et al. in [15,16,17]), Khirallah et al. [22], Moccari et al. in [25], among others.

The detailed description of the following concepts in complex dynamics can be found in [4,19]. Let $R:\widehat{\mathbb{C}}\to \widehat{\mathbb{C}}$ be a rational function on the Riemann sphere $\widehat{\mathbb{C}}=\mathbb{C}\cup$$\{\infty \}$, then we denote by $R(z)$ the operator $R(z)=\frac{P(z)}{Q(z)}$, where $P(z)$ and $Q(z)$ are complex polynomials with no common factors. Usually, R is obtained by applying an iterative scheme on a quadratic polynomial $p(z)=(z-a)(z-b)$.

Given $z_0\in \widehat{\mathbb{C}}$, we define the orbit of $z_0$ under R to be the sequence of points $z_0,z_1=R\left(z_0\right),z_2=R^2\left(z_0\right),\dots ,z_n=R^n\left(z_0\right),\dots$. Point $z_0$ is called the seed of the orbit. There are many different kind of orbits in a typical dynamical system. Undoubtedly, the most important kind of orbit is a fixed point, that is, a point $z_0$ that satisfies $R\left(z_0\right)=z_0$. Another important element is the periodic orbit or cycle. A point $z_0$ is n-periodic if $R^n\left(z_0\right)=z_0$ for some $n>0$, being $R^p(z_0)\ne z_0$ for any $p<n$. The least such n is called the prime period of the orbit.

Let us suppose that $z_0$ is a fixed point of R. Then $z_0$ is an attracting fixed point if $\left|R^{\prime }\left(z_0\right)\right|<1$. The point $z_0$ is a repelling fixed point if $\left|R^{\prime }\left(z_0\right)\right|>1$ and, finally, if $\left|R^{\prime }\left(z_0\right)\right|=1$, the fixed point is called neutral or indifferent.

Theorem 1.

[19] Attracting Fixed Point Theorem. Suppose $z_0$ is an attracting fixed point for R. Then there exists a domain D contained within the Riemann sphere such that $z_0\in D$ and in which the following condition is satisfied: if $z\in D$, then $R^n(z)\in D$ for all n and, moreover, $R^n(z)\to z_0$ as $n\to \infty$.

Theorem 2.

[19] Repelling Fixed Point Theorem. Suppose $z_0$ is a repelling fixed point for R. Then there exists a domain D contained within the Riemann sphere such that $z_0\in D$ and in which the following condition is satisfied: if $z\in D$ and $z\ne z_0$, then there is an integer $n>0$ such that $R^n(z)\notin D$.

Now, let us suppose $z_0$ is an attracting fixed point for R. The basin of attraction of $z_0$ is the set of all points whose orbits tend to $z_0$. The immediate basin of attraction of $z_0$ is the largest convex component containing $z_0$ that lies in the basin of attraction. The critical points of the operator are those $z^C$ that meet $R^{\prime }\left(z^C\right)=0$. A critical point is called free if it does not match with the roots of the polynomial. These two concepts are closely related with the next result.

Theorem 3.

[6,7] If R is a rational function, then the immediate basin of attraction of a periodic (or fixed) attracting point contains at least one critical point.

Fatou and Julia studied the iteration of rational maps $R:\widehat{\mathbb{C}}\to \widehat{\mathbb{C}}$ under the assumption that $deg(R)\ge 2$. They focused on a disjoint invariant decomposition of $\widehat{\mathbb{C}}$ into two sets. One of these sets is often called the Julia set. The other set we shall refer to it as the Fatou set [4]. The Fatou set is the union of all basins of attraction, and the Julia set is its boundary.

By Theorem 3, all the attracting behaviour of a rational function can be found by iterating the free critical points and classifying them by their asymptotical performance in the Fatou set. This is done by means of the parameter plane, in case the rational function R depends on a complex parameter $\gamma$. It is the graphical representation that provides information about the choice of one or another value of a parameter $\gamma$ within a family of iterative methods. This graphical representation directly relates each point in the complex plane to its corresponding parameter value that specifies each member of the family. Given a free critical point of R used as initial estimation, the parameter plane indicates which attracting periodic orbit or fixed point the orbit of the critical point converges to.

When all the free critical points are known, we calculate the parameter planes to determine which values of $\gamma$ relates the final state of all free critical points to any of the roots of $p(z)$. According to Theorem 3, this will ensure stability (on quadratic polynomials) of the associated elements of the class, as the only basins of attraction correspond with these roots.

Given a free critical point $cr(\gamma )$ as a seed of the rational function R, we define $D_{cr}$ as the set of values of $\gamma$ where the free critical point is in the basin of attraction of any of the roots of $p(z)$.

We also define the characteristic function of the set $D_{cr}$ as:

$$F_{D_{cr}}(\gamma )=\left\{\begin{array}{ccc}1& \mathrm{if}& \gamma \in D_{cr}\\ 0& \mathrm{if}& \gamma \notin D_{cr}.\end{array}\right.$$

Therefore, we consider the parameter plane by representing the characteristic function of $D_{cr}$, where the red color is assigned when $F_{D_{cr}}(\gamma )=1$ and black when $F_{D_{cr}}(\gamma )=0$.

Our interest is to determine the values of $\gamma$ in which the final state of the orbits of all free critical points is one of the polynomial roots. Since we construct a parameter plane for each free critical point, we proceed to construct the intersection of the parameter planes, that is called unified parameter plane [10]. In this plane, we graph $F_{{\bigcap }_{i=1}^n{D}_{c{r}_i}}(\gamma )$, where n is the number of free critical points of the rational operator. Each red-colored point corresponds to values of $\gamma$ for which all members of the family of iterative schemes are stable on quadratic polynomials.

One way to validate the information obtained through parameter planes is by using dynamical planes, which graphically represent the basins of attraction of a set of initial values for a specific member of the family, that is, for a value of $\gamma$. In other words, given a parameter value for which the associated method is stable, its dynamical plane is composed only of the basins of the polynomial roots. This indicates that, regardless of the initial estimate of the method, it will always be able to converge to one of the roots. On the other hand, if the parameter value is in the instability region on the parameter plane, in addition to the basins of the roots, other basins corresponding to attracting periodic or strange fixed points appear. Through the dynamical planes, we determine to which attracting fixed point the orbit of any initial estimate $z_0$ converges for a specific value of $\gamma$.

The basin of attraction for the root $z=a$ is represented in orange color, while the basin for $z=b$ is represented in blue. Additionally, we assign different colors such as green, red, etc., depending on the number of attracting fixed points associated with $\gamma$, and we use black to represent basins of periodic orbits. The decrease in brightness of each color is an indicator of the size of the orbit of each point, meaning that brighter colors represent points that require fewer iterations to reach the fixed or periodic point. The codes used to present the parameter and dynamical planes are defined by Chicharro et al. in [8].

Let us remark that all the information obtained depends on the quadratic polynomial $p(z)$ used to define the rational function R, by applying the class of iterative methods on it. A key tool to extend these results to any quadratic polynomial is the conjugation. Let f and g be two functions from the Riemann sphere to itself. An analytic conjugation between f and g is a diffeomorphism h of the Riemann sphere onto itself such that $h\circ f=g\circ h$, where h is called the conjugation. In [5], Blanchard introduced the conjugation map $\varphi (z)=\frac{z-a}{z-b}$, which is a Möbius transformation satisfying the following properties: $\varphi (\infty )=1$, $\varphi (a)=0$, and $\varphi (b)=\infty$ being a and b the roots of the arbitrary quadratic polynomial $p(z)=(z-a)(z-b)$.

Theorem 4.

Scaling Theorem [4,18] Let $f(z)$ be an analytic function and let $A(z)=\eta z+\beta$, with $\eta \ne 0$ an affine application. If $h(z)=\lambda (f\circ A)(z)$ with $\lambda \ne 0$, then the fixed point operator $R_f$ is analytically conjugate with $R_h$ through A, i.e. $\left(A\circ R_h\circ A^{-1}\right)(z)=R_f(z).$

If a family of iterative methods satisfies the Scaling Theorem, then the Möbius transformation can be applied. As demonstrated by Blanchard in [4], the asymptotic behavior of the system is qualitatively equivalent under conjugation. This allows us to perform dynamical analysis on an associated operator that does not depend on the roots a and b. Instead, the dynamical study is conducted based on their corresponding values in the new operator, which are $\varphi (a)=0$ and $\varphi (b)=\infty$. Additionally, the strange fixed point $z=1$ corresponds to the divergence of the original method. This makes the analysis of its stability particularly important.

In this investigation, using the weight function technique in Section 2, a parametric family of iterative methods of order four is designed. In Section 3, the rational function resulting from the fixed point operator applied on a quadratic polynomial $p(z)$ is analyzed. More specifically, we study the dynamics of the fixed points of this rational operator R, which allows us to determine its stability. Also, in Section 4, the numerical behavior of some of the most stable members on other non-polynomial functions is analyzed and compared with other known methods in the literature. Finally, in Section 5, some conclusions are presented, highlighting the most significant findings of the research.

2. Construction of a New Parametric Family of Iterative Methods

If $f(x)=0$ is a nonlinear equation, then an iterative method, under certain conditions, generates a sequence of real numbers that converges to a solution $\overline{x}$. However, the guarantee of convergence to a solution, the convergence rate, and the solution to which an iterative scheme converges all have a direct dependence on the initial estimate. The creation of a new class of iterative methods makes sense as long as there are members of the family that are competitive compared to existing schemes with good behavior. In our study, we present a new family of iterative methods using the weight function procedure, characterized by:

$$\begin{array}{cc}\hfill y_k& =x_k-\beta \frac{f\left(x_k\right)}{f^{\prime }\left(x_k\right)},\hfill \\ \hfill x_{k+1}& =x_k-H(\mu )\frac{f\left(x_k\right)}{f^{\prime }\left(x_k\right)},k=0,1,\dots ,\hfill \end{array}$$

(1)

where $\beta$ is an arbitrary parameter and H is a real function depending on $\mu (x)=\frac{f^{\prime }(x)}{f^{\prime }(y)}$. Now, we present the convergence result of this family, describing the conditions that function $H(\mu (x))$ must satisfy in order to achieve fourth-order convergence of (1).

Theorem 5.

Let $\overline{x}\in I$ a simple zero of a sufficiently differentiable function $f:I\subset \mathbb{R}⟶\mathbb{R}$ and $\mu (x)=\frac{f^{\prime }(x)}{f^{\prime }(y)}$. Let us assume an initial estimation $x_0$ close enough to a zero $\overline{x}$ of f. If a weight function $H\left(\mu (x_k)\right)$, $\mu (x_k)=\frac{f^{\prime }(x_k)}{f^{\prime }(y_k)}$, is chosen satisfying $H(1)=1$, $H^{\prime }(1)=\frac{3}{4}$, $H^{\prime \prime }(1)=\frac{3}{4}$ and $\left|H^{\prime \prime \prime }(1)\right|<+\infty$, then class (1) has fourth-order of convergence if and only if $\beta =\frac{2}{3}$, being the error equation is given by

$$e_{k+1}=\frac{1}{81}\left((117-32H^{\prime \prime \prime }(1))c_2^3-81c_3{c}_2+9c_4\right)e_k^4+O\left(e_k^5\right),$$

where $c_j=\frac{1}{j!}\frac{f^{(j)}(\overline{x})}{f^{\prime }(\overline{x})},\mathrm{for}j=2,3,\dots ,$ and $e_k=x_k-\overline{x}.$

Proof. 

By applying the Taylor expansion of $f(x_k)$ on $\overline{x}$, we have:

$$f(x_k)=f^{\prime }(\overline{x})\left[e_k+c_2{e}_k^2+c_3{e}_k^3+c_4{e}_k^4+O\left(e_k^5\right)\right],$$

(2)

and

$$f^{\prime }(x_k)=f^{\prime }(\overline{x})\left[1+2c_2{e}_k+3c_3{e}_k^2+4c_4{e}_k^3+O\left(e_k^4\right)\right],$$

(3)

dividing directly (2) by (3) we get

$$\begin{array}{c}\hfill \frac{f(x_k)}{f^{\prime }(x_k)}=e_k-c_2{e}_k^2+2\left(c_2^2-c_3\right)e_k^3+\left(-4c_2^3+7c_2{c}_3-3c_4\right)e_k^4+O\left(e_k^5\right).\end{array}$$

If in the first step of (1) we subtract $\overline{x}$ on both sides, then we have its error expression

$$\begin{array}{cc}\hfill y_k-\overline{x}& =e_k-\beta \frac{f(x_k)}{f^{\prime }(x_k)}\hfill \\ & =(1-\beta )e_k+\beta c_2{e}_k^2+2\beta \left(c_3-c_2^2\right)e_k^3+\beta \left(4c_2^3-7c_2{c}_c+3c_4\right)e_k^4+O\left(e_k^5\right).\hfill \end{array}$$

Then,

$$\begin{array}{c}\hfill f^{\prime }(y_k)=1-2(\beta -1)c_2{e}_k+B_1{e}_k^2+B_2{e}_k^3+O\left(e_k^4\right),\end{array}$$

(4)

being

$$\begin{array}{cc}\hfill & B_1=3{(\beta -1)}^2{c}_3+2\beta c_2^2,\hfill \\ \hfill & B_2=-4{(\beta -1)}^3{c}_4-4\beta c_2^3+2(5-3\beta )\beta c_2{c}_3.\hfill \end{array}$$

Dividing (3) by (4),

$$\begin{array}{c}\hfill \mu (x_k)=\frac{f^{\prime }(x_k)}{f^{\prime }(y_k)}=1+2\beta c_2{e}_k+D_1{e}_k^2+D_2{e}_k^3+O\left(e_k^4\right),\end{array}$$

where

$$\begin{array}{cc}\hfill & D_1=\left((4\beta -6)c_2^2-3(\beta -2)c_1{c}_3\right)\beta ,\hfill \\ \hfill & D_2=\left(2\left({\beta }^2-3\beta +2\right)c_2^3+\left(-3{\beta }^2+9\beta -7\right)c_2{c}_3+\left({\beta }^2-3\beta +3\right)c_4\right)4\beta .\hfill \end{array}$$

We expand $H\left(\mu (x_k)\right)$ around one because as k tends to infinity, variable $\mu$ approaches unity,

$$\begin{array}{cc}\hfill H\left(\mu (x_k)\right)=& H(1)+H^{\prime }(1){\nu }_k+\frac{H^{\prime \prime }(1)}{2!}{\nu }_k^2+\frac{H^{\prime \prime \prime }(1)}{3!}{\nu }_k^3+O\left({\nu }_k^4\right)\hfill \\ \hfill =& H(1)+2\beta H^{\prime }(1)c_2{e}_k+F_1{e}_k^2+\frac{2\beta }{3}(F_2+F_3)e_k^3+O\left(e_k^4\right),\hfill \end{array}$$

where

$$\begin{array}{cc}\hfill & {\nu }_k=\mu (x_k)-1,\hfill \\ \hfill & F_1=\left(2\left((2\beta -3)H^{\prime }(1)+\beta H^{\prime \prime }(1)\right)c_2^2-3(\beta -2)H^{\prime }(1)\right)\beta c_3,\hfill \\ \hfill & F_2=6\left({\beta }^2-3\beta +3\right)c_4{H}^{\prime }(1)-3\left(2\left(3{\beta }^2-9\beta +7\right)H^{\prime }(1)+3(\beta -2)\beta H^{\prime \prime }(1)\right)c_2{c}_3,\hfill \\ \hfill & F_3=2\left(6\left({\beta }^2-3\beta +2\right)H^{\prime }(1)+\beta \left((6\beta -9)H^{\prime \prime }(1)+\beta H^{\prime \prime \prime }(1)\right)\right)c_2^3.\hfill \end{array}$$

Subtracting $\overline{x}$ on both sides of the second step of (1),

$$\begin{array}{cc}\hfill e_{k+1}& =(x_k-\overline{x})-H\left(\frac{f^{\prime }(x_k)}{f^{\prime }(y_k)}\right)\frac{f\left(x_k\right)}{f^{\prime }\left(x_k\right)}\hfill \\ \hfill & =(1-H(1))e_k+\left(H(1)-2\beta H^{\prime }(1)\right)c_2{e}_k^2+(2H(1)+3(\beta -2)\beta H^{\prime }(1))c_3\hfill \\ \hfill & -2\left(H(1)+\beta (2(\beta -2)H^{\prime }(1)+\beta H^{\prime \prime }(1))c_2^2\right)e_k^3+O\left(e_k^4\right).\hfill \end{array}$$

(5)

Forcing the coefficients of $e_k$, $e_k^2$ and $e_k^3$ to be zero, we get $H(1)=1$, $H^{\prime }(1)=\frac{3}{4}$, $H^{\prime \prime }(1)=\frac{3}{4}$ and $\beta =\frac{2}{3}$. By replacing them in (5), we have

$$\begin{array}{c}\hfill e_{k+1}=\frac{1}{81}\left((117-32H^{\prime \prime \prime }(1))c_2^3-81c_2{c}_3+9c_4\right)e_k^4+O\left(e_k^5\right),\end{array}$$

and the optimal fourth-order of convergence is proven. □

2.1. Parametric Family

The proposed parametric family generalizes some existing classes of iterative methods. Let us recall the one-parameter families of fourth-order optimal iterative methods developed by Hueso et al. in [20],

$$\begin{array}{cc}\hfill y_k& =x_k-\frac{2}{3}\frac{f\left(x_k\right)}{f\left(x_k\right)},k=0,1,\dots ,\hfill \\ \hfill x_{k+1}& =x_k-\left[\frac{5-8\alpha }{8}+\alpha \left(\frac{f^{\prime }\left(x_k\right)}{f^{\prime }\left(y_k\right)}\right)+\frac{\alpha }{3}{\left(\frac{f^{\prime }\left(x_k\right)}{f^{\prime }\left(y_k\right)}\right)}^{-1}+\frac{9-8\alpha }{24}{\left(\frac{f^{\prime }\left(x_k\right)}{f^{\prime }\left(y_k\right)}\right)}^2\right]\frac{f\left(x_k\right)}{f^{\prime }\left(x_k\right)},\hfill \end{array}$$

(6)

and the parametric family defined by Cordero et al. in [12],

$$\begin{array}{cc}\hfill y_k& =x_k-\frac{2}{3}\frac{f\left(x_k\right)}{f^{\prime }\left(x_k\right)},k=0,1,\dots ,\hfill \\ \hfill x_{k+1}& =x_k-\left[1-\frac{3}{4}\left(\frac{f^{\prime }\left(y_k\right)}{f^{\prime }\left(x_k\right)}-1\right)+\frac{9}{8}{\left(\frac{f^{\prime }\left(y_k\right)}{f^{\prime }\left(x_k\right)}-1\right)}^2+\frac{\alpha }{6}{\left(\frac{f^{\prime }\left(y_k\right)}{f^{\prime }\left(x_k\right)}-1\right)}^3\right]\frac{f\left(x_k\right)}{f^{\prime }\left(x_k\right)}.\hfill \end{array}$$

(7)

Both families have the structure of the iterative function defined in equation (1) for $\beta =\frac{2}{3}$, and their weight functions satisfy Theorem 5. It is evident that the weight functions of these families share similar structures, suggesting the possibility of constructing a generalized family that encompasses both.

The following result shows that under certain conditions, the parametric weight function

$$\begin{array}{c}\hfill K(\mu (x))=a_0{\mu }^{k_0}(x)+a_1{\mu }^{k_1}(x)+a_2{\mu }^{k_2}(x)+a_3{\mu }^{k_3}(x),\end{array}$$

(8)

where $a_i\in \mathbb{R}$ and $k_i\in \mathbb{N}$ for $i=1,2,3$, satisfies the conditions of Theorem 5.

Theorem 6.

The parametric weight function $K(\mu (x))=a_0{\mu }^{k_0}(x)+a_1{\mu }^{k_1}(x)+a_2{\mu }^{k_2}(x)+a_3{\mu }^{k_3}(x)$ satisfies the conditions imposed to the weight function H of Theorem 5 if and only if

$$\begin{array}{cc}\hfill & k_i\ne k_j,\mathit{with}i\ne j,\hfill \\ & a_0=\frac{\left(k_1-k_2\right)\left(-3k_2+k_1\left(4k_2-3\right)+6\right)-\gamma {\varphi }_0}{4{\varphi }_3},\hfill \\ & a_1=\frac{\left(3\left(k_2-2\right)+k_0\left(3-4k_2\right)\right)\left(k_0-k_2\right)+\gamma {\varphi }_1}{4{\varphi }_3},\hfill \\ & a_2=\frac{\left(k_0-k_1\right)\left(-3k_1+k_0\left(4k_1-3\right)+6\right)-\gamma {\varphi }_2}{4{\varphi }_3},\hfill \\ & a_3=\gamma ,\hfill \end{array}$$

where

$$\begin{array}{cc}\hfill & {\varphi }_n=\prod _{\begin{array}{c}i<j\\ i,j\ne n\end{array}}(k_i-k_j),\mathit{with}i,j\mathit{and}n\mathit{in}\{0,1,2,3\}.\hfill \end{array}$$

Proof. 

Considering that

$$\begin{array}{ccc}\hfill K\left(\mu (x)\right)& =& a_0{\mu }^{k_0}(x)+a_1{\mu }^{k_1}(x)+a_2{\mu }^{k_2}(x)+a_3{\mu }^{k_3}(x),\hfill \\ \hfill K^{\prime }\left(\mu (x)\right)& =& a_0{k}_0{\mu }^{k_0-1}(x)+a_1{k}_1{\mu }^{k_1-1}(x)+a_2{k}_2{\mu }^{k_2-1}(x)+a_3{k}_3{\mu }^{k_3-1}(x),\hfill \\ \hfill K^{\prime \prime }\left(\mu (x)\right)& =& a_0\left(k_0-1\right)k_0{\mu }^{k_0-2}(x)+a_1\left(k_1-1\right)k_1{\mu }^{k_1-2}(x)+a_2\left(k_2-1\right)k_2{\mu }^{k_2-2}(x)\hfill \\ & & +a_3\left(k_3-1\right)k_3{\mu }^{k_3-2}(x),\hfill \end{array}$$

Then,

$$\begin{array}{cc}\hfill & K\left(1\right)=a_0+a_1+a_2+a_3,\hfill \\ & K^{\prime }\left(1\right)=a_0{k}_0+a_1{k}_1+a_2{k}_2+a_3{k}_3,\hfill \\ & K^{\prime \prime }\left(1\right)=a_0\left(k_0-1\right)k_0+a_1\left(k_1-1\right)k_1+a_2\left(k_2-1\right)k_2+a_3\left(k_3-1\right)k_3.\hfill \end{array}$$

If we equate These expressions to the conditions imposed in Theorem 5, then we have the following system

$$\begin{array}{ccc}\hfill a_0+a_1+a_2+a_3& =& 1,\hfill \\ \\ \hfill a_0{k}_0+a_1{k}_1+a_2{k}_2+a_3{k}_3& =& \frac{3}{4},\hfill \\ \\ \hfill a_0\left(k_0-1\right)k_0+a_1\left(k_1-1\right)k_1+a_2\left(k_2-1\right)k_2+a_3\left(k_3-1\right)k_3& =& \frac{3}{4},\hfill \end{array}$$

whose solution set for $a_3=\gamma$, and $k_i\ne k_j\mathrm{if}i\ne j$ is

$$\begin{array}{ccc}\hfill a_0& =& -\frac{4\gamma k_3^2-4\gamma k_1{k}_3-4\gamma k_2{k}_3+4\gamma k_1{k}_2+3k_1-4k_1{k}_2+3k_2-6}{4\left(k_0-k_1\right)\left(k_0-k_2\right)}=-\frac{3(k_2-2)+k_1\left(3-4k_2\right)}{4\left(k_0-k_1\right)\left(k_0-k_2\right)}+\frac{\gamma \left(k_1-k_3\right)\left(k_3-k_2\right)}{\left(k_0-k_1\right)\left(k_0-k_2\right)},\hfill \\ \\ \hfill a_1& =& -\frac{4\gamma k_3^2-4\gamma k_0{k}_3-4\gamma k_2{k}_3+4\gamma k_0{k}_2+3k_0-4k_0{k}_2+3k_2-6}{4\left(k_1-k_0\right)\left(k_1-k_2\right)}=\frac{3\left(k_2-2\right)+k_0\left(3-4k_2\right)}{4\left(k_0-k_1\right)\left(k_1-k_2\right)}+\frac{\gamma \left(k_0-k_3\right)\left(k_2-k_3\right)}{\left(k_0-k_1\right)\left(k_1-k_2\right)},\hfill \\ \\ \hfill a_2& =& -\frac{-4\gamma k_3^2+4\gamma k_0{k}_3+4\gamma k_1{k}_3-4\gamma k_0{k}_1-3k_0+4k_0{k}_1-3k_1+6}{4\left(k_1-k_2\right)\left(k_2-k_0\right)}=\frac{3\left(k_1-2\right)+k_0\left(3-4k_1\right)}{4\left(k_0-k_2\right)\left(k_2-k_1\right)}+\frac{\gamma \left(k_0-k_3\right)\left(k_1-k_3\right)}{\left(k_1-k_2\right)\left(k_2-k_0\right)},\hfill \\ \hfill a_3& =& \gamma .\hfill \end{array}$$

By defining

$$\begin{array}{cc}\hfill & {\varphi }_n=\prod _{\begin{array}{c}i<j\\ i,j\ne n\end{array}}(k_i-k_j),\mathrm{with}i,j\mathrm{and}n\mathrm{in}\{0,\cdots ,3\},\hfill \end{array}$$

and multiplying and dividing by ${\varphi }_3$, we have

$$\begin{array}{ccc}\hfill a_0& =& \frac{\left(k_1-k_2\right)\left(-3k_2+k_1\left(4k_2-3\right)+6\right)-\gamma \left(k_1-k_2\right)\left(k_1-k_3\right)\left(k_2-k_3\right)}{4{\varphi }_3}=\frac{\left(k_1-k_2\right)\left(-3k_2+k_1\left(4k_2-3\right)+6\right)-\gamma {\varphi }_0}{4{\varphi }_3},\hfill \\ \\ \hfill a_1& =& \frac{\left(3\left(k_2-2\right)+k_0\left(3-4k_2\right)\right)\left(k_0-k_2\right)+\gamma \left(k_0-k_2\right)\left(k_0-k_3\right)\left(k_2-k_3\right)}{4{\varphi }_3}=\frac{\left(3\left(k_2-2\right)+k_0\left(3-4k_2\right)\right)\left(k_0-k_2\right)+\gamma {\varphi }_1}{4{\varphi }_3},\hfill \\ \\ \hfill a_2& =& \frac{\left(k_0-k_1\right)\left(-3k_1+k_0\left(4k_1-3\right)+6\right)-\gamma \left(k_0-k_1\right)\left(k_0-k_3\right)\left(k_1-k_3\right)}{4{\varphi }_3}=\frac{\left(k_0-k_1\right)\left(-3k_1+k_0\left(4k_1-3\right)+6\right)-\gamma {\varphi }_2}{4{\varphi }_3},\hfill \\ \\ \hfill a_3& =& \gamma ,\hfill \end{array}$$

and the proof is finished. □

In Theorem 6, we observe that the parametric weight function $K(\mu (x))$ depends on the parameters $k_0,k_1,k_2,k_3,\mathrm{and}{a}_3$, so we denote it as $K_{(k_0,k_1,k_2,k_3,a_3)}(\mu (x))$. It is easy to see that the families of iterative methods from the iterative expressions (6) and (7) are subfamilies of the family formed by the iterative scheme (1) when $H(\mu )=K_{(-1,0,1,2,\frac{9-8\alpha }{24})}(\mu )$ and $H(\mu )=K_{(-3,-2,-1,0,\frac{1}{6}\alpha )}(\mu )$, respectively.

3. Complex Stability Analysis

In Section 2, we have set the conditions for the following iterative scheme to be a fourth-order optimal class of iterative methods.

$$\begin{array}{ccc}\hfill y_k& =& x_k-\frac{2}{3}\frac{f\left(x_k\right)}{f^{\prime }\left(x_k\right)},k=0,1,\dots ,\hfill \\ \hfill x_{k+1}& =& x_k-\left[a_0\mu {(x_k)}^{k_0}+a_1\mu {(x_k)}^{k_1}+a_2\mu {(x_k)}^{k_2}+a_3\mu {(x_k)}^{k_3}\right]\frac{f\left(x_k\right)}{f^{\prime }\left(x_k\right)},\hfill \end{array}$$

(9)

where $\mu (x_k)={\displaystyle \frac{f^{\prime }\left(x_k\right)}{f^{\prime }\left(y_k\right)}}$.

If we assign integer values to $k_0$, $k_1$, $k_2$, and $k_3$, all distinct from each other, and define $a_i$ for i in $\{0,1,2,3\}$ as in Theorem 6, then the following result demonstrates that the resulting one-parametric subfamily satisfies the scaling theorem. As a result, the Möbius transformation can be used as an analytical conjugation to analyze the stability of its members on quadratic polynomials by means of its associated rational operator.

Theorem 7.

Let $f(z)$ be an analytic function, let R be the rational operator associated with the iterative scheme (9), and let $A(z)={\eta }_{1}z+{\eta }_2$, with ${\eta }_1\ne 0$ an affine transformation. Let $h(z)=\lambda (f\circ A)(z)$, then the fixed point operators $R_f$ and $R_h$ are analytically conjugated by $A(z)$.

Proof. 

Let $N_f(z)=z-\frac{2}{3}\frac{f(z)}{f^{\prime }(z)}$ and considering that

$$\begin{array}{cc}\hfill A(x-y)& ={\eta }_1(x-y)+{\eta }_2=({\eta }_{1}x+{\eta }_2)-(\eta -1y+{\eta }_2)+{\eta }_2=A(x)-A(y)+{\eta }_2,\hfill \end{array}$$

and also $h(z)=\lambda (f\circ A)(z)=\lambda f(A(z)),\mathrm{then}{h}^{\prime }(z)=A^{\prime }(z)\lambda f^{\prime }(A(z))={\eta }_1\lambda (f^{\prime }\circ A)(z).$ Therefore,

$$\begin{array}{ccc}\hfill A\circ N_h\circ A^{-1}& =& A\left(A^{-1}(z)-\frac{2}{3}\frac{h\left(A^{-1}(z)\right)}{h^{\prime }\left(A^{-1}(z)\right)}\right)\hfill \\ & =& A\left(A^{-1}(z)-\frac{2}{3}\frac{f(z)}{{\eta }_1{f}^{\prime }(z)}\right)\hfill \\ & =& z-\frac{2}{3}\frac{f(z)}{f^{\prime }(z)}=N_f(z).\hfill \end{array}$$

So, $N_f$ and $N_h$ are analytically conjugated by $A(z)$. Then,

$$R_h=z-\left[\sum _{i=0}^3{a}_i{\left(\frac{h^{\prime }\left(z\right)}{h^{\prime }\left(N_h\right)}\right)}^{k_i}\right]\frac{h\left(z\right)}{h^{\prime }\left(z\right)}.$$

As

$$\begin{array}{cc}\hfill h(z)& =\lambda (f\circ A)(z)=\lambda f(A(z)),\hfill \\ \hfill h^{\prime }(z)& ={\eta }_1\lambda (f^{\prime }\circ A)(z),\hfill \end{array}$$

and

$$\begin{array}{cc}\hfill h^{\prime }(N_h)& ={\eta }_1\lambda (f^{\prime }\circ A\circ N_h)(z),\hfill \end{array}$$

then,

$$R_h=z-\left[\sum _{i=0}^3{a}_i{\left(\frac{f^{\prime }\circ A}{f^{\prime }\circ A\circ N_h}\right)}^{k_i}\right]\frac{1}{{\eta }_1}\frac{f\circ A}{f^{\prime }\circ A}.$$

Using that $A\circ N_h\circ A^{-1}=N_f$,

$$R_h\circ A^{-1}=A^{-1}(z)-\left[\sum _{i=0}^3{a}_i{\left(\frac{f^{\prime }(z)}{f^{\prime }\circ N_f}\right)}^{k_i}\right]\frac{1}{{\eta }_1}\frac{f(z)}{f^{\prime }(z)},$$

and therefore,

$$A\circ R_h\circ A^{-1}=A\left[A^{-1}(z)-\left[\sum _{i=0}^3{a}_i{\left(\frac{f^{\prime }(z)}{f^{\prime }\circ N_f}\right)}^{k_i}\right]\frac{1}{{\eta }_1}\frac{f(z)}{f^{\prime }(z)}\right],$$

as, $A(x-y)=A(x)-A(y)+{\eta }_2$, then,

$$A\circ R_h\circ A^{-1}=z-\left[\sum _{i=0}^3{a}_i{\left(\frac{f^{\prime }(z)}{f^{\prime }\circ N_f}\right)}^{k_i}\right]\frac{f(z)}{f^{\prime }(z)}=R_f,$$

so $R_f$ and $R_h$ are analytically conjugate by $A(z)$. □

3.1. Asymptotical Analysis of the Fixed Points of the Rational Operator

As a consequence of Theorem 7, we know that under the conditions established for the parameters $k_i$ and $a_i$ for $i\in \{0,1,2,3\}$ in class (9), the dynamical behavior of the associated rational operator remains invariant under the Möbius transformation through conjugation over quadratic polynomials.

Applying the family of iterative methods (9) on $p(z)=(z-a)(z-b)$, where $a,b\in \widehat{\mathbb{C}}$ we obtain the associated rational operator $R_K=\varphi \circ R_p\circ {\varphi }^{-1}$, after the conjugation with the Möbius map $\varphi (z)={\displaystyle \frac{z-a}{z-b}}$.

Now, for our dynamical analysis, we assign distinct integer values to $k_0,k_1,k_2$, and $k_3$, and we study the dynamics of the rational operator $R_K$ of the parametric family associated with these values. For the selection of the class, we use the following criteria:

1.

Since we aim to construct a parameter plane for each independent free critical point, we want the polynomial whose roots are critical points depending on $\gamma$ to have a degree no higher than 4. This is because it becomes extremely laborious to construct parameter planes when the critical points that depend on $\gamma$ are the roots of a high-degree polynomial.

2.

We want to find a radius disk where $|R_K^{\prime }(1,k_0,k_1,k_2,k_3,\gamma )|>1$ as big as possible. By making $z=1$ in $R_K^{\prime }(z,k_0,k_1,k_2,k_3,\gamma )$ it is easy to see that $R_K^{\prime }(1,k_0,k_1,k_2,k_3,\gamma )=\frac{A(k_0,k_1,k_2,k_3)}{B(k_0,k_1,k_2,k_3)+C(k_0,k_1,k_2,k_3)\gamma }$, where $A(k_0,k_1,k_2,k_3)=8(k_0-k_1)(k_0-k_2)(k_1-k_2),$

$$\begin{array}{cc}\hfill B(k_0,k_1,k_2,k_3)=& -{k_0}^2{k}_1{2}^{2-k_2}{3}^{k_2}+{k_0}^2{2}^{2-k_1}{3}^{k_1}{k}_2-{k_0}^2{2}^{-k_1}{3}^{k_1+1}+8{k_0}^2{k}_1\hfill \\ & +{k_0}^2{2}^{-k_2}{3}^{k_2+1}-8{k_0}^2{k}_2+k_0{k_1}^2{2}^{2-k_2}{3}^{k_2}-2^{2-k_0}{3}^{k_0}{k_1}^2{k}_2-8k_0{k_1}^2\hfill \\ & +2^{-k_0}{3}^{k_0+1}{k_1}^2-k_0{2}^{2-k_1}{3}^{k_1}{k_2}^2+2^{2-k_0}{3}^{k_0}{k}_1{k_2}^2\hfill \\ & +k_0{2}^{1-k_1}{3}^{k_1+1}{2}^{1-k_0}{3}^{k_0+1}{k}_1+8k_0{k_2}^2-2^{-k_0}{3}^{k_0+1}{k_2}^2-k_0{2}^{1-k_2}{3}^{k_2+1}\hfill \\ & +2^{1-k_0}{3}^{k_0+1}{k}_2-{k_1}^2{2}^{-k_2}{3}^{k_2+1}+8{k_1}^2{k}_2+2^{-k_1}{3}^{k_1+1}{k_2}^2-8k_1{k_2}^2\hfill \\ & +k_1{2}^{1-k_2}{3}^{k_2+1}-2^{1-k_1}{3}^{k_1+1}{k}_2,\hfill \end{array}$$

and

$$\begin{array}{cc}\hfill C(k_0,k_1,k_2,k_3)=& {k_0}^2{k}_1{2}^{2-k_2}{3}^{k_2}-{k_0}^2{2}^{2-k_1}{3}^{k_1}{k}_2-{k_0}^2{k}_1{2}^{2-k_3}{3}^{k_3}+{k_0}^2{2}^{2-k_1}{3}^{k_1}{k}_3+\hfill \\ \hfill & {k_0}^2{k}_2{2}^{2-k_3}{3}^{k_3}-{k_0}^2{2}^{2-k_2}{3}^{k_2}{k}_3-k_0{k_1}^2{2}^{2-k_2}{3}^{k_2}+2^{2-k_0}{3}^{k_0}{k_1}^2{k}_2+\hfill \\ & k_0{k_1}^2{2}^{2-k_3}{3}^{k_3}-2^{2-k_0}{3}^{k_0}{k_1}^2{k}_3+k_0{2}^{2-k_1}{3}^{k_1}{k_2}^2-2^{2-k_0}{3}^{k_0}{k}_1{k_2}^2-\hfill \\ & k_0{2}^{2-k_1}{3}^{k_1}{k_3}^2+2^{2-k_0}{3}^{k_0}{k}_1{k_3}^2-k_0{k_2}^2{2}^{2-k_3}{3}^{k_3}+2^{2-k_0}{3}^{k_0}{k_2}^2{k}_3+\hfill \\ & k_0{2}^{2-k_2}{3}^{k_2}{k_3}^2-2^{2-k_0}{3}^{k_0}{k}_2{k_3}^2-{k_1}^2{k}_2{2}^{2-k_3}{3}^{k_3}+{k_1}^2{2}^{2-k_2}{3}^{k_2}{k}_3+\hfill \\ & k_1{k_2}^2{2}^{2-k_3}{3}^{k_3}-2^{2-k_1}{3}^{k_1}{k_2}^2{k}_3-k_1{2}^{2-k_2}{3}^{k_2}{k_3}^2+2^{2-k_1}{3}^{k_1}{k}_2{k_3}^2.\hfill \end{array}$$

Then, $|\frac{A(k_0,k_1,k_2,k_3)}{C(k_0,k_1,k_2,k_3)}|>|\frac{B(k_0,k_1,k_2,k_3)}{C(k_0,k_1,k_2,k_3)}+\gamma |$ when $|R^{\prime }(1,k_0,k_1,k_2,k_3,\gamma )|>1$ and now we define $G(k_0,k_1,k_2,k_3)=\frac{A(k_0,k_1,k_2,k_3)}{C(k_0,k_1,k_2,k_3)}$ such that we can choose values of $k_0,k_1,k_2$ and $k_3$ where $|G(k_0,k_1,k_2,k_3)|$ can be made as large as desired.

3.

The coefficients of the rational operator must be simple in the sense that they do not exceed three digits.

In Table 1, the degrees of the polynomials whose zeros are free critical points of the rational operator associated with 330 one-parameter families constructed with values of $k_0$, $k_1$, $k_2$, and $k_3$ taken from the set of all unique permutations of integers from -5 to 5 are presented. Within this set, 9 families are identified that meet the criteria established in the first point. Of these, 8 families have an associated polynomial of degree 4, while one family has a polynomial of degree 2.

Table 2 presents the values of $k_j$ for $j=1,2,3,4$, radius $|G(k_0,k_1,k_2,k_3)|$ and rational function of those classes of iterative schemes such that the polynomials generating the free critical points have degree lower or equal to four. It is observed that the families by Hueso et al. (6) and Cordero et al. (7) meet our selection criteria. In the case of the first family, the radius of the disk where divergence is repelling, as a function of $\gamma$, is 24, and the polynomial that meets the first criterion is of degree 4. On the other hand, the second class has a radius as a function of $\gamma$ of 54, and it possesses the only second-degree polynomial that satisfies the first criterion among all possible formed families. In both families, the coefficients of the rational operators do not exceed three digits.

The family associated with $K_{(-4,-3,-2,-1,\gamma )}$ has a divergence radius greater than that of the two families mentioned earlier, and the polynomial degree is 4. However, the coefficients of the rational operator are not simple. On the other hand, the family associated with $K_{(-2,-1,0,1,\gamma )}$ has simple coefficients and a radius of 36, which is greater than the radius of the family studied by Hueso et al. (6). Therefore, we will focus our analysis on this family. It is worth mentioning that it would also be interesting to study the dynamical behavior of the family associated with $K_{(0,1,2,3,\gamma )}$ or $K_{(1,2,3,4,\gamma )}$, as well as other families if we do not consider the simple coefficient criterion.

From Table 2, we denote the rational operator of the family associated with $K_{(-2,-1,0,1,\gamma )}$ by

$$R(z,\gamma )=\frac{z^4\left(27z^4+126z^3+234z^2+198z+135-64\gamma \right)}{\left(135-64\gamma \right)z^4+198z^3+234z^2+126z+27}.$$

Let us now analyze some properties of the rational operator R.

Theorem 8.

Fixed points of the rational function $R(z,\gamma )$ are $z=0$, $z=\infty$ and also the following strange fixed points:

• $e{x}_1(\gamma )=1$ which corresponds to divergence from the original method,
• $e{x}_2(\gamma )=z_{1,1}=\frac{1}{2}\left(\frac{1}{9}S(\gamma )-\sqrt{{\left(-\frac{1}{9}S(\gamma )+\frac{17}{9S(\gamma )}+\frac{17}{9}\right)}^2-4}-\frac{17}{9S(\gamma )}-\frac{17}{9}\right),$
• $e{x}_3(\gamma )=z_{1,2}=\frac{1}{2}\left(\frac{1}{9}S(\gamma )+\sqrt{{\left(-\frac{1}{9}S(\gamma )+\frac{17}{9S(\gamma )}+\frac{17}{9}\right)}^2-4}-\frac{17}{9S(\gamma )}-\frac{17}{9}\right),$
• $e{x}_4(\gamma )=z_{2,1}=\frac{1}{2}\left(-\frac{1}{9}wS(\gamma )-\sqrt{-4+{\left(\frac{1}{9}wS(\gamma )-\frac{17\overline{w}}{9S(\gamma )}+\frac{17}{9}\right)}^2}+\frac{17\overline{w}}{9S(\gamma )}-\frac{17}{9}\right),$
• $e{x}_5(\gamma )=z_{2,2}=\frac{1}{2}\left(-\frac{1}{9}wS(\gamma )+\sqrt{-4+{\left(-\frac{1}{9}wS(\gamma )-\frac{17\overline{w}}{9S(\gamma )}+\frac{17}{9}\right)}^2}+\frac{17\overline{w}}{9S(\gamma )}-\frac{17}{9}\right),$
• $e{x}_6(\gamma )=z_{3,1}=\frac{1}{2}\left(-\frac{1}{9}\overline{w}S(\gamma )-\sqrt{-4+{\left(\frac{1}{9}\overline{w}S(\gamma )-\frac{17w}{9S(\gamma )}+\frac{17}{9}\right)}^2}+\frac{17w}{9S(\gamma )}-\frac{17}{9}\right),$
• $e{x}_7(\gamma )=z_{3,2}=\frac{1}{2}\left(-\frac{1}{9}\overline{w}S(\gamma )+\sqrt{-4+{\left(\frac{1}{9}\overline{w}S(\gamma )-\frac{17w}{9S(\gamma )}+\frac{17}{9}\right)}^2}+\frac{17w}{9S(\gamma )}-\frac{17}{9}\right),$

where $S(\gamma )=\sqrt[3]{3\sqrt{3}\sqrt{27648{\gamma }^2-60544\gamma +33327}-864\gamma +946}$ and $w=\frac{1}{2}\left(1-i\sqrt{3}\right).$

Proof. 

By fixed point definition, we have $R(z,\gamma )=z$. Then,

$$\frac{z^4\left(27z^4+126z^3+234z^2+198z+135-64\gamma \right)}{\left(135-64\gamma \right)z^4+198z^3+234z^2+126z+27}=z.$$

So,

$$z\left[27\left(z^7-1\right)+126z\left(z^5-1\right)+234z^2\left(z^3-1\right)+\left(63+64\gamma \right)z^3\left(z-1\right)\right]=0,$$

and

$$z\left(z-1\right)\left(27z^6+153z^5+387z^4+z^3(450+64\gamma )+387z^2+153z+27\right)=0.$$

Let $J(z,\gamma )$ denote the last factor of this equation. If we substitute $u=z+\frac{1}{z}$, then $u^2-2=z^2+\frac{1}{z^2}$ and $u^3-3u=z^3+\frac{1}{z^3}$. Thus, $J(z,\gamma )=z^3\left[27u^3+153u^2+306u+(144+64\gamma )\right]$. Therefore, the factor $27u^3+153u^2+306u+(144+64\gamma )$ must be zero. Consequently, we have

$$u_1(\gamma )=\frac{1}{9}S(\gamma )-\frac{17}{9S(\gamma )}-\frac{17}{9},$$

$$u_2(\gamma )=-\frac{1}{18}\left(1-i\sqrt{3}\right)S(\gamma )+\frac{17\left(1+i\sqrt{3}\right)}{18S(\gamma )}-\frac{17}{9},$$

$$u_3(\gamma )=-\frac{1}{18}\left(1+i\sqrt{3}\right)S(\gamma )+\frac{17\left(1-i\sqrt{3}\right)}{18S(\gamma )}-\frac{17}{9},$$

where

$$S(\gamma )=\sqrt[3]{3\sqrt{3}\sqrt{27648{\gamma }^2-60544\gamma +33327}-864\gamma +946}.$$

Therefore, $e{x}_{j+k}(\gamma )=\frac{u_i(\gamma )+{(-1)}^{j+k}\sqrt{u_i{(\gamma )}^2-4}}{2},$ for $i=1,2,3$, $j=1,3,5$ and $k=1,2.$ Let us notice that the strange fixed points $e{x}_{i+1}(\gamma )=z_{i,2}=\frac{1}{z_{i,1}}=\frac{1}{e{x}_i(\gamma )}$, since $z^2-u_i(\gamma )z+1=0$ is a symmetric polynomial.

• $e{x}_2(\gamma )=z_{1,1}=\frac{1}{2}\left(\frac{1}{9}S(\gamma )-\sqrt{{\left(-\frac{1}{9}S(\gamma )+\frac{17}{9S(\gamma )}+\frac{17}{9}\right)}^2-4}-\frac{17}{9S(\gamma )}-\frac{17}{9}\right),$
• $e{x}_3(\gamma )=z_{1,2}=\frac{1}{2}\left(\frac{1}{9}S(\gamma )+\sqrt{{\left(-\frac{1}{9}S(\gamma )+\frac{17}{9S(\gamma )}+\frac{17}{9}\right)}^2-4}-\frac{17}{9S(\gamma )}-\frac{17}{9}\right),$
• $e{x}_4(\gamma )=z_{2,1}=\frac{1}{2}\left(-\frac{1}{9}wS(\gamma )-\sqrt{-4+{\left(\frac{1}{9}wS(\gamma )-\frac{17\overline{w}}{9S(\gamma )}+\frac{17}{9}\right)}^2}+\frac{17\overline{w}}{9S(\gamma )}-\frac{17}{9}\right),$
• $e{x}_5(\gamma )=z_{2,2}=\frac{1}{2}\left(-\frac{1}{9}wS(\gamma )+\sqrt{-4+{\left(-\frac{1}{9}wS(\gamma )-\frac{17\overline{w}}{9S(\gamma )}+\frac{17}{9}\right)}^2}+\frac{17\overline{w}}{9S(\gamma )}-\frac{17}{9}\right),$
• $e{x}_6(\gamma )=z_{3,1}=\frac{1}{2}\left(-\frac{1}{9}\overline{w}S(\gamma )-\sqrt{-4+{\left(\frac{1}{9}\overline{w}S(\gamma )-\frac{17w}{9S(\gamma )}+\frac{17}{9}\right)}^2}+\frac{17w}{9S(\gamma )}-\frac{17}{9}\right),$
• $e{x}_7(\gamma )=z_{3,2}=\frac{1}{2}\left(-\frac{1}{9}\overline{w}S(\gamma )+\sqrt{-4+{\left(\frac{1}{9}\overline{w}S(\gamma )-\frac{17w}{9S(\gamma )}+\frac{17}{9}\right)}^2}+\frac{17w}{9S(\gamma )}-\frac{17}{9}\right),$

where $w=\frac{1}{2}\left(1-i\sqrt{3}\right).$

In a similar way, it can be easily checked that $z=0$ is a fixed point of $\frac{1}{R\left(\frac{1}{z},\gamma \right)}$, and thus, $z=\infty$ is a fixed point of R. □

It is clear that the stability of fixed points is influenced by the parameter of the family. This dependence can result in the presence of attracting strange fixed points, leading the corresponding iterative method to attracting elements that are not solutions to the problem being solved.

In order to perform the stability analysis of the fixed points, we evaluate the derivative operator in them,

$$R^{\prime }(z,\gamma )=-\frac{36z^3{(z+1)}^{4}N(z,\gamma )}{{\left(64\gamma z^4-135z^4-198z^3-234z^2-126z-27\right)}^2},$$

where $N(z,\gamma )=(192\gamma -405)z^4-(96\gamma +540)z^3+(64\gamma -990)z^2-(96\gamma +540)z+(192\gamma -405)$.

Fixed points $z=0$ and $z=\infty$ are always superattracting since $|R^{\prime }(z,\gamma )|=0$ for those values. However, the stability of the other fixed points, such as $z=1$, provides important numerical information, as it corresponds to the divergence of the original scheme. Hence, we will determine the stability of these fixed points in the forthcoming results.

Theorem 9.

The strange fixed point $e{x}_1(\gamma )=1$, if $\gamma \ne \frac{45}{4}$, has the following characteristics:

(i) 

It can not be superattracting.

(ii) 

When $\left|\frac{45}{4}-\gamma \right|<36$, $e{x}_1(\gamma )=1$ is a repelling fixed point.

(iii) 

If $\left|\frac{45}{4}-\gamma \right|=36$, $e{x}_1(\gamma )=1$ is a neutral or indifferent fixed point.

(iv) 

When $\left|\frac{45}{4}-\gamma \right|>36$, then $e{x}_1(\gamma )=1$ is an attracting fixed point.

Proof. 

(i)

The derivative $R_h^{\prime }(1,\gamma )=\frac{144}{45-4\gamma }$ is always different from zero, so $e{x}_1(\gamma )=1$ cannot be a superattracting. Moreover, it is not defined for $\gamma =\frac{45}{4}$, and $z=1$ is not a fixed point for this value of $\gamma$.

1.

(ii) According to the definition, $e{x}_1(\gamma )=1$ is a repelling fixed point if $\left|R_h^{\prime }(1,\gamma ))\right|=\left|\frac{144}{45-4\gamma }\right|>1$. Therefore, $\left|\frac{144}{45-4\gamma }\right|>1$ is equivalent to $\left|\frac{45}{4}-\gamma \right|<\frac{144}{4}=36$.

Let us consider an arbitrary complex number $\gamma =c+id$. Then, $\left|\frac{45}{4}-\gamma \right|<36$ is equivalent to ${(c-\frac{45}{4})}^2+d^2<{36}^2$, indicating that the divergence represents a repelling fixed point for values of $\gamma \ne \frac{45}{4}$ inside the disk $D_1$ centered at $C\left(\frac{45}{4},0\right)$ with a radius of 36. It is neutral at the circumference and attracting outside the disk. See Figure 1 for a visual representation.

The proof for items (iii) and (iv) follow a similar approach to item (ii).

□

We aim for the only attracting fixed points to be $z=0$ and $z=\infty$. Therefore, the stability function $\left|R^{\prime }(z,\gamma )\right|$ at the other fixed points should always be greater than one. Indeed, to prevent divergence, an stable method must be defined with a $\gamma$ value from within the disk $D_1$.

Figure 2 displays the stability function of $e{x}_2(\gamma )$ and $e{x}_3(\gamma )$, indicating that there are only attracting fixed points inside the orange cardioid contained within the rectangular region defined by

$$B(\gamma )=\{\gamma \in \mathbb{C}:\gamma =c+id,\frac{22}{5}<c<\frac{99}{8},-\frac{83}{20}<d<\frac{83}{20}\}.$$

Figure 3 displays the stability surfaces of the strange fixed points $e{x}_i(\gamma )$ for $i=4,5,6,7$, illustrating that they are always repelling for $\gamma$ values inside the disk $D_1$.

On the other hand, we know that $z=0$ and $z=\infty$ are critical points for all values of $\gamma$. The following result identifies the critical points of the rational operator, which allow us to determine for which values of $\gamma$ the zeros of $p(z)$ are the only attracting elements.

Theorem 10.

The rational operator R has the following free critical points:

i) 

$c{r}_1=-1,$

ii) 

$c{r}_2(\gamma )=z_{1,1}=\frac{1}{2}\left(v_1(\gamma )-\sqrt{v_1^2(\gamma )-4}\right),$

iii) 

$c{r}_3(\gamma )=z_{1,2}=\frac{1}{2}\left(v_1(\gamma )+\sqrt{v_1^2(\gamma )-4}\right),$

iv) 

$c{r}_4(\gamma )=z_{2,1}=\frac{1}{2}\left(v_2(\gamma )-\sqrt{v_2^2(\gamma )-4}\right),$

v) 

$c{r}_5(\gamma )=z_{2,2}=\frac{1}{2}\left(v_2(\gamma )+\sqrt{v_2^2(\gamma )-4}\right),$

where $v_1(\gamma )=\frac{2\left(-8\sqrt{3}\sqrt{83{\gamma }^2-90\gamma }+24\gamma +135\right)}{3(64\gamma -135)},$ and $v_2(\gamma )=\frac{2\left(8\sqrt{3}\sqrt{83{\gamma }^2-90\gamma }+24\gamma +135\right)}{3(64\gamma -135)}.$

Proof. 

It is straightforward that

$$R^{\prime }(z,\gamma )=-\frac{36z^3{(z+1)}^{4}M(z,\gamma )}{{\left(64\gamma z^4-135z^4-198z^3-234z^2-126z-27\right)}^2},$$

where $M(z,\gamma )=(192\gamma -405)z^4-(96\gamma +540)z^3+(64\gamma -990)z^2-(96\gamma +540)z+(192\gamma -405)$. Hence, the free critical points are $c{r}_1=-1$ and the solutions of the equation $M(z,\gamma )=0$. So, it can be expressed as

$$z^2\left[(192\gamma -405)z^2-(96\gamma +540)z+(64\gamma -990)-(96\gamma +540)\frac{1}{z}+(192\gamma -405)\frac{1}{z^2}\right]=0,$$

and therefore,

$$z^2\left[(192\gamma -405)(z^2+\frac{1}{z^2})-(96\gamma +540)(z+\frac{1}{z})+(64\gamma -990)\right]=0.$$

By making $v=z+\frac{1}{z}$, then $v^2-2=z^2+\frac{1}{z^2}$. Then,

$$z^2\left[(192\gamma -405)(v^2-2)-(96\gamma +540)v+(64\gamma -990)\right]=0,$$

and

$$z^2\left[(192\gamma -405)v^2-(96\gamma +540)v-(320\gamma +180)\right]=0.$$

Since $z\ne 0$, then

$$(192\gamma -405)v^2-(96\gamma +540)v-(320\gamma +180)=0.$$

Therefore, we have

$$v_1=\frac{2\left(-8\sqrt{3}\sqrt{83{\gamma }^2-90\gamma }+24\gamma +135\right)}{3(64\gamma -135)},$$

$$v_2=\frac{2\left(8\sqrt{3}\sqrt{83{\gamma }^2-90\gamma }+24\gamma +135\right)}{3(64\gamma -135)}.$$

So, the free critical points are $c{r}_{j+k}(\gamma )=\frac{v_i(\gamma )+{(-1)}^{j+k}\sqrt{v_i{(\gamma )}^2-4}}{2},$ for $i=1,2$, $j=1,3$ and $k=1,2.$

Let us notice that $c{r}_{i+1}(\gamma )=\frac{1}{c{r}_i(\gamma )}$ for $i=1,2$. Therefore, $c{r}_i(\gamma )=\frac{v_i(\gamma )\pm \sqrt{v_i{(\gamma )}^2-4}}{2},$

$$c{r}_2(\gamma )=\frac{1}{2}\left(v_1(\gamma )-\sqrt{v_1^2(\gamma )-4}\right),$$

$$c{r}_3(\gamma )=\frac{1}{2}\left(v_1(\gamma )+\sqrt{v_1^2-4}(\gamma )\right),$$

$$c{r}_4(\gamma )=\frac{1}{2}\left(v_2(\gamma )-\sqrt{v_2^2-4}(\gamma )\right),$$

$$c{r}_5(\gamma )=\frac{1}{2}\left(v_2(\gamma )+\sqrt{v_2^2-4}(\gamma )\right).$$

□

3.2. Parameter Plane

By means of Theorem 3, we employ critical points to locate the basins of attraction of periodic or fixed attracting points. To get this aim, we use the critical points as initial guesses to locate, through parameter planes, the values of $\gamma$ where the only attracting elements are the roots of polynomial $p(z)$.

We plot the parameter plane for values of $\gamma$ in the rectangular region defined by

$$\{\gamma \in \mathbb{C}:\gamma =c+id,-\frac{99}{4}<c<\frac{198}{4},-36<d<36\},$$

since this region contains the disk $\left|\frac{45}{4}-\gamma \right|=36$. To generate the plot, we use a grid of $2000\times 2000$ points, a maximum of 200 iterations, and a tolerance of ${10}^{-3}$. We recall that for the construction of the parameter and dynamical planes, we use the codes defined in Chicharro et al. [8]. Given a free critical point as the initial estimate for our iterative scheme, we define $D_{c{r}_i}$ as the loci of the complex plane defined by the values of $\gamma$ where the free critical point is in the basin of attraction of $z=0$ or $z=\infty$. The parameter plane is represented by the graph of the characteristic function of set $D_{c{r}_i}$, where the color red is assigned when $F_{D_{c{r}_i}}(\gamma )=1$ and black when $F_{D_{c{r}_i}}(\gamma )=0$. We plot the parameter plane for $\gamma$ values in the interior of the disk $D_1$ (defined in Theorem 9), as this region guarantees that the divergence since $e{x}_1(\gamma )=1$ is repelling.

Figure 4 displays the parameter planes corresponding to the critical points $c{r}_i(\gamma )$, where $i=2,3,4,5$. Let us remark that conjugate critical points have the same parameter plane, so only two of them are independent.

On the other hand, Figure 5 represents the unified parameter plane corresponding to the characteristic function $F_{{\bigcap }_{i=2}^5{D}_{c{r}_i}}(\gamma )$, where each point colored in red corresponds to values of $\gamma$ for which all members of our family of iterative schemes are stable on quadratic polynomials.

3.3. Dynamical Planes of Some Stable and Unstable Members

In Figure 5, we represent in red the values of $\gamma$ associated with stable members. To obtain additional complementary information that allows us to verify the stability characterization obtained through the parameter planes and the stability functions of the strange fixed points, we present in Figure 6 the corresponding dynamical planes for some stable and unstable values, as per the information obtained from the intersection of parameter planes.

To construct the dynamical planes in Figure 6, we define a grid in the complex plane with a $2000\times 2000$ point grid, where each one corresponds to a different value of the initial estimate $z_0$. Each dynamical plane shows the final state of the orbit of each point in the grid, considering a maximum of 200 iterations and a tolerance of ${10}^{-3}$. In each dynamical plane, we mark some special points. Repelling fixed points (RFP) with circles (∘), critical points (CP) with squares (□), and attracting fixed points (AFP) with asterisks (∗).

Through the dynamical planes, we determine to which attracting elements the orbit of any initial estimate $z_0$ will converge when using the rational operator R with a specific value of $\gamma$. The basin of attraction of $z=0$ is represented in orange, while the basin of $z=\infty$ is represented in blue. Additionally, we assign different colors such as green, red, etc., to other attracting strange fixed points, and we use black to represent basins of periodic orbits.

We have confirmed that the values of $\gamma$ in the red region of Figure 5 have only two basins of attraction, namely, $z=0$ and $z=\infty$. This guarantees that, regardless of the initial estimate $z_0$, the iterative scheme associated with these values of $\gamma$ will always converge to a root of $p(z)$. Similarly, it is validated that there are always pathologies that depend on the initial estimate $z_0$ for the methods associated with values of $\gamma$ in the black region of the parameter plane.

In Table 3, we present a classification of stability of the methods associated with the values of $\gamma$ presented in Figure 6, along with the number of special points and the period of the attracting orbits found of each $\gamma$.

4. Numerical Tests

In this section, we conduct numerical tests on some academical functions. Our goal is to determine whether methods that are stable in the previous analysis for quadratic polynomials perform better than unstable methods on non-polynomial functions. Subsequently, we compare the most effective methods on these functions with other well-known iterative methods in the field of numerical analysis, such as:

Newton’s method, whose iterative expression is

$$\begin{array}{cc}\hfill x_{k+1}& =x_k-\frac{f\left(x_k\right)}{f^{\prime }\left(x_k\right)},k=0,1,\dots ,\hfill \end{array}$$

Chun’s method, defined in [11] as

$$\begin{array}{cc}\hfill y_k& =x_k-\frac{f\left(x_k\right)}{f^{\prime }\left(x_k\right)},\hfill \\ \hfill x_{k+1}& =y_k-\frac{f\left(x_k\right)+2f\left(y_k\right)}{f\left(x_k\right)}\frac{f\left(y_k\right)}{f^{\prime }\left(x_k\right)},k=0,1,\dots ,\hfill \end{array}$$

Ostrowski’s method, proposed in [26] and expressed as

$$\begin{array}{cc}\hfill y_k& =x_k-\frac{f\left(x_k\right)}{f^{\prime }\left(x_k\right)},\hfill \\ \hfill x_{k+1}& =y_k-\frac{f\left(x_k\right)}{f\left(x_k\right)-2f\left(y_k\right)}\frac{f\left(y_k\right)}{f^{\prime }\left(x_k\right)},k=0,1,\dots ,\hfill \end{array}$$

Kung-Traub’s method, defined in [24] with the iterative expression

$$\begin{array}{cc}\hfill y_k& =x_k-\frac{f\left(x_k\right)}{f^{\prime }\left(x_k\right)},\hfill \\ \hfill x_{k+1}& =y_k-\frac{f{\left(x_k\right)}^2}{{\left(f\left(x_k\right)-f\left(y_k\right)\right)}^2}\frac{f\left(y_k\right)}{f^{\prime }\left(x_k\right)},k=0,1,\dots ,\hfill \end{array}$$

and Jarratt’s method, constructed in [21] as

$$\begin{array}{cc}\hfill y_k& =x_k-\frac{2}{3}\frac{f\left(x_k\right)}{f^{\prime }\left(x_k\right)},\hfill \\ \hfill x_{k+1}& =x_k-\frac{1}{2}\left(\frac{3f^{\prime }\left(y_k\right)+f^{\prime }\left(x_k\right)}{3f^{\prime }\left(y_k\right)-f^{\prime }\left(x_k\right)}\right)\frac{f\left(x_k\right)}{f^{\prime }\left(x_k\right)},k=0,1,\dots .\hfill \end{array}$$

For our comparative analysis, we consider as the best performing methods those iterative schemes that converge in the least number of iterations, the estimation of the error is as small as possible and the theoretical order is closely approximated by the computational convergence order (ACOC), defined by Cordero and Torregrosa in [13],

$$ACOC=\frac{ln|x_{k+1}-x_k|/|x_k-x_{k-1}|}{ln|x_k-x_{k-1}|/|x_{k-1}-x_{k-2}|}.$$

Test functions and their zeros to be estimated are the following,

$$\begin{array}{cc}{f}_1(x)=\sqrt{x^2+2x+5}-2sin(x)-x^2+3,\hfill & \overline{x}\approx 2.33196765588396\hfill \\ f_2(x)={(x-1)}^3-1+e^x,\hfill & \overline{x}\approx 0.297570665698774\hfill \\ f_3(x)=arctan(x)+e^{x^3},\hfill & \overline{x}\approx -0.757044010769526\hfill \\ f_4(x)=x{e}^{x^2}-cos(x),\hfill & \overline{x}\approx 0.588401776500996\hfill \end{array}$$

We present numerical results obtained using a Thinkpad T480s laptop equipped with an eighth-generation Core i7 processor (Intel(R) Core(TM) i7-8650U CPU @ 1.90GHz 2.11 GHz) and 16 GB of RAM. We have used MATLAB R2021b’s academic version, employing variable-precision arithmetics with 200 digits.

Table 4 displays the number of iterations required for each proposed method, using the parameter values specified in Table 3, to achieve the desired convergence for each of the functions. It can be observed that methods associated with values of $\gamma \ne \frac{45}{4}$, which are stable for real quadratic polynomials, demonstrate good performance on these test functions. However, it has been found that methods associated with stable complex $\gamma$ values, as well as methods associated with unstable $\gamma$ values, fail to converge in at least one of the test functions. For $\gamma =\frac{45}{4}$, the associated method does not converge for either $f_2$ or $f_4$. It is worth noting that the color intensity in the dynamical plane associated with this value of $\gamma$ indicates the number of iterations required to ensure convergence, the latter is more intense than those associated with the other stable real values of $\gamma$, indicating that those converging in fewer iterations on quadratic polynomials are more reliable when applied to non-polynomial test functions (see Figure 6).

As shown in Table 4, seven of the methods associated with real values guarantee converge in a maximum of nine iterations. Table 5, Table 6, Table 7 and Table 8 present the results obtained when applying Newton, Chun, Ostrowski, Kung and Traub, Jarratt, and the seven proposed methods. These tables provide information about the function to which each iterative method is applied, the initial estimate used ($x_0$), the target root ($\overline{x}$) to which the methods are expected to converge, the number of iterations (Iter) required for convergence, the approximate order of computational convergence (ACOC), execution time (Time), estimation of the error in the last iteration ($\left|x_{k+1}-x_k\right|$ and $\left|f\left(x_{k+1}\right)\right|$), and the approximate solution $\overline{x}$ found.

The results suggest that the methods associated with $\gamma$ values demonstrate competitive performance compared to Kung and Traub, Newton, Chun, Ostrowski and Jarratt methods based on our measurement parameters.

5. Conclusions

In this study, a set of conditions that weight functions must satisfy to ensure that the resulting iterative scheme has a convergence order of at least 4 has been established. Furthermore, a generalization of known uniparametric families of fourth-order optimal iterative methods has been presented. Criteria for selecting the best subfamilies of these uniparametric families have been established. Additionally, a dynamical analysis has been conducted, allowing us to identify stable members within the subfamily that exhibit good performance on quadratic polynomials. However, elements with pathological behavior have also been observed, which should be avoided in numerical applications.

Furthermore, the performance of stable members of the subfamily has been evaluated in comparison to unstable members when applied to non-polynomial functions. These results highlight the importance of stability in quadratic polynomials and its positive influence on the overall performance of iterative methods in various contexts.

Finally, stable methods within the subfamily that compete favorably with widely used iterative schemes in the literature have been identified. Through academical examples and detailed comparisons, it has been shown that the proposed methods represent viable and effective alternatives for solving nonlinear equations.