Next Article in Journal
Model Equivalence-Based Identification Algorithm for Equation-Error Systems with Colored Noise
Previous Article in Journal
On String Matching with Mismatches
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Dynamics and Fractal Dimension of Steffensen-Type Methods

1
Institute of Telecommunications and Multimedia Applications (iTEAM), Universitat Politècnica de Valencia, Camino de Vera, s/n, 46022-Valencia, Spain
2
Institute of Multidisciplinary Mathematics, Universitat Politècnica de València, Camino de Vera, s/n,46022-Valencia, Spain
*
Author to whom correspondence should be addressed.
Algorithms 2015, 8(2), 271-279; https://0-doi-org.brum.beds.ac.uk/10.3390/a8020271
Submission received: 20 March 2015 / Accepted: 25 May 2015 / Published: 1 June 2015

Abstract

:
In this paper, the dynamical behavior of different optimal iterative schemes for solving nonlinear equations with increasing order, is studied. The tendency of the complexity of the Julia set is analyzed and referred to the fractal dimension. In fact, this fractal dimension can be shown to be a powerful tool to compare iterative schemes that estimate the solution of a nonlinear equation. Based on the box-counting algorithm, several iterative derivative-free methods of different convergence orders are compared.

1. Introduction

A large number of problems in science and engineering require the solution of a nonlinear equation f(z) = 0. There are different techniques to tackle this problem, and one of the most commonly used is the numerical solution with iterative methods. The Newton’s method is a well-known iterative scheme to estimate the solution of nonlinear equations
z n + 1 = z n f ( z n ) f ( z n ) , n = 0 , 1 , 2 ,
However, in many situations it is not possible to evaluate the derivative. In Section 2, several optimal derivative-free methods, with increasing orders of convergence, are introduced, in order to avoid the evaluation of the derivative.
The efficiency index, introduced by Ostrowski (see [1]), feeds back the comparison between methods in terms of efficiency. It is defined as I = p1/d, where p is the order of convergence and d is the number of functional evaluations per step. A method is optimal if p = 2d−1, as Kung and Traub conjectured in [2].
The quality of an iterative method can be measured by distinct parameters, as the efficiency index, the order of convergence, …, but the stability is a capital element that needs to be analyzed. The dynamical planes of the iterative schemes supply this information graphically, and are developed in Section 3. Several authors have studied and compared the stability of different known iterative methods by means of their dynamical planes, firstly in the work of Varona [3] and later on developed by Amat et al. [4], Neta et al. [5,6], among others.
The dynamical behavior of Newton’s method on polynomials, as studied in [7,8] among others, shows that its stability is not guaranteed in the whole dynamical plane. Finding the fractal dimension of an iterative method is a technique to measure this Julia set and, therefore, a useful tool to characterize the stability of a method. In Section 4, the fractal dimension of the introduced methods is evaluated.
In this paper, we analyze the dynamical behavior of four derivative-free schemes, with orders of convergence 2, 4, 8 and 16, on different quadratic and cubic polynomials. From this analysis, some results can be conjectured. For example, the Julia sets of the corresponding rational functions have less complexity and the basins of attraction obtained by the different schemes are wider, becoming more similar to that of Newton’s method, as the order of convergence increases. From a numerical point of view, these facts are checked by the fractal dimension of each procedure, that gets increasingly close to Newton’s one.

2. Optimal Derivative-Free Methods

The well-known Steffensen’s method (see [9]) replaces the derivative of Equation (1) by the forward difference, so its iterative expression is
z n + 1 = z n [ f ( z n ) ] 2 f ( v n ) f ( z n )
where vn = zn + f(zn). The composition technique (described, for instance, in [10]) allows up the implementation of high-order methods. From two methods with orders of convergence p1 and p2 is possible to perform a method of p1 · p2 order. In [11], the authors describe a technique based on composition and Padé-like approximants to obtain optimal derivative-free methods.
We first compose Steffensen’s method with Newton’s scheme obtaining the fourth-order scheme
y n = z n f ( z n ) 2 f ( v n ) f ( z n ) z n + 1 = y n f ( y n ) f ( y n )
where vn = zn + f(zn). Now, in order to avoid the evaluation of f′(yn), we approximate it by the derivative m′(yn) of the following first degree Padé-like approximant
m ( t ) = a 1 + a 2 ( t y n ) 1 + a 3 ( t y n )
where a1, a2 and a3 are real parameters to be determined satisfying the following conditions:
m ( z n ) = f ( z n )
m ( y n ) = f ( y n )
and
m ( v n ) = f ( v n )
From Equation (5), it is immediately obtained
a 1 = f ( y n )
Conditions (4) and (6), give, respectively,
a 2 f ( z n ) a 3 = f [ z n , y n ]
and
a 2 f ( v n ) a 3 = f [ v n , y n ]
where f[zn, yn] denotes the divided difference f ( z n ) f ( y n ) z n y n. After some algebraic manipulations, the following values are obtained for the parameters:
a 2 = f [ y n , v n ] f ( v n ) f [ z n , y n , v n ] f [ v n , z n ]
and
a 3 = f [ z n , y n , v n ] f [ z n , v n ]
where f [ z n , y n , v n ] = f [ z n , y n ] f [ y n , v n ] z n v n is a second order divided difference.
The derivative of the Padé approximant evaluated in yn can be expressed as
m ( y n ) = f [ z n , y n ] f [ y n , v n ] f [ z n , v n ]
Substituting f′(yn) from Equation (3) by its approximation m′(yn), Cordero et al. in [11] obtained an optimal fourth-order iterative method denoted by M4, whose expression is:
y n = z n f ( z n ) 2 f ( v n ) f ( z n ) x n + 1 = y n f ( y n ) f [ z n , v n ] f [ z n , y n ] f [ y n , v n ]
By composing again the resulting scheme with Newton’s method, and estimating the last derivative by a Padé-like approximant of order two, the obtained optimal eighth-order method, denoted by M8, has the iterative scheme
{ y n = z n [ f ( z n ) ] 2 f ( v n ) f ( z n ) u n = y n f ( y n ) f [ z n , v n ] f [ z n , y n ] f [ y n , v n ] z n + 1 = u n f ( u n ) f [ z n , y n , v n ] f [ y n , v n , u n ] f [ u n , z n , y n ] ( u n z n ) + f [ z n , y n , v n ] f [ y n , u n ]
In a similar way, we can obtain a new optimal 16th-order method, denoted by M16, by composing the scheme M8 with Newton’s method and estimating the last derivative by a Padé-like approximant of third degree
m ( t ) = c 1 + c 2 ( t z n ) + c 3 ( t z n ) 2 + c 4 ( t z n ) 3 1 + + c 5 ( t z n )
obtaining the iterative expression
{ y n = z n [ f ( z n ) ] 2 f ( v n ) f ( z n ) u n = y n f ( y n ) f [ z n , v n ] f [ z n , y n ] f [ y n , v n ] w n = u n f ( u n ) f [ z n , y n , v n ] f [ y n , v n , u n ] f [ u n , z n , y n ] ( u n z n ) + f [ z n , y n , v n ] f [ y n , u n ] z n + 1 = w n f ( w n ) f [ z n , w n ] + ( z n w n ) { c 5 f [ z n , w n ] c 3 c 4 ( z n w n ) }
where c 5 = f [ z n , y n , v n , u n , w n ] f [ z n , y n , v n , u n ], c4 = f [zn; yn; vn;wn] + c5f [zn; yn; vn], c3 = f [zn; yn;wn] − c4(zn + yn − 2wn) + c5f [zn; yn]and f [·, ·, ·, ·] and f [·, ·, ·, ·, ·] are the divided differences of order three and four, respectively.
This procedure can be extended in order to obtain optimal derivative-free iterative schemes with convergence order 2k−1, k = 2, 3, 4, 5, …. All the methods designed in this way are optimal in the sense of Kung-Traub’s conjecture, since the methods of order 2k−1 need k functional evaluation per iteration, k = 2, 3, …

3. Complex Dynamics of Iterative Methods

The dynamical planes provides an at-a-glance status of the convergence and stability of a method. In this section, the dynamical planes of the introduced methods are shown. In order to get a better understanding, some dynamical concepts are recalled. For a more extensive and comprehensive review of these concepts, see for example [7,8].
Let R : Ĉ Ĉ be a rational function, where Ĉ is the Riemann sphere. The orbit of a point z0 ∈ Ĉ is defined as {z0, R(z0), …, Rn(z0), …}.
The dynamical behavior of the orbit of a point on the complex plane can be classified depending on its asymptotic behavior. In this way, a point z0 ∈ C is a fixed point of R if R(z0) = z0. A fixed point is attracting, repelling or neutral if |R′(z0)| is less than, greater than or equal to 1, respectively. Moreover, if |R′(z0)| = 0, the fixed point is superattracting.
Let z f * be an attracting fixed point of the rational function R. Its basin of attraction A ( z f * ) is defined as the set of pre-images of any order such that A ( z f * ) = { z 0 C ^ : R n ( z 0 ) z f * , n }. In our calculations, we usually consider the region [5, 5] × [5, 5] of the complex plane, with 600 × 600 points and we apply the corresponding iterative method starting in every z0 in this area. If the sequence generated by iterative method reaches a zero z* of the polynomial with a tolerance |zk – z*| < 103 and a maximum of 40 iterations, we decide that z0 is in the basin of attraction of these zero and we paint this point in a color previously selected for this root. In the same basin of attraction, the number of iterations needed to achieve the solution is showed in darker or brighter colors (the less iterations, the brighter color). Black color denotes lack of convergence to any of the roots (with the maximum of iterations established) or convergence to the infinity.
The set of points whose orbits tends to an attracting fixed point (or an attracting periodic orbit) z f * is defined as the Fatou set, F ( R ). The complementary set, the Julia set J ( R ), is the closure of the set consisting of its repelling fixed points, and establishes the boundaries between the basins of attraction.
If the rational function R is associated to the fixed point operator of the developed methods in Section 2 applied on a polynomial f(z), denoted generically as Of(z), it is possible to find its fixed and critical points. The fixed points zf verify Of(z) = z, and the critical points zc validate O f ( z ) = 0.
The expressions
S f ( z ) = z f 2 ( z ) f ( v ) f ( z )
F f ( z ) = y f ( y ) f [ z , v ] f [ z , y ] f [ y , v ]
E f ( z ) = u f ( u ) f [ z , y , v ] f [ y , v , u ] f [ u , z , y ] ( u z ) + f [ z , y , v ] f [ y , u ]
X f ( z ) = w f ( w ) f [ z , w ] + ( z w ) { c 5 f [ z , w ] c 3 c 4 ( z w ) }
where v = z +f(z), y = Sf(z), u = Ff(z) and w = Ef(z), are the fixed point operators of Steffensen’s, M4, M8 and M16 methods, respectively.
From the dynamical point of view, conjugacy classes play an important role in the understanding of the behavior of classes of maps in the following sense. Let us consider a map z → Mf(z) whose Mf is any iterative root-finding map. Since a conjugacy preserves fixed and periodic points, as well as their basins of attraction, the dynamical information concerning f is carried by the character of the fixed points of Mf. So, for polynomials p of degree greater than or equal to k, it is interesting to build up parametrized families of polynomial pµ, as simple as possible, so that there exist a conjugacy between M and Mp.
In order to study affine conjugacy classes of an iterative method Mf, we need to establish a result that is is called Scaling Theorem. As the authors proved in [12], Steffensen’s method does not satisfy the Scaling Theorem, and this statement can be proved in a similar way for M4, M8 and M16. Then, these schemes have not conjugacy classes and we must study their dynamics on specific polynomials. The behavior of each method is analyzed on four different polynomials: pc (z) = z2 + c and qc (z) = z3 + c, where c ∈ {1, i}. The introduced fixed point operators satisfy the symmetry property O f c ¯ ( z ¯ ) = O f c ( z ) ¯, ∀c, z ∈ C. Consequently, for polynomials with cR, the dynamical planes are symmetric respect to the abscissas axis. For polynomials with cC, the dynamical planes of O f c ¯ ( z ) are a vertical reflection of O f c ( z ). Therefore, the study of pc(z) and qc(z), with c ∈ {1, i}, directly involves the knowledge of c ∈ {1, −i}.
The dynamical planes of Sf(z), Ff(z), Ef(z) and Xf(z), when they are applied to p{1,i}(z) and q{1,i}(z), are displayed in Figures 14. The basins of attraction are colored in orange and blue—for quadratic polynomials—and also in green—for cubic polynomials. The black points are those that do not converge to any of the attracting fixed points. We observe that the basins of attraction get wider when the order of convergence is increased, obtaining fast convergence in regions of initial guesses where the original scheme is divergent. A reason for this behavior can be found in [12], where the authors prove that the infinity is an attracting fixed point of Steffensen’s method, and it is repulsive in case of M4 (in a similar way, it can be checked that it is also repulsive for M8 and M16).
As Figures 14 show, except for Steffensen’s method, the complexity of the dynamical planes gets smoother as the order of convergence increases.
The dynamical planes of the fixed point operator associated to Newton’s method, Nf(z), when it is applied to pc(z) = z2 + c and qc(z) = z3 + c, with c ∈ {1, i}, are shown in Figure 5.
If we focus on the evolution of M4 to M16, passing through M8, it is observed that they are closer to the Newton’s dynamical planes, for each polynomial.

4. Fractal Dimension of Iterative Methods

The fractal dimension of the Julia set is a useful tool to analyze how close is a dynamical plane to another. Usually, the fractal dimension is obtained by the box-counting algorithm, based on the Hausdorff dimension. The foundations of this algorithm (see, for instance, [13]) goes by covering the Julia set by boxes small and smaller, in order to find
D = lim ϵ 0 log ( N ( ϵ ) ) log ( ϵ )
where ϵ is the length of the boxes, N(ϵ) is the account of boxes that cover the Julia set, and D is the fractal dimension. Computationally, the value of D is obtained as the slope of the regression line of log ϵ vs. log N(ϵ).
Table 1 gathers the fractal dimension of operators (12), (13), (14) and Nf(z) when they are applied to the previously submitted polynomials in the top rows. The comparison of each method to Newton’s one is made by the percentage of their fractal dimension, shown in bottom rows.
Graphically, we observe that the dynamical plane of the different methods on an specific polynomial “tends” to the corresponding dynamical plane of Newton’s one. We can see, for example, in Figure 4, how the different pictures are looking increasingly to Figure 5d. Also, the Julia sets have less complexity as the order of convergence increases. In a numerical sense, taking into account the fractal dimension of each method, it gets close and closer from Newton’s one.
There are some studies about the fractal dimension of Newton’s dynamical plane. For instance, in [14] the fractal dimension of Nqc(z), with c = 1, is D = 1.44692 in the square [2.5, 2.5] × [2.5, 2.5], whereas our calculus in this region gets D = 1.4055. Also, in [15], in the square [1, 1] × [1, 1], D ≈ 1.42, while in our study, the value of the fractal dimension in this square is D = 1.4242. The exact value depends on the details of the algorithm, such as the thickness of the Julia set or the chosen sequence of ϵ. However, our purpose is the comparison of the methods, so finding the fractal dimensions with the same algorithm is enough.

Acknowledgments

The authors thank to the anonymous referees for their useful comments and suggestions to improve the final version of the manuscript.
This research was partially supported by Ministerio de Economía y Competitividad MTM2014-52016-C02-02 and FONDOCYT 2014-1C1-088, República Dominicana.

Author Contributions

The contributions of all of the authors have been similar. All of them have worked together to develop the present manuscript.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Ostrowski, A.M. Solutions of Equations and Systems of Equations; Academic Press: New York, NY, USA, 1966. [Google Scholar]
  2. Kung, H.T.; Traub, J.F. Optimal order of one-point and multipoint iteration. J. Assoc. Comput. Mach. 1974, 21, 643–651. [Google Scholar]
  3. Varona, J.L. Graphic and numerical comparison bteween iterative methods. Math. Intell. 2002, 24, 37–46. [Google Scholar]
  4. Amat, S.; Busquier, S.; Bermúdez, C.; Plaza, S. On two families of high order Newton type methods. Appl. Math. Lett. 2012, 25, 2209–2217. [Google Scholar]
  5. Neta, B.; Scott, M.; Chun, C. Basins of attraction for several methods to find simple roots of nonlinear equations. Appl. Math. Comput. 2012, 218, 10548–10556. [Google Scholar]
  6. Chun, C.; Neta, B. An analysis of a new family of eighth-order optimal methods. Appl. Math. Comput. 2014, 245, 86–107. [Google Scholar] [Green Version]
  7. Devaney, R.L. An Introduction to Chaotic Dynamical Systems; Addison-Wesley Publishing Company: New York, NY, USA, 1989. [Google Scholar]
  8. Blanchard, P. Complex Analytic Dynamics on the Riemann Sphere. Bull. AMS. 1984, 11, 85–141. [Google Scholar]
  9. Cordero, A.; Torregrosa, J.R. A class of Steffensen type methods with optimal order of convergence. Appl. Math. Comput. 2011, 217, 7653–7659. [Google Scholar]
  10. Ortega, J.M.; Rheinboldt, W.G. Iterative Solutions of Nonlinear Equations in Several Variables; Academic Press: New York, NY, USA, 1970. [Google Scholar]
  11. Cordero, A.; Hueso, J.L.; Martínez, E.; Torregrosa, J.R. A new technique to obtain derivative-free optimal iterative methods for solving nonlinear equations. J. Comput. Appl. Math. 2012, 252, 95–102. [Google Scholar]
  12. Chicharro, F.; Cordero, A.; Gutiérrez, J.M.; Torregrosa, J.R. Complex dynamics of derivative-free methods for nonlinear equations. Appl. Math. Comput. 2013, 219, 7023–7035. [Google Scholar]
  13. Peitgen, H.O.; Jürgens, H.; Saupe, D. Chaos and Fractals: New Frontiers of Science, 2nd ed; Springer-Verlag: New York, NY, USA, 2004. [Google Scholar]
  14. Gutiérrez, J.M.; Magreñán, Á.A.; Varona, J.L. The “Gauss-Seidelization” of iterative methods for solving nonlinear equations in the complex plane. Appl. Math. Comput. 2011, 218, 2467–2479. [Google Scholar]
  15. Epureanu, B.I.; Greenside, H.S. Fractal basins of attraction associated with a damped Newton’s method. SIAM Rev. 1998, 40, 102–109. [Google Scholar]
Figure 1. Dynamical planes of iterative methods on p1(z) = z2 + 1. (a) Sp1(z); (b) Fp1(z); (c) Ep1(z); (d) Xp1(z).
Figure 1. Dynamical planes of iterative methods on p1(z) = z2 + 1. (a) Sp1(z); (b) Fp1(z); (c) Ep1(z); (d) Xp1(z).
Algorithms 08 00271f1
Figure 2. Dynamical planes of iterative methods on pi(z) = z2 + i. (a) Spi(z); (b) Fpi(z); (c) Epi(z); (d) Xpi(z).
Figure 2. Dynamical planes of iterative methods on pi(z) = z2 + i. (a) Spi(z); (b) Fpi(z); (c) Epi(z); (d) Xpi(z).
Algorithms 08 00271f2
Figure 3. Dynamical planes of iterative methods on q1(z) = z3 + 1. (a) Sq1(z); (b) Fq1(z); (c) Eq1(z); (d) Xq1(z).
Figure 3. Dynamical planes of iterative methods on q1(z) = z3 + 1. (a) Sq1(z); (b) Fq1(z); (c) Eq1(z); (d) Xq1(z).
Algorithms 08 00271f3
Figure 4. Dynamical planes of iterative methods on q1(z) = z3 + i. (a) Sq1(z); (b) Fq1(z); (c) Eq1(z); (d) Xq1(z).
Figure 4. Dynamical planes of iterative methods on q1(z) = z3 + i. (a) Sq1(z); (b) Fq1(z); (c) Eq1(z); (d) Xq1(z).
Algorithms 08 00271f4
Figure 5. Dynamical planes of Newton’s scheme on pc(z) = z2 + c, and qc(z) = z3 + c, where c ∈ {1; i}. (a) Np1(z); (b) Npi(z); (c) Nq1(z); (d) Nqi(z).
Figure 5. Dynamical planes of Newton’s scheme on pc(z) = z2 + c, and qc(z) = z3 + c, where c ∈ {1; i}. (a) Np1(z); (b) Npi(z); (c) Nq1(z); (d) Nqi(z).
Algorithms 08 00271f5
Table 1. Fractal dimension and percentage of Newton’s method vs. M4, M8 and M16.
Table 1. Fractal dimension and percentage of Newton’s method vs. M4, M8 and M16.
Methodp1(z)pi(z)q1(z)qi(z)
Ff(z)Dimension1.71461.60341.70391.6314
Percentage58.32%66.73%81.30%84.91%
Ef(z)Dimension1.44771.19371.66161.5928
Percentage69.08%89.64%83.37%86.97%
Xf(z)Dimension1.36101.37901.58081.5363
Percentage73.48%77.59%87.63%90.17%
Nf(z)Dimension1.00001.07001.38531.3854
Percentage100.00%100.00%100.00%100.00%

Share and Cite

MDPI and ACS Style

Chicharro, F.I.; Cordero, A.; Torregrosa, J.R. Dynamics and Fractal Dimension of Steffensen-Type Methods. Algorithms 2015, 8, 271-279. https://0-doi-org.brum.beds.ac.uk/10.3390/a8020271

AMA Style

Chicharro FI, Cordero A, Torregrosa JR. Dynamics and Fractal Dimension of Steffensen-Type Methods. Algorithms. 2015; 8(2):271-279. https://0-doi-org.brum.beds.ac.uk/10.3390/a8020271

Chicago/Turabian Style

Chicharro, Francisco I., Alicia Cordero, and Juan R. Torregrosa. 2015. "Dynamics and Fractal Dimension of Steffensen-Type Methods" Algorithms 8, no. 2: 271-279. https://0-doi-org.brum.beds.ac.uk/10.3390/a8020271

Article Metrics

Back to TopTop