An Efficient Family of Traub-Steffensen-Type Methods for Solving Systems of Nonlinear Equations

Abstract

Based on Traub-Steffensen method, we present a derivative free three-step family of sixth-order methods for solving systems of nonlinear equations. The local convergence order of the family is determined using first-order divided difference operator for functions of several variables and direct computation by Taylor's expansion. Computational efficiency is discussed, and a comparison between the efficiencies of the proposed techniques with the existing ones is made. Numerical tests are performed to compare the methods of the proposed family with the existing methods and to confirm the theoretical results. It is shown that the new family is especially efficient in solving large systems.

1. Introduction

The problem of finding solution of the system of nonlinear equations 𝐹(𝑥)=0, where 𝐹∶𝐷→𝐷, 𝐷 is an open convex domain in 𝑅^𝑛, by iterative methods is an important and challenging task in numerical analysis and many applied scientific branches. One of the basic procedures for solving nonlinear equations is the quadratically convergent Newton method (see [1, 2]): 𝑥^(𝑘+1)=𝑥^(𝑘)−𝐹^′𝑥^(𝑘)−1𝐹𝑥^(𝑘),𝑘=0,1,2,…, (1) where [𝐹^′(𝑥)]⁻¹ is the inverse of the first Fréchet derivative 𝐹^′(𝑥) of the function 𝐹(𝑥).

In many practical situations, it is preferable to avoid the calculation of derivative 𝐹^′(𝑥) of the function 𝐹(𝑥). In such situations, it is preferable to use only the computed values of 𝐹(𝑥) and to approximate 𝐹^′(𝑥) by employing the values of 𝐹(𝑥) at suitable points. For example, a basic derivative free iterative method is the Traub-Steffensen method [3], which also converges quadratically and follows the scheme 𝑥^(𝑘+1)=𝐺_1,2𝑥^(𝑘)=𝑥^(𝑘)−𝑤^(𝑘),𝑥^(𝑘);𝐹⁻¹𝐹𝑥^(𝑘), (2) where [𝑤^(𝑘),𝑥^(𝑘);𝐹]⁻¹ is the inverse of the first-order divided difference [𝑤^(𝑘),𝑥^(𝑘);𝐹] of 𝐹 and 𝑤^(𝑘)=𝑥^(𝑘)+𝛾𝐹(𝑥^(𝑘));𝛾 is an arbitrary nonzero constant. Throughout this paper, 𝐺_𝑖,𝑝 is used to denote the 𝑖th iteration function of convergence order 𝑝. For 𝛾=1, the scheme (2) reduces to the well-known Steffensen method [4].

In recent years, many derivative free higher-order methods of great efficiency are developed for solving scalar equation 𝑓(𝑥)=0; see [5–14] and the references therein. For systems of nonlinear equations, however, the construction of efficient higher-order derivative free methods is a difficult task and therefore not many such methods can be found in the literature. Recently, based on Steffensen's scheme, that is, when 𝛾=1 in (2), a family of seventh-order methods has been proposed in [13]. Some important members of this family, as shown in [13], are given as follows: 𝑦^(𝑘)=𝐺_1,2𝑥^(𝑘),𝑧^(𝑘)=𝐺_1,4𝑥^(𝑘),𝑦^(𝑘)=𝑦^(𝑘)−𝑦^(𝑘),𝑥^(𝑘);𝐹+𝑦^(𝑘),𝑤^(𝑘)−𝑤;𝐹_(𝑘),𝑥^(𝑘);𝐹⁻¹𝐹𝑦^(𝑘),𝑥^(𝑘+1)=𝐺_1,7𝑥^(𝑘),𝑦^(𝑘),𝑧^(𝑘)=𝑧^(𝑘)−𝑧^(𝑘),𝑥^(𝑘);𝐹+𝑧^(𝑘),𝑦^(𝑘)−𝑦;𝐹_(𝑘),𝑥^(𝑘);𝐹⁻¹𝐹𝑧^(𝑘),𝑦^(𝑘)=𝐺_1,2𝑥^(𝑘),𝑧^(𝑘)=𝐺_2,4𝑥^(𝑘),𝑦^(𝑘)=𝑦^(𝑘)−𝑦^(𝑘),𝑥^(𝑘);𝐹⁻¹×𝑦^(𝑘),𝑥^(𝑘);𝐹−𝑦^(𝑘),𝑤^(𝑘)+𝑤;𝐹_(𝑘),𝑥^(𝑘);𝐹𝑦^(𝑘),𝑥^(𝑘);𝐹⁻¹𝐹𝑦^(𝑘),𝑥^(𝑘+1)=𝐺_2,7𝑥^(𝑘),𝑦^(𝑘),𝑧^(𝑘)=𝑧^(𝑘)−𝑧^(𝑘),𝑥^(𝑘);𝐹+𝑧^(𝑘),𝑦^(𝑘)−𝑦;𝐹_(𝑘),𝑥^(𝑘);𝐹⁻¹𝐹𝑧^(𝑘). (3)(4) Per iteration, both methods use four functions, five first-order divided differences, and three matrix inversions. The notable feature of these algorithms is their simple design which makes them easily implemented to solve systems of nonlinear equations. Here, the fourth-order method 𝐺_1,4(𝑥^(𝑘),𝑦^(𝑘)) is the generalization of the method proposed by Ren et al. in [5] and 𝐺_2,4(𝑥^(𝑘),𝑦^(𝑘)) is the generalization of the method by Liu et al. [6].

In this paper, our aim is to develop derivative free iterative methods that may satisfy the basic requirements of generating quality numerical algorithms, that is, the algorithms with (i) high convergence speed, (ii) minimum computational cost, and (iii) simple design. In this way, we here propose a derivative free family of sixth-order methods. The scheme is composed of three steps of which the first two steps consist of any derivative free fourth-order method with the base as the Traub-Steffensen iteration (2) whereas the third step is weighted Traub-Steffensen iteration. The algorithm of the present contribution is as simple as the methods (3) and (4) but with an additional advantage that it possesses high computational efficiency, especially when applied for solving large systems of equations.

The rest of the paper is summarized as follows. The sixth-order scheme with its convergence analysis is presented in Section 2. In Section 3, the computational efficiency of new methods is discussed and is compared with the methods which lie in the same category. Various numerical examples are considered in Section 4 to show the consistent convergence behavior of the methods and to verify the theoretical results. Section 5 contains the concluding remarks.

2. The Method and Its Convergence

Based on the above considerations of a quality numerical algorithm, we begin with the following three-step iteration scheme: 𝑦^(𝑘)=𝐺_1,2𝑥^(𝑘),𝑧^(𝑘)=𝐺₄𝑥^(𝑘),𝑦^(𝑘),𝑥^(𝑘+1)=𝐺₆𝑥^(𝑘),𝑦^(𝑘),𝑧^(𝑘)=𝑧^(𝑘)−𝑎𝐼+𝑤^(𝑘),𝑥^(𝑘);𝐹⁻¹×𝑏𝑦^(𝑘),𝑥^(𝑘);𝐹+𝑐𝑦^(𝑘),𝑤^(𝑘)×𝑤;𝐹_(𝑘),𝑥^(𝑘);𝐹⁻¹𝐹𝑧^(𝑘), (5) where 𝐺₄(𝑥^(𝑘),𝑦^(𝑘)) denotes any derivative free fourth-order scheme and 𝑎, 𝑏, 𝑐 are some parameters to be determined.

In order to find the convergence order of scheme (5), we first define divided difference operator for multivariable function 𝐹 (see [15]). The divided difference operator of 𝐹 is a mapping [⋅,⋅;𝐹]∶𝐷×𝐷⊂𝑅^𝑛×𝑅^𝑛→𝐿(𝑅^𝑛) defined by [𝑥+ℎ,𝑥;𝐹]=¹₀𝐹^′(𝑥+𝑡ℎ)𝑑𝑡,∀𝑥,ℎ∈𝑅^𝑛. (6) Expanding 𝐹^′(𝑥+𝑡ℎ) in Taylor series at the point 𝑥 and then integrating, we have [𝑥+ℎ,𝑥;𝐹]=¹₀𝐹^′(𝑥+𝑡ℎ)𝑑𝑡=𝐹^′(𝑥)+12𝐹^′′(𝑥)1ℎ+6𝐹^′′′(𝑥)ℎ²+𝑂ℎ³, (7) where ℎ^𝑖=(ℎ,ℎ,…^𝑖,ℎ),ℎ∈𝑅^𝑛.

Let 𝑒^(𝑘)=𝑥^(𝑘)−𝛼. Assuming that Γ=[𝐹^′(𝛼)]⁻¹ exists and then developing 𝐹(𝑥^(𝑘)) and its first three derivatives in a neighborhood of 𝛼, we have 𝐹𝑥^(𝑘)=𝐹^′(𝛼)𝑒^(𝑘)+𝐴₂𝑒^(𝑘)2+𝐴₃𝑒^(𝑘)3+𝐴₄𝑒^(𝑘)4+𝐴₅𝑒^(𝑘)5+𝑂𝑒^(𝑘)6,𝐹^′𝑥^(𝑘)=𝐹^′(𝛼)𝐼+2𝐴₂𝑒^(𝑘)+3𝐴₃𝑒^(𝑘)2+4𝐴₄𝑒^(𝑘)3+5𝐴₅𝑒^(𝑘)4+𝑂𝑒^(𝑘)5,𝐹^′′𝑥^(𝑘)=𝐹^′(𝛼)2𝐴₂+6𝐴₃𝑒^(𝑘)+12𝐴₄𝑒^(𝑘)2+20𝐴₅𝑒^(𝑘)3+𝑂𝑒^(𝑘)4,𝐹^′′′𝑥^(𝑘)=𝐹^′(𝛼)6𝐴₃+24𝐴₄𝑒^(𝑘)+60𝐴₅𝑒^(𝑘)2+𝑂𝑒^(𝑘)3, (8)(9)(10)(11) where 𝐴_𝑖=(1/𝑖!)Γ𝐹^(𝑖)(𝛼)∈𝐿_𝑖(𝑅^𝑛,𝑅^𝑛) and (𝑒^(𝑘))^𝑖=(𝑒^(𝑘),𝑒^(𝑘),…^𝑖,𝑒^(𝑘)),𝑒^(𝑘)∈𝑅^𝑛. 𝐴_𝑖 are symmetric operators that are used later on.

We can now analyze behavior of the scheme (5) through the following theorem.

Theorem 1. Let the function 𝐹∶𝐷⊂𝑅^𝑛→𝑅^𝑛 be sufficiently differentiable in an open neighborhood 𝐷 of its zero 𝛼 and 𝐺₄(𝑥^(𝑘),𝑦^(𝑘)) is a fourth-order iteration function which satisfies 𝑒_𝑧^(𝑘)=𝑧^(𝑘)−𝛼=𝐵₀𝑒^(𝑘)4+𝑂𝑒^(𝑘)5, (12) where 𝐵₀∈𝐿₄(𝑅^𝑛,𝑅^𝑛) and 𝑒^(𝑘)=𝑥^(𝑘)−𝛼. If an initial approximation 𝑥⁽⁰⁾ is sufficiently close to 𝛼, then the local order of convergence of method (5) is at least 6 provided 𝑎=3, 𝑏=−1, and 𝑐=−1.

Proof. Let 𝑒_𝑤^(𝑘)=𝑤^(𝑘)−𝛼=𝑥^(𝑘)+𝛾𝐹(𝑥^(𝑘))−𝛼. Then, using (8), it follows that 𝑒_𝑤^(𝑘)=𝐼+𝛾𝐹^′(𝛼)𝑒^(𝑘)+𝛾𝐹^′(𝛼)𝐴₂𝑒^(𝑘)2+𝛾𝐹^′(𝛼)𝐴₃𝑒^(𝑘)3+𝑂𝑒^(𝑘)4. (13) Employing (7) for 𝑥+ℎ=𝑤^(𝑘), 𝑥=𝑥^(𝑘), and ℎ=𝑒_𝑤^(𝑘)−𝑒^(𝑘) and then using (9)–(11), we write 𝑤^(𝑘),𝑥^(𝑘);𝐹=𝐹^′(𝛼)×𝐼+𝐴₂𝑒_𝑤^(𝑘)+𝑒^(𝑘)+𝐴₃𝑒_𝑤^(𝑘)2+𝑒^(𝑘)2+𝑒_𝑤^(𝑘)𝑒^(𝑘)+𝐴₄𝑒_𝑤^(𝑘)3+𝑒^(𝑘)3+𝑒_𝑤^(𝑘)2𝑒^(𝑘)+𝑒_𝑤^(𝑘)𝑒^(𝑘)2+𝑂𝑒^(𝑘)4. (14) Expanding in the formal power developments, the inverse of preceding operator can be written as 𝑤^(𝑘),𝑥^(𝑘);𝐹⁻¹=𝐼−𝐴₂𝑒_𝑤^(𝑘)+𝑒^(𝑘)+𝐴²₂−𝐴₃×𝑒_𝑤^(𝑘)2+𝑒^(𝑘)2+2𝐴²₂−𝐴₃𝑒_𝑤^(𝑘)𝑒^(𝑘)−𝐴³₂−𝐴₃𝐴₂−𝐴₂𝐴₃+𝐴₄×𝑒_𝑤^(𝑘)3+𝑒^(𝑘)3−3𝐴³₂−2𝐴₃𝐴₂−2𝐴₂𝐴₃+𝐴₄×𝑒_𝑤^(𝑘)2𝑒^(𝑘)+𝑒_𝑤^(𝑘)𝑒^(𝑘)2+𝑂𝑒^(𝑘)4Γ. (15) Using (8) and (15) to the required terms in the first step of (5), we find that 𝑒_𝑦^(𝑘)=𝑦^(𝑘)−𝛼=𝐴₂𝑒_𝑤^(𝑘)𝑒^(𝑘)−𝐴²₂−𝐴₃×𝑒_𝑤^(𝑘)2𝑒^(𝑘)+𝑒_𝑤^(𝑘)𝑒^(𝑘)2+𝑂𝑒^(𝑘)4. (16) Equation (7) for 𝑥+ℎ=𝑦^(𝑘), 𝑥=𝑥^(𝑘), and ℎ=𝑒_𝑦^(𝑘)−𝑒^(𝑘) yields 𝑦^(𝑘),𝑥^(𝑘);𝐹=𝐹^′(𝛼)𝐼+𝐴₂𝑒_𝑦^(𝑘)+𝑒^(𝑘)+𝐴₃𝑒_𝑦^(𝑘)𝑒^(𝑘)+𝑒^(𝑘)2+𝐴₄𝑒^(𝑘)3+𝑂𝑒^(𝑘)4. (17) Similarly, substituting 𝑥+ℎ=𝑦^(𝑘), 𝑥=𝑤^(𝑘), and ℎ=𝑒_𝑦^(𝑘)−𝑒_𝑤^(𝑘) in (7), we obtain 𝑦^(𝑘),𝑤^(𝑘);𝐹=𝐹^′(𝛼)𝐼+𝐴₂𝑒_𝑦^(𝑘)+𝑒_𝑤^(𝑘)+𝐴₃𝑒_𝑦^(𝑘)𝑒_𝑤^(𝑘)+𝑒_𝑤^(𝑘)2+𝐴₄𝑒^(𝑘)_w³+𝑂𝑒^(𝑘)4. (18) Then, using (15), (17), and (18) to the required terms, we find that 𝑎𝐼+𝑤^(𝑘),𝑥^(𝑘);𝐹⁻¹𝑏𝑦^(𝑘),𝑥^(𝑘);𝐹+𝑐𝑦^(𝑘),𝑤^(𝑘)=();𝐹𝑎+𝑏+𝑐𝐼−𝐴₂𝑏𝑒_𝑤^(𝑘)+𝑐𝑒^(𝑘)+𝐴²₂−𝐴₃𝑏𝑒_𝑤^(𝑘)2+𝑐𝑒^(𝑘)2+(𝑏+𝑐)𝐴²₂−𝐴₃𝑒_𝑤^(𝑘)𝑒^(𝑘)+(𝑏+𝑐)𝐴₂𝑒_𝑦^(𝑘)+𝑂𝑒^(𝑘)3. (19) Postmultiplying (19) by [𝑤^(𝑘),𝑥^(𝑘);𝐹]⁻¹ and simplifying 𝜃=𝑎𝐼+𝑤^(𝑘),𝑥^(𝑘);𝐹⁻¹×𝑏𝑦^(𝑘),𝑥^(𝑘);𝐹+𝑐𝑦^(𝑘),𝑤^(𝑘)×𝑤;𝐹_(𝑘),𝑥^(𝑘);𝐹⁻¹=(𝑎+𝑏+𝑐)𝐼−(𝑎+2𝑏+𝑐)𝐴₂𝑒_𝑤^(𝑘)−(𝑎+𝑏+2𝑐)𝐴₂𝑒^(𝑘)+(𝑎+3𝑏+𝑐)𝐴²₂−(𝑎+2𝑏+𝑐)𝐴₃𝑒_𝑤^(𝑘)2+(𝑎+𝑏+3𝑐)𝐴²₂−(𝑎+𝑏+2𝑐)𝐴₃𝑒^(𝑘)2+(𝑎+2(𝑏+𝑐))2𝐴²₂−𝐴₃𝑒_𝑤^(𝑘)𝑒^(𝑘)+(𝑏+𝑐)𝐴₂𝑒_𝑦^(𝑘)+𝑂𝑒^(𝑘)3Γ, (20) Taylor series of 𝐹(𝑧^(𝑘)) about 𝛼 yields 𝐹𝑧^(𝑘)=𝐹^′(𝛼)𝑒_𝑧^(𝑘)+𝐴₂𝑒_𝑧^(𝑘)2+𝑂𝑒_𝑧^(𝑘)3. (21) Then, using (20) and (21) in the third step of (5), we obtain the error equation as 𝑒^(𝑘+1)=𝑒_𝑧^(𝑘)−𝜃𝐹^′(𝛼)𝑒_𝑧^(𝑘)+𝑂𝑒_𝑧^(𝑘)2=−(𝑎+𝑏+𝑐−1)𝑒_𝑧^(𝑘)+(𝑎+2𝑏+𝑐)𝐴₂𝑒_𝑤^(𝑘)𝑒_𝑧^(𝑘)+(𝑎+𝑏+2𝑐)𝐴₂𝑒_𝑧^(𝑘)𝑒^(𝑘)−(𝑎+3𝑏+𝑐)𝐴²₂−(𝑎+2𝑏+𝑐)𝐴₃𝑒_𝑤^(𝑘)2𝑒_𝑧^(𝑘)−(𝑎+𝑏+3𝑐)𝐴²₂−(𝑎+𝑏+2𝑐)𝐴₃×𝑒_𝑧^(𝑘)𝑒^(𝑘)2−((×𝑎+2𝑏+𝑐))2𝐴²₂−𝐴₃𝑒_𝑤^(𝑘)𝑒_𝑧^(𝑘)𝑒^(𝑘)−(𝑏+𝑐)𝐴₂𝑒_𝑦^(𝑘)𝑒_𝑧^(𝑘)+𝑂𝑒^(𝑘)7. (22) In order to find 𝑎, 𝑏, and 𝑐, it will be sufficient to equate the factors 𝑎+𝑏+𝑐−1, 𝑎+2𝑏+𝑐, and 𝑎+𝑏+2𝑐 to 0 and then solving the resulting system of equations 𝑎+𝑏+𝑐=1,𝑎+2𝑏+𝑐=0,𝑎+𝑏+2𝑐=0, (23) we obtain 𝑎=3, 𝑏=−1, and 𝑐=−1.
Thus, for this set of values, the above error equation reduces to 𝑒^(𝑘+1)=𝐴²₂𝑒_𝑤^(𝑘)+𝑒^(𝑘)2𝑒_𝑧^(𝑘)+2𝐴₂𝑒_𝑦^(𝑘)𝑒_𝑧^(𝑘)−𝐴₃𝑒_𝑤^(𝑘)𝑒_𝑧^(𝑘)𝑒^(𝑘)+𝑂𝑒^(𝑘)7=𝐵₀𝐼+𝛾𝐹^′(𝛼)6𝐴²₂−𝐴₃+𝛾²𝐹^′(𝛼)²𝐴²₂𝑒^(𝑘)6+𝑂𝑒^(𝑘)7, (24) which shows the sixth order of convergence and hence the result follows.

Finally, the sixth-order family of methods is expressed by 𝑦^(𝑘)=𝐺_1,2𝑥^(𝑘),𝑧^(𝑘)=𝐺₄𝑥^(𝑘),𝑦^(𝑘),𝑥^(𝑘+1)=𝐺₆𝑥^(𝑘),𝑦^(𝑘),𝑧^(𝑘)=𝑧^(𝑘)−M𝑥^(𝑘),𝑤^(𝑘),𝑦^(𝑘)𝐹𝑧^(𝑘), (25) wherein M𝑥^(𝑘),𝑤^(𝑘),𝑦^(𝑘)=3𝐼−𝑤^(𝑘),𝑥^(𝑘);𝐹⁻¹𝑦^(𝑘),𝑥^(𝑘)+𝑦;𝐹_(𝑘),𝑤^(𝑘);𝐹𝑤^(𝑘),𝑥^(𝑘);𝐹⁻¹. (26)

Thus, the scheme (25) defines a new three-step family of derivative free sixth-order methods with the first two steps as any fourth-order scheme whose base is the Traub-Steffensen method (3). Some simple members of this family are as follows.

Method-I. The first method, which is denoted by 𝐺_1,6, is given by 𝑦^(𝑘)=𝐺_1,2𝑥^(𝑘),𝑧^(𝑘)=𝐺_1,4𝑥^(𝑘),𝑦^(𝑘),𝑥^(𝑘+1)=𝑧^(𝑘)−M𝑥^(𝑘),𝑤^(𝑘),𝑦^(𝑘)𝐹𝑧^(𝑘), (27) where 𝐺_1,4(𝑥^(𝑘),𝑦^(𝑘)) is the fourth-order method as given in the formula (3). It is clear that this formula uses four functions, three first-order divided differences, and two matrix inversions per iteration.

Method-II. The second method, that we denote by 𝐺_2,6, is given as 𝑦^(𝑘)=𝐺_1,2𝑥^(𝑘),𝑧^(𝑘)=𝐺_2,4𝑥^(𝑘),𝑦^(𝑘),𝑥^(𝑘+1)=𝑧^(𝑘)−M𝑥^(𝑘),𝑤^(𝑘),𝑦^(𝑘)𝐹𝑧^(𝑘), (28) where 𝐺_2,4(𝑥^(𝑘),𝑦^(𝑘)) is the fourth-order method shown in (4). This method also requires the same evaluations as in the above method.

3. Computational Efficiency

Here, we estimate the computational efficiency of the proposed methods and compare it with the existing methods. To do this, we will make use of efficiency index, according to which the efficiency of an iterative method is given by 𝐸=𝜌^1/𝐶, where 𝜌 is the order of convergence and 𝐶 is the computational cost per iteration. For a system of 𝑛 nonlinear equations in 𝑛 unknowns, the computational cost per iteration is given by (see [16]) 𝐶𝜇,𝑛,𝑙=𝑃₀(𝑛)𝜇+𝑃(𝑛,𝑙). (29) Here, 𝑃₀(𝑛) denotes the number of evaluations of scalar functions used in the evaluation of 𝐹 and [𝑥,𝑦;𝐹], and 𝑃(𝑛,𝑙) denotes the number of products needed per iteration. The divided difference [𝑥,𝑦;𝐹] of 𝐹 is an 𝑛×𝑛 matrix with elements given as (see [17, 18]) 𝑦,𝑥;𝐹_𝑖𝑗=𝑓_𝑖𝑦₁,…,𝑦_𝑗−1,𝑦_𝑗,𝑥_𝑗+1,…,𝑥_𝑛−𝑓_𝑖𝑦₁,…,𝑦_𝑗−1,𝑥_𝑗,𝑥_𝑗+1,…,𝑥_𝑛+𝑓_𝑖𝑥₁,…,𝑥_𝑗−1,𝑦_𝑗,𝑦_𝑗+1,…,𝑦_𝑛−𝑓_𝑖𝑥₁,…,𝑥_𝑗−1,𝑥_𝑗,𝑦_𝑗+1,…,𝑦_𝑛×2𝑦_𝑗−𝑥_𝑗⁻¹,1⩽𝑖,𝑗⩽𝑛. (30) In order to express the value of 𝐶(𝜇,𝑛,𝑙) in terms of products, a ratio 𝜇>0 between products and evaluations of scalar functions and a ratio 𝑙⩾1 between products and quotients are required.

To compute 𝐹 in any iterative function, we evaluate 𝑛 scalar functions (𝑓₁,𝑓₂,…,𝑓_𝑛), and if we compute a divided difference [𝑥,𝑦;𝐹], then we evaluate 2𝑛(𝑛−1) scalar functions, where 𝐹(𝑥) and 𝐹(𝑦) are computed separately. We must add 𝑛² quotients from any divided difference, 𝑛² products for multiplication of a matrix with a vector or of a matrix by a scalar, and 𝑛 products for multiplication of a vector by a scalar. In order to compute an inverse linear operator, we solve a linear system, where we have 𝑛(𝑛−1)(2𝑛−1)/6 products and 𝑛(𝑛−1)/2 quotients in the LU decomposition and 𝑛(𝑛−1) products and 𝑛 quotients in the resolution of two triangular linear systems.

The computational efficiency of the present sixth-order methods 𝐺_1,6 and 𝐺_2,6 is compared with the existing fourth-order methods 𝐺_1,4 and 𝐺_2,4 and with the seventh-order methods 𝐺_1,7 and 𝐺_2,7. In addition, we also compare the present methods with each other. Let us denote efficiency indices of 𝐺_𝑖,𝑝 by 𝐸_𝑖,𝑝 and computational cost by 𝐶_𝑖,𝑝. Then, taking into account the above and previous considerations, we have 𝐶_1,4=6𝑛²𝑛−3𝑛𝜇+32𝑛²(),𝐸+3𝑛−5+3𝑙4𝑛+1_1,4=4^1/𝐶1,4,𝐶_2,4=6𝑛²𝑛−3𝑛𝜇+32𝑛²(),𝐸+9𝑛−8+3𝑙4𝑛+2_2,4=4^1/𝐶2,4,𝐶_1,6=6𝑛²𝑛−2𝑛𝜇+32𝑛²(),𝐸+12𝑛−8+3𝑙4𝑛+3_1,6=6^1/𝐶1,6,𝐶_2,6=6𝑛²𝑛−2𝑛𝜇+32𝑛²(),𝐸+18𝑛−11+12𝑙𝑛+1_2,6=6^1/𝐶2,6,𝐶_1,7=10𝑛²𝑛−6𝑛𝜇+22𝑛²(),𝐸+3𝑛−5+𝑙13𝑛+3_1,7=7^1/𝐶1,7,𝐶_2,7=10𝑛²𝑛−6𝑛𝜇+22𝑛²(),𝐸+7𝑛−7+𝑙13𝑛+5_2,7=7^1/𝐶2,7. (31)(32)(33)(34)(35)(36)

3.1. Comparison between Efficiencies

To compare the computational efficiencies of the iterative methods, say 𝐺_𝑖,𝑝 against 𝐺_𝑗,𝑞, we consider the ratio 𝑅_{𝑖,𝑝;𝑗,𝑞}=log𝐸_𝑖,𝑝log𝐸_𝑗,𝑞=𝐶_𝑗,𝑞log𝑝𝐶_𝑖,𝑝𝑞.log (37) It is clear that if 𝑅_{𝑖,𝑝;𝑗,𝑞}>1, the iterative method 𝐺_𝑖,𝑝 is more efficient than 𝐺_𝑗,𝑞.

𝐺_1,4 versus 𝐺_1,6 Case. For this case, the ratio (37) is given by 𝑅_1,6;1,4=log6log42𝑛²+𝑛18𝜇+12𝑙+3−9𝜇+3𝑙−52𝑛²,+𝑛18𝜇+12𝑙+12−6𝜇+9𝑙−8 (38) which shows that 𝑅_1,6;1,4>1 for 𝜇>0, 𝑙⩾1, and 𝑛⩾9. Thus, we have 𝐸_1,6>𝐸_1,4 for 𝜇>0, 𝑙⩾1, and 𝑛⩾9.

𝐺_2,4 versus 𝐺_1,6 Case. In this case, the ratio (37) takes the following form 𝑅_1,6;2,4=log6log42𝑛²+𝑛18𝜇+12𝑙+9−9𝜇+6𝑙−82𝑛².+𝑛18𝜇+12𝑙+12−6𝜇+9𝑙−8 (39) In this case, it is easy to prove that 𝑅_1,6;2,4>1 for 𝜇>0, 𝑙⩾1, and 𝑛⩾2, which implies that 𝐸_1,6>𝐸_2,4.

𝐺_1,4 versus 𝐺_2,6 Case. The ratio (37) yields 𝑅_2,6;1,4=log6log42𝑛²+𝑛18𝜇+12𝑙+3−9𝜇+3𝑙−52𝑛²,+𝑛18𝜇+12𝑙+18−6𝜇+12𝑙−11 (40) which shows that 𝑅_2,6;1,4>1 for 𝜇>0, 𝑙⩾1, and 𝑛⩾19. Thus, we conclude that 𝐸_2,6>𝐸_1,4 for 𝜇>0, 𝑙⩾1, and 𝑛⩾19.

𝐺_2,4 versus 𝐺_2,6 Case. In this case, the ratio (37) is given by 𝑅_2,6;2,4=log6log42𝑛²+𝑛18𝜇+12𝑙+9−9𝜇+6𝑙−82𝑛².+𝑛18𝜇+12𝑙+18−6𝜇+12𝑙−11 (41) With the same range of 𝜇, 𝑙 as in the previous case and 𝑛⩾6, the ratio 𝑅_2,6;2,4>1, which implies that 𝐸_2,6>𝐸_2,4.

𝐺_1,6 versus 𝐺_2,6 Case. In this case, it is enough to compare the corresponding values of 𝐶_1,6 and 𝐶_2,6 from (33) and (34). Thus, we find that 𝐸_1,6>𝐸_2,6 for 𝜇>0, 𝑙⩾1, and 𝑛⩾2.

𝐺_1,6 versus 𝐺_1,7 Case. In this case, the ratio (37) is given by 𝑅_1,6;1,7=log63log72𝑛²+𝑛20𝜇+13𝑙+3−12𝜇+3𝑙−522𝑛².+𝑛18𝜇+12𝑙+12−6𝜇+9𝑙−8 (42) It is easy to show that 𝑅_1,6;1,7>1 for 𝜇>0, 𝑙⩾1, and 𝑛⩾4, which implies that 𝐸_1,6>𝐸_1,7 for this range of values of the parameters (𝜇,𝑛,𝑙).

𝐺_1,6 versus 𝐺_2,7 Case. The ratio (37) is given by 𝑅_1,6;2,7=log63log72𝑛²+𝑛20𝜇+13𝑙+7−12𝜇+5𝑙−722𝑛².+𝑛18𝜇+12𝑙+12−6𝜇+9𝑙−8 (43) With the same range of 𝜇, 𝑙 as in the previous cases and 𝑛⩾2, we have 𝑅_1,6;2,7>1, which implies that 𝐸_1,6>𝐸_2,7.

𝐺_2,6 versus 𝐺_1,7 Case. The ratio (37) yields 𝑅_2,6;1,7=log63log72𝑛²+𝑛20𝜇+13𝑙+3−12𝜇+3𝑙−522𝑛²,+𝑛18𝜇+12𝑙+18−6𝜇+12𝑙−11 (44) which shows that 𝑅_2,6;1,7>1 for 𝜇>0, 𝑙⩾1, and 𝑛⩾11, so it follows that 𝐸_2,6>𝐸_1,7.

𝐺_2,6 versus 𝐺_2,7 Case. For this case, the ratio (37) is given by 𝑅_2,6;2,7=log63log72𝑛²+𝑛20𝜇+13𝑙+7−12𝜇+5𝑙−722𝑛².+𝑛18𝜇+12𝑙+18−6𝜇+12𝑙−11 (45) In this case, also it is not difficult to prove that 𝑅_2,6;2,7>1 for 𝜇>0, 𝑙⩾1, and 𝑛⩾5, which implies that 𝐸_2,6>𝐸_2,7.

We summarize the above results in the following theorem.

Theorem 2. For 𝜇>0 and 𝑙⩾1, we have the following: (i)𝐸_1,6>𝐸_1,4 for 𝑛⩾9;(ii)𝐸_2,6>𝐸_1,4 for 𝑛⩾19;(iii)𝐸_2,6>𝐸_2,4 for 𝑛⩾6;(iv)𝐸_1,6>𝐸_1,7 for 𝑛⩾4;(v)𝐸_2,6>𝐸_1,7 for 𝑛⩾11;(vi)𝐸_2,6>𝐸_2,7 for 𝑛⩾5;(vii){𝐸_1,6>𝐸_2,4,𝐸_1,6>𝐸_2,6,𝐸_1,6>𝐸_2,7} for 𝑛⩾2. Otherwise, the comparison depends on 𝜇, 𝑙, and 𝑛.

4. Numerical Results

In this section, some numerical problems are considered to illustrate the convergence behavior and computational efficiency of the proposed methods. The performance is compared with the existing methods 𝐺_1,4, 𝐺_2,4, 𝐺_1,7, and 𝐺_2,7. All computations are performed using the programming package Mathematica [19] using multiple-precision arithmetic with 4096 digits. For every method, we analyze the number of iterations (𝑘) needed to converge to the solution such that ‖𝑥^(𝑘+1)−𝑥^(𝑘)‖+‖𝐹(𝑥^(𝑘))‖<10⁻²⁰⁰. In numerical results, we also include the CPU time utilized in the execution of program which is computed by the Mathematica command “TimeUsed[]”. In order to verify the theoretical order of convergence, we calculate the computational order of convergence (𝜌_𝑐) using the formula 𝜌_𝑐=log𝐹𝑥^(𝑘)/𝐹𝑥^(𝑘−1)log𝐹𝑥^(𝑘−1)/𝐹𝑥^(𝑘−2), (46) (see [20]) taking into consideration the last three approximations in the iterative process.

To connect the analysis of computational efficiency with numerical examples, the definition of the computational cost (29) is applied, according to which an estimation of the factor 𝜇 is claimed. For this, we express the cost of the evaluation of the elementary functions in terms of products, which depends on the computer, the software, and the arithmetics used (see [21, 22]). In Table 1, the elapsed CPU time (measured in milliseconds) in the computation of elementary functions and an estimation of the cost of the elementary functions in product units are displayed. The programs are performed in the processor, Intel (R) Core (TM) i5-480M CPU @ 2.67 GHz (64-bit Machine) Microsoft Windows 7 Home Basic 2009, and are compiled by Mathematica 7.0 using multiple-precision arithmetic. It can be observed from Table 1 that, for this hardware and the software, the computational cost of quotient with respect to product is 𝑙=2.8.

Table 1 

CPU time and estimation of computational cost of the elementary functions, where 𝑥=√3−1 and 𝑦=√5.

The present methods 𝐺_1,6 and 𝐺_2,6 are tested by using the values −0.01, 0.01, and 0.5 for the parameter 𝛾. The following problems are chosen for numerical tests.

Problem 1. Considering the system of two equations, 𝑥²₁−𝑥₂𝑥+1=0,₁−cos𝜋𝑥₂2=0. (47) In this problem, (𝑛,𝜇)=(2,52) are the values used in (31)–(36) for calculating computational costs and efficiency indices of the considered methods. The initial approximation chosen is 𝑥⁽⁰⁾={0.25,0.5}^𝑇 and the solution is 𝛼={0,1}^𝑇.

Problem 2. Consider the system of three equations: 10𝑥₁+sin𝑥₁+𝑥₂−1=0,8𝑥₂−cos²𝑥₃−𝑥₂−1=0,12𝑥₃+sin𝑥₃−1=0, (48) with initial value 𝑥⁽⁰⁾={0.8,0.5,0.125}^𝑇 towards the solution {}𝛼=0.0689783491726666…,0.2464424186091830…,0.0769289119875370…^𝑇. (49) For this problem, (𝑛,𝜇)=(3,103.33)are used in (31)–(36) to calculate computational costs and efficiency indices.

Problem 3. Next, consider the following boundary value problem (see [23]): 𝑦^′′+𝑦³=0,𝑦(0)=0,𝑦(1)=1. (50)Assume the following partitioning of the interval [0,1]: 𝑢₀=0<𝑢₁<𝑢₂<⋯<𝑢_𝑚−1<𝑢_𝑚𝑢=1,_𝑗+1=𝑢_𝑗1+ℎ,ℎ=𝑚. (51) Let us define 𝑦₀=𝑦(𝑢₀)=0,𝑦₁=𝑦(𝑢₁),…,𝑦_𝑚−1=𝑦(𝑢_𝑚−1),𝑦_𝑚=𝑦(𝑢_𝑚)=1. If we discretize the problem by using the numerical formula for second derivative, 𝑦_𝑘^′′=𝑦_𝑘−1−2𝑦_𝑘+𝑦_𝑘+1ℎ²,𝑘=1,2,3,…,𝑚−1, (52) we obtain a system of 𝑚−1 nonlinear equations in 𝑚−1 variables: 𝑦_𝑘−1−2𝑦_𝑘+𝑦_𝑘+1+ℎ²𝑦³_𝑘=0,𝑘=1,2,3,…,𝑚−1. (53) In particular, we solve this problem for 𝑚=5 so that 𝑛=4 by selecting 𝑦⁽⁰⁾={0.5,0.5,0.5,0.5}^𝑇 as the initial value. The solution of this problem is {}𝛼=0.21054188948074775…,0.42071046387616439…,0.62790045371805633…,0.82518822786851363…^𝑇 (54) and concrete values of the parameters are (𝑛,𝜇)=(4,4).

Problem 4. Consider the system of fifteen equations (see [16]): ¹⁵_{𝑗=1,𝑗≠𝑖}𝑥_𝑗−𝑒^−𝑥𝑖=0,1≤𝑖≤15. (55) In this problem, the concrete values of the parameters (𝑛,𝜇) are (15,81). The initial approximation assumed is 𝑥⁽⁰⁾={1,1,…,1}^𝑇 and the solution of this problem is {}𝛼=0.066812203179582582…,0.066812203179582582…,…,0.066812203179582582…^𝑇. (56)

Problem 5. Consider the system of fifty equations: 𝑥²_𝑖𝑥_𝑖+1()𝑥−1=0,𝑖=1,2,…,49₂₅₀𝑥₁−1=0. (57) In this problem, (𝑛,𝜇)=(50,2) are the values used in (31)–(36) for calculating computational costs and efficiency indices. The initial approximation assumed is 𝑥⁽⁰⁾={1.5,1.5,1.5,…,1.5}^𝑇 for obtaining the solution 𝛼={1,1,1,…,1}^𝑇.

Problem 6. Lastly, consider the nonlinear and nondifferentiable integral equation of mixed Hammerstein type (see [24]): 𝑥(𝑠)1=1+2¹₀𝐺(𝑠,𝑡)|𝑥(𝑡)|+𝑥(𝑡)²𝑑𝑡, (58) where 𝑥∈𝐶[0,1]; 𝑠,𝑡∈[0,1] and the kernel 𝐺 is given as follows: 𝐺()=()𝑠()𝑠,𝑡1−𝑠𝑡,𝑡⩽𝑠,1−𝑡,𝑠⩽𝑡. (59) We transform the above equation into a finite-dimensional problem by using Gauss-Legendre quadrature formula given as ¹₀𝑓(𝑡)𝑑𝑡≈^𝑛_𝑗=1𝜛_𝑗𝑓𝑡_𝑗, (60) where the abscissas 𝑡_𝑗 and the weights 𝜛_𝑗 are determined for 𝑛=8 by Gauss-Legendre quadrature formula. Denoting the approximation of 𝑥(𝑡_𝑖) by 𝑥_𝑖(𝑖=1,2,…,8), we obtain the system of nonlinear equations: 2𝑥_𝑖−2−⁸_𝑗=1𝑎_𝑖𝑗𝑥_𝑗+𝑥²_𝑗=0, (61) where 𝑎_𝑖𝑗=𝜛_𝑗𝑡_𝑗1−𝑡_𝑖𝜛if𝑗⩽𝑖,_𝑗𝑡_𝑖1−𝑡_𝑗if𝑖<𝑗. (62) In this problem, (𝑛,𝜇)=(8,10) are the values used in (31)–(36) for calculating computational costs and efficiency indices. The initial approximation assumed is 𝑥⁽⁰⁾={1,1,1,…,1}^𝑇 for obtaining the solution {}𝛼=1.0115010875012980…,1.0546781130247093…,1.1098925075633296…,1.1481439774759322,…,1.1481439774759322…,1.1098925075633296…,1.0546781130247093…,1.0115010875012980…^𝑇. (63)

Table 2 shows the numerical results obtained for the considered problems by various methods. Displayed in this table are the number of iterations (𝑘), the computational order of convergence (𝜌_𝑐), the computational costs (𝐶_𝑖,𝑝) in terms of products, the computational efficiencies (𝐸_𝑖,𝑝), and the mean elapsed CPU time (e-time). Computational cost and efficiency are calculated according to the corresponding expressions given by (31)–(36) by using the values of parameters 𝑛 and 𝜇 as calculated in each problem, while taking 𝑙=2.8 in each case. The mean elapsed CPU time is calculated by taking the mean of 50 performances of the program, where we use ‖𝑥^(𝑘+1)−𝑥^(𝑘)‖+‖𝐹(𝑥^(𝑘))‖<10⁻²⁰⁰ as the stopping criterion in single performance of the program.

Table 2 

Comparison of the performance of methods.

From the numerical results, we can observe that like the existing methods the present methods show consistent convergence behavior. It is seen that some methods do not preserve the theoretical order of convergence, especially when applied for solving some typical type of nonlinear systems. This can be observed in the last problem of nondifferentiable mixed Hammerstein integral equation where the seventh-order methods 𝐺_1,7 and 𝐺_2,7 yield the eighth order of convergence. However, for the present methods, the computational order of convergence overwhelmingly supports the theoretical order of convergence. Comparison of the numerical values of computational efficiencies exhibited in the second last column of Table 2 verifies the theoretical results of Theorem 2. As we know, the computational efficiency is proportional to the reciprocal value of the total CPU time necessary to complete running iterative process. This means that the method with high efficiency utilizes less CPU time than the method with low efficiency. The truthfulness of this fact can be judged from the numerical values of computational efficiency and elapsed CPU time displayed in the last two columns of Table 2, which are in complete agreement according to the notion.

5. Concluding Remarks

In the foregoing study, we have proposed iterative methods with the sixth order of convergence for solving systems of nonlinear equations. The schemes are totally derivative free and therefore particularly suited to those problems in which derivatives require lengthy computation. A development of first-order divided difference operator for functions of several variables and direct computation by Taylor's expansion are used to prove the local convergence order of new methods. Comparison of efficiencies of the new schemes with the existing schemes is shown. It is observed that the present methods have an edge over similar existing methods, especially when applied for solving large systems of equations. Six numerical examples have been presented and the relevant performances are compared with the existing methods. Computational results have confirmed the robust and efficient character of the proposed techniques. Similar numerical experimentations have been carried out for a number of problems and results are found to be on a par with those presented here.

We conclude the paper with the remark that in many numerical applications multiprecision in computations is required. The results of numerical experiments justify that the high-order efficient methods associated with a multiprecision arithmetic floating point are very useful, because they yield a clear reduction in the number of iterations to achieve the required solution.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.