Design and dynamical behavior of a fourth order family of iterative methods for solving nonlinear equations

1.   Introduction

An active topic in numerical analysis is the estimation of the solution of nonlinear equations and systems of equations by means of iterative processes. This activity results in a large variety of publications from different authors, such as Petković et al. ^[30] in the scalar case, Ahmad et al. ^[1], or Zhanlav et al. ^[36], among other researchers in vectorial problems.

The limited existence of analytical methods to find the solution compels us to use iterative methods to obtain an approximation. These schemes generate a sequence of real numbers that converges to a solution of the equation $f(x) = 0$, if it exists, being $f:I\subset\mathbb{R} \rightarrow \mathbb{R}$.

A scalar iterative method is defined by a fixed-point expression

where $g:\mathbb{R} \rightarrow \mathbb{R}$ is a fixed-point function. These methods can be classified as one-point or multipoint schemes. In the first ones, the $k+1$ iteration is obtained using only functional evaluations at the $k$-th iteration of $f$ and its derivatives. Meanwhile in multipoint methods, the $k+1$ iteration is obtained using functional evaluations at different points, which are taken as intermediate steps of the same iteration ^[34].

In this manuscript, by using the weight-function technique, we present a parametric family of optimal iterative schemes of fourth-order involving derivatives at different points. The stability of this class is studied in order to select its most stable elements and refuse the unstable ones.

There are problems in engineering, technology and in different fields of science that require the solution of nonlinear equations. For example, Qureshi et al. in ^[31] proposed a hybrid three-step iterative method for solving nonlinear equations in the field of medical science (blood rheology, population growth and neurophysiology). Yaseen et al. in ^[35] applied Jarratt-type methods to global positioning systems problems. Also, Chand et al. in ^[11] applied weight functions on several methods: Potra-Pták procedure and two iterative schemes (optimal and non-optimal, respectively) of fourth and sixth orders of convergence and a family of optimal eighth order methods, aiming to solve problems related to the effect of water flow, other factors in the flow of open channels (rivers or canals), and the determination of fluid flow through tubes and pipes ^[9].

Ostrowski in ^[29] defines the efficiency index of an iterative method as $I = p^{1/d}$, where $p$ is the convergence order and $d$ is the number of functional evaluations per iteration. According to Kung-Traub's conjecture ^[27]: The convergence order of an iterative method without memory cannot be greater than $2^{d-1}$ (called optimal order), where $d$ is the number of functional evaluations per iteration. As long as the condition established in the conjecture is satisfied, the method is considered optimal. For a deeper insight in optimal methods, we refer to ^[18].

Regarding the speed at which an iterative method reaches the zero of the equation $f(x) = 0$, it is possible to define a aproximated computational order of convergence (ACOC) ^[17], which considers the rate at which successive approximations approach the zero to estimate the order of convergence of the method:

On the other hand, the dynamical analysis has become a very useful tool for studying iterative methods (see the work of different authors and their colleagues as Behl in ^[6], Campos ^[10], Cordero ^[16], Khirallah ^[26], Kansal ^[25], Geum ^[22], Sayevand ^[32], or Scott in ^[33]) since it enables us to classify different iterative formulas, not only from the point of view of their order of convergence, but from the analysis of how these formulas behave depending on the chosen initial estimation. In addition, it provides valuable information about the stability and reliability of the iterative method. The dynamical analysis characterizes the stability of the fixed points of the rational function $R$ resulting from the application of the iterative expression on a low-degree polynomial $p(x)$. This study, joined with that of the existence of critical points, gives us a clear idea of the best elements of a class of iterative methods. Moreover, it is possible to obtain graphical representations that allow us to compare methods beyond mere numerical results. The dynamical plane allows us to visualize the stability of a method, the size of its convergence basins and the viability of certain initial estimations. Also, for uniparametric families of methods, the parameter plane will contribute to the choice of the most suitable family member.

Prior to performing the dynamical analysis of the family of iterative methods when applied over a low degree polynomial, we must remember a few ideas and definitions of complex dynamics ^[3,7]. Given a rational operator $R:\hat{\mathbb{C}}\rightarrow \hat{\mathbb{C}}$, where $\hat{ \mathbb{C}}$ is the Riemann sphere, the orbit of a point $x_{0}$ is defined as the sequence

A point $x^{*} \in \hat{ \mathbb{C}}$ is a $k$-periodic point of $R$ if it is kept invariant by the action of operator after $k$ iterations, that is, $R^k(x^{*}) = x^{*}$; if $k = 1$, it is called a fixed point. A periodic point is attracting if $\lvert R'(x^{*})\cdot (R^1(x^{*}))' \cdots (R^{k-1}(x^{*}))'\rvert < 1$, superattracting if $\lvert R'(x^{*})\cdot (R^1(x^{*}))' \cdots (R^{k-1}(x^{*}))' \rvert = 0$, repulsive if $\lvert R'(x^{*})\cdot (R^1(x^{*}))' \cdots (R^{k-1}(x^{*}))' \rvert > 1$ and neutral or parabolic if $\lvert R'(x^{*})\cdot (R^1(x^{*}))' \cdots (R^{k-1}(x^{*}))'\rvert = 1$. Furthermore, if $x^{*} \in \hat{ \mathbb{C}}$ is a fixed point, but not a zero of the polynomial $p(x)$, then it is called a strange fixed point.

On the other hand, a point $x_{c}$ is a critical point of $R$ if $R'(x_{c}) = 0$. Free critical points are those that do not coincide with the zeros of the polynomial.

The basin of attraction of an attracting fixed (or periodic) point $x^*$ is defined as

Moreover, the Fatou set of the rational function $R$ (see ^[21]), denoted by $\mathcal{F}(R)$, is the set of points $x^{*} \in \hat{ \mathbb{C}}$ whose orbits tend to an attractor of any kind. The Julia set ^[24], denoted by $\mathcal{J}(R)$, is the complementary of the Fatou set in $\hat{ \mathbb{C}}$, which includes the repulsive and neutral points and defines the boundaries between basins of attraction ^[12]. The key role of the critical points is stated in the following classical result (see, for example, ^[19]), which sets a relationship between the connected component of a basin of attraction including the attractive or periodic point (called immediate basin of attraction) and the critical point.

Theorem 1 (Fatou-Julia Theorem). Let $R$ be a rational function. The immediate basin of attraction of a periodic (or fixed) attractor point contains at least one critical point.

In order to perform the dynamical analysis, we use an arbitrary quadratic polynomial $p(x) = (x-a)(x-b)$, where $a, b \in \mathbb{C}$ and Möbius conjugacy map considered by Blanchard ^[8],

which satisfies the properties: $M(\infty) = 1, \; M(a) = 0, \; M(b) = \infty.$

It is easy to prove that the proposed class of iterative methods LCT, defined in Section 2, satisfies the scaling theorem (see, for example, ^[2]). So, the dynamical performance of the rational function obtained by applying the class of methods on quadratic polynomials is equivalent by conjugation. We then use the Möbius conjugacy map to get a conjugate rational function that does not depend on the roots of the polynomial $p(x)$ used. We discuss the discrete dynamical systems: We obtain the asymptotic performance of fixed points. In addition, we calculate the critical points and generate the parameter planes associated with the free critical points showing the stable and unstable regions, depending on the chosen parameter. These elements, therefore, allow us to determine which members of the class of iterative methods show more stable behavior.

This paper is organized as follows: In Section 2, we present the new class of multipoint methods and analyze its convergence. Using complex dynamics tools in Section 3, we study the dependence on initial estimates of the proposed family. These tools allow us to transform the analysis of the method into the study of a rational function obtained by applying the method over a generic quadratic polynomial. In Section 4, we perform some numerical tests, selecting parameter values within and outside the stability region of the parameter plane obtained in the previous section. Also, a comparison with other known methods is made. Finally, Section 5 is devoted to exhibit some conclusions.

2.   Development of a new class of iterative schemes

Chun in ^[15] stated that weight functions with a parameter can be successfully introduced into an iteration process to increase its order of the convergence. Also, Artidiello et al. proposed families of iterative schemes (see ^[4,5]) by using the weight function technique.

In this section, we present a new class of fourth-order multipoint Newton-like methods by means of a weight function $H_{m} (\mu)$, that we deonote by LCT:

where $H_{m} (\mu) = m_{1}+m_{2} \; \mu(x_{k}, y_{k}) +m_{3} \; \left[\mu(x_{k}, y_{k})\right]^{-1}+m_{4} \left[\mu(x_{k}, y_{k})\right] ^2$ is a real-valued weight function, $\mu(x_{k}, y_{k}) = \dfrac{f'(x_{k})}{f'(y_{k})}$, and $\theta$, $m_{i}$, $i = 1, 2, 3, 4$ are arbitrary parameters.

The next result states the conditions under which the class defined in (2.1) reaches order of convergence four. Hence, the methods that compose it are optimal in the sense of Kung-Traub's conjecture ^[27].

Theorem 2 (LCT family convergence order). Let $f:I \subset \mathbb{R} \rightarrow \mathbb{R}$ be a differentiable enough function defined in the open interval $I$. Let us also assume that $\xi$ is a simple zero of $f$ and $x_0$ is sufficiently close to $\xi$, then family (2.1), which we denote by LCT, has local order of convergence 4, as far as $m_{1} = \dfrac{5}{8}-m_{2}, \; m_{3} = \dfrac{m_{2}}{3}, \; m_{4} = \dfrac{3}{8}-\dfrac{m_{2}}{3}$, and $\theta = \frac{2}{3}$.

In this case, the error equation obtained is

where $C_{q} = \dfrac{1}{q!} \dfrac{f^{(q)}(\xi)}{f'(\xi)}$ for $q = 2, 3, ...$, and $e_{k} = x_{k}-\xi$.

Proof. Let us consider $\xi$ as a zero of $f(x)$, such that $f'(\xi)\neq 0$ and let $x_{k} = \xi+e_{k}$. Using a Taylor series expansion of $f(x_k)$ and $f'(x_k)$ around $x = \xi$, we get

and

respectively.

By direct division of (2.2) and (2.3), we obtain

Since $x_{k} = \xi+e_{k}$, the error at the first step is

Thus,

and

Now, let us consider the $\mu = \mu(x_k, y_k)$ function as

We propose the real-valued weight function depending of parameters $m_{1}$, $m_{2}$, $m_{3}$ and $m_{4}$ as

Combining (2.1) and (2.9), we obtain the error equation

In order to attain a class of iterative schemes with fourth-order of convergence, the coefficients of $e_{k}$, $e_{k}^2$ and $e_{k}^3$ in (2.10) must be simultaneously zero. We then solve the resulting system of equations to obtain the parametric solution:

By replacing (2.11) in (2.10), we have a family of fourth-order optimal methods, with error equation

□

We have proven that the convergence order of the proposed class is 4 and it requires three functional evaluations per iteration: $f(x_k), f'(x_k)$ and $f'(y_k)$. Therefore, $d = 3$ and the proposed family of methods is optimal under Kung-Traub's conjecture ^[27], with an efficiency index of $I = p^{1/d} = 4^{1/3} \approx 1.5874$.

In the next section, we analyze the behavior of the class of iterative schemes using complex dynamics tools, which transforms the study of the family into the analysis of its associated rational function when the class is applied on a generic quadratic polynomial.

3.   Stability of the parametric family

The proposed class of iterative schemes (2.11) holds a great amount of elements, as the only restrictions are those required in Theorem 2. Selecting a member requires to choose a specific weight function satisfying those conditions and, among all the infinite members that this weight function defines, all the schemes have the same order of convergence. The next question to be answered is how to choose the final members in order to be applied on a problem. A key characteristic of iterative methods, as important as the order of convergence, is the stability. This fact allows us to select the most appropriate members to be used, as the notion of stability involves the dependence on the initial estimation.

The usual reasons for the lack of convergence of an iterative method are divergence (convergence to infinity), a bad initial estimation and the presence of attracting elements (different to the solutions) that makes the iterate generate a periodic orbit that never converges. So, fixed (and periodic points) are the object of this study. As the most simple nonlinear function is a quadratic polynomial, we make the stability analysis on a generic one. The obtained results can not be directly extrapolated to any nonlinear function, but in practice the bad performance on quadratic polynomials yields to the worst performance on other functions. Moreover, the number of attracting elements is defined by the amount of free critical points (see Theorem 1). Therefore, the quantity of critical points and the basins of attraction that they hold are necessary to completely understand the performance of the different elements of the class.

In the family of methods proposed in (2.1), let us denote $m_{2} = \alpha$. Using conditions defined in (2.11), we have

Applying the LCT family on the quadratic polynomial $p(x) = (x-a)(x-b)$, where $a, b \in \mathbb{C}$ and applying the Möbius conjugacy map, we obtain the rational operator $Lg_{p, \alpha}(x)$,

with $\alpha \in \mathbb{C}$ as arbitrary parameter. The rational operator associated to that family of methods, using Möbius conjugacy map, is

depending on $\alpha$.

Factor $x^4$ obtained in (3.3) shows that the family of methods, for any parameter $\alpha$, has order of convergence 4 (at least over quadratic polynomials) and $R_{\alpha}(x)$ does not depend on parameters $a$ and $b$, as it satisfies the scaling theorem.

3.1.   Fixed points: calculation and stability analysis

Previously, we obtained the rational operator and with it, we now get the fixed points of operator $R_\alpha(x)$. A fixed point is one that is kept invariant by the operator, therefore, all the roots of $R_{\alpha}(x) = x$ are naturally fixed points. The existence of strange fixed points is not desirable (from the numerical point of view) or, if they exist, the best situation is that they are repulsive. The performance of an iterative method may be compromised if there are strange fixed points or periodic orbits of its associate rational function that are attracting.

Campos et al. ^[10] studied and characterized the stability of $z = 1$ as strange fixed point of any multipoint Newton-like methods (those where the intermediate evaluations are variations of Newton's scheme), as well as the character of point $z = -1$. They proved the hypothesis that makes this point a strange fixed or critical point, being pre-image of $z = 1$. These results involve the first item of the next result (the existence of $z = 1$ as strange fixed point), and also its stability, of which we provide in Theorem 4 an alternative proof. The rest of the items in Theorem 3 or in the rest of the dynamical analysis cannot be deduced from any result appearing in ^[10].

Theorem 3. The zeros of the equation $R_{\alpha}(x) = x$ are the fixed points of the rational function $R_{\alpha}(x)$, that is, $x = 0$, $x = \infty$, and also the (strange) fixed points:

i). $ex_{0}(\alpha) = 1$, whether $\alpha \neq -\frac{153}{8}$.

ii). $ex_{i}(\alpha)$ (for $i = 1, 2, \ldots, 6$), which corresponds to the 6 zeros of the polynomial

Thus, in general there are 9 fixed points, of which $x = 0$ and $x = \infty$ are related to the zeros of $p(x)$. However, this amount of fixed points is reduced for $\alpha = \frac{423}{8}$, as it can be checked that $ex_{1}\left(\frac{423}{8}\right) = ex_{2}\left(\frac{423}{8}\right) = 1 = ex_0\left(\frac{423}{8}\right)$, and when $\alpha = -\frac{153}{8}$, the rational function is simplified and $x = 1$ is not a fixed point.

Proof. Using the definition of fixed point, we solve

Simplifying and factoring the previous equation, we get

It can be checked in a similar way that $x = 0$ is a fixed point of $1/R_\alpha(1/x)$, so $x = \infty$ is a fixed point of $R_\alpha$. □

The derivative of the rational operator associated to the LCT family of methods, which is used to analyze the stability of fixed points and to obtain critical points, is:

with a free parameter $\alpha \in \mathbb{C}$.

In order to perform the stability analysis of the fixed points, we evaluate the derivative operator in each of them. As $\lvert R_{\alpha} '(0) \rvert = \left| \frac{1}{R_{\alpha} '(1/0)} \right| = 0$, $x = 0$ and $x = \infty$ are superattracting fixed points. Nevertheless, the stability of strange fixed points changes according to the values of parameter $\alpha$.

Theorem 4 (Stability of $ex_{0}\left(\alpha \right) = 1$). The asymptotical behavior of the strange fixed point of $R_\alpha(x)$, $ex_{0}\left(\alpha \right) = 1$ (for $\alpha \neq -\dfrac{153}{8}$), obtained by evaluating $\lvert R_{\alpha}' (1) \rvert$, is:

i). It is an attractor for $\Bigg\lvert \alpha +\dfrac{153}{8} \Bigg\rvert < 72$, but it cannot be superattracting.

ii). It is parabolic for $\Bigg\lvert \alpha +\dfrac{153}{8} \Bigg\rvert = 72$ and repulsive in other cases.

Proof. The stability function of $ex_0(\alpha)$ is

Let us remark that $\Bigg\lvert R'_{\alpha}(1)\Bigg\rvert\neq 0$ for all $\alpha\in \mathbb{C}$, so it cannot be superattracting. Therefore,

Let us consider $\alpha = a+ib$ is an arbitrary complex number, then we have the inequation

By simplifying, we get

Thus,

□

It is possible to analyze the stability of a strange fixed point numerically by graphing its stability surface, which is a three-dimensional representation of the stability function. On the $z$-axis, the stability function of the fixed point $\lvert R_{\alpha} '(ex_i(\alpha)) \rvert$ is represented, while on the complex plane we have the real and imaginary parts of parameter $\alpha$. These diagrams allow us to identify the regions of the plane where a strange fixed point is an attractor or a repulsor. Figures 1 and 2 display the stability surfaces of $ex_i(\alpha)$ for $i = 0, 1, 2, 3, 4, 5, 6$ and show that the strange fixed points are repulsive in the grey complex area, while the colored surface represents the attraction zone, that is, the zone of complex values of $\alpha$ where the strange fixed points are attracting. In the boundary between those areas, the strange fixed points are parabolic or neutral.

By using this technique, it can be stated that

$({\rm{i}})$ $ex_{1}\left(\alpha \right)$ is parabolic for $\alpha = \dfrac{423}{8}$, where $ex_{1}\left(\dfrac{423}{8}\right) = ex_{0}\left(\dfrac{423}{8}\right) = 1$; in other cases, $ex_{1}\left(\alpha \right)$ for $\alpha \neq \dfrac{423}{8}$ is a repulsor (see Figure 2(a)).

$({\rm{ii}})$ $ex_{2}\left(\alpha \right)$ is parabolic for $\alpha = \dfrac{423}{8}$, where $ex_{1}\left(\dfrac{423}{8}\right) = ex_{0}\left(\dfrac{423}{8}\right) = 1$; it is also parabolic for $\alpha = \frac{5}{288}(20\pm \sqrt{5}i)$. For any other value of $\alpha$, it is a repulsor (see Figure 2(b)).

$({\rm{iii}})$ $ex_{3}\left(\alpha \right)$ and $ex_{4}\left(\alpha \right)$ have the same stability: They are superattracting for $\alpha \approx 0.346163 \pm 0.0326984 i$, attracting in the colored area of Figure 2(c) and (d). For any other value of $\alpha$, they are repulsive.

$({\rm{iv}})$ $ex_{5}\left(\alpha \right)$ is superattracting for $\alpha = 30.2862$, attracting for any $\alpha$ in the region of the complex plane delimited by $\left\lbrace z \in \mathbb{C}: z = a+ib, 20 < a < 55, -18 < b < 18\right\rbrace$ (see Figure 2(e)). It is parabolic in its boundary and, for any other value of $\alpha$, is repulsive.

$({\rm{v}})$ $ex_{6}\left(\alpha \right)$ is superattracting for $\alpha = 30.2862$. For any other value of $\alpha$, it is a repulsor (see Figure 2(f)).

Now, we know some areas of the complex plane where choosing the value of $\alpha$ must be avoided, as the rational function $R_\alpha(x)$ has an attracting strange fixed point, and our iterative methods can fail. However, attracting fixed points are not the only pathological elements that can be found. In order to get them, we use the critical points of the rational function.

3.2.   Critical points: calculation and analysis

The critical points are those making null the derivative of the rational operator. Note that superattracting fixed points are simultaneously critical points. Those critical points that do not coincide with the zeros of the polynomial $p(x)$ ($0$ and $\infty$, after Möbius map) are called free critical points, which are related to the different basins of attraction, according to Julia-Fatou Theorem 1.

Theorem 5. Critical points of the rational operator $R_\alpha(x)$ are $x = 0$ and $x = \infty$, directly related to the zeros of the polynomial and free critical points:

i). $cr_{1}(\alpha) = -1;$

ii). $cr_{i}(\alpha) = \frac{1}{3} \left(-1\pm2 i \sqrt{2}\right), i = 2, 3;$

iii). $cr_{i}\left(\alpha\right) = \frac{1}{24}\left(1-\frac{2565}{A}-\frac{8B}{\sqrt{A^2}}\mp 16\sqrt{\frac{153 \left(-5265+D\right)-4\alpha\left(E+D\right)}{A^3}}\right), i = 4, 5;$

iv). $cr_{i}\left(\alpha\right) = \frac{1}{24}\left(1-\frac{2565}{A}+\frac{8B}{\sqrt{A^2}}\mp 16\sqrt{\frac{-153\left(5265+D\right)-4\alpha\left(E-D\right)}{A^3}}\right), i = 6, 7;$

where $A$$=$$117$$+$$64 \alpha$, $B$$=$$\sqrt{2025+16 \alpha (-981+1348 \alpha)}$, $D = {\sqrt{A^2}} \cdot B$, $E$$=$$1246347+8 \alpha$$\left(91611+7616 \alpha \right)$.

Proof. Using the derivative of the rational operator,

being $r(x) = 351+1224 x+1746 x^2+1224 x^3+351 x^4$+$192 \alpha -32 x \alpha -64 x^2 \alpha -32 x^3 \alpha +192 x^4 \alpha$. The zeros of the previous equation are the critical points of the rational operator. However, the amount of critical points can be reduced for values of $\alpha$ satisfying that any root of $r(x)$, $cr_i(\alpha)$, $i = 4, 5, 6, 7$ coincides with $cr_j(\alpha)$, $j = 1, 2, 3$. □

Moreover, the critical points $cr_1(\alpha) = -1$, $cr_2(\alpha) = \frac{1}{3} \left(-1+2 i \sqrt{2}\right)$, and $cr_3(\alpha) = \frac{1}{3} \left(-1-2 i \sqrt{2}\right)$ are pre-images of $ex_{0}(\alpha) = 1$. Thus, their stability corresponds to the stability of the fixed point $ex_0(\alpha) = 1$.

Remark 1. Conjugate free critical points are those that satisfy $cr_{i}(\alpha) = \frac{1}{cr_{j}(\alpha)}\,$ for $\, i\neq\ j$. From the proof of Theorem 5, it is straightforward that $cr_{3}(\alpha)$, $cr_{5}(\alpha)$, and $cr_{7}(\alpha)$ are conjugate critical points of $cr_{2}(\alpha)$, $cr_{4}(\alpha)$, and $cr_{6}(\alpha)$, respectively.

So, there are four independent free critical points. Furthermore, for $\alpha = 9 (109 - 24 \sqrt{6})/2696\,$ and $\, \alpha = 9 (109 + 24 \sqrt{6})/2696$, critical points $cr_{4}(\alpha)$ and $cr_{6}(\alpha)$ match in pairs for $\alpha$$\approx$$- 0.795073 - 0.606513 \; i$ and $\alpha\approx- 0.657558 - 0.753404 \; i$, respectively, and critical points $cr_{5}(\alpha)$ and $cr_{7}(\alpha)$ match in pairs for $\alpha\approx- 0.795073 + 0.606513 \; i$ and $\alpha\approx- 0.657558 + 0.753404 \; i$, respectively.

Remark 2. Moreover, some strange fixed points coincide with critical point $ex_0(\alpha) = 1$ or any of its pre-images at several values of the parameter: $\alpha = -\frac{153}{8}$, $\alpha = 0$, and $\alpha = \frac{9}{8}$. At these values of $\alpha$, $x = 1$ is repulsive, so they do not generate their own basin of attraction.

3.3.   Bifurcation diagram

The bifurcation diagram is a graphical tool that represents the real values of the fixed and critical points of the rational operator depending on the parameter $\alpha$ ^[28]. The parameter is represented on the abscissa axis and the real value of the fixed or critical point is represented on the ordinate axis. It is a very useful tool for understanding the dynamical behavior of a dynamical system, confirming the stability analysis of the fixed points.

Figure 3 shows that $ex_{1}(\alpha)$, $ex_{2}(\alpha)$, $cr_{4}(\alpha)$, and $cr_{5}(\alpha)$ coincide with $cr_{1}(\alpha)$ for $\alpha = 0$, $cr_{6}(\alpha)$ and $cr_{7}(\alpha)$ coincide with each other for $\alpha = -\frac{153}{8}$. Thus, for $\alpha = 0$, $ex_{1}(\alpha)$ and $ex_{2}(\alpha)$ are superattracting fixed points. This confirms the results presented in Section 3.1.

3.4.   Parameter spaces

The main goal for drawing dynamical and parameter planes is to have a visual interpretation of the behavior of the iterative methods in terms of the initial estimation taken, or the selected member of the class, respectively. A parameter space is constructed by means of a mesh where each node is related to a different complex value of parameter $\alpha$. Once the value of $\alpha$ is selected, a free independent critical point is used as initial estimation. The point of the mesh is colored in red if the critical point converges to $0$ or $\infty$ (belongs to their basin of attraction) and it is represented in black color in other cases. In the parameter plane, the red values of the parameter are related with a good performance (convergence to the roots) and black values correspond to unstable performance.

The long-term performance of the $R_{\alpha}(x)$ operator depends on parameter $\alpha$. Therefore, it is important to find regions of the parameter plane as stable as possible and those corresponding values of $\alpha$ will give us the best members of the family (in terms of numerical stability) ^[12]. In order to obtain these parameter spaces, we must associate each point of the parameter plane with a complex value of $\alpha$, that is, with an element of (3.2). Each $\alpha$ belonging to the same connected component of the parameter plane produces sets of schemes (3.2) with similar asymptotic performance.

Parameter spaces belonging to the family of methods are generated on a 16 MB RAM computer with an Intel Core i7 processor using a MATLAB R2021b programming package.

The graph has been drawn for values of $\alpha$ in $[-100, 60] \times [-80, 80]$, with a mesh of $1000 \times 1000$ points, a maximum of $200$ iterations, and a $10^{-3}$ tolerance. It shows the behavior of a method of the LCT class associated with $\alpha$ by means of a free critical point $cr(\alpha)$ appearing in Theorem 5 as the initial estimate. In the parameter spaces displayed in this section, the values of $\alpha$ for which the method converges to zero or infinity are plotted in red, and those values that converge to any other attracting element are plotted in black.

$R_\alpha(x)$ has at most $7$ free critical points. Of them, $cr_{1}(\alpha)$, $cr_{2}(\alpha)$, and $cr_{3}(\alpha)$ have the same parameter plane corresponding to the stability of $ex_{1}$, as they are its pre-images, so their parameter planes are not necessary. The remaining free critical points are conjugated in pairs: $cr_{4}(\alpha)$ to $cr_{5}(\alpha)$, $cr_{6}(\alpha)$ to $cr_{7}(\alpha)$ (see Remark 1), which results in $2$ different parameter spaces, namely $P_{1}$ (for $x = cr_{4}(\alpha), cr_{5}(\alpha)$), and $P_{2}$ (for $x = cr_{6}(\alpha), cr_{7}(\alpha)$), shown in Figure 4.

Taking into account the wideness of the red areas in both parameter planes, the members of the proposed class can be considered as mainly stable. However, those elements in the black regions show any kind of unstable behavior. In Figure 5, the unified parameter plane ^[13] is presented, where the white color represents those values of the parameter that are simultaneously red in all parameter planes (including those from $cr_{2}(\alpha)$ and $cr_{3}(\alpha)$) meanwhile, in the black color appears those that are also black in any of the parameter planes.

Let us remark that the outer black circle corresponds to the area where the strange fixed point $ex_0(\alpha) = 1$ is attracting; meanwhile, the big circle at the right side is related to the values of $\alpha$ where $ex_6(\alpha)$ is attracting. The rest of the black areas correspond to values of the parameter that define iterative methods where there exist attracting periodic orbits of different periods.

3.5.   Dynamical planes

The stability of specific members of the family of iterative methods, given by a fixed value of $\alpha$, is studied using dynamical planes ^[12]. This graphical tool is constructed by defining a mesh where each complex node corresponds to a value of $x_{0}$. The convergence of the method to any of the attracting fixed points starting from $x_{0}$ with a maximum of $200$ iterations and a tolerance of $10^{-3}$ is shown in it. The fixed points are presented with a white circle ($\circ$), the critical points with a white square ($\Box$), and the attracting points with an asterisk ($\ast$).

To study the stability of some elements of LCT class, several values of parameter $\alpha$ are chosen inside and outside the white stability regions observed in the unified parameter space, as shown in Figure 5. We choose $\alpha = 0$, $\alpha = -\frac{1}{2}$, $\alpha = \frac{9}{8}$ and $\alpha = \frac{-153}{8}$ as values within the stability zone.

For $\alpha = -\frac{153}{8}$, $x = 1$ is not a fixed point (Theorem 3). Figure 6(a) shows 2 basins of attraction: the zero basin (orange) and the $\infty$ basin (blue). All free critical points are located in the infinity basin and the strange fixed points lay at the Julia set, as they are repulsive.

Let us notice in Figures 6(b) and 6(c) that, for parameters $\alpha = -\frac{1}{2}$ and $\alpha = 0$, respectively, there are also 2 basins of attraction of zero (orange) and $\infty$ (blue) that are superattracting. All free critical points are located in the zero basin and the behavior is stable. $cr_{1}(\alpha) = -1$ is pre-image of $x = 1$, that is, it is repulsive and belongs to the Julia set. For $\alpha = \frac{9}{8}$, Figure 6(d) shows that there are 2 basins of attraction: the zero basin (orange) and the $\infty$ basin (blue). All critical points are located in the $\infty$ basin, which is a sign of stable behavior.

As values out of stability zone, $\alpha = -55-55i$, $\alpha = \frac{423}{8}$, and $\alpha = \pm 100$ are chosen.

Figure 7(a) shows that for parameters $\alpha = -55-55i$, there are 3 basins of attraction: the zero basin (orange), the $\infty$ basin (blue), and another basin (black) that corresponds to a periodic orbit of period 3. This is an important finding because, according to Sharkovsky's theorem (^[19,20]), if there are orbits of period 3, we can assert that orbits of any period can be obtained.

For $\alpha = \mp100$ (Figure 7(b) and 7(c), respectively), there are 3 basins of attraction: the zero basin (orange), the $\infty$ basin (blue), and another one that is not related to the zeros (green). The superattracting fixed points are $0$ and $\infty$. There is an attracting fixed strange point: $ex_{0} = 1$. All free critical points are located in the green region, so they do not converge to any of the zeros of the polynomial, as they converge to $ex_{0}(\alpha) = 1$. That is a sign of unstable behavior.

For parameter $\alpha = \frac{423}{8}$, Figure 7(d) shows 3 basins of attraction: the zero $\infty$ basins and another one in black. There is an attracting strange fixed point: $ex_{0}(\alpha) = 1$. All free critical points are located in the black region, so they do not converge to any of the zeros of the polynomial, as they converge to $ex_{0}(\alpha) = 1$. It is an unstable behavior.

In the next section, we choose some members of the proposed family of methods in order to perform numerical tests, and also compare them to other known iterative schemes.

4.   Numerical results

Now, numerical tests are performed in order to corroborate the convergence and stability of the proposed family (3.1) of methods whose stability has been previously studied. The process consists of selecting values of parameter $\alpha$ that lie inside and outside the stability region of Scetion 3, and then applying them to six nonlinear test functions, with the following expressions and zeros

We conduct the study in two stages:

$(1)$ The selection of three values of $\alpha$ within the stability zone ($\alpha = -\frac{1}{2}, \; \alpha = 0, \; \alpha = \frac{9}{8}$) and three values outside the stability area ($\alpha = -100, \; \alpha = \frac{423}{8}, \; \alpha = 100$).

$(2)$ A comparative study between one selected member with best stability properties of the family of methods and three other iterative methods with four order of convergence: Chun ^[14], Jarratt ^[23], and Ostrowski ^[29].

To perform numerical tests, we choose different values for the initial estimates relative to the zero $\xi :$ very near $\left(x_{0}\approx 1.1 \; \xi \right)$, near $\left(x_{0}\approx 3 \; \xi \right)$, far $\left(x_{0}\approx 10 \; \xi \right)$, and very far $\left(x_{0}\approx 100 \; \xi \right)$.

Calculations are performed on a 16 MB of RAM computer with an Intel Core i7 processor using MATLAB R2021b programming software, taking variable precision arithmetics with 2000 digits of mantissa and an error tolerance of $\epsilon = 10^{-100}$. For each method, the stopping criterium is $\lvert x_{k+1}-x_{k}\rvert + \lvert f(x_{k+1})\rvert < \epsilon$, where $\lvert x_{k+1}-x_{k}\rvert$ is an error estimate between two consecutive iterations and $\lvert f(x_{k+1})\rvert$ being the residual error. The number of iterations (iter) required to converge to the solution is also shown.

The ACOC is calculated in order to check the theoretical order of convergence $p$. In the numerical tests performed in this section, if ACOC does not stabilize along the successive iterations, the result will display "-", and if any scheme does not converge in 50 iterations, the result will display "nc".

The execution time (ex-time) needed to converge to the solution is shown (in seconds) and it is calculated as the average of 10 consecutive runs of the scheme.

4.1.   First stage

Now, we consider values of $\alpha$ within the stability areas of the parameter planes $\left(\alpha = -\frac{1}{2}, 0, \frac{9}{8} \right)$ and also outside $\left(\alpha = -100, \frac{423}{8}, 100 \right)$.

In Tables 1–3, the numerical performance of iterative methods associated with values of $\alpha$ inside the stability region for very near, near, far, and also very far from initial guesses are shown. From these tables we notice that the schemes associated with $\alpha = -\frac{1}{2}, 0, \frac{9}{8}$ always reach the solution (except in $f_{4}(x), f_{5}(x)$, and $f_{6}(x)$ for far and very far values from $\xi$). The amount of iterates needed to converge varies for each seed and each test function, but when the initial estimation is very close to the zero, the methods converge to $\xi$ in 5 iterations. When the initial estimation is close to the zero, they converge to $\xi$ between a minimum of 5 and a maximum of 11 iterations. In cases where the estimates are far and very far from the zero, the iterative schemes reach $\xi$ between 7 and 26 iterations. For $f_{1}(x)$, when the initial estimate is very far from the zero, it takes more than 1600 iterations for the iterative process to reach $\xi$.

Tables 4–6 show that methods related to $\alpha = -100, \frac{423}{8}, 100$ do not always reach the solution (except, in a few cases, for initial estimates very near $\xi$), confirming that these values are indeed in the unstable area, even when the function is not a polynomial. The convergence depends, to a great extent, on the nonlinear test function used and the initial estimate. Thus, for estimates near, far, and very far from the zero, these methods do not reach the solution.

4.2.   Second stage

In this stage, we perform a comparative analysis of an LCT($\alpha$) class considering three methods with the same order of convergence: Chun ^[14], Jarratt ^[23], and Ostrowski ^[29] in order to contrast their performance on nonlinear equations. A stable member of LCT($\alpha$) family is chosen, the scheme related to $\alpha = -\frac{1}{2}$, that is, LCT$(-\frac{1}{2})$, based on its performance at the first stage.

Tables 7–10 show the numerical results of LCT $\left(-\frac{1}{2} \right)$ and three known methods, considering near, very near, far, and very far seeds. Therefore, from the results shown in these tables, we conclude that the LCT $\left(-\frac{1}{2} \right)$ scheme has similar or best numerical performance to comparison methods, considering a stable element of the family $\left(\alpha = -\frac{1}{2} \right)$.

For initial estimations close to the zero, the Chun's scheme does not converge to the zero in $f_{5}(x)$ and the Ostrowski's scheme does not converge to the zero in $f_{6}(x)$.

For initial estimations far from the zero, none of the schemes converge in $f_{5}(x)$. For estimations very far away from the zero, the schemes do not converge in $f_{5}(x)$ and $f_{6}(x)$.

5.   Conclusions

A new family of optimal fourth-order multipoint methods has been proposed for solving nonlinear equations $f(x) = 0$. As the order of convergence is not the only key fact to be considered, the stability of the class has been analyzed by using complex dynamical tools. This study allows us to select the elements of the family with wider sets of converging initial estimations, as well as the member with chaotic behavior. These results are confirmed by the numerical tests, in which we also compare the proposed methods with other known ones, with the same order of convergence.

Future works based on these results are focused on the extension to vectorial problems and the semilocal analysis of convergence, using majorizing sequences.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

The authors would like to thank the anonymous reviewers for their comments and suggestions that have improved the final version of this manuscript.

Conflict of interest

Alicia Cordero and Juan R. Torregrosa are the Guest Editor of special issue “Recent advances in iterative procedures for solving nonlinear problems” for AIMS Mathematics. Alicia Cordero and Juan R. Torregrosa were not involved in the editorial review and the decision to publish this article.