A New Nonlinear Ninth-Order Root-Finding Method with Error Analysis and Basins of Attraction

: Nonlinear phenomena occur in various ﬁelds of science, business, and engineering. Research in the area of computational science is constantly growing, with the development of new numerical schemes or with the modiﬁcation of existing ones. However, such numerical schemes, objectively need to be computationally inexpensive with a higher order of convergence. Taking into account these demanding features, this article attempted to develop a new three-step numerical scheme to solve nonlinear scalar and vector equations. The scheme was shown to have ninth order convergence and requires six function evaluations per iteration. The efﬁciency index is approximately 1.4422, which is higher than the Newton’s scheme and several other known optimal schemes. Its dependence on the initial estimates was studied by using real multidimensional dynamical schemes, showing its stable behavior when tested upon some nonlinear models. Based on absolute errors, the number of iterations, the number of function evaluations, preassigned tolerance, convergence speed, and CPU time (sec), comparisons with well-known optimal schemes available in the literature showed a better performance of the proposed scheme. Practical models under consideration include open-channel ﬂow in civil engineering, Planck’s radiation law in physics, the van der Waals equation in chemistry, and the steady-state of the Lorenz system in meteorology.


Introduction
Applied mathematics, biological sciences, and engineering disciplines are full of systems and processes whose response is not proportional to the input, thereby demanding numerical schemes with a high order of convergence and a reasonable computational cost. The literature has provided us with various numerical schemes to fulfill this demand, while still leaving significant room for continuous improvements in this vibrant area of numerical analysis. Fortunately, the need to solve nonlinear equations does not end with mathematical studies since many ordinary, partial, integral, and integro-differential equations often give rise to systems of nonlinear equations after discretization, so they must be handled with numerical schemes. Generally, nonlinear phenomena designed in terms of models with nonlinear equations are seen not only in academia and the physical industry, but also in the fields of medical sciences and other real-world situations. For instance, the triangulation of GPS signals, combustion, heat transfer, and fluid flow amongst others are some of those fields [1][2][3][4].

Previous Materials and Methods
In this section, we will briefly review some well-known numerical schemes that are frequently used to determine approximate solutions to nonlinear models. Indeed, when it comes to solving such nonlinear models, there is objectively no analytical technique to yield exact results. Thus, we are obliged to use numerical schemes.

Brief Overview of Some Existing Schemes
When someone thinks of solving a nonlinear equation with a numerical scheme, the first scheme that comes to mind is the classical Newton-Raphson [5] which has optimal second-order convergence. In short, we used the symbol NR2 for this scheme which uses two function evaluations-one function and one first-order derivative-per iteration. The scheme is given by In [20], one can find the classical third-order Halley's numerical scheme (denoted by Halley3). This scheme needs three function evaluations (one function, one first-order derivative and one second-order derivative) per iteration and is given by By using NR2 (1) as a predictor, Halley's as corrector, and the finite difference approximation of the second derivative, a new numerical scheme was suggested in [21] named the modified Halley method (MHM5). It possesses fifth-order convergence and needs four function evaluations (two function and two first-order derivatives) per iteration. The scheme is given by With the use of five function evaluations (three function, and two first-order derivatives), a three-step numerical scheme was discussed in [22] whose theoretical order of convergence, as proved by the authors, is 9. The scheme is given by A two-step numerical scheme of order nine proposed in [23] uses seven function evaluations per iteration (two function, two first-order derivatives, two second-order derivatives and one third-order derivative). The scheme is given by Most recently, the authors in [24] devised, using a variational iteration approach, a new three-step numerical scheme to solve nonlinear scalar equations. Their proposed scheme is ninth-order convergent with six function evaluations (three function, two first-order derivatives and one second-order derivative) per iteration. The scheme is as follows: The numerical schemes discussed above will be used for comparison purposes in the present study.

Proposed Numerical Scheme
After the thorough investigation of the existing literature, we observed that several authors have proposed new hybrid type schemes using the second-order NR scheme as their starting step, aiming to increase the order of convergence and reduce the number of function evaluations. However, we took a step forward to use the third-order Halley's scheme in the first step: to be blended with the third-order modified Newton scheme given below: The aim was to obtain a new nonlinear scheme with a higher order of convergence while the number of function evaluations per iteration was reduced as much as possible. Such an idea was carried out in [13] and later used by many authors. Hence, after blending both schemes above (7) and (8), we obtain the following new nonlinear scheme: The proposed scheme will be referred to as PS9 for the rest of the paper, wherein its theoretical analysis, including convergence theory, stability, and bifurcation phenomenon are discussed in detail.

Convergence Analysis for
In the following theorem, we prove the ninth-order of convergence for the proposed nonlinear three-step scheme given in (9). Theorem 1. Let γ be the root of a sufficiently differentiable function f : Π ⊂ R → R on an open interval Π. Then, the three-step scheme PS9 presented in (9) has ninth-order convergence, and the error equation is: where A = 1 128 Proof. Let γ be the root of f (x) = 0, x j be the jth approximation to the root provided by PS9 and ε j = x j − γ be the error term after the jth iteration. Using Taylor's expansion for f (x j ) about γ, we have: Expanding f (x j ) via Taylor's series about γ, we have: Multiplying (11) and (12), we obtain: . (13) Now, the Taylor expansion of the denominator in the fractional term in the first step of (9) multiplied by its numerator results in: Using (14) in the first step of (9), we obtain: Using Taylor's expansion for f (y j ) about γ, we have: Using the Taylor's expansion for 1/ f (y j ) about γ, we have Multiplying (16) by (17), we obtain: Using (18) in the second step of (9), we have: Now, using the Taylor expansion for f (z j ) about γ, we have: Using (16), (17), and (20) in the third step of (9), we have (note that, for brevity, we suppressed the γ in the derivatives): Using (19) in (21), we have: Finally, using (15) in (22), we obtain: which proves the ninth-order convergence of the proposed scheme PS9 in (9).
2.2.2. Convergence Analysis for F(x) = 0 Now, let us consider a nonlinear system of N equations and N variables, F(x) = 0, where F : Π ⊂ R N → R N is a sufficiently Frechet differentiable function. In this case, the PS9 method adopts the form: where F is the Jacobian matrix of F.
We present the following results to show the error equation, and therefore, the order of convergence of the proposed scheme in the case of a system of nonlinear equations.
Then, for any x and h ∈ R N the following expression holds where ||R r || ≤ 1 r! sup 0<t<1 ||F (r) (x + ht)||||h|| r and the notation Theorem 2. Let the function F : Π ⊂ R N → R N be sufficiently differentiable in a convex set Π containing a simple zero γ of F(x). Let us consider that F (x) is continuous and nonsingular in γ. If the initial guess x 0 is close to γ then the sequence x j obtained with the proposed three-step scheme PS9 converges to γ with ninth-order convergence.
Proof. Let γ be the root of F(x), x j be the jth approximation to the root by PS9 and ε j = x j − γ be the error vector after the jth iteration. Using Taylor's expansion of F(x j ) about γ, we have: Using Taylor's expansion of F (x j ) about γ, we have: Multiplying (26) and (27), we obtain: Using Taylor's expansion for the inverse matrix 2F ( we obtain: Applying (29) to 2F(x j )F (x j ) after using (28), we have: Using (30) in the first step of (24), we obtain: Likewise, using Taylor's expansion for F(y n ) and F (x n ) about γ; and F (y n ) in the second step of (24), we obtain: Finally, we use the Taylor series for F(z n ) about γ and earlier found series to obtain the following: Putting (31) in the above equation, we obtain the required result as follows: Hence, the ninth-order convergence of the proposed scheme PS9 has also been proven in the case of nonlinear systems.

Numerical Comparison
A few numerical schemes were briefly introduced in Section 2.1, whose practical comparison can be made with the proposed scheme. The following abbreviations will be used hereinafter:

NR2
Newton-Raphson scheme with second-order convergence given in (1); Halley3 Halley's scheme with third-order convergence given in (2); MHM5 Modified Halley's scheme with fifth-order convergence given in (3); HD9 A numerical scheme with ninth-order convergence given in (4); NZ9 A numerical scheme with ninth-order convergence given in (5); NA9 A numerical scheme with ninth-order convergence given in (6); PS9 Proposed numerical scheme with ninth-order convergence given in (9); FE Number of function evaluations; EI Efficiency index; γ Exact root of f (x) = 0; IG Initial guess x 0 ; COC Computational order of convergence.
Based upon orders of convergence, the number of function evaluations per iteration in both the scalar and vector versions of the given nonlinear problem, as well as efficiency indices, we present in Table 1 a practical comparison of the schemes under consideration. It can be observed that the number of function evaluations per iteration is either equivalent or smaller than those having ninth-order convergence, whereas the efficiency index of PS9 possesses a similar type of character as HD9, being the only scheme with the highest efficiency index. The last column of Table 1 shows the number of new function evaluations with the nonlinear problem being a system with a dimension n > 1. In such a case, the coefficient of n shows the number of function evaluations, the coefficient of n 2 shows the number of first-order derivatives, the coefficient of n 3 shows the number of second-order derivatives, and finally the coefficient of n 4 is for the number of third-order derivatives; one per each iteration step. In this situation, PS9 stands with Halley3 and NA9 in terms of efficiency index, whereas NZ9 takes an extra function evaluation with the fourth-order derivative.

Basins of Attraction
One of the dynamical concepts concerning the stability and convergence of roots for f (z) = 0 by some iterative methods is the concept of a basin of attraction which is a region of 2D space within which iterations will be iterated into the attractor irrespective of the choice of an initial guess z 0 . This concept is found to have several applications in scientific fields including applied mathematics, arts, design, architecture, and many more. Interesting applications of basins of attraction can be seen in the references [25][26][27][28][29] and most of the references cited therein. Nonlinear models under both scalar and vector cases have aesthetically interesting behavior that is not usually observed in the linear models. The basins of attraction will demonstrate various dynamical characteristics under the proposed iterative method PS9. Such regions, in the present study, were obtained through graphical utilities provided by some routines of MATLAB R2020a. To serve the purpose, we selected a rectangular region R on [−2, 2] × [−2, 2] under discretization of R into 1000 by 1000 grid points containing the solution of the complex-valued functions under consideration with considerably higher resolution. Some complex valued functions, as given below, are taken to illustrate the dynamical aspects of the behavior of PS9 via aesthetically amazing basins of attraction and their corresponding iteration coloring: The above functions were selected, randomly containing both complex polynomials and transcendental functions. Graphs were obtained with MATLAB R2020a within R = [−2, 2] × [−2, 2] taking = 10 −2 and k = 60, where k refers to the maximum number of iterations chosen. It should be noted that the region R is the phase-space wherein the basins of attractions are plotted under PS9 for the complex-valued functions described in Example (1). As seen in Figures 1-6, the highest number of iterations is required by those initial guesses that are at or near the boundaries (black lines in the figures). When the initial guess z 0 is taken within a small neighborhood of the exact solution, then the minimum number of iterations is needed. While running the simulations to obtain the fractal structures shown in Figures 1-6, we computed the total average time (seconds) required by PS9 for the complex functions under consideration, as given in Table 2. As observed, complex functions with higher degrees (in the case of polynomials) and hyperbolic functions are computationally more expensive. Table 2. Execution time (seconds) required by PS9 for each P i (z), i = 1, 2, . . . , 6.

Numerical Experiments
In this section, we performed some numerical simulations to illustrate the performance of the proposed scheme with ninth-order convergence, i.e., PS9, as given in Equation (9). Different types of nonlinear functions of single and multi-variable character are taken into consideration. To further strengthen the results, some practical nonlinear models were also touched upon, including those used in meteorology, civil, environmental, and chemical engineering. We use the stopping criterion | f (x n )| ≤ ε, where the tolerance is chosen as ε = 10 −200 and the required precision is set to be as large as 16,000. The numerical results obtained have been tabulated, where each table contains initial guesses, the number of iterations (N), the number of function evaluations (FEs), the computational order of convergence (COC), the absolute error (|δ|)), the absolute value of the underlying function at the final approximated solution (| f (x N )|), and the CPU time consumed by each scheme in seconds. The COC has been estimated using the usual formula [30,31]: The numerical results were performed with the Mathematica 12.1 system running under Windows 10 OS with an installed memory of 24.0 GB having a processor Intel(R) Core(TM) i7-1065G7 CPU @ 1.30 GHz speed.
The PS9 performance shows the smallest absolute error and therefore the smallest functional value along with comparable CPU times, as can be seen in Tables 3-7 for the five test functions listed in problem 2. Under different initial guesses taken near and far from the approximate roots, it can be seen from these tables that the proposed strategy takes reasonable and sometimes smaller iterations to meet the stopping criteria. The COC turns out to be 9, confirming the theoretical ninth-order of convergence derived in Section 2.2.1. The functional absolute values at the last approximation considered and the corresponding absolute errors are smaller with PS9 compared to other schemes, while consuming a reasonable amount of CPU time. Moreover, schemes given as NR2, NZ9 and NA9 in [5,23,24], respectively, could not produce results for the test functions f 4 (x) under the initial guess x 0 = −5, and f 5 (x) under the initial guess x 0 = 0.8 including the scheme MHM5.
where m 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 W and the depth of water in the channel h, it is known that A = Wh, and: With these values, (37) becomes: Table 6. Numerical results for problem 2 ( f 4 (x)) with two different initial guesses. Here, "div" stands for "divergence". Table 7. Numerical results for problem 2 ( f 5 (x)) with two different initial guesses. As per the requirement wherein one needs to determine the depth of water h in the channel for a given quantity of water, the above equation can be rewritten as follows:

Scheme IG
The depth of water h in the channel has been estimated while assuming the rest of the parameters as Q = 14.15 m 3 /s, W = 4.572 m, n = 0.017 and m = 0.0015. The initial guess is taken to be h 0 = 8.5 m. The numerical results under different schemes are shown in Table 8. It can be observed from this table that the functional value and the absolute error are smallest with HD9, whereas PS9 competes well enough and consumes less time while maintaining ninth-order convergence, as given by COC in the fourth column of Table 8. Consider the following Planck's radiation law problem: which calculates the energy density within an isothermal blackbody. Here, λ stands for the wavelength of the radiation, T is the absolute temperature of the blackbody, k is Boltzmann's constant, h is the Planck's constant, and c shows the speed of light. Here, we compute the wavelength λ which corresponds to the maximum energy density φ(λ). Using (41), we obtain: or: It is easy to see that a maximum for φ occurs when M = 0, that is, when: Substituting x = ch/λkT, the above equation can be rewritten as The numerical results from the application problem 4 regarding Planck's radiations are listed in Table 9 in which it is observed that the functional value and the corresponding absolute error are smaller under PS9 with a reasonable amount of CPU time while taking only five iterations to meet the stopping criterion. It is worth noting that the schemes NR2, Halley3, MHM5, and HD9 converge to a solution other than the required one. Table 9. Numerical results for the problem 4 with initial guess x 0 = 1.5 m. A well-known model from chemistry used for determining the difference between ideal and real gases is known as the van der Waals nonlinear equation which is given as follows: where k and h are the positive real parameters also known as the van der Waals constants which depend on the type of gas under consideration. Moreover, constants P, V, T are the pressure, volume, and temperature of the gas. R ≈ 0.0820578 is the gas constant and n represents the number of moles. The equation given above is simplified as follows: The following third degree polynomial can be obtained under k = 18, h = 0.1154, n = 1.4, P = 40 atm, and T = 500 • C: The approximate solution of (48) correct to 40 dp is as follows: 2.1544346900318837217592935665193504952593.
We solved the Equation (48) with the schemes under consideration taking V 0 = 0.0, resulting in the data shown in Table 10. It can be observed that PS9 produces the smallest absolute error with the least CPU time, followed by the ninth-order scheme NA9. Each scheme, for this particular problem, needs many more iterations than PS9. This applied problem shows the highest performance of PS9 in terms of number of iterations, number of function evaluations, absolute error, absolute functional value and CPU time. The schemes NR2 and PS9 have been used for solving two nonlinear three-dimensional systems given by problems 6 and 7. Here, we chose only NR2 and PS9 because the simula-tions obtained under other schemes behave dynamically in the same way as that shown in the univariate cases discussed above. For the steady-state of the Lorenz Equations (49) and the system (51), Tables 11 and 12 show that, at little high computational cost, PS9 yields the smallest error norm δ n+1 = ||x n+1 − x n || ∞ and the smallest absolute functional values under different initial conditions. Additionally, PS9 takes very few iterations to meet the stopping criterion. The numerical results of these two systems are in favor of PS9, showing that PS9 is applicable for both univariate (single nonlinear equation) and multivariate (system of nonlinear equations) cases, while maintaining its ninth-order convergence. Example 6. Lorenz system in steady-state [32] (p. 816) The approximate solution for the system (49) correct to 50 dp is given by where its approximate solution up to 50 dp is: The numerical results for the system (51) are shown in Table 12 showing the good performance of the proposed method.

Concluding Remarks
A new ninth-order convergent scheme, with six function evaluations per iteration, to solve nonlinear models in both univariate and multivariate cases is proposed by blending a third-order modified Newton's and the third-order Halley's schemes. The theoretical order of convergence with the asymptotic error term was obtained via Taylor's expansion and the efficiency index was calculated to be approximately 1.4422. The proposed scheme was found to be accurate and reliable enough in terms of iteration counts, number of function evaluations, absolute errors, and CPU time. This was further confirmed when the schemes under consideration were used to solve several applied models in different fields. Two-dimensional phase diagrams, using the proposed scheme, were obtained to show basins of attraction for the stability of PS9.
The presented scheme (9), which comes in the category of the schemes without memory, requires the evaluation of three inverse operators per iteration. In addition, such schemes also require many multiplication and division operations while computing derivatives and inverse operators. All such evaluations will increase, to some extent, the amount of computational work of the scheme. Thus, future work would include the modification of PS9 as a numerical scheme with memory, keeping in mind the time efficiency and the computational complexity. Moreover, first-and second-order derivatives involved in PS9 can be replaced using appropriate finite-difference approximations to avoid complexity with the analytical expressions. Such analyses will be carried out in our future research studies.