Development of a Family of Jarratt-Like Sixth-Order Iterative Methods for Solving Nonlinear Systems with Their Basins of Attraction

We develop a family of three-step sixth order methods with generic weight functions employed in the second and third sub-steps for solving nonlinear systems. Theoretical and computational studies are of major concern for the convergence behavior with applications to special cases of rational weight functions. A number of numerical examples are illustrated to confirm the convergence behavior of local as well as global character of the proposed and existing methods viewed through the basins of attraction.

This paper is devoted to devise a class of sixth-order iterative root-finders for nonlinear systems by employing a three-step weighted Jarratt-like method below: where s = f (x n ) −1 f (y n ), γ is a parameter to be determined later and T f , L f : C → C are weight functions being analytic [16][17][18] in a neighborhood of 1. Note that Scheme (1) uses two functional values as well as two derivatives. We are certainly able to introduce generic weight functions using one derivative and three functional values to develop general optimal eighth-order methods that covers the existing ones for the zero of a given scalar function. However, expanding such approach to a nonlinear system requires different weight functions. For unified analysis to be performed in both scalar and vector functions, we aim to develop a family of Jarratt-like sixth-order iterative methods by maintaining the same form of weight functions with two derivatives as well as two functional values. This extension to nonlinear systems is the main strength of this paper.
The robustness of the current analysis presented here covers most existing studies on higher-order root-finders using two derivatives and two function values for both scalar and vector equations. The results of Theorem 1 give us not only fairly generic scalar function solvers, but also some advantage of extending to a nonlinear system with any finite dimension. Such an extension is evidently characterized by Theorem 2 to be studied in this analysis.
Our major aim is not only to design a class of sixth-order methods by fully specifying the algebraic structure of generic weight functions T f (s) and L f (s), but also to investigate their basins of attraction behind the extraneous fixed points [19] when applied to polynomials. The last sub-step of (1) in the form of weighted Newton's method is clearly more convenient in dealing with extraneous fixed points which are the roots of the weight function T f (s) + L f (s) · f (z) f (x) . The extraneous fixed points may lead us to attractive, indifferent, repulsive and chaotic orbits via the related basins of attraction.
Section 2 investigates the main theorem regarding the convergence behavior with the desired forms of weight functions, while Section 3 deals with special cases of weight functions that can cover many of the existing studies using two derivatives and two functional evaluations. Section 4 discusses the computational and long-term orbit behavior of the proposed iterative methods regarding scalar functions. Section 5 presents numerical experiments in a d-dimensional Euclidean space by solving a system of nonlinear vector equations f : R d → R d encountered in a real life with d ∈ {3, 4, 9, 10}. In addition, computational efficiency is addressed with issues related to the accuracy and applicability of the proposed methods. Concluding remarks are stated in Section 6.

Main Theorem
The main theorem for a nonlinear scalar equation will be pursued and extended later in Section 5 to a system of nonlinear vector equations: Theorem 1. Assume that f : D ⊂ C → C has a simple root α and is analytic in a region D containing α.
Let θ j = f (j) (α) j! f (α) for j = 2, 3, · · · . Let x 0 be an initial guess chosen in a sufficiently small neighborhood of α. Let T f , L f : C → C be analytic in a neighborhood of 1.
In view of (20) and (21), we can select a total of 154 special pairs of second-order rational weight functions (T f (s), L f (s)). Excluding known studies, the following pairs of weight functions (T f (s), L f (s)) may be of great interest to us. The corresponding methods to such pairs of weight functions (T f (s), L f (s)) are denoted by LK i for 1 ≤ i ≤ 10, respectively, and indicated on the right of (22) and (23).
One should be aware that Method LK 2 can be found by taking α = β = 2 from (16).

Computational Experiments on Local and Global Convergence
For computational experiments, we first deal with local convergence of methods (1) for a variety of test functions along with the existing methods EM1-EM7; then we discuss global convergence underlying extraneous fixed points via basins of attraction. Numerical experiments for 17 methods EM1-EM7 and LK1-LK10 were implemented with Mathematica with 300 and 140 digits of minimum number of precision for scalar and vector equations, respectively.

Definition 1. (Computational Convergence Order (COC) and Approximated Computational Convergence
Order (ACOC)) Assume that theoretical asymptotic error constant η = lim n→∞ |e n | |e n−1 | p and convergence order p ≥ 1 are known. Define p n = log |e n /η| log |e n−1 | as the computational convergence order. Note that lim n→∞ p n = p. Approximated computational convergence orderp n is defined asp n = log(|x n −x n−1 |/|x n−1 −x n−2 |) log(|x n−1 −x n−2 |/|x n−2 −x n−3 |) requiring knowledge at four points x n , x n−1 , x n−2 , x n−3 . Table 1 lists test scalar functions to measure the convergence behavior of proposed scheme (1). Computed values of x n is listed with up to 15 significant digits for proper readability. The error bound = 1 2 × 10 −180 is assigned for scalar equations. Table 1. Test functions f i (x) with zeros α and initial guesses x 0 .

Local Convergence
Here log z(z ∈ C) represents a principal analytic branch with − π ≤ Im(log z) < π.
According to Table 2, sixth-order convergence is clearly seen. The values of computational asymptotic error constant agree up to 10 significant digits with η. It appears that the computational convergence order well approaches 6. In Table 3, we compare numerical errors |x n − α| of proposed Methods LK1-LK10 with those of existing Methods EM1-EM7. The least errors within the prescribed error bound are highlighted in bold face. According to the comparison, Methods LK1 and LK8 display slightly better convergence for most test functions, while other remaining methods exhibit similar convergence. A close inspection of the asymptotic error constant η(α, θ i , T f , L f ) ≈ |x n+1 −α| |x n −α| 6 addresses one's attention to the local convergence dependent on f (x), x 0 , α, T f and L f . Accordingly, the convergence of one method is hardly expected to be always better than the others.

Global Convergence
We usually locate a zero α of f (x) by means of a fixed point ξ of iterative maps: In general, W f might possess other extraneous fixed points ξ = α. It is well known that extraneous fixed points may result in attractive, indifferent or repulsive cycles as well as other periodic orbits influencing global convergence. Combining proposed methods (1) with maps (24), we find: where can be regarded as a weight function of the classical Newton's method. We are interested in the dynamics [24][25][26][27][28][29] of maps (1) underlying their extraneous fixed points of associated with their basins of attraction.
Good initial guesses for the numerical solutions of methods (1) can be determined from the basins of attraction which exhibits convergence of global character. Table 4 features statistical data including the average number of iterations per point, CPU time (in seconds), and number of points requiring 40 iterations. In the following examples, we take a 6 × 6 square centered at the origin containing all the zeros of the given test functions. We then take 601 × 601 equally spaced points in the square as initial points for methods (1). We color the point based on the root it converged to. This way we can figure out if the method converged within the maximum number of iteration and if it converged to the root closer to the initial point. Figures 1-4 present basins of attraction of 17 iterative maps when applied to various polynomials.
Example 1. As a first example, we have taken a quadratic polynomial with all real roots: Clearly the roots are ±1. Basins of attraction for all the listed methods are given in Figure 1. Consulting Table 4, we find that Methods LK1 and LK8 use the fewest iterations per point on average (AvgCon), and they also have the fewest black points. Other remaining methods have AvgCon ranging from 3.43 to 4.83. The fastest methods are EM2 with 49.406 s and EM6 with 48.032 s. Observe that Methods LK5 and LK6 exhibit more chaotic nature along the imaginary axis than others. Example 2. In our second example, we have taken a cubic polynomial: Basins of attraction are given in Figure 2. We now consult Table 4 to find that the methods with the fewest AvgCon are EM1 with 4.2698 and EM6 with 3.9447 iterations. All the others require between 3.97 and 8.44. In terms of CPU time in seconds, the fastest are EM2 (129.172 s) and EM6 (128.797 s) and the slowest are LK5 (306.672 s) and EM5 (172.548 s). The methods having the most black points are LK5 with 75,147 and LK7 with 16,293, while the methods having the fewest are LK1 with 418 points and EM6 with no points. Method LK5 displays most chaotic nature near the basin boundaries axis, followed by Method LK4. Methods LK1, LK8, LK9 and EM6 present better stability than others. Example 3. As a third example, we have taken another cubic polynomial: Now, all the roots are real. The basins for this example are plotted in Figure 3. Based on The basins are given in Figure 4. In terms of AvgCon, EM1 with 5.1189 and EM4 with 4.4454 are the best, while LK5 with 6.6932 and EM1 with 5.4899 are the worst. The fastest are EM2 with 110.578 s and EM6 with 100.078 s, while the slowest are LK4 with 541.970 s and EM1 with 285.126 s. The methods having the fewest black points are LK5 with 4 and EM2 with 0, while the methods having the most black points are LK4 with 8444 and LK1 with 5762. Methods EM1, EM3, EM5, LK3, LK4 and LK7 are more chaotic than the others, while Methods EM6 and LK6 are more stable. Method LK1 is of somewhat different stability character.  (a)EM1 Basins of attraction of the listed methods, for the roots of the polynomial z 3 + z 2 + z + 1. Figure 3. Basins of attraction of the listed methods, for the roots of the polynomial z(z 3 + 1).   Conv, convergent behavior; CPU, processing CPU time in seconds; NConpts, number of points whose each orbit is non-convergent but bounded; Conpts, number of points whose each orbit is convergent; AvgCon, average number of iterations for convergence per point; Bkpts, number of points whose each orbit tends to infinity within 40 iterations.
Theorem 1 suggests us to use T f and L f as at most third-degree matrix polynomials in (s − I).
j=4 F j e j + O(||e|| 7 ) with E j = F j (γ, c 2 , c 3 , c 4 , c 5 , c 6 , T 0 , T 1 , T 2 ). Equating the coefficients of the first and second-order terms in e z yields: We obtain: j=3 H j e j + O(||e 7 ||) with H j = H j (γ, c 2 , c 3 , c 4 , c 5 , c 6 , T 0 , T 1 , T 2 , L 0 , L 1 , L 2 ). Now, we annihilate the first five coefficients of the terms up to the fifth-order terms of e (n+1) with the use of (34) by taking the set of coefficients below: The discussions thus far lead us to the following theorem for nonlinear systems of equations.
Theorem 2. Let f : Ω ⊂ C d → C d with d ∈ N have a simple root α and be sufficiently Frétchet differentiable in Ω containing α. Let x (0) be an initial guess chosen close to α. Let T f , L f : C d×d → C d×d be matrix functions sufficiently Frétchet differentiable in a neighborhood of I, being defined by: given, then iterative scheme (1) reduces to a family of sixth-order methods satisfying the error equation below. For n = 0, 1, 2, · · · , (37) Equation (37) clearly reduces to (2) for a scalar function by identifying c i with θ i . In what follows, we employ several test examples for the zeros of vector-valued functions to verify the convergence behavior claimed here. In terms of Euclidean norm || • ||, we display the error sizes for e k = ||x (k+1) − x (k) ||, residual error || f (x (k+1) )|| and ACOC using the error criterion of ||x (k+1) − x (k) || < 10 −140 within 20 iterations.

Test Example 1
We consider a nonlinear algebraic vector equation f : R 3 → R 3 defined by f (x) = 0 with x = (x 1 , x 2 , x 3 ) T as follows: The exact solution is given by x = (1, 2, π) T . We try to solve (38) with an initial guess vector x (0) = (0.8, 1.8, 3.0) T by method (1), and find the results in Table 5. We observe that ACOC approaches up to 6, which is the theoretical order of convergence.

Test Example 2
We consider a nonlinear ODE boundary-value problem given below: The exact solution is found to be y = (sin x) 2 . With the use of the central finite difference method, the first and second derivatives are approximated by: where y n = y(x n ), h = 1 N ( π 2 − π 6 ) = π 3N , N is the number of divisions of the interval [ π 6 , π 2 ]. It can be shown that y(x + h) = y(x) + O(h 3 ), y (x) = O(h 2 ) and y (x) = O(h 2 ) in view of Taylor expansion of y(x + h) about x. This discretization yields the algebraic equations with 6 unknowns y 0 , y 1 , y 2 , y 3 , y 4 , y 5 : for j = 0, 1, · · · , N − 1, with boundary conditions y 0 = y(x 0 ) = y( π 6 ) = 1 4 , and y N = y(x N ) = y( π 2 ) = 1. Further computation after selecting N = 5 gives us a nonlinear algebraic vector equation f : R 4 → R 4 defined by f (y) = 0 with y = (y 1 , y 2 , y 3 , y 4 ) T of the form: After solving (42) with an initial guess vector y (0) = (0.6, 0.7, 0.8, 0.9) T by a typical method LK1, we find the results in Table 6 and Figure 5. It is seen that ACOC approaches up to 6, which is the theoretical order of convergence. The errors (sin x) 2 − y i for 1 ≤ i ≤ 4) at the internal nodes are, respectively, given by: Ç 0.0045808024173838738376216319522, 0.00731167760933327109201071704851, 0.0073683157995758596940525729079, 0.00474700015993753712029421111822. å As a remark, we should note that the numerical solution by the central finite-difference methods is accurate within the range of ∆y(x) = O(h 3 ) with h = 5π 3N = 0.00918704.

Test Example 3
A two-dimensional nonlinear reaction-diffusion equation for a concentration u(x, t) of the substance under consideration in a bounded domain Ω ⊂ R 2 with continuous boundary ∂Ω is represented by an initial boundary value problem: where d > 0 is a diffusion coefficient, a is a positive constant, g is continuous on ∂Ω, and ∆ is the Laplacian operator. For brevity of analysis, let d = 1, a = 1, and Ω = [0, 1] × [0, 1] (i.e., unit square region). We are interested in steady state solutions to (43), which lead us to elliptic partial differential equations with Dirichlet boundary conditions as follows: By using central divided differences with step h = 1/4 in each component of the space vector, we discretize (44) into a nonlinear system of equations with 25 nodes, 9 of which constitute interior nodal variables x 1 , x 2 , · · · , x 9 in Ω, while the remaining 16 nodes are boundary nodes. As a result, we obtain a nonlinear algebraic vector equation f : R 9 → R 9 defined by: Ψ(x) = (x 2 1 , x 2 2 , · · · , x 2 9 ) T and b = ( 29 16 , 7 8 , 29 16 , 7 8 , 0, 7 8 , 29 16 , 7 8 , 29 16 ) T . We solve (38) with an initial guess vector x (0) = (1, 1, 1, 1, 1, 1, 1, 1, 1) T by a typical method LK1, and find the results in Table 7. It is evident that ACOC reaches up to 6, being the theoretical order of convergence. As can be seen in Table 7, the methods with γ = 2/3 appear to converge more quickly and better than those with γ = 1.
Interior 16 nodal values of the steady-state solution of u(x, t) are illustrated with adjacent nodal points connected by straight lines in Figure 6.
The above nonlinear system is described in [4]. Selecting n = 10, we find d = 10 and solve (46) in R 10 with an initial guess vector x (0) = (0.75, 0.75, · · · , 0.75) T for the desired root x = (0.5149332, 0.5149332, · · · , 0.5149332) T ∈ R 10 in Table 8. It is evident that ACOC reaches up to 6, being the theoretical order of convergence. As can be seen in Table 8, the methods with γ = 2/3 appear to converge more quickly and better than those with γ = 1.

Computational Efficiency
The computational efficiency of an iterative method is defined by an efficiency index E = ρ 1/d [30], with ρ as the order of convergence and d as the number of functional evaluations per iteration. We require n scalar functions for each f and n 2 for each f . The concept of the efficiency index E applied to a nonlinear system of vector equations has been extended to treat the concept of computational efficiency by using CE = ρ 1/(d+op) [4], where op is the number of operations associated with products and quotients. Suppose that n is the size of the matrix needed in the nonlinear system of vector equations. Matrix inversion requires n 3 −n 3 product-quotient operations and LU-decomposition technique for solving linear systems requires n 2 product-quotient operations, including the n 2 product-quotient operations related to matrix multiplication by a vector. Note that each method treated here follows the three set of linear systems (1) and has one matrix inverse f (x n ) −1 . Consequently, the number of functional evaluations plus product-quotient operations d + op becomes 2n + 2n 2 + n 3 −n 3 + 6n 2 = n 3 +18n 2 +5n 3 , which gives us the computational efficiency CE = 6 3 n 3 +18n 2 +5n for each listed method.
Many real-life application problems include ones related to: interval arithmetic benchmark, neurophysiology, chemical equilibrium, kinematic application, combustion application, and economics modeling, whose studies are described in [31]. The methods used therein are based on second-order Newton-like approach which may be more efficient in real-life problems in terms of speed and computational cost. On the other hand, our proposed family of sixth-order methods (1) is much more accurate than Newton-like methods, but has more complexities owing to the high-order formulation and require more CPU time to get the desired solution.
One certainly has to acknowledge that determining a better method than the other one should be avoided through solving a function with a randomly chosen initial guess vector and comparing the number of convergent iterations.

Conclusions
A family of Jarratt-like iterative methods for scalar and vector equations is developed and its convergence properties are theoretically established through Theorems 1 and 2. Computational aspects applied to various test equations agree well with the convergence behavior claimed in the theory developed. Global convergence behavior of the listed methods is illustrated for typical polynomials based upon their basins of attraction. The basins of attraction suggest selecting members of the iterative methods (1) give better convergence.
We will focus our future study on extending the current approach with different weight functions to the development of higher-order iterative root-finders.