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

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.


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  ( () ) ,  = 0, 1, 2, . . ., (1) where [  ()] −1 is the inverse of the first Fréchet derivative   () of the function ().
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 sixthorder 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.
We can now analyze behavior of the scheme (5) through the following theorem.
Method-I.The first method, which is denoted by  1,6 , is given by 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 where  2,4 ( () ,  () ) is the fourth-order method shown in (4).This method also requires the same evaluations as in the above method.

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])  (, , ) =  0 ()  +  (, ) . (29) Here,  0 () 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]) 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 ( 1 ,  2 , . . .,   ), and if we compute a divided difference [, ; ], then we evaluate 2( − 1) scalar functions, where () and () are computed separately.We must add  2 quotients from any divided difference,  2 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 fourthorder 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 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: Otherwise, the comparison depends on , , and .

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 −200 .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 (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  1 that, for this hardware and the software, the computational cost of quotient with respect to product is  = 2.8.
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 4. Consider the system of fifteen equations (see [16]): In this problem, the concrete values of the parameters (, ) are (15,81).The initial approximation assumed is  (0) = {1, 1, . . ., 1}  and the solution of this problem is  = {0.066812203179582582 . . ., 0.066812203179582582 . . ., . . ., 0.066812203179582582 . ..}  . ( Problem 5. Consider the system of fifty equations: In this problem, (, ) = (50, 2) are the values used in (31)-(36) for calculating computational costs and efficiency indices.The initial approximation assumed is  (0) = {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]): where  ∈ [0, 1]; ,  ∈ [0, 1] and the kernel  is given as follows: We transform the above equation into a finite-dimensional problem by using Gauss-Legendre quadrature formula given as 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: where 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 −200 as the stopping criterion in single performance of the program.
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.

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.

Table 1 :
CPU time and estimation of computational cost of the elementary functions, where  = √ 3 − 1 and  = √5.
5, 0.5}  as the initial value.The solution of this problem is

Table 2 :
Comparison of the performance of methods.