Dynamics and Fractal Dimension of Steffensen-Type Methods

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-\frac{f(z_n)}{f^{\prime }(z_n)},n=0,1,2,\dots$$

(1)

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 = p^1/d, where p is the order of convergence and d is the number of functional evaluations per step. A method is optimal if p = 2^d−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-\frac{{[f(z_n)]}^2}{f(v_n)-f(z_n)}$$

(2)

where v_n = z_n + f(z_n). The composition technique (described, for instance, in [10]) allows up the implementation of high-order methods. From two methods with orders of convergence p₁ and p₂ is possible to perform a method of p₁ · p₂ 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

$$\begin{array}{l}{y}_n=z_n-\frac{f{(z_n)}^2}{f(v_n)-f(z_n)}\\ z_{n+1}=y_n-\frac{f(y_n)}{f^{\prime }(y_n)}\end{array}$$

(3)

where v_n = z_n + f(z_n). Now, in order to avoid the evaluation of f′(y_n), we approximate it by the derivative m′(y_n) of the following first degree Padé-like approximant

$$m(t)=\frac{a_1+a_2(t-y_n)}{1+a_3(t-y_n)}$$

where a₁, a₂ and a₃ are real parameters to be determined satisfying the following conditions:

$$m(z_n)=f(z_n)$$

(4)

$$m(y_n)=f(y_n)$$

(5)

and

$$m(v_n)=f(v_n)$$

(6)

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[z_n, y_n] denotes the divided difference $\frac{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]-\frac{f(v_n)f[z_n,y_n,v_n]}{f[v_n,z_n]}$$

and

$$a_3=-\frac{f[z_n,y_n,v_n]}{f[z_n,v_n]}$$

where $f[z_n,y_n,v_n]=\frac{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 y_n can be expressed as

$$m^{\prime }(y_n)=\frac{f[z_n,y_n]f[y_n,v_n]}{f[z_n,v_n]}$$

(7)

Substituting f′(y_n) from Equation (3) by its approximation m′(y_n), Cordero et al. in [11] obtained an optimal fourth-order iterative method denoted by M4, whose expression is:

$$\begin{array}{l}{y}_n=z_n-\frac{f{(z_n)}^2}{f(v_n)-f(z_n)}\\ x_{n+1}=y_n-\frac{f(y_n)f[z_n,v_n]}{f[z_n,y_n]f[y_n,v_n]}\end{array}$$

(8)

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

$$\{\begin{array}{l}{y}_n=z_n-\frac{{[f(z_n)]}^2}{f(v_n)-f(z_n)}\\ u_n=y_n-\frac{f(y_n)f[z_n,v_n]}{f[z_n,y_n]f[y_n,v_n]}\\ z_{n+1}=u_n-\frac{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]}\end{array}$$

(9)

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)=\frac{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

$$\{\begin{array}{l}{y}_n=z_n-\frac{{[f(z_n)]}^2}{f(v_n)-f(z_n)}\\ u_n=y_n-\frac{f(y_n)f[z_n,v_n]}{f[z_n,y_n]f[y_n,v_n]}\\ w_n=u_n-\frac{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-\frac{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)\}}\end{array}$$

(10)

where $c_5=-\frac{f[z_n,y_n,v_n,u_n,w_n]}{f[z_n,y_n,v_n,u_n]}$, c₄ = f [z_n; y_n; v_n;w_n] + c₅f [z_n; y_n; v_n], c₃ = f [z_n; y_n;w_n] − c₄(z_n + y_n − 2w_n) + c₅f [z_n; y_n]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 2^k−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 2^k−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 z₀ ∈ Ĉ is defined as {z₀, R(z₀), …, Rⁿ(z₀), …}.

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 z₀ ∈ C is a fixed point of R if R(z₀) = z₀. A fixed point is attracting, repelling or neutral if |R′(z₀)| is less than, greater than or equal to 1, respectively. Moreover, if |R′(z₀)| = 0, the fixed point is superattracting.

Let $z_f^{*}$ be an attracting fixed point of the rational function R. Its basin of attraction $\mathcal{A}(z_f^{*})$ is defined as the set of pre-images of any order such that $\mathcal{A}(z_f^{*})=\left\{z_0\in \widehat{\mathrm{C}}:R^n(z_0)\to z_f^{*},n\to \infty \right\}$. 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 z₀ in this area. If the sequence generated by iterative method reaches a zero z^* of the polynomial with a tolerance |z_k – z^*| < 10⁻³ and a maximum of 40 iterations, we decide that z₀ 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, $\mathcal{F}(R)$. The complementary set, the Julia set $\mathcal{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 O_f(z), it is possible to find its fixed and critical points. The fixed points z_f verify O_f(z) = z, and the critical points z_c validate $O_f^{\prime }(z)=0$.

The expressions

$$S_f(z)=z-\frac{f^2(z)}{f(v)-f(z)}$$

(11)

$$F_f(z)=y-\frac{f(y)f[z,v]}{f[z,y]-f[y,v]}$$

(12)

$$E_f(z)=u-\frac{f(u)f[z,y,v]}{f[y,v,u]f[u,z,y](u-z)+f[z,y,v]f[y,u]}$$

(13)

$$X_f(z)=w-\frac{f(w)}{f[z,w]+(z-w)\{c_{5}f[z,w]-c_3-c_4(z-w)\}}$$

(14)

where v = z +f(z), y = S_f(z), u = F_f(z) and w = E_f(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 → M_f(z) whose M_f 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 M_f. 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_pµ and M_p.

In order to study affine conjugacy classes of an iterative method M_f, 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: p_c (z) = z² + c and q_c (z) = z³ + c, where c ∈ {1, i}. The introduced fixed point operators satisfy the symmetry property $O_{f_{\overline{c}}}(\overline{z})=\overline{O_{f_c}(z)}$, ∀c, z ∈ C. Consequently, for polynomials with c ∈ R, the dynamical planes are symmetric respect to the abscissas axis. For polynomials with c ∈ C, the dynamical planes of $O_{f_{\overline{c}}}(z)$ are a vertical reflection of $O_{f_c}(z)$. Therefore, the study of p_c(z) and q_c(z), with c ∈ {1, i}, directly involves the knowledge of c ∈ {−1, −i}.

The dynamical planes of S_f(z), F_f(z), E_f(z) and X_f(z), when they are applied to p_{1,i}(z) and q_{1,i}(z), are displayed in Figures 1–4. 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 1–4 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, N_f(z), when it is applied to p_c(z) = z² + c and q_c(z) = z³ + 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=-\underset{ϵ\to 0}{\lim }\frac{\log (N(ϵ))}{\log (ϵ)}$$

(15)

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. Fractal dimension and percentage of Newton’s method vs. M4, M8 and M16.

Method		p1(z)	pi(z)	q1(z)	qi(z)
Ff(z)	Dimension	1.7146	1.6034	1.7039	1.6314
	Percentage	58.32%	66.73%	81.30%	84.91%
Ef(z)	Dimension	1.4477	1.1937	1.6616	1.5928
	Percentage	69.08%	89.64%	83.37%	86.97%
Xf(z)	Dimension	1.3610	1.3790	1.5808	1.5363
	Percentage	73.48%	77.59%	87.63%	90.17%
Nf(z)	Dimension	1.0000	1.0700	1.3853	1.3854
	Percentage	100.00%	100.00%	100.00%	100.00%

gathers the fractal dimension of operators (12), (13), (14) and N_f(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 N_qc(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.