# King-Type Derivative-Free Iterative Families: Real and Memory Dynamics

Authors: F. I. Chicharro; A. Cordero; J. R. Torregrosa; M. P. Vassileva
Journal: Complexity
 (2017)

---

## Abstract

A biparametric family of derivative-free optimal iterative methods of order four, for solving nonlinear equations, is presented. From the error equation of this class, different families of iterative schemes with memory can be designed increasing the order of convergence up to six. The real stability analysis of the biparametric family without memory is made on quadratic polynomials, finding areas in the parametric plane with good performance. Moreover, in order to study the real behavior of the parametric class with memory, we associate it with a discrete multidimensional dynamical system. By analyzing the fixed and critical points of its vectorial rational function, we can select those methods with best stability properties.


---

## Body

## 1. Introduction

The solution of a nonlinear equationf(x)=0 is a common problem in topics related to several disciplines of science. However, most of these equations can not be solved analytically. In this case, iterative methods acquire prominence, because they can reach an approximate solution of f(x). Newton’s scheme is the most-applied iterative method, but the evaluation of the derivative of f is its most remarkable inconvenience.Several authors have made great efforts to study the complex dynamics of some derivative-free iterative methods [1–3], but the analysis of the real dynamics is not as usual as the complex one [4, 5].During the last years, many researchers have focused on the construction of optimal multipoint methods without memory based on the Kung and Traub’s conjecture [6], which claims that to create an optimal method without memory with optimal order 2n, n+1 functional evaluations are required. Construction and development of multipoint methods with memory have more advantages in comparison with methods without memory. In other words, they are able to increase the convergence order without any new functional evaluation. Due to this basic fact, they can also get higher efficiency index (the efficiency index was defined by Ostrowski [7] as I=p1/d, where p is the order of convergence and d is the number of functional evaluations required per iteration). As another plus, they generate more stable computations than optimal methods without memory, even with higher convergence orders.The basic idea for the construction of multipoint methods with memory was introduced by Traub [8], who presented a version with memory from Steffensen’s method. Recently, based on this method, some schemes with memory have been developed by several authors. We can see interesting overviews in [9, 10].In this paper, we deeply analyze the real dynamics of a version with memory of a King-type derivative-free iterative family. Two new tools are introduced in the real dynamics section: the unified parameter plane and the dynamical line. The former gathers the information of different parameters planes into one plane. The latter becomes a comfortable tool to visualize the dynamical behavior of a method for a set of initial points in the real line.Starting from King’s family of fourth-order schemes(1)yk=xk-fxkf′xk,xk+1=yk-fxk+αfykfxk+α-2fykfykf′xk,where α is a disposable parameter, we design a derivative-free family of fourth-order methods. The order of convergence can be increased introducing memory, as can be deduced from the error equation.In the manuscript, we use symbols~ and O in the following way: if limn→∞an/bn=C, with C being a nonzero constant, we write an=O(bn) or an~bn. The notation and techniques used in the proofs of the results can be found in [9–11].The rest of this paper is organized as follows. In Section2, we introduce the derivative-free parametric family and its error equation. Section 3 addresses the increase of the order of convergence with the introduction of the memory in the iterative schemes. In Section 4, the real dynamics of the biparametric family is studied for different quadratic polynomials, while Section 5 covers the dynamical study of the methods with memory. Finally, Section 6 collects some conclusions about the obtained results.


## 2. Parametric Families of Iterative Schemes

Following the structure of King’s family and replacing the derivative by a first-order divided difference we present the following schemes:(2)yk=xk-fxkfxk,wk,xk+1=yk-fxk+αfykμfxk+βfykfykfyk,wk,where wk=xk+γf(xk), with α, β, γ, and μ being real parameters, γ≠0.The order of convergence of methods (2) is established in the following result. The proof only requires the development in Taylor series of the elements of the iterative expression and some algebraic manipulations.Theorem 1.
Let us suppose thatf:I⊆R→R is a sufficiently differentiable function in an open interval I and x∗∈I is a simple root of f(x)=0. If the initial approximation x0 is close enough to x∗, then the iterative scheme (2) has optimal fourth-order convergence when μ=1 and β=α-1 for all nonzero γ, being in this case the error equation (3)ek+1=1+γf′x∗2c22+α+αγf′x∗c22-c3ek4+Oek5,where ek=xk-x∗, k=0,1,…, and cj=(1/j!)f(k)(x∗)/f′(x∗),k≥2.Proof.
By expanding in Taylor seriesf(xk), we get(4)fxk=f′x∗ek+c2ek2+c3ek3+c4ek4+c5ek5+c6ek6+Oek7,fxk,wk=f′x∗1+2+γf′x∗c2ek+γf′x∗c22+3+3γf′x∗+γ2f′x∗2c3ek2+2+γf′x∗2γf′x∗c2c3+2+2γf′x∗+γ2f′x∗2c4ek3+γ2f′x∗2c22c3+γf′x∗3+2γf′x∗c32+γf′x∗7+8γf′x∗+3γ2f′x∗2c2c4+5+10γf′x∗+10γ2f′x∗2+5γ3f′x∗3+γ4f′x∗4c5ek4+γ2f′x∗24+3γf′x∗c22c4+γf′x∗9+10γf′x∗+3γ2f′x∗2c3c4+γf′x∗c22γf′x∗c32+11+20γf′x∗+15γ2f′x∗2+4γ3f′x∗3c5+6+15γf′x∗+20γ2f′x∗2+15γ3f′x∗3+6γ4f′x∗4+γ5f′x∗5c6ek5+Oek6.Therefore, the error at the first step is (5)yk-x∗=1+f′x∗γc2ek2+-2+2f′x∗γ+f′x∗2γ2c22+2+3f′x∗γ+f′x∗2γ2c3ek3+4+5f′x∗γ+3f′x∗2γ2+f′x∗3γ3c23-7+10f′x∗γ+7f′x∗2γ2+2f′x∗3γ3c2c3+3+6f′x∗γ+4f′x∗2γ2+f′x∗3γ3c4ek4+Oek5.By obtaining the Taylor expansion of f(yk) and f[yk,wk], we get the following expression for the error equation: (6)ek+1=1+f′x∗γ-1+μc2μek2+1μ2β1+f′x∗γ2-μ-3+α1+f′x∗γ2+2f′x∗γ-2+μ+f′x∗2γ2-2+μ+2μc22+2+3f′x∗γ+f′x∗2γ2-1+μμc3ek3+Oek4. It is clear that, imposing μ=1, third-order is achieved and, then, the error equation is (7)ek+1=--1+α-β1+f′x∗γ2c22ek3+α1+f′x∗γ6+7f′x∗γ+3f′x∗2γ2+β1+f′x∗γ2-1+f′x∗γ3+3f′x∗γ+2f′x∗2γ2+β+f′x∗βγ2+β6+7f′x∗γ+3f′x∗2γ2c23-1+f′x∗γ2-3-2f′x∗γ+2α2+f′x∗γ-2β2+f′x∗γc2c3ek4+Oek5. Again, setting β=α-1, we get order four and(8)ek+1=1+f′x∗γ2c22+α+f′x∗αγc22-c3ek4-1+f′x∗γ10+11f′x∗γ+5f′x∗2γ2+α21+f′x∗γ3+α7+13f′x∗γ+8f′x∗2γ2+2f′x∗3γ3c24-14+19f′x∗γ+7f′x∗2γ2+3α1+f′x∗γ22+f′x∗γc22c3+2+3f′x∗γ+f′x∗2γ2c32+2+3f′x∗γ+f′x∗2γ2c2c4ek5+Oek6.Let us remark not only that the fourth-order has been proven but also that factor (1+f′(x∗)γ) appears also in the term of order five.Let us observe that this family of derivative-free methods (2) supports Kung-Traub’s conjecture [6], having optimal efficiency index I=41/3≈1.587. On the other hand, by observing the expressions of the error equation (8) we can choose different values of the free disposable parameters in order to obtain iterative methods with memory, increasing the order of convergence.


## 3. Iterative Methods with Memory

We are going to design derivative-free schemes with memory based on the proposed methods (2).From (8) we can assure that the order of convergence of family (2) increases up to six if γ=-1/f′(x∗), but the value of f′(x∗) is not available in practice and such acceleration is not possible. However, we can use an approximation f′¯(x∗)≈f′(x∗), calculated by using known information. Therefore, by setting γ=-1/f′¯(x∗) we can increase the convergence order without using new functional evaluations. The main idea in constructing methods with memory consists of the calculation of the parameter γ=γk as the iteration proceeds by γk=-1/f′¯(x∗), k=1,2,…. We are going to consider different approximations of f′(x∗).(1) Let N1(t)=N1(t,xk,xk-1) be Newton’s interpolation polynomial of first degree through two available approximations xk, xk-1; that is N1(t)=f(xk)+f[xk,xk-1](t-xk), so(9)γk=-1N1′xk=-1fxk,xk-1=-xk-xk-1fxk-fxk-1, and the algorithm denoted by MM1(α) can be presented in the following way: (i)
x0,γ0 are given.(ii)
wk=xk+γkf(xk),k=0,1,2,….(iii)
yk=xk-fxk/fxk,wk.(iv)
xk+1=yk-f(xk)+αf(yk)/f(xk)+(α-1)f(yk)f(yk)/f[yk,wk],whereγk=-1/f[xk,xk-1] and α is a free parameter.The order of convergence of this algorithm is analyzed in the following result.Theorem 2.
Letx∗ be a simple zero of a sufficiently differentiable function f:I⊆R→R in an open interval I. If x0 is close enough to x∗ and γ0 is given, then the R-order of family MM1(α) is at least 2+6≈4.45 that corresponds to the positive root of polynomial p2-4p-2.Proof.
By expanding in Taylor seriesf(xk) and f(xk-1), we get(10)fxk=f′x∗ek+c2ek2+c3ek3+c4ek4+c5ek5+c6ek6+Oek7,fxk-1=f′x∗ek-1+c2ek-12+c3ek-13+c4ek-14+c5ek-15+c6ek-16+Oek-17. Then,(11)γk=-1fxk,xk-1=1f′x∗-1+c2ek+ek-1+-2c22+c3ekek-1+-c22+c3ek2+ek-12+3c23-4c2c3+c4ek2ek-1+ekek-12+c23-2c2c3+c4ek3+ek-13+O4ek,ek-1,wk-x∗=ekek-1c25-4c3c23+3c4c22+3c32-2c5c2-2c3c4+c6ek-14-c24-3c3c22+2c4c2+c32-c5ek-13+c23-2c3c2+c4ek-12+c3-c22ek-1+c2+ek2ek-14c25-13c3c23+7c4c22+8c32-3c5c2-4c3c4+c6ek-13+-3c24+7c3c22-3c4c2-2c32+c5ek-12+2c23-3c3c2+c4ek-1-c22+c3+O4ek,ek-1,where O4(ek,ek-1) indicates that the sum of the exponents of ek-1 and ek in the rejected terms of the development is at least 4.
The error at the next step is(12)yk-x∗=c2ek2ek-1c2+c3-c22ek-1+c23-2c3c2+c4ek-12-c24-3c3c22+2c4c2+c32-c5ek-13+ek3ek-1-2c23+2c2c3+2c24-4c22c3+c32+c2c4ek-1+-2c25+6c23c3-4c22c4+c3c4+c2-2c32+c5ek-12+O5ek,ek-1and, then,(13)fxkfyk,wk=c22ek-1+-c23+c2c3ek-12+c2c23-2c2c3+c4ek-13ek2+-2c23+2c2c3ek-1+c24-4c22c3+c32+c2c4ek-12ek3+O5ek,ek-1.Moreover,(14)fxk+αfykfxk+α-1fyk=c22ek-1+-c23+c2c3ek-12+c2c23-2c2c3+c4ek-13ek+-3c23+2c2c3ek-1+--4+αc24-5c22c3+c32+c2c4ek-12ek2+6c24-9c22c3+c32+2c2c4ek-13ek3+O5ek,ek-1.Finally, by combining the previous expansions from (12)–(14), (15)ek+1=2c25-c23c3ek4ek-12+O6ek,ek-1. Since the lower term of the error equation is (2c25-c23c3)ek-12ek4, applying Theorem  9.2.9 of [12], the powers of ek and ek-1 are 4 and 2, respectively, so the polynomial whose real roots give the R-order of the method is p2-4p+2 and the order of the method is, at least, 2+6.(2) Let N2(t)=N2(t,xk,xk-1,yk-1) (can be also used as N2(t)=N2(t,xk,xk-1,wk-1)) be Newton’s interpolation polynomial of second degree; that is, N2(t)=f(xk)+f[xk,xk-1](t-xk)+f[xk,xk-1,yk-1](t-xk)(t-xk-1), and therefore(16)γk=-1N2′xk=-1fxk,xk-1+fxk,xk-1,yk-1xk-xk-1, and a similar algorithm to the previous one, denoted by MM2(α), can be presented. The convergence of the method is established in the following result, whose proof is similar to that of the last proposed method. This is the reason why it is avoided.Theorem 3.
Letx∗ be a simple zero of a sufficiently differentiable function f:I⊆R→R in an open interval I. If x0 is close enough to x∗ and γ0 is given, then the R-order of family MM2 is at least 5.(3) Let N3(t)=N3(t,xk,xk-1,yk-1,wk-1) be Newton’s interpolation polynomial of third degree: N3(t)=f(xk)+f[xk,xk-1](t-xk)+f[xk,xk-1,yk-1](t-xk)(t-xk-1)+f[xk,xk-1,yk-1,wk-1](t-xk)(t-xk-1)(t-yk-1), and then (17)N3′xk=fxk,xk-1+fxk,xk-1,yk-1xk-xk-1+fxk,xk-1,yk-1,wk-1xk-xk-1xk-yk-1,γk=-1N3′xk. The algorithm denoted by MM3(α) can be presented in the following way: (i)
x0, γ0 are given.(ii)
wk=xk+γkf(xk),k=0,1,2,….(iii)
yk=xk-f(xk)/f[xk,wk].(iv)
xk+1=yk-f(xk)+αf(yk)/f(xk)+(α-1)f(yk)f(xk)/f[yk,wk],whereγk=-1/N3′(xk) and α is a free parameter.Theorem 4.
Letx∗ be a simple zero of a sufficiently differentiable function f:I⊆R→R in an open interval I. If x0 is close enough to x∗ and γ0 is given, then the R-order of family MM3 is at least 6.Proof.
By using Taylor expansion aroundx∗ of the different elements of N3′(xk), we prove that(18)1+γkf′x∗=1-f′x∗N3′x~c4ek-1ey,k-1ew,k-1~ek-1ey,k-1ew,k-1,where ek=xk-x∗, ey,k=yk-x∗, and ew,k=wk-x∗.
Let{xk} be the sequence generated by the algorithm MM3. If {xk} converges to x∗ with R-order r, then(19)ek+1~Dk,rekr, such that limk→∞Dk,r=Dr, with this limit being the asymptotic error constant of the method. So,(20)ek+1~Dk,rDk-1,rek-1rr=Dk,rDk-1,rrek-1r2.
Let us suppose that{wk} has R-order r1, and then(21)ew,k~Dk,r1ekr1~Dk,r1Dk-1,rrek-1rr1.
Analogously, if sequence{yk} has R-order r2, we obtain(22)ey,k~Dk,r2Dk-1,rrek-1rr2.
On the other hand, we obtain the error equation of the different sequences without memory(23)ek+1=1+γf′x∗2Mek4+Oek5,where M=c2(2+α(1+γf′(x∗))c22-c3) and(24)ew,k=1+γf′x∗ek+Oek2,ey,k=c21+γf′x∗ek2+Oek3.
Therefore, the corresponding error relations with memory are(25)ek+1~1+γf′x∗2Mkek4,ew,k~1+γf′x∗ek,ey,k~c21+γf′x∗ek2, and, by using expression (18) and simple algebraic manipulations, we obtain(26)ek+1~Dk-1,r12Dk-1,r22Dk-1,r4ek-14r+2r1+2r2+2,(27)ew,k~Dk-1,r1Dk-1,r2Dk-1,rek-1r+r1+r2+1,(28)ey,k~c2Dk-1,r1Dk-1,r2Dk-1,r2ek-12r+r1+r2+1.If we compare the exponents of ek-1 between “(20), (21), (22), and (26)), ((21), (22), (26), and (27)” and “(22), (26), (27), and (28),” we obtain the following nonlinear system:(29)r2=4r+2r1+2r2+2,rr1=r+r1+r2+1,rr2=2r+r1+r2+1.
It is easy to obtain the unique real solution of this system:r1=2, r2=3, and r=6. Therefore, the R-order of algorithm MM3 is 6 for any value of parameter α.Many other approximations off′(x∗) are possible, but they either are of sixth order (with more computational cost) or are lower than six.


## 4. Real Dynamics

In this section, the real dynamics of the iterative family are analyzed. After an introduction of real dynamics fundamentals [13, 14], the features of the family are described when it is applied over the nonlinear polynomials x2-1, x2, and x2+1.


### 4.1. Fundamentals on Real Dynamics

LetM:R→R be a rational function. The orbit of a point x0∈R is defined as the set (30)x0,Mx0,M2x0,…,Mnx0,…. A point x0 is a fixed point, x0F, of M if M(x0F)=x0F. The multiplier M′x0F classifies the fixed points as attracting, repelling, or neutral if its value is lower than, greater than, or equal to 1, respectively. When the multiplier is zero, the fixed point is superattracting. Both attracting and superattracting fixed points are denoted as x∗.The character of the infinity is set by the study ofI(x)=1/M1/x [15]. The point ∞ is a fixed point of M if and only if x=0 is a fixed point of I. In this case, the multiplier of ∞ is I′0. When this value is zero, the point ∞ is a superattracting point.The basin of attraction of an attracting fixed pointx∗, A(x∗), is defined as the set of preimages of any order such that (31)Ax∗=x0∈R:Mnx0⟶x∗,n⟶∞.The Fatou set,F(M), includes the points whose orbits tend to an attracting point x∗. The Julia set, J(M), is its complementary. It covers the repelling points and sets the borders between the basins of attraction.The fixed point operator of the biparametric family (2) is(32)Mfx=y-fx+αfyfx+α-1fyfyfy,w,where w=x+γf(x) and y=x-f(x)/f[x,w].


### 4.2. Application on Quadratic Polynomials

We are going to analyze the dynamical behavior of the rational function obtained when family (32) is applied on f(x)∈{x2-1,x2,x2+1}.


#### 4.2.1.f(x)=x2-1

In this case, expression (32) is Mf(x)=N11(x)/D10(x), where N11(x) and D10(x) are polynomials of degrees 11 and 10, respectively, depending on x and on parameters α and γ.There are 10 fixed points ofMf(x). x1,2F=x1,2∗=±1 are superattracting points, while x3-10F(α,γ) are the roots of an 8th-degree polynomial: x3,4F∈R and x5-10F∈C. Since our purpose is the study of the real dynamics, x5-10F are rejected. The evaluation of Mf′x3,4F establishes the behavior of these points. For different values of α and γ, x3,4F have different dynamical features, as the stability plane of Figure 1 represents. A mesh of 100×100 points covers the values of α∈-2,2andγ∈[-2,2]. The white area represents where the multiplier is lower than 1, while the corresponding black represents where the multiplier is greater than 1. Its dynamical meaning is immediate, since white and black regions represent attracting and repelling behavior, respectively.Figure 1
Stability plane of (a)x3F and (b) x4F.
(a)
(b)Therefore,x3F has both behaviors, depending on the value of α and γ. For assuring the convergence of the method to x1,2∗, white areas must be avoided.The point∞ is a parabolic fixed point, since I(0)=0 and I′(0)=1.ComputingMf′(x)=0, 16 critical points can be found. As expected, x1,2C=x1,2∗. x3,4C=-1±γ/γ are preimages of x1,2∗.  x5-16C are the roots of a polynomial of degree 12, that depend on the parameters α and γ. Six of these roots are real (x5-10C) while the other ones are complex (x11-16C).A classical result (attributed to Fatou and Julia) establishes that there is at least one critical point associated with each invariant Fatou component, so it is interesting to analyze the behavior of each free critical point (critical point different to the roots of the polynomial) used as initial estimation for all the elements of the family.Just as the convergence plane [4] represents the final state of the orbits of every starting point in a single-parametric family of methods, we use the parameters plane to guess the final state of the orbit of one critical point in a biparametric family of methods, in a similar manner as [16]. Therefore, the parameters plane represents the family of methods Mf(x), where each point (α,γ) stands for an individual method. If the orbit of the free critical point xC tends to -1,1, or ∞, the point is colored in red, green, or blue, respectively. If xC does not converge to any of those points, the point is colored in black. Figure 2 shows the parameter planes of x5-10C.Figure 2
Parameter planes of (a)x5C, (b) x6C, (c) x7C, (d) x8C, (e) x9C, and (f) x10C (red: x1∗, green: x2∗, blue: x∞∗, and black: others).
(a)
(b)
(c)
(d)
(e)
(f)The unified parameter plane represented in Figure3 gathers in one image the main information of Figure 2. It is composed of the superposition of the black regions. In this way, the election of a point in the white region guarantees that the corresponding scheme tends to ±1 or ∞, while the set of methods in the black region do not.Figure 3
Unified parameter plane ofx5-10C.In order to visualize the behavior of each method, a new representation concept is introduced: the dynamical line. Analogous to the dynamical plane for complex variable, the dynamical line represents the basins of attraction, plotting in different colors where the orbit of each initial estimation tends. In coherence to Figure2, the red color is assigned to those initial approximations that tend to x1∗=-1. Analogously, the green color is assigned to x2∗=1, the blue color is assigned to x∗=∞, and the black one is assigned to those initial guesses that tend to a point different than the previous three. The final orbit of 1001 initial values of x∈[-2,2] has been calculated, and the convergence to a point is set when its distance is lower than 10-3.Figure4 represents the dynamical line of a set of methods whose orbits tend to -1,1, or ∞, that is, methods that belong to the white area of the unified parameter plane. On the contrary, Figure 5 shows methods whose points also tend to another fixed point.Figure 4
Dynamical lines of methods with (a)(α,γ)=(0.5,0.05), (b) (α,γ)=(0.5,1), (c) (α,γ)=(0.5,2), (d) (α,γ)=(1,0.05), (e) (α,γ)=(1,1), (f) (α,γ)=(1,2), (g) (α,γ)=(2,0.05), (h) (α,γ)=(2,1), and (i) (α,γ)=(2,2) (red: x1∗, green: x2∗, and blue: x∞∗).
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)Figure 5
Dynamical lines of methods with (a)(α,γ)=(-0.5,-0.5), (b) (α,γ)=(-0.5,0.5), (c) (α,γ)=(-1,-0.5), and (d) (α,γ)=(-1,0.5) (red: x1∗, green: x2∗, blue: x∞∗, and black: others).
(a)
(b)
(c)
(d)Note that every point in Figure4 tends to one of the attracting points. However, when a value of (α,γ) of the black region of Figure 3 is chosen, some points tend to a different point. For instance, these points can be found in Figures 5(a)–5(d) for x0=-0.2, x0=0.2, x0=-0.28, and x0=0.28, respectively.As expected from the attracting area ofx3F and the unified parameter plane, the election of methods with α>0 is capital, as proved in Figures 4 and 5.


#### 4.2.2.f(x)=x2

Applying (32) over f(x)=x2, the fixed point operator is(33)Mx=x1+γx5+5γx+γ2x2+α2+γx1+γx22+γx3+γx3+2γx+α1+γx2.The real fixed points of the operator are x1F=0 and x2F(α,γ), while the other two fixed points are complex. Both real fixed points have a different dynamical behavior depending on the parameter α. x1F is attracting for every value of α but the interval -3.25<α<-2.875. x2F is repelling for every value of α∈(-3.25,-0.098). The stability planes of Figure 6 illustrate this behavior. Consequently, for (i)
α∈(-∞,-3.25)∪(-0.098,+∞), x1F=x1∗ is attracting, and x2F is repelling,(ii)
α∈(-3.25,-2.875), x1F is repelling, and x2F=x2∗ is attracting,(iii)
α∈(-2.875,-0.098), x1F=x1∗, and x2F=x2∗ are attracting.Figure 6
Stability plane of (a)x1F and (b) x2F.
(a)
(b)From a stability point of view, just the first set of values ofα cause the family to have a good performance, since the only attracting point is the root of f(x).In a similar fashion as the casef(x)=x2-1, the infinity is a parabolic fixed point.The zero values of the derivative of (33) are the critical points, which correspond to the roots of polynomial (34)x8α2γ8+x714α2γ7+4αγ7+x677α2γ6+50αγ6+2γ6+x5224α2γ5+252αγ5+26γ5+x4383α2γ4+666αγ4+139γ4+x3398α2γ3+996αγ3+376γ3+x2247α2γ2+842αγ2+529γ2+x84α2γ+372αγ+360γ+12α2+66α+90. The polynomial has four real roots and four complex ones. No one of the critical points agrees with the roots of the function x2, so every one is free. Figure 7 shows the parameters planes of each real critical point. When the critical point tends to x1∗, x2∗, or x∞∗, the points are plotted in red, green, or blue, respectively, while black is kept for convergence to any other point. Note in Figures 7(b), 7(c), and 7(d) there is a thin black region, where the critical points have periodical orbits, so they do not converge to any one of the fixed points.Figure 7
Parameter planes of (a)x1C, (b) x2C, (c) x3C, and (d) x4C (red: x1∗, green: x2∗, blue: x∞∗, and black: others).
(a)
(b)
(c)
(d)For gathering the information of the parameters planes, the unified parameter plane is shown in Figure8.Figure 8
Unified parameter plane ofx1-4C.Finally, the dynamical behavior is analyzed with the dynamical line. The representation of the methods that converge to any of the roots is covered in Figure9, while Figure 10 illustrates the methods that also converge to any other point or have a periodical orbit. Color map is held.Figure 9
Dynamical lines of methods with (a)(α,γ)=(0.1124,-1), (b) (α,γ)=(-2,0.5), (c) (α,γ)=(1,1), and (d) (α,γ)=(-0.23,-1.5) (red: x1∗, green: x2∗, blue: x∞∗, and black: others).
(a)
(b)
(c)
(d)Figure 10
Dynamical lines of methods with (a)(α,γ)=(-2.876,1) and (b) (α,γ)=(-2.876,-1) (red: x1∗, green: x2∗, blue: x∞∗, and black: others).
(a)
(b)


#### 4.2.3.f(x)=x2+1

Finally, the application of (32) on f(x)=x2+1 results in the fixed point operator M=N11(x)/D10(x), where N11(x) and D10(x) are polynomials of degrees 11 and 10, respectively.The fixed points are the roots off(x) and the roots of an eighth-degree polynomial. Both groups of roots are complex, so they are left out of the purpose of the current study.Regarding the critical points, an analogous comment can be made. The roots off(x) are critical points. The rest of the critical points are x3,4C=-1/γ±i and the roots of a polynomial of degree 12. No critical point is real, so the real dynamics can not be analyzed.


## 4.1. Fundamentals on Real Dynamics

LetM:R→R be a rational function. The orbit of a point x0∈R is defined as the set (30)x0,Mx0,M2x0,…,Mnx0,…. A point x0 is a fixed point, x0F, of M if M(x0F)=x0F. The multiplier M′x0F classifies the fixed points as attracting, repelling, or neutral if its value is lower than, greater than, or equal to 1, respectively. When the multiplier is zero, the fixed point is superattracting. Both attracting and superattracting fixed points are denoted as x∗.The character of the infinity is set by the study ofI(x)=1/M1/x [15]. The point ∞ is a fixed point of M if and only if x=0 is a fixed point of I. In this case, the multiplier of ∞ is I′0. When this value is zero, the point ∞ is a superattracting point.The basin of attraction of an attracting fixed pointx∗, A(x∗), is defined as the set of preimages of any order such that (31)Ax∗=x0∈R:Mnx0⟶x∗,n⟶∞.The Fatou set,F(M), includes the points whose orbits tend to an attracting point x∗. The Julia set, J(M), is its complementary. It covers the repelling points and sets the borders between the basins of attraction.The fixed point operator of the biparametric family (2) is(32)Mfx=y-fx+αfyfx+α-1fyfyfy,w,where w=x+γf(x) and y=x-f(x)/f[x,w].


## 4.2. Application on Quadratic Polynomials

We are going to analyze the dynamical behavior of the rational function obtained when family (32) is applied on f(x)∈{x2-1,x2,x2+1}.


### 4.2.1.f(x)=x2-1

In this case, expression (32) is Mf(x)=N11(x)/D10(x), where N11(x) and D10(x) are polynomials of degrees 11 and 10, respectively, depending on x and on parameters α and γ.There are 10 fixed points ofMf(x). x1,2F=x1,2∗=±1 are superattracting points, while x3-10F(α,γ) are the roots of an 8th-degree polynomial: x3,4F∈R and x5-10F∈C. Since our purpose is the study of the real dynamics, x5-10F are rejected. The evaluation of Mf′x3,4F establishes the behavior of these points. For different values of α and γ, x3,4F have different dynamical features, as the stability plane of Figure 1 represents. A mesh of 100×100 points covers the values of α∈-2,2andγ∈[-2,2]. The white area represents where the multiplier is lower than 1, while the corresponding black represents where the multiplier is greater than 1. Its dynamical meaning is immediate, since white and black regions represent attracting and repelling behavior, respectively.Figure 1
Stability plane of (a)x3F and (b) x4F.
(a)
(b)Therefore,x3F has both behaviors, depending on the value of α and γ. For assuring the convergence of the method to x1,2∗, white areas must be avoided.The point∞ is a parabolic fixed point, since I(0)=0 and I′(0)=1.ComputingMf′(x)=0, 16 critical points can be found. As expected, x1,2C=x1,2∗. x3,4C=-1±γ/γ are preimages of x1,2∗.  x5-16C are the roots of a polynomial of degree 12, that depend on the parameters α and γ. Six of these roots are real (x5-10C) while the other ones are complex (x11-16C).A classical result (attributed to Fatou and Julia) establishes that there is at least one critical point associated with each invariant Fatou component, so it is interesting to analyze the behavior of each free critical point (critical point different to the roots of the polynomial) used as initial estimation for all the elements of the family.Just as the convergence plane [4] represents the final state of the orbits of every starting point in a single-parametric family of methods, we use the parameters plane to guess the final state of the orbit of one critical point in a biparametric family of methods, in a similar manner as [16]. Therefore, the parameters plane represents the family of methods Mf(x), where each point (α,γ) stands for an individual method. If the orbit of the free critical point xC tends to -1,1, or ∞, the point is colored in red, green, or blue, respectively. If xC does not converge to any of those points, the point is colored in black. Figure 2 shows the parameter planes of x5-10C.Figure 2
Parameter planes of (a)x5C, (b) x6C, (c) x7C, (d) x8C, (e) x9C, and (f) x10C (red: x1∗, green: x2∗, blue: x∞∗, and black: others).
(a)
(b)
(c)
(d)
(e)
(f)The unified parameter plane represented in Figure3 gathers in one image the main information of Figure 2. It is composed of the superposition of the black regions. In this way, the election of a point in the white region guarantees that the corresponding scheme tends to ±1 or ∞, while the set of methods in the black region do not.Figure 3
Unified parameter plane ofx5-10C.In order to visualize the behavior of each method, a new representation concept is introduced: the dynamical line. Analogous to the dynamical plane for complex variable, the dynamical line represents the basins of attraction, plotting in different colors where the orbit of each initial estimation tends. In coherence to Figure2, the red color is assigned to those initial approximations that tend to x1∗=-1. Analogously, the green color is assigned to x2∗=1, the blue color is assigned to x∗=∞, and the black one is assigned to those initial guesses that tend to a point different than the previous three. The final orbit of 1001 initial values of x∈[-2,2] has been calculated, and the convergence to a point is set when its distance is lower than 10-3.Figure4 represents the dynamical line of a set of methods whose orbits tend to -1,1, or ∞, that is, methods that belong to the white area of the unified parameter plane. On the contrary, Figure 5 shows methods whose points also tend to another fixed point.Figure 4
Dynamical lines of methods with (a)(α,γ)=(0.5,0.05), (b) (α,γ)=(0.5,1), (c) (α,γ)=(0.5,2), (d) (α,γ)=(1,0.05), (e) (α,γ)=(1,1), (f) (α,γ)=(1,2), (g) (α,γ)=(2,0.05), (h) (α,γ)=(2,1), and (i) (α,γ)=(2,2) (red: x1∗, green: x2∗, and blue: x∞∗).
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)Figure 5
Dynamical lines of methods with (a)(α,γ)=(-0.5,-0.5), (b) (α,γ)=(-0.5,0.5), (c) (α,γ)=(-1,-0.5), and (d) (α,γ)=(-1,0.5) (red: x1∗, green: x2∗, blue: x∞∗, and black: others).
(a)
(b)
(c)
(d)Note that every point in Figure4 tends to one of the attracting points. However, when a value of (α,γ) of the black region of Figure 3 is chosen, some points tend to a different point. For instance, these points can be found in Figures 5(a)–5(d) for x0=-0.2, x0=0.2, x0=-0.28, and x0=0.28, respectively.As expected from the attracting area ofx3F and the unified parameter plane, the election of methods with α>0 is capital, as proved in Figures 4 and 5.


### 4.2.2.f(x)=x2

Applying (32) over f(x)=x2, the fixed point operator is(33)Mx=x1+γx5+5γx+γ2x2+α2+γx1+γx22+γx3+γx3+2γx+α1+γx2.The real fixed points of the operator are x1F=0 and x2F(α,γ), while the other two fixed points are complex. Both real fixed points have a different dynamical behavior depending on the parameter α. x1F is attracting for every value of α but the interval -3.25<α<-2.875. x2F is repelling for every value of α∈(-3.25,-0.098). The stability planes of Figure 6 illustrate this behavior. Consequently, for (i)
α∈(-∞,-3.25)∪(-0.098,+∞), x1F=x1∗ is attracting, and x2F is repelling,(ii)
α∈(-3.25,-2.875), x1F is repelling, and x2F=x2∗ is attracting,(iii)
α∈(-2.875,-0.098), x1F=x1∗, and x2F=x2∗ are attracting.Figure 6
Stability plane of (a)x1F and (b) x2F.
(a)
(b)From a stability point of view, just the first set of values ofα cause the family to have a good performance, since the only attracting point is the root of f(x).In a similar fashion as the casef(x)=x2-1, the infinity is a parabolic fixed point.The zero values of the derivative of (33) are the critical points, which correspond to the roots of polynomial (34)x8α2γ8+x714α2γ7+4αγ7+x677α2γ6+50αγ6+2γ6+x5224α2γ5+252αγ5+26γ5+x4383α2γ4+666αγ4+139γ4+x3398α2γ3+996αγ3+376γ3+x2247α2γ2+842αγ2+529γ2+x84α2γ+372αγ+360γ+12α2+66α+90. The polynomial has four real roots and four complex ones. No one of the critical points agrees with the roots of the function x2, so every one is free. Figure 7 shows the parameters planes of each real critical point. When the critical point tends to x1∗, x2∗, or x∞∗, the points are plotted in red, green, or blue, respectively, while black is kept for convergence to any other point. Note in Figures 7(b), 7(c), and 7(d) there is a thin black region, where the critical points have periodical orbits, so they do not converge to any one of the fixed points.Figure 7
Parameter planes of (a)x1C, (b) x2C, (c) x3C, and (d) x4C (red: x1∗, green: x2∗, blue: x∞∗, and black: others).
(a)
(b)
(c)
(d)For gathering the information of the parameters planes, the unified parameter plane is shown in Figure8.Figure 8
Unified parameter plane ofx1-4C.Finally, the dynamical behavior is analyzed with the dynamical line. The representation of the methods that converge to any of the roots is covered in Figure9, while Figure 10 illustrates the methods that also converge to any other point or have a periodical orbit. Color map is held.Figure 9
Dynamical lines of methods with (a)(α,γ)=(0.1124,-1), (b) (α,γ)=(-2,0.5), (c) (α,γ)=(1,1), and (d) (α,γ)=(-0.23,-1.5) (red: x1∗, green: x2∗, blue: x∞∗, and black: others).
(a)
(b)
(c)
(d)Figure 10
Dynamical lines of methods with (a)(α,γ)=(-2.876,1) and (b) (α,γ)=(-2.876,-1) (red: x1∗, green: x2∗, blue: x∞∗, and black: others).
(a)
(b)


### 4.2.3.f(x)=x2+1

Finally, the application of (32) on f(x)=x2+1 results in the fixed point operator M=N11(x)/D10(x), where N11(x) and D10(x) are polynomials of degrees 11 and 10, respectively.The fixed points are the roots off(x) and the roots of an eighth-degree polynomial. Both groups of roots are complex, so they are left out of the purpose of the current study.Regarding the critical points, an analogous comment can be made. The roots off(x) are critical points. The rest of the critical points are x3,4C=-1/γ±i and the roots of a polynomial of degree 12. No critical point is real, so the real dynamics can not be analyzed.


## 4.2.1.f(x)=x2-1

In this case, expression (32) is Mf(x)=N11(x)/D10(x), where N11(x) and D10(x) are polynomials of degrees 11 and 10, respectively, depending on x and on parameters α and γ.There are 10 fixed points ofMf(x). x1,2F=x1,2∗=±1 are superattracting points, while x3-10F(α,γ) are the roots of an 8th-degree polynomial: x3,4F∈R and x5-10F∈C. Since our purpose is the study of the real dynamics, x5-10F are rejected. The evaluation of Mf′x3,4F establishes the behavior of these points. For different values of α and γ, x3,4F have different dynamical features, as the stability plane of Figure 1 represents. A mesh of 100×100 points covers the values of α∈-2,2andγ∈[-2,2]. The white area represents where the multiplier is lower than 1, while the corresponding black represents where the multiplier is greater than 1. Its dynamical meaning is immediate, since white and black regions represent attracting and repelling behavior, respectively.Figure 1
Stability plane of (a)x3F and (b) x4F.
(a)
(b)Therefore,x3F has both behaviors, depending on the value of α and γ. For assuring the convergence of the method to x1,2∗, white areas must be avoided.The point∞ is a parabolic fixed point, since I(0)=0 and I′(0)=1.ComputingMf′(x)=0, 16 critical points can be found. As expected, x1,2C=x1,2∗. x3,4C=-1±γ/γ are preimages of x1,2∗.  x5-16C are the roots of a polynomial of degree 12, that depend on the parameters α and γ. Six of these roots are real (x5-10C) while the other ones are complex (x11-16C).A classical result (attributed to Fatou and Julia) establishes that there is at least one critical point associated with each invariant Fatou component, so it is interesting to analyze the behavior of each free critical point (critical point different to the roots of the polynomial) used as initial estimation for all the elements of the family.Just as the convergence plane [4] represents the final state of the orbits of every starting point in a single-parametric family of methods, we use the parameters plane to guess the final state of the orbit of one critical point in a biparametric family of methods, in a similar manner as [16]. Therefore, the parameters plane represents the family of methods Mf(x), where each point (α,γ) stands for an individual method. If the orbit of the free critical point xC tends to -1,1, or ∞, the point is colored in red, green, or blue, respectively. If xC does not converge to any of those points, the point is colored in black. Figure 2 shows the parameter planes of x5-10C.Figure 2
Parameter planes of (a)x5C, (b) x6C, (c) x7C, (d) x8C, (e) x9C, and (f) x10C (red: x1∗, green: x2∗, blue: x∞∗, and black: others).
(a)
(b)
(c)
(d)
(e)
(f)The unified parameter plane represented in Figure3 gathers in one image the main information of Figure 2. It is composed of the superposition of the black regions. In this way, the election of a point in the white region guarantees that the corresponding scheme tends to ±1 or ∞, while the set of methods in the black region do not.Figure 3
Unified parameter plane ofx5-10C.In order to visualize the behavior of each method, a new representation concept is introduced: the dynamical line. Analogous to the dynamical plane for complex variable, the dynamical line represents the basins of attraction, plotting in different colors where the orbit of each initial estimation tends. In coherence to Figure2, the red color is assigned to those initial approximations that tend to x1∗=-1. Analogously, the green color is assigned to x2∗=1, the blue color is assigned to x∗=∞, and the black one is assigned to those initial guesses that tend to a point different than the previous three. The final orbit of 1001 initial values of x∈[-2,2] has been calculated, and the convergence to a point is set when its distance is lower than 10-3.Figure4 represents the dynamical line of a set of methods whose orbits tend to -1,1, or ∞, that is, methods that belong to the white area of the unified parameter plane. On the contrary, Figure 5 shows methods whose points also tend to another fixed point.Figure 4
Dynamical lines of methods with (a)(α,γ)=(0.5,0.05), (b) (α,γ)=(0.5,1), (c) (α,γ)=(0.5,2), (d) (α,γ)=(1,0.05), (e) (α,γ)=(1,1), (f) (α,γ)=(1,2), (g) (α,γ)=(2,0.05), (h) (α,γ)=(2,1), and (i) (α,γ)=(2,2) (red: x1∗, green: x2∗, and blue: x∞∗).
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)Figure 5
Dynamical lines of methods with (a)(α,γ)=(-0.5,-0.5), (b) (α,γ)=(-0.5,0.5), (c) (α,γ)=(-1,-0.5), and (d) (α,γ)=(-1,0.5) (red: x1∗, green: x2∗, blue: x∞∗, and black: others).
(a)
(b)
(c)
(d)Note that every point in Figure4 tends to one of the attracting points. However, when a value of (α,γ) of the black region of Figure 3 is chosen, some points tend to a different point. For instance, these points can be found in Figures 5(a)–5(d) for x0=-0.2, x0=0.2, x0=-0.28, and x0=0.28, respectively.As expected from the attracting area ofx3F and the unified parameter plane, the election of methods with α>0 is capital, as proved in Figures 4 and 5.


## 4.2.2.f(x)=x2

Applying (32) over f(x)=x2, the fixed point operator is(33)Mx=x1+γx5+5γx+γ2x2+α2+γx1+γx22+γx3+γx3+2γx+α1+γx2.The real fixed points of the operator are x1F=0 and x2F(α,γ), while the other two fixed points are complex. Both real fixed points have a different dynamical behavior depending on the parameter α. x1F is attracting for every value of α but the interval -3.25<α<-2.875. x2F is repelling for every value of α∈(-3.25,-0.098). The stability planes of Figure 6 illustrate this behavior. Consequently, for (i)
α∈(-∞,-3.25)∪(-0.098,+∞), x1F=x1∗ is attracting, and x2F is repelling,(ii)
α∈(-3.25,-2.875), x1F is repelling, and x2F=x2∗ is attracting,(iii)
α∈(-2.875,-0.098), x1F=x1∗, and x2F=x2∗ are attracting.Figure 6
Stability plane of (a)x1F and (b) x2F.
(a)
(b)From a stability point of view, just the first set of values ofα cause the family to have a good performance, since the only attracting point is the root of f(x).In a similar fashion as the casef(x)=x2-1, the infinity is a parabolic fixed point.The zero values of the derivative of (33) are the critical points, which correspond to the roots of polynomial (34)x8α2γ8+x714α2γ7+4αγ7+x677α2γ6+50αγ6+2γ6+x5224α2γ5+252αγ5+26γ5+x4383α2γ4+666αγ4+139γ4+x3398α2γ3+996αγ3+376γ3+x2247α2γ2+842αγ2+529γ2+x84α2γ+372αγ+360γ+12α2+66α+90. The polynomial has four real roots and four complex ones. No one of the critical points agrees with the roots of the function x2, so every one is free. Figure 7 shows the parameters planes of each real critical point. When the critical point tends to x1∗, x2∗, or x∞∗, the points are plotted in red, green, or blue, respectively, while black is kept for convergence to any other point. Note in Figures 7(b), 7(c), and 7(d) there is a thin black region, where the critical points have periodical orbits, so they do not converge to any one of the fixed points.Figure 7
Parameter planes of (a)x1C, (b) x2C, (c) x3C, and (d) x4C (red: x1∗, green: x2∗, blue: x∞∗, and black: others).
(a)
(b)
(c)
(d)For gathering the information of the parameters planes, the unified parameter plane is shown in Figure8.Figure 8
Unified parameter plane ofx1-4C.Finally, the dynamical behavior is analyzed with the dynamical line. The representation of the methods that converge to any of the roots is covered in Figure9, while Figure 10 illustrates the methods that also converge to any other point or have a periodical orbit. Color map is held.Figure 9
Dynamical lines of methods with (a)(α,γ)=(0.1124,-1), (b) (α,γ)=(-2,0.5), (c) (α,γ)=(1,1), and (d) (α,γ)=(-0.23,-1.5) (red: x1∗, green: x2∗, blue: x∞∗, and black: others).
(a)
(b)
(c)
(d)Figure 10
Dynamical lines of methods with (a)(α,γ)=(-2.876,1) and (b) (α,γ)=(-2.876,-1) (red: x1∗, green: x2∗, blue: x∞∗, and black: others).
(a)
(b)


## 4.2.3.f(x)=x2+1

Finally, the application of (32) on f(x)=x2+1 results in the fixed point operator M=N11(x)/D10(x), where N11(x) and D10(x) are polynomials of degrees 11 and 10, respectively.The fixed points are the roots off(x) and the roots of an eighth-degree polynomial. Both groups of roots are complex, so they are left out of the purpose of the current study.Regarding the critical points, an analogous comment can be made. The roots off(x) are critical points. The rest of the critical points are x3,4C=-1/γ±i and the roots of a polynomial of degree 12. No critical point is real, so the real dynamics can not be analyzed.


## 5. Memory Dynamics

Introducing memory in the iterative expression of a method can improve its order of convergence. Below some fundamentals about methods with memory are introduced [11, 17].


### 5.1. Basics on Iterative Methods with Memory

The standard form of an iterative method with memory of first order is(35)xk+1=ϕxk,xk-1,k≥1,where x0 and x1 are the initial estimations. Fixed points must satisfy two conditions: xk+1=xk and xk-1=xk, so the fixed point function Φ can be defined by means of (36)Φxk-1,xk=xk,xk+1=xk,ϕxk-1,xk,k=1,2,….So, (xk-1,xk) is a fixed point of Φ if Φ(xk-1,xk)=(xk-1,xk) and, consequently, xk+1=xk and xk=xk-1. Therefore, the discrete dynamical system Φ:R2→R2, such that (37)Φx→=Φz,x=x,ϕz,x,where ϕ is the operator of the iterative scheme with memory, has been defined. Below some definitions related to iterative methods with memory are recalled.Theorbit of a pointx→0 is defined as the set {x→,Φ(x→),Φ2(x→),…,Φn(x→),…}. A point x→0=(z,x)0 is afixed pointx→0F=(z,x)0F of Φ if z=x and x=ϕ(z,x). If a fixed point x→F of operator Φ is different from (r,r), where f(r)=0, it is calledstrange fixed point. A point x→T is T-periodic if ΦT(x→T)=x→T and Φt(x→T)≠x→T, for t<T.The dynamical behavior of a pointx→ is defined from its asymptotical behavior. For this purpose, Theorem  5.1 of [18] is recalled.Theorem 5.
LetΦ from Rn to Rn be C2. Assume x→T is T-periodic. Let λ1,λ2,…,λn be the eigenvalues of Φ′(x→T). Then (1)
if all the eigenvaluesλj have λj<1, then x→T is attracting;(2)
if one eigenvalueλj0 has λj0>1, then x→T is unstable, that is, repelling or saddle;(3)
if all the eigenvaluesλj have λj>1, then x→T is repelling.In addition,(1)
if all the eigenvaluesλj have λj≠1, the fixed point ishyperbolic;(2)
if there exist an eigenvalueλi:λi<1 and another eigenvalue λj:λj>1, then the hyperbolic point is recalled assaddle point.A critical pointx→C satisfies det(Φ′(x→C))=0. The basin of attraction of a fixed attracting point, A(x→∗), is defined as the set of preimages of any order such that (38)Ax→∗=x→0∈Rn:Φmx→0⟶x→∗,m⟶∞.The drawing tool to represent the basins of attraction is the dynamical plane [16, 19]. This tool is based on the application to complex dynamics, where horizontal and vertical axes stand for the real and imaginary parts of a complex value, respectively. For real dynamics with memory, the horizontal axis is devoted to the current iteration xk while the vertical one represents the last iteration xk-1.


### 5.2. Stability of Family MM1(α)

As we have seen in Section3, the memory is introduced in the iterative expression by replacing parameter γ by -1/f′¯(x∗), where f′¯(x∗) denotes any approximation of f′(x∗). The first family of iterative methods with memory designed in Section 3 is expressed as(39)wk=xk-fxkfxk,xk-1,yk=xk-fxkfxk,wk,xk+1=yk-fxk+αfykfxk+α-1fykfykfyk,wk,where α is a free parameter. We are going to analyze the dynamical behavior of this class on the quadratic polynomials x2-1, x2, and x2+1.


#### 5.2.1.f(x)=x2-1

When (39) is applied over f(x)=x2-1, the fixed point operator is (40)Mαz,x=-x2-12z2-1x+zax2-1z2-1+x2+2xz+12/2x3z+3x2z2+1+6xz+z2+1ax2-1z2-1+x+zxx2+3xz+3+z+x2z+2x+zx2+2xz+1-1.The fixed points of the operator are(z,x)1,2,3F={(-1,-1),(0,0),(1,1)}. The fixed points (z,x)1F=(z,x)1∗=(-1,-1) and (z,x)3F=(z,x)3∗=(1,1) are superattracting, while the behavior of (z,x)2F depends on the parameter α. In this way, (z,x)2F=(z,x)2∗ is attracting just in case -1/2<α≤-1/17(5+22).Figure11 represents the dynamical plane of the method when (z,x)2F is attracting (a) and when it is not (b). White stars plot the attracting points. Orange points represent the set of initial points that converge to (z,x)1∗=(-1,-1), while green points are the (z,x)3∗=(1,1) corresponding ones. Blue is devoted to initial estimations that finally converge to (z,x)2∗=(0,0). The initial points are a mesh of 501×501 points of values of x and z from −2 to 2. Every initial point has converged to one of the attracting fixed points in less than 100 iterations. Further information about the representation can be found in [16].Figure 11
Dynamical planes of methodMα(z,x) for f(x)=x2-1.
(a)
α
=
-
0.47, and z2F is attracting
(b)
α
=
0.5,   and z2F is repelling


#### 5.2.2.f(x)=x2

The application of (39) over f(x)=x2 results in the following fixed point operator:(41)Mαz,x=xzx3+6x2z+αxz2+10xz2+2αz3+5z3x+2z2x+3zx2+4xz+αz2+3z2.The fixed points are (z,x)1F=(0,0) and (z,x)∞F=∞. Some singularities can be found when the denominator of (41) is zero, for z=-x/2,z=-2x/3, and α=-8. This behavior is illustrated in Figure 12, where two dynamical planes are represented for α≠-8 (a) and α=-8 (b).Figure 12
Dynamical planes of methodMα(z,x) for f(x)=x2.
(a)
α
=
-
1
(b)
α
=
-
8


#### 5.2.3.f(x)=x2+1

Finally, we apply (39) on f(x)=x2+1. The fixed point operator is a rational function whose numerator has degree 9 and its denominator has degree 8. Note the degrees are the higher sum of the powers of x and z. Let r1-8 be the roots of p(x)=x8(6α+49)+x6(4α-108)+x4(86-8α)+x2(-4α-12)+2α+1. Depending on the region of α, the amount of fixed points changes. Table 1 summarizes the fixed points depending on the value of α.Table 1
Fixed points depending onα.
α
Fixed points
<
-
8.1667
r
1
,
r
1
,
(
r
2
,
r
2
)
,
(
r
3
,
r
3
)
,
(
r
4
,
r
4
)
,
(
0,0
)
∈
(
-
8.1667
,
-
2.61803
)
(
r
1
,
r
1
)
,
(
r
2
,
r
2
)
,
(
r
3
,
r
3
)
,
(
r
4
,
r
4
)
,
(
r
5
,
r
5
)
,
(
r
6
,
r
6
)
,
(
0,0
)
=
-
2.61803
(
±
1.26518
,
±
1.26518
)
,
(
±
0.488916
,
±
0.488916
)
,
(
0,0
)
∈
(
-
2.61803
,
-
2.57392
)
(
r
1
,
r
1
)
,
(
r
2
,
r
2
)
,
(
r
3
,
r
3
)
,
(
r
4
,
r
4
)
,
(
r
5
,
r
5
)
,
(
r
6
,
r
6
)
,
(
0,0
)
=
-
2.57392
(
±
1.31541
,
±
1.31541
)
,
(
±
0.487946
,
±
0.487946
)
,
(
0,0
)
∈
(
-
2.57392
,
-
0.5
)
(
r
1
,
r
1
)
,
(
r
2
,
r
2
)
,
(
0,0
)
=
-
0.5
(
±
0.361815
,
±
0.361815
)
,
(
0,0
)
∈
(
-
0.5
,
-
0.381966
)
(
r
1
,
r
1
)
,
(
r
2
,
r
2
)
,
(
r
3
,
r
3
)
,
(
r
4
,
r
4
)
,
(
0,0
)
=
-
0.381966
(
±
0.172745
,
±
0.172745
)
,
(
0,0
)
∈
(
-
0.381966
,
-
0.324219
)
(
r
1
,
r
1
)
,
(
r
2
,
r
2
)
,
(
r
3
,
r
3
)
,
(
r
4
,
r
4
)
,
(
0,0
)
=
-
0.324219
(
±
0.262345
,
±
0.262345
)
,
(
0,0
)
∈
(
-
0.324219,0
)
(
0,0
)
∈
(
0,4.73341
)
(
0,0
)
=
4.73341
(
±
0.788301
,
±
0.788301
)
,
(
0,0
)
>
4.73341
(
r
1
,
r
1
)
,
(
r
2
,
r
2
)
,
(
r
3
,
r
3
)
,
(
r
4
,
r
4
)
,
(
0,0
)Every fixed point whose value is known is repelling, except the cases collected in Table2. Related to Theorem 5, only the point (0,0) is an attracting fixed point for α=-0.381966; the rest of the points are nonhyperbolic.Table 2
Fixed points without repelling behavior.
α
Fixed points
Eigenvalues
-
0.5
(
0,0
)
(
0,1
)
-
0.381966
(
0,0
)
0.190983
+
0.7626
i
(
-
1,1
)
,
λ
1,2
=
0.786151
-
0.324219
(
±
0.262345
,
±
0.262345
)
(
0.438608,1
)
-
0.324219
(
0,0
)
-
0.0421652
+
1.04046
i
(
-
1,1
)
,
λ
1,2
=
1.04131
4.73341
(
±
0.788301
,
±
0.788301
)
(
-
0.613745,1
)The dynamical planes of the cases in Table2 are represented in Figures 13(a)–13(d). In this case, white stars represent the fixed points. Only for the particular case of α=-0.381966 is there an attracting fixed point that attracts a few points of the memory plane, as can be observed in Figure 13.Figure 13
((a)–(d)) Dynamical planes of methodMα(z,x) for f(x)=x2+1. ((e)–(h)) Orbit of points xk for 5000 iterations; inset: orbit at iterations 4000 to 4020.
(a)
α
=
-
0.5
(b)
α
=
-
0.381966
(c)
α
=
-
0.324219
(d)
α
=
4.73341
(e)
α
=
-
0.5
(f)
α
=
-
0.381966
(g)
α
=
-
0.324219
(h)
α
=
4.73341Since the usual behavior is the nonconvergence to any of the roots, the orbit of different initial points have been obtained and represented in Figures13(e)–13(h).Forα∈{-0.5,-0.324219,4.73341}, whose corresponding Figures are (e), (g), and (h), the orbit of the analyzed points falls into a periodical orbit, as the magnification in their insets shows.


## 5.1. Basics on Iterative Methods with Memory

The standard form of an iterative method with memory of first order is(35)xk+1=ϕxk,xk-1,k≥1,where x0 and x1 are the initial estimations. Fixed points must satisfy two conditions: xk+1=xk and xk-1=xk, so the fixed point function Φ can be defined by means of (36)Φxk-1,xk=xk,xk+1=xk,ϕxk-1,xk,k=1,2,….So, (xk-1,xk) is a fixed point of Φ if Φ(xk-1,xk)=(xk-1,xk) and, consequently, xk+1=xk and xk=xk-1. Therefore, the discrete dynamical system Φ:R2→R2, such that (37)Φx→=Φz,x=x,ϕz,x,where ϕ is the operator of the iterative scheme with memory, has been defined. Below some definitions related to iterative methods with memory are recalled.Theorbit of a pointx→0 is defined as the set {x→,Φ(x→),Φ2(x→),…,Φn(x→),…}. A point x→0=(z,x)0 is afixed pointx→0F=(z,x)0F of Φ if z=x and x=ϕ(z,x). If a fixed point x→F of operator Φ is different from (r,r), where f(r)=0, it is calledstrange fixed point. A point x→T is T-periodic if ΦT(x→T)=x→T and Φt(x→T)≠x→T, for t<T.The dynamical behavior of a pointx→ is defined from its asymptotical behavior. For this purpose, Theorem  5.1 of [18] is recalled.Theorem 5.
LetΦ from Rn to Rn be C2. Assume x→T is T-periodic. Let λ1,λ2,…,λn be the eigenvalues of Φ′(x→T). Then (1)
if all the eigenvaluesλj have λj<1, then x→T is attracting;(2)
if one eigenvalueλj0 has λj0>1, then x→T is unstable, that is, repelling or saddle;(3)
if all the eigenvaluesλj have λj>1, then x→T is repelling.In addition,(1)
if all the eigenvaluesλj have λj≠1, the fixed point ishyperbolic;(2)
if there exist an eigenvalueλi:λi<1 and another eigenvalue λj:λj>1, then the hyperbolic point is recalled assaddle point.A critical pointx→C satisfies det(Φ′(x→C))=0. The basin of attraction of a fixed attracting point, A(x→∗), is defined as the set of preimages of any order such that (38)Ax→∗=x→0∈Rn:Φmx→0⟶x→∗,m⟶∞.The drawing tool to represent the basins of attraction is the dynamical plane [16, 19]. This tool is based on the application to complex dynamics, where horizontal and vertical axes stand for the real and imaginary parts of a complex value, respectively. For real dynamics with memory, the horizontal axis is devoted to the current iteration xk while the vertical one represents the last iteration xk-1.


## 5.2. Stability of Family MM1(α)

As we have seen in Section3, the memory is introduced in the iterative expression by replacing parameter γ by -1/f′¯(x∗), where f′¯(x∗) denotes any approximation of f′(x∗). The first family of iterative methods with memory designed in Section 3 is expressed as(39)wk=xk-fxkfxk,xk-1,yk=xk-fxkfxk,wk,xk+1=yk-fxk+αfykfxk+α-1fykfykfyk,wk,where α is a free parameter. We are going to analyze the dynamical behavior of this class on the quadratic polynomials x2-1, x2, and x2+1.


### 5.2.1.f(x)=x2-1

When (39) is applied over f(x)=x2-1, the fixed point operator is (40)Mαz,x=-x2-12z2-1x+zax2-1z2-1+x2+2xz+12/2x3z+3x2z2+1+6xz+z2+1ax2-1z2-1+x+zxx2+3xz+3+z+x2z+2x+zx2+2xz+1-1.The fixed points of the operator are(z,x)1,2,3F={(-1,-1),(0,0),(1,1)}. The fixed points (z,x)1F=(z,x)1∗=(-1,-1) and (z,x)3F=(z,x)3∗=(1,1) are superattracting, while the behavior of (z,x)2F depends on the parameter α. In this way, (z,x)2F=(z,x)2∗ is attracting just in case -1/2<α≤-1/17(5+22).Figure11 represents the dynamical plane of the method when (z,x)2F is attracting (a) and when it is not (b). White stars plot the attracting points. Orange points represent the set of initial points that converge to (z,x)1∗=(-1,-1), while green points are the (z,x)3∗=(1,1) corresponding ones. Blue is devoted to initial estimations that finally converge to (z,x)2∗=(0,0). The initial points are a mesh of 501×501 points of values of x and z from −2 to 2. Every initial point has converged to one of the attracting fixed points in less than 100 iterations. Further information about the representation can be found in [16].Figure 11
Dynamical planes of methodMα(z,x) for f(x)=x2-1.
(a)
α
=
-
0.47, and z2F is attracting
(b)
α
=
0.5,   and z2F is repelling


### 5.2.2.f(x)=x2

The application of (39) over f(x)=x2 results in the following fixed point operator:(41)Mαz,x=xzx3+6x2z+αxz2+10xz2+2αz3+5z3x+2z2x+3zx2+4xz+αz2+3z2.The fixed points are (z,x)1F=(0,0) and (z,x)∞F=∞. Some singularities can be found when the denominator of (41) is zero, for z=-x/2,z=-2x/3, and α=-8. This behavior is illustrated in Figure 12, where two dynamical planes are represented for α≠-8 (a) and α=-8 (b).Figure 12
Dynamical planes of methodMα(z,x) for f(x)=x2.
(a)
α
=
-
1
(b)
α
=
-
8


### 5.2.3.f(x)=x2+1

Finally, we apply (39) on f(x)=x2+1. The fixed point operator is a rational function whose numerator has degree 9 and its denominator has degree 8. Note the degrees are the higher sum of the powers of x and z. Let r1-8 be the roots of p(x)=x8(6α+49)+x6(4α-108)+x4(86-8α)+x2(-4α-12)+2α+1. Depending on the region of α, the amount of fixed points changes. Table 1 summarizes the fixed points depending on the value of α.Table 1
Fixed points depending onα.
α
Fixed points
<
-
8.1667
r
1
,
r
1
,
(
r
2
,
r
2
)
,
(
r
3
,
r
3
)
,
(
r
4
,
r
4
)
,
(
0,0
)
∈
(
-
8.1667
,
-
2.61803
)
(
r
1
,
r
1
)
,
(
r
2
,
r
2
)
,
(
r
3
,
r
3
)
,
(
r
4
,
r
4
)
,
(
r
5
,
r
5
)
,
(
r
6
,
r
6
)
,
(
0,0
)
=
-
2.61803
(
±
1.26518
,
±
1.26518
)
,
(
±
0.488916
,
±
0.488916
)
,
(
0,0
)
∈
(
-
2.61803
,
-
2.57392
)
(
r
1
,
r
1
)
,
(
r
2
,
r
2
)
,
(
r
3
,
r
3
)
,
(
r
4
,
r
4
)
,
(
r
5
,
r
5
)
,
(
r
6
,
r
6
)
,
(
0,0
)
=
-
2.57392
(
±
1.31541
,
±
1.31541
)
,
(
±
0.487946
,
±
0.487946
)
,
(
0,0
)
∈
(
-
2.57392
,
-
0.5
)
(
r
1
,
r
1
)
,
(
r
2
,
r
2
)
,
(
0,0
)
=
-
0.5
(
±
0.361815
,
±
0.361815
)
,
(
0,0
)
∈
(
-
0.5
,
-
0.381966
)
(
r
1
,
r
1
)
,
(
r
2
,
r
2
)
,
(
r
3
,
r
3
)
,
(
r
4
,
r
4
)
,
(
0,0
)
=
-
0.381966
(
±
0.172745
,
±
0.172745
)
,
(
0,0
)
∈
(
-
0.381966
,
-
0.324219
)
(
r
1
,
r
1
)
,
(
r
2
,
r
2
)
,
(
r
3
,
r
3
)
,
(
r
4
,
r
4
)
,
(
0,0
)
=
-
0.324219
(
±
0.262345
,
±
0.262345
)
,
(
0,0
)
∈
(
-
0.324219,0
)
(
0,0
)
∈
(
0,4.73341
)
(
0,0
)
=
4.73341
(
±
0.788301
,
±
0.788301
)
,
(
0,0
)
>
4.73341
(
r
1
,
r
1
)
,
(
r
2
,
r
2
)
,
(
r
3
,
r
3
)
,
(
r
4
,
r
4
)
,
(
0,0
)Every fixed point whose value is known is repelling, except the cases collected in Table2. Related to Theorem 5, only the point (0,0) is an attracting fixed point for α=-0.381966; the rest of the points are nonhyperbolic.Table 2
Fixed points without repelling behavior.
α
Fixed points
Eigenvalues
-
0.5
(
0,0
)
(
0,1
)
-
0.381966
(
0,0
)
0.190983
+
0.7626
i
(
-
1,1
)
,
λ
1,2
=
0.786151
-
0.324219
(
±
0.262345
,
±
0.262345
)
(
0.438608,1
)
-
0.324219
(
0,0
)
-
0.0421652
+
1.04046
i
(
-
1,1
)
,
λ
1,2
=
1.04131
4.73341
(
±
0.788301
,
±
0.788301
)
(
-
0.613745,1
)The dynamical planes of the cases in Table2 are represented in Figures 13(a)–13(d). In this case, white stars represent the fixed points. Only for the particular case of α=-0.381966 is there an attracting fixed point that attracts a few points of the memory plane, as can be observed in Figure 13.Figure 13
((a)–(d)) Dynamical planes of methodMα(z,x) for f(x)=x2+1. ((e)–(h)) Orbit of points xk for 5000 iterations; inset: orbit at iterations 4000 to 4020.
(a)
α
=
-
0.5
(b)
α
=
-
0.381966
(c)
α
=
-
0.324219
(d)
α
=
4.73341
(e)
α
=
-
0.5
(f)
α
=
-
0.381966
(g)
α
=
-
0.324219
(h)
α
=
4.73341Since the usual behavior is the nonconvergence to any of the roots, the orbit of different initial points have been obtained and represented in Figures13(e)–13(h).Forα∈{-0.5,-0.324219,4.73341}, whose corresponding Figures are (e), (g), and (h), the orbit of the analyzed points falls into a periodical orbit, as the magnification in their insets shows.


## 5.2.1.f(x)=x2-1

When (39) is applied over f(x)=x2-1, the fixed point operator is (40)Mαz,x=-x2-12z2-1x+zax2-1z2-1+x2+2xz+12/2x3z+3x2z2+1+6xz+z2+1ax2-1z2-1+x+zxx2+3xz+3+z+x2z+2x+zx2+2xz+1-1.The fixed points of the operator are(z,x)1,2,3F={(-1,-1),(0,0),(1,1)}. The fixed points (z,x)1F=(z,x)1∗=(-1,-1) and (z,x)3F=(z,x)3∗=(1,1) are superattracting, while the behavior of (z,x)2F depends on the parameter α. In this way, (z,x)2F=(z,x)2∗ is attracting just in case -1/2<α≤-1/17(5+22).Figure11 represents the dynamical plane of the method when (z,x)2F is attracting (a) and when it is not (b). White stars plot the attracting points. Orange points represent the set of initial points that converge to (z,x)1∗=(-1,-1), while green points are the (z,x)3∗=(1,1) corresponding ones. Blue is devoted to initial estimations that finally converge to (z,x)2∗=(0,0). The initial points are a mesh of 501×501 points of values of x and z from −2 to 2. Every initial point has converged to one of the attracting fixed points in less than 100 iterations. Further information about the representation can be found in [16].Figure 11
Dynamical planes of methodMα(z,x) for f(x)=x2-1.
(a)
α
=
-
0.47, and z2F is attracting
(b)
α
=
0.5,   and z2F is repelling


## 5.2.2.f(x)=x2

The application of (39) over f(x)=x2 results in the following fixed point operator:(41)Mαz,x=xzx3+6x2z+αxz2+10xz2+2αz3+5z3x+2z2x+3zx2+4xz+αz2+3z2.The fixed points are (z,x)1F=(0,0) and (z,x)∞F=∞. Some singularities can be found when the denominator of (41) is zero, for z=-x/2,z=-2x/3, and α=-8. This behavior is illustrated in Figure 12, where two dynamical planes are represented for α≠-8 (a) and α=-8 (b).Figure 12
Dynamical planes of methodMα(z,x) for f(x)=x2.
(a)
α
=
-
1
(b)
α
=
-
8


## 5.2.3.f(x)=x2+1

Finally, we apply (39) on f(x)=x2+1. The fixed point operator is a rational function whose numerator has degree 9 and its denominator has degree 8. Note the degrees are the higher sum of the powers of x and z. Let r1-8 be the roots of p(x)=x8(6α+49)+x6(4α-108)+x4(86-8α)+x2(-4α-12)+2α+1. Depending on the region of α, the amount of fixed points changes. Table 1 summarizes the fixed points depending on the value of α.Table 1
Fixed points depending onα.
α
Fixed points
<
-
8.1667
r
1
,
r
1
,
(
r
2
,
r
2
)
,
(
r
3
,
r
3
)
,
(
r
4
,
r
4
)
,
(
0,0
)
∈
(
-
8.1667
,
-
2.61803
)
(
r
1
,
r
1
)
,
(
r
2
,
r
2
)
,
(
r
3
,
r
3
)
,
(
r
4
,
r
4
)
,
(
r
5
,
r
5
)
,
(
r
6
,
r
6
)
,
(
0,0
)
=
-
2.61803
(
±
1.26518
,
±
1.26518
)
,
(
±
0.488916
,
±
0.488916
)
,
(
0,0
)
∈
(
-
2.61803
,
-
2.57392
)
(
r
1
,
r
1
)
,
(
r
2
,
r
2
)
,
(
r
3
,
r
3
)
,
(
r
4
,
r
4
)
,
(
r
5
,
r
5
)
,
(
r
6
,
r
6
)
,
(
0,0
)
=
-
2.57392
(
±
1.31541
,
±
1.31541
)
,
(
±
0.487946
,
±
0.487946
)
,
(
0,0
)
∈
(
-
2.57392
,
-
0.5
)
(
r
1
,
r
1
)
,
(
r
2
,
r
2
)
,
(
0,0
)
=
-
0.5
(
±
0.361815
,
±
0.361815
)
,
(
0,0
)
∈
(
-
0.5
,
-
0.381966
)
(
r
1
,
r
1
)
,
(
r
2
,
r
2
)
,
(
r
3
,
r
3
)
,
(
r
4
,
r
4
)
,
(
0,0
)
=
-
0.381966
(
±
0.172745
,
±
0.172745
)
,
(
0,0
)
∈
(
-
0.381966
,
-
0.324219
)
(
r
1
,
r
1
)
,
(
r
2
,
r
2
)
,
(
r
3
,
r
3
)
,
(
r
4
,
r
4
)
,
(
0,0
)
=
-
0.324219
(
±
0.262345
,
±
0.262345
)
,
(
0,0
)
∈
(
-
0.324219,0
)
(
0,0
)
∈
(
0,4.73341
)
(
0,0
)
=
4.73341
(
±
0.788301
,
±
0.788301
)
,
(
0,0
)
>
4.73341
(
r
1
,
r
1
)
,
(
r
2
,
r
2
)
,
(
r
3
,
r
3
)
,
(
r
4
,
r
4
)
,
(
0,0
)Every fixed point whose value is known is repelling, except the cases collected in Table2. Related to Theorem 5, only the point (0,0) is an attracting fixed point for α=-0.381966; the rest of the points are nonhyperbolic.Table 2
Fixed points without repelling behavior.
α
Fixed points
Eigenvalues
-
0.5
(
0,0
)
(
0,1
)
-
0.381966
(
0,0
)
0.190983
+
0.7626
i
(
-
1,1
)
,
λ
1,2
=
0.786151
-
0.324219
(
±
0.262345
,
±
0.262345
)
(
0.438608,1
)
-
0.324219
(
0,0
)
-
0.0421652
+
1.04046
i
(
-
1,1
)
,
λ
1,2
=
1.04131
4.73341
(
±
0.788301
,
±
0.788301
)
(
-
0.613745,1
)The dynamical planes of the cases in Table2 are represented in Figures 13(a)–13(d). In this case, white stars represent the fixed points. Only for the particular case of α=-0.381966 is there an attracting fixed point that attracts a few points of the memory plane, as can be observed in Figure 13.Figure 13
((a)–(d)) Dynamical planes of methodMα(z,x) for f(x)=x2+1. ((e)–(h)) Orbit of points xk for 5000 iterations; inset: orbit at iterations 4000 to 4020.
(a)
α
=
-
0.5
(b)
α
=
-
0.381966
(c)
α
=
-
0.324219
(d)
α
=
4.73341
(e)
α
=
-
0.5
(f)
α
=
-
0.381966
(g)
α
=
-
0.324219
(h)
α
=
4.73341Since the usual behavior is the nonconvergence to any of the roots, the orbit of different initial points have been obtained and represented in Figures13(e)–13(h).Forα∈{-0.5,-0.324219,4.73341}, whose corresponding Figures are (e), (g), and (h), the orbit of the analyzed points falls into a periodical orbit, as the magnification in their insets shows.


## 6. Conclusions

We have designed a biparametric family of fourth-order optimal schemes for solving nonlinear equations. Its error equation allows us to introduce memory increasing the order of convergence until six, without need of more functional evaluations. We have considered in depth the dynamical study of this family. Several tools related to real and memory dynamics have been applied, supplying a very important visual information about the different behaviors. Moreover, the dynamical line and the unified parameter plane have been introduced for the first time to the best of our knowledge.The real dynamics of the biparametric family has revealed its good performance in terms of stability. For polynomials with complex roots, the iterative family does not fail, while for the rest of polynomials there are wide regions for both parameters where the stability of the family is guaranteed.The memory dynamics have improved the order of convergence of the original method. For polynomials with real roots, the iterative family of methods has a broad stable behavior. When the family is applied over polynomials with complex roots, many fixed points with repelling behavior appear, and some periodical orbits can be found.



---

*Source: 2713145-2017-10-31.xml*