Design and Complex Dynamics of Potra–Pták-Type Optimal Methods for Solving Nonlinear Equations and Its Applications

: In this paper, using the idea of weight functions on the Potra–Pták method, an optimal fourth order method, a non optimal sixth order method, and a family of optimal eighth order methods are proposed. These methods are tested on some numerical examples, and the results are compared with some known methods of the corresponding order. It is proved that the results obtained from the proposed methods are compatible with other methods. The proposed methods are tested on some problems related to engineering and science. Furthermore, applying these methods on quadratic and cubic polynomials, their stability is analyzed by means of their basins of attraction.


Introduction
For solving nonlinear equations iteratively, the Newton's method given by is one of the most commonly used methods. The efficiency index as defined by Ostroswki in [1], which relates the order of convergence of a method p with the number of function evaluations per iteration d, is given by the expression p 1/d . Newton's method is quadratically convergent and requires two function evaluations per iteration and, thereby, has the efficiency index value of 2 1/2 ≈ 1.414. Numerous methods have appeared giving higher order of convergence or better efficiency. One of the recent strategies to increase the order of the methods is the use of weight functions [2][3][4][5]. In this regard, Sharma and Behl [6] presented the fourth order method: Similarly, Sharifi et al. [7] used weight functions on the third order Heun's method and proposed the fourth order method f (x n ) , According to Kung and Traub [8], an iterative method is said to be optimal if its order is 2 d−1 , where d is the number of function evaluations per iteration. Notice that Newton's method as well as (1) and (2) are all optimal. Potra and Pták [9], as an attempt to improve Newton's method, gave the method ( This method is cubically convergent but is not optimal, as it requires three function evaluations per iteration. The aim, in the present paper, is to further investigate the method (3). Precisely, we use weight functions and improve the order of convergence of (3). We do it in three ways which correspond to the methods of orders 4, 6 and 8. Out of these, the methods with orders 4 and 8 are optimal.
Dynamics of a rational operator give important information about the convergence, efficiency and stability of the iterative methods. During the last few decades, many researchers, e.g., [10][11][12][13][14][15][16] and references therein, study the dynamical behavior of rational operators associated with iterative methods. Furthermore, there is an extensive literature [17][18][19][20][21] to understand and implement further results on the dynamics of rational functions. In this paper, we also analyze the dynamical behavior of the methods that we have developed in this paper. Furthermore, at the end of this work, the basins of attraction are also presented and compared among the proposed and other methods of the corresponding order.
The remaining part of the paper is organized as follows. In Section 2, the development of the methods and their convergence analysis are given. In Section 3, the proposed methods are tested on some functions, and the results are compared with other methods in the head of Numerical Examples. In Section 4, the proposed methods are tested on some engineering and science related designs. Section 5 is devoted to analyze the stability of the introduced methods by means of complex dynamics. In this sense, the study of the rational function resulting from the application of the methods to several nonlinear functions is developed, and their basins of attraction are represented. Finally, Section 6 covers the conclusions of the research.

Development of Methods and Their Convergence Analysis
In this section, the methods of order four, six and eight are introduced, and its convergence is analyzed.

Optimal Fourth Order Method
Based on the Potra-Pták method (3), we propose the following two-step method using a weight function, whose iterative expression is where w(t n ) = a 1 + a 2 t n + a 3 t 2 n and t n = f (y n ) f (x n ) . The convergence of (4) is proved in the following theorem.
Theorem 1. Let f be a real or complex valued function defined in the interval I having a sufficient number of smooth derivatives. Let α be a simple root of the equation f (x) = 0 and the initial point x 0 is close enough to α. Then, the method (4) is fourth order of convergence if a 1 = 1, a 2 = 0 and a 3 = 2.
Let d n = y n − α, then, from the first equation of (4), we get so that, using Taylor's series expansion of f (y n ) about α, we get Now, from (5) and (7), we get Therefore, using the results obtained above in the second equation of (4), we get In order to obtain fourth order of convergence, in view of (9), we must have In order to obtain sixth order of convergence, the coefficients of e 4 n and e 5 n must vanish in (14), i.e., b 1 = 1 and b 2 = 2. Therefore, the error equation of the method (11) becomes e n+1 = c 2 18c 4 2 − 9c 2 2 c 3 + c 2 3 e 6 n + O e 7 n , and the assertion follows.
In view of Theorem 2, the following is the sixth order method

Optimal Eighth Order Method
Notice that the method (15) is not optimal as it requires four function evaluation per iteration to achieve sixth order of convergence. Its efficiency index is 1.5651, which is less than that of the fourth order method (10). However, an eighth order method is obtained by (10) using an additional Newton step. The resulting iterative scheme is Nevertheless, this method requires five function evaluation per iteration, so that its efficiency index reduces to 1.5157, and, moreover, it is not optimal. Towards making the method (16) more efficient and optimal, we approximate f (z) as where Here, J and G are some appropriate weight functions of two variables and one variable, respectively. This type of approximations was done by Matthies et al. in [22]. Accordingly, we propose the following method: where t n , u n , and s n , are as in (17). For the method (18), we take the functions J and G as J(t n , u n ) = 1 + 2t n + (β + 2)u n + 3t 2 and where β and λ belong to C. We prove the following result.
Theorem 3. Let f be a real or complex valued function defined on some interval I having a sufficient number of smooth derivatives. Let α be a simple root of the equation f (x) = 0 and the initial point x 0 is close enough to α. Then, (18) is an eighth order of convergence for the functions J and G given by (19) and (20), respectively.

Numerical Examples
In this section, we test the performance of the methods proposed in Section 2 with the help of some numerical examples. We compare the results obtained with the known methods of the corresponding order. We consider the following nonlinear equations and initial guesses: In the previous section, we have proved the theoretical order of convergence of various methods. For practical purposes, we can test numerically the order of convergence of these methods by using Approximated Computational Order of Convergence (or ACOC), defined by Cordero and Torregrosa [23]. They defined the ACOC of a sequence {x k }, k ≥ 0 as The use of ACOC, given by (23), serves as a practical check on the theoretical error calculations. We apply our proposed methods and other existing methods as discussed in the following subsections on each of the test functions. Various results of up to four iterations are observed, and we compare the results obtained at the 4th iteration among different methods of the corresponding order and shown in Tables 1-3. For a particular test function, we take the same initial guess x 0 for each of the methods under consideration. We compare the approximate error ∆x n ≡ |x n − x n−1 |, the approximate solution x n , the absolute value of corresponding functional value | f (x n )|, and approximated computational order of convergence (ACOC) at n = 4. In the tables, "NC" stands for no convergence of the method. We use Mathematica 9.0 for the calculations.

Comparison of the Fourth Order Method
Let us denote our method (10)  Jarratt's method [24], denoted by M 44 and given by Kung-Traub [8] method, denoted by M 45 , and given by All the methods M 4i , i = 1, 2, 3, 4, 5 are optimal. Table 1 records the performance of all these methods.

Comparison of Sixth Order Methods
We denote our sixth order method (15) by M 61 . We shall compare this method with • M 62 : Method of Neta [25] with a = 1, given by M 63 : Method of Grau et al. [26] given by M 64 : Method of Sharma and Guha [27] with a = 2, given by M 65 : Method of Chun and Neta [28] given by The comparison of the methods M 6i , i = 1, 2, 3, 4, 5 is tabulated in Table 2. From the table, we observe that the proposed method M 61 is compatible with the other existing methods. We can see that method M 63 gives different results for the test functions f 2 and f 5 with given initial guesses.

Comparison of Eighth Order Methods
Consider the eighth order method (18), which involves the parameter pair (β, λ). We denote • M 81 the case where (β, λ) = (0, 0), whose iterative expression results in M 82 for (β, λ) = (1, 1), resulting in the iterative scheme given by M 81 : M 83 for (β, λ) = (0, 1), whose iterative method is Along with these, we take the following methods for the comparison of numerical results: • Matthies et al. in [22] presented an optimal class of 8th order method from the Kung-Traub method [8]. For some particular values of the parameters, one of the methods denoted by M 84 is given by 2+t n +5u n +4t 2 n +4t 3 n 2−3t n +u n +2t 2 n · 2+s n 2−s n .
• Babajee et al. in [11] presented a family of eighth order methods. For some fixed values of parameters, the method denoted by M 85 is given by Chun and Lee in [29] presented a family of optimal eighth order methods. For some particular values of parameters, the method denoted by M 86 is given by In all the above methods, t n , u n and s n are as given in (17). The performance of the methods M 8i , i = 1, 2, 3, 4, 5, 6 are recorded in Table 3.   From Tables 1-3, we observe that the proposed methods are compatible with other existing methods (and sometimes perform better than other methods) of the corresponding order. Not any particular method is superior to others for all examples. Among the family of eighth order methods (18), from Table 3, we observe that the method M 81 performs better than other two. For more understanding about the iterative methods, we study the dynamics of these methods in the next section.

Applications
The applications discussed in Sections 4.1-4.3 are based on standard engineering examples, and we refer to [30]. We use the proposed methods M 41 , M 61 , and M 8i , i = 1, 2, 3 to obtain the various results from the first three iterations of these examples. In particular, we compute the value of the unknowns x n−1 and x n , absolute value of the function f (x n ) and absolute value of the difference d of unknown in two consecutive iterations, i.e., d = |x n − x n−1 |, n = 1, 2, 3.

Pipe Friction Problem
Determining fluid flow through pipes and tubes has great relevance in many areas of engineering and science. In engineering, typical applications include the flow of liquids and gases through pipelines and cooling systems. Scientists are interested in topics ranging from flow in blood vessels to nutrient transmission through a plant's vascular system. The resistance to flow in such conduits is parameterized by a dimensionless number called the friction factor f . For a flow with turbulence, the Colebrook equation [31] provides a means to calculate the friction factor: where is the roughness (m), D is the diameter (m) and Re is the Reynolds number Here, ρ denotes the fluid density (kg/m 3 The results obtained by the various methods are presented in Table 4.

Open-Channel Flow
An open problem in civil engineering is to relate the flow of water with other factors affecting the flow in open channels such as rivers or canals. The flow rate is determined as the volume of water passing a particular point in a channel per unit time. A further concern is related to what happens when the channel is slopping.
Under uniform flow conditions, the flow of water in an open channel is given by Manning's equation where S is the slope of the channel, A is the cross-sectional area of the channel, R is the hydraulic radius of the channel and n is the Manning's roughness coefficient. For a rectangular channel having the width B and the defth of water in the channel y, it is known that With these values, (26) becomes Now, if it is required to determine the depth of water in the channel for a given quantity of water, (27) can be rearranged as In our work, we estimate y when the remaining parameters are assumed to be given as Q = 14.15 m 3 /s, B = 4.572 m, n = 0.017 and S = 0.0015. We choose as an initial guess y 0 = 4.5 m. The results obtained by the various methods are presented in Table 5.

Ideal and Non-Ideal Gas Laws
The ideal gas law is where P is the absolute pressure, V is the volume, n is the number of moles, R is the universal gas constant and T is the absolute temperature. Due to its limited use in engineering, an alternative equation of state for gases is the given van der Waals equation [32][33][34][35] where v = V n is the molal volume and a, b are empirical constants that depend on the particular gas. The computation of the molal volume is done by solving We take the remaining parameters as R = 0.082054 L atm/(mol K), for carbon dioxide a = 3.592, b = 0.04267, T = 300 K, p = 1 atm, and the initial guess for the molal volume is taken as v 0 = 3. The results obtained by the various methods are presented in Table 6. In this table, I ND stands for indeterminate form.

Dynamical Analysis
The stability analysis of the methods M 41 , M 61 and M 8i , i = 1, 2, 3, is performed in this section. The dynamics of the proposed methods on a generic quadratic polynomial will be studied, analyzing the associated rational operator for each method. This analysis shows their performance depending on the initial estimations. In addition, method M 41 is analyzed for cubic polynomials. First, we recall some basics on complex dynamics.

Basics on Complex Dynamics
Let R :Ĉ −→Ĉ be a rational function defined on the Riemann sphere. Let us recall that every holomorphic function from the Riemann sphere to itself is in fact a rational function R(z) = P(z) Q(z) , where P and Q are complex polynomials (see [36]). For older work on dynamics on the Riemann sphere, see, e.g., [37].
The orbit of a point z 0 ∈Ĉ is composed by the set of its images by R, i.e., Note that the roots z * of an equation f (z) = 0 are fixed points of the associated operator of the iterative method. Fixed points that do not agree with a root of f (x) = 0 are strange fixed points. The asymptotical behavior of a fixed point z F is determined by the value of its multiplier µ = |R (z F )|. Then, z F is attracting, repelling or neutral if µ is lower, greater or equal to 1, respectively. In addition, it is superattracting when µ = 0. For an attracting fixed point z F , its basin of attraction is defined as the set of its pre-images of any order: A(z F ) = {z 0 ∈Ĉ : R n (z 0 ) −→ z F , n → ∞}.
The dynamical plane represents the basins of attraction of a method. By iterating a set of initial guesses, their convergence is analyzed and represented. The points z C ∈Ĉ that satisfy R (z C ) = 0 are called critical points of R. When a critical point does not agree with a solution of f (x) = 0, it is a free critical point. A classical result [21] establishes that there is at least one critical point associated with each immediate invariant Fatou component.

Rational Operators
Let p(z) be a polynomial defined onĈ. Corresponding to the methods developed in this paper, i.e., methods (10), (15) and family (18), we define the operators R 4 (z), R 6 (z) and R 8 (z), respectively, inĈ as follows: . Thus, it proves that (R Theorem 4 allows for generalizing the dynamical study of a specific polynomial to a generic family of polynomials by using an affine map. Analogous to the way we proved the Scaling Theorem for the operator R 4 , it also follows that the fixed point operators R 6 and R 8 obey the Scaling Theorem.
Regarding the critical points of S 6 , the roots of p(z) are critical points, and S 6 has the free critical points z C 3 = −1 and the roots of a tenth-degree polynomial, z C 4 , . . . , z C 11 . The dynamical planes are a useful tool in order to analyze the stability of an iterative method. Taking each point of the plane as initial estimation to start the iterative process, they represent the convergence of the method depending on the initial guess. In this sense, the dynamical planes show the basins of attraction of the attracting points. Figure 1 represents the dynamical planes of the methods S 4 and S 6 . The generation of the dynamical planes follows the guidelines established in [38]. A mesh of 500 × 500 complex values has been set as initial guesses in the intervals −5 < {z} < 5, −5 < {z} < 5. The roots z * 1 = 0 and z * 2 = ∞ are mapped with orange and blue colors, respectively. The regions where the colors are darker represent that more iterations are necessary to converge than with the lighter colors, with a maximum of 40 iterations of the methods and a stopping criteria of a difference between two consecutive iterations lower than 10 −6 .
As Figure 1 illustrates, there is convergence to the roots for every initial guess. Let us remark that, when the order of the method increases, the basin of attraction of z * 1 = 0 becomes more intricate. Finally, for the fixed point operators associated with family M 8 , the solutions of S 8i (z) = z for i = 1, 2, 3 give the superattracting fixed points z F 1,2 = z * 1,2 and the repelling point z F 3 = 1. In addition, S 81 has 28 repelling points. S 82 and S 83 have 38 repelling points, corresponding to the roots of polynomials of 28 and 38 degree, respectively, and the strange fixed points z F 4,5 = 1 2 (−1 ± √ 5). The number of critical points of the fixed point operators S 8i are collected in Table 7. In addition, the number of strange fixed points and free critical points are also included in the table for all of the methods.  Figure 2 represents the dynamical planes of the methods S 81 , S 82 and S 83 . Since the original methods satisfy the Scaling Theorem, the generation of one dynamical plane involves the study of every quadratic polynomial. There is an intricate region around z = −1 in Figure 2a, becoming wider in Figure 2b,c around z = −1.5. However, for every initial guess in the three dynamical planes of Figure 2, there is convergence to the roots.

Dynamics on Cubic Polynomials
The stability of method M 41 on cubic polynomials is analyzed below. As stated by the authors in [39], the Scaling Theorem reduces the dynamical analysis on cubic polynomials to the study of dynamics on the cubic polynomials p 0 (z) = z 3 , p + (z) = z 3 + z, p − (z) = z 3 − z and the family of polynomials p γ (z) = z 3 + γz + 1. Let us recall that the first one only has the root z * 1 = 0, while p + (z) and p − (z) have three simple roots: z * 1 = 0 and z * 2,3 = ∓i or z * 2,3 = ∓1, respectively. For each γ ∈ C, the polynomial p γ (z) also has three simple roots that depend on the value of γ. They will be denoted by z * 1,2,3 (γ).
By applying method M 41 to polynomials p 0 (z), p + (z) and p − (z), the fixed point operators obtained are, respectively, Moreover, there is the presence of free critical points with values z C 4,5 = ±i 5 23 for S 4,+ (z) and z C 4,5 = ± 5 23 for S 4,− (z). As for quadratic polynomials, the dynamical planes of method M 41 when it is applied to the cubic polynomials have been represented in Figure 3. Depending on the roots of each polynomial, the convergence to z * 1 = 0 is represented in orange, while the convergence to z * 2 and z * 3 is represented in blue and green, respectively. It can be see in Figure 3 that there is full convergence to a root in the three cases. However, there are regions with darker colors that indicate a higher number of iterations until the convergence is achieved. When method M 41 is applied on p γ (z), the fixed point function turns into The fixed points of S 4,γ (z) are the roots of the polynomial z * 1,2,3 (γ), being superattracting, and the strange fixed points z F 4−9 (γ) that are the roots of the sixth-degree polynomial q(z, γ) = 35z 6 + 37γz 4 + 7z 3 + 11γ 2 z 2 + γz + γ 3 − 1.
As the asymptotical behavior of z F 4 (γ), . . . , z F 9 (γ) depends on the value of γ, the stability planes corresponding to these points are represented in Figure 4. For each strange fixed point, a mesh of 100 × 100 points covers the values of (γ) ∈ [−5, 5] and (γ) ∈ [−5, 5]. The stability plane shows the values for the parameter where |S 4,γ (z F )| is lower or greater than 1, represented in red or green, respectively. The solutions of S 4,γ (z) = 0 are the critical points z C 1,2,3 (γ) = z * 1,2,3 (γ) and the free critical points z C 4 = 0 and When the fixed point function has dependence on a parameter, another useful representation is the parameters' plane. This plot is generated in a similar way to the dynamical planes, but, in this case, by iterating the method taking as an initial guess a free critical point and varying the value of γ in a complex mesh of values, so each point in the plane represents a method of the family. The parameters' plane helps to select the values for the parameter that give rise to the methods of the family with more stability. The parameters' planes of the four free critical points are shown in Figure 5. Parameter γ takes the values of 500 × 500 points in a complex mesh in the square [−5, 5] × [ −5, 5]. Each point is represented in orange, green or blue when the corresponding method converges to an attracting fixed point. The iterative process ends when the maximum number of 40 iterations is reached, in which case the point is represented in black, or when the method converges as soon as, by the stopping criteria, a difference between two consecutive iterations lower than 10 −6 is reached.
For the parameters' planes in Figure 5, there is not any black region. This guarantees that the corresponding iterative schemes converge to a root of p γ (z) for all the values of γ.
In order to visualize the basins of attraction of the fixed points, several values of γ have been chosen to perform the dynamical planes of method M 41 . These values have been selected from the different regions of convergence observed in the parameters planes. Figure 6, following the same code of colours and stopping criteria as in the other representations, shows the dynamical planes obtained when these values of γ are fixed. As Figure 6 shows, there is not any initial guess that tends to a point different than the roots. This fact guarantees the stability of these methods on the specific case of any cubic polynomial.

Conclusions
Two iterative schemes of orders of convergence four and six, and a family of methods of order eight have been introduced. The method of order four and the family of order eight are optimal in the sense of Kung-Traub's conjecture. The development of the order of convergence of every method has been performed. For every method, we have made a numerical experiment, over both test functions and real engineering problems. In order to analyze the stability of the introduced methods, the dynamical behavior of them has been studied. The results confirm that the methods have wide basins of attraction, guaranteeing the stability over some nonlinear problems.