Fixed Point Root-Finding Methods of Fourth-Order of Convergence

Abstract

In this manuscript, by using the weight-function technique, a new class of iterative methods for solving nonlinear problems is constructed, which includes many known schemes that can be obtained by choosing different weight functions. This weight function, depending on two different evaluations of the derivative, is the unique difference between the two steps of each method, which is unusual. As it is proven that all the members of the class are optimal methods in the sense of Kung-Traub’s conjecture, the dynamical analysis is a good tool to determine the best elements of the family in terms of stability. Therefore, the dynamical behavior of this class on quadratic polynomials is studied in this work. We analyze the stability of the presented family from the multipliers of the fixed points and critical points, along with their associated parameter planes. In addition, this study enables us to select the members of the class with good stability properties.

1. Introduction

There is a large number of problems in Science and Engineering modeled by nonlinear equations $f\left(z\right)=0$, where $f:D\subseteq {\mathbb{R}}^n\to {\mathbb{R}}^n$ is a vectorial real function defined in a convex set D, with $n\ge 1$. Usually these problems cannot be analytically solved, so we use iterative schemes. The best known is Newton’s method, the characteristics of which have been improved in much of the research (see, for instance, the texts [1,2] and the references therein). Different techniques have been employed to design these new methods being, in the majority of cases, multistep schemes, due to the limitations of one-step methods.

However, in many applications in geometric modeling, it is necessary to find all the roots of a nonlinear problem within a domain. In these cases, subdivision solvers are the most usual approach when looking for roots globally (see, for example [3] and the references therein). As subdivision solvers usually employ Newton’s scheme when the initial estimation is close to the root, higher order methods as the proposed in this manuscript may be useful to improve the local behavior of these solvers.

The weight function procedure is used to achieve this goal, which is able to increase the order of convergence of known methods (see [1]), achieving optimal schemes under the point of view of Kung-Traub’s conjecture [4]. This conjecture claims that an iterative method without memory which uses d functional evaluations per iteration can reach, at most, order of convergence $2^{d-1}$, being optimal in this bound. Although the aim of many researches in this area is to design optimal high-order methods, it is also known that the higher the order is, the more sensitive the scheme to initial estimations will be.

In our study, through weight function procedure and Newton’s scheme, a new family of iterative methods is presented, with the following iterative expression:

$$\begin{array}{ccc}{y}_n\hfill & =\hfill & z_n-\gamma {\displaystyle \frac{f\left(z_n\right)}{f^{\prime }\left(z_n\right)}},\hfill \\ z_{n+1}\hfill & =\hfill & z_n-H\left(t_n\right){\displaystyle \frac{f\left(z_n\right)}{f^{\prime }\left(z_n\right)}},n=0,1,\dots ,\hfill \end{array}$$

(1)

where $gamma$ is a real parameter and the variable of the weight function $H\left(t\right)$ is $t=\frac{f^{\prime }\left(y\right)}{f^{\prime }\left(z\right)}$.

In the next section, we present the convergence result of this family, describing the conditions that function $H\left(t\right)$ must satisfy for reaching order four. We also analyze how to extend this structure (1) for designing derivative-free methods and for solving nonlinear systems, $n>1$. Many known methods, for equations or systems, are particular cases of our class. In Section 3, we study the dynamical behavior of a sub-class obtained by using a particular weight function. The aim of this section is to check how, using complex dynamics tools, we can determine which elements of a given family have good stability properties and which have chaotic behavior.

2. Convergence and Stability

The order of convergence of family (1) is analyzed in the next result, for $n=1$.

Theorem 1.

Let $f:D\subseteq \mathbb{R}\to \mathbb{R}$ be a sufficiently differentiable function in an open interval D and let $\overline{z}\in D$ be a simple solution of the nonlinear equation $f\left(z\right)=0$. Starting from a known initial estimation $z_0$ close enough to $\overline{z}$, if $\gamma =2/3$ and function H satisfies $H\left(1\right)=1$, $H^{\prime }\left(1\right)=-3/4$, $H^{″}\left(1\right)=9/4$ and $|H^{‴}\left(1\right)|<+\infty$, then sequence ${\left\{z_n\right\}}_{n\ge 0}$ obtained from (1) converges to $\overline{z}$ with order of convergence four, being the error equation

$$e_{n+1}=\left(\left(5+\frac{38}{81}{H}^{‴}\left(1\right)\right)c_2^3-c_2{c}_3+c_4/9\right)e_n^4+O\left(e_n^5\right),$$

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

Proof. 

It is known that $f\left(z_n\right)$ and $f^{\prime }\left(z_n\right)$ can be expressed as

$$f\left(z_n\right)=f^{\prime }\left(\overline{z}\right)[e_n+c_2{e}_n^2+c_3{e}_n^3+c_4{e}_n^4+O\left(e_n^5\right)]$$

and

$$f^{\prime }\left(z_n\right)=f^{\prime }\left(\overline{z}\right)[1+2c_2{e}_n^2+3c_3{e}_n^3+4c_4{e}_n^4+O\left(e_n^5\right)].$$

By direct division,

$$\frac{f\left(z_n\right)}{f^{\prime }\left(z_n\right)}=e_n-c_2{e}_n^2+(2c_2^2-2c_3)e_n^3+[-8c_2^3+c_2(4c_2^2-3c_3)+10c_2{c}_3-3c_4]e_n^4+O\left(e_n^5\right).$$

Then,

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

Again, by expanding in Taylor series,

$$\begin{array}{cc}\hfill f^{\prime }\left(y_n\right)& =f^{\prime }\left(\overline{z}\right)\left[1-2\left[(-1+\gamma )c_2\right]e_n+[2\gamma c_2^2+3{(-1+\gamma )}^2{c}_3]e_n^2\right.\\ & \left.+2[-2\gamma c_2^3+(5-3\gamma )\gamma c_2{c}_3-2{(-1+\gamma )}^3{c}_4]e_n^3\right]+O\left(e_n^4\right).\end{array}$$

Therefore,

$$\begin{array}{ccc}\hfill t_n={\displaystyle \frac{f^{\prime }\left(y_n\right)}{f^{\prime }\left(z_n\right)}}& =& 1-2\gamma c_2{e}_n+3\gamma [2c_2^2+(-2+\gamma )c_3]e_n^2\hfill \\ & & -4\left[\gamma (4c_2^3+(-7+3\gamma )c_2{c}_3+(3-3\gamma +{\gamma }^2)c_4)\right]e_n^3+O\left(e_n^4\right)\hfill \end{array}$$

and, taking into account that t tends to 1 when n tends to infinity, the weight function can be expressed as

$$\begin{array}{ccc}\hfill H\left(t_n\right)& =& h_0+h_1(t_n-1)+\frac{1}{2}{h}_2{(t_n-1)}^2+\frac{1}{6}{h}_3{(t_n-1)}^3+O\left({(t_n-1)}^4\right)\hfill \\ & =& h_0-2\gamma h_1{c}_2{e}_n+[2{\gamma }^2{h}_2{c}_2^2+3\gamma h_1(2c_2^2+(-2+\gamma )c_3)]]e_n^2\hfill \\ & & +\left[-\frac{4}{3}{\gamma }^3{h}_3{c}_2^3-6{\gamma }^2{h}_2{c}_2(2c_2^2+(-2+\gamma )c_3)-4\gamma h_1(4c_2^3+(-7+3\gamma )c_2{c}_3\right.\hfill \\ & & \left.+(3-3\gamma +{\gamma }^2)c_4)\right]e_n^3+O\left(e_n^4\right),\hfill \end{array}$$

where $h_0$, $h_1$, $h_2$ and $h_3$ denotes $H\left(1\right)$, $H^{\prime }\left(1\right)$, $H^{″}\left(1\right)$ and $H^{‴}\left(1\right)$, respectively. Finally,

$$e_{n+1}=z_{n+1}-\overline{z}=z_n-\overline{z}-H\left(t_n\right)\frac{f\left(z_n\right)}{f^{\prime }\left(z_n\right)}=$$

$$=(1-h_0)e_n+(h_0+2\gamma h_1)c_2{e}_n^2+[-2(h_0+\gamma (4h_1+\gamma h_2))c_2^2+(2h_0-3(-2+\gamma )\gamma h_1)c_3]e_n^3$$

$$+[(4h_0+26\gamma h_1+14{\gamma }^2{h}_2+\frac{4}{3}{\gamma }^3{h}_3)c_2^3+(-7h_0+\gamma ((-38+15\gamma )h_1+6(-2+\gamma )\gamma h_2))c_2{c}_3$$

$$+(3h_0+4\gamma (3-3\gamma +{\gamma }^2)h_1)c_4]e_n^4+O\left(e_n^5\right).$$

Therefore, in order to eliminate first order error,

$$1-h_0=0,\phantom{aa}\mathrm{which} \mathrm{means} \mathrm{that}\phantom{aa}{h}_0=H\left(1\right)=1.$$

For the second order error:

$$h_0+2\gamma h_1=0,$$

being $h_0=1$

$$h_1=H^{\prime }\left(1\right)=\frac{-1}{2\gamma }.$$

With these values of $h_0$ and $h_1$, the error equation is now

$$\begin{array}{ccc}\hfill e_{n+1}& =& [(2-2{\gamma }^2{h}_2)c_2^2+\frac{1}{2}(-2+3\gamma )c_3]e_n^3\hfill \\ & & +[(-9+14{\gamma }^2{h}_2+\frac{4{\gamma }^3{h}_3}{3})c_2^3\frac{3}{2}(8-5\gamma -8{\gamma }^2{h}_2+4{\gamma }^3{h}_2)c_2{c}_3\hfill \\ & & +(-3+6\gamma -2{\gamma }^2)c_4]e_n^4+O(e_n^5).\hfill \end{array}$$

So, we can eliminate the third order error if

$$2-2{\gamma }^2{h}_2=0\phantom{aa}\mathrm{and}\phantom{aa}-2+3\gamma =0,$$

which means that

$$h_2=H^{″}\left(1\right)=\frac{1}{{\gamma }^2}\phantom{aa}\mathrm{and}\phantom{aa}\gamma =\frac{2}{3}.$$

Hence it is proven that if $\gamma =2/3$ and function H satisfies $H\left(1\right)=1$, $H^{\prime }\left(1\right)=-3/4$ and $H^{″}\left(1\right)=9/4$, the order of convergence is four. □

If we want to design a derivative-free scheme with a similar structure as (1), we change $f^{\prime }\left(z_n\right)$ by the divided difference $f[z_n,w_n]$, where $w_n=z_n+\rho f\left(z_n\right)$ with $\rho$ a real parameter, obtaining the following iterative expression:

$$\begin{array}{ccc}{y}_n\hfill & =\hfill & z_n-\gamma {\displaystyle \frac{f\left(z_n\right)}{f[z_n,w_n]}},\hfill \\ z_{n+1}\hfill & =\hfill & z_n-H\left(t_n\right){\displaystyle \frac{f\left(z_n\right)}{f[z_n,w_n]}},n=0,1,\dots ,\hfill \end{array}$$

(2)

being in this case $t={\displaystyle \frac{f[y,w]}{f[z,w]}}$ the variable of the weight function. Unfortunately, family (2) does not reach order four. However, if we use as variable of the weight function $t={\displaystyle \frac{f[y,z]}{f[z,w]}}$, being $w_n=z_n+\rho f{\left(z_n\right)}^2$, the class

$$\begin{array}{ccc}{y}_n\hfill & =\hfill & z_n-\gamma {\displaystyle \frac{f\left(z_n\right)}{f[z_n,w_n]}},\hfill \\ z_{n+1}\hfill & =\hfill & z_n-H\left(t_n\right){\displaystyle \frac{f\left(z_n\right)}{f[z_n,w_n]}},n=0,1,\dots ,\hfill \end{array}$$

(3)

reaches order four under certain conditions as we establish in the following result

Theorem 2.

Let $f:D\subseteq \mathbb{R}\to \mathbb{R}$ be a sufficiently differentiable function in an open interval D and let $\overline{z}\in D$ be a simple solution of the nonlinear equation $f\left(z\right)=0$. Starting from a known initial estimation $z_0$ close enough to $\overline{z}$, if $\gamma =1$ and function H satisfies $H\left(1\right)=1$, $H^{\prime }\left(1\right)=-1$, $H^{″}\left(1\right)=4$ and $|H^{‴}\left(1\right)|<+\infty$, then sequence ${\left\{z_n\right\}}_{n\ge 0}$ obtained from (3) converges to $\overline{z}$ with order of convergence four, being the error equation

$$e_{n+1}=(-f^{\prime }{\left(\overline{z}\right)}^2\rho c_2^2+\frac{1}{6}(30+h_3)c_2^3-c_2{c}_3)e_n^4+O\left(e_n^5\right),$$

where $c_j=\frac{f^{\left(j\right)}\left(\overline{z}\right)}{j!f^{\prime }\left(\overline{z}\right)}$, $j=2,3,\dots$, $e_n=z_n-\overline{z}$ and $h_3=H^{‴}\left(1\right)$.

On the other hand, family (1) can be extended to the multidimensional case, $n>1$, so we accomplish a family of iterative methods for solving nonlinear systems $F\left(z\right)=0$, where $F:D\subseteq {\mathbb{R}}^n\to {\mathbb{R}}^n$. In this case, the iterative expression is

$$\begin{array}{ccc}{y}^{\left(n\right)}\hfill & =\hfill & z^{\left(n\right)}-\gamma {\left[F^{\prime }\left(z^{\left(n\right)}\right)\right]}^{-1}F\left(z^{\left(n\right)}\right),\hfill \\ z^{(n+1)}\hfill & =\hfill & z^{\left(n\right)}-H\left(t^{\left(n\right)}\right){\left[F^{\prime }\left(z^{\left(n\right)}\right)\right]}^{-1}F\left(z^{\left(n\right)}\right),n=0,1,\dots ,\hfill \end{array}$$

(4)

where $F^{\prime }\left(z^{\left(n\right)}\right)$ is the Jacobian matrix of F evaluated in the iteration $z^{\left(n\right)}$. The variable of weight is $t={\left[F^{\prime }\left(z\right)\right]}^{-1}{F}^{\prime }\left(y\right)$ and $H\left(t\right)$ is a matrix function $H:X\to X$, where $X={\mathbb{R}}^{n\times n}$, such that

(i)

$H^{\prime }\left(u\right)\left(v\right)=h_{1}uv$, being $H^{\prime }$ the first derivative of H, $H^{\prime }:X\to \mathcal{L}\left(X\right)$, $h_1\in \mathbb{R}$ and $\mathcal{L}\left(X\right)$ denotes the space of linear mappings from X to itself.

(ii)

$H^{″}(u,v)\left(w\right)=h_{2}uvw$, being $H^{″}$ the second derivative of H, $H^{″}:X\times X\to \mathcal{L}\left(X\right)$ and $h_2\in \mathbb{R}$.

(iii)

$H^{‴}(u,vw)\left(s\right)=h_{3}uvws$, being $H^{‴}$ the third derivative of H, $H^{‴}:X\times X\times X\to \mathcal{L}\left(X\right)$ and $h_3\in \mathbb{R}$.

Then, the Taylor expansion of H around the identity matrix I gives

$$H(t^{(n)})\approx H(I)+h_1(t^{(n)}-I)+\frac{1}{2}{h}_2{(t^{(n)}-I)}^2+\frac{1}{6}{h}_3{(t^{(n)}-I)}^3.$$

A similar result to Theorem 1 can be established for class (4), obtaining conditions for function $H\left(t\right)$ to reach order four. In the proof we use some tools and notations introduced in [5].

Theorem 3.

Let $F:D\subseteq {\mathbb{R}}^n\to {\mathbb{R}}^n$ be a sufficiently differentiable function in an convex set D and let $\overline{z}\in D$ be a solution of the nonlinear system $F\left(z\right)=0$, such that $F^{\prime }\left(\overline{z}\right)$ is nonsingular. Starting from a known initial estimation $z^{\left(0\right)}$ close enough to $\overline{z}$, if $\gamma =2/3$ and function H satisfies $H\left(I\right)=I$, $h_1=-3/4$, $h_2=9/4$ and $|h_3|<+\infty$, then sequence ${\left\{z_n\right\}}_{n\ge 0}$ obtained from (4) converges to $\overline{z}$ with order of convergence four, being the error equation

$$e_{n+1}=\left(\left(5+\frac{38}{81}{h}_3\right)C_2^3-C_2{C}_3+C_4/9\right)e_n^4+O\left(e_n^5\right),$$

where $C_j=\frac{1}{j!}{\left[F^{\prime }\left(\overline{z}\right)\right]}^{-1}{F}^{\left(j\right)}\left(\overline{z}\right)$, $j=2,3,\dots$ and $e_n=z^{\left(n\right)}-\overline{z}$.

Proof. 

By using Taylor expansion of $F\left(z^{\left(n\right)}\right)$ and $F^{\prime }\left(z^{\left(n\right)}\right)$ around $\overline{z}$,

$$F\left(z^{\left(n\right)}\right)=F^{\prime }\left(\overline{z}\right)[e_n+C_2{e}_n^2+C_3{e}_n^3+C_4{e}_n^4]+O\left(e_n^5\right),$$

$$F^{\prime }\left(z^{\left(n\right)}\right)=F^{\prime }\left(\overline{z}\right)[I+2C_2{e}_n+3C_3{e}_n^2+4C_4{e}_n^3]+O\left(e_n^4\right).$$

From the above expression, we conjecture

$$F^{\prime }{\left(z^{\left(n\right)}\right)}^{-1}=[I+X_2{e}_n+X_3{e}_n^2+X_4{e}_n^3]F^{\prime }{\left(\overline{z}\right)}^{-1}+O\left(e_n^4\right)$$

and, from $F^{\prime }{\left(z^{\left(n\right)}\right)}^{-1}{F}^{\prime }\left(z^{\left(n\right)}\right)=F^{\prime }\left(z^{\left(n\right)}\right)F^{\prime }{\left(z^{\left(n\right)}\right)}^{-1}=I$, we have

$$F^{\prime }{\left(z^{\left(n\right)}\right)}^{-1}=[I-2C_2{e}_n+(4C_2^2-3C_3)e_n^2+(-4C_4+6C_2{C}_3+6C_3{C}_2-8C_2^3)e_n^3]F^{\prime }{\left(\overline{z}\right)}^{-1}+O\left(e_n^4\right).$$

Then,

$$F^{\prime }{\left(z^{\left(n\right)}\right)}^{-1}F\left(z^{\left(n\right)}\right)=e_n-C_2{e}_n^2+2(C_2^2-C_3)e_n^3+(-3C_4+4C_2{C}_3+3C_3{C}_2-4C_2^3)e_n^4+O\left(e_n^5\right).$$

Similarly, we calculate

$$y^{\left(n\right)}-\overline{z}=(1-\gamma )e_n+\gamma C_2{e}_n^2-2\gamma (C_2^2-C_3)e_n^3-\gamma (-3C_4+4C_2{C}_3+3C_3{C}_2-4C_2^3)e_n^4+O\left(e_n^5\right).$$

So,

$$\begin{array}{ccc}\hfill F^{\prime }\left(y^{\left(n\right)}\right)& =& F^{\prime }\left(\overline{z}\right)\left[I+2C_2(y^{\left(n\right)}-\overline{z})+3C_3{(y^{\left(n\right)}-\overline{z})}^2+4C_4{(y^{\left(n\right)}-\overline{z})}^3\right]+O{(y^{\left(n\right)}-\overline{z})}^4)\hfill \\ & =& F^{\prime }\left(\overline{z}\right)\left[I+2C_2(1-\gamma )e_n+(2\gamma C_2^2+3C_3{(1-\gamma )}^2)e_n^2+(-4\gamma (C_2^2-C_2{C}_3)+6\gamma (1-\gamma )C_3{C}_2\right.\hfill \\ & & \left.+4{(1-\gamma )}^3{C}_4)e_n^3\right]+O\left(e_n^4\right).\hfill \end{array}$$

Therefore, variable t of the weight function is described as

$$\begin{array}{ccc}\hfill t=F^{\prime }{\left(z_n\right)}^{-1}{F}^{\prime }\left(y_n\right)& =& I-2\gamma C_2{e}_n+(6\gamma C_2^2+3C_3({\gamma }^2-2\gamma ))e_n^2\hfill \\ & & +(-16\gamma C_2^3+(6{\gamma }^2-16\gamma )C_2{C}_3+(12\gamma -6{\gamma }^2)C_3{C}_2\hfill \\ & & +(-4{\gamma }^3+12{\gamma }^2-12\gamma )C_4)e_n^3+O\left(e_n^4\right)\hfill \end{array}$$

and

$$\begin{array}{ccc}\hfill H\left(t\right)& =& H\left(I\right)+h_1(t-I)+\frac{1}{2}{h}_2{(t-I)}^2+\frac{1}{6}{h}_3{(t-I)}^3\hfill \\ & =& H\left(I\right)-2\gamma C_2{h}_1{e}_n+(6\gamma h_1{C}_2^2+3h_1{C}_3({\gamma }^2-2\gamma )+2h_2{\gamma }^2{C}_2^2)e_n^2\hfill \\ & & +\left[-4/3\gamma (12h_1+\gamma (9h_2+h_3\gamma ))C_2^3-4h_1\gamma ({\gamma }^2-3\gamma +3)C_4\right.\hfill \\ & & \left.\gamma (h_1(16-6\gamma )-3h_2\gamma (\gamma -2))C_2{C}_3-3\gamma (\gamma -2)(2h_1+h_2\gamma )C_3{C}_2\right]e_n^3+O\left(e_n^4\right).\hfill \end{array}$$

Therefore,

$$\begin{array}{ccc}\hfill e_{n+1}& =& (I-H\left(I\right))e_n+(H\left(I\right)C_2+2\gamma h_1{C}_2)e_n^2+\left(-2H\left(I\right)(C_2^2-C_3)-2\gamma h_1{C}_2^2\right.\hfill \\ & & \left.-8\gamma h_1{C}_2^2-3h_1{C}_3({\gamma }^2-2\gamma )-2h_2{\gamma }^2{C}_2^2\right)e_n^3+O\left(e_n^4\right).\hfill \end{array}$$

In order to reach order four, the coefficients of $e_n$, $e_n^2$ and $e_n^3$ must be zero, so

$$H\left(I\right)=I,\phantom{aaa}{h}_1=\frac{-1}{2\gamma },\phantom{aaa}{h}_2=\frac{9}{4}\phantom{aaa}and\phantom{aaa}\gamma =\frac{2}{3}.$$

With these values, the error equation is

$$e_{n+1}=\left(\left(5+\frac{38}{81}{h}_3\right)C_2^3-C_2{C}_3+C_4/9\right)e_n^4+O\left(e_n^5\right)$$

and the proof is finished. □

Many known schemes designed for solving nonlinear equations or nonlinear systems can be obtained as particular cases of (4) by using different weight functions satisfying the conditions of Theorem 3. For example, the classical Jarratt’s scheme [6]

$$\begin{array}{ccc}{y}_n\hfill & =\hfill & z_n-\frac{2}{3}{\displaystyle \frac{f\left(z_n\right)}{f^{\prime }\left(z_n\right)}},\hfill \\ z_{n+1}\hfill & =\hfill & z_n-\frac{3f^{\prime }\left(y_n\right)+f^{\prime }\left(z_n\right)}{6f^{\prime }\left(y_n\right)-2f^{\prime }\left(z_n\right)}{\displaystyle \frac{f\left(z_n\right)}{f^{\prime }\left(z_n\right)}},n=0,1,\dots ,\hfill \end{array}$$

(5)

is obtained from (1) for equations and (4) for systems, by using

$$H\left(t\right)={(6t-2I)}^{-1}(3t+I).$$

In a similar way, the method constructed by Khattri and Abbasbandy in [7] is an element of (4), $n\ge 1$, by using the weight function

$$H\left(t\right)=I+\frac{21}{8}t-\frac{9}{2}{t}^2+\frac{15}{8}{t}^3.$$

The parametric family of iterative methods for solving nonlinear equations or systems, designed by Kanwar et al. in [8], is a particular case of (1) or (4) using the weight function

$$H\left(t\right)=\frac{1}{2}\frac{{\alpha }_1^2-22{\alpha }_1{\alpha }_2-27{\alpha }_2^2+3({\alpha }_1^2+10{\alpha }_1{\alpha }_2+5{\alpha }_2^2)t}{({\alpha }_1+3{\alpha }_{2}t)(3({\alpha }_1+{\alpha }_2)t-{\alpha }_1-5{\alpha }_2)},$$

where ${\alpha }_1$ and ${\alpha }_2$ are free real parameters such that ${\alpha }_1\ne {\alpha }_2$ and ${\alpha }_1\ne 3{\alpha }_2$.

Sharma and Arora published in [9] a method for solving nonlinear systems which we can obtain from (4) by using

$$H\left(t\right)=\frac{23}{8}I+t\left(-3I+\frac{9}{8}t\right).$$

Hueso et al. presented in [10] a parametric family of iterative schemes, which is obtained from the weight function

$$H\left(t\right)=\frac{5-8\alpha }{8}I+\alpha t^{-1}+\frac{\alpha }{3}t+\frac{9-8\alpha }{24}{t}^{-2},$$

where $\alpha$ is a real free parameter.

Finally, Ghorbanzadeh and Soleymani presented in [11] an iterative method solving nonlinear equations or nonlinear systems, which is a particular case of our scheme using the weight function

$$H\left(t\right)=4{(I+3t)}^{-1}\left[I+\frac{9}{16}{(t-I)}^2\right].$$

One of the most simple sub-classes of family (1) is obtained when $H\left(t\right)$ is the Taylor polynomial of third degree, satisfying $\gamma =2/3$, $H\left(1\right)=1$, $H^{\prime }\left(1\right)=-3/4$ and $H^{″}\left(1\right)=9/4$, that is, the weight function is:

$$H\left(t\right)=1-\frac{3}{4}(t-1)+\frac{9}{8}{(t-1)}^2+\frac{1}{6}\alpha {(t-1)}^3,$$

being $\alpha =H^{‴}\left(1\right)$ the free parameter. In this way, we obtain a parametric family of iterative methods of order four, denoted by CGT$\alpha$. In the next sections, we are going to analyze de dynamical behavior of this class in terms of parameter $\alpha$, for finding the methods with good stability properties and to avoid the elements of the family with chaotic behavior.

3. Stability of Family Cgt$\mathit{\alpha }$

We are going to introduce some concepts of complex dynamics that we need to analyze the stability of the proposed family on quadratic polynomials. These are the most simple nonlinear functions and, although the stability conclusions for quadratic polynomials can not be extended to any nonlinear function, in practice a bad method for quadratic polynomials does not perform adequately for other functions. Moreover, cubic polynomials and high-order iterative methods usually yields to rational operators of very high degree; this makes the analysis unaffordable in the most of cases.

After that, we will calculate the fixed and critical points of the rational operator associated with the class on polynomials. The stability of the fixed points, the behavior of the critical points as initial estimation of the iterative schemes and the basins of attraction of the periodic points give us relevant information about the stability of the methods.

3.1. Basic Concepts On Dynamics

A more complete description of the following concepts can be found in [12]). Given a rational function $R:\widehat{\mathbb{C}}\to \widehat{\mathbb{C}}$, where $\widehat{\mathbb{C}}$ is the Riemann sphere, the orbit of a point $z_0\in \widehat{\mathbb{C}}$ is defined as:

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

We analyze the phase plane of operator R by classifying the starting points from the asymptotic behavior of their orbits.

A $z_0\in \widehat{\mathbb{C}}$ is called a fixed point if $R\left(z_0\right)=z_0$. A periodic point $z_0$ of period $p>1$ is a point such that $R^p\left(z_0\right)=z_0$ and $R^k\left(z_0\right)\ne z_0$, for $k<p$. A pre-periodic point is a point $z_0$ that is not periodic but there exists a $k>0$ such that $R^k\left(z_0\right)$ is periodic. A critical point $z_0$ is a point where the derivative of the rational function vanishes, $R^{\prime }\left(z_0\right)=0$. Moreover, a fixed point $z_0$ is called attractor if $|R^{\prime }\left(z_0\right)|<1$, superattractor if $|R^{\prime }\left(z_0\right)|=0$, repulsor if $|R^{\prime }\left(z_0\right)|>1$ and parabolic if $|R^{\prime }\left(z_0\right)|=1$.

The basin of attraction of an attractor $\overline{z}$ is defined as:

$$\mathcal{A}\left(\overline{z}\right)=\{z_0\in \widehat{\mathbb{C}}:R^n\left(z_0\right)\to \overline{z},n\to \infty \}.$$

The immediate basin of attraction of an attractor is the connected component of its basin of attraction that holds the attractor.

The Fatou set of the rational function R, $\mathcal{F}\left(R\right)$, is the set of points $z\in \widehat{\mathbb{C}}$ whose orbits tend to an attractor (fixed point, periodic orbit or infinity). Its complement in $\widehat{\mathbb{C}}$ is the Julia set, $\mathcal{J}\left(R\right)$. That means that the basin of attraction of any fixed point belongs to the Fatou set and the boundaries of these basins of attraction belong to the Julia set.

The following Theorem establishes a classical result of Fatou and Julia that we use in the study of parameter space associated with the family.

Theorem 4

([13]). Let R be a rational function. The immediate basin of attraction of an attracting fixed or periodic point holds, at least, a critical point.

By using this result, one can be sure to find all the stable behaviors associated with a rational function R, by analyzing the performance of R on the set of critical points.

Taking into account that the proposed family of iterative methods satisfies the Scaling Theorem, it can be seen in the literature ([12]) that the roots of a polynomial can be changed by an affine map without qualitative changes on the dynamics of the family. Therefore, we can use a generic quadratic polynomial $p\left(z\right)=(z-a)(z-b)$.

By replacing $f\left(x\right)$ by $p\left(z\right)$ in the iterative expression of the proposed family CGT$\alpha$, the following fixed point operator is obtained:

$$\begin{array}{c}\hfill T_{p,\alpha ,a,b}\left(z\right)=z+\frac{(-a+z)(-b+z)(1+\frac{(a-z)(b-z)}{{(a+b-2z)}^2}+\frac{2{(a-z)}^2{(b-z)}^2}{{(a+b-2z)}^4}+\frac{32\alpha {(a-z)}^3{(-b+z)}^3}{81{(a+b-2z)}^6})}{a+b-2z},\end{array}$$

which depends on parameters $\alpha$, a and b.

In [13], Blanchard considered the conjugacy map ${\displaystyle \varphi \left(z\right)=\frac{z-a}{z-b}}$, a Möbius transformation with these properties:

$$\left(\mathrm{i}\right)\varphi \left(\infty \right)=1,\left(\mathrm{ii}\right)\varphi \left(a\right)=0,\left(\mathrm{iii}\right)\varphi \left(b\right)=\infty .$$

It was proven that Newton’s operator, for quadratic polynomials, is conjugate to the rational map $z^2$. In an similar way, operator $T_{p,\alpha ,a,b}$ on quadratic polynomials is conjugated to operator $O_{\alpha }\left(z\right)$,

$$O_{\alpha }\left(z\right)=\left(\varphi \circ T_{p,\alpha ,a,b}\circ {\varphi }^{-1}\right)\left(z\right)=z^4\frac{405+32\alpha +1134z+1134z^2+486z^3+81z^4}{81+486z+1134z^2+1134z^3+405z^4+32\alpha z^4}.$$

(6)

In this new operator $O_{\alpha }$, parameters a and b have been obviated.

3.2. Study of the Fixed Points

Stability and reliability of the members of the family are analyzed in the rest of the paper, regarding the properties of its associate rational function when the class is applied on polynomial $p\left(z\right)$. Firstly, fixed points of the rational function $O_{\alpha }\left(z\right)$ are calculated. Specifically, we focus on the points that are not related to the original roots of polynomial $p\left(z\right)$, which are named strange fixed points.

In the next result, some properties of the strange fixed points are described.

Theorem 5.

Fixed points of the rational function $O_{\alpha }\left(x\right)$ are the roots of equation $O_{\alpha }\left(z\right)=z$. Therefore, we obtain $z=0$, $z=\infty$ (corresponding to the roots of $p\left(z\right)$), and the following strange fixed points:

• $e{x}_1\left(\alpha \right)=1$,
• the roots of symmetric sixth-degree polynomial $r\left(t\right)=81+567t+1701t^2+(2430-32\alpha )t^3+1701t^4+567t^5+81t^6$, or analogously,

  $$\begin{array}{ccc}\hfill e{x}_{2,3}\left(\alpha \right)& =& \frac{s_1\left(\alpha \right)\pm \sqrt{s_1{\left(\alpha \right)}^2-4}}{2},e{x}_2\left(\alpha \right)=\frac{1}{e{x}_3\left(\alpha \right)},\hfill \\ \hfill e{x}_{4,5}\left(\alpha \right)& =& \frac{s_2\left(\alpha \right)\pm \sqrt{s_2{\left(\alpha \right)}^2-4}}{2},e{x}_4\left(\alpha \right)=\frac{1}{e{x}_5\left(\alpha \right)},\hfill \\ \hfill e{x}_{6,7}\left(\alpha \right)& =& \frac{s_3\left(\alpha \right)\pm \sqrt{s_3{\left(\alpha \right)}^2-4}}{2},e{x}_6\left(\alpha \right)=\frac{1}{e{x}_7\left(\alpha \right)},\hfill \end{array}$$

where $s_1\left(\alpha \right)$, $s_2\left(\alpha \right)$ and $s_3\left(\alpha \right)$ are the roots of the third-degree polynomial $s\left(t\right)=81t^3+567t^2+1458t+1296-32\alpha$, that is,

$$\begin{array}{ccc}\hfill s_1\left(\alpha \right)& =& \frac{1}{9}\left(-21-\frac{153^{1/3}}{{\left(\mathsf{\Phi }\left(\alpha \right)\right)}^{1/3}}+{\left(9\mathsf{\Phi }\left(\alpha \right)\right)}^{1/3}\right),\hfill \\ \hfill s_2\left(\alpha \right)& =& \frac{1}{18}\left(-42+\frac{153^{1/3}\left(1+i\sqrt{3}\right)}{{\left(\mathsf{\Phi }\left(\alpha \right)\right)}^{1/3}}+i\left(i+\sqrt{3}\right){\left(9\mathsf{\Phi }\left(\alpha \right)\right)}^{1/3}\right),\hfill \\ \hfill s_3\left(\alpha \right)& =& \frac{1}{18}\left(-42+\frac{153^{1/3}\left(1-i\sqrt{3}\right)}{{\left(\mathsf{\Phi }\left(\alpha \right)\right)}^{1/3}}-\left(1+i\sqrt{3}\right){\left(9\mathsf{\Phi }\left(\alpha \right)\right)}^{1/3}\right),\hfill \end{array}$$

being $\mathsf{\Phi }\left(\alpha \right)=24+16\alpha +\sqrt{1701+768\alpha +256{\alpha }^2}$.

It is possible to find values of parameter α where two or more strange fixed points coincide. Consequently, operator $O_{\alpha }\left(x\right)$ provides seven strange fixed points, except in the following cases:

(i) 

If $\alpha =0$ or $\alpha =\frac{891}{4}$, there are only five strange fixed points.

(ii) 

If $\alpha =\frac{3}{16}i\left(8i+5\sqrt{5}\right)$ or $\alpha =-\frac{3}{16}i\left(-8i+5\sqrt{5}\right)$ the strange fixed points are $e{x}_1=1$ and the roots of a sixth-degree polynomial, which are two simple and two double roots.

3.3. Stability of the Fixed Points

In this section, we observe that not only the number but also the stability of the fixed points depend on the parameter of the family. This fact can lead to the existence of attracting strange fixed points, which can make the iterative scheme converge to a false solution.

It is known that $z=0$ and $z=\infty$ are always superattracting fixed points, as the order of convergence of the class is greater than 2. However, relevant numerical information is provided by the stability of the other fixed points (for example, $z=1$ corresponds to the divergence of the original method). Therefore, we determine this stability in the next results.

Theorem 6.

Strange fixed point $e{x}_1\left(\alpha \right)=1$, with $\alpha \ne -{\displaystyle \frac{405}{4}}$, has the following character:

(i) 

When $\left|\alpha +\frac{405}{4}\right|<324$ , $e{x}_1\left(\alpha \right)=1$ is a repulsor.

(ii) 

If $\left|\alpha +\frac{405}{4}\right|=324$, $e{x}_1\left(\alpha \right)=1$ is a parabolic point.

(iii) 

When $\left|\alpha +\frac{405}{4}\right|>324$, then $e{x}_1\left(\alpha \right)=1$ is an attractor.

Proof. 

It is easy to prove that

$$O_{\alpha }^{\prime }\left(1\right)=\frac{1296}{405+4\alpha }.$$

So,

$$\left|\frac{1296}{405+4\alpha }\right|\le 1\phantom{aaa}\mathrm{is} \mathrm{equivalent} \mathrm{to}\phantom{aaa}1296\le \left|405+4\alpha \right|.$$

Let us consider $\alpha =c+id$ an arbitrary complex number. Then,

$${1296}^2\le {405}^2+16c^2+16d^2+3240c.$$

By simplifying

$$1515591-3240c-16c^2-16d^2\le 0,$$

that is,

$${\left(c+\frac{405}{4}\right)}^2+d^2\le {324}^2.$$

Therefore,

$$\left|O_{\alpha }^{\prime }\left(1\right)\right|\le 1\phantom{aaa}\mathrm{if} \mathrm{and} \mathrm{only} \mathrm{if}\phantom{aaa}\left|\alpha +\frac{405}{4}\right|\le 324.$$

□

Stability regions of every strange fixed point ($e{x}_i\left(\alpha \right)$, $i=1,2,\dots ,7$) are shown in Figure 1.

We can observe from the stability of $e{x}_i\left(\alpha \right),i=4,5,6,7$ that these strange fixed points are repulsive for any value of parameter $\alpha$.

3.4. Analysis of the Critical Points

Regarding the dynamics of the family, it is significant to analyze the critical points of the rational function $O_{\alpha }\left(z\right)$ different from 0 and ∞, which are called free critical points.

With the purpose of calculating the critical points, we study which points make null the derivative of $O_{\alpha }\left(z\right)$.

$$O_{\alpha }^{\prime }\left(z\right)=324z^3{(1+z)}^6\frac{405{(1+z)}^2+16\alpha \left(2-3z+2z^2\right)}{{\left(81+486z+1134z^2+1134z^3+(405+32\alpha )z^4\right)}^2}.$$

We know that $z=0$ and $z=\infty$, which are related to the roots of the polynomial by means of Möbius map, are critical points and give rise to their respective Fatou components. Nevertheless, several free critical points can be obtained, some of them depending on the value of the parameter $\alpha$.

Proposition 1.

The number of critical points of the rational function $O_{\alpha }\left(z\right)$ depend on the value of the parameter α:

(a) 

If $\alpha =0$, there exist one only free critical point, $z=-1$, which is a pre-image of the fixed point $z=1$.

(b) 

In all other cases, the free critical points are $c{r}_1\left(\alpha \right)=-1$,

$$\begin{array}{ccc}\hfill c{r}_2\left(\alpha \right)& =& \frac{-405+24\alpha -4\sqrt{7}\sqrt{-\alpha (405+4\alpha )}}{405+32\alpha }=\frac{1}{c{r}_3},\hfill \end{array}$$

which means that there exists only one independent free critical point.

Proof. 

We obtain the previous critical points, since

$$O_{\alpha }^{\prime }\left(z\right)=324z^3{(1+z)}^6\frac{405{(1+z)}^2+16\alpha \left(2-3z+2z^2\right)}{{\left(81+486z+1134z^2+1134z^3+(405+32\alpha )z^4\right)}^2}=$$

$$=324(405+32\alpha )z^3{(1+z)}^6\frac{(z-c{r}_2\left(\alpha \right))(z-c{r}_3\left(\alpha \right))}{{(81+486z+1134z^2+1134z^3+(405+32\alpha )z^4)}^2}.$$

□

In the case of strange fixed point $e{x}_i\left(\alpha \right),i=2,3$, we can see that there exist a ball in which the strange fixed points are attractors, whereas outside this ball, $e{x}_i\left(\alpha \right),i=2,3$ are repulsive.

3.5. The Parameter Plane

We have observed that the dynamical behavior of operator $O_{\alpha }\left(x\right)$ depends on the values of the parameter $\alpha$. Taking into account Theorem 4, we are interested in knowing what happens with the free critical points, and if any of them gives rise to a basin of attraction different from those of zero and infinity. In order to have knowledge of this, we obtain the parameter plane associated with the family.

We depict the parameter space associated with a free critical point of operator (6) by associating each point of the parameter plane with a complex value of $\alpha$, i.e., with an element of family (1). Every value of $\alpha$ belonging to the same connected component of the parameter space gives rise to subsets of schemes of family (1) with similar dynamical behavior. Therefore, our objective is to find regions of the parameter plane as much stable as possible, due to the fact that the values of $\alpha$ in these regions will provide the best members of the family in terms of numerical stability.

Since $c{r}_2\left(\alpha \right)={\displaystyle \frac{1}{c{r}_3\left(\alpha \right)}}$, we have at most one free independent critical point, consequently, there exist only one parameter plane of the family. If we consider the independent free critical point as a starting point of the iterative schemes of the family associated with each complex value of $\alpha$, this point of the complex plane is painted in red if the method converges to any of the roots (zero and infinity) and they are black in other cases. The lower the number of iterations, the brighter the color used. Following this procedure, we obtain the parameter plane presented in Figure 2, by using the processes described in [14]. We have used a mesh of $2000\times 2000$ points, 500 maximum iterations and ${10}^{-3}$ as the tolerance used in the stopping criterium. We also show a zoom of this parameter plane in Figure 3a, focusing on the biggest red area, where the members of family (1) are, in general, very stable. These would be the best values, a priori, to apply the iterative methods in terms of stability.

Nevertheless, we can see black regions that inform us about different pathological behavior of the elements of the family. The black ball represented in Figure 3b, correspond to values of parameter $\alpha$ for which $e{x}_i\left(\alpha \right),i=2,3$ are attracting. Besides, the big black ball that surrounds the parameter plane in Figure 2 correspond to values of the parameter for which $e{x}_1\left(\alpha \right)$ is attracting.

In addition, we can analyze the remaining black regions by drawing dynamical planes with the parameter $\alpha$ corresponding to values inside these black regions.

4. Dynamical Planes

Now, we will analyze the qualitative behavior of the different elements of family through the dynamical planes. These elements are selected taking into account the conclusions obtained by studying the parameter plane of the family. Dynamical planes with stable behavior are depicted by using points in the red region of the parameter plane, whereas dynamical planes with unstable performance are calculated with points in the black region of the parameter plane of the family.

Every dynamical plane presented here has been generated by using the routines appearing in [14]. The dynamical plane related to a value of parameter $\alpha$ is obtained by iterating an element of family (1). Each point of the complex plane is used as initial estimation. The points whose orbit converges to infinity are plotted in blue, in orange those points converging to zero, in green the points whose orbit converges to one of the strange fixed points, which appear marked as a white star in the figures if they are attracting fixed points, and in black if it reaches the maximum number of 50 iterations without converging to any of the fixed points. A mesh of $600\times 600$ points has been used and a tolerance of ${10}^{-3}$.

Some dynamical planes are shown in Figure 4. These planes correspond to values of the parameter $\alpha$ which from the parameter plane, give us elements of the family with stable behavior. Therefore, we can see only two basins of attraction, that correspond to zero (orange basin) and infinity (blue basin).

On the other hand, as it has been said, unstable behavior is found when we choose values of the parameter in the black region of parameter plane. These dynamical planes are shown in Figure 5.

In Figure 5a,c,e,f the dynamical planes of the iterative method corresponding to $\alpha =65$, $\alpha =85+35i$, $\alpha =-140$ and $\alpha =-200$, respectively, are presented, with regions of slow convergence. Figure 5b, corresponding to $\alpha =90$, shows the existence of four different basins of attraction, two of them of the superattractors 0 and ∞ and the other ones corresponding to strange fixed points.

Finally, in Figure 5d, corresponding to $\alpha =\frac{891}{4}$, we can observe an orbit of period two.

5. Conclusions

In this paper, a class of optimal iterative methods for solving nonlinear equations is presented, which holds many known methods as particular elements. This family is extended in two different directions: a class of derivative-free methods for solving nonlinear equations and a multidimensional family for nonlinear systems. Both classes preserve the order of convergence of the initial family. The dynamical analysis in these areas will be object of study in future works.

Moreover, a complex dynamical study for a parametric sub-family on quadratic polynomials has been presented. We have been able to prove, by means of the parameter plane, that there are many values of $\alpha$ with good stability properties, which means that there exist plenty of elements of the family with stable behavior. Nevertheless, there are other values of $\alpha$ with convergence anomalies that must be avoided in practical applications.