An Efficient Class of Traub-Steffensen-Like Seventh Order Multiple-Root Solvers with Applications

Abstract: Many higher order multiple-root solvers that require derivative evaluations are available in literature. Contrary to this, higher order multiple-root solvers without derivatives are difficult to obtain, and therefore, such techniques are yet to be achieved. Motivated by this fact, we focus on developing a new family of higher order derivative-free solvers for computing multiple zeros by using a simple approach. The stability of the techniques is checked through complex geometry shown by drawing basins of attraction. Applicability is demonstrated on practical problems, which illustrates the efficient convergence behavior. Moreover, the comparison of numerical results shows that the proposed derivative-free techniques are good competitors of the existing techniques that require derivative evaluations in the iteration.


Introduction
Solving nonlinear equations is an important task in numerical analysis and has numerous applications in engineering, mathematical biology, physics, chemistry, medicine, economics, and other disciplines of applied sciences [1][2][3].Due to advances in computer hardware and software, the problem of solving the nonlinear equations by computational techniques has acquired an additional advantage of handling the lengthy and cumbersome calculations.In the present paper, we consider iterative techniques for computing multiple roots, say α, with multiplicity m of a nonlinear equation f (x) = 0, that is f (j) (α) = 0, j = 0, 1, 2, . . ., m − 1 and f (m) (α) = 0.The solution α can be calculated as a fixed point of some function M : D ⊂ C → C by means of the fixed point iteration: where x ∈ D is a scalar.Many higher order techniques, based on the quadratically-convergent modified Newton's scheme (see [4]): have been proposed in the literature; see, for example, [5][6][7][8][9][10][11][12][13][14][15][16][17][18][19][20] and the references therein.The techniques based on Newton's or the Newton-like method require the evaluations of derivatives of first order.There is another class of multiple-root techniques involving derivatives of both the first and second order; see [5,21].However, higher order derivative-free techniques to handle the case of multiple roots are yet to be explored.The main problem of developing such techniques is the difficulty in finding their convergence order.Derivative-free techniques are important in the situations when the derivative of the function f is difficult to compute or is expensive to obtain.One such derivative-free technique is the classical Traub-Steffensen method [22].The Traub-Steffensen method actually replaces the derivative in the classical Newton's method with a suitable approximation based on the difference quotient, or writing more concisely: where is a first order divided difference.In this way, the modified Newton's scheme (2) assumes the form of the modified Traub-Steffensen scheme: The Traub-Steffensen scheme (3) is a noticeable improvement of Newton's scheme, since it maintains the quadratic convergence without using any derivative.
The aim of the present contribution is to develop derivative-free multiple-root iterative techniques with high computational efficiency, which means the techniques that may attain a high convergence order using as small a number of function evaluations as possible.Consequently, we develop a family of derivative-free iterative methods of seventh order convergence that requires only four function evaluations per full iteration.The scheme is composed of three steps, out of which the first step is the classical Traub-Steffensen iteration (3) and the last two steps are Traub-Steffensen-like iterations.The methodology is based on the simple approach of using weight functions in the scheme.Many special cases of the family can be generated depending on the different forms of weight functions.The efficacy of the proposed methods is tested on various numerical problems of different natures.In the comparison with existing techniques requiring derivative evaluations, the new derivative-free methods are observed to be computationally more efficient.
We summarize the contents of the rest of paper.In Section 2, the scheme of the seventh order multiple-root solvers is developed and its order of convergence is determined.In Section 3, the basins of attractors are presented to check the stability of new methods.To demonstrate the performance and comparison with existing techniques, the new techniques are applied to solve some practical problems in Section 4. Concluding remarks are given in Section 5.

Development of the Family of Methods
Given a known multiplicity m ≥ 1, we consider a three-step iterative scheme with the first step as the Traub-Steffensen iteration (3) as follows: in a neighborhood of 0, and the function G(u, w) : C × C → C is holomorphic in a neighborhood of (0, 0).Note that second and third steps are weighted by the factors H(u) and G(u, w), so these factors are called weight factors or weight functions.We shall find conditions under which the scheme (4) achieves convergence order as high as possible.In order to do this, let us prove the following theorem: Theorem 1.Let f : C → C be an analytic function in a region enclosing a multiple zero α with multiplicity m.Assume that initial guess x 0 is sufficiently close to α, then the iteration scheme defined by (4) possesses the seventh order of convergence, provided that the following conditions are satisfied: Proof.Let the error at the n th iteration be e n = x n − α.Using the Taylor expansion of f (x n ) about α, we have that: where Using (5) in t n = x n + β f (x n ), we obtain that: Taylor's expansion of f (t n ) about α is given as: By using Equations ( 5)- (7) in the first step of (4), after some simple calculations, it follows that: where The expressions of ω i are very lengthy, so they are not written explicitly.
Taylor's expansion of f (y n ) about α is given by: where ωi = ωi (m, C 1 , C 2 , . . ., C 7 ).By using ( 5) and ( 9), we get the expression of u as: where . We expand the weight function H(u) in the neighborhood of 0 by the Taylor series, then we have that: Inserting Equations ( 5), (9), and (11) in the second step of Scheme (4) and simplifying, where In order to attain higher order convergence, the coefficients of e 2 n and e 3 n should be simultaneously equal to zero.That is possible only for the following values of H(0) and H (0): By using the above values in (12), we obtain that: Expansion of f (z n ) about α leads us to the expression: From ( 5), (9), and (15), we get the expressions of v and w as: and: where τ i and ς i are some expressions of m, H (0), H (0), C 1 , C 2 , . . ., C 7 .
Expanding the function G(u, w) in the neighborhood of origin (0, 0) by Taylor series: where G ij = ∂ i+j ∂u i ∂w j G(u, w)| (0,0) .Then by substituting (5), ( 16), (17), and (18) into the last step of Scheme (4), we obtain the error equation: where It is clear from Equation ( 19) that we will obtain at least fifth order convergence if we have G 00 (0, 0) = 1.Moreover, we can use this value in ξ 1 = 0 to obtain: By using G 00 = 1 and (20) in ξ 2 = 0, the following equation is obtained: Using the above values in (19), the final error equation is given by: Hence, the seventh order convergence is established.

Forms of the Weight Function
Numerous special cases of the family (4) are generated based on the forms of weight functions H(u) and G(u, w) that satisfy the conditions of Theorem 1.However, we restrict ourselves to simple forms, which are given as follows:

I. Some particular forms of H(u)
Case I(a).When H(u) is a polynomial weight function, e.g., By using the conditions of Theorem 1, we get A 0 = 1, A 1 = 2 and A 2 = −1.Then, H(u) becomes:

Case I(b).
When H(u) is a rational weight function, e.g., Using the conditions of Theorem 1, we get that . Therefore,

Case I(c).
When H(u) is a rational weight function, e.g., Using the conditions of Theorem 1, then we get A 0 = 3, A 1 = 2, and A 2 = 1.H(u) becomes:

Case I(d).
When H(u) is a rational function of the form: Using the conditions of Theorem 1, we get

II. Some particular forms of G(u, w)
Case II(a).When G(u, w) is a polynomial weight function, e.g., Using the conditions of Theorem 1, then we get where A 4 and A 5 are free parameters.

Case II(b).
When G(u, w) is a rational weight function, e.g., Using the conditions of Theorem 1, we have where A 3 and B 3 are free parameters.

Case II(c).
When G(u, w) is the sum of two weight functions H 1 (u) and By using the conditions of Theorem 1, we get: where B 2 is a free parameter.

Case II(d).
When G(u, w) is the sum of two rational weight functions, that is: By using the conditions of Theorem 1, we obtain that: Case II(e).When G(u, w) is product of two weight functions, that is: Using the conditions of Theorem 1, then we get:

Complex Dynamics of Methods
Here, our aim is to analyze the complex dynamics of the proposed methods based on the visual display of the basins of attraction of the zeros of a polynomial p(z) in the complex plane.Analysis of the complex dynamical behavior gives important information about the convergence and stability of an iterative scheme.Initially, Vrscay and Gilbert [23] introduced the idea of analyzing complex dynamics.Later on, many researchers used this concept in their work, for example (see [24][25][26] and the references therein).We choose some of the special cases (corresponding to the above forms of H(u) and G(u, w)) of family ( 4 We take the initial point as z 0 ∈ D, where D is a rectangular region in C containing all the roots of p(z) = 0.The iterative methods starting at a point z 0 in a rectangle either converge to the zero of the function p(z) or eventually diverge.The stopping criterion considered for convergence is 10 −3 up to a maximum of 25 iterations.If the desired tolerance is not achieved in 25 iterations, we do not continue and declare that the iterative method starting at point z 0 does not converge to any root.The strategy adopted is the following: A color is assigned to each starting point z 0 in the basin of attraction of a zero.If the iteration starting from the initial point z 0 converges, then it represents the basins of attraction with that particular color assigned to it, and if it fails to converge in 25 iterations, then it shows the black color.
To view complex geometry, we analyze the attraction basins for the methods NM-1(a-d) and NM-2(a-d) on the following two polynomials: Problem 1.In the first example, we consider the polynomial p 1 (z) = (z 2 − 1) 2 , which has zeros {±1} with multiplicity two.In this case, we use a grid of 400 × 400 points in a rectangle and assign the color green to each initial point in the basin of attraction of zero ' − 1 and the color red to each point in the basin of attraction of zero '1 .Basins obtained for the methods NM-1(a-d) and NM-2(a-d) are shown in Figures 1-4 corresponding to β = 0.01, 0.002.Observing the behavior of the methods, we see that the method NM-2(d) possesses a lesser number of divergent points and therefore has better stability than the remaining methods.Notice that there is a small difference in the basins for the rest of the methods with the same value of β.Notice also that the basins are becoming wider as parameter β assumes smaller values.From the graphics, we can easily observe the behavior and applicability of any method.If we choose an initial guess z 0 in a region wherein different basins of attraction touch each other, it is difficult to predict which root is going to be attained by the iterative method that starts from z 0 .Therefore, the choice of z 0 in such a region is not a good one.Both black regions and the regions with different colors are not suitable to assume the initial guess as z 0 when we are required to achieve a particular root.The most intricate geometry is between the basins of attraction, and this corresponds to the cases where the method is more demanding with respect to the initial point.We conclude this section with a remark that the convergence behavior of the proposed techniques depends on the value of parameter β.The smaller the value of β is, the better the convergence of the method.

Numerical Examples and Discussion
In this section, we implement the special cases NM-1(a-d) and NM-2(a-d) that we have considered in the previous section of the family (4), to obtain zeros of nonlinear functions.This will not only illustrate the methods' practically, but also serve to test the validity of theoretical results.The theoretical order of convergence is also confirmed by calculating the computational order of convergence (COC) using the formula (see [27]): The performance is compared with some well-known higher order multiple-root solvers such as the sixth order methods by Geum et al. [8,9], which are expressed below: First method by Geum et al. [8]: where function in the neighborhood of origin (0, 0).The authors have also studied various forms of the function Q f leading to sixth order convergence of (21).We consider the following four special cases of function Q f (u, s) in the formula (21) and denote the corresponding methods by GKN-1(j), j = a, b, c, d: , where , Second method by Geum et al. [9]: where of 0, and K f : C 2 → C is holomorphic in a neighborhood of (0, 0).Numerous cases of G f and K f have been proposed in [9].We consider the following four special cases and denote the corresponding methods by GKN-2(j), j = a, b, c, d: Computations were carried out in the programming package Mathematica with multiple-precision arithmetic.Numerical results shown in Tables 1-4 include: (i) the number of iterations (n) required to converge to the solution, (ii) the values of the last three consecutive errors e n = |x n+1 − x n |, (iii) the computational order of convergence (COC), and (iv) the elapsed CPU time (CPU-time).The necessary iteration number (n) and elapsed CPU time are calculated by considering The convergence behavior of the family of iterative methods ( 4) is tested on the following problems: Example 1. (Eigenvalue problem).Finding eigenvalues of a large square matrix is one of the difficult tasks in applied mathematics and engineering.Finding even the roots of the characteristic equation of a square matrix of order greater than four is a big challenge.Here, we consider the following 6× 6 matrix: The characteristic equation of the above matrix (M) is given as follows: This function has one multiple zero at α = 1 of multiplicity three.We choose initial approximation x 0 = 0.25.Numerical results are shown in Table 1.
A numerical study, for different values of the parameters α and K, has been performed in [28].As a particular example, let us take α = 1 4 and K = π 5 .Consider this particular case four times with same values of the parameters, then the required nonlinear function is: This function has one multiple zero at α = 0.80926328 . . . of multiplicity four.The required zero is calculated using initial approximation x 0 = 1.Numerical results are displayed in Table 2.The relationship between the Mach number before the corner (i.e., M 1 ) and after the corner (i.e., M 2 ) is given by (see [3]): where b = γ + 1 γ − 1 and γ is the specific heat ratio of the gas.
As a particular case study, the equation is solved for M 2 given that M 1 = 1.5, γ = 1.4,and δ = 10 0 .Then, we have that: Consider this particular case three times with the same values of the parameters, then the required nonlinear function is: This function has one multiple zero at α = 1.8411027704 . . . of multiplicity three.The required zero is calculated using initial approximation x 0 = 1.5.Numerical results are shown in Table 3.
which has a multiple zero at α = −0.72855964390156 . . . of multiplicity four.Numerical results are shown in Table 4 with initial guess x 0 = −0.5.Example 5. Consider the standard function, which is given as (see [8]): The multiple zero of function f 5 is α = 2 with multiplicity five.We choose the initial approximation x 0 = 1.5 for obtaining the zero of the function.Numerical results are exhibited in Table 5.
which has a zero α = 3 of multiplicity three.Let us choose the initial approximation x 0 = 3.5 for obtaining the zero of the function.Numerical results are shown in Table 6.
The zero of function f 7 is α = 1.4632625480850 . . .with multiplicity four.We choose the initial approximation x 0 = 1.3 to find the zero of this function.Numerical results are displayed in Table 7.It is clear from the numerical results shown in Tables 1-7 that the accuracy in the successive approximations increases as the iterations proceed.This shows the stable nature of the methods.Moreover, the present methods like that of existing methods show consistent convergence behavior.We display the value zero of |e n | in the iteration at which |x n+1 − x n | + |F(x n )| < 10 −350 .The values of the computational order of convergence exhibited in the penultimate column in each table verify the theoretical order of convergence.However, this is not true for the existing methods GKN-1(a-d) and GKN-2(a) in Example 2. The entries in the last column in each table show that the new methods use less computing time than the time used by existing methods.This verifies the computationally-efficient nature of the new methods.Similar numerical tests, performed for many problems of different types, have confirmed the aforementioned conclusions to a large extent.
We conclude the analysis with an important problem regarding the choice of initial approximation x 0 in the practical application of iterative methods.The required convergence speed of iterative methods can be achieved in practice if the selected initial approximation is sufficiently close to the root.Therefore, when applying the methods for solving nonlinear equations, special care must be given for guessing close initial approximations.Recently, an efficient procedure for obtaining sufficiently close initial approximation has been proposed in [29].For example, the procedure when applied to the function of Example 1 in the interval [0, 1.5] using the statements: in programming package Mathematica yields a close initial approximation x 0 = 1.04957 to the root α = 1.

Conclusions
In the present work, we have designed a class of seventh order derivative-free iterative techniques for computing multiple zeros of nonlinear functions, with known multiplicity.The analysis of convergence shows the seventh order convergence under standard assumptions for the nonlinear function, the zeros of which we have searched.Some special cases of the class were stated.They were applied to solve some nonlinear equations and also compared with existing techniques.Comparison of the numerical results showed that the presented derivative-free methods are good competitors of the existing sixth order techniques that require derivative evaluations.The paper is concluded with the remark that unlike the methods with derivatives, the methods without derivatives are rare in the literature.Moreover, such algorithms are good options to Newton-like iterations in the situation when derivatives are difficult to compute or expensive to obtain.
) to analyze the basins.Let us choose the combinations of special Cases II(c) (for B 2 = 1) and II(d) with I(a), I(b), I(c), and I(d) in (4) and denote the corresponding new methods by NM-i(j), i = 1, 2 and j = a, b, c, d.

Figure 4 .Problem 2 .
Figure 4. Basins of attraction for NM-2 (for β = 0.002) in polynomial p 1 (z).Problem 2. Let us take the polynomial p 2 (z) = (z 3 + z) 3 having zeros {0, ±i} with multiplicity three.To see the dynamical view, we consider a rectangle D = [−2, 2] × [−2, 2] ∈ C with 400 × 400 grid points and allocate the colors red, green, and blue to each point in the basin of attraction of −i, 0, and i, respectively.Basins for this problem are shown in Figures 5-8 corresponding to parameter choices β = 0.01, 0.002 in the proposed methods.Observing the behavior, we see that again, the method NM-2(d) has better convergence behavior due to a lesser number of divergent points.Furthermore, observe that in each case, the basins are becoming larger with the smaller values of β.The basins in the remaining methods other than NM-2(d) are almost the same.

Table 1 .
Comparison of the performance of methods for Example 1.

Table 6 .
Comparison of the performance of methods for Example 6.