An Optimal Fourth Order Derivative-Free Numerical Algorithm for Multiple Roots

: A plethora of higher order iterative methods, involving derivatives in algorithms, are available in the literature for ﬁnding multiple roots. Contrary to this fact, the higher order methods without derivatives in the iteration are difﬁcult to construct, and hence, such methods are almost non-existent. This motivated us to explore a derivative-free iterative scheme with optimal fourth order convergence. The applicability of the new scheme is shown by testing on different functions, which illustrates the excellent convergence. Moreover, the comparison of the performance shows that the new technique is a good competitor to existing optimal fourth order Newton-like techniques.

Numerous higher order methods, either independent or based on the modified Newton's method [4]: have been studied in the literature (see, for example, [5][6][7][8][9][10][11][12][13][14][15][16][17]). Such techniques require the information of either first order derivatives or both first and second order derivatives. On the contrary, higher order derivative-free methods to handle the case of multiple zeros are seldom explored in the literature. These methods are important in cases where the derivative Φ of Φ is expensive to compute. One such method is the Traub-Steffensen iteration [18], which uses: for the derivative Φ in Newton Formula (1).
is the first order divided difference. Then, the Newton method (1) takes the form of the modified Traub-Steffensen scheme: Very recently, Kumar et al. [19] and Sharma et al. in [20,21] proposed one-point, two-point, and three-point derivative-free methods with second, fourth, and eighth order convergence, respectively, to compute the multiple zeros. The number of function evaluations corresponding to second, fourth, and eighth order methods is two, three, and four, and so, as per the Kung-Traub hypothesis, these methods possess optimal convergence [22]. The main goal of this work is to construct derivative-free multiple root numerical methods of good computational efficiency, which means the iterative methods of high convergence order with the minimum number of function evaluations. This leads us to develop a two-step derivative-free scheme with fourth order convergence. The proposed scheme consumes only three function evaluations per full iteration, and hence, it is optimal in the sense of the Kung-Traub hypothesis [22]. The procedure is based on the classical Traub-Steffensen iteration (2) in the first step and the Traub-Steffensen-type iteration in the second step.

Development of the Scheme
In what follows, we develop an iterative scheme to compute a multiple root of equation Φ(u) = 0 with multiplicity m > 1. Let us consider the following two-step scheme based on (2): where α 1 , α 2 , α 3 , α 4 are unknown parameters and s k = m Φ(w k ) Φ(u k ) . Observe that this scheme uses the first step as the Traub-Steffensen iteration (2) and the next step as the Traub-Steffensen-like iteration.
We consider principal analytic branches of s k since it is a one-to-m multi-valued function. Hence, it is convenient to treat it as the principal root. For example, the principal root is given by this convention of Arg(p) for p ∈ C agrees with that of the Log[p] command of Mathematica [23] to be employed later in the sections of basins of attraction and numerical experiments.
For the sake of clarity, we prove the results separately for different cases depending on the multiplicity m. Firstly, we consider the case m = 2 and prove the following result: Theorem 1. Let the mapping f : C → C be analytic in a domain containing a multiple zero (say, α) having multiplicity m = 2. Suppose that the starter u 0 is sufficiently close to α, then the formula (3) has the convergence order four, provided that α 1 = 1 4α 3 , α 2 = − 1 2α 3 , and α 4 = 2α 3 .
We can obtain at least fourth order convergence by setting the coefficients of σ 2 k and σ 3 k simultaneously equal to zero. The resulting equations imply that: Consequently, the above error equation reduces to: This establishes the fourth order convergence. Proof. Taking into account that Φ(α) = 0, Φ (α) = 0, Φ (α) = 0 and Φ (3) (α) = 0, the Taylor series development of Φ(u k ) about α gives: where P n = 3!
for n ∈ N.
We state the following theorems for the cases m = 4, 5, 6 (without proof). . Moreover, the scheme satisfies the error equation: for n ∈ N. . Moreover, the scheme satisfies the error equation: for n ∈ N. . Moreover, the scheme satisfies the error equation: for n ∈ N.

Remark 1.
Observing the error equation of each of the above cases, we see that the parameter β does not appear in the equations for m = 4, 5, 6. This leads to the notion: when m ≥ 4, the error equation in each such case would not contain the β term. We shall prove this fact in the next section.

Generalization of the Method
Based on the previous ideas for the case m ≥ 4, we prove the fourth order convergence of Formula (3) by the following theorem: Theorem 6. Using the conditions of Theorem 1, the order of convergence of Formula (3) for the case m ≥ 4 is at least four, if , α 4 = −mα 3 . Moreover, the error equation in the scheme is given by: for n ∈ N.
Similarly, the expansion of Φ(v k ) about α is: where From the first step of Equation (3): . (20) The expansion of Φ(w k ) around α yields: Using (18) and (21) in the expressions of s k , we have that: . (22) Inserting (18)-(22) in the second step of (3), then we have: where Make the coefficients of σ 2 k and σ 3 k simultaneously equal to zero to obtain fourth order convergence. The resulting equations yield: As a result, the error equation is given by: Thus, the theorem is proven.

Remark 2.
The new scheme (3) attains the convergence rate of order four using the conditions of Theorems 1, 2, and 6. This rate is achieved by utilizing only three functional evaluations viz. Φ(u k ), Φ(v k ), and Φ(w k ) per iteration. Therefore, the scheme (3) is optimal by the Kung-Traub hypothesis [22].

Remark 3.
The parameter β used in the expression of v k appears only in the error expressions of the cases m = 2, 3 and not for m ≥ 4 (see Equation (25)). However, for the case m ≥ 4, we have seen that this parameter is included in the coefficient of σ 5 k and in the higher order terms. These terms are lengthy to evaluate in general. Moreover, we do not need these in order to show the desired fourth order convergence.

Remark 4.
For future reference, the proposed method (3) is written as: wherein we have considered α 4 = 2α 3 . This iterative scheme satisfies the common conditions of Theorems 1, 2, and 6.

Basins of Attraction
Our point here is to analyze the new technique by a geometrical tool, namely the basins of attraction of the zeros of a polynomial function P(z) in the Argand plane. The examination of the basins of attraction of roots by the iterative scheme gives a significant idea about the convergence of the scheme. This thought was presented in [24][25][26][27][28][29][30][31]. In order to assess the basins, we take 10 −3 as the stopping condition for convergence, but to a maximum of 25 iterations. If this tolerance is not achieved in the required iterations, the procedure is dismissed with the result showing the divergence of the iteration function starting from z 0 . While drawing the basins, the following criterion is adopted: a color is allotted to every initial guess z 0 in the attraction basin of a zero. If the iterative formula beginning at the point z 0 converges, then it forms the basins of attraction with that assigned color, and if the formula fails to converge in the required number of iterations, then it is painted with black color.
We discuss the basins of attraction by applying the method (26) (choosing β = 10 −2 , 10 −4 , 10 −6 ) on the following two polynomials: Problem 1. In this example, we consider the polynomial P 1 (z) = (z 2 + 5z + 6) 2 , which has zeros {−3, −2} with multiplicity two. For this situation, we utilize a square shape D ∈ C of size [−4, 4] × [−4, 4] and allot the shading of yellow and blue to each underlying point in the attraction basin of zero "−2" and "−3". The basins obtained for the method appear in Figure 1 for the parameter values β = 10 −2 , 10 −4 , 10 −6 . Observing the behavior of Method (3), we see that the basins become more qualitative as parameter β accepts small values. It is clear that the attraction basins show the convergence behavior and suitability of an iterative scheme relying on the conditions. In the event that we pick a starter z 0 in a zone where different basins of attraction meet one another, it is difficult to foresee which root will be attained by the method that starts in z 0 . Thus, the selection of z 0 in such a zone is not a wise decision. The graphics show that the dark zones and the zones with various hues are not appropriate to speculate about z 0 when we need to accomplish a specific root. The most alluring pictures show up when we have extremely perplexing wildernesses between basins of attraction. We close this segment with a comment that the nature of the proposed technique relies on the estimation of parameter β. The smaller the estimation of β, the better the convergence of the method.

Numerical Results
In order to check the performance and validity of the theoretical results, we apply the new method (NM) to solve some nonlinear problems. The theoretical fourth order of convergence is verified by using the formula of the approximate computational order of convergence (ACOC; see [32]): The performance of NM is compared with some existing well-known optimal fourth order methods with and without derivative evaluations.

Sharma-Sharma method (SSM):
Zhou-Chen-Song method (ZM): Soleymani-Babajee-Lotfi method (SM): Kansal-Kanwar-Bhatia method (KM): where p = m m+2 . Sharma-Kumar-Jäntschi method (SM-1): Sharma-Kumar-Jäntschi method (SM-2): . The calculations were executed in the programmable package of the Mathematica software [23] with higher precision. The results of the new method were obtained by choosing a value of 0.01 for the parameter β. The numerical results shown in Tables 1-4 include: (i) the number of iterations (k) that are needed to converge to the required solution such that |u k+1 − u k | + |Φ(u k )| < 10 −100 , (ii) the estimated error |u k+1 − u k | in the first three iterations, (iii) the approximated computational order of convergence (ACOC) using Formula (27), and (iv) the elapsed CPU time to run the program measured by the Mathematica command "TimeUsed[ ]".
The following numerical examples were selected for testing: Example 1. Consider the van der Waals equation: which shows the nature of a real gas by the inclusion of parameters a 1 and a 2 in the ideal gas equation. The volume V in terms of the remaining parameters can be found by the equation: For a given a set of values of a 1 and a 2 of a particular gas, one can find values of n, P, and T, so that the equation has three roots. For a particular set of values, we have: This function has three zeros: one is a simple zero α = 1.72 and the other is a repeated zero α = 1.75 of multiplicity two. The methods were tested for initial guess u 0 = 2.3 to find the desired zero α = 1.75.
The computed results are given in Table 1.

Example 2.
Let λ, c, T, k, and h be wavelength of the radiation, the speed of light, the absolute temperature of the black body, Boltzmann's constant, and Planck's constant, respectively. Then, the Planck law of radiation [33] to find the energy density in a black body is given as: One wants to determine the wavelength λ corresponding to the maximum energy density φ(λ). From Equation (28), it follows that:
Setting u = ch/λkT, the above equation assumes the form: The following function can be obtained by considering the above case three times: The trivial zero u = 0 is not taken for discussion. The non-trivial zero can be guessed somewhere near u = 5, since for u = 5, the left-hand side of (29) is zero and the right-hand side is e −5 ≈ 6.74 × 10 −3 . In fact the expected root α ≈ 4.96511423174427630369 is obtained with the starter u 0 = 5.4. Thereby, the wavelength (λ) for the maximum energy density is: .
The numerical results are displayed in Table 2.
Example 3. The problem of isentropic supersonic flow around a sharp expansion corner is considered. The mathematical expression between the Mach number before the corner (i.e., M 1 ) and after the corner (i.e., M 2 ) is defined by (see [34]): where b = γ+1 γ−1 , γ being the specific heat ratio of the gas. Given that M 1 = 1.5, γ = 1.4 and δ = 10 0 , we solve the equation for M 2 . This case results in: Let us consider this equation four times, and so, the required function is: This function possesses zero α ≈ 1.8411027704926161 with multiplicity four. The required zero is calculated using starter u 0 = 1.5. The results obtained by the methods are displayed in Table 3.

Example 4.
Lastly, the standard test function defined by: is considered. The function Φ 4 has multiple imaginary zero α = i of multiplicity five. We select the initial value u 0 = 1.3i to obtain the required zero of the function. The numerical results so produced are shown in Table 4. It can be seen from the numerical results displayed in Tables 1-4 that the proposed method supported the theoretical results proven in Sections 2 and 3 and, like existing methods, possessed consistent convergence. Moreover, the CPU time consumed by the methods as shown in the tables proved computationally efficient nature of the new technique as compared to the CPU time of the considered existing methods of the same order. The main aim of applying the new derivative-free method on different nonlinear equations was purely to show its accuracy and consistency. Similar numerical testing, carried out for a number of problems of different types, confirmed the above conclusions to a good extent.

Conclusions
In the study, we proposed a derivative-free numerical method with optimal fourth order convergence for approximating the repeated roots of nonlinear equations. The analysis of the convergence under standard hypotheses proved the convergence order four. The method was employed to solve nonlinear problems including those arising from real-world applications. The performance was compared with existing techniques (with and without derivatives) of the same order. The testing of the numerical results showed the presented derivative-free method as a good competitor of the already established fourth order techniques that use derivative information in the algorithm. We conclude this work with a remark: the derivative-free method presented here can be a better choice compared to existing Newton-type schemes in the cases where derivatives are difficult to obtain or expensive to compute.