A New Perspective on the Convergence of Mean-Based Methods for Nonlinear Equations

1. Introduction

In mathematics and engineering, many real-world and physical phenomena are modeled using nonlinear equations or systems. Solving such problems has led to the development of numerous numerical methods, which serve as fundamental tools for approximating solutions that cannot be found exactly or analytically.

Among these, iterative methods play a crucial role in finding the roots of nonlinear equations. Since a variety of iterative strategies exist, it becomes essential to evaluate their convergence order, stability, and computational efficiency. These aspects allow us to compare methods and choose the most appropriate one for a given problem.

Iterative methods are typically classified based on whether they are single-step or multi-step, with or without memory, and whether they require derivatives or not . In this context, we present new multipoint iterative methods aimed at approximating the zeros of nonlinear functions. These methods are inspired by the work in [1], which investigates and enhances Newton-type methods by incorporating convex combinations of classical means, achieving third-order convergence.

With the aim of improving this order of convergence, we introduce new multipoint iterative schemes based on the composition and in the use of weight function. These functions combine data from multiple evaluations of the function and its derivatives throughout previous iterations to improve both accuracy and efficiency. The proposed methods achieve fourth-order convergence and satisfy the optimality condition postulated by the Kung-Traub conjecture [2], which states:

The order of convergence p of a without memory iterative method cannot exceed $2^{d-1}$, where d is the number of functional evaluations per iteration.

A method that reaches this bound is called optimal. This and other criteria described below will allow us to define which method is more effective than another.

One of the most widely used and foundational approaches in nonlinear root-finding is Newton’s method. It is defined as:

$$x_{n+1}=x_n-{\displaystyle \frac{f\left(x_n\right)}{f^{\prime }\left(x_n\right)}},n=0,1,2,\dots ,$$

(1)

provided $f^{\prime }\left(x_n\right)\ne 0$. Under appropriate smoothness conditions and for simple roots, the method exhibits quadratic convergence, meaning:

$$|x_{n+1}-\alpha |\approx C|x_n{-\alpha |}^2,$$

for some constant $C>0$ and being $\alpha$ a root of $f\left(x\right)=0$.

In 2000, Weerakoon and Fernando [3] proposed a third-order variant of Newton’s method. This method replaces the rectangular approximation in the integral form of Newton’s method with a trapezoidal approximation, reducing truncation error and improving convergence. Their method, known as the trapezoidal or arithmetic mean method, is defined as:

$$y_n=x_n-{\displaystyle \frac{f\left(x_n\right)}{f^{\prime }\left(x_n\right)}},x_{n+1}=x_n-{\displaystyle \frac{2f\left(x_n\right)}{f^{\prime }\left(x_n\right)+f^{\prime }\left(y_n\right)}},n=0,1,2,\dots .$$

(2)

This method laid the foundation for subsequent generalizations using other types of means. Researchers such as Chicharro et al. [4] and Cordero et al. [1] expanded this idea by incorporating various mathematical means to construct families of third-order methods:

$$x_{n+1}=x_n-{\displaystyle \frac{f\left(x_n\right)}{M_m\left[f^{\prime }\left(x_n\right),f^{\prime }\left(y_n\right)\right]}},n=0,1,2,\dots ,$$

(3)

where $y_n$ denotes the Newton step and $M_m(x,y)$ represents the chosen mean applied to the values x and y.

1.1. Types of Means

Below are the different types of convex averages used in the literature for the design and analysis of various iterative methods. These concepts constitute a fundamental reference and will serve as a methodological basis for the development of our own iterative procedures.

Arithmetic Mean $M_A$: The arithmetic mean of two real numbers x and y is given by:

$$M_A(x,y)={\displaystyle \frac{x+y}{2}}.$$

This mean appears in the trapezoidal scheme (2).

Harmonic Mean $M_H$: The harmonic mean of two positive real numbers x and y is defined as:

$$M_H(x,y)={\displaystyle \frac{2xy}{x+y}}.$$

This mean is particularly sensitive to small values and is known for its use in rates and resistances. In the context of iterative methods, its reciprocal nature often yields improved stability under specific conditions. The following scheme arises by replacing in (2) the arithmetic mean by this mean, as done in [5]:

$$x_{n+1}=x_n-{\displaystyle \frac{f\left(x_n\right)(f^{\prime }\left(x_n\right)+f^{\prime }\left(y_n\right))}{2f^{\prime }\left(x_n\right)f^{\prime }\left(y_n\right)}},n=0,1,2,\dots$$

Contraharmonic Mean $M_C$: The contraharmonic mean is given by:

$$M_C(x,y)={\displaystyle \frac{x^2+y^2}{x+y}},$$

and is always greater than or equal to the arithmetic mean. It accentuates larger values, making it suitable when higher magnitudes dominate the behavior of the function. When this mean is incorporated into iterative schemes, the obtained method presented in [1,6] is:

$$x_{n+1}=x_n-{\displaystyle \frac{(f^{\prime }\left(x_n\right)+f^{\prime }\left(y_n\right))f\left(x_n\right)}{{\left(f^{\prime }\left(x_n\right)\right)}^2+{\left(f^{\prime }\left(y_n\right)\right)}^2}},n=0,1,2,\dots$$

Different authors have proved that all these schemes have order of convergence three.

1.2. Some Characteristics of the Iterative Methods

To analyze an iterative method in greater depth, it is essential to understand certain concepts related to the mathematical notation used, the order of convergence, the efficiency index, the computational order of convergence, as well as the fundamental theorems and conjectures that support the correct formulation of the proposed new multipoint methods. Each of these aspects is detailed below.

Order of convergence

The speed at which a sequence $\left\{x_n\right\}$ approaches a solution $\alpha$ is quantified by the order of convergence p. Formally, the sequence $\left\{x_n\right\}$ is said to converge to $\alpha$ with order $p\ge 1$ and asymptotic error constant $C>0$ if:

$$\underset{n\to \infty }{lim}{\displaystyle \frac{|x_{n+1}-\alpha |}{|x_n{-\alpha |}^p}}=C.$$

This limit establishes an asymptotic relation that describes how rapidly the errors $e_n=x_n-\alpha$ decay as the number of iterations increases. Specifically:

• If $p=1$ and $0<C<1$, the convergence is linear,
• If $p=2$, the convergence is quadratic,
• If $p=3$, the convergence is cubic,
• For $p>3$, the method has higher-order convergence.

In practice, a high value of p implies faster convergence toward the root of $f\left(x\right)=0$, assuming the constant C remains reasonably small. However, higher-order methods often require more function or derivative evaluations per iteration, increasing computational cost.

The error equation of an iterative method can be expressed as:

$$e_{n+1}=C{e}_n^p+\mathcal{O}\left(e_n^{p+1}\right),C\in \mathbb{R},$$

where C is the asymptotic error constant and $\mathcal{O}\left(e_n^{p+1}\right)$ denotes higher-order terms that become negligible as n increases. This expression is central to the local convergence analysis of iterative schemes.

Numerical estimation of the order of convergence

Since the exact root $\alpha$ is typically unknown, practical estimation of the convergence order relies on approximate values of the iterates. Two widely used techniques are:

• The computational order of convergence (COC) defined in [3],
• The approximate computational order of convergence (ACOC) defined in [7].

These tools are commonly used in numerical experimentation to assess the performance of iterative schemes.

The classical estimate (COC) [3], assuming knowledge of the root $\alpha$, is given by:

$$p\approx \mathrm{COC}={\displaystyle \frac{ln\left|\frac{x_{n+1}-\alpha }{x_n-\alpha }\right|}{ln\left|\frac{x_n-\alpha }{x_{n-1}-\alpha }\right|}},n\ge 2.$$

When $\alpha$ is unknown, the ACOC [7] formula provides a root-free estimate using only the iterates:

$$p\approx \mathrm{ACOC}={\displaystyle \frac{ln\left|\frac{x_{n+1}-x_n}{x_n-x_{n-1}}\right|}{ln\left|\frac{x_n-x_{n-1}}{x_{n-1}-x_{n-2}}\right|}},n\ge 3.$$

(4)

This tool allows us to approximate the theoretical order of convergence p without requiring knowledge of the exact solution. Its reliability increases as the iterates approach the root, provided the errors remain sufficiently small to avoid numerical cancellation or round-off issues.

Efficiency Index of an Iterative Method

To assess the computational efficiency of an iterative method, one must consider its convergence order and the number of function evaluations per iteration. Ostrowski (1966) introduced the efficiency index I, defined as:

$$I=p^{1/d},$$

(5)

where p is the order of convergence and d is the total number of functional evaluations per iteration (including derivatives, if applicable). This index provides a comparative efficiency measure across methods with varying orders and computational demands.

More recently, the concept has been extended to the computational efficiency index (CEI) [8], which includes not only functional evaluations but also the products/quotients of the iterative method.

$$\mathrm{CEI}=p^{1/(d+op)},$$

(6)

where $op$ refers to the number of products/quotients produced in each iteration.

These indicators allow different iterative methods to be compared regarding convergence speed and total computational cost.

In this manuscript, Section 2 presents some known fourth-order iterative methods used in the numerical section for comparison. In Section 3, the new schemes are presented and their order of convergence is proven. Section 4 deals with the dynamical analysis of one of the proposed families of iterative schemes. The best method in terms of stability is compared in the numerical section with known schemes. Two academic examples and two applied problems confirm the theoretical results. With some conclusions and references, we conclude the manuscript.

2. Some Fourth-Order Methods in the Literature

In recent decades, there has been an urgent need to develop iterative methods with high orders of convergence that do not require new functional evaluations or derivatives. Since Traub’s initial contributions [9] with his method known as the frozen derivative, various approaches have been proposed to address this challenge. The following iterative expression defines Traub’s method:

$$y_n=x_n-{\displaystyle \frac{f\left(x_n\right)}{f^{\prime }\left(x_n\right)}},x_{n+1}=y_n-{\displaystyle \frac{f\left(y_n\right)}{f^{\prime }\left(x_n\right)}},n=0,1,2,\dots$$

This scheme achieves cubic convergence without the need to evaluate the second derivative. Similarly, Jarratt [10] introduced a two-step iterative scheme:

$$y_n=x_n-{\displaystyle \frac{2}{3}}{\displaystyle \frac{f\left(x_n\right)}{f^{\prime }\left(x_n\right)}},x_{n+1}=x_n-{\displaystyle \frac{1}{2}}\left({\displaystyle \frac{3f^{\prime }\left(y_n\right)+f^{\prime }\left(x_n\right)}{3f^{\prime }\left(y_n\right)-f^{\prime }\left(x_n\right)}}\right){\displaystyle \frac{f\left(x_n\right)}{f^{\prime }\left(x_n\right)}},n=0,1,2,\dots ,$$

This also avoids evaluating second derivatives and achieves fourth-order convergence.

Based on these ideas, many multipoint methods have been developed to achieve even higher convergence orders. The specialized literature, including the works of Chun [11], Ostrowski [12], King [13], among others, offers a wide range of fourth-order schemes based on the adjustment of parameters and weighting functions. These schemes will serve as benchmarks for comparison with the methods proposed in this work.

Among them, it is worth highlighting the family of fourth-order methods introduced by Artidiello [14]. It is based on a weighting function $H\left(\mu \right)$ that generalizes several known schemes. Its formulation follows a two-step scheme:

$$\begin{array}{cc}\hfill y_n& =x_n-{\displaystyle \frac{f\left(x_n\right)}{f^{\prime }\left(x_n\right)}},\hfill \\ \hfill x_{n+1}& =y_n-H\left({\mu }_n\right){\displaystyle \frac{f\left(y_n\right)}{f^{\prime }\left(x_n\right)}},n=0,1,2,\dots ,\hfill \end{array}$$

(7)

where $\mu ={\displaystyle \frac{f\left(y\right)}{f\left(x\right)}}$.

Theorem 1 

([14]). Let $f:I\subset \mathbb{R}\to \mathbb{R}$ be a sufficiently differentiable function on an open interval I containing a simple root α of $f\left(x\right)=0$, and $H:\mathbb{R}\to \mathbb{R}$ any sufficiently differentiable function satisfying:

$$H\left(0\right)=1,H^{\prime }\left(0\right)=2,\left|H^{\prime \prime }\left(0\right)\right|<\infty .$$

Then, for an initial estimate $x_0$ close enough to α, method (7) converges with order at least 4 and its error equation is:

$$e_{n+1}=\left(\left({\displaystyle \frac{5-H^{\prime \prime }\left(0\right)}{2}}\right)c_3^2-c_2{c}_3\right)c_2^2{e}_n^4+\mathcal{O}\left(e_n^5\right),$$

where $c_j={\displaystyle \frac{f^{\left(j\right)}\left(\alpha \right)}{j!f^{\prime }\left(\alpha \right)}},j=2,3,\dots ,$ and $e_n=x_n-\alpha$.

This theoretical framework not only unifies classical methods (such as those of Ostrowski or Chun) through specific choices of $H\left(\mu \right)$, but also provides a rigorous basis for designing new schemes. Table 1 summarizes several iterative methods that achieve fourth order thanks to the correct consideration of weighting functions.

Each of these weight functions iterative methods satisfies the conditions $H\left(0\right)=1,H^{\prime }\left(0\right)=2,\left|H^{\prime \prime }\left(0\right)\right|<\infty$, achieving the corresponding scheme a fouth-order convergence.

3. General Framework of Mean-Based Iterative Methods

As mentioned in the SubSection 1.1, all iterative methods exhibit cubic convergence. Consequently, this motivates the creation of a new family of iterative methods in which weight functions play a key role in increasing the order of convergence from three to four, without adding new functional evaluations.

We consider the following multipoint method:

$$x_{n+1}=x_n-H\left({\mu }_n\right){\displaystyle \frac{M_T(f\left(x_n\right),f\left(y_n\right))}{f^{\prime }\left(x_n\right)}},$$

(8)

where $y_n$ represents the Newton step, $M_T$ is an arbitrary mean applied to $f\left(x_n\right)$ and $f\left(y_n\right)$, and $H\left(\mu \right)$ is a weight function, with $\mu ={\displaystyle \frac{f\left(y\right)}{f\left(x\right)}}$.

By choosing different symmetric means $M_T$, such as: the arithmetic mean ($M_A$) and ($M_{A_y}$); the harmonic mean ($M_{H_y}$); the counterharmonic mean ($M_C$) and ($M_{c_y}$), we obtain several iterative schemes with different correction steps, as summarized in Table 2.

Theorem 2 

Let $f:I\subset \mathbb{R}\to \mathbb{R}$ be a sufficiently differentiable function on an open interval I, holding its simple root $\alpha \in I$, that is, $f\left(\alpha \right)=0$ and $f^{\prime }\left(\alpha \right)\ne 0$. The multipoint iterative method defined by a Newton step:

$$y_n=x_n-{\displaystyle \frac{f\left(x_n\right)}{f^{\prime }\left(x_n\right)}},$$

(9)

followed by a corrector step

$$x_{n+1}=x_n-H\left({\mu }_n\right){\displaystyle \frac{M_T(f\left(x_n\right),f\left(y_n\right))}{f^{\prime }\left(x_n\right)}},n=0,1,\dots ,$$

(10)

where:

• $M_T(f\left(x_n\right),f\left(y_n\right))$ is a symmetric bivariate function representing a mean ($M_A$, $M_{A_y}$, $M_{H_y}$, $M_C$ or $M_{C_y}$),
• $\mu ={\displaystyle \frac{f\left(y\right)}{f\left(x\right)}}$, is the variable of the weight funtion H,
• $H\left(\mu \right)$ is a weight function with a Taylor expansion around $\mu =0$:

  $$H\left(\mu \right)=H\left(0\right)+H^{\prime }\left(0\right)\mu +{\displaystyle \frac{1}{2}}{H}^{\prime \prime }\left(0\right){\mu }^2+{\displaystyle \frac{1}{6}}{H}^{\prime \prime \prime }\left(0\right){\mu }^3+O\left({\mu }^4\right).$$

If the coefficients $H\left(0\right),H^{\prime }\left(0\right),H^{\prime \prime }\left(0\right),H^{\prime \prime \prime }\left(0\right)$ satisfy specific conditions given in Table 3 then the method achieves fourth-order convergence.

All methods achieve fourth-order convergence, with the generalized error equation:

$$e_{n+1}=\left({\gamma }_1{c}_2^3+{\gamma }_2{c}_2{c}_3\right)e_n^4+O\left(e_n^5\right),$$

$$c_j={\displaystyle \frac{1}{j!}}{\displaystyle \frac{f^{\left(j\right)}\left(\alpha \right)}{f^{\prime }\left(\alpha \right)}},j=2,3,\dots ,$$

(11)

and $e_n=x_n-\alpha$ represents the iteration error. The constants ${\gamma }_1$ and ${\gamma }_2$ depend on the chosen mean and on the coefficients of the weight function $H\left(\mu \right)$.

Proof. 

Let $e_n=x_n-\alpha$ denote the error at the n-th iteration. Since f is sufficiently differentiable and $\alpha$ is a simple root, we can use Taylor expansions of $f\left(x_n\right)$ and $f^{\prime }\left(x_n\right)$ around $\alpha$. In terms of $e_n$, we have:

$$f\left(x_n\right)=f\left(\alpha \right)+f^{\prime }\left(\alpha \right)e_n+{\displaystyle \frac{f^{\prime \prime }\left(\alpha \right)}{2!}}{e}_n^2+{\displaystyle \frac{f^{\prime \prime \prime }\left(\alpha \right)}{3!}}{e}_n^3+{\displaystyle \frac{f^{\left(4\right)}\left(\alpha \right)}{4!}}{e}_n^4+O\left(e_n^5\right).$$

Since $f\left(\alpha \right)=0$, it follows that

$$f\left(x_n\right)=f^{\prime }\left(\alpha \right)e_n+{\displaystyle \frac{f^{\prime \prime }\left(\alpha \right)}{2!}}{e}_n^2+{\displaystyle \frac{f^{\prime \prime \prime }\left(\alpha \right)}{3!}}{e}_n^3+{\displaystyle \frac{f^{\left(4\right)}\left(\alpha \right)}{4!}}{e}_n^4+O\left(e_n^5\right).$$

We simplify the calculations using the expression of the constants (11)

$$f\left(x_n\right)=f^{\prime }\left(\alpha \right)\left(e_n+c_2{e}_n^2+c_3{e}_n^3+c_4{e}_n^4+O\left(e_n^5\right)\right).$$

In a similar way,

$$\begin{array}{cc}\hfill f^{\prime }\left(x_n\right)& =f^{\prime }\left(\alpha \right)+f^{\prime \prime }\left(\alpha \right)e_n+{\displaystyle \frac{f^{\prime \prime \prime }\left(\alpha \right)}{2}}{e}_n^2+{\displaystyle \frac{f^{\left(4\right)}\left(\alpha \right)}{6}}{e}_n^3+O\left(e_n^4\right),\hfill \\ \hfill & =f^{\prime }\left(\alpha \right)\left(1+2c_2{e}_n+3c_3{e}_n^2+4c_4{e}_n^3+O\left(e_n^5\right)\right).\hfill \end{array}$$

Therefore,

$${\displaystyle \frac{f\left(x_n\right)}{f^{\prime }\left(x_n\right)}}=e_n-c_2{e}_n^2+(2c_2^2-2c_3)e_n^3+(-4c_2^3+7c_2{c}_3-3c_4)e_n^4+O\left(e_n^5\right),$$

$$\begin{array}{cc}\hfill y_n& =x_n-\left(e_n-c_2{e}_n^2+(2c_2^2-2c_3)e_n^3+(-4c_2^3+7c_2{c}_3-3c_4)e_n^4+O\left(e_n^5\right)\right)\hfill \\ \hfill & =e_n+\alpha -e_n+c_2{e}_n^2-(2c_2^2-2c_3)e_n^3-(-4c_2^3+7c_2{c}_3-3c_4)e_n^4+O\left(e_n^5\right).\hfill \end{array}$$

So,

$$y_n-\alpha =c_2{e}_n^2-(2c_2^2-2c_3)e_n^3+(4c_2^3-7c_2{c}_3+3c_4)e_n^4+O\left(e_n^5\right).$$

Hence, the error in the approximation of $y_n$ becomes:

$$e_n^y=c_2{e}_n^2-(2c_2^2-2c_3)e_n^3+(4c_2^3-7c_2{c}_3+3c_4)e_n^4+O\left(e_n^5\right).$$

(12)

In a similar way, we have:

$$f\left(y_n\right)=f^{\prime }\left(\alpha \right)e_n^y+{\displaystyle \frac{f^{\prime \prime }\left(\alpha \right)}{2!}}{\left(e_n^y\right)}^2+{\displaystyle \frac{f^{\prime \prime \prime }\left(\alpha \right)}{3!}}{\left(e_n^y\right)}^3+{\displaystyle \frac{f^4\left(\alpha \right)}{4!}}{\left(e_n^y\right)}^4+O{\left(e_n^y\right)}^5,$$

we rewrite the expansion in terms of normalized coefficients (11):

$$f\left(y_n\right)=f^{\prime }\left(\alpha \right)\left(e_n^y+c_2{\left(e_n^y\right)}^2+c_3{\left(e_n^y\right)}^3+c_4{\left(e_n^y\right)}^4+O{\left(e_n^y\right)}^5\right).$$

In terms of $e_n$, we obtain:

$$f\left(y_n\right)=f^{\prime }\left(\alpha \right)\left(c_2{e}_n^2+(2c_3-2c_2^2)e_n^3+(5c_2^3-7c_2{c}_3+3c_4)e_n^4+O\left(e_n^5\right)\right).$$

(13)

Therefore,

$$f\left(x_n\right)+f\left(y_n\right)=f^{\prime }\left(\alpha \right)\left(e_n+2c_2{e}_n^2+\left(-2c_2^2+3c_3\right)e_n^3+\left(5c_2^3-7c_2{c}_3+4c_4\right)e_n^4+O{\left(e_n\right)}^5\right).$$

(14)

We calculate, by direct division, ${\mu }_n={\displaystyle \frac{f\left(y_n\right)}{f\left(x_n\right)}}$

$${\mu }_n={\displaystyle \frac{f\left(y_n\right)}{f\left(x_n\right)}}=c_2{e}_n+(-3c_2^2+2c_3)e_n^2+(8c_2^3-10c_2{c}_3+3c_4)e_n^3+O\left(e_n^4\right).$$

Since ${\mu }_n={\displaystyle \frac{f\left(y_n\right)}{f\left(x_n\right)}}\to 0$ as $n\to \infty$, we expand the weight function $H\left(\mu \right)$ in a Taylor series around $\mu =0$:

$$H\left(\mu \right)=H\left(0\right)+H^{\prime }\left(0\right)\mu +{\displaystyle \frac{1}{2}}{H}^{\prime \prime }\left(0\right){\mu }^2+{\displaystyle \frac{1}{6}}{H}^{\prime \prime \prime }\left(0\right){\mu }^3+O\left({\mu }^4\right).$$

Substituting the expression in terms of $e_n$, we obtain:

$$\begin{array}{cc}\hfill H\left({\mu }_n\right)=& H\left(0\right)+H^{\prime }\left(0\right)c_2{e}_n+\left({\displaystyle \frac{1}{2}}{c}_2^2(H^{\prime \prime }\left(0\right)-6H^{\prime }\left(0\right))+2c_3{H}^{\prime }\left(0\right)\right)e_n^2\hfill \\ \hfill & +\left(c_2^3\left(8H^{\prime }\left(0\right)-3H^{\prime \prime }\left(0\right)+{\displaystyle \frac{H^{\prime \prime \prime }\left(0\right)}{6}}\right)+2c_3{c}_2(H^{\prime \prime }\left(0\right)-5H^{\prime }\left(0\right))+3c_4{H}^{\prime }\left(0\right)\right)e_n^3+O\left(e_n^4\right)\hfill \end{array}$$

(15)

Now, we detail the expansion of $M_T(f\left(x_n\right),f\left(y_n\right))$ using the arithmetic mean $M_A$ of Table 2. For the rest of the means, all calculations are analogous.

We return to (14) and (13), to substitute them into the expression:

$${\displaystyle \frac{f\left(x_n\right)+f\left(y_n\right)}{2f^{\prime }\left(x_n\right)}}={\displaystyle \frac{e_n}{2}}-c_2^2{e}_n^3+\left({\displaystyle \frac{9c_2^3}{2}}-{\displaystyle \frac{7c_2{c}_3}{2}}\right)e_n^4+O\left(e_n^5\right),$$

multiplied by (15):

$$\begin{array}{cc}\hfill H\left({\mu }_n\right){\displaystyle \frac{f\left(x_n\right)+f\left(y_n\right)}{2f^{\prime }\left(x_n\right)}}=& {\displaystyle \frac{e_{n}H\left(0\right)}{2}}+{\displaystyle \frac{1}{2}}{c}_2{e}_n^2{H}^{\prime }\left(0\right)\hfill \\ \hfill & +\left({\displaystyle \frac{1}{4}}{c}_2^2(-4H\left(0\right)-6H^{\prime }\left(0\right)+H^{\prime \prime }\left(0\right))+c_3{H}^{\prime }\left(0\right)\right)e_n^3\hfill \\ \hfill & +{\displaystyle \frac{1}{12}}\left(c_2^3(54H\left(0\right)+36H^{\prime }\left(0\right)-18H^{\prime \prime }\left(0\right)+H^{\prime \prime \prime }\left(0\right))\right.\hfill \\ \hfill & \left.-6c_3{c}_2(7H\left(0\right)+10H^{\prime }\left(0\right)-2H^{\prime \prime }\left(0\right))+18c_4{H}^{\prime }\left(0\right)\right)e_n^4\hfill \\ \hfill & +O\left(e_n^5\right)\hfill \end{array}$$

(16)

Thus, equation $M_A$ of Table 2 expands to:

$$\begin{array}{cc}\hfill x_{n+1}=& e_n+\alpha -{\displaystyle \frac{e_{n}H\left(0\right)}{2}}-{\displaystyle \frac{1}{2}}{c}_2{e}_n^2{H}^{\prime }\left(0\right)\hfill \\ \hfill & +\left({\displaystyle \frac{1}{4}}{c}_2^2(-4H\left(0\right)-6H^{\prime }\left(0\right)+H^{\prime \prime }\left(0\right))+c_3{H}^{\prime }\left(0\right)\right)e_n^3\hfill \\ \hfill & -{\displaystyle \frac{1}{12}}\left(c_2^3(54H\left(0\right)+36H^{\prime }\left(0\right)-18H^{\prime \prime }\left(0\right)+H^{\prime \prime \prime }\left(0\right))-6c_3{c}_2(7H\left(0\right)+10H^{\prime }\left(0\right)-2H^{\prime \prime }\left(0\right))\right.\hfill \\ \hfill & \left.+18c_4{H}^{\prime }\left(0\right)\right)e_n^4+O\left(e_n^5\right),\hfill \end{array}$$

(17)

that is,

$$\begin{array}{cc}\hfill x_{n+1}-\alpha =& \left(1-{\displaystyle \frac{H\left(0\right)}{2}}\right)e_n-{\displaystyle \frac{1}{2}}{c}_2{H}^{\prime }\left(0\right)e_n^2\hfill \\ \hfill & +\left(c_2^2\left(H\left(0\right)+{\displaystyle \frac{3}{2}}{H}^{\prime }\left(0\right)-{\displaystyle \frac{1}{4}}{H}^{\prime \prime }\left(0\right)\right)-c_3{H}^{\prime }\left(0\right)\right)e_n^3\hfill \\ \hfill & +{\displaystyle \frac{1}{12}}\left(-c_2^3(54H\left(0\right)+36H^{\prime }\left(0\right)-18H^{\prime \prime }\left(0\right)+H^{\prime \prime \prime }\left(0\right))\right.\hfill \\ \hfill & \left.+6c_3{c}_2(7H\left(0\right)+10H^{\prime }\left(0\right)-2H^{\prime \prime }\left(0\right))-18c_4{H}^{\prime }\left(0\right)\right)e_n^4\hfill & \hfill +O\left(e_n^5\right).\end{array}$$

(18)

By solving the system obtained from eliminating the first-, second, and third-order error terms, ensuring that:

$$H\left(0\right)=2,H^{\prime }\left(0\right)=0,H^{\prime \prime }\left(0\right)=8,\left|H^{\prime \prime \prime }\left(0\right)\right|<\infty ,$$

the error of the method based on the arithmetic mean $M_A$ becomes:

$$e_{n+1}={\displaystyle \frac{1}{12}}\left((36-H^{\prime \prime \prime }\left(0\right))c_2^3-12c_2{c}_3\right)e_n^4+O\left(e_n^5\right).$$

This finishes the proof for the case of the arithmetic mean. The order of convergence of the remaining methods is obtained in a similar way, replacing the mean function $M_T$ and using the corresponding values of the coefficients of $H\left(\mu \right)$ presented in Table 3. Proceeding in this manner, the methods indicated in Table 2 achieve a fourth-order convergence. □

4. Dynamical Analysis

The order of convergence of an iterative method is not the only relevant criterion when evaluating its performance. In fact, the dynamical of the method, that is, the behavior of its orbits under different initial estimations, plays a fundamental role in its overall analysis. For it, tools from complex analysis are used, which represent the evolution of the methods in the Riemann sphere $\mathbb{C}\cup \{\infty \}$ [17,18,19].

We start from a rational function resulting from the application of an iterative method to a polynomial of low degree, denoted by $R:\widehat{\mathbb{C}}\to \widehat{\mathbb{C}}$, where $\widehat{\mathbb{C}}=\mathbb{C}\cup \{\infty \}$ denotes the Riemann sphere. The orbit of a point $z_0\in \widehat{\mathbb{C}}$ is given by the sequence:

$$\{z_0,R\left(z_0\right),R^2\left(z_0\right),\dots ,R^n\left(z_0\right),\dots \}.$$

We are interested in studying the asymptotic behavior of the orbits, so we must classify the different points of the rational operator R. A point $\widehat{z}$ is k-periodic $k\ge 1$, if

$$R^k\left(\widehat{z}\right)=\widehat{z},\mathrm{and}{R}^p\left(\widehat{z}\right)\ne \widehat{z},\forall 1\le p<k.$$

We say that it is a fixed point of R if

$$R\left(\widehat{z}\right)=\widehat{z}.$$

If this fixed point is not a solution of the polynomial, it is called a strange fixed point, as well as being numerically undesirable, since the iterative method can converge on them under certain initial guesses [18].

The dynamical behavior of these fixed points $z^{*}$ is classified according to the modulus of the derivative $R^{\prime }\left(z^{*}\right)$. If $|R^{\prime }\left(z^{*}\right)|<1$, the fixed point is an attractor; if $|R^{\prime }\left(z^{*}\right)|=1$, it is called parabolic or indifferent; if $|R^{\prime }\left(z^{*}\right)|>1$, it is repellent; and if $|R^{\prime }\left(z^{*}\right)|=0$, it is a superattractor.

On the other hand, the study of the basin of attraction of an attractor $z^{*}$ is defined as the set of points that converge to it

$$A\left(z^{*}\right)=\{z_0\in \widehat{\mathbb{C}}:\underset{n\to \infty }{lim}{R}^n\left(z_0\right)=z^{*}\}.$$

The Fatou set F is the union of the basins of attraction. The Julia set J is its topological complement in the Riemann sphere and represents the union of the frontiers of the basins of attraction.

The following classic result, given by Fatou [20] and Julia [21], includes both periodic points (of any period) and fixed points, considered as periodic points of unit period.

Theorem 3 

Let R be a rational function. The immediate basins of attraction of each attractive periodic point contain at least one critical point.

Using this key result, the entire attraction behavior can be found using the critical points as seeds of the iterative process [22]. A point $z^{*}\in \widehat{\mathbb{C}}$ is critical for R if:

$$R^{\prime }\left(z^{*}\right)=0.$$

We are going to start with a rational function resulting from the application of an iterative scheme on a quadratic polynomial. In order to obtain global results for the class of quadratic polynomial we prove a Scaling Theorem for the corresponding iterative method.

4.1. Conjugacy Classes

Let f and g be two analytic functions defined on the Riemann sphere. An analytic conjugacy between f and g is a diffeomorphism h on the Riemann sphere such that $h\circ f\circ h^{-1}=g$.

We now state a general theorem that applies to all types of symmetric means described in Table 2.

Theorem 4 

Let $f:\widehat{\mathbb{C}}\to \mathbb{C}$ be an analytic function on the Riemann sphere, and let $h\left(z\right)=\alpha z+\beta$ be an affine transformation with $\alpha \ne 0$, and $g\left(z\right)=\lambda f\left(h\right(z\left)\right)$, with $\lambda \ne 0$. Let us consider the iterative scheme defined by

$$G_f\left(z\right)=y_f-H\left({\mu }_f\right){\displaystyle \frac{M_T\left(f\left(z\right),f\left(y_f\right)\right)}{f^{\prime }\left(z\right)}},$$

(19)

where $y_f=z-{\displaystyle \frac{f\left(z\right)}{f^{\prime }\left(z\right)}}$ is Newton’s method, being $M_T$ one of the means that provide the schemes $M_{A_y}$, $M_A$, $M_{H_y}$ or $M_C$, $M_{C_y}$, Here, $H\left(\mu \right)$ satisfies the conditions indicated in Table 3.

Then, (19) is analytically conjugate to the analogous method applied to g, that is:

$$(h\circ G_g\circ h^{-1})\left(z\right)=G_f\left(z\right),$$

(20)

sharing the same essential dynamics.

Proof. 

To prove the general result, we consider a particular case of the mean. For the rest of the methods, the proof is analogous. We choose the case of $M_{A_y}$, whose scheme is given by:

$$x_{n+1}=y_n-H\left({\mu }_n\right){\displaystyle \frac{f\left(x_n\right)+f\left(y_n\right)}{2f^{\prime }\left(x_n\right)}}.$$

(21)

As can be seen, its structure is representative of the methods included in Table 2. We know that the affine function $h\left(z\right)=\alpha z+\beta$ has an inverse given by $h^{-1}\left(z\right)={\displaystyle \frac{z-\beta }{\alpha }}$.

By hypothesis, $g\left(z\right)=\lambda f\left(h\right(z\left)\right)=\lambda f(\alpha z+\beta )$. By the chain rule, we obtain:

$$g\left(h^{-1}\left(z\right)\right)=\lambda f\left(z\right),g^{\prime }\left(h^{-1}\left(z\right)\right)=\lambda \alpha f^{\prime }\left(z\right).$$

(22)

Defining the operator $G_g\left(z\right)$ as

$$G_g\left(z\right)=y_g-H\left({\mu }_g\right){\displaystyle \frac{M_T\left(g\left(z\right),g\left(y_g\right)\right)}{g^{\prime }\left(z\right)}},$$

(23)

and evaluated at $h^{-1}\left(z\right)$, we obtain:

$$\begin{array}{cc}\hfill G_g\left(h^{-1}\left(z\right)\right)& =h^{-1}\left(z\right)-{\displaystyle \frac{g\left(h^{-1}\left(z\right)\right)}{g^{\prime }\left(h^{-1}\left(z\right)\right)}}\hfill \\ \hfill & -H\left({\mu }_g\left(h^{-1}\left(z\right)\right)\right){\displaystyle \frac{g\left(h^{-1}\left(z\right)\right)+g\left(h^{-1}\left(z\right)-\frac{g\left(h^{-1}\left(z\right)\right)}{g^{\prime }\left(h^{-1}\left(z\right)\right)}\right)}{2g^{\prime }\left(h^{-1}\left(z\right)\right)}}.\hfill \end{array}$$

(24)

Using the identities from (22) and taking into account that

$$h\left(h^{-1}\left(z\right)\right)=z,h\left(h^{-1}\left(z\right)-{\displaystyle \frac{g\left(h^{-1}\left(z\right)\right)}{g^{\prime }\left(h^{-1}\left(z\right)\right)}}\right)=z-{\displaystyle \frac{f\left(z\right)}{f^{\prime }\left(z\right)}},$$

we deduce that:

$$g\left(h^{-1}\left(z\right)-{\displaystyle \frac{g\left(h^{-1}\left(z\right)\right)}{g^{\prime }\left(h^{-1}\left(z\right)\right)}}\right)=\lambda f\left(z-{\displaystyle \frac{f\left(z\right)}{f^{\prime }\left(z\right)}}\right).$$

(25)

Substituting (22) and (25) into (24), we obtain:

$$\begin{array}{ccc}\hfill G_g\left(h^{-1}\left(z\right)\right)& =& {\displaystyle \frac{z-\beta }{\alpha }}-{\displaystyle \frac{\lambda f\left(z\right)}{\lambda \alpha f^{\prime }\left(z\right)}}-H\left({\mu }_g\left(h^{-1}\left(z\right)\right)\right){\displaystyle \frac{\lambda f\left(z\right)+\lambda f\left(z-{\displaystyle \frac{f\left(z\right)}{f^{\prime }\left(z\right)}}\right)}{2\lambda \alpha f^{\prime }\left(z\right)}}\hfill \\ \hfill & =& {\displaystyle \frac{z-\beta }{\alpha }}-{\displaystyle \frac{f\left(z\right)}{\alpha f^{\prime }\left(z\right)}}-H\left({\mu }_g\left(h^{-1}\left(z\right)\right)\right){\displaystyle \frac{f\left(z\right)+f\left(z-{\displaystyle \frac{f\left(z\right)}{f^{\prime }\left(z\right)}}\right)}{2\alpha f^{\prime }\left(z\right)}}\hfill \end{array}$$

Now, we apply the transformation h:

$$\begin{array}{ccc}\hfill h\left(G_g\left(h^{-1}\left(z\right)\right)\right)& =& \alpha \left[{\displaystyle \frac{z-\beta }{\alpha }}-{\displaystyle \frac{f\left(z\right)}{\alpha f^{\prime }\left(z\right)}}-H\left({\mu }_g\left(h^{-1}\left(z\right)\right)\right){\displaystyle \frac{f\left(z\right)+f\left(z-{\displaystyle \frac{f\left(z\right)}{f^{\prime }\left(z\right)}}\right)}{2\alpha f^{\prime }\left(z\right)}}\right]+\beta \hfill \\ \hfill & =& z-{\displaystyle \frac{f\left(z\right)}{f^{\prime }\left(z\right)}}-H\left({\mu }_g\left(h^{-1}\left(z\right)\right)\right){\displaystyle \frac{f\left(z\right)+f\left(z-{\displaystyle \frac{f\left(z\right)}{f^{\prime }\left(z\right)}}\right)}{2f^{\prime }\left(z\right)}}\hfill \end{array}$$

We observe that the term ${\mu }_g\left(h^{-1}\left(z\right)\right)$ transforms as:

$${\mu }_g\left(h^{-1}\left(z\right)\right)={\displaystyle \frac{g\left(h^{-1}\left(z\right)-\frac{g\left(h^{-1}\left(z\right)\right)}{g^{\prime }\left(h^{-1}\left(z\right)\right)}\right)}{g\left(h^{-1}\left(z\right)\right)}}={\displaystyle \frac{f\left(z-\frac{f\left(z\right)}{f^{\prime }\left(z\right)}\right)}{f\left(z\right)}}={\mu }_f\left(z\right).$$

(26)

Substituting (26) into (), we finally deduce:

$$h\left(G_g\left(h^{-1}\left(z\right)\right)\right)=z-{\displaystyle \frac{f\left(z\right)}{f^{\prime }\left(z\right)}}-H\left({\mu }_f\left(z\right)\right){\displaystyle \frac{f\left(z\right)+f\left(z-\frac{f\left(z\right)}{f^{\prime }\left(z\right)}\right)}{2f^{\prime }\left(z\right)}}=G_f\left(z\right),$$

(27)

which proves the desired identity (20), and confirms that $G_f$ and $G_g$ are analytically conjugate through the affine transformation $h\left(z\right)$. □

The same reasoning extends directly to the methods in Table 2, since:

• In each case, the correction term maintains the form $M_T\left(f\left(z\right),f\left(y_f\right)\right)$, with symmetric combinations based on means.
• The affine transformation h acts compatibly on both $f\left(z\right)$ and $f\left(y_f\right)$, preserving the functional structure of the correction.
• The identity ${\mu }_g\left(h^{-1}\left(z\right)\right)={\mu }_f\left(z\right)$ holds, since it depends only on the ratio $f\left(y_f\right)/f\left(z\right)$, scaled by $\lambda$.

Therefore, the result holds for the entire family of iterative methods based on symmetric means, as described in Table 2.

4.2. Dynamics of Fourth-Order Methods

As shown in Table 2, five different parametric families of iterative methods are identified, each of them associated with specific conditions on the function $H\left(\mu \right)$. To satisfy these conditions, polynomial weight functions have been selected. However, other types of functions could also be considered, as long as they meet the imposed constraints. This also allows the introduction of an additional parameter $\beta$ in cases where $H\left(\mu \right)$ is bounded.

Table 4 presents the polynomials chosen by the authors for the development of this paper.

Here, ${\mu }_n={\displaystyle \frac{f\left(y_n\right)}{f\left(x_n\right)}}$ and $\beta$ is a free parameter.

To observe the dynamics of these iterative methods, we will take as an example the arithmetic mean family $M_{A_y}$, which is defined as:

$$y_n=x_n-{\displaystyle \frac{f\left(x_n\right)}{f^{\prime }\left(x_n\right)}},$$

$$x_{n+1}=y_n-\left(2{\displaystyle \frac{f\left(x_n\right)}{f\left(y_n\right)}}-2{\left({\displaystyle \frac{f\left(x_n\right)}{f\left(y_n\right)}}\right)}^2-\beta {\left({\displaystyle \frac{f\left(x_n\right)}{f\left(y_n\right)}}\right)}^3\right){\displaystyle \frac{f\left(x_n\right)+f\left(y_n\right)}{2f^{\prime }\left(x_n\right)}}.$$

(28)

The other cases can be analyzed in a similar way.

4.3. Rational Operator

Proposition 1. 

Let us consider the generic quadratic polynomial $p\left(z\right)=(z-a)(z-b)$, of roots a and b. The rational operator related to family $M_{A_y}$ given in (28) on $p\left(z\right)$ is:

$$R_p(z,\beta )=-{\displaystyle \frac{z^4\left(-\beta +2z^6+16z^5+54z^4+96z^3-\beta z^2+92z^2-3\beta z+44z+8\right)}{\beta z^6-8z^6+3\beta z^5-44z^5+\beta z^4-92z^4-96z^3-54z^2-16z-2}},$$

(29)

being $\beta \in \mathbb{C}$ an arbitrary parameter.

Proof. 

We apply the iterative scheme $M_{A_y}$ to $p\left(z\right)$ and obtain a rational function $A_p(z,\beta )$ that depends on the roots $a,b$ and the parameter $\beta \in \mathbb{C}$. Then, we apply a Möbius transformation [19,23,24] on $A_p(z,\beta )$ with

$$h\left(z\right)={\displaystyle \frac{z-a}{z-b}},$$

which satisfies $h\left(a\right)=0,h\left(b\right)=\infty$ and $h(\infty )=1$. This transformation maps the roots a and b to the points 0 and ∞, respectively, whose nature is attractive, and the divergence of the method to 1. Thus, the new conjugate rational operator is defined as:

$$\begin{array}{cc}\hfill R_p(z,\beta )& :=h\circ A_p(z,\beta )\circ h^{-1}\left(z\right),\hfill \end{array}$$

(30)

$$\begin{array}{cc}\hfill & =-{\displaystyle \frac{z^4\left(-\beta +2z^6+16z^5+54z^4+96z^3-beta{z}^2+92z^2-3\beta z+44z+8\right)}{\beta z^6-8z^6+3\beta z^5-44z^5+\beta z^4-92z^4-96z^3-54z^2-16z-2}}.\hfill \end{array}$$

(31)

which no longer depends on the parameters a and b. □

Thus, this transformation facilitates the analysis of the dynamics of iterative methods by allowing the standardization of roots and the structural study of dynamic planes and their stability regions [8].

4.4. Fixed Points of the Operator

Now, we calculate all the fixed points of $R_p(z,\beta )$, to subsequently analyze their character (attractive, repulsive, neutral, or parabolic). Taking into account that the method has order four, the points $z=0$ and $z=\infty$ are always superattractor fixed points, since they come from the roots of the polynomial.

The fixed points of $R_p(z,\beta )$ are $z=0$, $z=\infty$, and nine strange fixed points:

• $z=1$ is a strange fixed point if $\beta \ne {\displaystyle \frac{312}{5}}$,
• The roots of the polynomial:

  $$P_{\beta }\left(t\right)=2+18t+72t^2+(160+\beta )t^3+(208+3\beta )t^4+(160+\beta )t^5+72t^6+18t^7+2t^8,$$

  (32)

  denoted by $e{x}_i\left(\beta \right)$, $i=1,\dots ,8$, for any $\beta \in \mathbb{C}$.

Now, we study the stability of the strange fixed point $z=1$.

Proposition 2. 

The strange fixed point $z=1$ has the following character:

• If $\beta ={\displaystyle \frac{312}{5}}$, $z=1$ is not a fixed point.
• If $\left|\beta -{\displaystyle \frac{312}{5}}\right|>{\displaystyle \frac{1024}{5}}$, $z=1$ is attractor.
• If $\left|\beta -{\displaystyle \frac{312}{5}}\right|={\displaystyle \frac{1024}{5}}$, $z=1$ is parabolic.
• If $\left|\beta -{\displaystyle \frac{312}{5}}\right|<{\displaystyle \frac{1024}{5}}$, $z=1$ is repulsive.

Proof. 

As seen in the previous section, the behavior of the fixed point can be determined according to the value of the stability function: it will be an attractor if $|R_p^{\prime }(z,\beta )|<1$, a repulsor if $|R_p^{\prime }(z,\beta )|>1$, superattractor if $|R_p^{\prime }(z,\beta )|=0$ and parabolic if $|R_p^{\prime }(z,\beta )|=1$. The expression of operator $|R_p^{\prime }(z,\beta )|$ is

$$R_p^{\prime }(z,\beta )=-{\displaystyle \frac{2z^3{(z+1)}^8\left(4\beta +4\beta z^4-32z^4+7\beta z^3-156z^3-12\beta z^2-248z^2+7\beta z-156z-32\right)}{{\left(\beta z^6-8z^6+3\beta z^5-44z^5+\beta z^4-92z^4-96z^3-54z^2-16z-2\right)}^2}}.$$

Therefore,

$$R_p^{\prime }(1,\beta )={\displaystyle \frac{1024}{312-5\beta }}.$$

If $\beta ={\displaystyle \frac{312}{5}}$, then $z=1$ is not a fixed point. To determine whether it is attractive or repulsive, we solve:

$$\left|{\displaystyle \frac{1024}{312-5\beta }}\right|\le 1⟺{1024}^2\le {|312-5\beta |}^2.$$

Expressing the right side in terms of $\Re \left(\beta \right)$ and $\Im \left(\beta \right)$:

$${|312-5(\Re \left(\beta \right)+i\Im \left(\beta \right))|}^2={(312-5\Re \left(\beta \right))}^2+25\Im {\left(\beta \right)}^2.$$

So,

$${1024}^2\le {312}^2-3120\Re \left(\beta \right)+25\Re {\left(\beta \right)}^2+25\Im {\left(\beta \right)}^2$$

By simplifying, we get

$${\left(\Re \left(\beta \right)-{\displaystyle \frac{312}{5}}\right)}^2+\Im {\left(\beta \right)}^2\ge {\left({\displaystyle \frac{1024}{5}}\right)}^2,$$

thus,

$$\left|\beta -{\displaystyle \frac{312}{5}}\right|\ge {\displaystyle \frac{1024}{5}}.$$

□

Graphically, the behavior of the fixed point $z=1$ is visualized in Mathematica using the graph of the function $\Phi \left(\beta \right)=\left|{\displaystyle \frac{1024}{312-5\beta }}\right|$.

In Figure 1, the attraction zones are the yellow area and the repulsion zone corresponds to the gray area. For values of $\beta$ within the disk, $z=1$ is repulsive; while for values of $\beta$ outside the gray disk, $z=1$ becomes attractive. Therefore, it is natural to select values within the gray disk, since repulsive divergence improves the performance of the iterative scheme.

For the eighth roots $e{x}_i\left(\beta \right)$, $i=1,\dots ,8$, of polynomial $P_{\beta }\left(t\right)$, we obtain the following results:

• $|R_p^{\prime }\left(e{x}_1\left(\beta \right)\right)|=0$, there is no $\beta$ value,
• $|R_p^{\prime }\left(e{x}_2\left(\beta \right)\right)|=0$, for ${\beta }_1=0.408822$, and ${\beta }_2=18.0802$,
• $|R_p^{\prime }\left(e{x}_3\left(\beta \right)\right)|=0$, for ${\beta }_1=0.408822$, and ${\beta }_2=18.0802$,
• $|R_p^{\prime }\left(e{x}_4\left(\beta \right)\right)|=0$, for ${\beta }_{3,4}=0.782768\pm 1.1103i$,
• $|R_p^{\prime }\left(e{x}_5\left(\beta \right)\right)|=0$, there is no $\beta$ value,
• $|R_p^{\prime }\left(e{x}_6\left(\beta \right)\right)|=0$, for ${\beta }_{3,4}=0.782768\pm 1.1103i$,
• $|R_p^{\prime }\left(e{x}_7\left(\beta \right)\right)|=0$, for ${\beta }_5=-78.5858$,
• $|R_p^{\prime }\left(e{x}_8\left(\beta \right)\right)|=0$, for ${\beta }_5=-78.5858$.

In the following figure we represent the stability functions of the strange fixed points points $e{x}_i\left(\beta \right)$, $i=1,2,\dots ,8$.

For each root $e{x}_i\left(\beta \right)$ evaluated in the rational derivative operator (), its stability surfaces are constructed. In this context, the graphical representation distinguishes the orange regions as zones of attraction $|R_p^{\prime }\left(e{x}_i\left(\beta \right)\right)|<1$, the gray regions as zones of repulsion $|R_p^{\prime }\left(e{x}_i\left(\beta \right)\right)|>1$, superattraction zones when it is the vertex of the cone $|R_p^{\prime }\left(e{x}_i\left(\beta \right)\right)|=0$, and parabolic zones when it is at the boundary $|R_p^{\prime }\left(e{x}_i\left(\beta \right)\right)|=1$.

From Figure 2, the following conclusions are drawn:

• As the derivative operator associated with the strange fixed points $e{x}_{1,5}\left(\beta \right)$ can not be zero, it can be seen in Figure 2 that the resulting surface has only one gray region. This indicates that these fixed points are repulsive throughout the analyzed range, which is desirable, as it prevents convergence to a strange fixed point.
• Furthermore, at points $e{x}_{2,3}\left(\beta \right)$, we obtain we obtain ${\beta }_1=0.408822$ and ${\beta }_2=18.0802$. Figure 2 shows an inverted cone-shaped surface (normally yellow), representing an attractor inside the cone and a superattractor at its vertex ${\beta }_2$ (that of ${\beta }_1$ is similar, so it is omitted). The associated unstable domain is approximately $[0,0.2]\times [0,0.2]$, indicating a small but localized region. Similarly, by setting the derivative operator associated with the roots $e{x}_{4,6}\left(\beta \right)$ to zero, we obtain ${\beta }_3=0.782768+1.1103i$ and ${\beta }_4=0.782768-1.1103i$. Figure 2 and Figure 2 show behavior qualitatively similar to that of ${\beta }_2$, with a comparable domain.
• By setting the derivative operator associated with the roots $e{x}_{7,8}\left(\beta \right)$ to zero, we obtain ${\beta }_5=-78.5858$. As illustrated in Figure 2, a considerably wider region of attraction appears, approximately $(0,100)\times (0,100)$, indicating that the method shows marked instability for these values of $\beta$.

Therefore, to ensure the robustness of the method, values of $\beta$ where some root $e{x}_i\left(\beta \right)$ is an attractor or superattractor should be avoided. In contrast, values such as $\beta$ where all strange fixed points are repellers, are preferable to ensure stable numerical behavior.

Just as we have studied strange fixed points, we must also analyze critical points, since, recalling Theorem 3, it turns out that each attraction basin of an attractive periodic point (any period) contains at least one critical point.

4.5. Critical points of the operator

Proposition 3. 

The critical points of the rational operator $R_p(z,\beta )$ are $z=0$, $z=\infty$, directly related to the zeros of the polynomial, and the following free critical points:

$$z=-1,$$

$$z_{1,2}\left(\beta \right)=\pm {\displaystyle \frac{1}{2}}\sqrt{-a\left(\beta \right)+b\left(\beta \right)-c\left(\beta \right)}-d\left(\beta \right)-e\left(\beta \right),$$

$$z_{3,4}\left(\beta \right)=\pm {\displaystyle \frac{1}{2}}\sqrt{a\left(\beta \right)+b\left(\beta \right)-c\left(\beta \right)}+d\left(\beta \right)-e\left(\beta \right),$$

where the auxiliary functions $a\left(\beta \right)$, $b\left(\beta \right)$, $c\left(\beta \right)$, $d\left(\beta \right)$, and $e\left(\beta \right)$ are algebraic simplifications used for easy of notation:

$$\begin{array}{cc}\hfill a\left(\beta \right)& ={\displaystyle \frac{2(\beta -8)}{\sqrt{369{\beta }^2-1800\beta +784}}}\left({\displaystyle \frac{{(7\beta -156)}^3}{64{(\beta -8)}^3}}-{\displaystyle \frac{(-3\beta -62)(7\beta -156)}{{(\beta -8)}^2}}+{\displaystyle \frac{2(7\beta -156)}{\beta -8}}\right),\hfill \\ \hfill b\left(\beta \right)& ={\displaystyle \frac{{(7\beta -156)}^2}{32{(\beta -8)}^2}}-{\displaystyle \frac{-3\beta -62}{\beta -8}}-{\displaystyle \frac{8\beta -64}{4(\beta -8)}},\hfill \\ \hfill c\left(\beta \right)& ={\displaystyle \frac{\sqrt{369{\beta }^2-1800\beta +784}}{16(\beta -8)}},\hfill \\ \hfill d\left(\beta \right)& ={\displaystyle \frac{\sqrt{369{\beta }^2-1800\beta +784}}{16(\beta -8)}},\hfill \\ \hfill e\left(\beta \right)& ={\displaystyle \frac{7\beta -156}{16(\beta -8)}}.\hfill \end{array}$$

Thus, there are five free critical points, except for $\beta =0$, $\beta ={\displaystyle \frac{312}{5}}$, and $\beta =8$, where only three free critical points exist.

Proof. 

To prove the result, we recall that the derivative of the rational operator () is:

$$R_p^{\prime }(z,\beta )=-{\displaystyle \frac{2z^3{(z+1)}^8\left(4\beta +4\beta z^4-32z^4+7\beta z^3-156z^3-12\beta z^2-248z^2+7\beta z-156z-32\right)}{{\left(\beta z^6-8z^6+3\beta z^5-44z^5+\beta z^4-92z^4-96z^3-54z^2-16z-2\right)}^2}}$$

It is easily observed that its roots are $z=0$, $z=\infty$, $z=-1$, $z_{1,2}\left(\beta \right)$ and $z_{3,4}\left(\beta \right)$. These last four correspond to the roots of the polynomial of degree 4 in the numerator.

Now, let us observe that for certain values of $\beta$, only three free critical points exist. One such case is $\beta =0$, where the derivative of the operator simplifies to:

$${\displaystyle \frac{2z^3{(z+1)}^6\left(8z^2+23z+8\right)}{{(2z+1)}^2{\left(2z^3+6z^2+4z+1\right)}^2}}.$$

Here, the strange critical points are $z=-1$, and the conjugate pair

$$z={\displaystyle \frac{1}{16}}\left(-23\pm \sqrt{273}\right).$$

When $\beta ={\displaystyle \frac{312}{5}}$, the derivative operator becomes:

$$-{\displaystyle \frac{10x^3{(x+1)}^8\left(272x^2+895x+272\right)}{{\left(136x^5+494x^4+420x^3+180x^2+45x+5\right)}^2}}.$$

In this scenario, the strange critical points are $z=-1$, and the conjugate pair

$$z={\displaystyle \frac{1}{544}}\left(-895\pm 3\sqrt{56121}\right).$$

And finally, when $\beta =8$, the derivative operator becomes:

$${\displaystyle \frac{2x^4{(x+1)}^8\left(25x^2+86x+25\right)}{{\left(10x^5+42x^4+48x^3+27x^2+8x+1\right)}^2}},$$

whose zeros are $z=-1$ and the conjugate pair

$$z={\displaystyle \frac{1}{25}}\left(-43\pm 6\sqrt{34}\right).$$

□

For the free critical point $z=-1$, we have $R_p(-1,\beta )=1$, which is a strange fixed point. Therefore, the parameter plane associated with this critical point is not of much interest, since we already know the stability of $z=1$.

To visualize the behavior of the free critical points that depend on $\beta$, we plot the parameter planes. In each parameter plane, we use each free critical point as an initial estimation. A mesh of $2000\times 2000$ points is defined in the complex plane. Each point of the mesh corresponds to a value of $\beta$, that is an iterative method member to the family, and for each one of them, we iterate the rational function $R_p(z,\beta )$. If the orbit of the critical point converges to $z=0$ or $z=\infty$ in a maximum of 100 iterations, the point is colored red; otherwise, it is colored black.

As a first step, we graph the parameter plane of the conjugate pair $z_{1,2}\left(\beta \right)$, both in the domain $D_1=[-150,-25]\times [-225,225]$, which represents a broad stable performance region around the origin, and $D_2=[-150,25]\times [-60,60]$, which represents a divergence zone related to the $e{x}_i$.

It is observed that there are many values of the parameter $\beta$ for which the free critical points converge to the roots $z=0$ or $z=\infty$, visually, showing convergence in an approximate domain of $[-50,100]\times [-100,100]$. On the other hand, Figure 3 presents a divergence detail and refers to the strange fixed points $e{x}_7\left(\beta \right)$ and $e{x}_8\left(\beta \right)$, which, when computing their derivative operator set to zero, yielded a value of $\beta =-78.5858$, which precisely aligns with the divergence observed in Figure 2.

Likewise, the parameter plane of the conjugate pair $z_{3,4}\left(\beta \right)$ is shown, both in the domain $D_3=[-150,275]\times [-210,210]$, and a detail in the domain $D_4=[70,270]\times [-100,100]$, which represents the region that has not converged to any of the roots.

Figure 4 shows very stable behavior in an approximate domain of $[-50,50]\times [-100,100]$, while when the domain $D_4$ corresponding to Figure 4 is viewed in more detail, the mostly black region demonstrates a divergent behavior of the studied method (28).

4.6. Dynamical Planes

In the case of dynamical planes, each point in the complex plane is considered as a starting point $z_0$ of the iterative scheme and is represented with different colors depending on the point it converges to. In this case, points that converge to $z=\infty$ are colored blue, and those that converge to $z=0$ are colored orange. These dynamical planes have been generated using a grid of $800\times 800$ points and a maximum of 100 iterations per point. In these planes, the fixed points are represented by a white circle, the critical points by a white square, and the attracting points by an asterisk.

Next, the dynamical planes are plotted based on the values for $\beta$ obtained from the strange fixed points of the operator () and from the observations in the parameter plane.

In Figure 5 and Figure 6, good behavior of the method can be observed when choosing $\beta =0$ and $\beta =1$. The colors of the plots indicate convergence of the method to $z=0$ and $z=\infty$, which are the roots.

A notable case is $\beta ={\displaystyle \frac{312}{5}}$, since in Proposition Section 4.4 it was established that when $\beta$ takes that value, $z=1$ is not a fixed point, as observed in Figure 7.

In Figure 8, it can be clearly seen that $z=1$ is no longer characterized as a strange fixed point of the method. Moreover, we recall that when $\left|\beta -{\displaystyle \frac{312}{5}}\right|<{\displaystyle \frac{1024}{5}}$, $z=1$ is repulsive, as shown in Figure 5 and Figure 6, and when $\left|\beta -{\displaystyle \frac{312}{5}}\right|>{\displaystyle \frac{1024}{5}}$, $z=1$ is an attractor, as shown in Figure 9 and Figure 10.

Observe that when evaluating the dynamical plane in Figure 9, the method converges to $z=1$; a strange fixed point of the rational operator, even though the root $z=0$ is relatively closer. Furthermore, note the notable instability of choosing such a high value of $\beta$.

Based on the previous study (see Figure 2), when considering the value ${\beta }_5=-78.5858$, complex dynamical behavior was observed. In that figure, it is seen that the associated attraction cone covers a significantly larger area compared to other values of $\beta$. Additionally, this basin of attraction is related to the parameter plane of the conjugate critical point $z_{1,2}$. In this scenario, the method converges to basins different from the roots 0 and ∞, indicating that it is not suitable for root-finding. Therefore, values such as $\beta =-78.5858$ should be avoided when applying this method. Likewise, another example of divergence $\beta =200$.

Figure 11. $\beta =-78.5858$.

Figure 11. $\beta =-78.5858$.

Figure 12 represents the divergence observed in the parameter plane of Figure 4.

The analysis of the remaining iterative methods presented in Table 4 has been carried out similarly. The same qualitative information obtained in the previous analysis was also found in the rest of the families. The union of parameter planes is the same for all families. Therefore, their qualitative performance is the same, with a wide area of complex values for parameter $\beta$ giving rise to stable methods.

5. Numerical Examples

The iterative methods used in this section are presented in Table 4. This table considers all the conditions established in Table 3 with respect to $H_{\mu }$, aiming to guarantee fourth-order convergence. In particular, the parameter $\beta =0$ is selected, given its favorable behavior observed in the dynamical analysis.

To assess the efficiency of the newly proposed iterative methods, a comparison is made with classical algorithms from the literature, as listed in Table 1.

This section evaluates the performance of the newly proposed multipoint iterative methods compared to several well-established fourth-order methods. The comparison includes the following performance indicators:

• Number of iterations required for convergence (Iter),
• Approximate Computational Order of Convergence (ACOC),
• Bounds of the errors (Incr, Incr2), where

  $$Incr=|x_{n+1}-x_n|,Incr2=|f\left(x_{n+1}\right)|,$$
• Efficiency Index (EI),
• Execution time (Time), measured in seconds.

5.1. Academic Example 1: $f_1\left(x\right)=cos\left(x\right)-x$

The target root of the problem considered is $x^{*}\approx 0.7391$, obtained from the proposed nonlinear function. For the iterative process, $x_0=0$ was taken as the initial estimate, with a convergence tolerance of ${10}^{-50}$ and a maximum number of 100 iterations. The implemented algorithm uses as a stopping criterion.

$$|x_{n+1}-x_n|<{10}^{-50},or\left|f\left(x_{n+1}\right)\right|<{10}^{-50}.$$

If neither criterion is met, the procedure ends when the maximum number of iterations is reached. Table 5 presents the results obtained under these conditions.

The results reported in Table 5 underscore the strong competitiveness of the newly developed methods when compared with classical schemes. The proposed approaches exhibit comparable, and in several cases superior, performance in terms of both accuracy and computational efficiency. This is reflected not only in the reduced number of iterations required to approximate the root, but also in the exact convergence order of four achieved by several of the new methods. In particular, tiny residual errors were obtained, such as $|f\left(x_{n+1}\right)|\approx {10}^{-197}$ for method $M_{H_y}$ and exactly zero for method $M_{A_y}$.

Nevertheless, it is important to acknowledge the effective performance of the classical iterative schemes. For instance, the Jarratt’s method demonstrated remarkable robustness and efficiency, reaching errors of the order of ${10}^{-143}$ with relatively low computational cost.

On the other hand, the Artidiello IV method, despite its high estimated order of convergence (ACOC $=4.25$), produced significantly larger final errors ($\sim {10}^{-18}$). This behavior may indicate the presence of numerical instabilities or sensitivity to the transcendental nature of the problem.

5.2. Academic Example 2: $f_2\left(x\right)=e^x-cos\left(x\right)-2$

The real root of this function is $x^{*}\approx 0.9488$, with an initial estimate $x_0=0.5$, the same tolerance ${10}^{-50}$, the maximum number of 100 iterations, and the same stopping criterium as before. The results are given in Table 6.

In this nonlinear and rapidly varying function, $M_{A_y}$ exhibits the most outstanding precision, with final error on the order of ${10}^{-208}$, confirming its high robustness. The classical methods Jarratt, Ostrowski, and Artidiello IV maintain excellent convergence behavior with very low errors ($\sim {10}^{-176}$) and lower execution times.

Overall, methods $M_C$ and $M_{C_y}$ achieve superlinear convergence with low number of iterations, and their performance may improve with adaptive strategies. In contrast, methods such as $M_{A_y}$ offer a balance between precision and convergence, suggesting their suitability for problems demanding extremely high accuracy.

The proposed methods $M_C$, $M_{C_y}$, and $M_{H_y}$ demonstrate solid fourth-order behavior, but exhibit larger errors compared to Example 1, indicating sensitivity to the exponential component of $f_2\left(x\right)$. Among all methods, Kung-Traub’s performance is clearly limited in convergence order $ACOC\approx 3.0$, highlighting its theoretical constraint, though its execution time is the lowest.

Applied Problems

Problem 1: Chemical Equilibrium in Ammonia Synthesis

The analysis of chemical equilibrium systems using numerical methods has been widely addressed in the scientific literature. Solving complex nonlinear equations that model fractional conversions in reactive processes such as ammonia synthesis requires robust and efficient techniques.

This work analyzes the chemical equilibrium corresponding to the ammonia synthesis reaction from nitrogen and hydrogen, in a molar ratio of 1:3 [25], under standard industrial conditions (500 °C and 250 atm). The equation that describes this reaction is the following:

$$f\left(x\right)=x^4-7.79075x^3+14.7445x^2+2.511x-1.674,$$

(33)

where x is the fractional conversion. Of the four real roots of this equation, only one $x^{*}\approx 0.27776$ lies within the physical interval $[0,1]$, and therefore it is the only one with chemical significance.

In Table 7, it can be seen that all methods converge rapidly to the physically meaningful root, demonstrating high efficiency for this type of problem. For this results, we can used $x_0=0.5$, and the same stopping criterium as in the previous examples.

Method $M_{A_y}$ exhibits the best performance in terms of accuracy, reaching an Incr2 on the order of ${10}^{-208}$, positioning it as the most precise among the set, at the cost of one additional iteration.

Among the classical methods, Ostrowski and Zhao et al. stand out due to their excellent accuracy, with very small errors ($\sim {10}^{-138}$ and $\sim {10}^{-145}$, respectively) and theoretically ideal ACOC. The Artidiello IV method also shows good performance, achieving an Incr2 as low as ${10}^{-161}$, although its practical order of convergence falls slightly below the expected value.

It is worth noting that the Kung-Traub method is the only one in the set whose practical order of convergence approaches 3 (ACOC $\approx 2.994$), a limitation reflected in the relative error and intrinsic efficiency (EI). Nevertheless, it exhibits competitive computational time, being faster than some more precise methods.

Regarding the proposed methods $M_A$, $M_C$, $M_{C_y}$, and $M_{H_y}$, a consistent stability in the order of convergence and a good approximation in the errors can be observed, with results close to the classical methods, though without systematically surpassing them. Their computational performance is acceptable, albeit slightly inferior in terms of time.

In summary, for this chemical equilibrium problem, the methods with the best overall performance considering accuracy, efficiency, and stability are $M_{A_y}$, Ostrowski, and Zhao et al., all offering an ideal combination of minimal errors, fulfilled theoretical order of convergence, and low execution times.

Table 8. Comparison of iterative methods with initial point $x_0=0$.

Method	Root	Iter	Incr	Incr2	ACOC	Time (s)
$M_A$	-0.3841	10	$3.656\times {10}^{-31}$	$4.79\times {10}^{-120}$	3.997	0.1931
$M_{A_y}$	-0.3841	9	$5.865\times {10}^{-157}$	$3.893\times {10}^{-208}$	4.000	0.1948
$M_{H_y}$	-0.3841	9	$6.609\times {10}^{-35}$	$1.24\times {10}^{-134}$	3.998	0.1766
$M_C$	-0.3841	11	$4.265\times {10}^{-16}$	$1.52\times {10}^{-59}$	3.928	0.2086
$M_{C_y}$	-0.3841	13	$1.081\times {10}^{-36}$	$7.571\times {10}^{-142}$	3.998	0.2341
Jarratt	0.2778	4	$3.431\times {10}^{-16}$	$2.049\times {10}^{-61}$	3.906	0.08466
Ostrowski	0.2778	4	$3.157\times {10}^{-16}$	$1.459\times {10}^{-61}$	3.997	0.0515
King ($\beta =-1$)	-0.3841	19	$1.563\times {10}^{-33}$	$1.602\times {10}^{-129}$	3.998	0.293
Kung-Traub	0.2778	6	$1.707\times {10}^{-51}$	$8.41\times {10}^{-152}$	3.000	0.09162
Zhao et al.	0.2778	5	$1.077\times {10}^{-34}$	$8.707\times {10}^{-136}$	4.002	0.09587
Chun	-0.3841	73	$1.426\times {10}^{-42}$	$1.898\times {10}^{-165}$	3.999	0.8783
Artidiello IV	0.2778	5	$1.516\times {10}^{-39}$	$9.091\times {10}^{-156}$	3.960	0.08355

Table 8. Comparison of iterative methods with initial point $x_0=0$.

Method	Root	Iter	Incr	Incr2	ACOC	Time (s)
$M_A$	-0.3841	10	$3.656\times {10}^{-31}$	$4.79\times {10}^{-120}$	3.997	0.1931
$M_{A_y}$	-0.3841	9	$5.865\times {10}^{-157}$	$3.893\times {10}^{-208}$	4.000	0.1948
$M_{H_y}$	-0.3841	9	$6.609\times {10}^{-35}$	$1.24\times {10}^{-134}$	3.998	0.1766
$M_C$	-0.3841	11	$4.265\times {10}^{-16}$	$1.52\times {10}^{-59}$	3.928	0.2086
$M_{C_y}$	-0.3841	13	$1.081\times {10}^{-36}$	$7.571\times {10}^{-142}$	3.998	0.2341
Jarratt	0.2778	4	$3.431\times {10}^{-16}$	$2.049\times {10}^{-61}$	3.906	0.08466
Ostrowski	0.2778	4	$3.157\times {10}^{-16}$	$1.459\times {10}^{-61}$	3.997	0.0515
King ($\beta =-1$)	-0.3841	19	$1.563\times {10}^{-33}$	$1.602\times {10}^{-129}$	3.998	0.293
Kung-Traub	0.2778	6	$1.707\times {10}^{-51}$	$8.41\times {10}^{-152}$	3.000	0.09162
Zhao et al.	0.2778	5	$1.077\times {10}^{-34}$	$8.707\times {10}^{-136}$	4.002	0.09587
Chun	-0.3841	73	$1.426\times {10}^{-42}$	$1.898\times {10}^{-165}$	3.999	0.8783
Artidiello IV	0.2778	5	$1.516\times {10}^{-39}$	$9.091\times {10}^{-156}$	3.960	0.08355

Using $x_0=0$ as initial estimation, a bifurcation in the convergence of the methods is observed. In particular, all mean-base d methods $M_A$, $M_{A_y}$, $M_{H_y}$, $M_C$, and $M_{C_y}$ converge to the root $x^{*}\approx -0.384094$, which lacks physical interpretation within the context of the chemical equilibrium problem. Although this root is mathematically valid, it does not represent the desired solution, indicating a strong dependence of these methods on the initial condition and their sensitivity to the local behavior of the function.

The case of the King and Chun methods is particularly notable: although they also converge to the same non-physical root, they do so after a large number of iterations (19 and 73, respectively), reflecting lower numerical efficiency compared to the other methods.

On the other hand, classical methods such as Jarratt, Ostrowski, Zhao et al., Kung-Traub, and Artidiello IV correctly converge to the physically meaningful root $x^{*}\approx 0.27776$, requiring approximately 4 to 6 iterations and yielding extremely low final errors.

This highlights the importance of appropriately selecting the initial condition in problems with multiple solutions.

Problem 2: Determination of the Maximum in Planck’s Radiation Law

The study of blackbody radiation through numerical methods has been fundamental in the development of quantum physics. As noted in [26] in their work “Didactic Procedure for Teaching Planck’s Formula”, determining the spectral maximum in Planck’s distribution requires advanced techniques to solve nonlinear transcendental equations.

We analyze the equation derived from Planck’s radiation law that determines the wavelength corresponding to the maximum energy density:

$$f_5\left(x\right)=e^{-x}-1+{\displaystyle \frac{x}{5}},$$

where $x={\displaystyle \frac{hc}{\lambda kT}}$. Among the possible solutions, only $x\approx 4.9651142317$ has physical meaning in this context [27].

Table 9 shows that all proposed methods $M_A,M_{A_y},M_{H_y},M_C,M_{C_y}$ converge correctly to the physically valid root $x^{*}\approx 4.9651142317$, even starting from a distant initial condition ($x_0=1$).

Method $M_{A_y}$ stands out for its quartic convergence (ACOC = 4.0) and a final error of exactly 0 in only 5 iterations, making it the most accurate method. Although it is slightly more computationally expensive, it achieves the highest relative efficiency (EI ≈ 1.587) among the proposed methods.

Methods $M_A,M_C,M_{C_y}$ show quadratic convergence $ACOC\approx 2$, with errors of order ${10}^{-67}$ in only 4 iterations, with low execution times. The method $M_{H_y}$ slightly improves the order of convergence $ACOC\approx 2.67$, which represents a good compromise between efficiency and robustness.

In contrast, several classical methods converge to values that do not represent the physically significant root. For example, Chun converges to a negative value ($x\approx -158.3$) with a clear divergence (Incr2 $\sim {10}^{68}$). Meanwhile, the other methods diverge completely, presenting values with no physical meaning.

This highlights that, under unfavorable initial conditions, the proposed methods remain stable, while the classical ones are sensitive. Therefore, for problems such as Planck’s, methods based on the mean offer a more robust and reliable alternative.

Table 10 shows that all evaluated methods successfully converged to the physically meaningful root of Planck’s equation, despite starting from a distant initial condition ($x_0=10$).

The $M_{A_y}$ method maintains its outstanding performance, with quartic convergence order (ACOC = 4.0), zero error in the second increment (Incr2 = 0), and maximum intrinsic efficiency ($\mathrm{EI}=1.587$). This confirms its high robustness and accuracy even under unfavorable scenarios.

The methods $M_A$, $M_{H_y}$, $M_C$, and $M_{C_y}$ exhibited similar behavior, with $ACOC\approx 3.45$ and convergence achieved in just 3 iterations, maintaining errors on the order of ${10}^{-17}$ and very low execution times (< 0.04 s).

Among the classical methods, King, Ostrowski, Zhao et al., and Jarratt also showed good efficiency and stable convergence. The Kung-Traub method, although achieving higher precision $Incr2\approx {10}^{-122}$, resulted in a lower ACOC (3.0) and reduced efficiency due to its higher cost per iteration.

In summary, all the analyzed methods proved to be efficient from a distant initial point, with $M_{A_y}$ standing out for its stability and performance.

6. Conclusions

This manuscrpt presents a new perspective on the design, analysis, and dynamical behavior of fourth-order multipoint iterative methods, constructed through convex combinations of classical means and parameterized weight functions. By extending the Newton-type scheme and incorporating arithmetic, harmonic, and contra-harmonic, a versatile family of optimal methods is developed, complying with the Kung–Traub conjecture by achieving order four with a minimal number of functional evaluations. Other means such as Heronian, and centroidal have been used without positive results. The resulting iterative methods do not reach order four, so the are not optimal schemes.

The general formulation, grounded in a solid theoretical framework (Taylor expan- sions, affine conjugation, and local error analysis), enabled the derivation of explicit conditions on the weight functions to ensure fourth-order convergence. This was complemented by a rigorous dynamical systems analysis using tools such as conjugated rational operators, stability surfaces, parameter planes, and dynamical planes on the Riemann sphere.

The results reveal that the proposed parametric families, particularly those associated with the modified arithmetic mean ($M{A}_y$), exhibit stable convergence behavior over large regions of the complex parameter space $\beta$. Nevertheless, critical regions of divergence were identified, including attraction basins unrelated to the roots or convergence toward strange fixed points. Such zones often visualized as black or green regions in the parameter and dynamical planes must be avoided in practice. In this context, the detailed study of free critical points proved fundamental, as each attractive basin must contain at least one critical point. Their behavior provided early insights into the method’s stability and convergence characteristics. Noteworthy cases such as $\beta =-78.5858$ or $\beta =400$ illustrated convergence toward undesirable fixed points (e.g., $z=1$), despite proximity to the root $z=0$, emphasizing the importance of well-informed parameter selection.

Finally, numerical experiments confirmed the competitiveness of the proposed schemes. In various tests, both academic and applied, methods based on convex means showed efficient performance comparable to classical methods. In most of the problems considered, they managed to converge to the root in approximately 3 to 5 iterations. The only exception was the first applied example, which required a greater number of iterations; however, it also reached the root with errors of the order of $\sim {10}^{-208}$. Likewise, in the second applied problem, their robustness was evident, as they were the only methods that converged to the root, even starting from initial conditions far from the solution, unlike classical methods.

In summary, the proposed class of parametric iterative methods based on convex means not only achieves high computational efficiency but, when coupled with dynamical analysis, offers a robust and predictive framework for the stable and accurate solution of nonlinear equations.