New Predictor-Corrector Methods with High Efficiency for Solving Nonlinear Systems

Abstract

A new set of predictor-corrector iterative methods with increasing order of convergence is proposed in order to estimate the solution of nonlinear systems. Our aim is to achieve high order of convergence with few Jacobian and/or functional evaluations. Moreover, we pay special attention to the number of linear systems to be solved in the process, with different matrices of coefficients. On the other hand, by applying the pseudocomposition technique on each proposed scheme we get to increase their order of convergence, obtaining new efficient high-order methods. We use the classical efficiency index to compare the obtained procedures and make some numerical test, that allow us to confirm the theoretical results.

1. Introduction

Many relationships in nature are inherently nonlinear, which according to these effects are not in direct proportion to their cause. Approximating a solution 𝜉 of a nonlinear system, 𝐹(𝑥)=0, is a classical problem that appears in different branches of science and engineering (see, e.g. [1]). In particular, the numerical solution of nonlinear equations and systems is needed in the study of dynamical models of chemical reactors [2] or in radioactive transfer [3]. Moreover, many of numerical applications use high precision in their computations; in [4], high-precision calculations are used to solve interpolation problems in astronomy; in [5] the authors describe the use of arbitrary precision computations to improve the results obtained in climate simulations; the results of these numerical experiments show that the high-order methods associated with a multiprecision arithmetic floating point are very useful, because it yields a clear reduction in iterations. A motivation for an arbitrary precision in interval methods can be found in [6], in particular for the calculation of zeros of nonlinear functions.

Recently, many robust and efficient methods with high convergence order have been proposed to solve nonlinear equations, but in most of cases the schemes cannot be extended to multivariate problems. Few papers for the multidimensional case introduce methods with high order of convergence. The authors design in [7] a modified Newton-Jarrat scheme of sixth order; in [8] a third-order method is presented for computing real and complex roots of nonlinear systems; Shin et al. compare in [9] Newton-Krylov methods and Newton-like schemes for solving big-sized nonlinear systems; the authors in [10] and A. Iliev and I. Iliev in [11] show general procedures to design high-order methods by using frozen Jacobian and Taylor expansion, respectively. Special case of sparse Jacobian matrices is studied in [12].

Dayton et al. in [13] formulate the multiplicity for the general nonlinear system at an isolated zero. They present an algorithm for computing the multiplicity structure, propose a depth-deflation method for accurate computation of multiple zeros, and introduce the basic algebraic theory of the multiplicity.

In this paper, we present three new Newton-like schemes, of order of convergence four, six, and eight, respectively. After the analysis of convergence of the new methods, we apply the pseudocomposition technique in order to get higher-order procedures. This technique (see [14]) consists of the following: we consider a method of order of convergence 𝑝 as a predictor, whose penultimate step is of order 𝑞, and then we use a corrector step based on the Gaussian quadrature. So, we obtain a family of iterative schemes whose order of convergence is min{𝑞+𝑝,3𝑞}. This is a general procedure to improve the order of convergence of known methods.

To analyze and compare the efficiency of the proposed methods we use the classic efficiency index 𝐼=𝑝^1/𝑑 due to Ostrowski [15], where 𝑝 is the order of convergence and 𝑑 is the number of functional evaluations at each iteration.

The convergence theorem in Section 2 is demonstrated by means of the 𝑛-dimensional Taylor expansion of the functions involved. Let 𝐹∶𝐷⊆𝑅^𝑛→𝑅^𝑛 be sufficiently Frechet differentiable in 𝐷. By using the notation introduced in [7], the 𝑞th derivative of 𝐹 at 𝑢∈𝑅^𝑛, 𝑞≥1, is the 𝑞-linear function 𝐹^(𝑞)(𝑢)∶𝑅^𝑛×⋯×𝑅^𝑛→𝑅^𝑛 such that 𝐹^(𝑞)(𝑢)(𝑣₁,…,𝑣_𝑞)∈𝑅^𝑛. It is easy to observe that(1)𝐹^(𝑞)(𝑢)(𝑣₁,…,𝑣_𝑞−1,⋅)∈ℒ(𝑅^𝑛), (2)𝐹^(𝑞)(𝑢)(𝑣_𝜎(1),…,𝑣_𝜎(𝑞))=𝐹^(𝑞)(𝑢)(𝑣₁,…,𝑣_𝑞), for all permutation 𝜎 of {1,2,…,𝑞}.

So, in the following we will denote:(a)𝐹^(𝑞)(𝑢)(𝑣₁,…,𝑣_𝑞)=𝐹^(𝑞)(𝑢)𝑣₁⋯𝑣_𝑞, (b)𝐹^(𝑞)(𝑢)𝑣^𝑞−1𝐹^(𝑝)𝑣^𝑝=𝐹^(𝑞)(𝑢)𝐹^(𝑝)(𝑢)𝑣^{𝑞+𝑝−1}.

It is well known that, for 𝜉+ℎ∈𝑅^𝑛 lying in a neighborhood of a solution 𝜉 of the nonlinear system 𝐹(𝑥)=0, Taylor’s expansion can be applied (assuming that the jacobian matrix 𝐹^′(𝜉) is nonsingular), and 𝐹(𝜉+ℎ)=𝐹^′(𝜉)ℎ+^𝑝−1_𝑞=2𝐶_𝑞ℎ^𝑞+𝑂ℎ^𝑝, (1.1) where 𝐶_𝑞=(1/𝑞!)[𝐹^′(𝜉)]⁻¹𝐹^(𝑞)(𝜉), 𝑞≥2. We observe that 𝐶_𝑞ℎ^𝑞∈𝑅^𝑛 since 𝐹^(𝑞)(𝜉)∈ℒ(𝑅^𝑛×⋯×𝑅^𝑛,𝑅^𝑛) and [𝐹^′(𝜉)]⁻¹∈ℒ(𝑅^𝑛).

In addition, we can express the Jacobian matrix of 𝐹, 𝐹^′ as 𝐹^′(𝜉+ℎ)=𝐹^′(𝜉)𝐼+^𝑝−1_𝑞=2𝑞𝐶_𝑞ℎ^𝑞−1+𝑂ℎ^𝑝, (1.2) where 𝐼 is the identity matrix. Therefore, 𝑞𝐶_𝑞ℎ^𝑞−1∈ℒ(𝑅^𝑛). From (1.2), we obtain 𝐹^′(𝜉+ℎ)⁻¹=𝐼+𝑋₂ℎ+𝑋₃ℎ²+𝑋₄ℎ³+⋯𝐹^′(𝜉)⁻¹+𝑂ℎ^𝑝, (1.3) where 𝑋₂=−2𝐶₂,𝑋₃=4𝐶²₂−3𝐶₃,….

We denote 𝑒_𝑘=𝑥^(𝑘)−𝜉 the error in the 𝑘th iteration. The equation 𝑒_(𝑘+1)=𝐿𝑒_𝑘^𝑝+𝑂[𝑒_𝑘^𝑝+1], where 𝐿 is a 𝑝-linear function 𝐿∈ℒ(𝑅^𝑛×⋯×𝑅^𝑛,𝑅^𝑛), is called the error equation and 𝑝 is the order of convergence.

The rest of the paper is organized as follows: in the next section, we present the new methods of order four, six, and eight, respectively. Moreover, the convergence order is increased when the pseudocomposition technique is applied. Section 3 is devoted to the comparison of the different methods by means of several numerical tests.

2. Design and Convergence Analysis of the New Methods

Let us introduce now a new Jarratt-type scheme of five steps which we will denote as M8. We will prove that its first three steps define a fourth-order scheme, denoted by M4, and its four first steps become a sixth-order method that will be denoted by M6. The coefficients involved have been obtained optimizing the order of the convergence, and the whole scheme requires three functional evaluations of 𝐹 and two of 𝐹^′ to attain eighth order of convergence. Let us also note that the linear systems to be solved in first, second, and last step have the same matrix and also have the third and fourth steps, so the number of operations involved is not as high as it can seem.

Theorem 2.1. Let 𝐹∶Ω⊆𝑅^𝑛→𝑅^𝑛 be a sufficiently differentiable in a neighborhood of 𝜉∈Ω which is a solution of the nonlinear system 𝐹(𝑥)=0, and let 𝑥⁽⁰⁾ be an initial estimation close enough to the solution 𝜉. One also supposes that 𝐹^′(𝑥) is continuous and nonsingular at 𝜉. Then, the sequence {𝑥^(𝑘)}_𝑘≥0 obtained by 𝑦^(𝑘)=𝑥^(𝑘)−23𝐹^′𝑥^(𝑘)−1𝐹𝑥^(𝑘),𝑧^(𝑘)=𝑦^(𝑘)+16𝐹^′𝑥^(𝑘)−1𝐹𝑥^(𝑘),𝑢^(𝑘)=𝑧^(𝑘)+𝐹^′𝑥^(𝑘)−3𝐹^′𝑦^(𝑘)−1𝐹𝑥^(𝑘),𝑣^(𝑘)=𝑧^(𝑘)+𝐹^′𝑥^(𝑘)−3𝐹^′𝑦^(𝑘)−1𝐹𝑥^(𝑘)+2𝐹𝑢^(𝑘)𝑥^(𝑘+1)=𝑣^(𝑘)−12𝐹^′𝑥^(𝑘)−15𝐹^′𝑥^(𝑘)−3𝐹^′𝑦^(𝑘)𝐹^′𝑥^(𝑘)−1𝐹𝑣^(𝑘) (2.1) converges to 𝜉 with order of convergence eight. The error equation is 𝑒_𝑘+1=𝐶²₂−12𝐶₃2𝐶³₂+2𝐶₃𝐶₂−2𝐶₂𝐶₃−209𝐶₄𝑒⁸_𝑘+𝑂𝑒⁹_𝑘. (2.2)

Proof. From (1.1) and (1.2) we obtain 𝐹𝑥^(𝑘)=𝐹^′(𝜉)𝑒_𝑘+𝐶₂𝑒²_𝑘+𝐶₃𝑒³_𝑘+𝐶₄𝑒⁴_𝑘+𝐶₅𝑒⁵_𝑘+𝐶₆𝑒⁶_𝑘+𝐶₇𝑒⁷_𝑘+𝐶₈𝑒⁸_𝑘+𝑂𝑒⁹_𝑘,𝐹^′𝑥^(𝑘)=𝐹^′(𝜉)𝐼+2𝐶₂𝑒_𝑘+3𝐶₃𝑒²_𝑘+4𝐶₄𝑒³_𝑘+5𝐶₅𝑒⁴_𝑘+6𝐶₆𝑒⁵_𝑘+7𝐶₇𝑒⁶_𝑘+8𝐶₈𝑒⁷_𝑘+𝑂𝑒⁸_𝑘. (2.3) As [𝐹^′(𝑥^(𝑘))]⁻¹𝐹^′(𝑥^(𝑘))=𝐼, we calculate 𝐹^′𝑥^(𝑘)−1=𝐼+𝑋₂𝑒_𝑘+𝑋₃𝑒²_𝑘+𝑋₄𝑒³_𝑘+𝑋₅𝑒⁴_𝑘+𝑋₆𝑒⁵_𝑘+𝑋₇𝑒⁶_𝑘+𝑋₈𝑒⁷_𝑘𝐹^′(𝜉)⁻¹+𝑂𝑒⁸_𝑘, (2.4) where 𝑋₁=𝐼 and 𝑋_𝑠=−∑^𝑠_𝑗=2𝑗𝑋_{𝑠−𝑗+1}𝐶_𝑗, for 𝑠=2,3,… So, 𝐹^′𝑥^(𝑘)−1𝐹𝑥^(𝑘)=𝑒_𝑘+𝑀₂𝑒²_𝑘+𝑀₃𝑒³_𝑘+𝑀₄𝑒⁴_𝑘+𝑀₅𝑒⁵_𝑘+𝑀₆𝑒⁶_𝑘+𝑀₇𝑒⁷_𝑘+𝑀₈𝑒⁸_𝑘+𝑂𝑒⁹_𝑘, (2.5) where 𝑀₂=𝐶₂+𝑋₄ and 𝑀_𝑠=𝐶_𝑠+∑^𝑠_𝑗=3𝑋_{𝑠−𝑗+2}𝐶_𝑗−1+𝑋_𝑠+2,𝑠=3,4,….
Then, 𝑦^(𝑘)=𝜉+(1/3)𝑒_𝑘−(2/3)𝑀 and 𝑧^(𝑘)=𝜉+(1/2)𝑒_𝑘−(1/2)𝑀, where 𝑀=𝑀₂𝑒²_𝑘+𝑀₃𝑒³_𝑘+𝑀₄𝑒⁴_𝑘+𝑀₅𝑒⁵_𝑘+𝑀₆𝑒⁶_𝑘+𝑀₇𝑒⁷_𝑘+𝑀₈𝑒⁸_𝑘+𝑂[𝑒⁹_𝑘].
The Taylor expansion of 𝐹^′(𝑦^(𝑘)) is 𝐹^′𝑦^(𝑘)=𝐹^′(𝜉)𝐼+𝑄₁𝑒_𝑘+𝑄₂𝑒²_𝑘+𝑄₃𝑒³_𝑘+𝑄₄𝑒⁴_𝑘+𝑄₅𝑒⁵_𝑘+𝑄₆𝑒⁶_𝑘+𝑄₇𝑒⁷_𝑘+𝑄₈𝑒⁸_𝑘+𝑂𝑒⁹_𝑘, (2.6) where 𝑄₁=23𝐶₂,𝑄₂=13𝐶₃−43𝐶₂𝑀₂,𝑄₃=4𝐶27⁴−43𝐶₃𝑀₂−43𝐶₂𝑀₃,𝑄₄=5𝐶81⁵−89𝐶₄𝑀₂+43𝐶₃𝑀²₂−2𝑀₃−43𝐶₂𝑀₄,𝑄₅=2𝐶81⁶−40𝐶81⁵𝑀₂+49𝐶₄𝑀²₂−𝑀₃+43𝐶₃𝑀₂𝑀₃+𝑀₃𝑀₂−𝑀₄−43𝐶₂𝑀₅,𝑄₆=7𝐶729⁷−20𝐶81⁶𝑀₂+40𝐶81⁵3𝑀²₂−𝑀₃+8𝐶27⁴6𝑀₂𝑀₃+6𝑀₃𝑀₂−3𝑀₄−4𝑀³₂−𝐶₂𝑀₆,𝑄₇=43𝐶₂𝑀₇+43𝐶₃𝑀₂𝑀₅+𝑀₃𝑀₄+𝑀₄𝑀₃+𝑀₅𝑀₂+𝑀₆𝑀₂−𝑀₇+8𝐶27⁴6𝑀₂𝑀₄+6𝑀²₃+6𝑀₄𝑀₂−3𝑀₅−4𝑀²₂𝑀₃−4𝑀₂𝑀₃𝑀₂−4𝑀₃𝑀²₂+40𝐶81⁵3𝑀₂𝑀₃+3𝑀₃𝑀₂−𝑀₄−4𝑀³₂+20𝐶81⁶4𝑀²₂−𝑀₃−28𝐶243⁷𝑀₂+8𝐶2187⁸,𝑄₈4=−3𝐶₂𝑀₈+43𝐶₃𝑀₂𝑀₆+𝑀₃𝑀₅+𝑀²₄+𝑀₅𝑀₃+𝑀₆𝑀₂−𝑀₇+169𝐶₄𝑀₂𝑀₅+𝑀₃𝑀₄+𝑀₄𝑀₃+𝑀₅𝑀₂−𝑀₆+32𝐶27⁴𝑀²₂𝑀₄+𝑀₂𝑀²₃𝑀₄𝑀₂+𝑀₃𝑀₂𝑀₃+𝑀²₃𝑀₂+𝑀₄𝑀²₂+40𝐶81⁵−𝑀₅−4𝑀²₂𝑀₃−4𝑀₂𝑀₃𝑀₂−4𝑀₃𝑀²₂+3𝑀₂𝑀₄+3𝑀²₃+3𝑀₄𝑀₂+2𝑀₄+20𝐶41⁶4𝑀₂𝑀₃+4𝑀₃𝑀₂−𝑀₄−8𝑀³₂+28𝐶243⁷5𝑀²₂−𝑀₃−112𝐶2187⁸𝑀₂+1𝐶729⁹. (2.7) We also obtain the Taylor expansion of 𝐹^′(𝑥^(𝑘))−3𝐹^′(𝑦^(𝑘)): 𝐹^′𝑥^(𝑘)−3𝐹^′𝑦^(𝑘)=𝐹^′(𝜉)−2𝐼+𝐴₁𝑒_𝑘+𝐴₂𝑒²_𝑘+𝐴₃𝑒³_𝑘+𝐴₄𝑒⁴_𝑘+𝐴₅𝑒⁵_𝑘+𝐴₆𝑒⁶_𝑘+𝐴₇𝑒⁷_𝑘+𝐴₈𝑒⁸_𝑘+𝑂𝑒⁹_𝑘, (2.8) where 𝐴_𝑠=(𝑠+1)𝐶_𝑠+1−3𝑄_𝑠,𝑠=1,2,…. As [𝐹^′(𝑥^(𝑘))−3𝐹^′(𝑦^(𝑘))]⁻¹[𝐹^′(𝑥^(𝑘))−3𝐹^′(𝑦^(𝑘))]=𝐼, we obtain 𝐹^′𝑥^(𝑘)−3𝐹^′𝑦^(𝑘)−1=−12𝐼+𝑌₂𝑒_𝑘+𝑌₃𝑒²_𝑘+𝑌₄𝑒³_𝑘+𝑌₅𝑒⁴_𝑘+𝑌₆𝑒⁵_𝑘+𝑌₇𝑒⁶_𝑘+𝑌₈𝑒⁷_𝑘𝐹^′(𝜉)⁻¹+𝑂𝑒⁸_𝑘, (2.9) where 𝑌₂=0 and 𝑌_𝑠=(1/2)∑^𝑠_𝑗=3𝑌_{𝑠−𝑗+2}𝐴_𝑗−2−(1/4)𝐴_𝑠−1,𝑠=3,4,….
So, 𝐹^′𝑥^(𝑘)−3𝐹^′𝑦^(𝑘)−1𝐹𝑥^(𝑘)1=−2𝑒_𝑘+𝑅₂𝑒²_𝑘+𝑅₃𝑒³_𝑘+𝑅₄𝑒⁴_𝑘+𝑅₅𝑒⁵_𝑘+𝑅₆𝑒⁶_𝑘+𝑅₇𝑒⁷_𝑘+𝑅₈𝑒⁸_𝑘+𝑂𝑒⁹_𝑘, (2.10) where 𝑅₂=𝑌₂−(1/2)𝐶₂ and 𝑅_𝑠=𝑌_𝑠+∑^𝑠_𝑗=3𝑌_{𝑠−𝑗+2}𝐶_𝑗−1−(1/2)𝐶_𝑠,𝑠=3,4,….
We now calculate 𝑢^(𝑘)=𝑧^(𝑘)+[𝐹^′(𝑥^(𝑘))−3𝐹^′(𝑦^(𝑘))]⁻¹𝐹(𝑥^(𝑘)), and the error equation of the method at this step is 𝑒_𝑢(𝑘)=12𝑒_𝑘−12𝑀₂𝑒²_𝑘+𝑀₃𝑒³_𝑘+𝑀₄𝑒⁴_𝑘+𝑀₅𝑒⁵_𝑘+𝑀₆𝑒⁶_𝑘+𝑀₇𝑒⁷_𝑘+𝑀₈𝑒⁸_𝑘−12𝑒_𝑘+𝑅₂𝑒²_𝑘+𝑅₃𝑒³_𝑘+𝑅₄𝑒⁴_𝑘+𝑅₅𝑒⁵_𝑘+𝑅₆𝑒⁶_𝑘+𝑅₇𝑒⁷_𝑘+𝑅₈𝑒⁸_𝑘+𝑂𝑒⁹_𝑘=𝑃₄𝑒⁴_𝑘+𝑃₅𝑒⁵_𝑘+𝑃₆𝑒⁶_𝑘+𝑃₇𝑒⁷_𝑘+𝑃₈𝑒⁸_𝑘+𝑂𝑒⁹_𝑘, (2.11) where 𝑃_𝑠=−(1/2)𝑀_𝑠+𝑅_𝑠, 𝑠=4,5,… Then the first three steps define a fourth-order procedure, and 𝐹𝑢^(𝑘)=𝐹^′(𝜉)𝑃₄𝑒⁴_𝑘+𝑃₅𝑒⁵_𝑘+𝑃₆𝑒⁶_𝑘+𝑃₇𝑒⁷_𝑘+𝑃₈+𝑃²₄𝑒⁸_𝑘+𝑂𝑒⁹_𝑘,𝐹𝑥^(𝑘)+2𝐹𝑢^(𝑘)=𝐹^′(𝜉)𝑒_𝑘+𝐶₂𝑒²_𝑘+𝐶₃𝑒³_𝑘+𝐶₄+2𝑃₄𝑒⁴_𝑘+𝐶₅+2𝑃₅𝑒⁵_𝑘+𝐶₆+2𝑃₆𝑒⁶_𝑘+𝐶₇+2𝑃₇𝑒⁷_𝑘+𝐶₈+2𝑃₈𝑒⁸_𝑘+𝑂𝑒⁹_𝑘. (2.12) So, 𝐹^′𝑥^(𝑘)−3𝐹^′𝑦^(𝑘)−1𝐹𝑥^(𝑘)+2𝐹𝑢^(𝑘)1=−2𝑒_𝑘+𝐿₂𝑒²_𝑘+𝐿₃𝑒³_𝑘+𝐿₄𝑒⁴_𝑘+𝐿₅𝑒⁵_𝑘+𝐿₆𝑒⁶_𝑘+𝐿₇𝑒⁷_𝑘+𝐿₈𝑒⁸_𝑘+𝑂𝑒⁹_𝑘, (2.13) where 𝐿₂=𝑌₂−(1/2)(𝐶₂+𝑃₂) and 𝐿_𝑠=𝑌_𝑠+∑^𝑠_𝑗=3𝑌_{𝑠−𝑗+2}𝐶_𝑗−1−(1/2)(𝐶_𝑠+𝑃_𝑠),𝑠=3,4,….
At the fourth step, we calculate 𝑒_𝑣(𝑘)=𝑧_𝑘−𝜉+𝐹^′𝑥^(𝑘)−3𝐹^′𝑦^(𝑘)−1𝐹𝑥^(𝑘)+2𝐹𝑢^(𝑘)=𝑁₆𝑒⁶_𝑘+𝑁₇𝑒⁷_𝑘+𝑁₈𝑒⁸_𝑘+𝑂𝑒⁹_𝑘, (2.14) where 𝑁_𝑠=𝐿_𝑠−(1/2)𝑀_𝑠 and 𝑠=6,7,… This shows that the first four steps of the method define a sixth-order scheme. Indeed, 𝐹𝑣^(𝑘)=𝐹^′(𝜉)𝑁₆𝑒⁶_𝑘+𝑁₇𝑒⁷_𝑘+𝑁₈𝑒⁸_𝑘+𝑂𝑒⁹_𝑘,5𝐹^′𝑥^(𝑘)−3𝐹^′𝑦^(𝑘)=𝐹^′(𝜉)2𝐼+𝐵₁𝑒_𝑘+𝐵₂𝑒²_𝑘+𝐵₃𝑒³_𝑘+𝐵₄𝑒⁴_𝑘+𝐵₅𝑒⁵_𝑘+𝐵₆𝑒⁶_𝑘+𝐵₇𝑒⁷_𝑘+𝐵₈𝑒⁸_𝑘+𝑂𝑒⁹_𝑘, (2.15)(2.16) where 𝐵_𝑠=5(𝑠+1)𝐶_𝑠+1−3𝑄_𝑠 and 𝑠=1,2… Then, 𝐹^′𝑥^(𝑘)−15𝐹^′𝑥^(𝑘)−3𝐹^′𝑦^(𝑘)=2𝐼+𝑆₁𝑒_𝑘+𝑆₂𝑒²_𝑘+𝑆₃𝑒³_𝑘+𝑆₄𝑒⁴_𝑘+𝑆₅𝑒⁵_𝑘+𝑆₆𝑒⁶_𝑘+𝑆₇𝑒⁷_𝑘+𝑆₈𝑒⁸_𝑘+𝑂𝑒⁹_𝑘,𝐹^′𝑥^(𝑘)−1𝐹𝑣^(𝑘)=2𝑁₆𝑒⁶_𝑘+𝑁₇+𝑋₂𝑁₆𝑒⁷_𝑘+𝑁₈+𝑋₂𝑁₇+𝑋₃𝑁₆𝑒⁸_𝑘+𝑂𝑒⁹_𝑘, (2.17)(2.18) where 𝑆₁=𝐵₁+2𝑋₂ and 𝑆_𝑠=𝐵_𝑠+∑^𝑠_𝑗=2𝑋_{𝑠−𝑗+2}𝐵_𝑗−1+2𝑋_𝑠+1, 𝑠=2,3,…. By multiplying (2.17) and (2.18), 𝐹^′𝑥^(𝑘)−15𝐹^′𝑥^(𝑘)−3𝐹^′𝑦^(𝑘)𝐹^′𝑥^(𝑘)−1𝐹𝑣^(𝑘)=2𝑁₆𝑒⁶_𝑘+𝐾₇𝑒⁷_𝑘+𝐾₈𝑒⁸_𝑘+𝑂𝑒⁹_𝑘, (2.19) where 𝐾₇=2(𝑁₇+𝑋₂𝑁₆)+𝑆₁𝑁₆ and 𝐾₈=2(𝑁₈+𝑋₂𝑁₇+𝑋₃𝑁₆)+𝑆₁(𝑁₇+𝑋₂𝑁₆)+𝑆₂𝑁₆.
Finally, 𝑥^(𝑘+1)=𝜉+𝐿₈−12𝑀₈−122𝑁₈+𝑋₂𝑁₇+𝑋₃𝑁₆+𝑆₁𝑁₇+𝑋₂𝑁₆+𝑆₂𝑁₆𝑒⁸_𝑘+𝑂𝑒⁹_𝑘, (2.20) and the error equation is 𝑒_𝑘+1=𝐶²₂−12𝐶₃2𝐶³₂+2𝐶₃𝐶₂−2𝐶₂𝐶₃−209𝐶₄𝑒⁸_𝑘+𝑂𝑒⁹_𝑘, (2.21) proving that the order of convergence of the analyzed method is eight.

In [14] the authors presented a new procedure to design higher-order schemes. This technique, called pseudocomposition, uses the two last steps of the predictor method to obtain a corrected scheme with higher order of convergence.

Theorem 2.2 (see [14]). Let 𝐹∶Ω⊆𝑅^𝑛→𝑅^𝑛 be differentiable enough Ω, let 𝜉∈Ω be a solution of the nonlinear system 𝐹(𝑥)=0, and let 𝑥⁽⁰⁾ be an initial estimation close enough to the solution 𝜉. We suppose that 𝐹^′(𝑥) is continuous and nonsingular at 𝜉. Let 𝑦^(𝑘) and 𝑧^(𝑘) be the penultimate and final steps of orders 𝑞 and 𝑝, respectively, of a certain iterative method. Taking this scheme as a predictor we get a new approximation 𝑥^(𝑘+1) of 𝜉 given by 𝑥^(𝑘+1)=𝑦^(𝑘)−2^𝑚_𝑖=1𝜔_𝑖𝐹^′𝜂_𝑖^(𝑘)−1𝐹𝑦^(𝑘), (2.22) where 𝜂_𝑖^(𝑘)=(1/2)[(1+𝜏_𝑖)𝑧^(𝑘)+(1−𝜏_𝑖)𝑦^(𝑘)] and 𝜏_𝑖, 𝜔_𝑖,𝑖=1,2,…,𝑚 are the nodes and weights of the orthogonal polynomial corresponding to the Gaussian quadrature used. Then,(1)the obtained set of families will have an order of convergence at least 𝑞; (2)if 𝜎=2 is satisfied, then the order of convergence will be at least 2𝑞; (3)if, also, 𝜎₁=0, the order of convergence will be min{𝑝+𝑞,3𝑞},
where ∑^𝑛_𝑖=1𝜔_𝑖=𝜎 and ∑^𝑛_𝑖=1𝜔_𝑖𝜏^𝑗_𝑖/𝜎=𝜎_𝑗with𝑗=1,2.

Depending on the orthogonal polynomial corresponding to the Gaussian quadrature used in the corrector step, this procedure will determine a family of schemes. Furthermore, it is possible to obtain different methods in these families by using distinct number of nodes corresponding to the orthogonal polynomial used (see Table 1). However, according to the proof of Theorem 2.2 the order of convergence of the obtained methods does not depend on the number of nodes used.

Table 1 

Values of the parameters for the different quadratures used.

Let us note that these methods, obtained by means of Gaussian quadratures, seem to be known interpolation quadrature schemes such as midpoint, trapezoidal, or Simpson’s method (see [16]). It is only a similitude, as they are not applied on the last iteration 𝑥^(𝑘), and the last step of the predictor 𝑧^(𝑘), but on the two last steps of the predictor. In the following, we will use a midpoint-like as a corrector step, which corresponds to a Gauss-Legendre quadrature with one node; for this scheme the order of convergence will be at least min{𝑞+𝑝,3𝑞}, by applying Theorem 2.2.

The pseudocomposition can be applied to the proposed scheme M8 with iterative expression (2.1), but also to M6. By pseudocomposing on M6 and M8 there can be obtained two procedures of order of convergence 10 and 14 (denoted by PsM10 and PsM14), respectively. Let us note that it is also possible to pseudocompose on M4, but the resulting scheme would be of third order of convergence, which is worse than the original M4, so it will not be considered.

Following the notation used in (2.1), the last step of PsM10 is 𝑥^(𝑘+1)=𝑢^(𝑘)−2𝐹^′𝑣^(𝑘)+𝑢^(𝑘)2⁻¹𝐹𝑢^(𝑘), (2.23) and the last three steps of psM14 can be expressed as 𝑣^(𝑘)=𝑧^(𝑘)+𝐹^′𝑥^(𝑘)−3𝐹^′𝑦^(𝑘)−1𝐹𝑥^(𝑘)+2𝐹𝑢^(𝑘),𝑤^(𝑘+1)=𝑣^(𝑘)−12𝐹^′𝑥^(𝑘)−15𝐹^′𝑥^(𝑘)−3𝐹^′𝑦^(𝑘)𝐹^′𝑥^(𝑘)−1𝐹𝑣^(𝑘),𝑥^(𝑘+1)=𝑣^(𝑘)−2𝐹^′𝑤^(𝑘)+𝑣^(𝑘)2⁻¹𝐹𝑣^(𝑘). (2.24)

In Figure 1, we analyze the efficiency indices of the proposed methods, compared with Newton and Jarrat’s schemes and between themselves. There can be deduced the following conclusions: the new methods M4, M6, and M8 (and also the pseudocomposed PsM10 and PsM14) improve Newton and Jarratt’s schemes (in fact, the indices of M4 and Jarratt’s are equal). Indeed, for 𝑛≥3 the best index is that of M8. Nevertheless, none of the pseudocomposed methods improve the efficiency index of their original partners. Nevertheless, as we will see in the following section, the pseudocomposed schemes will show a very stable behavior that makes them worth.

Figure 1 

Efficiency index of the different methods for different sizes of the system.

3. Numerical Results

In order to illustrate the effectiveness of the proposed methods, we will compare them with other known schemes. Numerical computations have been performed in MATLAB R2011a by using variable-precision arithmetic, which uses floating-point representation of 2000 decimal digits of mantissa. The computer specifications are Intel(R) Core(TM) i5-2500 CPU @ 3.30 GHz with 16.00 GB of RAM. Each iteration is obtained from the former by means of an iterative expression 𝑥^(𝑘+1)=𝑥^(𝑘)−𝐴⁻¹𝑏, where 𝑥^(𝑘)∈ℝ^𝑛, 𝐴 is a real matrix 𝑛×𝑛 and 𝑏∈ℝ^𝑛. The matrix 𝐴 and vector 𝑏 are different according to the method used, but in any case, we calculate 𝐴⁻¹𝑏 as the solution of the linear system 𝐴𝑦=𝑏, with Gaussian elimination with partial pivoting. The stopping criterion used is ||𝑥^(𝑘+1)−𝑥^(𝑘)||<10⁻²⁰⁰ or ||𝐹(𝑥^(𝑘))||<10⁻²⁰⁰.

Firstly, let us consider the following nonlinear systems of different sizes:(1)𝐹₁=(𝑓₁(𝑥),𝑓₂(𝑥),…,𝑓_𝑛(𝑥)), where 𝑥=(𝑥₁,𝑥₂,…,𝑥_𝑛)^𝑇 and 𝑓_𝑖∶ℝ^𝑛→ℝ,𝑖=1,2,…,𝑛, such that 𝑓_𝑖(𝑥)=𝑥_𝑖𝑥_𝑖+1𝑓−1,𝑖=1,2,…,𝑛−1,_𝑛(𝑥)=𝑥_𝑛𝑥₁−1. (3.1) When 𝑛 is odd, the exact zeros of 𝐹₁(𝑥) are: 𝜉₁=(1,1,…,1)^𝑇 and 𝜉₂=(−1,−1,…,−1)^𝑇. (2)𝐹₂(𝑥₁,𝑥₂)=(𝑥²₁−𝑥₁−𝑥²₂−1,−sin(𝑥₁)+𝑥₂) and the solutions are 𝜉₁≈(−0.845257,−0.748141)^𝑇 and 𝜉₂≈(1.952913,0.927877)^𝑇. (3)𝐹₃(𝑥₁,𝑥₂)=(𝑥²₁+𝑥²₂−4,−exp(𝑥₁)+𝑥₂−1), being the solutions 𝜉₁≈(1.004168,−1.729637)^𝑇 and 𝜉₂≈(−1.816264,0.837368)^𝑇. (4)𝐹₄(𝑥₁,𝑥₂,𝑥₃)=(𝑥²₁+𝑥²₂+𝑥²₃−9,𝑥₁𝑥₂𝑥₃−1,𝑥₁+𝑥₂−𝑥²₃) with roots 𝜉₁≈(2.14025,−2.09029,−0.223525)^𝑇, 𝜉₂≈(2.491376,0.242746,1.653518)^𝑇 and 𝜉₁≈(0.242746,2.491376,1.653518)^𝑇.

Table 2 presents results showing the following information: the different iterative methods employed (Newton (NC), Jarratt (JT), the new methods M4, M6, and M8, and the pseudocomposed PsM10 and PsM14), the number of iterations 𝐼𝑡𝑒𝑟 needed to converge to the solution Sol, the value of the stopping factors at the last step, and the computational order of convergence 𝜌 (see [17]) approximated by the formula: 𝑥𝜌≈ln^(𝑘+1)−𝑥^(𝑘)/𝑥^(𝑘)−𝑥^(𝑘−1)ln𝑥^(𝑘)−𝑥^(𝑘−1)/𝑥^(𝑘−1)−𝑥^(𝑘−2). (3.2) The value of 𝜌 which appears in Table 2 is the last coordinate of the vector 𝜌 when the variation between their coordinates is small. Also the elapsed time, in seconds, appears in Table 2, being the mean execution time for 100 performances of the method (the command cputime of MATLAB has been used).

Table 2 

Numerical results for functions 𝐹₁ to 𝐹₄.

We observe from Table 2 that not only the order of convergence and the number of new functional evaluations and operations are important in order to obtain new efficient iterative methods to solve nonlinear systems of equations. A key factor is the range of applicability of the methods. Although they are slower than the original methods when the initial estimation is quite good, when we are far from the solution or inside a region of instability, the original schemes do not converge or do it more slowly, the corresponding pseudocomposed procedures usually still converge or do it faster.

The advantage of pseudocomposition can be observed in Figures 2(a) and 2(b) (methods M6 and PsM10) and Figures 3(a) and 3(b) (methods M8 and PsM14) where the dynamical plane on 𝑅² is represented: we consider the system of two equations and two unknowns 𝐹₂(𝑥)=0, for any initial estimation in 𝑅² represented by its position in the plane, a different color (blue or orange, as there exist only two solutions in this region) is used for the different solutions found (marked by a white point in the figure). Black color represents an initial point in which the method converges to infinity, and the green one means that no convergence is found (usually because any linear system cannot be solved). It is clear that when many initial estimations tend to infinity (see Figure 3(a)), the pseudocomposition “cleans” the dynamical plane, making the method more stable as it can find one of the solutions by using starting points that do not allow convergence with the original scheme (see Figure 3(b)).

(a) M6

(b) PsM10

Figure 2 

Real dynamical planes for system 𝐹2(𝑥)=0 and methods M6 and PsM10.

(a)  M8

(b) PsM14

Figure 3 

Real dynamical planes for system 𝐹2(𝑥)=0 and methods M8 and PsM14.

If an analogous study is made on system 𝐹₃(𝑥)=0, similar conclusions can be obtained, as the effect of smoothness is clear when the real dynamical plane of a method and its pseudocomposed partner are compared. So, in Figure 4 the amount of points in the lower half of the plane that converge to one of the roots is higher after the pseudocomposition, and, in Figure 5, there is a big green region of no convergence for method M8 that shows to be convergent when pseudocomposition is applied in PsM14.

(a) M6

(b) PsM10

Figure 4 

Real dynamical planes for system 𝐹3(𝑥)=0 and methods M6 and PsM10.

(a) M8

(b) PsM14

Figure 5 

Real dynamical planes for system 𝐹3(𝑥)=0 and methods M8 and PsM14.

We conclude that the presented schemes M4, M6, and M8 show to be excellent, in terms of order of convergence and efficiency, but also that the pseudocomposition technique achieves to transform them in competent and more robust new schemes.

Acknowledgments

The authors would like to thank the referees for the valuable comments and for the suggestions to improve the readability of the paper. This research was supported by Ministerio de Ciencia y Tecnología MTM2011-28636-C02-02 and FONDOCYT República Dominicana.