InTechOpen uses cookies to offer you the best online experience. By continuing to use our site, you agree to our Privacy Policy.

Mathematics » "Nonlinear Systems - Design, Analysis, Estimation and Control", book edited by Dongbin Lee, Tim Burg and Christos Volos, ISBN 978-953-51-2715-4, Print ISBN 978-953-51-2714-7, Published: October 19, 2016 under CC BY 3.0 license. © The Author(s).

Chapter 5

Design, Analysis, and Applications of Iterative Methods for Solving Nonlinear Systems

By Alicia Cordero, Juan R. Torregrosa and Maria P. Vassileva
DOI: 10.5772/64106

Article top


Stability function of z = 1 as fixed point of operator Op(z, a1, b2).
Figure 1. Stability function of z = 1 as fixed point of operator Op(z, a1, b2).
Parameter plane of operator Op(z, a1, b2).
Figure 2. Parameter plane of operator Op(z, a1, b2).
Dynamical plane of operator Op(z, a1, b2).
Figure 3. Dynamical plane of operator Op(z, a1, b2).
Dynamical plane of multivariate HMT method on p(x) for θ ≈ −0.3838109.
Figure 4. Dynamical plane of multivariate HMT method on p(x) for θ ≈ −0.3838109.
Different parameter plots of HMT family on p(x).
Figure 5. Different parameter plots of HMT family on p(x).
4-Periodic orbits in HMT family.
Figure 6. 4-Periodic orbits in HMT family.

Design, Analysis, and Applications of Iterative Methods for Solving Nonlinear Systems

Alicia Cordero1, Juan R. Torregrosa1 and Maria P. Vassileva2
Show details


In this chapter, we present an overview of some multipoint iterative methods for solving nonlinear systems obtained by using different techniques such as composition of known methods, weight function procedure, and pseudo-composition, etc. The dynamical study of these iterative schemes provides us valuable information about their stability and reliability. A numerical test on a specific problem coming from chemistry is performed to compare the described methods with classical ones and to confirm the theoretical results.

Keywords: system of nonlinear equations, iterative methods, order of convergence, weight function procedure, stability, basin of attraction

1. Introduction

The problem of solving equations and systems of nonlinear equations is among the most important in theory and practice, not only of applied mathematics, but also in many branches of science, engineering, physics, computer science, astronomy, finance, etc. A glance at the literature shows a high level of contemporary interest.

The search for solutions of systems of nonlinear equations is an old, frequent, and important problem for many applications in mathematics and engineering (for example, see Refs. [15]).

The main goal of this chapter is to describe different methods for approximating a solution ξ of a system of nonlinear equations F(x) = 0, where F : Ω ⊆ ℝn → ℝn is a sufficiently differentiable function on the convex set Ω ⊆ ℝn. The most commonly used techniques are iterative methods, which, from an initial guess a sequence of iterates is built, were converging to the solution of the problem under some conditions. Although not as many as in the case of scalar equations, some publications have appeared in the recent years, proposing different iterative methods for solving nonlinear systems (see, for example, Refs. [68], among others). They have made several modifications to the classical methods to accelerate the convergence and to reduce the number of operations and functional evaluations per step of the iterative method. Newton’s method is the most used iterative technique for solving these kind of problems (see Ref. [9]), whose iterative expression is


where F′ (x(k)) denotes the Jacobian matrix associated to function F on x(k).

We remember the concepts of order of convergence and efficiency index of an iterative scheme.

Definition 1.1. Let {x(k)}k≥0 be a sequence inn convergent to ξ. Then, the convergence is said to be linear, if there exist M, 0 < M < 1, and k0 ∈ ℕ such that x(k+1)ξMx(k)ξ,kk0, and of order p, p > 1, if there exist M, M > 0, and k0 ∈ ℕ such that x(k+1)ξMx(k)ξp,kko.

On the other hand, Ostrowski [10] introduced the efficiency index of an iterative method as p1/d, where p is the order of convergence and d is the number of functional evaluations per iteration. In the multidimensional case, it is more useful for the efficiency index to be defined as p1/(d+op), where op is the number of products-quotients per iteration.

The most direct technique is to adapt the methods designed for solving nonlinear equations to the multidimensional case. This process is easy only if, in the denominators of the iterative expression, does not appear any evaluation of the nonlinear function that describes the system. This is the case of Newton’s schemes and Newton-type methods coming from quadrature formulas, such as those described in Refs. [1118]. In Refs. [19] and [20], the authors designed a general procedure called pseudo-composition that allows to obtain predictor-corrector methods with high order of convergence. These multipoint schemes use any method as a predictor and a corrector step, where Gaussian quadrature is used.

On the other hand, other multipoint schemes have been developed by using different techniques: Adomian decomposition [2123], the replacement of the second derivative in one-point schemes by some approximation that yields multipoint iterative methods [8, 9], the Steffensen-type methods adapted to multidimensional case (see, for instance, Refs. [17, 24], among others). In all these papers, the references therein are also important.

A generally used technique for constructing iterative schemes is the composition of known methods. This technique was introduced by Traub [9]: if two iterative methods with orders of convergence p1 and p2, respectively, are composed, the resulting scheme has order p1p2. The use of this technique greatly increases the number of functional evaluations of F and F′. Therefore, it is necessary to avoid as much as possible such new evaluations by means of approximation techniques. In Ref. [25], the authors composed twice Newton’s scheme with itself ‘frozening’ the derivatives and, by means of undetermined coefficients method, obtained the following fifth-order iterative scheme


where y(k) is a Newton’s step. Let us observe that it only needs three functional evaluations and one Jacobian evaluation. Moreover, all the linear systems involved have the same matrix of coefficients. As a consequence, the efficiency index of this method is the best one, as far as we know. A similar procedure is used by Cordero et al. [26], getting Newton-Jarratt type methods of fifth and sixth orders. In addition, the following method described by Arroyo et al. [27] belongs to the class of Jarratt-type methods, but it has order of convergence five,


where y(k) is a Newton’s step.

In recent years, the technique of weight functions has also been developed, mainly for scalar equations. Weight functions are introduced in the iterative scheme to increase the order of convergence without increasing the number of functional evaluations. Among others, Sharma et al. [28] constructed a fourth-order scheme by using this procedure and, more recently, Artidiello et al. [7, 29] presented different families of high-order iterative methods by using matrix weight functions.

As we have previously mentioned, most of the iterative methods for nonlinear equations are not directly extendable to systems. However, in Refs. [6, 30], the authors present a general procedure to transform any scalar iterative method to the multidimensional case.

On the other hand, the dynamical analysis of an iterative method is becoming a trend in recent publications on iterative methods for scalar equations because it allows us to classify the different iterative formulas, not only from the point of view of its order of convergence, but also analyzing how these formulas behave as a function of the initial estimate that is taken. Another advantage of this analysis is to select the more stable elements of a parametric family whose members have the same order of convergence (see, for example, Ref. [31]). A first step in this direction on nonlinear systems was given by Cordero et al. [32], and a deeper analysis was made by Cordero et al. [33], studying the behavior of several methods on particular polynomial systems. In any case, the dynamical study in the multidimensional case is an emerging research topic with a promising future.

The rest of the chapter is organized as follows. In Section 2, we describe different techniques for designing iterative methods for nonlinear systems. In Section 3, we give some touches about the dynamic study of an iterative method for the scalar case and its extension to the multidimensional case. In order to check the introduced methods and compare them with other classical ones, in Section 4, we apply them on different test problems.

2. Design of the methods

In this section, we introduce three techniques for designing iterative methods for solving nonlinear systems of equations: pseudo-composition, weight function procedure, and a technique for extending scalar methods to the multidimensional case, in a non-trivial way.

The convergence results are going to be demonstrated by means of the n-dimensional Taylor expansion of the functions involved. Let F: Ωnn be a sufficiently Frèchet differentiable function in Ω. By using the notation introduced by Cordero et al. [26], the qth derivative of F at un,q1, is the q-linear function F(q)(u):n×n×nn such that F(q)(u)(v1,,vq)n. It is easy to observe that: and F(q)(u)(v1,…,vq−1,·) ∈,media/equl.png(ℝn) and F(q)(u)(vσ1,,vσq)=F(q)(u)(v1,,vq), for all permutation σ of {1, 2, …, q}. We will use the notation:


For ξ+hn lying in a neighborhood of a solution ξ of the nonlinear system F(x) = 0 and assuming that the Jacobian matrix F′(ξ) is nonsingular, Taylor’s expansion can be applied, obtaining


where Cq=(1q!)[F'(ξ)]1F(q)(ξ),q ≥ 2 We observe that CqFqn since F(q)(ξ) ∈media/equl.png(ℝn × ⋯ × ℝn → ℝn). In addition, we can express the Jacobian matrix of F, F′ as


where I is the identity matrix. Therefore, qCqhq1media/equl.png(ℝn). From this expansion we can conjecture that


and taking into account that [F'(ξ+h)]1F'(ξ+h)=F'(ξ+h)[F'(ξ+h)]1=I, we obtain


We denote ek=x(k)ξ the error in the kth iteration. The equation ek+1=Lekp+O[ekp+1], where L is a p-linear function Lmedia/equl.png(ℝn×…×ℝn, ℝn), is called error equation, and p is the order of convergence.

2.1. Pseudo-composition technique

We use the generic formulas of the Gaussian quadrature and develop families of predictor-corrector iterative methods, variants of Newton’s scheme, for solving nonlinear systems. Starting with any method of order p as a predictor and correcting over Gaussian quadrature, we will show that the final order of the obtained method will depend, among other things, on the order of the last two steps of the predictor. Let


be the penultimate and last steps of any iterative method with orders of convergence p and q, respectively. Taking this scheme as a predictor, we introduce the Gaussian quadrature as a corrector and we get four cases with the following iterative formulas:


where for all cases K=i=1mωiF(ηi(k)),ηi(k)=12[(1+τi)z(k)+(1τi)y(k)], and ωi and τi are the weights and nodes, respectively, of the orthogonal polynomial of degree m, which defines the corresponding Gaussian quadrature. Then, ηi(k) is calculated by using the points obtained in the last two steps of the predictor.

To simplify the calculations, we use the following notation: j=q5q1Mjekj=A1(q) and j=q5q1Njekj=A2(p), where the subscripts in parentheses denote the value of the smallest power assumed by j in the sum. By using this notation, ηi(k) can be expressed as
ηi(k)=ξ+12(R+τiS)(q), where R=A2(p)+A1(q) and S=A2(p)A1(q). By expansion of F(y(k)),F(z(k)) and F'(ηi(k)) in the Taylor series around ξ, we obtain

whereB=C1R+34C3R2+12C4R3,C=C2S+34C3(RS+SR)+12C4(R2S+RSR+SR2), D=34C3S2+12C4(SRS+S2R+R2S) and E=12C4S3. We also introduce the following notation: i=1mωi=σ and 1σi=1mωiτij=σj with j = 1, 2, …, which will allow us to simplify the analysis of the convergence conditions of the described methods.

Now, we develop the expression K=i=1mωiF(ηi(k)) appearing in Eq. (2) and we obtain K=σF(ξ)[I+K1(q)+K2(2q)+K3(3q)]+O[ek4q],whereK1(q)=C2(R+σ1S)(q), K2(2q)=34C3(R2+σ1(RS+SR)+σ2S2)(2q)andK3(3q)=[12C4(R3+σ1(R2S+RSR+SR2)+σ2(RS2+SRS+S2R)+σ2S3)(3q). Recalling that K1K=I, we get K1=σ1[I+K1(q)+K2(2q)+K3(3q)][F(ξ)]1+O[ek4q], where K1(q)=K1(q),K2(2q)=(K12K2)(2q), and K3(3q)=(K13+K2K1+K1K2K3)(3q). Therefore, considering case (a), we obtain for L=2K1F(y(k)) the following expression:


Since x(k+1)=y(k)L, the error equation can be expressed as


We note that if σ = 2 we obtain order of convergence at least 2q. The possibility of obtaining a convergence order greater than 2q depends on the expression (C2A1+K1)A1. We develop it and get (C2A1+K1)A1=σ1C2(A12)(2q)(1+σ1)C2(A2A1)(p+q). Then, σ1 = 0, the error equation is:


It is clear that the higher order of convergence in case (a) is min{p+q, 3q}. Finally, it is easy to prove that in case (b), the order of convergence of the resulting method is at least p. If σ = 2 the convergence order is p + q unless σ1 = 2 or C2 = 0, being in case the order of convergence is at least 2q + p. In cases (c) and (d) the convergence order is q, which is lower than the order of the predictor. It should be noticed that in case (b), two further functional evaluations are required and a new linear system per step must be solved. This causes the obtained method to be inefficient, from the standpoint of computational efficiency. In addition, we would like to remark that the order of convergence can be greater than p + q depending on expressions A1 and A2, which represent the errors of penultimate and last steps of the predictor method, and also of σ1 in case (b). The comments above allow us to state the following result, whose proof can be found by Cordero et al. [19], which establishes the order of convergence of the schemes that are obtained by using any method as a predictor of order p and eventually correcting it by the use of Gaussian quadrature, in case (a).

Theorem 2.1. Let F:Ωnn be sufficiently differentiable function in Ω and F′(x) continuous and nonsingular at ξΩ, solution of the nonlinear system. Let y(k) and z(k) be the penultimate and last steps of orders q and p, respectively, of a certain iterative method. Taking this scheme as a predictor we get a new approximation x(k + 1) of ξ given by Eq. (2). Then,

  1. the methods of the obtained families have an order of convergence at least q,

  2. if σ = 2 is satisfied, then the order of convergence is at least 2q,

  3. if, also, σ1 = 0 the order of convergence is min{p + q, 3q}.

In Table 1, we show the values of σ and σ1 for some Gaussian quadrature.

In terms of computational efficiency, the most efficient methods are those which use fewer nodes and few functional evaluations, so we only consider case (a). Also, Theorem 2.1 shows that the order of convergence does not depend on the number of nodes, it only depends on the order of convergence of the penultimate and last step of the predictor method. Therefore, it is computationally more efficient to use one or two nodes. We note that the Gauss-Chebyshev quadrature does not fulfill the second condition of Theorem 2.1. Then, its order of convergence is q. The method obtained by using Gauss-Radau quadrature of one node does not fulfill the third condition but it verifies the second one; hence, its order of convergence is at least 2q. The remaining quadrature with nodes 1, 2, and 3 satisfies the conditions of Theorem 2.1 and the order of convergence is at least min{p + q, 3q}.

Number of nodesGaussian Quadrature

Table 1.

Quadrature formula used.

If we use case (a) and the Gauss-Legendre quadrature with 1 node or Gauss-Lobatto quadrature with one node such as corrector, we obtain the midpoint method, where K=2F(y(k)+z(k)2). In case of using Gauss-Radau quadrature with one node, we obtain Newton’s method (K=F(y(k))). Finally, if we use Gauss-Lobatto quadrature with two nodes or Gauss-Radau quadrature with two nodes such as corrector, we obtain trapezoidal method with K=2F(y(k)+z(k)) and Noor’s scheme where K=8[F(y(k))+3F(y(k)+2z(k)3)], respectively.

2.2. Weight function procedure

The different methods obtained in the previous section are not optimal in the sense of Kung-Traub conjecture [34], when they are applied to scalar equations. By using the weight functions technique, we can increase the order of convergence of the designed methods without adding new functional evaluations.

We denote by X=n×n the Banach space of all n × n real square matrices. The weight function in this context is a Frèchet differentiable function H:XX satisfying:

  1. (i) H(u)(v)=H1uv,H being the first derivative of H,H:Xmedia/equl.png(X),H1 and media/equl.png(X) denotes the space of linear mapping from X to itself.

  2. (ii) H(u,v)(w)=H2uvw,H being the second derivative of H,H:X×Xmedia/equl.png(X),H2.

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


Now, by using a relaxed Newton’s method as a predictor and a weight function procedure in the corrector step, we design the following families of two-point schemes:


where β is a non-zero real parameter and u(k)=1i=1mωi[F(x(k))]1i=1mωiF(ηi(k)). Here, we have used the same notations previously introduced. The following result establishes the order of convergence of this iterative scheme under some conditions of function H.

Theorem 2.2. Let ξ ∈ Ω be a zero of a sufficiently differentiable function F:Ωnn. Let us also suppose that the initial estimation x(0) is close enough to the solution ξ and F′(ξ) is nonsingular. The iterative methods (Eq. 4) have order of convergence four if


and a sufficiently differentiable function H is chosen satisfying the conditions


and H‴ is a bounded operator, where σ and σj j = 1, 2, depend on the Gaussian quadrature used (see Ref. [35]).

Table 2 shows the values of σ1, σ2, and β, depending on the weights and the nodes of the orthogonal polynomials used in the Gaussian quadrature and also the corresponding values of the weight function and its derivatives, which are used in the iterative scheme.

QuadratureNo. of nodesσσ1σ2βH0(I)H1H2
Gauss-Chebyshev 1π004/3(π/2)Iπ/83π/8
Gauss-Legendre 12004/3I1/43/4
Gauss-Lobatto 12004/3I1/43/4
Gauss-Radau 12−100(1/2)I
2 2 0 1/3 1 I 0 2 

Table 2.

Nodes, H0(I), H1 and H2, for different Gaussian quadrature.

Let us remark that the class of methods designed allows the use of any Gaussian quadrature. In fact, when the technique of pseudo-composition was defined in Section 2.1 based also in Gaussian quadrature, the Chebyshev orthogonal polynomials could not be employed, as they did not verify the hypothesis of the main theorem. Nevertheless, with this new design, all the orthogonal polynomials can be applied and all of them derive optimal methods in the scalar case. In practice, we only use quadrature formulas with one and two nodes because, according to Theorem 2.2, the order of convergence is independent of the number of nodes.

For the one-dimensional case, according to the Kung-Traub’s conjecture [34], the obtained fourth-order methods are optimal (in case of Gauss-Radau with one node, classical Newton’s method is obtained). Other new methods are obtained in the rest of cases.

QuadratureNo. of nodesNameH(t)

Table 3.

Notation and weight functions for different Gaussian quadrature formulas.

In Table 3, we show the weight functions used for each iterative scheme coming from the respective Gaussian quadrature rules. Let us note that the iterative method coming from Gauss-Lobatto with one node is the same as the one resulting from the application of Gauss-Legendre, also with one node. In this case, both coincide with the fourth-order procedure recently published by Sharma et al. [28]. Let us note that other weight functions should derive in other new schemes. For example, the expression of the method obtained by using Gauss-Legendre quadrature with one node is


2.3. Divided difference operator

In this section, we present a design published by Cordero et al. [30] of some families of parametric iterative methods for solving nonlinear equations by means of some known schemes and afterwards extend one of them to systems of nonlinear equations. For this purpose, we use Ostrowski’s [10] and Chun’s [36] methods with iterative schemes xk+1=ykf(xk)f(xk)2f(yk)f(yk)f(xk) and xk+1=ykf(xk)+2f(yk)f(xk)f(yk)f(xk), respectively, where yk is a Newton’s step. We propose a new family as a generalization of the previous methods in the following form:


where α, a1, a2, b1, and b2 are real parameters. In the following result, we show which values of these parameters are necessary to guarantee at least order of convergence 4.

Theorem 2.3. Let f:I be a sufficiently differentiable function in an open interval I, such that ξ ∈ I is a simple root of the nonlinear equation f(x) = 0. If α=1,a2=a12(b22),b1=11a1 and for all a1 and b2 with, a10 then sequence {xk}k0 obtained from (Eq. 5) converges to ξ with local order of convergence at least four. In this case, the error equation is


where ek=xkξ and cq=(1q!)f(q)(ξ)f'(ξ),q2.

It is easy to prove this result by using any symbolic software as Wolfram Mathematica. To extend family (Eq. 5) to multivariate case, we need to rewrite its iterative expression in such a way that no functional evaluations of f remain at the denominator, as they will become vectors in the multidimensional case. So, let us consider that the first step of (Eq. 5) can be rewritten as f(xk)=1α(xkyk)f(xk). By using this, we can rewrite quotient f(yk)f(xk) as


where f[xk,yk]=f(xk)f(yk)xkyk is the first-order divided difference. By using this transformation, the proposed family (Eq. 5) is fully extensible to several variables in the following way


where [x(k),y(k);F] denotes the divided difference operator of F on x(k) and y(k), I is the identity matrix, and F(x(k)) is the Jacobian matrix of the system.

Since the analysis of the local convergence is based on Taylor series expansion around the solution, we need to obtain the corresponding development of the divided difference operator. Let us denote by [x(k),y(k);F] the divided difference operator defined by Ortega and Rheinboldt [37] as the function [,;F]:Ω×Ωn×nmedia/equl.png(n) that satisfies [x,y;F](xy)=F(x)F(y),x,yΩ. To achieve this, we use the Genocchi-Hermite formula (see [37])


and, by developing F(x+th) in Taylor series around x, we obtain


Assuming that


and replacing these developments in the formula of Genocchi-Hermite, denoting the second point of the divided difference by y = x + h and the error at the first step by ek,y=yξ, we have [x,y;F]=F'(ξ)[I+C2(ek,yek)+C3ek2]+O[ek3]. In particular, if y is the approximation of the solution provided by Newton's method, i.e., h=xy=[F'(x)]1F(x), we obtain [x,y;F]=F'(ξ)[I+C2ek+(C22+C3)ek2]+O[ek3]. These tools allow us to prove the following result.

Theorem 2.4. Let f:Ωnn be a sufficiently differentiable function in a convex set Ω, and ξΩ be a solution of the nonlinear system of equations F(x) = 0. Then, the sequence {x(k)}k0 obtained by using expression (6) converges to ξ with local order of convergence at least four if α=1,a2=a12(b22),b1=11a1 and for all a1 and b2 with a1 ≠ 0. The error equation is


where ek=x(k)ξ and Cq=(1q!)[F'(ξ)]1F(q)(ξ),q2.

Proof: By using Taylor expansion around ξ, we obtain:


Let us consider [F'(x(k))]1=(I+X2ek+X3ek2+X4ek3)[F'(ξ)]1+O[ek4]. Forcing [F'(x(k))]1F'(x(k))=I, we get X2 = −2C2, X3 = 2C22 −3C3 and X4=4C4+6C3C24C22+6C2C3. These expressions allow us to obtain for first step of iterative formula (6)


where A2=C2X2,A3=C3C2X2X3 and A4=C4C3X2C2X3X4. By using these results and the Taylor series expansion around ξ, we obtain


where B1=β, B2=(α + β2)C2, B3=−αA3 + 2αβC2A3 + 3αβ2C3C2 + β4C4, B4=−αA4 + α2C23−2αβC2A3 + 3αβ2C3C2 + β4C4 and β = 1 − α. We calculate the Taylor expansion of [x(k), y(k);F] by using Eq. (9), [ x(k), y(k);F]= F'(ξ)(I + D2ek + D3ek2 + D4ek3) + O[ek4], where D2= (2−α)C2, D3= αC22 + (3−3α + α2)C3 and D4= 2αC2C3 + α(3−2α)C3C2−(4−6α + 4α2α3)C4.



where E2 = αa2C2, E3 = αC23 + α(α−3)C3 and E4 = 6αC2C3−2αC23−4C4 + 5α(2−α)C3C2 + (4−6α + 4α2α3)C3C4.

Thus, we obtain G1(x(k), y(k)) as the inverse of matrix M: G1(x(k), y(k))=I + Y2ek+Y3ek2+Y4ek3+O[ek4], where Y2=αa2a1C2, Y3=αa2a12[(αa23)C22+(α3)C3], Y4=αa2a13(8a1+3αa1a2+3αa2α2a23)C23 and G2(x(k), y(k)) = b1+F2ek+F3ek2+F4ek3+F5ek4+O[ek5] where F2=αb2C2, F3=−αb2[2C22−(α−3)C3] and F4=b2[α(6−4α+α2)C4−6α(2−α)C3C2+4(α+1)C23−6(α+1)C2C3].

Finally, we obtain the error equation of the proposed method


where H1=1a1(1+a1(b11))(α1). If α = 1, then H1 = 0 and the error equation takes the form: ek+1=H2'ek2+H3'ek3+H4'ek4+O[ek5] where H2'=1a1(1+a1(b11))C2. We note that if b1=11/a1, then H2'=0 We introduce this value of b1 and obtain the new form of the error equation ek+1=H3''ek3+H4''ek4+O[ek5] where H3''=(1/a12)[a2a12(b22)]C22. Finally, if a2=a12(b22), the error equation is:


and the proof is finished.□

By using the same hypothesis as in Theorem 2.4, the iterative scheme of the class (Eq. 6) takes the form:


In the following we propose some particular cases:

1. When a1 = 1, the iterative expression, being G(x(k),y(k))=G1(x(k),y(k))+G2(x(k),y(k)) takes the form:


and we have a family of schemes with interesting particular cases, among others,

a) If b2 = 2, Chun's method transferred to systems is obtained


b) If b2 = 0 we get Ostrowski's scheme transferred to systems


2. When b2 = 0 for any a1 ≠ 0, the iterative expression of the parametric family is


If we express −2(a1 − 1) = β we get the King's family transferred to systems


3. When b2 = 1 for any a1 ≠ 0 and a1 ≠ 1 the iterative form, the last step of the method (Eq. 10) takes the form:


3. Dynamic studies of some methods

In the last years, a new branch of the analysis of iterative methods for solving nonlinear equations or systems has taken relevance: the dynamical analysis of the rational functions associated with the fixed point operator associated with the iterative scheme on polynomials. By using complex or real dynamics techniques, the stability and reliability of a method can be checked. Indeed, if a parametric family of iterative procedures is considered, this kind of analysis allows selecting those elements of the class with better stability and also to know which ones behave chaotically even on the most simple functions, as low degree polynomials.

The use of these dynamical tools is very frequent on scalar iterative methods; see, for example, Refs. [3845] and the references therein, but in the multidimensional case, it is a starting area of research.

Now, let us recall some basic concepts on complex dynamics. Given a rational function media/uueqn1.png, where media/uueqn2.png is the Riemann sphere, the orbit of a point z0media/uueqn2.png is defined as {z0,R(z0),R2(z0),,Rn(z0),}. A point z*media/uueqn2.png is called a fixed point of R(z) if it verifies that R(z*)=z*. Moreover, z* is called a periodic point of period p > 1 if it is a point such that Rp(z*)=z* but Rk(z*)z0, for each k < p. Moreover, a point z* is called pre-periodic if it is not periodic but there exists a k > 0 such that Rk(z*) is periodic.

There exist different types of fixed points depending on their associated multiplier media/uequ1R.png Taking the associated multiplier into account, a fixed point z* is called: superattracting if media/uueqn3.png, attracting if media/uueqn4.png, repulsive if media/uueqn5.png, and parabolic if media/uueqn6.png. The fixed point operator of any iterative method on an arbitrary polynomial p(z) is a rational function. The fixed points of this rational function that do not correspond to the roots of the polynomial p(z) are called strange fixed points. On the other hand, a critical point z* is a point satisfying |R'(z*)|=0.

The basin of attraction of an attractor α is defined as media/First equation.png The Fatou set of the rational function R, media/uueqn8.png is the set of points zmedia/uueqn2.png whose orbits tend to an attractor (fixed point, periodic orbit or infinity). Its complement in media/uueqn9.png is the Julia set, media/uueqn10.png. That means that the basin of attraction of any fixed point belongs to the Fatou set and the boundaries of these basins of attraction belong to the Julia set.

Some other concepts are a key fact in this kind of analysis, as the immediate basin of attraction of an attracting fixed point α (considered as a periodic point of period 1), as the connected component of the basin containing α. This concept is directly related with the existence of critical points, as can be seen in the following classical result, due to Fatou and Julia.

Theorem 3.1. Let R be a rational function. The immediate basins of attraction of any attracting periodic point hold, at least, a critical point.

If a scaling theorem can be established, the qualitative behavior of the class of iterative schemes is analyzed on a generic quadratic polynomial p(z)=z2c, its dynamics being analytically conjugated by affine transformations. Then, the fixed points of their associated rational function (being or not roots of the polynomial) are calculated and it is studied if they are stable (that is, if the successive iterations of the method converge to them) or if they are unstable, repelling the iterations near them. Finally, the calculation of critical points and their use as starting points of the iterative process will allow us, by applying Theorem 3.1, to find all the values of the parameter (that is, methods of the family) that do not converge to the roots of the polynomial, resulting in unstable behavior (attracting periodic orbits, chaos, etc.).

This kind of analysis has been developed by Cordero et al. [46] on the fourth-order biparametric Ostrowski-Chun family, whose iterative expression was shown in Eq. (5). In this manuscript, the strange fixed points and free independent critical points of the fixed point rational operator


have been identified. As the behavior of the elements of the family depends on two parameters, only real values have been considered in the plane (a1, b2) in order to calculate the multipliers of the strange fixed points: z = 1 and


where A, B, C, and D depend on both parameters a1 and b2. In Figure 1, the stability function of strange fixed point z = 1 is represented, |Op'(1,a1,b2)|. The orange region corresponds to real values of parameters a1 and b2 where z = 1 is attracting, meanwhile gray region includes values of a1 and b2 where this strange fixed point is repulsive, and therefore, the corresponding iterative methods have better behavior. Similar studies can be made on the rest of strange fixed points, searching for a more stable position of plane (a1, b2).


Figure 1.

Stability function of z = 1 as fixed point of operator Op(z, a1, b2).

On the other hand, as the dynamics of critical points could lead to a Fatou component, their associated parameter planes are represented to see the behavior of the method when the initial estimate is a critical point. In this case, z = −1 and cri,i=1,2,3,4, are free critical points of Op(z,a1, b2) such that cr1=1cr2 and cr3=1cr4, that will be used to calculate the different parameter planes associated with the fixed point operator. As z = −1 is a pre-image of z = 1, only two critical points can be considered as free and independent: cr2 and cr4.

The parameter space of a free critical point is obtained by associating each point of the parameter plane with real values of free parameters a1 and b2, so every point on the plane represents a different member of the iterative family. These parameter planes (one of them can be seen in Figure 2) have been created using a vectorized version of the MATLAB code programs presented by Chicharro et al. [47], with 800 × 800 different combinations of a1 and b2. Black points correspond to the parameter values for which the associated iterative method does not converge to the conjugated values of the zeros of polynomial p(z) with a tolerance of 10−3 after 500 iterations, taking the same free independent critical point as the initial estimate.


Figure 2.

Parameter plane of operator Op(z, a1, b2).

Points shown in red in Figure 2 (and also simultaneously red in the rest of parameter planes corresponding to other free independent critical points of the rational function) correspond to the most stable methods of the family. In these terms, Figure 3 shows the dynamical plane associated with one of these stable elements of the family: each point with 800 × 800 mesh in the complex plane is an initial estimation for the iterative method. This point has an associated color depending on the convergence of the method after 80 iterations: if it converges to a fixed point that corresponds to a root of p(z), then a color (orange and blue) is assigned (see Figure 3a). If not, then it diverges or converges to any other element (attracting strange fixed point, periodic orbit, etc.), and it is painted in black (see Figure 3b). In the former case, the orbit of an initial point in a black region is marked in yellow.


Figure 3.

Dynamical plane of operator Op(z, a1, b2).

Selected values of parameters a1 and b2 in these red regions were tested numerically, not only on scalar equations, but also on multidimensional problems. These tests showed that good performance of the member of the family observed in the dynamical study could also be noticed in the multidimensional case.

However, very recently real multivariate dynamical tools have revealed also to be a reliable implement to analyze the stability of iterative methods, specially designed for solving nonlinear systems. In order to study the stability of the fixed points of vectorial rational functions associated with iterative schemes for solving nonlinear systems, a new tool was presented by Cordero et al. [33]. In it, the authors checked the consistence of this tool by applying it on known methods as Newton's and Traub's schemes and also on a family of parametric iterative procedures


Presented by Hueso et al. [48] and denoted by HMT.

In order to analyze, under a multidimensional point of view, the qualitative behavior of this family, the Cordero et al. [33] introduced some concepts of multidimensional real discrete dynamics that we will recall in the following; some of them are the natural extensions of the defined concepts in complex dynamics, but others will be defined, as being different from those.

Let us denote by G(x) the vectorial fixed-point function associated with the iterative method or family on polynomial p(x). The dynamical behavior of the orbit of a point of n can be classified depending on its asymptotic behavior. In this way, a point x* is a fixed point of G if G(x*) = x*.

Theorem 3.2. [43, page 558] Let G:nn be a function in media/uequ2C.png. Assume x* is a period-k point. Let λ1,λ2,,λn be the eigenvalues of G'(x*).

  1. a) If all the eigenvalues λj have |λj|<1, then x* is attracting.

  2. b) If one eigenvalue λj0 has |λj0|>1, then x* is unstable, that is, repelling or saddle.

  3. c) If all the eigenvalues λj have |λj|>1, then x* is repelling.

Moreover, a fixed point is called hyperbolic if all the eigenvalues λj of G′(x*) have |λj|1. Indeed, if there is an eigenvalue λi such that |λi|<1 and also there exists an eigenvalue λj such that |λj|>1, the hyperbolic point is called saddle point.

To avoid the calculation of spectrum of G′(x*), the authors proposed by Cordero et al. [33] a result that, being consistent with the previous theorem, provided a practical tool to classify the stability of fixed points in many multidimensional cases.

Proposition 3.1. Let x* be a fixed point of G. then,

  1. a) if media/uueqn11.png, for all {1,2,,n}, then x* is attracting.

  2. b) if media/uueqn12.png, for all i,jϵ{1,2,,n}, then x* is superattracting.

  3. c) if media/uueqn13.png, for all jϵ{1,2,,n}, then x* is unstable and lies at the Julia set.

gi(x),i=1,2,,nbeing the coordinate functions of the fixed point multivariate function G.

Let us remark that if the order of convergence of the iterative method is at least two, then the roots of the nonlinear function are superattracting fixed points of the vectorial rational function corresponding to the iterative scheme. If a fixed point is not a root of the nonlinear function, it is called strange fixed point and its character can be analyzed in the same manner.

The concept of critical point can be defined following the idea of multivariate convergence of iterative methods.

Definition 3.1. A fixed point x* is a critical point of G if its coordinate functions gi(x) satisfy media/uueqn14.png, for all i,j ϵ{1,2,,n}.

From this definition, a superattracting fixed point will also be a critical point of the operator and, from numerical point of view, the iterative scheme will have, at least, quadratic order of convergence. A critical point that is not root of p(x) will be called free critical point.

The stability of this family is studied on two systems of real quadratic polynomials, p(x) = 0 and q(x) = 0, where


The components of the iteration function G are, in case of p(x),


It was proved by Cordero et al. [33] that the number of fixed points of the vectorial rational function G(x) associated with HMT iterative method on p(x) is 64, four of them (those corresponding to the roots of p(x)) are superattractive for any value of θ. Moreover, there are 48 unstable strange fixed points and the character of the final 12 real strange fixed points is simultaneously attractive for two ranges of values of θ, [−0.3847551, −0.3838109] and [0.3838109, 0.3847551], being superattractive if θ0.3838109 or θ0.3838109.

In Figure 4, we can see the dynamical plane of HMT method for θ0.3838109, where those 12 attracting strange fixed points appear as white circles, with small basins of attraction. The roots of p(x) appear as white stars, in their own basins of attraction.


Figure 4.

Dynamical plane of multivariate HMT method on p(x) for θ ≈ −0.3838109.

The parametric plots, as the extension of parameter planes in multivariable case, were revealed to be an interesting procedure that allowed to detect the most stable and unstable elements of a family of iterative methods. By using the information that gives us the iteration of the elements of the family on the different free critical points, we can know the global behavior of the family depending on the value of the parameter. Orbits of each free critical point are showed in Figure 5. In each one of these pictures, a different free critical point is used as initial guess of each member of the class of iterative schemes, taking values of parameter θ in [−5,5]. The values of θ corresponding to the members of the family are placed in the abscissa axis and the ordinate axis corresponds to 0.1, 0.2, 0.3, or 0.4 if the iterative process has converged to each one of the solutions of the quadratic polynomial real system, p(x), respectively. Moreover, the ordinate of a point is −0.1 if the process diverges and it is null in other cases.


Figure 5.

Different parameter plots of HMT family on p(x).


Figure 6.

4-Periodic orbits in HMT family.

The goal of these graphics is to show the elements of the family HMT (that is, the values of parameter θ) that present unstable behavior (attracting strange fixed points, periodic orbits, …). These plots were obtained by using 20,000 subintervals, a maximum of 40 iterations and an error estimation of 10−3, when iterates tend to a fixed point.

Thanks to this analysis, some non-desirable behavior was detected, as attracting periodic orbits of different periods, as can be observed in Figure 6, where two periodic orbits of period four are showed in yellow.

4. Numerical results

In this section, we are going to check the numerical performance of some elements of family (Eq. 10) compared with the multidimensional version of some other known methods of the same order and also with Newton's one.

Definition 4.1. Let ξ be a zero of function F and suppose that x(k−1), x(k) and x(k+1) are three consecutive iterations close to ξ. Then, the computational order of convergence p can be approximated using the formula (see Ref. [12]).


Numerical computations have been carried out in MATLAB, with variable precision arithmetic that uses floating point representation of 1000 decimal digits of mantissa. If n > 1, every iterate x(k+1) is obtained from the previous one, x(k), by adding one term of the form A−1b. Matrix A and vector b are different according to the method used, but in any case the inverse calculation −A−1b is computed by solving the linear system Ay = −b, by using Gaussian elimination with partial pivoting. Nevertheless, if several systems with the same matrix A must be solved in iteration, LU factorization is made once and the corresponding triangular systems are solved by substitution. The stopping criterion used is F(x(k))<10700 or x(k+1)x(k)<10700.

4.1. Molecular interaction problem

We are going to solve the equation of molecular interaction,


with the boundary conditions u(x,0)=2x2x+1,u(x,1)=2,u(0,y)=2y2y+1 and u(1, y) = 2, for all (x, y) in their domain.

For approximating its solution, we are going to use central divided differences procedure to transform the problem in a nonlinear system of equations. This system is solved by means of the proposed methods of order four and another known ones.

The discretization process yields the nonlinear system of equations,


where ui,j denotes the estimation of the unknown u(xi,yj), where xi = ih with i = 1,2,…, nx and yj = jk with j = 1,2, …, ny are the nodes in both variables, with h=1nx,k=1ny and nx = ny.

For checking purposes, we consider nx = ny = 4, getting a mesh of 5 × 5 points. Applying the boundary conditions we have only nine unknowns, which we rename as:


Then, the nonlinear system associated with the partial differential equation, including the boundary conditions, can be expressed as




I being the 3 × 3 identity matrix, φ(x)=h2(x12,x22,,x92)T and b=(74,1,278,1,0,2,278,2,4)T.

In this case,


For solving this system, we apply the extension for systems of Ostrowski's method (OM), Chun's scheme (CM), Jarratt's method (JM), Newton's method (NM), and three elements of family (Eq. 10) obtained by selecting stable values of the parameters.


In Table 4, the approximated computational order of convergence ACOC, the number of iterations, the difference between the two last iterations and the residual of the function at the last iteration are shown, for each one of the methods. In all cases, the initial estimation is a null vector and the Euclidean norm is used in the calculation of the residuals.


Table 4.

Numerical results for molecular interaction problem.

All the checked schemes provide the same solution of the nonlinear system. In Table 4, we can observe that all the fourth-order methods have a similar performance, but we note that the lowest error corresponds to method MA, duplicating the number of exact digits with respect to the other ones.

5. Conclusions

Many problems in science and engineering are modeled in such a way that, for their solution, it is necessary to solve systems of nonlinear equations. Therefore, designing iterative methods for solving these types of problems is an important task and it is a fruitful area of research. In this chapter, a review of the different techniques for constructing iterative methods is presented. Moreover, it is shown that real discrete dynamics tools are useful for analyzing the stability of the designed methods, selecting those with good dynamical behavior. In the numerical section, a chemical problem is used for testing the presented methods and the theoretical results are confirmed.


This research was partially supported by Ministerio de Economía y Competitividad of Spain, MTM2014-52016-C2-2-P, and by Ministry of Higher Education Science and Technology of Dominican Republic, FONDOCYT 2014-1C1-088.


1 - Bruns D.D., Bailey J.E. Nonlinear feedback control for operating a nonisothermal CSTR near and unstable steady state. Chemical Engineering Science. 1977;32:257–264.
2 - Ezquerro J.A., Gutiérrez J.M., Hernández M.A., Salanova M.A. Chebyshev-like methods and quadratic equations. Revue d’analyse Numérique et de Théorie de l’approximation. 2000;28:23–35.
3 - He Y., Ding C. Using accurate arithmetics to improve numerical reproducibility and stability in parallel applications. Journal of Supercomputing. 2001;18:259–277.
4 - Iliev A., Kyurkchev N. Nontrivial methods in numerical analysis: select topics in numerical analysis. Saarbcken, Germany: Lap Lambert Academic Publishing; 2010.
5 - Zhang Y., Huang P. High-precision time-interval measurement techniques and methods. Progress in Astronomy. 2006;24(1):1–15.
6 - Abad M., Cordero A., Torregrosa J.R. A family of seventh-order schemes for solving nonlinear systems. Bulletin Mathématique Societe des Sciences Mathematiques de Roumanie. 2014;57(105):133–145.
7 - Artidiello S., Cordero A., Torregrosa J.R., Vassileva M.P. Design of high-order iterative methods for nonlinear systems by using weight-function procedure. Abstract and Applied Analysis. 2015;2015:12. DOI: 10.1155/2015/289029
8 - Babajee D.K.R., Dauhoo M.Z., Darvishi M.T., Karami A., Barati A. Analysis if two Chebyshev-like three-order methods free from second derivative for solving systems of nonlinear equation. Journal of Computational and Applied Mathematics. 2010;233(8):2002–2012.
9 - Traub J.F. Iterative methods for the solution of equations. New York: Prentice Hall; 1964.
10 - Ostrowski A.M. Solutions of equations and systems of equations. Academic Press; 1966.
11 - Cordero A., Torregrosa J.R. Variants of Newton’s method for functions of several variables. Applied Mathematics and Computation. 2006;183:199–208.
12 - Cordero A., Torregrosa J.R. Variants of Newton’s method using fifth-order quadrature formulas. Applied Mathematics and Computation. 2007;190:686–698.
13 - Cordero A., Torregrosa J.R. On interpolation variants of Newton’s method for functions of several variables. Journal of Computational and Applied Mathematics. 2010;234:34–43.
14 - Frontini M., Sormani E. Third-order methods from quadrature formulae for solving systems of nonlinear equations. Applied Mathematics and Computation. 2004;149:771–782.
15 - Gerlach J. Accelerated convergence in Newton’s method. SIAM Review. 1994;36(2):272–276.
16 - Ozban A.Y. Some new variants of Newton’s method. Applied Mathematics Letters. 2014;17:677–682.
17 - Wang X., Zhang T., Qian W., Teng M. Seventh-order derivative-free iterative method for solving nonlinear systems. Numerical Algorithms. 2015;70:545–558.
18 - Weerakoon S., Fernando T.G.Y. A variant of Newton's method with accelerated third-order convergence. Applied Mathematics Letters. 2000;13:87–93.
19 - Cordero A., Torregrosa J.R., Vassileva M.P. Pseudo-composition: a technique to design predictor-corrector methods for systems of nonlinear equations. Applied Mathematics and Computation. 2012;218:11496–11508.
20 - Vassileva M.P. Métodos iterativos eficientes para la resolución de sistemas no lineales [thesis]. Valencia: Universitat Politècnica de València; 2011. 217 p. Available from:
21 - Adomian G. Solving frontier problems of Physics: the decomposition method. The Netherlands: Kluwer Academic ed; 1994.
22 - Babajee D.K.R., Dauhoo M.Z., Darvishi M.T., Barati A. A note on the local convergence of iterative methods based on Adomian decomposition method and three-node quadrature root. Applied Mathematics and Computation. 2008;200(1):452–458.
23 - Cordero A., Martínez E., Torregrosa J.R. Iterative methods of order four and five for systems of nonlinear equations. Journal of Computational and Applied Mathematics. 2009;231(2):541–551.
24 - Liu Z., Zheng Q., Zhao P. A variant of Steffensen’s method of fourth-order of convergence and its applications. Applied Mathematics and Computation. 2010;216:1978–1983.
25 - Arroyo V., Cordero A., Torregrosa J.R. Approximation of artificial satellites preliminary orbits: the efficiency challenge. Journal of Mathematical and Computer Modelling. 2011;54(7–8):1802–1807.
26 - Cordero A., Hueso J.L., Martínez E., Torregrosa J.R. A modified Newton-Jarratt’s composition. Numerical Algorithms. 2010;55:87–99.
27 - Arroyo V., Cordero A., Torregrosa J.R., Vassileva M.P. Artificial satellites preliminary orbit determination by modified high-order Gauss methods. International Journal of Computer Mathematics. 2012;89(3):347–356.
28 - Sharma J.R., Guha R.K., Sharma R. An efficient fourth-order weighted-Newton method for systems of nonlinear equations. Numerical Algorithms. 2013;62:307–323.
29 - Artidiello S., Cordero A., Torregrosa J.R., Vassileva M.P. Multidimensional generalization of iterative methods for solving nonlinear problems by means of weight-functions procedure. Applied Mathematics and Computation. 2015;268:1064–1071.
30 - Cordero A., García-Maimó J., Torregrosa J.R., Vassileva M.P. Solving nonlinear problems by Ostrowski-Chun type parametric families. Journal of Mathematical Chemistry. 2015;53:430–449.
31 - Cordero A., Feng L., Magreñán A.A., Torregrosa J.R. A new fourth-order family for solving nonlinear problems and its dynamics. Journal of Mathematical Chemistry. 2015;53:893–910.
32 - Cordero A., Torregrosa J.R., Vassileva M.P. Increasing the order of convergence of iterative schemes for solving nonlinear systems. Journal of Computational and Applied Mathematics. 2013;252:86–94.
33 - Cordero A., Soleymani F., Torregrosa J.R. Dynamical analysis of iterative methods for nonlinear systems or how to deal with the dimension? Applied Mathematics and Computation. 2014;244:398–412.
34 - Kung H.T., Traub F.F. Optimal order of one-point and multipoint iteration. Journal ACM. 1974;21:643–651.
35 - Cordero A., Torregrosa J.R., Vassileva M.P. Weighted-Gaussian correction of Newton-type methods for solving nonlinear systems. Bulletin Mathematique de la Societe des Sciences Mathematiques de Roumanie. 2016;59(107):23–38.
36 - Chun C. Construction of Newton-like iterative methods for solving nonlinear equations. Numerical Mathematics. 2006;104:297–315.
37 - Ortega J.M., Rheinboldt W.C. Iterative solutions of nonlinear equations in several variables. Academic Press Inc.; 1970.
38 - Amat S., Busquier S., Bermúdez C., Plaza S. On two families of high order Newton type methods. Applied Mathematic Letters. 2012;25:2209–2217.
39 - Amat S., Busquier S., Plaza S. Chaotic dynamics of a third-order Newton-type method. Journal of Computational and Applied Mathematics. 2010;366:24–32.
40 - Babajee D.K.R., Cordero A., Soleymani F., Torregrosa J.R. On improved three-step schemes with high efficiency index and their dynamics. Numerical Algorithms. 2014;65(1):153–169.
41 - Campos B., Cordero A., Magreñán A.A., Torregrosa J.R., Vindel P. Study of a biparametric family of iterative methods. Abstract and Applied Analysis. 2014;2014:12.
42 - Chun C., Neta B., Kozdon J., Scott M. Choosing weight functions in iterative methods for simple roots. Applied Mathematics and Computation. 2014;227:788–800.
43 - Cordero A., Maimó J.G., Torregrosa J.R., Vassileva M.P., Vindel P. Chaos in King’s iterative family. Applied Mathematics Letters. 2013;26:842–848.
44 - Cordero A., Torregrosa J.R., Vindel P. Dynamics of a family of Chebyshev-Halley type methods. Applied Mathematics and Computation. 2013;219:8568–8583.
45 - Magreñán A.A. Different anomalies in a Jarratt family of iterative root-finding methods. Applied Mathematics and Computation. 2014;233:29–38.
46 - Cordero A., Maimó J.G., Torregrosa J.R., Vassileva M.P. Stability of a fourth order bi-parametric family of iterative methods. Journal of Computational and Applied Mathematics. DOI: 10.1016/
47 - Chicharro F.I., Cordero A., Torregrosa J.R. Drawing dynamical and parameter planes of iterative families and methods. The Scientific World Journal. 2013;2013:11.
48 - Hueso J.L., Martínes E., Torregrosa J.R. New modifications of Potra-Pták’s method with optimal fourth and eighth orders of convergence. Journal of Computational and Applied Mathematics. 2010;234:2969–2976.
49 - Robinson R.C. An introduction to dynamical systems, continuous and discrete. Providence: American Mathematical Society; 2012.