Family of fourth-order optimal classes for solving multiple-root nonlinear equations

We present a new iterative procedure for solving nonlinear equations with multiple roots with high efficiency. Starting from the arithmetic mean of Newton’s and Chebysev’s methods, we generate a two-step scheme using weight functions, resulting in a family of iterative methods that satisfies the Kung and Traub conjecture, yielding an optimal family for different choices of weight function. We have performed an in-depth analysis of the stability of the family members, in order to select those members with the highest stability for application in solving mathematical chemistry problems. We show the good characteristics of the selected methods by applying them on four relevant chemical problems.


Introduction
Throughout history, physical phenomena have been modeled through mathematical expressions. On many occasions, it is necessary to obtain the solution of the nonlinear equation f (x) = 0, f : D ⊆ R → R. However, it is not always possible to obtain an analytical solution to such a problem and we have to resort to approximate solutions. Iterative methods obtain approximations to the solution as accurate as we need. The best known and most widely used scheme is Newton's method Under the convergence conditions -initial estimate close to the root and sufficiently differentiable function -Newton's method converges quadratically as long as we look for a simple root. But there are many physical phenomena in which we have to look for a multiple root, and there the convergence is impaired. Nonlinear equations f (x) = 0 with multiple roots of multiplicity m > 1 satisfy f (r m ) = f (r m ) = · · · = f (m−1) (r m ) = 0 and f (m) (r m ) = 0. Many authors have designed iterative methods for this purpose [1][2][3][4][5][6]. Kansal et al. [7] designed two-step optimal methods of fourth order, starting from the arithmetic mean between second order methods for multiple roots, and including weight functions and accelerating parameters in the final scheme. Behl et al. [8,9] followed a similar strategy to reach order of convergence four.
In this paper we propose the design and analysis of an iterative class taking the arithmetic mean of the accelerated Newton's method [10] x n+1 = x n − m f (x n ) f (x n ) , n = 0, 1, 2, . . . and Chebyshev's iterative scheme [11] x 3 , n = 0, 1, 2, . . . , obtaining the following iterative method: where Our aim is to design an optimal method in the sense of the Kung-Traub conjecture [11] which states that the order of convergence p of an iterative scheme without memory is at most 2 d−1 , where d is the number of different functional evaluations performed by the iterative algorithm, reaching the optimality when p = 2 d−1 . Moreover, Traub showed in [11] that in order to design a one-step method with order of convergence p, its iterative expression must contain derivatives at least up to order p − 1.
Let us notice that (1) is a one-step iterative method that uses three different functional evaluations on each iterative step, being one of them a second-order derivative of the function. With these features, its order of convergence is at most three, so this method cannot reach the optimality. For this reason, we add an extra step in (1) to obtain a two-step iterative scheme that does not require derivatives of order two and uses only three functional evaluations so that according to the Kung-Traub conjecture it could be possible to reach the optimal order of convergence 4.
The paper is structured as follows. Section 2 analyzes the convergence of the introduced family. Section 3 studies the stability of the family, in order to find its best members in terms of initial guesses. In Sect. 4 a numerical analysis is performed, showing the validity of the method in chemistry and academic problems. Finally, Sect. 5 collects the main conclusions of the manuscript.

Design of the fourth-order optimal family
Consider the Newton-type iterative method for multiple roots Using Taylor series developments around x = x n , we have: From (2), the second derivative of the function can be approximated as f (x n ) . In addition, we can write Replacing (3) in (1), we obtain It can be proved the linear convergence of (4). Also, although the derivative of order two is not used, it is not an optimal method because three functional evaluations are still carried out. Based on its iterative structure we design another method with higher order of convergence. In this sense, we add two free parameters a and b in (4). Then the biparametric family is The convergence of family (5) is analyzed below.
Theorem 1 Let r m be a multiple zero with a muliplicity m ≥ 1 of a sufficiently differentiable function f : I ⊆ R −→ R defined in an open interval I such that r m ∈ I . If the initial estimation x 0 is close enought to r m , then the iterative scheme defined by (5) has order of convergence three when parameters a and b are defined as following: In this case, the error equation of the resulting method is given by , and e n = x n − r m , ∀n ∈ N.
Proof Using Taylor series developments around x = r m and taking into account and its derivative From (6) and (7), we obtain Then we have and then can be written as where β = 3m 4 (−48 − 24m + 36m 2 + 30m 3 + 9m 4 + m 5 )c 3 . Now, from (7) and (9) we have Using (8) and (10) on the second step of the iterative scheme we obtain the error equation From the error equation, the iterative class (5) can achieve cubic order when it holds Equivalently, the order is cubic when parameters a and b take the value and then the error equation of the method results into In order to construct a fourth-order optimal method with less than three different functional evaluations, we propose the introduction of a weight function on (5) obtaining n = 0, 1, 2, . . . , (11) where t n = and Q ∈ C 2 (R) is a real variable weight function. The following result shows the conditions that the weight function Q must satisfy in (11) to reach the optimal order of convergence four. If the initial estimation x 0 is close enought to r m , then the iterative scheme defined by (11) has order of convergence four when the weight function Q satisfies: Under these conditions, the following error equation is hold: Proof From (6) and (7), we can write Now we can , so we consider the Taylor series expansion of the Using (8), (10) and (13) in the iterative expression (11), the error equation holds where To obtain an optimal family of iterative methods of order 4, the coefficients of e n , e 2 n and e 3 n in the error Eq. (14) must be zero simultaneously. Solving for K 1 = 0, K 2 = 0 and K 3 = 0 we obtain the following conditions to cancel the terms up to order three of Eq. (14) Using the previous conditions, the error equation of family (11) shows that the order of convergence is four: According to Theorem 2, we have designed a two-step family of iterative methods that only performs three functional evaluations on each iteration and reaches the optimal order of convergence p = 4.
In the following, we will try to simplify the iterative expression of (11). First, using the notation u n = f (y n ) f (x n ) we can write and then On the other hand, from Theorem 2, the weight function can be written as family (11) can be simplified as being α ∈ C a free parameter. A new family of optimal iterative methods belonging to (11) is designed. Let us denote the iterative class (16) by NC4 family. Next section is devoted to perform a dynamical analysis of NC4 family in order to choose the best members in terms of stability.

Stability of the family of methods
The stability analysis of a family of iterative methods studies the behaviour of the schemes for a wide set of initial guesses. In this sense, this analysis discriminates whether the methods are useful to solve the nonlinear problems or not. In section 3.1 the scaling theorem is introduced, and the Möbius transformation is applied to reduce the amount of parameters involved. Section 3.2 analyzes fixed and critical points. Finally, section 3.3 is devoted to select the members of the family with better stability behaviour.
Complex dynamics is used to study the dynamical behaviour of the rational operator associated to family NC4 applied on polynomials. The concepts can be reviewed in more detail at [12][13][14][15].

Operator simplification through Möbius transformation
The iterative family is applied for solving the general nonlinear multiple-root polyno- . Theorem 3 points that family NC4 satisfies the scaling theorem, in order to simplify the stability analysis.

Theorem 3 (Scaling theorem for family NC4) Let f be an analytic function inĈ, and let T
that is, O g and O f are affine conjugated by T .
Since family (16) satisfies the scaling theorem, the fixed point operators associated to family NC4 applied to analytic functions are affine conjugated by an affine map.
The rational operator of family NC4 applied on where 3 , Note that (17) depends on the variable z, the parameter α and the roots δ 1 and δ 2 . In order to overcome the dependence on the roots, we consider the Möbius transformation [12] M obtaining its affined conjugated operator Since M(δ 1 ) = 0, M(δ 2 ) = ∞ and M(∞) = 1, Möbius transformation maps roots δ 1 and δ 2 with 0 and ∞, respectively, while the divergent behaviour will be at 1.

Selection of the best members of the family
Selecting a value of α in the black region of Figure 2 or 4 results in a set of initial guesses that converge to a strange attracting fixed point. Figure 7shows two samples. Figure 7(a) represents the dynamical plane of the method α = −400 + 200i. A wide set of initial estimations converge to the root 0. However, some initial guesses converge to the strange attracting fixed point t 1 ≈ 11.4574 represented in white star. A similar behaviour is represented in Figure 7(b), for the method corresponding to α = −1500, where some initial estimations converge to the strange attracting fixed point s 1 ≈ 7.9661, represented in white star.
Free critical points can have their own basin of attraction, different from the roots of the polynomial. Selecting a value of α in the black region of Figure 6 results in iterative methods whose free critical points converge to their own root. Figure 8illustrates this behaviour. Figure 8(a) represents the dynamical plane of the method α = 0. In this case, u 1 ≈ 26.8431 -represented in white square-is inside a basin of attraction of a twoperiodic orbit. Figure 8(b) is the corresponding dynamical plane to α = −119 + 157i.

Numerical tests
In this section we compare the introduced method NC4 with relevant methods of the literature. The selected methods can be found in Sharma and Sharma [17] denoted as SS, Behl et al [9] denoted as B4, and Sharifi et al. [18] denoted as SH. The comparison is performed on different chemical nonlinear equations that involve multiple roots. For The numerical results of each problem are displayed in a table. In particular, each method is applied until | f (x k+1 )| < tol or |x k+1 − x k | < tol. If it is not stated differently, tol = 10 −300 is used along the following examples. Both values are shown in the result tables for the last iterate computed. If the targeted solution is available, we introduce the value of |x k+1 − α|. Furthermore, we include the number of iterates (it) required to reach those tolerances and the approximated computational order of convergence ACOC [19].

Non-ideal gas model
The Van der Waals gas equation describes the evolution of a non-ideal gas from its idealised version. It uses two parameters a 1 and a 2 to study the nonideality of the gas.
The equations are described as follows: where P is the pressure of the gas, V is its volume, n are the moles of gas, R is the universal constant of ideal gases and T is the absolute temperature of the gas.
Taking proper values for n, P y T and the constants a 1 and a 2 we end up with equation This is a polynomial with a double root at 1.75 and a simple root at 1.72. Table 1shows the results for the Non-ideal gas model.
We observe that for the parameter −1500 the NC4 method does not converge, this can be expected since the parameter is in the region where at least one of the strange fixed points is attracting.

Stirred-tank reactor
Let us consider a stirred tank reactor in which an isothermal fluid is stirred continuously, given in Constantinides and Mostoufi [20].
Two components A and R are injected in the reactor at ratios Q and q − Q respectively. The following reactions occurs in the reactor [21] proposes a feedback control system in order to control the velocity of the reaction. After his analysis the equation of the transfer function of the reaction is where K c is the proportional gain of the controller in the control system. The control is stable for those values of K c with negative real part of the transfer function. In control theory, the locatoin the poles of the transfer function improves the knowledge of the behaviour of the controller. Therefore, we are interested finding the roots of equation f 2 (x) = x 4 + 11.5x 3 + 47.49x 2 + 83.06325x + 51.23266875 = 0.
Here, x = −2.85 is a root of multiplicity 2. There are two other roots located at x = −1.45 and x = −4.35. The results are displayed in Table 2. Table 1 Numerical results for the Non-ideal gas model   Table 2 Numerical results for the Stirred-tank reactor model There is no convergence of the method. * *

Converges to a simple root
Two different initial conditions have been considered. The second one is far from the root, and therefore, several of the methods do not converge. Three of the parameters tested are able to converge to the desired solution. We observe that the SS method converges to a simple root, yielding an undesired result.

Oceanic acidity
The C O 2 concentration of the ocean is modelled according to McHugh et al. [22] and developed by Babajee [23] and Kansal et al [24].
The model computes the acidity levels of the ocean by computing the roots of a fourth order polynomial. The hypothesis considered by Babajee [23] in order to simplify the problem are The values K 0 , K 1 , K 2 , K W and K B are the equilibrium constants of the reactions that involve the acidification process. The parameter A represents the alcalinity of the ocean water, P t is the partial pressure of C O 2 .
In [24] the values A = 2.050 and B = 0.409 are taken from [25,26] and obtained from [23]: P t = 200, K 0 = 3.347(−5), K 1 = 9.747(−4), K 2 = 8.501(−7), K W = 6.46(−9), K B − > 1.881(−6), The values of the constants K i can be extracted from [23]. The value of P t is a variable that we set in P t = 148.508. Under these conditions, we obtain the polynomial that has a double root at x = 178.977 and two simple roots at x = −360.003 y x = 11.2859. Table 3shows the numerical results to find the zero of multiplicity 2 for taht equation. In that case the tolerance required is 10 −200 .

Fixed points in a bi-electron model
A classical problem in chemistry is study of the movement of an electron in the hydrogen atom with a circularly polarised microwave field [27].
The dynamics of this model are given by the Hamiltonian where K is the intensity of the microwave field. If a a negative charged nucleus is considered instead of positive one, the resulting model is the bi-electron model One is interested in searching in the fixed points of that model. Once computed the dynamical equations it is clear that y = 0 and x must verify: Considering K = 2 this equation has one double root for x = 1 and one single root at x = 0. Table 4 shows the results for the bi-electron model.

Conclusion
In the present paper we have introduced a parametric family of order 4 to find roots of multiplicity greater than 2 of nonlinear equations. In order to select the most stable members of the class, a complete study on the stability and the parameters of the Table 3 Numerical results for the ocean acidity model          There is no convergence of the method Table 4 Numerical results for the ocean acidity model                     There is no convergence of the method family has been performed. The application of these schemes is supported with several numerical experiments in the field of mathematical chemistry, showing the good properties of the introduced methods.