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

Abstract

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.

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. [1–5]).

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 : Ω ⊆ ℝⁿ → ℝⁿ is a sufficiently differentiable function on the convex set Ω ⊆ ℝⁿ. 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. [6–8], 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

$x^{\left(k+1\right)}=x^{\left(k\right)}-{\left[F^{\text{'}}\left(x^{\left(k\right)}\right)\right]}^{-1}F\left(x^{\left(k\right)}\right),k=0,1,\dots$E1

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 in ℝⁿ convergent to ξ. Then, the convergence is said to be linear, if there exist M, 0 < M < 1, and k₀ ∈ ℕ such that $\Vert x^{(k+1)}-\xi \Vert \le M\Vert x^{(k)}-\xi \Vert ,\forall k\ge k_0$, and of order p, p > 1, if there exist M, M > 0, and k₀ ∈ ℕ such that $\Vert x^{(k+1)}-\xi \Vert \le M{\Vert x^{(k)}-\xi \Vert }^p,\forall k\ge k_o$.

On the other hand, Ostrowski [10] introduced the efficiency index of an iterative method as p^1/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 $p^{1/(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. [11–18]. 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 [21–23], 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 p₁ and p₂, respectively, are composed, the resulting scheme has order p₁p₂. 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

$z^{\left(k\right)}=y^{\left(k\right)}-5{\left[F^{\text{'}}\left(x^{\left(k\right)}\right)\right]}^{-1}F\left(y^{\left(k\right)}\right),x^{(k+1)}=z^{\left(k\right)}+\frac{1}{5}{\left[F^{\text{'}}\left(x^{\left(k\right)}\right)\right]}^{-1}\left[15F\left(y^{\left(k\right)}\right)-F(z^{\left(k\right)})\right],$

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,

$x^{(k+1)}=y^{(k)}+{\left[F^{\text{'}}\left(x^{\left(k\right)}\right)-5F^{\text{'}}\left(y^{\left(k\right)}\right)\right]}^{-1}\left[F^{\text{'}}\left(x^{\left(k\right)}\right)-5F^{\text{'}}\left(y^{\left(k\right)}\right)\right]{\left[F^{\text{'}}\left(x^{\left(k\right)}\right)\right]}^{-1}F\left(y^{\left(k\right)}\right),$

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: \Omega \subseteq {ℝ}^n\to {ℝ}^n$ be a sufficiently Frèchet differentiable function in Ω. By using the notation introduced by Cordero et al. [26], the qth derivative of F at $u\in {ℝ}^n,q\ge 1$, is the q-linear function $F^{\left(q\right)}\left(u\right):{ℝ}^n\times {ℝ}^n\cdots \times {ℝ}^n\to {ℝ}^n$ such that $F^{\left(q\right)}\left(u\right)\left(v_1,\dots ,v_q\right)\in {ℝ}^n$. It is easy to observe that: and F^(q)(u)(v₁,…,v_q−1,·) ∈,

(ℝⁿ) and $F^{\left(q\right)}\left(u\right)\left(v_{\sigma 1},\dots ,v_{\sigma q}\right)=F^{\left(q\right)}\left(u\right)\left(v_1,\dots ,v_q\right)$, for all permutation σ of {1, 2, …, q}. We will use the notation:

$F^{\left(q\right)}\left(u\right)\left(v_1,v_2,\dots ,v_q\right)=F^{\left(q\right)}\left(u\right)v_1{v}_2\dots v_q\text{and}{F}^{\left(q\right)}\left(u\right)v^{q-1}{F}^{\left(p\right)}{v}^{\left(p\right)}=F^{(q)}(u)F^{(p)}(u)v^{q+p-1}.$E96321

For $\xi +h\in {ℝ}^n$ 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

$F\left(\xi +h\right)=F^{\text{'}}\left(\xi \right)\left[h+{\displaystyle \sum }_{q=1}^{p-1}{C}_q{h}^q\right]+O\left[h^p\right],$

where $C_q=\left(\frac{1}{q!}\right){\left[F^{\text{'}}\left(\xi \right)\right]}^{-1}{F}^{(q)}\left(\xi \right),$q ≥ 2 We observe that $C_q{F}^q\in {ℝ}^n$ since F(q)(ξ) ∈

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

$F^{\text{'}}\left(\xi +h\right)=F^{\text{'}}\left(\xi \right)\left[I+{\displaystyle \sum }_{q=1}^{p-1}q{C}_q{h}^{q-1}\right]+O\left[h^p\right],$

where I is the identity matrix. Therefore, $q{C}_q{h}^{q-1}\in$

(ℝⁿ). From this expansion we can conjecture that

${\left[F^{\prime }(\xi +h)\right]}^{-1}=\left[I+X_{2}h+X_3{h}^2+X_4{h}^3+\cdots \right]{[F^{\prime }(\xi )]}^{-1}+O[h^p]$

and taking into account that ${\left[F\text{'}(\xi +h)\right]}^{-1}{F}^{\text{'}}\left(\xi +h\right)=F^{\text{'}}\left(\xi +h\right){\left[F\text{'}(\xi +h)\right]}^{-1}=I$, we obtain

$X_2=-2C_2,X_3=4C_2^2-3C_3,X_4=-8C_2^3+6C_2{C}_3+6C_3{C}_2-4C_4,\dots$E5258525

We denote $e_k=x^{(k)}-\xi$ the error in the kth iteration. The equation $e_{k+1}=L{e}_k^p+O\left[e_k^{p+1}\right]$, where L is a p-linear function $L\in$

(ℝⁿ×…×ℝⁿ, ℝⁿ), 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

$y^{(k)}=\xi +{\displaystyle \sum }_{j=q}^{5q-1}{M}_j{e}_k^j+O\left[e_k^{5q}\right],z^{(k)}=\xi +{\displaystyle \sum }_{j=p}^{5p-1}{N}_j{e}_k^j+O\left[e_k^{5p}\right],$E800

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:

$\begin{array}{l}(\text{a})x^{(k+1)}=y^{(k)}-2K^{-1}F(y^{(k)}),\text{(b)}{x}^{(k+1)}=z^{(k)}-2K^{-1}F(z^{(k)}),\\ \text{(c)}{x}^{(k+1)}=y^{(k)}-2K^{-1}F(z^{(k)}),\text{(d)}{x}^{(k+1)}=z^{(k)}-2K^{-1}F(y^{(k)}),\end{array}$E2

where for all cases $K={\displaystyle \sum {}_{i=1}^m}{\omega }_i{F}^{\prime }({\eta }_i^{(k)}),{\eta }_i^{(k)}=\frac{1}{2}[(1+{\tau }_i)z^{(k)}+(1-{\tau }_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, ${\eta }_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: ${\displaystyle \mathit{\sum }{}_{\mathit{j}\mathit{=}\mathit{q}}^{5\mathit{q}\mathit{-}\mathit{1}}}{\mathit{M}}_{\mathit{j}}{\mathit{e}}_{\mathit{k}}^{\mathit{j}}\mathit{=}{\mathit{A}}_{\mathit{1}(\mathit{q})}$ and ${\displaystyle \mathit{\sum }{}_{\mathit{j}\mathit{=}\mathit{q}}^{5\mathit{q}\mathit{-}\mathit{1}}}{\mathit{N}}_{\mathit{j}}{\mathit{e}}_{\mathit{k}}^{\mathit{j}}\mathit{=}{\mathit{A}}_{2(\mathit{p})}$, where the subscripts in parentheses denote the value of the smallest power assumed by j in the sum. By using this notation, ${\eta }_i^{(k)}$ can be expressed as
${\eta }_i^{(k)}=\xi +\frac{1}{2}{(R+{\tau }_{i}S)}_{(q)}$, where $R=A_{2(p)}+A_{1(q)}$ and $S=A_{2(p)}-A_{1(q)}$. By expansion of $F(y^{(k)}),F(z^{(k)})$ and $F\text{'}({\eta }_i^{\left(k\right)})$ in the Taylor series around ξ, we obtain

$F(y^{(k)})=F^{\prime }(\xi )[A_{1(q)}+C_2{A}_1^2{}_{(2q)}+C_3{A}_1^3{}_{(3q)}+C_4{A}_1^4{}_{(4q)}]+O[e_k^{5q}],F(z^{(k)})=F^{\prime }(\xi )[A_{1(p)}+C_2{A}_1^2{}_{(2p)}+C_3{A}_1^3{}_{(3p)}+C_4{A}_1^4{}_{(4p)}]+O[e_k^{5p}],F({\eta }_i^{(k)})=F^{\prime }(\xi )[I+B+C{\tau }_i+D{\tau }_i^2+E{\tau }_i^3]+O[e_k^{4q}],$E4563217

where$B=C_{1}R+\frac{3}{4}{C}_3{R}^2+\frac{1}{2}{C}_4{R}^3,C=C_{2}S+\frac{3}{4}{C}_3(RS+SR)+\frac{1}{2}{C}_4(R^{2}S+RSR+S{R}^2)$, $D=\frac{3}{4}{C}_3{S}^2+\frac{1}{2}{C}_4(SRS+S^{2}R+R^{2}S)$ and $E=\frac{1}{2}{C}_4{S}^3$. We also introduce the following notation: ${\displaystyle \sum {}_{i=1}^m}{\omega }_i=\sigma$ and $\frac{1}{\sigma }{\displaystyle \sum {}_{i=1}^m}{\omega }_i{\tau }_i^j={\sigma }_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={\displaystyle \sum {}_{i=1}^m}{\omega }_i{F}^{\prime }({\eta }_i^{(k)})$ appearing in Eq. (2) and we obtain $K=\sigma F^{\prime }\left(\xi \right)\left[I+K_{1(q)}+K_{2(2q)}+K_{3(3q)}\right]+O[e_k^{4q}]$,where$K_{1(q)}=C_2{(R+{\sigma }_{1}S)}_{(q)}$, $K_{2(2q)}=\frac{3}{4}{C}_3{(R^2+{\sigma }_1\left(RS+SR\right)+{\sigma }_2{S}^2)}_{(2q)}$and$K_{3(3q)}=[\frac{1}{2}{C}_4{(R^3+{\sigma }_1\left(R^{2}S+RSR+S{R}^2\right)+{\sigma }_2\left(R{S}^2+SRS+S^{2}R\right)+{\sigma }_2{S}^3)}_{(3q)}$. Recalling that $K^{-1}K=I$, we get $K^{-1}={\sigma }^{-1}\left[I+{K^{\prime }}_{1(q)}+{K^{\prime }}_{2(2q)}+{K^{\prime }}_{3(3q)}\right]{\left[F^{\prime }\left(\xi \right)\right]}^{-1}+O[e_k^{4q}]$, where ${K^{\prime }}_{1(q)}=-K_{1(q)},{K^{\prime }}_{2(2q)}={(K_1^2-K_2)}_{(2q)}$, and ${K^{\prime }}_{3(3q)}={\left(-K_1^3+K_2{K}_1+K_1{K}_2-K_3\right)}_{(3q)}$. Therefore, considering case (a), we obtain for $L=2K^{-1}F(y^{(k)})$ the following expression:

$L=\frac{2}{\sigma }{A}_{1(q)}+\frac{2}{\sigma }{\left[\left(C_2{A}_1+{K^{\prime }}_1\right)A_1\right]}_{(2q)}+\frac{2}{\sigma }{\left[\left(C_2{A}_1^2+{K^{\prime }}_1{C}_2{A}_1+{K^{\prime }}_2\right)A_1\right]}_{(3q)}+\frac{2}{\sigma }{\left[\left(C_4{A}_1^3+{K^{\prime }}_1{C}_3{A}_1^2+{K^{\prime }}_2{C}_2{A}_1+{K^{\prime }}_3\right)A_1\right]}_{(4q)}+O[e_k^{4q}].$E1478521

Since $x^{(k+1)}=y^{(k)}-L,$ the error equation can be expressed as

$\begin{array}{l}{e}_{k+1}=A_{1(q)}-\frac{2}{\sigma }{A}_{1\left(q\right)}-\frac{2}{\sigma }{\left[\left(C_2{A}_1+{K^{\prime }}_1\right)A_1\right]}_{\left(2q\right)}\\ -\frac{2}{\sigma }{\left[\left(C_2{A}_1^2+{K^{\prime }}_1{C}_2{A}_1+{K^{\prime }}_2\right)A_1\right]}_{(3q)}\end{array}-\frac{2}{\sigma }{\left[\left(C_4{A}_1^3+{K^{\prime }}_1{C}_3{A}_1^2+{K^{\prime }}_2{C}_2{A}_1+{K^{\prime }}_3\right)A_1\right]}_{(4q)}+O[e_k^{4q}]$

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 $(C_2{A}_1+{K^{\prime }}_1)A_1$. We develop it and get $\left(C_2{A}_1+{K^{\prime }}_1\right)A_1={\sigma }_1{C}_2{\left(A_1^2\right)}_{(2q)}-(1+{\sigma }_1)C_2{\left(A_2{A}_1\right)}_{(p+q)}$. Then, σ₁ = 0, the error equation is:

$e_{k+1}=C_2{\left(A_2{A}_1\right)}_{\left(p+q\right)}-\frac{2}{\sigma }{\left[\left(C_2{A}_1^2+{K^{\prime }}_1{C}_2{A}_1+{K^{\prime }}_2\right)A_1\right]}_{\left(3q\right)}-\frac{2}{\sigma }{\left[\left(C_4{A}_1^3+{K^{\prime }}_1{C}_3{A}_1^2+{K^{\prime }}_2{C}_2{A}_1+{K^{\prime }}_3\right)A_1\right]}_{\left(4q\right)}+O\left[e_k^{4q}\right].$E3

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 σ₁ = 2 or C₂ = 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 A₁ and A₂, which represent the errors of penultimate and last steps of the predictor method, and also of σ₁ 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:\text{Ω}\subset {ℝ}^n\to {ℝ}^n$ 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, σ₁ = 0 the order of convergence is min{p + q, 3q}.

In Table 1, we show the values of σ and σ₁ 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 nodes	Gaussian Quadrature
	Chebyshev		Legendre		Lobatto		Radau
	σ	σ1	σ	σ1	σ	σ1	σ	σ1
1	π	0	2	0	2	0	2	−1
2	π	0	2	0	2	0	2	0
3	π	0	2	0	2	0	2	0

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^{\prime }\left(\frac{y^{\left(k\right)}+z^{\left(k\right)}}{2}\right)$. In case of using Gauss-Radau quadrature with one node, we obtain Newton’s method $(K=F^{\prime }(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^{\prime }\left(y^{(k)}+z^{(k)}\right)$ and Noor’s scheme where $K=8\left[F^{\prime }\left(y^{\left(k\right)}\right)+3F^{\prime }\left(\frac{y^{\left(k\right)}+2z^{\left(k\right)}}{3}\right)\right]$, 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\times n}$ the Banach space of all n × n real square matrices. The weight function in this context is a Frèchet differentiable function $H:X\to X$ satisfying:

1. (i)

   $H^{\prime }(u)(v)=H_{1}uv,H^{\prime }$ being the first derivative of $H,H^{\prime }:X\to$

   
   $(X),H_1\in ℝ$ and $$
   
   $(X)$ denotes the space of linear mapping from X to itself.

2. (ii)

   $H^{″}(u,v)(w)=H_{2}uvw,H^{″}$ being the second derivative of $H,H^{″}:X\times X\to$

   
   $(X),H_2\in ℝ$.

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

$H\left({\psi }^{(k)}\right)\approx H_0\left(I\right)+H_1\left({\psi }^{(k)}-I\right)+\frac{1}{2}{H}_2{({\psi }^{(k)}-I)}^2$

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:

$z^{(k)}=x^{(k)}-\beta {[F^{\prime }(x^k)]}^{-1}F(x^k),x^{(k+1)}=z^{(k)}-2H\left(u^{\left(k\right)}\right){\left[{\displaystyle \sum }_{i=1}^n{\omega }_i{F}^{\text{'}}\left({\eta }_i^{\left(k\right)}\right)\right]}^{-1}F\left(x^{\left(k\right)}\right)$E4

where β is a non-zero real parameter and $u^{(k)}=\frac{1}{{{\displaystyle \sum }}_{i=1}^m{\omega }_i}{\left[F^{\prime }(x^{(k)})\right]}^{-1}{\displaystyle \sum {}_{i=1}^m}{\omega }_i{F}^{\prime }({\eta }_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:\Omega \subseteq {ℝ}^n\to {ℝ}^n$. Let us also suppose that the initial estimation x⁽⁰⁾ is close enough to the solution ξ and F′(ξ) is nonsingular. The iterative methods (Eq. 4) have order of convergence four if

$\beta =\frac{4(1+{\sigma }_1)}{3(1+2{\sigma }_1+{\sigma }_2)}$

and a sufficiently differentiable function H is chosen satisfying the conditions

$H_0\left(I\right)=\frac{\sigma }{2}I,H_1=\frac{\sigma (1+2{\sigma }_1+4{\sigma }_1^2+3{\sigma }_2)}{8{(1+{\sigma }_1)}^2}{H}_2=-\frac{3\sigma (-1-4{\sigma }_1-2{\sigma }_1^2+4{\sigma }_1^3-4{\sigma }_2-8{\sigma }_1{\sigma }_2+2{\sigma }_1^2{\sigma }_2-3{\sigma }_2^2)}{8{\left(1+{\sigma }_1\right)}^4}$

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 σ₁, σ₂, 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.

Quadrature	No. of nodes	σ	σ1	σ2	β	H0(I)	H1	H2
Gauss-Chebyshev	1	π	0	0	4/3	(π/2)I	π/8	3π/8
Gauss-Legendre	1	2	0	0	4/3	I	1/4	3/4
Gauss-Lobatto	1	2	0	0	4/3	I	1/4	3/4
	2	2	0	1	2/3	I	−1/2	6
Gauss-Radau	1	2	−1	0	0	(1/2)I	−	−
	2	2	0	1/3	1	I	0	2

Table 2.

Nodes, H₀(I), H₁ and H₂, 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.

Quadrature	No. of nodes	Name	H(t)
Gauss-Chebyshev	1	GC1	$\frac{\pi }{16}\frac{5-12t+15t^2}{t^2}$
Gauss-Legendre	1	GLe1	$\frac{1}{8}(9-4t+3t^2)$
Gauss-Lobatto	2	GLo2	$\frac{9}{2}-\frac{13}{2}t+3t^2$
Gauss-Radau	2	GR2	$t^2-2t+2$

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

$y^{\left(k\right)}=x^{\left(k\right)}-\frac{4}{3}{\left[F^{\text{'}}\left(x^{\left(k\right)}\right)\right]}^{-1}F\left(x^{\left(k\right)}\right),z^{\left(k\right)}=\frac{1}{2}\left(x^{\left(k\right)}+y^{\left(k\right)}\right),\begin{array}{l}{x}^{(k+1)}=x^{(k)}-\frac{9}{8}{\left[F^{\text{'}}\left(x^{\left(k\right)}\right)\right]}^{-1}F\left(x^{\left(k\right)}\right)\\ +\left(\frac{1}{2}I-\frac{3}{8}{\left[F^{\text{'}}\left(x^{\left(k\right)}\right)\right]}^{-1}{\left[F^{\text{'}}\left(z^{\left(k\right)}\right)\right]}^{-1}\right){\left[F^{\text{'}}\left(x^{\left(k\right)}\right)\right]}^{-1}F\left(x^{\left(k\right)}\right).\end{array}$

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 $x_{k+1}=y_k-\frac{f(x_k)}{f\left(x_k\right)-2f(y_k)}\frac{f(y_k)}{f^{\prime }(x_k)}$ and $x_{k+1}=y_k-\frac{f\left(x_k\right)+2f(y_k)}{f(x_k)}\frac{f(y_k)}{f^{\prime }(x_k)}$, respectively, where y_k is a Newton’s step. We propose a new family as a generalization of the previous methods in the following form:

$y_k=x_k-\alpha \frac{f(x_k)}{f^{\prime }(x_k)}{x}_{k+1}=y_k-\left[\frac{f\left(x_k\right)}{a_{1}f\left(x_k\right)+a_{2}f\left(y_k\right)}+\frac{b_{1}f\left(x_k\right)+b_{2}f\left(y_k\right)}{f\left(x_k\right)}\right]\frac{f\left(y_k\right)}{f^{\prime }\left(x_k\right)},$E5

where α, a₁, a₂, b₁, and b₂ 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\subseteq ℝ\to ℝ$ 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 $\alpha =1,a_2=a_1^2\left(b_2-2\right),b_1=1-\frac{1}{a_1}$ and for all a₁ and $b_2\in ℝ$ with, $a_1\ne 0$ then sequence ${\left\{x_k\right\}}_{k\ge 0}$ obtained from (Eq. 5) converges to ξ with local order of convergence at least four. In this case, the error equation is

$e_{k+1}=\left(\left(5-a_1{\left(b_2-2\right)}^2\right)c_2^3-c_2{c}_3\right)e_k^4+O\left[e_k^5\right],$E56324456

where $e_k=x_k-\xi$ and $c_q=\left(\frac{1}{q!}\right)\frac{f^{(q)}\left(\xi \right)}{f^{\text{'}}\left(\xi \right)},q\ge 2$.

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\left(x_k\right)=\frac{1}{\alpha }\left(x_k-y_k\right)f^{\prime }\left(x_k\right)$. By using this, we can rewrite quotient $\frac{f\left(y_k\right)}{f\left(x_k\right)}$ as

$\frac{f\left(y_k\right)}{f\left(x_k\right)}=1-\alpha \frac{f\left[x_k,y_k\right]}{f^{\prime }\left(x_k\right)},$E9000

where $f\left[x_k,y_k\right]=\frac{f\left(x_k\right)-f\left(y_k\right)}{x_k-y_k}$ 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

$y^{(k+1)}=x^{(k)}-\alpha {\left[F^{\text{'}}\left(x^{\left(k\right)}\right)\right]}^{-1}F\left(x^{\left(k\right)}\right),x^{(k+1)}=y^{(k)}-\left(G_1\left(x^{\left(k\right)},y^{\left(k\right)}\right)+G_2\left(x^{\left(k\right)},y^{\left(k\right)}\right)\right){\left[F^{\text{'}}\left(x^{\left(k\right)}\right)\right]}^{-1}F\left(y^{(k)}\right),G_1\left(x^{(k)},y^{(k)}\right)={\left[\left(a_1+a_2\right)I-\alpha a_2{\left[F^{\text{'}}\left(x^{\left(k\right)}\right)\right]}^{-1}\left[x^{(k)},y^{(k)};F\right]\right]}^{-1},G_2\left(x^{(k)},y^{(k)}\right)=\left(b_1+b_2\right)I-\alpha b_2{\left[F^{\text{'}}\left(x^{\left(k\right)}\right)\right]}^{-1}\left[x^{\left(k\right)},y^{\left(k\right)};F\right],$E6

where $\left[x^{\left(k\right)},y^{\left(k\right)};F\right]$ denotes the divided difference operator of F on x^(k) and y^(k), I is the identity matrix, and $F^{\prime }\left(x^{\left(k\right)}\right)$ 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^{\left(k\right)},y^{\left(k\right)};F]$ the divided difference operator defined by Ortega and Rheinboldt [37] as the function $\left[\cdot ,\cdot ;F\right]:\Omega \times \Omega \in {ℝ}^n\times {ℝ}^n\to$

$({ℝ}^n)$ that satisfies $\left[x,y;F\right]\left(x-y\right)=F\left(x\right)-F\left(y\right),\forall x,y\in \Omega$. To achieve this, we use the Genocchi-Hermite formula (see [37])

$\left[x,x+h;F\right]=\underset{0}{\overset{1}{{\displaystyle \int }}}{F}^{\prime }(x+th)dt$E963214

and, by developing $F^{\prime }(x+th)$ in Taylor series around x, we obtain

$\underset{0}{\overset{1}{{\displaystyle \int }}}{F}^{\prime }\left(x+th\right)dt=F^{\prime }\left(x\right)+\frac{1}{2}{F}^{″}\left(x\right)h+\frac{1}{6}{F}^{‴}\left(x\right)h^2+O\left[h^3\right].$E7

Assuming that

$F(x)=F^{\text{'}}(\xi )\left(e_k+C_2{e}_k^2+C_3{e}_k^3+C_4{e}_k^4\right)+O\left[e_k^5\right],F\text{'}(x)=F^{\text{'}}(\xi )\left(I+2C_2{e}_k+3C_3{e}_k^2+4C_4{e}_k^3\right)+O\left[e_k^4\right],F\text{'}\text{'}(x)=F^{\text{'}}(\xi )\left(2C_2+6C_3{e}_k+12C_4{e}_k^2\right)+O\left[e_k^3\right],F\text{'}\text{'}\text{'}(x)=F^{\text{'}}(\xi )\left(6C_3+24C_4{e}_k\right)+O\left[e_k^2\right],$E8

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 $e_{k,y}=y-\xi$, we have $\left[x,y;F\right]=F^{\text{'}}\left(\xi \right)\left[I+C_2(e_{k,y}-e_k)+C_3{e}_k^2\right]+O\left[e_k^3\right]$. In particular, if y is the approximation of the solution provided by Newton's method, i.e., $h=x-y={\left[F^{\text{'}}(x)\right]}^{-1}F(x)$, we obtain $\left[x,y;F\right]=F^{\text{'}}(\xi )\left[I+C_2{e}_k+\left(C_2^2+C_3\right)e_k^2\right]+O\left[e_k^3\right]$. These tools allow us to prove the following result.

Theorem 2.4. Let $f:\Omega \subseteq {ℝ}^n\to {ℝ}^n$ be a sufficiently differentiable function in a convex set Ω, and $\xi \in \Omega$ be a solution of the nonlinear system of equations F(x) = 0. Then, the sequence ${\left\{x^{(k)}\right\}}_{k\ge 0}$ obtained by using expression (6) converges to ξ with local order of convergence at least four if $\alpha =1,a_2=a_1^2(b_2-2),b_1=1-\frac{1}{a_1}$ and for all a₁ and $b_2\in ℝ$ with a₁ ≠ 0. The error equation is

$e_{k+1}=\left(\left(5-a_1{(b_2-2)}^2\right)C_2^3-C_2{C}_3\right)e_k^4+O\left[e_k^5\right],$E456321

where $e_k=x^{(k)}-\xi$ and $C_q=\left(\frac{1}{q!}\right){\left[F^{\text{'}}(\xi )\right]}^{-1}{F}^{(q)}(\xi ),q\ge 2$.

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

$F\left(x^{(k)}\right)=F^{\text{'}}(\xi )\left(e_k+C_2{e}_k^2+C_3{e}_k^3+C_4{e}_k^4\right)+O\left[e_k^5\right],F^{\text{'}\left(x^{\left(k\right)}\right)}=F^{\text{'}}(\xi )\left(I+2C_2{e}_k+3C_3{e}_k^2+4C_4{e}_k^3\right)+O\left[e_k^4\right].$

Let us consider ${\left[\text{F'}\left({\text{x}}^{(\text{k})}\right)\right]}^{-1}=\left(\text{I}+{\text{X}}_2{\text{e}}_{\text{k}}+{\text{X}}_3{\text{e}}_{\text{k}}^2+{\text{X}}_4{\text{e}}_{\text{k}}^3\right){\left[\text{F'(ξ)}\right]}^{-1}+\text{O}\left[{\text{e}}_{\text{k}}^4\right]$. Forcing ${\left[\text{F'}\left({\text{x}}^{(\text{k})}\right)\right]}^{-1}{\text{F}}^{\text{'}}\left({\text{x}}^{(\text{k})}\right)=\text{I}$, we get X₂ = −2C₂, X₃ = 2$C_2^2$ −3C₃ and ${\text{X}}_4=-4{\text{C}}_4+6{\text{C}}_3{\text{C}}_2-4{\text{C}}_2^2+6{\text{C}}_2{\text{C}}_3$. These expressions allow us to obtain for first step of iterative formula (6)

$y^{(k)}=\xi +(1-\alpha )e_k-\alpha \left(A_2{e}_k^2+A_3{e}_k^3+A_4{e}_k^4\right)+O\left[e_k^5\right],$E9

where $A_2=-C_2-X_2,A_3=-C_3-C_2{X}_2-X_3$ and $A_4=-C_4-C_3{X}_2-C_2{X}_3-X_4$. By using these results and the Taylor series expansion around ξ, we obtain

$F\left(y^{(k)}\right)=F^{\text{'}}(\xi )\left(B_1{e}_k+B_2{e}_k^2+B_3{e}_k^3+B_4{e}_k^4\right)+O\left[e_k^5\right],$

where B₁=β, B₂=(α + β²)C₂, B₃=−αA₃ + 2αβC₂A₃ + 3αβ²C₃C₂ + β⁴C₄, B₄=−αA₄ + α²C₂³−2αβC₂A₃ + 3αβ²C₃C₂ + β⁴C₄ and β = 1 − α. We calculate the Taylor expansion of [x^(k), y^(k);F] by using Eq. (9), [ x^(k), y^(k);F]= F^'(ξ)(I + D₂e_k + D₃e_k² + D₄e_k³) + O[e_k⁴], where D₂= (2−α)C₂, D₃= αC₂² + (3−3α + α²)C₃ and D₄= 2αC₂C₃ + α(3−2α)C₃C₂−(4−6α + 4α²−α³)C₄.

Then,

$\begin{array}{l}M=(a_1+a_2)I-\alpha a_2{\left[F^{\text{'}}\left(x^{(k)}\right)\right]}^{-1}\left[x^{(k)},y^{(k)};F\right]\\ =a_1+E_2{e}_k+E_3{e}_k^2+E_4{e}_k^3+O\left[e_k^4\right],\end{array}$

where E₂ = αa₂C₂, E₃ = αC₂³ + α(α−3)C₃ and E₄ = 6αC₂C₃−2αC₂³−4C₄ + 5α(2−α)C₃C₂ + (4−6α + 4α²−α³)C₃C₄.

Thus, we obtain G₁(x^(k), y^(k)) as the inverse of matrix M: G₁(x^(k), y^(k))=I + Y₂e_k+Y₃e_k²+Y₄e_k³+O[e_k⁴], where $Y_2=\frac{\alpha a_2}{a_1}{C}_2$, $Y_3=\frac{\alpha a_2}{a_1^2}\left[\left(\alpha a_2-3\right)C_2^2+\left(\alpha -3\right)C_3\right]$, $Y_4=\frac{\alpha a_2}{a_1^3}\left(8a_1+3\alpha a_1{a}_2+3\alpha a_2-{\alpha }^2{a}_2^3\right)C_2^3$ and G₂(x^(k), y^(k)) = b₁+F₂e_k+F₃e_k²+F₄e_k³+F₅e_k⁴+O[e_k⁵] where F₂=αb₂C₂, F₃=−αb₂[2C₂²−(α−3)C₃] and F₄=b₂[α(6−4α+α²)C₄−6α(2−α)C₃C₂+4(α+1)C₂³−6(α+1)C₂C₃].

Finally, we obtain the error equation of the proposed method

$e_{k+1}=H_1{e}_k+H_2{e}_k^2+H_3{e}_k^3+H_4{e}_k^4+O\left[e_k^5\right],$E1003

where $H_1=\frac{1}{a_1}\left(1+a_1\left(b_1-1\right)\right)(\alpha -1)$. If α = 1, then H₁ = 0 and the error equation takes the form: $e_{k+1}=H_2^{\text{'}}{e}_k^2+H_3^{\text{'}}{e}_k^3+H_4^{\text{'}}{e}_k^4+O\left[e_k^5\right]$ where $H_2^{\text{'}}=-\frac{1}{a_1}\left(1+a_1(b_1-1)\right)C_2$. We note that if $b_1=1-1/a_1$, then $H_2^{\text{'}}=0$ We introduce this value of b₁ and obtain the new form of the error equation $e_{k+1}=H_3^{\text{'}\text{'}}{e}_k^3+H_4^{\text{'}\text{'}}{e}_k^4+O\left[e_k^5\right]$ where $H_3^{\text{'}\text{'}}=\left(1/a_1^2\right)\left[a_2-a_1^2(b_2-2)\right]C_2^2$. Finally, if $a_2=a_1^2(b_2-2)$, the error equation is:

$e_{k+1}=-\left[\left(a_1{\left(b_2-2\right)}^2-5\right)C_2^3+C_2{C}_3\right]e_k^4+O\left[e_k^5\right]$

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:

$y^{(k+1)}=x^{(k)}-{\left[F^{\text{'}}\left(x^{\left(k\right)}\right)\right]}^{-1}F\left(x^{\left(k\right)}\right),x^{(k+1)}=y^{(k)}-\left(G_1\left(x^{\left(k\right)},y^{\left(k\right)}\right)+G_2\left(x^{\left(k\right)},y^{\left(k\right)}\right)\right){\left[F^{\text{'}}\left(x^{\left(k\right)}\right)\right]}^{-1}F\left(y^{(k)}\right),G_1\left(x^{(k)},y^{(k)}\right)=\frac{1}{a_1}{\left[\left(1+a_1{b}_2-2a_1\right)I-a_1\left(b_2-2\right){\left[F^{\text{'}}\left(x^{\left(k\right)}\right)\right]}^{-1}\left[x^{(k)},y^{(k)};F\right]\right]}^{-1}{G}_2\left(x^{(k)},y^{(k)}\right)=\frac{1}{a_1}\left[\left(a_1+a_1{b}_2-1\right)I-b_2{\left[F^{\text{'}}\left(x^{(k)}\right)\right]}^{-1}\left[x^{(k)},y^{(k)};F\right]\right].$E10

In the following we propose some particular cases:

1. When a₁ = 1, the iterative expression, being $G\left(x^{(k)},y^{(k)}\right)=G_1\left(x^{(k)},y^{(k)}\right)+G_2\left(x^{(k)},y^{(k)}\right)$ takes the form:

$\begin{array}{l}G\left(x^{(k)},y^{(k)}\right)={\left[\left(b_2-1\right)I-\left(b_2-2\right){\left[F^{\text{'}}\left(x^{(k)}\right)\right]}^{-1}\left[x^{(k)},y^{(k)};F\right]\right]}^{-1}\\ +b_{2}I-b_2{\left[F^{\text{'}}\left(x^{(k)}\right)\right]}^{-1}\left[x^{(k)},y^{(k)};F\right]\end{array}$

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

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

$x^{(k+1)}=y^{(k)}-\left(I-2{\left[F^{\text{'}\left(x^{(k)}\right)}\right]}^{-1}\left[x^{(k)},y^{(k)};F\right]\right){\left[F^{\text{'}\left(x^{(k)}\right)}\right]}^{-1}F\left(y^{(k)}\right)$E1000

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

$x^{(k+1)}=y^{(k)}-{\left(-I+2{\left[F^{\text{'}\left(x^{(k)}\right)}\right]}^{-1}\left[x^{(k)},y^{(k)};F\right]\right)}^{-1}{\left[F^{\text{'}\left(x^{(k)}\right)}\right]}^{-1}F\left(y^{(k)}\right).$E1001

2. When b₂ = 0 for any a₁ ≠ 0, the iterative expression of the parametric family is

$G\left(x^{(k)},y^{(k)}\right)=\frac{\left(3-2a_1\right)I+2\left(a_1-1\right){\left[F\text{'}\left(x^{(k)}\right)\right]}^{-1}\left[x^{(k)},y^{(k)};F\right]}{\left(1-2a_1\right)I+2a_1{\left[F\text{'}\left(x^{(k)}\right)\right]}^{-1}\left[x^{(k)},y^{(k)};F\right]}.$

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

$\begin{array}{l}{x}^{(k+1)}\\ =y^{(k)}\\ -\frac{\left(1+\beta \right)I+\beta {\left[F^{\text{'}\left(x^{(k)}\right)}\right]}^{-1}\left[x^{(k)},y^{(k)};F\right]}{\left(\beta -1\right)I+\left(\beta -2\right){\left[F^{\text{'}\left(x^{(k)}\right)}\right]}^{-1}\left[x^{(k)},y^{(k)};F\right]}{\left[F^{\text{'}\left(x^{(k)}\right)}\right]}^{-1}F\left(y^{(k)}\right).\end{array}$

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

$\begin{array}{l}{x}^{(k+1)}=y^{(k)}-(\frac{1}{a_1}{\left[\left(1-a_1\right)I-a_1{\left[F^{\text{'}}\left(x^{\left(k\right)}\right)\right]}^{-1}\left[x^{\left(k\right)},y^{\left(k\right)};F\right]\right]}^{-1}\\ +\frac{1}{a_1}[\left(2a_1-1\right)I\\ -{\left[F^{\text{'}}\left(x^{\left(k\right)}\right)\right]}^{-1}\left[x^{\left(k\right)},y^{\left(k\right)};F\right]]){\left[F^{\text{'}}\left(x^{\left(k\right)}\right)\right]}^{-1}F\left(y^{\left(k\right)}\right).\end{array}$

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. [38–45] 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 point $z_0\in$

is defined as $\left\{z_0,R(z_0),R^2(z_0),\dots ,R^n(z_0),\dots \right\}$. A point $z^{*}\in$

is called a fixed point of $R(z)$ if it verifies that $R\left(z^{*}\right)=z^{*}$. Moreover, z^* is called a periodic point of period p > 1 if it is a point such that $R^p\left(z^{*}\right)=z^{*}$ but $R^k\left(z^{*}\right)\ne z_0$, 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 $R^k\left(z^{*}\right)$ 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: superattracting if

, attracting if

, repulsive if

, and parabolic if

. 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 $\left|R^{\text{'}}\left(z^{*}\right)\right|=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\in$

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)=z^2-c$, 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

$\begin{array}{l}{O}_p\left(z,a_1,b_2\right)\\ =-z^4\frac{{\left(z+1\right)}^2\left(z^2+4z+5\right)-a_1\left(b_2^2-b_2\left(z^3+4z^2+5z+4\right)+2{\left(z+1\right)}^2\left(z+2\right)\right)}{z^4\left(a_1{\left(b_2-2\right)}^2-5\right)+z^3\left(-5a_1\left(b_2-2\right)-14\right)-2z^2\left(2a_1\left(b_2-2\right)+7\right)-z\left(a_1\left(b_2-2\right)+6\right)-1}\end{array}$E4568956

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 (a₁, b₂) in order to calculate the multipliers of the strange fixed points: z = 1 and

$e{x}_1=\frac{1}{12}\left(A-\sqrt{3B}-\sqrt{6(C+D)}\right),e{x}_2=\frac{1}{12}\left(A-\sqrt{3B}+\sqrt{6(C+D)}\right),e{x}_3=\frac{1}{12}\left(A+\sqrt{3B}-\sqrt{6(C+D)}\right),e{x}_4=\frac{1}{12}\left(A+\sqrt{3B}+\sqrt{6(C+D)}\right),$

where A, B, C, and D depend on both parameters a₁ and b₂. In Figure 1, the stability function of strange fixed point z = 1 is represented, $\left|O_p^{\text{'}}\left(1,a_1,b_2\right)\right|$. The orange region corresponds to real values of parameters a₁ and b₂ where z = 1 is attracting, meanwhile gray region includes values of a₁ and b₂ 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 (a₁, b₂).

Figure 1.

Stability function of z = 1 as fixed point of operator O_p(z, a₁, b₂).

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 $c{r}_i,i=1,2,3,4,$ are free critical points of O_p(z,a₁, b₂) such that $c{r}_1=\frac{1}{c{r}_2}$ and $c{r}_3=\frac{1}{c{r}_4}$, 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: cr₂ and cr₄.

The parameter space of a free critical point is obtained by associating each point of the parameter plane with real values of free parameters a₁ and b₂, 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 a₁ and b₂. 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⁻³ after 500 iterations, taking the same free independent critical point as the initial estimate.

Figure 2.

Parameter plane of operator O_p(z, a₁, b₂).

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 O_p(z, a₁, b₂).

Selected values of parameters a₁ and b₂ 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

$y^{(k)}=x^{(k)}-\theta {\left[F^{\text{'}}\left(x^{(k)}\right)\right]}^{-1}F\left(x^{(k)}\right),z^{(k)}=x^{(k)}-{\left[F^{\text{'}}\left(x^{(k)}\right)\right]}^{-1}\left(\theta F\left(x^{(k)}\right)+F\left(y^{(k)}\right)\right),x^{(k+1)}=x^{(k)}-{\left[F^{\text{'}}\left(x^{(k)}\right)\right]}^{-1}\left(\theta F\left(x^{(k)}\right)+F\left(y^{(k)}\right)+F\left(z^{(k)}\right)\right),$

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:{ℝ}^n\to {ℝ}^n$ be $afunctionin$

. Assume x^* is a period-k point. Let ${\lambda }_1,{\lambda }_2,\dots ,{\lambda }_n$ be the eigenvalues of G'(x^*).

1. a)

   If all the eigenvalues ${\lambda }_j$ have $\left|{\lambda }_j\right|<1$, then x^* is attracting.

2. b)

   If one eigenvalue ${\lambda }_j{}_{{}_0}$ has $\left|{\lambda }_{j_0}\right|>1$, then x^* is unstable, that is, repelling or saddle.

3. c)

   If all the eigenvalues ${\lambda }_j$ have $\left|{\lambda }_j\right|>1$, then x^* is repelling.

Moreover, a fixed point is called hyperbolic if all the eigenvalues ${\lambda }_j$ of G′(x^*) have $\left|{\lambda }_j\right|\ne 1$. Indeed, if there is an eigenvalue ${\lambda }_i$ such that $\left|{\lambda }_i\right|<1$ and also there exists an eigenvalue ${\lambda }_j$ such that $\left|{\lambda }_j\right|>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

   
   , for all $\left\{1,2,\dots ,n\right\}$, then x^* is attracting.

2. b)

   if

   
   , for all $i,jϵ\left\{1,2,\dots ,n\right\}$, then x^* is superattracting.

3. c)

   if

   
   , for all $jϵ\left\{1,2,\dots ,n\right\}$, then x^* is unstable and lies at the Julia set.

$g_i(x),i=1,2,\dots ,n$being 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 $g_i(x)$ satisfy

, for all $i,jϵ\left\{1,2,\dots ,n\right\}$.

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

$\begin{array}{cc}p\left(x_1,x_2\right)=\{\begin{array}{c}{x}_1^2-1,\\ x_2^2-1,\end{array}and& q\left(x_1,x_2\right)=\{\begin{array}{c}{x}_1^2-1,\\ x_2^2-1.\end{array}\end{array}$

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

$g_1\left(x_1,x_2\right)=-\frac{8{\theta }^2{x}_1^2{\left(x_1^2-1\right)}^3+{\theta }^4{\left(x_1^2-1\right)}^4-16x_1^4\left(3x_1^4+6x_1^2-1\right)}{128x_1^7},$

$g_2\left(x_1,x_2\right)=-\frac{8{\theta }^2{x}_2^2{\left(x_2^2-1\right)}^3+{\theta }^4{\left(x_2^2-1\right)}^4-16x_2^4\left(3x_2^4+6x_2^2-1\right)}{128x_2^7}.$

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 $\theta \approx -0.3838109$ or $\theta \approx 0.3838109$.

In Figure 4, we can see the dynamical plane of HMT method for $\theta \approx -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⁻³, 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]).

$p\approx \frac{\ln \left(\Vert x^{\left(k+1\right)}-x^{\left(k\right)}\Vert /\Vert x^{(k)}-x^{(k-1)}\Vert \right)}{\ln \left(\Vert x^{(k)}-x^{(k-1)}\Vert /\Vert x^{(k-1)}-x^{(k-2)}\Vert \right)}$

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⁻¹b. Matrix A and vector b are different according to the method used, but in any case the inverse calculation −A⁻¹b 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 $\Vert F\left(x^{(k)}\right)\Vert <{10}^{-700}$ or $\Vert x^{(k+1)}-x^{(k)}\Vert <{10}^{-700}$.

4.1. Molecular interaction problem

We are going to solve the equation of molecular interaction,

$u_{xx}+u_{yy}=u^2,(x,y)\in [0,1]\times [0,1],$

with the boundary conditions $u(x,0)=2x^2-x+1,u(x,1)=2,u(0,y)=2y^2-y+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,

$u_{i+1,j}-4u_{i,j}+u_{i-1,j}+u_{i,j+1}+u_{i,j-1}-h^2{u}_{i,j}^2=0,i=1,2,\dots ,nx,j=1,2,\dots ,ny,$

where u_i,j denotes the estimation of the unknown u(x_i,y_j), where x_i = ih with i = 1,2,…, nx and y_j = jk with j = 1,2, …, ny are the nodes in both variables, with $h=\frac{1}{nx},k=\frac{1}{ny}$ 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:

$x_1=u_{1,1},x_2=u_{2,1},x_3=u_{3,1},x_4=u_{1,2},x_5=u_{2,2},x_6=u_{3,2},x_7=u_{1,3},x_8=u_{2,3},x_9=u_{3,3}.$E4563

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

$F(x)=Ax+\phi (x)-b=0,$E4564

where

$A=\left(\begin{array}{ccc}M& -I& 0\\ -I& M& -I\\ 0& -I& M\end{array}\right),M=\left(\begin{array}{ccc}4& -1& 0\\ -1& 4& -1\\ 0& -1& 4\end{array}\right)$E4565

I being the 3 × 3 identity matrix, $\phi (x)=h^2{\left(x_1^2,x_2^2,\dots ,x_9^2\right)}^T$ and $b={\left(\frac{7}{4},1,\frac{27}{8},1,0,2,\frac{27}{8},2,4\right)}^T$.

In this case,

$F^{\text{'}}(x)=A+2h^{2}diag\left(x_1,x_2,\dots x_9\right).$E4566

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.

$MA:a_1=\frac{5}{4},b_2=0,MB:a_1=1,b_2=1,MC:a_1=1,b_2=3.$E4567

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.

Method	ACOC	Iteration	$\Vert x^{(k+1)}-x^{(k)}\Vert$	$\Vert F(x^{(k)})\Vert$
NM	1.9999	9	1.482e-413	6.448e-828
JM	3.9954	5	1.482e-413	1.976e-1007
OM	3.9964	5	1.482e-413	1.618e-1007
CM	3.9959	5	1.998e-353	1.618e-1007
MA	4.0519	5	5.362e-510	1.707e-2007
MB	3.9960	5	7.123e-362	1.409e-1449
MC	3.9960	5	3.110e-362	3.811e-1451

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.

Acknowledgments

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.