On a Reduced Cost Higher Order Traub-Steffensen-Like Method for Nonlinear Systems

: We propose a derivative-free iterative method with ﬁfth order of convergence for solving systems of nonlinear equations. The scheme is composed of three steps, of which the ﬁrst two steps are that of third order Traub-Steffensen-type method and the last is derivative-free modiﬁcation of Chebyshev’s method. Computational efﬁciency is examined and comparison between the efﬁciencies of presented technique with existing techniques is performed. It is proved that, in general, the new method is more efﬁcient. Numerical problems, including those resulting from practical problems viz. integral equations and boundary value problems, are considered to compare the performance of the proposed method with existing methods. Calculation of computational order of convergence shows that the order of convergence of the new method is preserved in all the numerical examples, which is not so in the case of some of the existing higher order methods. Moreover, the numerical results, including the CPU-time consumed in the execution of program, conﬁrm the accurate and efﬁcient behavior of the new technique.


Introduction
We are concerned with the problem of solving a system of nonlinear equations (1) This problem can precisely be stated as to find a solution vector α = (α 1 , α 2 , ..., α m ) T such that F(α) = 0, where F(x) : D ⊂ R m −→ R m is the given nonlinear vector function F(x) = ( f 1 (x), f 2 (x), ..., f m (x)) T and x = (x 1 , x 2 , ..., x m ) T . The vector α can be computed as a fixed point of some function M : D ⊂ R m → R m by means of fixed point iteration x (k+1) = M(x (k) ), k ≥ 0. ( Many applied problems in Science and Engineering are reduced to solve numerically the system F(x) = 0 of nonlinear equations (see, for example [1-6]). A plethora of iterative methods are developed in literature for solving such equations. A classical method is cubically convergent Chebyshev's method (see [7]) where L F (x (k) ) = F (x (k) ) −1 F (x (k) )F (x (k) ) −1 F(x k ). This one-point iterative scheme depends explicitly on the first two derivatives of F. In [7], Ezquerro and Hernández present modification in Chebyshev's method that avoids the computation of second derivative F while maintaining third-order of convergence. It has the following form: There is an interest in constructing derivative free iterative processes obtained by considering an approximation of the first derivative of F from a divided difference of first order. One class of such methods is called the class of Secant-type methods which is obtained by replacing F with the divided difference operator [x (k−1) , x (k) ; F]. Using this operator a family of derivative free methods is given in [8]. The authors call this family the Chebyshev-Secant-type method and it is defined as where a, b and c are non-negative parameters. Another class of derivative free methods is the class of Steffensen-type processes that replaces F with operator [w(x (k) ), x (k) ; F], wherein w : R m → R m . The work presented in [9] analyzes Steffensen-type iterative method which is given as (k) ), x (k) ; F] −1 F(x (k) ), (k) ), x (k) ; F] −1 b F(x (k) ) + c F(y (k) ) , k ≥ 0.
The two-step third order Traub-Steffensen-type method, i.e., the case of (6) for a = b = c = 1, can be written as x (0) ∈ D, w(x (k) ) = x (k) + βF(x (k) ), y (k) = M 2,1 (x (k) ), x (k+1) = M 3,1 (x (k) , y (k) ) = y (k) − [w(x (k) ), x (k) ; F] −1 F(y (k) ), k ≥ 0, where M 2,1 (x (k) ) = x (k) − [w(x (k) ), x (k) ; F] −1 F(x (k) ) is the quadratically convergent Traub-Steffensen scheme. Here and in the sequel, the symbol M p,i is used for denoting an i-th iteration function of convergence order p. It can be observed that the third order scheme (7) is computationally more efficient than quadratically convergent Traub-Steffensen scheme. The reason is that the convergence order is increased from two to three at the cost of only one function evaluation without adding extra inverse operator. We discuss computational efficiency in later sections. Researchers have always been trying to develop the iterative method with increasing efficiency since different methods converge to the solution with different convergence speed. This can be done either by increasing the convergence order or by decreasing the computational cost or both. In [11], Ren et al. have derived a fourth order derivative-free method that uses three F, three divided differences and two matrix inversions per iteration. Zheng et al. [12] have constructed two families of fourth order derivative-free methods for scalar nonlinear equations, that are extendable to solve systems of nonlinear equations. First family requires to evaluate three F, three divided differences and two matrix inversions, whereas the second family needs three F, three divided differences and three matrix inversions. Grau et al. presented a fourth order derivative-free method in [13] utilizing four F, two divided differences and two matrix inversions. Sharma and Arora [14] presented a fourth order derivative-free method that uses the evaluations of three F, three divided differences and one matrix inversion per each step.
In search of more fast techniques, researchers have also introduced sixth and seventh order derivative-free methods in [13,[15][16][17][18]. The sixth order method in [13] proposed by Grau et al. requires five F, two divided differences and two matrix inverses. Sharma and Arora [17] also developed a method of at least sixth order which requires evaluation of five functions, two divided difference and one matrix inversion per iteration. The seventh order method proposed by Sharma and Arora [15] utilizes four F, five divided differences and two matrix inversions per iteration. The seventh order methods presented by Wang and Zhang [16] use four F, five divided differences and three matrix inversions. Ahmad et al. [18] proposed eighth order derivative free method without memory which uses six functions evaluations, three divided difference and one matrix inversion.
The main goal in this study is to develop a derivative-free method of high computational efficiency, that means a method with high convergence speed and low computational cost. Consequently, we present a Traub-Steffensen-type method of fifth order of convergence which requires the evaluations four F, two divided differences and only one matrix inversion per step. The scheme of the present contribution is simple and consists of three steps. Of the three steps, the first two are that of cubically convergent Traub-Steffensen-type scheme (7) whereas the third is derivative-free modification of Chebyshev's scheme (3). We show that the proposed method is more efficient than existing methods of similar nature.
The content of the rest of the paper is summarized as follows. Basic definitions relevant to the present work are stated in Section 2. In Section 3, the scheme of fifth order method is introduced and its convergence behavior is studied. In Section 4, the computational efficiency of the new method is examined and also compared with the existing derivative-free methods. In Section 5, the basins of attractors are presented to check the stability and convergence of the new method. Numerical tests are performed in Section 6 to verify the theoretical results as proved in Sections 3 and 4. Section 7 contains the concluding remarks.

Computational Order of Convergence
Let α be a solution of the function F(x) = 0 and x (k−2) , x (k−1) , x (k) and x (k+1) be the four consecutive iterations close to α. Then, the computational order of convergence (say, p c ) can be calculated using the formula (see [19])

Divided Difference
Divided difference operator for multivariable function F (see [4,5,20]) is a mapping [·, · ; F] : If F is differentiable, we can also define first order divided difference as (see [4,20]) This also implies that It can be seen that the divided difference operator [x, y ; F] is an m × m matrix and the definitions (9) and (10) are equivalent (for details see [20]). For computational purpose the following definition (see [5]), is used

Computational Efficiency
Computational efficiency of an iterative method for solving F(x) = 0 is calculated by the efficiency index E = p 1/C , (for detail see [21,22]), where p is the order of convergence and C is the total cost of computation. The cost of computation C is measured in terms of the total number of function evaluations per iteration and the number of operations (that means products and quotients) per iteration.
Note that this is a scheme whose first two steps are that of third order Traub-Steffensen-type method (7) whereas third step is based on Chebyshev's method (3). The scheme requires first and second derivatives of F at y (k) . To make this a derivative-free method, we describe an approach as follows: Consider the Taylor expansion of F(z (k) ) about y (k) , Then, it follows that Using the fact that F(z (k) ) − F(y (k) ) = [z (k) , y (k) ; F](z (k) − y (k) ), (see, for example [4,5]), we can write (15) as Then, using the second step of (13) in the above equation, it follows that Let us assume F (y (k) ) ≈ [w (k) , x (k) ; F], then (17) implies In addition, we have that Using (18) in (19), we obtain that Now, we can write the third-step of (13) in modified form as Thus, we define the following new method: wherein Since the scheme (22) is composed of Traub-Steffensen like steps, we call it the Traub-Steffensen-like method.
In order to explore the convergence properties of Traub-Steffensen-like method, we recall some important results from the theory of iteration functions. First, we state the following well-known result (see [3,23]): then α is a point of attraction for the iteration x (k+1) = M(x (k) ), where ρ is a spectral radius of M (α).
Next, we state a result which has been proven in [24] by Madhu et al. and that shows α is a point of attraction for a general iteration function of the form M(x) = P(x) − Q(x)R(x).

Lemma 2.
Let F : D ⊂ R m → R m be sufficiently Fréchet differentiable at each point of an open convex set D of α ∈ D, which is a solution of the nonlinear system F(x) = 0. Suppose that P, Q, R : D ⊂ R m → R m are sufficiently Fréchet differentiable functions (depending on F) at each point in the set D with the properties P(α) = α, Q(α) = 0, R(α) = 0. Then, there exists a ball Let us also recall the definition (10) of divided difference operator. Then, expanding F (x + th) in (10) by Taylor series at the point x and thereafter integrating, we have that where Assuming that Γ = F (α) −1 exists, then expanding F(x (k) ) and its first three derivatives in a neighborhood of α by Taylor's series, we have that and where i−times · · · · , e (k) ), e (k) ∈ R m . We are in a situation to analyze the behavior of Traub-Steffensen-like method. Thus, the following theorem is proved: Assume that x ∈ S =S(α, ) and F (x) is continuous and nonsingular at α, and x (0) close to α. Then, α is a point of attraction of the sequence {x (k) } generated by the Traub-Steffensen-like method (22). Furthermore, the sequence so developed converges locally to α with order at least 5.
Proof. First we show that α is a point of attraction of Traub-Steffensen-like iteration. In this case, we have that so that ρ(M (α)) = 0 < 1 and by Lemma 1, α is a point of attraction of (22). Let e (k) . Then using (25), it follows that e (k) w − e (k) in Equation (24) and then using (26)- (29), we can write where λ = βF (α), Expansion of the inverse of preceding divided difference operator is given as By using (25) and (31) in the first step of method (22), we get e (k) Taylor expansion of F(y k ) about α yields, From the second step of (22), on using (31) and (33), it follows that By Taylor expansion of F(z k ) about α, Equation (24), for From (31) and (36), we have Equations (31) and (37) yield Applying Equations (34), (35) and (38) in the last step of method (22) and then simplifying, we get the error equation This completes the proof of Theorem 1.
Thus, the Traub-Steffensen-like method (22) defines a one-parameter (β) family of derivative-free fifth order methods. Now onwards we denote it by M 5,1 . In terms of computational cost M 5,1 utilizes four functions, two divided difference and one matrix inversion per each step. In the next section we will compare the computational efficiency of the new method with the existing derivative-free methods.

Computational Efficiency
In order to find the computational efficiency we will use the definition given in Section 2.3. The various evaluations and arithmetic operations that contribute towards the cost of computation are described as follows. For the computation of F in any iterative function we evaluate m scalar functions f i , (1 ≤ i ≤ m) and when computing a divided difference [x, y ; F] (see, Section 2.2) we evaluate m(m − 1) scalar functions, wherein F(x) and F(y) are evaluated separately. Furthermore, one has to add m 2 divisions from any divided difference. For the computation of inverse linear operator, a linear system can be solved that requires m(m − 1)(2m − 1)/6 products and m(m − 1)/2 divisions in the LU decomposition process, and m(m − 1) products and m divisions in the resolution of two triangular linear systems. Moreover, we add m products for the multiplication of a vector by a scalar and m 2 products for multiplying a matrix by a vector or of a matrix by a scalar.
Fifth order method by Kumar et al. (M 5,2 ): Sixth order method by Grau et al. (M 6,1 ): Wang-Zhang seventh order method (M 7,1 ): Sharma-Arora seventh order method (M 7,2 ): Let us denote efficiency indices of the methods M p,i by E p,i and their computational costs by C p,i . Then, using the definition of the Section 2.3 taking into account the above considerations of evaluations and operations, we have that m and E 2,1 = 2 1/C 2,1 .
To compare the efficiency of considered iterative methods, say M p,i against M q,j , we consider the ratio It is clear that when R p,i;q,j > 1, the iterative method M p,i is more efficient than M q,j . M 3,1 versus M 2,1 case: For this case the ratio (50) is given by It can be easily shown that R 3,1;2,1 > 1 for m ≥ 2. This implies that E 3,1 > E 2,1 for m ≥ 2. Thus, M 3,1 is more efficient than M 2,1 as we have stated in the introduction section. The ratio 50 is given by .
The above results are summarized in the following theorem: E 5,1 > E 6,1 f or m ≥ 8.

Complex Dynamics of Methods
Our aim is to analyze the complex dynamics of the new method based on graphical tool 'basins of attraction' of the zeros of polynomial P(z) in complex plane. Visual display of the basins gives important information about the stability and convergence of iterative methods. This idea was introduced initially by Vrscay and Gilbert [26]. In recent times, many authors have used this concept in their work, see, for example [27,28] and references therein. We consider the method (22) to analyze the basins of attraction.
To start with we take the initial point z 0 in a rectangular region R ∈ C that contains all the zeros of a polynomial P(z). The iterative method, when starting from point z 0 in a rectangle, either converges to the zero P(z) or eventually diverges. Stopping condition for convergence is considered as 10 −3 to a maximum of 25 iterations. If the required tolerance is not achieved in 25 iterations, we conclude that the iterative scheme starting at point z 0 does not converge to any root. The strategy adopted is as follows: A color is allocated to each initial point z 0 in the basin of attraction of a zero. If the iteration initiating at z 0 converges, then it represents the attraction basin with that assigned color to it, otherwise in the failing (divergence) situation in 25 iterations the iteration represents the black color.
We analyze the basins of attraction of the new method (for the choices β = 10 −2 , 10 −4 , 10 −8 ) on following three polynomials: Example 1. In the first case, consider the polynomial P 1 (z) = z 2 − 1 which has zeros {±1}. A grid of 400 × 400 points in a rectangle D ∈ C of size [−2, 2] × [−2, 2] is used for drawing the graphics. We assign the color red to each initial point in the basin of attraction of zero '1' and the color green to the points in the basin of attraction of zero '−1'. The graphics are shown in Figure 1 corresponding to β = 10 −2 , 10 −4 , 10 −8 . Observing the behavior of the basins of the new method, we conclude that the convergence domain becoming wider as parameter β assumes smaller values since black zones (divergent points) are getting smaller in size.

Numerical Tests
In this section, some numerical tests on different problems are performed to demonstrate the convergence behavior and computational efficiency of the method M 5,1 . A comparison between the performance of M 5,1 with the existing methods M 2,1 , M 3,1 , M 4,j (j = 1, 2, 3), M 5,2 , M 6,1 , M 7,1 and M 7,2 is also drawn. The programs are performed in the processor with specifications Intel (R) Core (TM) i5-4210U CPU @ 1.70 GHz 2.40 GHz (64-bit Operating System) Microsoft Windows 10 Professional and are complied by Mathematica 10.0 using multiple-precision arithmetic. We record the number of iterations (k) required to converge to the solution such that the stopping condition is satisfied. In order to verify the theoretical order of convergence, the computational order of convergence (p c ) is obtained by using the Formula (8). In the comparison of performance of considered methods, we also include the real CPU time elapsed during the execution of program computed by the Mathematica command "TimeUsed[ ]".
The methods M 2,1 , M 3,1 , M 4,3 , M 5,1 and M 7,2 are tested by using the value 0.01 for the parameter β. In numerical experiments we consider the following five problems: Example 4. Let us consider the system of two equations (selected from [29]): The initial guess assumed is x (0) = {−1, −2} T for obtaining the solution Example 5. Now considering the mixed Hammerstein integral equation (see [4]): wherein x ∈ C[0, 1]; s, t ∈ [0, 1] and the kernel G is The above equation is transformed to a finite-dimensional problem by using the Gauss-Legendre quadrature formula where the weights j and abscissas t j are obtained for m = 8 by Gauss-Legendre quadrature formula. Then, setting x(t i ) = x i , i = 1, 2, ....., 8, we obtain the following system of nonlinear equations wherein the abscissas t j and the weights j are known and produced in Table 1 for m = 8. The initial approximation assumed is
Methods ||x (2) − x (1) || ||x (3) − x (2) || ||x (4) − x (3) || k p c C p,i E p,i e-Time From the numerical results displayed in Tables 3-7, it can be observed that like that of the existing methods the proposed new method shows consistent convergence behavior. Seventh order methods produce approximations with large accuracy due to their higher order of convergence, but they are less efficient. In Example 6, M 4,1 and M 7,1 do not converge to the required solution α 1 . Instead, they converge to solution α 2 which is far off from initial approximation chosen. Calculation of computational order of convergence shows that the order of convergence of the new method is preserved in all the numerical examples. However, this is not true for some existing methods, e.g., M 4,j (j = 1, 2, 3), M 5,2 , M 6,1 , M 7,1 and M 7,2 , in Example 8. Values of the efficiency index shown in the penultimate column of each table also verify the theoretical results stated in Theorem 2. The efficiency results are also in complete agreement with the CPU time utilized in the execution of the program since the method with large efficiency uses less computing time than the method with small efficiency. Moreover, the proposed method utilizes less CPU time than existing higher order methods which points to the dominance of the method. In fact, the new method is especially more efficient for large systems of nonlinear equations.

Conclusions
In the foregoing study, we have developed a fifth order iterative method for approximating solution of systems of nonlinear equations. The methodology is based on third order Traub-Steffensen method and further developed by using derivative free modification of classical Chebyshev's method. The iterative scheme is totally derivative-free and so particularly suitable to those problems where derivatives are lengthy to compute. To prove the local fifth order of convergence for the new method, a development of first-order divided difference operator and direct computation by Taylor's expansion are used.
We have examined the computational efficiency of the new method. A comparison of efficiencies with that of the existing most efficient methods is also performed. It is proved that, in general, the new algorithm is more efficient. Numerical experiments are performed and the performance is compared with existing derivative-free methods. From numerical results it has been observed that the proposed method has equal or better convergence compared to existing methods. Theoretical results related to convergence order and computational efficiency have also been verified in the considered numerical problems. Similar numerical tests, performed for a variety of other different problems, have confirmed the above drawn conclusions to a good extent.

Conflicts of Interest:
The authors declare no conflict of interest.