# Dynamical Techniques for Analyzing Iterative Schemes with Memory

Authors: Neha Choubey; A. Cordero; J. P. Jaiswal; J. R. Torregrosa
Journal: Complexity
 (2018)

---

## Abstract

We construct a new biparametric three-point method with memory to highly improve the computational efficiency of its original partner, without adding functional evaluations. In this way, through different estimations of self-accelerating parameters, we have modified an existing seventh-order method. The parameters have been defined by Hermite interpolating polynomial that allows the accelerating effect. In particular, theR-order of the proposed iterative method with memory is increased from seven to ten. A real multidimensional analysis of the stability of this method with memory is made, in order to study its dependence on the initial estimations. Taking into account that usually iterative methods with memory are more stable than their derivative-free partners and the obtained results in this study, the behavior of this scheme shows to be excellent, but for a small domain. Numerical examples and comparison are also provided, confirming the theoretical results.


---

## Body

## 1. Introduction

In recent times, solving nonlinear problems described byfx=0 is burning difficulty in real-world phenomena. In this direction, numerous iterative methods have been projected (see, e.g., [1–7]). These iterative methods have a noteworthy area of research in numerical analysis, as they can be applied in several areas of pure or applied sciences. Out of them, the most eminent one-point iterative method without memory is the Newton-Raphson scheme, which is given by
(1)xn+1=xn−fxnf′xn,n=0,1,…and has second-order convergence. One drawback of this method is the clause f′xn≠0, which confines its applications.The first objective and inspiration to design iterative methods for solving this kind of problem are to get the highest order of convergence with the least computational cost. Therefore, a lot of researchers have paid much interest to construct optimal multipoint methods without memory, in the sense of Kung-Traub’s conjecture [8], using n+1 functional evaluation which can reach the optimal 2n. The paper follows three main goals: primarily, to avoid the restriction f′x≠0 in the practice; along with the subsequent goal, to develop the proposed families to methods with memory in an approach that achieves convergence R-order 10 without any new functional evaluation; and to analyze the dependence of the resulting scheme from the set of initial estimations used.Multipoint schemes have enormous practical meaning, as they conquer theoretical limits of any point-to-point scheme about order and efficiency. They also generate approximations of higher precision and highly improved computer arithmetics, and symbolic calculation has allowed efficient execution of multipoint methods.On the other hand, iterative schemes with memory utilize information from the recent and preceding iterations. Traub, in [9], designed the first method with memory by a small change in the well-known Steffensen’s scheme as follows:
(2)un=xn+ρnfxn,ρn≠0,xn+1=xn−fxnfxn,un,n≥0,where x0 and ρ0 are given; ρn is a self-accelerating parameter defined as
(3)ρn+1=−1N1′xn,N1x=fxn+x−xnfxn,un,n≥0,where N1x is Newton’s interpolating polynomial of the first degree. The order of convergence of (2) was 2.41.More recently, some authors have constructed iterative schemes with memory from optimal procedures of different orders, mainly four (see, e.g., [10–12]), eight ([13–15], among others), or even general n-point schemes [16, 17]. Some good reviews regarding the acceleration of convergence order by using memory are [18, 19].The convergence order of a new method is improved regarding its without memory partner; it is derived by adding self-accelerating parameters but holding the derivatives in the iterative expression. The accelerated convergence rate has been obtained without additional functional evaluations, which results in higher computational efficiency. However, another important aspect of an iterative method to be considered is the numerical stability, that is, the analysis that tells us how dependent the scheme of the initial guesses used is. The dynamical performance of the rational functions associated to iterative schemes is a very useful element to study their dependence on initial estimations. In recent years, complex discrete dynamics has been widely used on iterative methods without memory (see, e.g., [20–23]). Nevertheless, it is known that iterative schemes with memory cannot be analyzed by means of these techniques. This is the reason why the authors focused in [24–26] their qualitative study by transforming them into multidimensional dynamical systems, as their qualitative properties can be analyzed by using standard tools as the classical Hénon map (see [27, 28]).In this manuscript, we have obtained a multipoint method with memory, to solve the nonlinear equations followed by its convergence and stability analysis. Sections2, 3, and 4 can be summarized as follows: in Section 2, we design a biparametric three-point iterative scheme with memory. The method has been obtained by employing two self-accelerating parameters, which use current and previous data from available information. These parameters are recalculated at each iteration by using Hermite interpolating polynomials, increasing from 7 to 10 the order of convergence. Section 3 is related to the qualitative analysis of the proposed scheme in terms of the real discrete dynamical system involved. After introducing some basic concepts, the stability of the method on a low degree generic polynomial is studied. In general, it presents a very stable behavior, but also a chaotic one for the values of the parameter in a small domain. In Section 4, some numerical examples are presented confirming the results proven in Section 2.


## 2. Convergence Analysis

In this section, we demonstrate the improvement of the convergence order of the method given by Khattri and Argyros [29] by adding up two different parameters in the first and the third steps. Firstly, we consider the seventh order without memory scheme presented in [29].
(4)yn=xn−fxnf′xn,zn=yn−fyn−f′xn+2fyn−fxn/yn−xn,xn+1=zn−fznα1fyn−fxn/yn−xn+α2fzn−fyn/zn−yn+1−α1−α2fzn−fxn/zn−xn,where α1=−1 and α2=1. The error equation for the last step of (4) is
(5)en+1=c3−c22c22c3en7+Oen8,where cj=fjα/j!f′α, for j=2,3,…,n.By using two parametersγ and λ in (4) and replacing α1=−1 and α2=1, we obtain the biparametric family.
(6)yn=xn−fxnf′xn−γfxn,zn=yn−fyn−f′xn+2fyn−fxn/yn−xn,xn+1=zn−fzn−fyn−fxnyn−xn+fzn−fynzn−yn+fzn−fxnzn−xn+λzn−xnzn−yn.The error equations of each step of (6) are
(7)en,y=yn−α=c2−γen2+−2c22−γ2+2c2γ+2c3en3+⋯+Oen8,(8)en,z=zn−α=γ−c2γ−c2c2+c3en4+8γc23−4c24+3γc3−2c32+c22−6γ2+8c3+2c2γ3−5γc3−c4+2γc4en5+⋯+Oen8,and(9)en+1=γ−c22γ−c2c2+c3−λ+f′αc3f′αen7+Oen8.From (9), we can assure that the order of convergence of (6) is still seven, with the independence of parameters γ and λ. This order can be improved from seven to ten, by taking γ=c2 and λ=f′αc3, but root α is not known.So, to improve the rate of convergence of (6), we recalculate the value of parameters γ and λ in each iteration, by taking γ≈c2 and λ≈f′αc3, as f′α, f′′α, and f′′′α are not provided. We denote by γn and λn these estimations, and they are computed by using the current and previous iteration satisfying limn→∞γn=c2=f′′α/2f′α and limn→∞λn=f′αc3=f′′′α/6. We consider the following formulas for γn and λn, using Hermite interpolating polynomials of different degrees:
(10)γn=H5″xn2f′xn,λn=H6‴yn6.Theorem 1.
Let f:I⊆ℝ→ℝ be a sufficiently differentiable function and Hmx be the interpolating polynomial of Hermite with mth degree that interpolates fx at nodes yn,xn,xn,t0,…,tm−3 belonging to interval I, and the derivative fm+1 is continuous in I and the Hermite interpolating polynomial Hmx satisfies conditions Hmyn=fyn, Hmxn=fxn, Hm′xn=f′xn, and Hmtj=ftj,j=0,1,…,m−3.
Let us define et,j=tj−α and j=0,1,…,m−2 and
(i)
All nodes yn,xn,t0,…,tm−3 are close enough to zero α.
(ii)
Conditions en=xn−α=Oet,0…et,m−3 and en,y=yn−α=Oen2,et,0…et,m−3 hold. Then,
(11)λn−f′αc3~−1m−3f′αcm+1∏j=0m−3et,j.Proof.
It is known that the error expression ofHmx can be obtained by
(12)fx−Hmx=fm+1ξm+1!x−ynx−xn2∏j=0m−3x−tj,ξ∈I.After differentiating (12) at x=yn, we get
(13)f‴yn−Hm‴yn=6fm+1ξm+1!∏j=0m−3yn−tj+2yn−xn∑k=0m−3∏j≠k,j=0m−3yn−tj+yn−xn22∑k=0m−3∑i≠k,i=0m−3∏j≠k,j≠i,j=0m−3yn−tj,ξ∈I,Hm‴yn=f‴yn−6fm+1ξm+1!∏j=0m−3yn−tj+2yn−xn∑k=0m−3∏j≠k,j=0m−3yn−tj+yn−xn22∑k=0m−3∑i≠k,i=0m−3∏j≠k,j≠i,j=0m−3yn−tj,ξ∈I.Taylor’s development of the derivative off at xn,yn∈I and ξ∈I about α provides
(14)f′′′yn=f′α6c3+24c4en,y+Oen,y2and(15)fm+1ξ=f′αm+1!cm+1+m+2!cm+2eξ+Oeξ2,where eξ=ξ−α. Putting (14) and (15) in (13), we obtain
(16)Hm‴yn=6f′αc3+4c4en,y−cm+1∏j=0m−3yn−tj+2yn−xn∑k=0m−3∏j≠k,j=0m−3yn−tj+yn−xn22∑k=0m−3∑i≠k,i=0m−3∏j≠k,j≠i,j=0m−3yn−tj,which implies
(17)λn=Hm‴yn6~f′αc3+4c4en,y−cm+1∏j=0m−3yn−tj+2yn−xn∑k=0m−3∏j≠k,j=0m−3yn−tj+yn−xn22∑k=0m−3∑i≠k,i=0m−3∏j≠k,j≠i,j=0m−3yn−tj,that is,
(18)λn−f′αc3~f′α−cm+1∏j=0m−3yn−tj+2yn−xn∑k=0m−3∏j≠k,j=0m−3yn−tj+yn−xn22∑k=0m−3∑i≠k,i=0m−3∏j≠k,j≠i,j=0m−3yn−tj+4c4en,yor
(19)λn−f′αc3∼−1m−3f′αcm+1∏j=0m−3et,j.The concept of theR-order of convergence [30] and the subsequent declaration (see [4], p.287) is used to approximate the order of convergence of (6).Then, we can prove the following result.Theorem 2.
Let f:I⊆ℝ→ℝ be a sufficiently differentiable function and let γn and λn be the varying parameters in iterative scheme (6) obtained by means of (10). If the initial estimate x0 is close enough to a simple root α of fx, then, the R-order of the iterative method is at least 10. We denote this scheme by OM10.Proof.
Letxn be the generated sequence by scheme OM10. As it converges to the root α of fx, with R-order r, we write
(20)en+1~Dn,renr,where en=xn−α. If n→∞, then, Dn,r→Dr. Therefore,
(21)en+1~Dn,rDn−1,ren−1rr=Dn,rDn−1,rren−1r2.Let us consider now sequencesyn and zn with R-order p and q, respectively; then,
(22)en,y~Dn,penp~Dn,pDn−1,ren−1rp=Dn,pDn−1,rpen−1rp,and(23)en,z~Dn,qenq~Dn,qDn−1,ren−1rq=Dn,qDn−1,rqen−1rq.The following error expression of the method with memory (6) can be obtained by (7), (8), and (9) and the varying parameter γn and λn.
(24)en,y=yn−α~c2−γnen2,(25)en,z=zn−α~γn−c2γn−c2c2+c3en4,and(26)en+1=xn+1−α~γ−c22γ−c2c2+c3−λ+f′αc3f′αen7.Here, we excluded higher order terms in (24), (25), and (26).In method OM10,γn is calculated by (10): Hermite interpolating polynomial H5x satisfies conditions H5xn=fxn, H5′xn=f′xn, H5zn−1=fzn−1, H5yn−1=fyn−1, H5xn−1=fxn−1, and H5′xn−1=f′xn−1. After differentiating twice the error expression of H5x at point x=xn, we get
(27)H′5′xn=f′′xn−2f6ξ6!n−xn−12xn−zn−1xn−yn−1,ξ∈I.By using Taylor’s expansion off′ and f′′ at xn∈I and ξ∈I about α,
(28)f′xn=f′α1+2c2en+3c3en2+Oen3,(29)f′′xn=f′α2c2+6c3en+Oen2,and(30)fm+1ξ=f′αm+1!cm+1+m+2!cm+2eξ+Oeξ2.From (30), we have
(31)f6ξ=f′α6!c6+7!c7eξ+Oeξ2,where eξ=ξ−α. By replacing (29) and (31) into (27), we obtain
(32)H′5′xn∼2f′αc2−c6en−1,yen−1,zen−12+3c3en.Now, dividing (32) by (28),
(33)γn=H5''xn2f'xn~c2−c6en−1,yen−1,zen−12+3c3−2c22en,or(34)γn−c2~−c6en−1,yen−1,zen−12.Asλn is calculated by (10), the third derivative of Hermite interpolating polynomial H6x satisfies
(35)H′′6′yn=f′′′yn−f7ξ7!6yn−xn−12yn−yn−1yn−zn−1+16yn−xn−12yn−xnyn−zn−1+16yn−xn−12yn−xnyn−yn−1+16yn−xn2yn−xn−1yn−zn−1+16yn−xn2yn−xn−1yn−yn−1+24yn−xn−1yn−xnyn−yn−1yn−zn−1+12yn−xn−12yn−xn2+6yn−xn2yn−yn−1yn−zn−1,ξ∈I.From (30), we have
(36)f7ξ=f′α7!c7+8!c8eξ+Oeξ2.Substituting (14) and (36) into (35), we obtain
(37)H′′6′yn∼6f′αc3−c7en−1,yen−1,zen−12,and hence,(38)λn=H′′6′yn6∼f′αc3−c7en−1,yen−1,zen−12,or(39)λn−f′αc3~−f′αc7en−1,yen−1,zen−12.According to (24), (25), (26), (34), and (39), we get
(40)en,y~−c6en−1,yen−1,zen−12Dn−1,ren−1r2~−c6Dn−1,pDn−1,qDn−1,r2en−1pen−1qen−12ren−12~−c6Dn−1,pDn−1,qDn−1,r2en−1p+q+2r+2,(41)en,z~c6en−1,yen−1,zen−12Dn−1,ren−1r4~c6Dn−1,pDn−1,qDn−1,r4en−1pen−1qen−14ren−12~c6Dn−1,pDn−1,qDn−1,r4en−1p+q+4r+2and(42)en+1~c62en−1,yen−1,zen−122−c7f′αen−1,yen−1,zen−12Dn−1,ren−1r7~−c62c7f′αen−1,y3en−1,z3en−16Dn−1,r7en−17r~−f′αc62c7Dn−1,p3Dn−1,q3Dn−1,r7en−13pen−13qen−17ren−16~−f′αc62c7Dn−1,p3Dn−1,q3Dn−1,r7en−13p+3q+7r+6.By comparing the exponents ofen−1 featuring in the three pairs of relation (22)-(40), (23)-(41), and (21)-(42), the subsequent system is provided.
(43)2r+p+q+2=rp,4r+p+q+2=rq,7r+3p+3q+6=r2.The unique positive solution of (43) is r=10, q=5, and p=3. So R-order of method (6), for γn and λn being calculated by (10), is at least 10.


## 3. Multidimensional Dynamical Study

In this section, we build the discrete dynamical system associated to OM10, for overcoming its qualitative analysis. The general iterative expression of a fixed-point procedure that employs two previous iterations in order to get the following one is(44)xn+1=hxn−1,xn,n≥1,where x0 and x1 are initial approximations. But if the scheme is a three-step one, as OM10, not only xn−1 but also the subsequent intermediate points yn−1 and zn−1 must be used to calculate the following iterate xn with a high order of convergence. So the fixed point iteration is expressed as
(45)xn+1=gxn−1,yn−1,zn−1,xn,n≥1.To obtain the fixed points, we define thefixed point multidimensional function G:ℝ4→ℝ4 associated to g by means of
(46)Gxn−1,yn−1,zn−1,xn=xn,yn,zn,xn+1=xn,yn,zn,gxn−1,yn−1,zn−1,xn,n=1,2,…Therefore, a fixed point ofG satisfies xn+1=xn and xn−1=yn−1=zn−1=xn.We define a discrete dynamical system inℝ4 from function G:ℝ4→ℝ4, given by
(47)Gz,zy,zz,x=x,xy,xz,gz,xy,xz,x,where the steps of nth and (n−1)th iterations are denoted by z=xn−1, zy=yn−1, zz=zn−1, xy=yn, xz=zn, and x=xn. Fixed points z,zy,zz,x of G satisfy z=zy=zz=x and x=gz,xy,xz,x. By imposing these conditions to the rational operator G, it can be reduced to a real-valued function g~x, to study the asymptotic stability of these points.If all the components of a fixed point ofG are different from r, r being a zero of the nonlinear scalar function f, then, it is called strange fixed point.LetG:ℝ4→ℝ4 be a vectorial function. The orbit of a point x¯∈ℝ4 is defined as x¯,Gx¯,…,Gnx¯,….A pointx∗∈ℝ4 is a k-periodic point if Gkx∗=x∗ and Gpx∗≠x∗, for p=1,2,…,k−1.The stability of fixed points in multivariate operators (see, e.g., [31]) is analyzed as follows.Theorem 3.
Let G from ℝm to ℝm be of class C2. Assume x∗ is a k-periodic point. Let λ1,λ2,…,λm be the eigenvalues of G′x∗.
(a)
If all the eigenvaluesλj have ∣λj∣<1, then, x∗ is attracting.
(b)
If one eigenvalueλj0 has ∣λj0∣>1, then, x∗ is unstable, that is, repelling or saddle.
(c)
If all the eigenvaluesλj have ∣λj∣>1, then, x∗ is repelling.In addition, a fixed point ishyperbolic if λj satisfies λj≠1∀j. In particular, if λi exists such that ∣λi∣<1 and λj such that ∣λj∣>1, then, the hyperbolic point is a saddle point.Let us remark the difference between the stability of a fixed pointx∗ in one-dimensional dynamics and that in multidimensional dynamics: meanwhile, in the scalar case, the stability of the fixed point depends on the value of the derivative operator at the point: (∣R′x∗∣<1 means that x∗ is attracting, superattracting if ∣R′x∗∣=0 and repulsive if ∣R′x∗∣>1, R being the rational function associated to the iterative method applied on a polynomial px); in the multidimensional case, the eigenvalues of the Jacobian matrix associated to the fixed point operator are the elements that determine the character of the fixed points.Moreover, a pointx is called critical if the associated Jacobian matrix G′x satisfies detG′x=0 (in the one-dimensional case, a critical point is that what makes the derivative fixed point operator vanish). One way to calculate critical points, for iterative methods of convergence order higher than two, is to find those fixed points with null eigenvalues λj=0, ∀j. In the case that they are also fixed points, they are called superattracting, as an extension of the scalar case.So,x∗ being an attracting fixed point of G, its basin of attraction Ax∗ is defined as
(48)Ax∗=x0∈ℝn:Gmx0→x∗,m→∞.A key result from Julia and Fatou [32] proves the relationship between the existence of free critical points (they are called in this way if they are different from the roots of px) and the set of initial points that converge to an attracting periodic point: there is at least one critical point associated with each invariant Fatou component. So observing the orbit of the critical points, all the attracting behaviors can be found. This result is valid as for complex as for real iterative functions. The main drawback is that often the analytic expression of the critical points cannot be found in high-order methods, because of the elevated degree of the polynomials involved in the rational function.The union of all the basins of attraction defines theFatou set, usually represented as a dynamical plane. It is numerically constructed starting with a mesh of initial guesses, iterating the method on them and assigning different colors depending on the basin they converge to.


### 3.1. Dynamical Analysis of OM10

Now, we study the behavior of the operator associated to scheme OM10 on quadratic polynomials. In order to get this aim, the asymptotic stability of the fixed points of the associate rational vector applied onpx=x2−c, with c∈ℝ, must be analyzed. It is clear that, if c>0, px has two simple real roots, it has one double root at x=0 if c=0 and px does not have real roots if c<0.We need to calculate the fixed points of the associate fixed point operator on method OM10 onpx=x2−c, Mz,zy,zz,x,c, that is, a high-degree rational function depending on the parameter of the polynomial c. However, it is not viable to analyze directly the fixed points of this operator, as usually indetermination appears when we calculate the fixed points due to cancelations of factors as z−x, zy−x, ansd zz−x in the denominator of Mz,zy,zz,x,c. In order to avoid this and using that all the fixed points have equal components, we force firstly zy=zz=x and, after simplifying, z=x. The resulting one-dimensional reduced operator is
(49)Mz,zy,zz,x,c∣zy=zz=x∣z=x=m~x,c=x2x−12+c6x2−6x+1+c222x−1c+xx−1.Moreover, there are two values ofc which reduce the rational function:
(i)
Ifc=0, then,
(50)m~x,0=−1+xx−2+4x.
(ii)
Ifc=1/4, then,
(51)m~x,14=5−4x+4x2−8+16x.Proposition 1.
The number of real fixed points (and their stability) of m~x,c depends on the value of parameter c:
(a)
If c<0, px has no real roots, so fixed points s1c=2/3−1/31−3c and s2c=2/3+1/31−3c are strange, being both repulsive.
(b)
If c=0, the only (double) root of px, x=0 is an attracting fixed point of m~x,0 and also x=1/3 is a repulsive strange fixed point.
(c)
When 0<c<1/4 or c>1/4, fixed points c,c and −c,−c are superattracting (corresponding to the roots of px). Moreover, s1c and s2c are strange fixed points if 1/4<c≤1/3; s1c is always repulsive, and s2c is attracting if c∈7/25,1/3 and superattracting if c=5/16.
For c=1/4, there is only one (superattracting) fixed point that is a root of px, x=−1/2 and one strange fixed point, x=5/6, that is repulsive.Proof.
It is straightforward that identitym~x,c=x leads us to the fixed points. Specifically,
(52)m~x,c=c2+c6x2−6x+1+x−12x222x−1c+x−1x=x⇔c−x21+c−4x+3x2=0,
provided that2x−1c+x−1x≠0. That is, the fixed points of m~x,c are the roots of px or the real roots of polynomial 1+c−4x+3x2. Then, real strange fixed points s1c and s2c appear when c≤1/3.Regarding their stability, it can be checked that the eigenvalues of the reduced Jacobian matrix coincide with the derivative ofm~xx,c: .
(53)EigM′z,zy,zz,x,czy=zz=xz=x=m~xx,c=x−12−cx2−c−2c+2x2−2x+121−2x2c+x−1x2,and then, their stability is easily stated.(54)m~xs1c,c=−24c−21−3c+7−6c+1−3c+1>1,∀c∈ℝ,and(55)m~xs2c,c=24c−21−3c−76c+1−3c−1<1,ifc∈725,13,where m~xs2c,c=0 for c=5/16. This is graphically checked in Figure 1.Figure 1
Absolute value of the eigenvaluesλi of Msic,sic,sic,sic,c,i=1,2.
(a)
λ1=m~s1c,c
(b)
λ2=m~s2c,cAn example of the behavior stated at Proposition1 is presented at Figure 2, where the basins of attraction for different values of c are showed. These pictures have been generated by means of a real mesh x,z of 800×800 points, following the routines appearing in [33]. For each initial pair x,z, intermediate values zy and zz have been generated from z by adding a small random value (in MATLAB code, zy=z+0.01∗rand1 and zz=z+0.01∗rand1). When the method is executed on each initial estimation x,z under these conditions, the point of the mesh is represented in orange or blue color if the process has converged to one of the roots c or −c at a distance lower than 10−3; other colors correspond to convergence to strange fixed points, and black points to initial guesses that make the method diverge or do not converge after a maximum number of iterations, 200.Figure 2
Dynamical planes ofOM10 on px=x2−c.
(a)
c=5/16
(b)
c=5/16
(c)
c=5/16
(d)
c=0.1
(e)
c=1
(f)
c=0Let us remark thatc=5/16 is the value of parameter c that makes the strange fixed point s2c superattractor; even in this situation, it is difficult to find initial estimations converging to it (a converging orbit to s25/16 is plotted in Figure 2(b)), as the most of them converge to the roots (orange and blue areas of Figure 2(a)) but some of them converge to 2-periodic orbits (black area appearing in Figure 2(c) where the initial estimation of the orbit is selected). Also, c=0.1 is used and convergence is observed (Figure 2(d)) to both roots with black areas of slow convergence. No other attractors appear. In Figure 2(e), the most frequent behavior is the convergence to the roots, but a black line of convergence to a 2-periodic orbit near 1,1 also appears. In the case of double root, corresponding to c=0 (see Figure 2(f)), convergence to x=z=0 is very slow.On the other hand, critical points play a key role in the stability of the system, as they appear always in any basin of attraction. Then, it is necessary, not only to calculate them but also to analyze their asymptotic behavior, in order to detect all the attracting elements: attracting fixed points, periodic orbits, strange attractors, and so on. This is the reason why we calculate the free critical points ofm~xx,c, with equal components, in the following result. Maybe other critical points of the multivariate rational functions exist, but their calculation is not feasible due to the high degree of the polynomials (of several variables) involved.Proposition 2.
The number of real critical points of m~x,c depends on parameter c, ±c, c>0 are always critical points,  and also,
(a)
If c>0, c1c=1−c and c2c=1+c are free critical points.
(b)
If c>1/4, there also exist two more free critical points: c3c=1/2−1/2−1+4c and c4c=1/2+1/2−1+4c.The proof of this result is direct from solvingm~xx,c=x−12−cx2−c−2c+2x2−2x+1/21−2x2c+x−1x2=0 or, equivalently, forcing the eigenvalues of the reduced Jacobian matrix to be zero.Let us remark that, although the dynamics of the method seems to be stable from the dynamical planes plotted and the analysis of the strange fixed points, the presence of any of these critical points can make them bifurcate into periodic orbits or to generate chaos for specific values ofc. In what follows, we analyze how it changes depending on c and if there exists any other kind of bifurcations that leads to repulsive or attractive points or to periodic or strange attractors.


### 3.2. Feigenbaum Diagrams

In order to analyze bifurcations, we employ Feigenbaum diagrams of the multidimensional rational function related to OM10 onpx, by using as initial estimation real critical points described in Proposition 2 and searching those ranges where changes of behavior happen. The critical points are not independent but conjugated two-by-two: c1c=−1/c2c and c3c=−1/c4c. Nevertheless, it is necessary to analyze the bifurcation diagram corresponding each one of them, as their performance is different.By using the critical pointc1c as the starting guess, a bifurcation diagram is observed in Figure 3(a). In this case, it is observed that the critical point is in the basin of attraction of the roots −c, as converges to it for c>0. A symmetric behavior is obtained for c2c as the initial estimation, belonging to the basin of attraction of the root c. In Figure 3(b), critical point c3c shows also convergence to one of the roots, but c4c shows a chaotic behavior for 1/4<c<1/3, with period-doubling bifurcation cascades as can be observed in Figures 3(c) and 3(d).Figure 3
Bifurcation diagrams of family OM10 onpx for free independent critical points.
(a)
c1c
(b)
c3c
(c)
c4c
(d)
c4cWe plot the iteration ofm~x,c for the values of c in one of the blue regions of Figure 3(d), in the x,z space. In the sequence shown in Figure 4, several strange attractors appear, separated or unified depending on very close values of c. These plots have been generated by fixing a value of c, with an initial estimation at x=0.8. OM10 which has been applied on it, plotting one point per iteration (blue for the first 2000 iterations, green for iteration 2001 to 4000, and magenta for iteration 4001 to 10,000). The resulting images show that several attracting points joint into wandering areas that are unified and separated depending on c.Figure 4
Bifurcation diagrams of family OM10 onpx for free independent critical points.
(a)
c=0.262
(b)
c=0.264
(c)
c=0.2644
(d)
c=0.2646
(e)
c=0.2648
(f)
c=0.265


## 3.1. Dynamical Analysis of OM10

Now, we study the behavior of the operator associated to scheme OM10 on quadratic polynomials. In order to get this aim, the asymptotic stability of the fixed points of the associate rational vector applied onpx=x2−c, with c∈ℝ, must be analyzed. It is clear that, if c>0, px has two simple real roots, it has one double root at x=0 if c=0 and px does not have real roots if c<0.We need to calculate the fixed points of the associate fixed point operator on method OM10 onpx=x2−c, Mz,zy,zz,x,c, that is, a high-degree rational function depending on the parameter of the polynomial c. However, it is not viable to analyze directly the fixed points of this operator, as usually indetermination appears when we calculate the fixed points due to cancelations of factors as z−x, zy−x, ansd zz−x in the denominator of Mz,zy,zz,x,c. In order to avoid this and using that all the fixed points have equal components, we force firstly zy=zz=x and, after simplifying, z=x. The resulting one-dimensional reduced operator is
(49)Mz,zy,zz,x,c∣zy=zz=x∣z=x=m~x,c=x2x−12+c6x2−6x+1+c222x−1c+xx−1.Moreover, there are two values ofc which reduce the rational function:
(i)
Ifc=0, then,
(50)m~x,0=−1+xx−2+4x.
(ii)
Ifc=1/4, then,
(51)m~x,14=5−4x+4x2−8+16x.Proposition 1.
The number of real fixed points (and their stability) of m~x,c depends on the value of parameter c:
(a)
If c<0, px has no real roots, so fixed points s1c=2/3−1/31−3c and s2c=2/3+1/31−3c are strange, being both repulsive.
(b)
If c=0, the only (double) root of px, x=0 is an attracting fixed point of m~x,0 and also x=1/3 is a repulsive strange fixed point.
(c)
When 0<c<1/4 or c>1/4, fixed points c,c and −c,−c are superattracting (corresponding to the roots of px). Moreover, s1c and s2c are strange fixed points if 1/4<c≤1/3; s1c is always repulsive, and s2c is attracting if c∈7/25,1/3 and superattracting if c=5/16.
For c=1/4, there is only one (superattracting) fixed point that is a root of px, x=−1/2 and one strange fixed point, x=5/6, that is repulsive.Proof.
It is straightforward that identitym~x,c=x leads us to the fixed points. Specifically,
(52)m~x,c=c2+c6x2−6x+1+x−12x222x−1c+x−1x=x⇔c−x21+c−4x+3x2=0,
provided that2x−1c+x−1x≠0. That is, the fixed points of m~x,c are the roots of px or the real roots of polynomial 1+c−4x+3x2. Then, real strange fixed points s1c and s2c appear when c≤1/3.Regarding their stability, it can be checked that the eigenvalues of the reduced Jacobian matrix coincide with the derivative ofm~xx,c: .
(53)EigM′z,zy,zz,x,czy=zz=xz=x=m~xx,c=x−12−cx2−c−2c+2x2−2x+121−2x2c+x−1x2,and then, their stability is easily stated.(54)m~xs1c,c=−24c−21−3c+7−6c+1−3c+1>1,∀c∈ℝ,and(55)m~xs2c,c=24c−21−3c−76c+1−3c−1<1,ifc∈725,13,where m~xs2c,c=0 for c=5/16. This is graphically checked in Figure 1.Figure 1
Absolute value of the eigenvaluesλi of Msic,sic,sic,sic,c,i=1,2.
(a)
λ1=m~s1c,c
(b)
λ2=m~s2c,cAn example of the behavior stated at Proposition1 is presented at Figure 2, where the basins of attraction for different values of c are showed. These pictures have been generated by means of a real mesh x,z of 800×800 points, following the routines appearing in [33]. For each initial pair x,z, intermediate values zy and zz have been generated from z by adding a small random value (in MATLAB code, zy=z+0.01∗rand1 and zz=z+0.01∗rand1). When the method is executed on each initial estimation x,z under these conditions, the point of the mesh is represented in orange or blue color if the process has converged to one of the roots c or −c at a distance lower than 10−3; other colors correspond to convergence to strange fixed points, and black points to initial guesses that make the method diverge or do not converge after a maximum number of iterations, 200.Figure 2
Dynamical planes ofOM10 on px=x2−c.
(a)
c=5/16
(b)
c=5/16
(c)
c=5/16
(d)
c=0.1
(e)
c=1
(f)
c=0Let us remark thatc=5/16 is the value of parameter c that makes the strange fixed point s2c superattractor; even in this situation, it is difficult to find initial estimations converging to it (a converging orbit to s25/16 is plotted in Figure 2(b)), as the most of them converge to the roots (orange and blue areas of Figure 2(a)) but some of them converge to 2-periodic orbits (black area appearing in Figure 2(c) where the initial estimation of the orbit is selected). Also, c=0.1 is used and convergence is observed (Figure 2(d)) to both roots with black areas of slow convergence. No other attractors appear. In Figure 2(e), the most frequent behavior is the convergence to the roots, but a black line of convergence to a 2-periodic orbit near 1,1 also appears. In the case of double root, corresponding to c=0 (see Figure 2(f)), convergence to x=z=0 is very slow.On the other hand, critical points play a key role in the stability of the system, as they appear always in any basin of attraction. Then, it is necessary, not only to calculate them but also to analyze their asymptotic behavior, in order to detect all the attracting elements: attracting fixed points, periodic orbits, strange attractors, and so on. This is the reason why we calculate the free critical points ofm~xx,c, with equal components, in the following result. Maybe other critical points of the multivariate rational functions exist, but their calculation is not feasible due to the high degree of the polynomials (of several variables) involved.Proposition 2.
The number of real critical points of m~x,c depends on parameter c, ±c, c>0 are always critical points,  and also,
(a)
If c>0, c1c=1−c and c2c=1+c are free critical points.
(b)
If c>1/4, there also exist two more free critical points: c3c=1/2−1/2−1+4c and c4c=1/2+1/2−1+4c.The proof of this result is direct from solvingm~xx,c=x−12−cx2−c−2c+2x2−2x+1/21−2x2c+x−1x2=0 or, equivalently, forcing the eigenvalues of the reduced Jacobian matrix to be zero.Let us remark that, although the dynamics of the method seems to be stable from the dynamical planes plotted and the analysis of the strange fixed points, the presence of any of these critical points can make them bifurcate into periodic orbits or to generate chaos for specific values ofc. In what follows, we analyze how it changes depending on c and if there exists any other kind of bifurcations that leads to repulsive or attractive points or to periodic or strange attractors.


## 3.2. Feigenbaum Diagrams

In order to analyze bifurcations, we employ Feigenbaum diagrams of the multidimensional rational function related to OM10 onpx, by using as initial estimation real critical points described in Proposition 2 and searching those ranges where changes of behavior happen. The critical points are not independent but conjugated two-by-two: c1c=−1/c2c and c3c=−1/c4c. Nevertheless, it is necessary to analyze the bifurcation diagram corresponding each one of them, as their performance is different.By using the critical pointc1c as the starting guess, a bifurcation diagram is observed in Figure 3(a). In this case, it is observed that the critical point is in the basin of attraction of the roots −c, as converges to it for c>0. A symmetric behavior is obtained for c2c as the initial estimation, belonging to the basin of attraction of the root c. In Figure 3(b), critical point c3c shows also convergence to one of the roots, but c4c shows a chaotic behavior for 1/4<c<1/3, with period-doubling bifurcation cascades as can be observed in Figures 3(c) and 3(d).Figure 3
Bifurcation diagrams of family OM10 onpx for free independent critical points.
(a)
c1c
(b)
c3c
(c)
c4c
(d)
c4cWe plot the iteration ofm~x,c for the values of c in one of the blue regions of Figure 3(d), in the x,z space. In the sequence shown in Figure 4, several strange attractors appear, separated or unified depending on very close values of c. These plots have been generated by fixing a value of c, with an initial estimation at x=0.8. OM10 which has been applied on it, plotting one point per iteration (blue for the first 2000 iterations, green for iteration 2001 to 4000, and magenta for iteration 4001 to 10,000). The resulting images show that several attracting points joint into wandering areas that are unified and separated depending on c.Figure 4
Bifurcation diagrams of family OM10 onpx for free independent critical points.
(a)
c=0.262
(b)
c=0.264
(c)
c=0.2644
(d)
c=0.2646
(e)
c=0.2648
(f)
c=0.265


## 4. Numerical Comparison

In this section, after a review of two and three steps with memory methods, we have tried some schemes to be compared with proposed biparametric three-step method OM10, (6). This method is compared with methods XW41 (16–18), XW42 (16–19), XW43 (16–20), XW44 (17-18), XW45 (17–19), and XW46 (17–20) presented in [5], XW5(18–24) given in [1], method SK7 in [29], and methods XW81 (37–34), XW82 (37–35), XW83 (37–36), XW84 (38–34), XW85 (38–35), and XW86 (38–36) given in [6] by using the test functions shown in Table 1.Table 1
Nonlinear test functions along with their zeros.
Nonlinear function
Root
f1=x3−x2−1
1.4655…
f2=x−13−1
2.0000…
f3=x3−10
2.1544…A numerical test of methods with memory is usually made by using starting pointx0 and also starting values of the accelerating parameters: γ0 to calculate y1 and z1 and thereafter λ0 in order to estimate x1. We check the performance of the methods by using this strategy, in order to see if the unstable behavior is avoided or, on the contrary, increased. For this, we use γ0=0.01 and λ0=2. It can be observed in Figure 5 that the stability of the method is highly improved when appropriate initial values are considered, not only x0 but also the steps y1 and z1 and, subsequently, the second iterate x1. The areas of convergence around the searched roots are wide in spite of the complexity of the nonlinear functions involved, and the behavior of the method is stable.Figure 5
Dynamical planes of OM10 on test functions.
(a)
f1x
(b)
f2x
(c)
f3xTable1 is furnished with the considered nonlinear test functions with their root (α). In the same table, there are infinite numbers of digits after the decimal but we have mentioned only four digits (the nonlinear functions are taken from [6, 7]). In Table 2, the absolute errors ∣xk−α∣ are given for presented method OM10. The computational order of convergence COC approximated using expression (see [34])
(56)COC≈ln∣fxn+1/fxn∣ln∣fxn/fxn−1∣,Table 2
Numerical comparison of bi-parameter with memory method.
Example
Method
Guess
∣x1−α∣
∣x2−α∣
∣x3−α∣
∣x4−α∣
COC
f1
XW41 (16–18),T0=1.5,a=8
1.3
NC
—
—
—
—
XW41 (16–19),T0=1.5,a=8
NC
—
—
—
—
XW41 (16–20),T0=1.5,a=8
NC
—
—
—
—
XW41 (17-18),T0=1.5,b=−2
0.0437e − 0
0.1582e − 6
0.7900e − 33
0.3499e − 148
4.5591
XW41 (17–19),T0=1.5,b=−2
0.0437e − 0
0.4206e − 8
0.1566e − 37
0.1579e − 154
4.3383
XW41 (17–20),T0=1.5,b=−2
0.0437e − 0
0.4206e − 8
0.3135e − 35
0.1015e − 152
4.3310
XW81 (37–34),T0=1.5,L=0
0.0240e − 0
0.2017e − 17
0.3801e − 153
0.1089e − 1374
9.0001
XW82 (37–35),T0=1.5,L=0
0.0240e − 0
0.5752e − 18
0.3290e − 174
0.1235e − 1736
10.000
XW83 (37–36),T0=1.5,L=0
0.0240e − 0
0.5752e − 18
0.3290e − 174
0.1250e − 1736
10.000
XW84 (38–34),T0=1.5,K=6
0.0828e − 0
0.2901e − 10
0.3076e − 80
0.1528e − 647
8.1216
XW85 (38–35),T0=1.5,K=6
0.0528e − 0
0.4136e − 13
0.1322e − 124
0.1472e − 1239
10.000
XW86 (38–36),T0=1.5,K=6
0.0528e − 0
0.4136e − 13
0.1322e − 124
0.1472e − 1239
10.000
XW5 (18–24),T0=−1.5,K0=−1.5
0.8364e − 3
0.1055e − 13
0.7715e − 66
0.1276e − 315
4.7910
SK7,α1=−1,α2=1
0.2002e − 1
0.4807e − 41
0.1535e − 289
0.3491e − 2029
7.0000
OM10 (4)–(8),γ=1.5,λ=1.2
0.7454e − 9
0.1755e − 72
0.9235e − 729
0.1496e − 7291
10.000
f2
XW41 (16–18),T0=1.5,a=8
2.5
NC
—
—
—
—
XW41 (16–19),T0=1.5,a=8
NC
—
—
—
—
XW41 (16–20),T0=1.5,a=8
NC
—
—
—
—
XW41 (17-18),T0=1.5,b=−2
NC
—
—
—
—
XW41 (17–19),T0=1.5,b=−2
NC
—
—
—
—
XW41 (17–20),T0=1.5,b=−2
NC
—
—
—
—
XW81 (37–34),T0=1.5,L=0
3.2783e − 4
5.6640e − 33
7.5622e − 293
9.9847e − 2632
9.0000
XW82 (37–35),T0=1.5,L=0
3.2783e − 3
1.7686e − 36
3.6959e − 354
5.8712e − 3586
10.000
XW83 (37–36),T0=1.5,L=0
3.2783e − 3
1.7686e − 36
3.6959e − 354
5.8712e − 3586
10.000
XW84 (38–34),T0=1.5,K=6
NC
—
—
—
—
XW85 (38–35),T0=1.5,K=6
NC
—
—
—
—
XW86 (38–36),T0=1.5,K=6
NC
—
—
—
—
XW5 (18–24),T0=−1.5,K0=−1.5
1.2283e − 2
4.9823e − 11
5.6852e − 53
1.0497e − 262
5.000
SK7,α1=−1,α2=1
1.0478e − 4
3.0825e − 29
5.8764e − 201
5.3772e − 1403
7.0000
OM10 (4)–(8),γ=1.5,λ=1.2
1.5515e − 3
3.9470e − 30
4.5314e − 296
1.8024e − 2955
10.000
f3
XW41 (16–18),T0=−1.5,a=8
2.0
2.5035e − 3
5.5559e − 15
6.3146e − 68
2.3358e − 309
4.8601
XW42 (16–19),T0=−1.5,a=8
2.5035e − 3
1.0056e − 15
8.3197e − 70
7.0470e − 305
4.3466
XW43 (16–20),T0=−1.5,a=8
2.5035e − 3
1.0056e − 15
1.6639e − 69
4.5101e − 303
4.3429
XW44 (17-18),T0=−1.5,b=−2
1.6966e − 3
1.0769e − 15
3.750e − 71
4.4929e − 324
4.5604
XW45 (17–19),T0=−1.5,b=−2
1.6966e − 3
1.4415e − 16
1.6141e − 73
4.4565e − 321
4.3469
XW46 (17–20),T0=−1.5,b=−2
1.6966e − 3
1.4415e − 16
3.2282e − 73
2.8521e − 319
4.3434
XW81 (37–34),T0=−1.5,L=0
1.4745e − 7
1.2463e − 66
1.9291e − 598
9.8391e − 5285
9.0000
XW82 (37–35),T0=−1.5,L=0
1.4745e − 7
5.9981e − 73
7.4408e − 727
6.4227e − 7266
10.000
XW83 (37–36),T0=−1.5,L=0
1.4745e − 7
5.9981e − 73
7.4408e − 727
6.4227e − 7266
10.000
XW84 (38–34),T0=−1.5,K=6
7.7322e − 5
5.9725e − 39
3.3296e − 315
2.3998e − 2359
8.1235
XW85 (38–35),T0=−1.5,K=6
1.5456e − 4
1.6975e − 41
4.3428e − 411
5.2174e − 4107
10.000
XW86 (38–36),T0=−1.5,K=6
1.5456e − 4
1.6975e − 41
4.3428e − 411
5.2174e − 4107
10.000
XW5 (18–24),T0=−1.5,K0=−1.5
4.6468e − 4
1.8636e − 19
1.9319e − 96
2.3132e − 481
5.000
SK7,α1=−1,α2=1
7.0030e − 9
1.8357e − 60
1.5607e − 421
5.0127e − 2949
7.0000
OM10 (4)–(8),γ0=−1.5, λ0=−1.5
3.0734e − 7
3.7135e − 70
2.4628e − 699
4.0535e − 6991
10.000to check the theoretical order. All the numerical results revealed in Table2 agree with previous results; for this, we have considered up to 2000 significant digits by using “Set Accuracy” command in Mathematica 8. OM10 has been used to solve the nonlinear functions, and the calculated results are compared with the other methods XW41 (16–18), XW42 (16–19), XW43 (16–20), XW44 (17-18), XW45 (17–19), XW46 (17–20),XW5(18–24), SK7, XW81 (37–34), XW82 (37–35), XW83 (37–36), XW84 (38–34), XW85 (38–35), and XW86 (38–36). In Table 2, “NC” means not convergent. From Table 2, it is very easy to identify that the results obtained by the proposed method are quite superior to the other two- and three-step schemes.


## 5. Conclusion

In this work, we have presented a new family of biparametric three-step schemes with memory to solve nonlinear equations. As a result, we have used two self-correcting parameters, designed by Hermite interpolating polynomials in the seventh-order method to achieve higher order convergence without any additional calculation. TheR-order of OM10 increases up to 10. The stability analysis of this family has been made by transforming it in a discrete dynamical system and studying the asymptotical behavior of the fixed and critical points. The method has been shown to be very stable except in few cases. These results have been checked by using dynamical planes, and the proposed method has been compared in performance and computational efficiency with a few existing methods by numerical examples. This confirms the validity of theoretical results.



---

*Source: 1232341-2018-06-28.xml*