An optimal eighth order derivative free multiple root finding numerical method and applications to chemistry

In this paper, we present an optimal eighth order derivative-free family of methods for multiple roots which is based on the first order divided difference and weight functions. This iterative method is a three step method with the first step as Traub–Steffensen iteration and the next two taken as Traub–Steffensen-like iteration with four functional evaluations per iteration. We compare our proposed method with the recent derivative-free methods using some chemical engineering problems modelled as nonlinear equations with simple and multiple roots. Stability of the presented family of methods is demonstrated by using the graphical tool known as basins of attraction.


Introduction
Many iterative methods have been developed for the solution of problems modelled as nonlinear equations with simple or multiple roots but it is very difficult to construct an optimal iterative method for finding multiple roots. It is noteworthy that most of the higher order optimal methods that were obtained during the past years require derivative evaluation of the involved function, as in [1], while higher order derivativefree methods exist rarely in literature because it is a difficult task to preserve the order of convergence when we replace derivatives by differences. Kansal et al. [7] and Kumar et al. [8] gave second order iterative schemes to find repeated roots of nonlinear equations. Sharma et al. [17], Kumar et al. [9,10], Behl et al. [2] and Rani and Kansal [13] proposed a fourth-order root finding methods for multiple roots. Qudsi et al. [11,12] presented three step sixth order iterative methods for finding the multiple roots. Sharma et al. [15] proposed seventh order convergent iterative scheme for multiple roots. Sharma et al. in [14,16] presented eight order scheme for computing multiple root of nonlinear equations. nevertheless, other derivative-free methods for multiple roots have been generated by using different approaches, such as [4].
Motivated by the need to find the more efficient iterative method for computing the root of nonlinear equations, we propose an optimal eighth order derivative-free method for computing multiple root of nonlinear equations with multiplicity η ≥ 1 in the next section. This method is derivative-free with four functional evaluations with a univariate and a multivariate weight functions. In Sect. 3, we discuss some particular cases of weight functions and choose four of them for further investigation. We compare our methods with two of the recent derivative-free schemes of seventh [15] and eighth order [14,16] using physical applications from chemical engineering and dynamical behaviour in Sect. 4. Our proposed method has wider convergence region as compared to recent root-finding iterative method.

Proposed optimal eighth order scheme
Let τ = w be a multiple root with multiplicity η ≥ 1 of the function g. Let K : C → C and L : C 3 → C be analytic functions in the neighborhood of 0 and (0, 0, 0) respectively. We propose a family of derivative-free methods of eighth-order involving first-order divided differences for computing multiple roots with multiplicity η > 1 given by: where p k = ( g(μ k ) g(τ k ) ) 1 η , q k = ( g(υ k ) g(τ k ) ) 1 η and r k = ( g(υ k ) g(μ k ) ) 1 η . Following result can be used to investigate the convergence order of presented method (1) and the conditions to be imposed on weight functions K and L. Theorem 1 Let g : C → C be analytic in a region enclosing multiple root of g with known multiplicity η. Suppose that the initial guess τ 0 be sufficiently close to the multiple zero w. Then, the class of iterative schemes defined by (1) has eighth-order of convergence, when the following conditions are satisfied: where The error equation for the proposed schemes is: Proof Let w be a multiple root of g such that g ( j) (w) = 0, j = 0, 1, 2, . . . , η − 1 with g (η) (w) = 0. Let e k = τ k − w be error in the kth iteration. Considering Taylor expansion of g(τ k ) about w, we have: which can be expressed as: where and j ∈ N. For ρ k = τ k + γ g (τ k ), we obtain: or Using Taylor expansion of g (ρ k ) about w, By employing (3) and (6) in the first step of (1) and after some algebraic manipulations, we get: wherez i =z i (η, d 1 , . . . , d 8 ) for i = 1, 2, 3, 4. Now, the Taylor expansion of g (μ k ) about w is given by; Using (5) and (8) By expanding the weight function K ( p k ) in the neighborhood of 0, Then, for we get the expansion of υ k around w using (5), (8), (9) and (10): If we choose the values of K (0), K (0) and K (0) in (12), given as: we obtain Next, from Taylor expansion of g(υ k ) about w: By using (5) and (14), becomes: Similarly, from (8) and (14) The expansion of the weight function L ( p k , q k , r k ) in the neighborhood of (0, 0, 0) is given by, L ( p k , q k , r k ) = L 000 + p k L 100 + q k L 010 + r k L 001 + p k q k L 110 + p k r k L 101 +q k r k L 011 + r 2 k 2 L 002 , and therefore, e k+1 = υ k − p k L 000 + p k L 100 + q k L 010 + r k L 001 + p k q k L 110 + p k r k L 101 After substituting the values of p k , q k and r k from (9), (15) and (16), the error equation becomes: As L i jk = ∂ i+ j+k ∂ p i ∂q j ∂r k L ( p, q, r ) | (0,0,0) , in order to achieve higher order of convergence, we use the values: By using (18) in (17), Furthermore, from (19), and from conditions (20), Now, by using L 011 = 4 in (21), we get Clearly, for K (0) = 36 in (21), the error equation indicates eighth order of convergence as an optimal order of convergence: From Theorem 1, we can get several new root finding methods for multiple roots by using different cases for K ( p k ) and L( p k , q k , r k ) in the proposed scheme (1). We discuss few particular cases of the proposed class. It is important to note that the choice of specific values of parameter γ can be made under the point of view of an improvement of the stability and a widening of the set of converging initial estimations. However, as in the iterative expression of the new class, g[τ k , ρ k ] is an estimation of g (τ k ), it is clear that the approximation of the derivative is better when γ is close to zero.

Some special cases of weight functions
Many special cases of the proposed scheme (1) can be generated by using different forms of weight functions K ( p) and L( p, q, r ). These weight functions satisfy the conditions given in Theorem 1. We discuss some simple cases as follows: By using the condition set (2), the weight function K ( p) becomes: By taking the other weight function L( p, q, r ) as a second-degree polynomial, and applying the conditions stated in Theorem 1, L( p, q, r ) becomes L( p, q, r ) = 2q + 4qr + r + r 2 .
Therefore, the resulting scheme is denoted by FZ1 and has the iterative expression: Case 2 By considering the rational form of weight function K ( p) as: and imposing conditions (2), we get K ( p) in the form, Let us take L( p, q, r ) as the second-degree polynomial which becomes, Thus, the resulting member of class (1) is denoted by FZ2 and has the following expression: Case 3 Let us consider the weight function K ( p) in the form of a proper rational function: By using the conditions of Theorem 1, K ( p) becomes; Taking again the second weight function L( p, q, r ) as the second-degree polynomial L( p, q, r ) = q + pr + 4qr + r + r 2 , the obtained element of family (1) is denoted by FZ3 and can be expressed as Case 4 Let us consider the weight function K ( p) in the form: After applying conditions (2), the weight function K ( p) becomes: and L( p, q, r ) is the same as in previous cases, So that, the resulting scheme is denoted by FZ4 and has the iterative expression

Numerical and dynamical analysis
In this section, we check the convergence behavior and performance of Case 1 to Case 4 of our proposed eighth order scheme denoted respectively by FZ1, FZ2, FZ3 and FZ4 by carrying out some nonlinear equations from real life applications of chemical engineering. We compare the methods with the recent derivative free methods of seventh order (see [15], Case I(a), Case II(c)) denoted by SH1, SH2 and eighth order (see [16], M-4 and [14]) denoted as SH3 and SH4. The recent seventh order methods are defined by, in case of SH1 and SH2 is The above mentioned eight order scheme (see [16]) denoted by SH3 is defined as where, The most recent eighth order method (see [14]), denoted by SH4, is defined as: For testing those methods numerically, we have performed all computations in computer algebra software Maple 16 using 300 significant digits of precision. We take the value of γ = 0.001. Tables show per step numerical errors of approximating real root |τ k − τ k−1 | of first three iterations, the absolute residual error of the test function at the third iteration, |g(τ 3 )|, and the computational order of convergence (see [6]).
For dynamical analysis, we use the graphical tool known as basins of attraction. For a polynomial function, basins of attraction is also called a polynomiograph. The dynamical analysis gives us information about the stability and convergence regions of iterative schemes. The choice of initial guess decide the convergence of nonlinear function towards the exact root. Many researchers have used basins of attraction to show the convergence behaviour of their schemes. Consider a function g k (τ ) where τ ∈ C and w k be the root of g k (τ ). We have chosen the parameter β = 0.001, 10 −5 as tolerance and we work out the schemes with maximum 15 number of iterations. A grid of 1000 × 1000 points in the complex plane [a, b] × [c, d] is taken into consideration where the values of a, b, c and d depends upon the zero of the nonlinear function. We select ' hot' color map and assign black color to the divergence region with few exceptions. Different shades of colors are decided for different number of iterations used by the iterative method to converge to the root with required accuracy. We compare our presented schemes namely F Z1−F Z4 with recent methods named as S H1−S H4.

Applications to chemical engineering problems
Chemical Engineering deals with the applications of chemistry for industrial purposes. It effectively analyze the chemical methods to convert materials into more applicable and useful materials. Let us consider some examples from chemical engineering.

Example 1 Continuous Stirred Tank Reactors (CSTR)
Let us consider an isothermal continuous stirred tank reactor. The components E and M are fed to the reactor at the rates of Q and q-Q respectively. The following reaction scheme develops in the reactor (see [3]): Table 1 Computational comparison of iterative schemes for g 1 (τ ) We take multiple root w = −2.85 with multiplicity η = 2. By taking the initial guess τ 0 = −3.13, the numerical results are given in Table 1. . When we observe these basins of attraction, we see that our newly proposed scheme have better convergence as compared to the methods S H1 − S H4. We notice that S H1 − S H4 take minimum 2 and maximum 12 iterations, F Z1 − F Z4 take minimum 2 and maximum 10 iterations to converge to the root.

Equations of state
An equation of state (EOS) relates the molar volume, temperature and the pressure of a substance. Its simplest form is the ideal gas equation of state: where P, V and T represents pressure, volume and temperature of a gas respectively. Here n is the number of moles of gas and R is the universal gas constant. There are many improved variants of EOS but we will not discuss all of them. We start with the virial equation of state which is the first improved form of EOS. We also discuss four improved cubic EOS: the Van Der Waals (VW), the Redlich-Kwong (RK), Soave Redlich-Kwong (SRK) and Peng-Robinson (PR) forms. All these equations can be written in the following form (see [18]): with different values of constants c 1 and c 2 . Here a and b are the parameters defined as:    Table 2.   It can be observed that the computational estimation of the order of convergence COC is the best at FZ1-FZ4, meanwhile the error at the third iteration is the lowest at SH1-SH2 and FZ2.

Example 2 Virial Equation of State Let us consider the virial equation of state:
In Example 2, the basins of attraction obtained for the methods SH1-SH4 and  conclude that SH1-SH4 take minimum two and maximum 10 numbers of iterations to converge to the root. On the other hand, FZ1-FZ4 take minimum two and maximum 8 iterations to converge to the root. In this example, black color shows the divergence region only for SH3 and in all other cases, it indicates the hue of two iterations.

Example 3 Van Der Waals Equation of
State By using the values c 1 = c 2 = 0, λ b = 1 3 , λ a = 9 8 , z = 3 8 and α(T r ) = 1, we get the reduced Eq. (22) in the form: For finding the volume V of a particular gas, we have to use the values of remaining parameters. For the given parameters a and b of a particular gas, we can get values of n (number of moles), P (pressure) and T (absolute temperature), such that the equation (23) has multiple root with multiplicity 3. By using these values, we get the nonlinear function g 3 (τ ) = τ 3 − 5.22τ 2 + 9.0825τ − 5.2675, Table 3 Computational comparison of iterative schemes for g 3 (τ ) Methods  where τ = V . This equation has the multiple root w = 1.75 with multiplicity η = 2. By taking the initial guess τ 0 = 2.00, the numerical results are presented in Table 3.
Again proposed methods reach better estimations of the theoretical order of convergence in the first three iterations. Moreover, the best error is obtained by SH4, FZ3 and FZ4.
Let us solve the above equation for Carbon Dioxide (C O 2 ) at T = T c = 304.2, P = P c = 72.85. By using these values, we get the nonlinear equation: Methods where R = 0.082057366080960, and solving for τ = V , gives us the multiple root w = 0.1142157436 with multiplicity η = 3. By using the initial guess τ 0 = 0.8, the numerical results are shown in Table 4, with excellent results in the first three iterations.
In Example 4, the basins of attraction obtained for the methods SH1-SH4 and FZ1-FZ4 are shown in Figs. 7, 8. We draw basins of attraction in the region [−2, 2]×[−2, 2] and conclude that SH1-SH4 and FZ1-FZ4 use minimum two and maximum three numbers of iterations to converge to the root. In this example, black color shows the hue of two iterations for all iterative schemes. FZ1-FZ4 has globally wider black region as compared to SH1-SH4.

Example 5 Soave Redlich-Kwong Equation of State
w represents the acentric factor of a specific gas. We solve this equation for Ammonia at T = 302.15, P = 229.9×10 3 . We use T c = 405.55, P c = 1.128×10 7 , w = 0.250 and R = 8.314 and get the nonlinear equation: where τ = V . It has a simple root w = 0.0001547767475. We choose the initial guess τ 0 = 0.0003885 and the numerical results are presented in Table 5.  Table 5 shows the best performance of FZ1-FZ4 schemes, even for simple roots.
For g 5 (τ ), the basins of attraction obtained for the methods S H1 − S H4 and F Z1 − F Z4 are shown in Figs. 9, 10. For this example, we consider the region [−1, 1] × [−1, 1] for dynamical analysis. The basins of attraction of S H1 − S H4 and F Z1 − F Z4 indicate that all iterative schemes take minimum two and maximum six iterations to converge to the root except S H3 which take minimum 2 and maximum 7 numbers of iterations for convergence towards the root. For all iterative schemes, hue of black color shows the divergence region except F Z3 and S H3 in which it indicates two iterations.

Example 6 Peng-Robinson Equation of State
We also consider PR equation of state for finding multiple root of V at critical isotherm by using the values c 1 = 2, c 2 = −1, λ b = 0.2531, λ a = 1.4874, z = 0.3074 and α(T r ) = 1. Equation (22) takes the form:   Let us solve Peng-Robinson equation of state for Oxygen O 2 at T = T c = 153, P = P c = 50 and using these values we get following nonlinear equation: where R = 0.0831446261815324. Solving for τ = V , gives us the multiple root w = 0.07820926381 with multiplicity η = 3. By using the initial guess τ 0 = 0.9, the numerical results are presented in Table 6.
In this last example, we solve Peng-Robinson equation of state for oxygen O 2 and the basins of attraction obtained for the methods SH1-SH4 and F Z1− F Z4 are shown in Figs. 11, 12 We use the region [−2, 2] × [−2, 2] to draw the basins of attraction and examine that SH1-SH4 and FZ1-FZ4 use minimum two and maximum three numbers of iterations to converge to the root. In this example, the black color indicates the hue of two iterations for all iterative schemes. We observe that newly proposed scheme has globally wider black region as compared to recent methods.

Conclusion
There exist many numerical root-solvers with high order of convergence that require derivative evaluations for finding the multiple root of a function. The derivative free higher-order techniques for computing multiple root are rare and yet to be explored. It is not an easy task to attain an optimal derivative-free root-finder family of methods. We have presented a derivative-free optimal eighth order iterative scheme for finding the zeros of a nonlinear equations with known multiplicity of zeros with a univariate and a multivariate weight functions. The number of functional evaluations required for our presented scheme is four. We have compared the numerical results of our presented scheme with the most recent methods in the literature. We conclude that our method gives eight order of convergence also for the case when function has root of multiplicity one. The dynamical comparison of our presented method with recent seventh and eighth order methods shows that newly proposed method has, in general, wider convergence region. It is apparent from the numerical results and construction that our presented family of method is optimal and efficient in terms of small residual errors. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.