Open access peer-reviewed chapter

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

By Alicia Cordero, Juan R. Torregrosa and Maria P. Vassileva

Submitted: October 22nd 2015Reviewed: May 5th 2016Published: October 19th 2016

DOI: 10.5772/64106

Downloaded: 1383


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.


  • 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→ ℝnis 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 Fon 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 innconvergent to ξ. Then, the convergence is said to be linear, if there exist M, 0 < M< 1, and k0 ∈ ℕ such thatx(k+1)ξMx(k)ξ,kk0, and of order p, p> 1, if there exist M, M> 0, and k0 ∈ ℕ such thatx(k+1)ξMx(k)ξp,kko.

On the other hand, Ostrowski [10] introduced the efficiency indexof an iterative method as p1/d, where pis the order of convergence and dis 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 opis 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-compositionthat 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 Fand 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: Ωnnbe a sufficiently Frèchet differentiable function in Ω. By using the notation introduced by Cordero et al. [26], the qth derivative of Fat un,q1, is the q-linear function F(q)(u):n×n×nnsuch that F(q)(u)(v1,,vq)n. It is easy to observe that: and F(q)(u)(v1,…,vq−1,·) ∈,

(ℝ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 ξ+hnlying 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 CqFqnsince F(q)(ξ) ∈

(ℝn× ⋯ × ℝn→ ℝn). In addition, we can express the Jacobian matrix of F, F′ as


where Iis the identity matrix. Therefore, qCqhq1

(ℝ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

X2=2C2,X3=4C223C3,X4=8C23+6C2C3+6C3C24C4, E5258525

We denote ek=x(k)ξthe error in the kth iteration. The equation ek+1=Lekp+O[ekp+1], where Lis a p-linear function L

(ℝn×…×ℝn, ℝn), is called error equation, and pis 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 pas 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 pand 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 ωiand τiare 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 jin 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=σjwith 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 2qdepends 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+ qunless σ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+ qdepending 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 pand eventually correcting it by the use of Gaussian quadrature, in case (a).

Theorem 2.1.LetF:Ωnnbe 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 byEq. (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 ismin{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×nthe Banach space of all n× nreal square matrices. The weight function in this context is a Frèchet differentiable function H:XXsatisfying:

  1. H(u)(v)=H1uv,Hbeing the first derivative of H,H:X

    (X)denotes the space of linear mapping from Xto itself.

  2. H(u,v)(w)=H2uvw,Hbeing the second derivative of H,H:X×X


Then, the Taylor expansion of Haround 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 functionF:Ω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 ykis 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.Letf:Ibe 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=11a1and for all a1 andb2with, a10then sequence{xk}k0obtained from (Eq. 5) converges to ξ with local order of convergence at least four. In this case, the error equation is



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 fremain 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)xkykis 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 Fon x(k) and y(k), Iis 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×n

(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 + hand 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 yis 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.Letf:Ωnnbe 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)}k0obtained by using expression (6) converges to ξ with local order of convergence at least four ifα=1,a2=a12(b22),b1=11a1and for alla1 andb2witha1 ≠ 0. The error equation is



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=C3C2X2X3and 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)C23and 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'=0We 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

, where
is the Riemann sphere, the orbit of a pointz0
is defined as {z0,R(z0),R2(z0),,Rn(z0),}. A point z*
is called a fixed pointof R(z)if it verifies that R(z*)=z*. Moreover, z*is called a periodic pointof 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-periodicif 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

Taking the associated multiplier into account, a fixed point z* is called: superattractingif
, attractingif
, repulsiveif
, and parabolicif
. 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 pointz*is a point satisfying |R'(z*)|=0.

The basin of attraction of an attractor α is defined as

The Fatou set of the rational function R,
is the set of points z
whose orbits tend to an attractor (fixed point, periodic orbit or infinity). Its complement in
is the Julia set,
. 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 = 1and


where A, B, C,and Ddepend 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 ofz= 1 as fixed point of operatorOp(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= −1and cri,i=1,2,3,4,are free critical points of Op(z,a1,b2) such that cr1=1cr2and cr3=1cr4, that will be used to calculate the different parameter planes associated with the fixed point operator. As z= −1is 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 a1and 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 operatorOp(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 operatorOp(z,a1,b2).

Selected values of parameters a1and b2in 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 ncan be classified depending on its asymptotic behavior. In this way, a point x*is a fixed point of Gif G(x*) = x*.

Theorem 3.2.[43, page 558] LetG:nnbe a function in

. Assume x*is a period-k point. Letλ1,λ2,,λnbe the eigenvalues of G'(x*).
  1. If all the eigenvaluesλjhave|λj|<1, then x*is attracting.

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

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

Moreover, a fixed point is called hyperbolicif all the eigenvalues λjof G′(x*)have |λj|1. Indeed, if there is an eigenvalue λisuch that |λi|<1and also there exists an eigenvalue λjsuch 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. if

    , for all{1,2,,n}, then x*is attracting.

  2. if

    , for alli,jϵ{1,2,,n}, then x*is superattracting.

  3. if

    , for alljϵ{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 functionsgi(x)satisfy

, for alli,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) = 0and q(x) = 0, where


The components of the iteration function Gare, 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.3838109or θ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 onp(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 onp(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 Aand vector bare different according to the method used, but in any case the inverse calculation −A−1bis computed by solving the linear system Ay= −b, by using Gaussian elimination with partial pivoting. Nevertheless, if several systems with the same matrix Amust be solved in iteration, LUfactorization is made once and the corresponding triangular systems are solved by substitution. The stopping criterion used is F(x(k))<10700or 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+1and 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,jdenotes the estimation of the unknown u(xi,yj), where xi= ihwith i= 1,2,…, nxand yj= jkwith j= 1,2, …, nyare the nodes in both variables, with h=1nx,k=1nyand 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)Tand 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.

© 2016 The Author(s). Licensee IntechOpen. This chapter is distributed under the terms of the Creative Commons Attribution 3.0 License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

How to cite and reference

Link to this chapter Copy to clipboard

Cite this chapter Copy to clipboard

Alicia Cordero, Juan R. Torregrosa and Maria P. Vassileva (October 19th 2016). Design, Analysis, and Applications of Iterative Methods for Solving Nonlinear Systems, Nonlinear Systems - Design, Analysis, Estimation and Control, Dongbin Lee, Tim Burg and Christos Volos, IntechOpen, DOI: 10.5772/64106. Available from:

chapter statistics

1383total chapter downloads

More statistics for editors and authors

Login to your personal dashboard for more detailed statistics on your publications.

Access personal reporting

Related Content

This Book

Next chapter

Nonlinear State and Parameter Estimation Using Iterated Sigma Point Kalman Filter: Comparative Studies

By Marwa Chaabane, Imen Baklouti, Majdi Mansouri, Nouha Jaoua, Hazem Nounou, Mohamed Nounou, Ahmed Ben Hamida and Marie‐France Destain

Related Book

First chapter

Recent Advances in Fragment Molecular Orbital-Based Molecular Dynamics (FMO-MD) Simulations

By Yuto Komeiji, Yuji Mochizuki, Tatsuya Nakano and Hirotoshi Mori

We are IntechOpen, the world's leading publisher of Open Access books. Built by scientists, for scientists. Our readership spans scientists, professors, researchers, librarians, and students, as well as business professionals. We share our knowledge and peer-reveiwed research papers with libraries, scientific and engineering societies, and also work with corporate R&D departments and government entities.

More About Us