Development of a New Multi-step Iteration Scheme for Solving Non-Linear Models with Complex Polynomiography

The appearance of nonlinear equations in science, engineering, economics, and medicine cannot be denied. Solving such equations requires numerical methods having higher-order convergence with cost-eﬀectiveness, for the equations do not have exact solutions. In the pursuit of eﬃcient numerical methods, an attempt is made to devise a modiﬁed strategy for approximating the solution of nonlinear models in either scalar or vector versions. Two numerical methods of second-and sixth-order convergence are carefully merged to obtain a hybrid multi-step numerical method with twelfth-order convergence while using seven function evaluations per iteration, resulting in the eﬃciency index of about 1.4262. The convergence is also ascertained theoretically, and the asymptotic error constant is computed. Furthermore, various numerical methods of varying orders are used to compare the numerical results obtained with the proposed hybrid method in approximate solution, number of iterations, absolute error, absolute functional value, and the machine time measured in seconds. The obtained results outperformed the chosen methods when applied models arising from physical and natural ﬁelds were solved. Finally, to observe the convergence graphically, some complex polynomials are plotted as polynomiographs, wherein the rapid convergence of the proposed method is guaranteed.


Introduction
Computing the approximate zeros of the nonlinear scalar and vector functions is one of the most important and interesting research areas in the modern age.
ere are many applications of the root-finding methods in different disciplines of science as well as in arts and economics. With the help of several mathematical techniques, a variety of complex problems in different applied sciences can be modulated in the form of nonlinear equations and then can be readily solved via different root-finding techniques. A root-finding method in mathematics and computer technology is a method for finding zeros, commonly known as "roots" of continuous functions. A zero of a real-valued or a complex-valued function f, is a value r such that f(r) � 0. Mostly, the root-finding techniques give approximations to zeros, expressed either as floating-point integers or as tiny isolating intervals, or discs for real or complex roots, because the zeros of a function cannot be calculated accurately with available analytical techniques.
Most numerical approaches for root-finding rely on the iteration process, which generates a series of discrete points that should converge towards the root as a limit. ese iteration schemes start with one or more estimations of the root as initial inputs, and every iteration of the process generates a more accurate estimation of the exact root [1][2][3][4].
Since iterations must be ended at some point, these approaches yield estimation to the root rather than a precise solution.
Many approaches compute successive approximations by considering an auxiliary function on the values that came before them. erefore, the limit is a fixed point of the function f, which is selected to have the solution of the original equation as fixed points and to converge to these fixed points quickly. In pursuing accurate and efficient roof-finding techniques, several techniques have been proposed in the past and current literature from the Newton method through the techniques proposed in the ongoing year. Perhaps, the most commonly used root-finding method is the Newton Rahpson method N2 [5] with quadratic convergence. Its computational step is shown below that uses two function evaluations: 1 for the function f(x) itself and 1 for the first-order derivative f ′ (x): where f ′ (x i ) ≠ 0. e researchers in [6] devised a modified version of an existing algorithm aiming at the removal of first-order derivatives.
ey came up with a two-step method with fourth-order convergence abbreviated by N4. One of the advantages of the algorithm was the use of only four function evaluations per iteration, as depicted in the following computational scheme: where (3) e researchers in [7] devised a three-step method having eighth-order convergence as denoted by W 8 . Despite being eighth-order convergent, one of the advantages of the algorithm was the use of only four function evaluations per iteration, as depicted in the following computational scheme:

Formulation of the Proposed Method
It has been observed in the current literature that new modified root-finding techniques are being proposed because of increasing the efficiency of the existing ones. In search of such algorithms, some researchers have merged two existing iterative methods of convergence order m and n to obtain a hybrid method with convergence order mn. In this respect, the convergence order is improved. Nonetheless, the computational aspect was ignored, resulting in an increased number of function evaluations in most newly modified hybrid approaches. For example, authors in [13] proposed an iterative third-order method with five function evaluations required per iteration, including another algorithm presented in the same research work with a fourthorder three-step method that requires eight function evaluations. Likewise, authors in [14,15] have employed an excessive number of first-order derivatives, leading to high computational effort and machine time. e primary concern of the present work is to propose a hybrid method with possible higher-order convergence with the minimum number of function evaluations so that the computational cost in terms of arithmetic operations and CPU time be reduced. e proposed method comes from Newton's method (m � 2) and a three-step method (n � 6) in [16,17], leading to produce the proposed approach with convergence order mn � 12 while using just seven function evaluations per iteration. e four-step proposed method results in the following computational steps, whose flowchart is depicted in Figure 1: for i � 0, 1, 2, . . .. e methods as mentioned earlier, including the one proposed as the four-step hybrid method in the present article, are compared in Figure 2 with each other based on efficiency index (p 1/n ), order of convergence (p), and the number of function evaluations (n) used by each method per iteration.

Convergence Analysis of the Proposed Method
is section has been divided into two sections wherein the convergence of the proposed four-step hybrid method in both scalar and vector form has been discussed in detail. It is noted that the twelfth-order convergence is theoretically verified in each case with the aid of Taylor's expansion.

Convergence eory with Scalar Form.
In this subsection, we theoretically prove the local order of convergence for the proposed method given in (9). Theorem 1. Suppose that α ∈ Q be the required simple root for a sufficiently differentiable function f: Q⊆R ⟶ R within an open real interval Q. en, the proposed four-step numerical method (9) possesses twelfth-order convergence with the error equation given by: where

Proof.
Suppose α be a simple root of f(x) i � 0, where x i be the i th approximation for the root by the proposed method (9), and ζ i � x i − α be the error term in variable x at the i th iteration step. Employing the single real variable Taylor's series in [9] for f(x i ) around α, we obtain Similarly, using the Taylor's series for 1/f ′ (x i ) around α, we obtain Multiplying (11) and (12), we obtain Figure 1: Flow chart for the proposed four-step hybrid method given in (9).

N2
N4 W8 N9 P10 N11 N12 P12   Complexity Substituting (13) in the first step of (9), we obtain where η i � w i − α. Using the Taylor's series for f(w i ) around α, we obtain Similarly, using the Taylor's series for Multiplying (15) and (16), we obtain Substituting (17) in the second step of (9), we obtain where ] i � y i − α. Using the Taylor's series for f(y i ) around α, we obtain Similarly, using the Taylor's series for 1/f ′ (y i ) around α, we obtain Multiplying (19) and (20), we obtain Substituting (21) in the first step of (9), we obtain where Substituting (19), (20), and (23) in the fourth step of (9), one obtains Finally, using , and (14) and (18), (22) for the above equation, we obtain Hence, the twelfth-order convergence of the proposed method P12 given by (9) for the nonlinear functions in single variable (f(x) � 0) is proved.
Proof. Suppose α be a simple root of f(x i ) � 0, where x i be the i th approximation for the root by the proposed method (9), and ζ i � x i − α be the error term in variable x at the i th iteration step. Employing the multi variable Taylor's series given in the theorem [9] for f(x i ) around α, we obtain Again, using the Taylor's expansion for the inverted Jacobian matrix f ′ (x i ) − 1 around α, we obtain 6 Complexity Multiplying (27) and (28), we obtain Substituting (29) in the first step of (9), we obtain Using the Taylor's series for f(w i ) around α, we obtain Again, using the Taylor's expansion for the inverted Jacobian Multiplying (31) and (32), we obtain Substituting (33) in the second step of (9), we obtain where ] i � y i − α. Using the Taylor's expansion for f(y i ) around α, we obtain Again, using the Taylor's expansion for the inverted Jacobian matrix f ′ (y i ) − 1 around α, we obtain Multiplying (35) and (36), we obtain Substituting (37) in the first step of (9), we obtain Complexity where ρ i � z i − α. Using the Taylor's series for f(z i ) around α, we obtain Substituting (35), (36), and (39) in the fourth step of (9), one obtains Finally, using , and (30) and (34), (38) for the above equation, we obtain Hence, the twelfth-order convergence of the proposed multi-step (four-step) hybrid method P12 mentioned by (9) for the nonlinear functions in multi-variable (f(x) � 0) case is proved. □ □

Polynomiography
Polynomiography is a process that integrates mathematics and art to create a new type of visual art. e produced graphics result from algorithmic visualization of iterative approaches for solving a polynomial equation. is term was first introduced by Dr. Bahman Kalantari at the start of the 21st century [18]. Dr. Bahman Kalantari's study on polynomial root-finding, which is an old and traditional discipline that continues to find new implications with each generation of mathematicians and scientists, inspired the concepts of Polynomiography. Dr. Kalantari invented the term "polynomiography," which is a mixture of the word "polynomial" with the suffix "graphy." A "polynomiograph" is a separately produced picture resulting from Polynomiography. It is defined as "An iterative procedure for producing two-dimensional colored pictures (polynomiographs) that represent polynomials." In recent years, researchers worked in the field of Polynomiography along with its implementations in other fields. In [19], the authors introduced a new mathematical art with the help of Newton-Ellipsoid method. Gdawiec et al. in [20], presented the visual analysis of Newton's method with fractional-order derivatives. e authors employed the techniques of coloring by roots and coloring via iterations to study the convergence and dynamical aspects of the processes visualized by polynomiographs.
Naseem et al. in [21] presented some new graphical art with the help of newly suggested ninth-order iteration schemes. Scot et al. in [22], presented the basin of attraction for several methods and examined its dependence on their convergence orders. In [23], the authors introduced a new family of eighth-order methods and then drew their basins of attractions by assuming different polynomials. Finally, in [24], the authors generated some new fractal patterns by combing two root-finding methods. e obtained fractal patterns were diverse and had many applications in the textile and ceramic industries.
We use a rectangle R ∈ C along with the dimension [− 2, 2] × [− 2, 2], accuracy ε � 1 × 10 − 3 and the max. no. of iterations T � 20 to create the polynomiographs over the complex plane C through the computer software by taking multiple complex polynomials. e color black is allocated to the spots where the method failed to converge. e partitioning of R determines the pixel density of the created visual representations; for example, if we partition the rectangle R into a grid of 2000 × 2000, the plotting polynomiographs will then have better resolution.
For drawing graphical objects in the complex plane, we use the four complex polynomials listed below: For coloring the iterations, we employ the colormap given in Figure 3.

Problem 1. Polynomiographs for the Polynomial q 1 (t)
rough Various Methods In this example, we investigate and compare the dynamical results obtained through different iteration schemes with our presented method by considering the cubic 8 Complexity polynomial t 3 − 1 which possesses three distinct simple zeros: We executed all the methods to achieve the simple zeros of the considered polynomials and the results can be visualized in Figure 4.

Problem 2. Polynomiographs for the Polynomial q 2 (t)
rough Various Methods In this example, we investigate and compare the dynamical results obtained through different iteration schemes with our presented method by considering the 3rd-degree polynomial 3t 3 + 2t 2 − t + 1 which possesses three distinct simple zeros: We executed all methods to achieve the simple zeros of the considered polynomials and the results can be seen in Figure 5.

Problem 3. Polynomiographs for the Polynomial q 3 (t)
rough Various Methods In this example, we investigate and compare the dynamical results obtained through different methods with our presented method by considering the 4th-degree polynomial t 4 − 1 which possesses the following simple zeros:1, − 1, i, and − i. We executed all the methods to achieve the simple zeros of the considered polynomials and the results are given in Figure 6. rough Various Methods In this example, we investigate and compare the dynamical results obtained through different iteration methods with our presented method by considering 4th-degree polynomial (t 2 − 1)(t 2 − 2) which possesses the simple zeros:1, − 1, 2, and − 2. We executed all methods to achieve the simple zeros of the considered polynomials and the results can be visualized in Figure 7.
In all given examples, a detailed graphical analysis of the designed algorithm has been provided via polynomigraphs. For plotting polynomiographs on the complex plane, we take two cubic-degree polynomials namely: q 1 (t), and q 2 (t) and two quatric-degree polynomials represented by q 3 (t), and q 4 (t) respectively. e plotted graphs tell us about the convergence speed and the iterations performed by the method for drawing these objects. e second characteristic is the dynamics of the iteration scheme. In each polynomiographs, the individual root has been denoted by the blue colored dot. e black colored zones denote the divergence area or deficiency of the method through which the polynomiographs has been plotted. e darker or brighter zones in the provided polynomiographs showing less iterations performed to approximate the solution. One can easily observe the superiority of the proposed method over the others by examining the more darker or brighter zones of the polynomiographs drawn by the suggested method.

Numerical Simulations: Realworld Scenarios
is section of the paper discusses the real-life applications by applying the newly proposed hybrid method. We also present a numerical comparison with other existing most frequently used methods: N2, N4, W8, N9, P10, N11, and N12, whose computational steps are shown in the introductory section above. In each applied model, we set the tolerance to be ε � 10 − 100 as the stopping criterion of the iterative process of every method under consideration: |x N − x e | < ε. Two additional methods with fifth-and sixthorder of convergence are also included for the simulations of a six-dimensional model chosen from the field of neurophysiology based on the reason that some of the methods under discussion in the above introductory section did not prove to be valuable candidates when it comes to the simulations of the nonlinear models presented in the system or vector form.

Problem 5.
e Plank's radiation law in physics explains the spectral density of radiation emitted by a black body in thermal equilibrium at the temperature T and the condition that there must not be a flow of energy between the body and its surroundings. In other words, the law is introduced to determine the amount of energy density in a black body based on isothermal properties. Moreover, it is sometimes used to estimate the maximal radiation's wavelength. As described in [25], the maximal wavelength of the radiation may be written in the form of the following nonlinear equation in scalar version: where x stands for the maximal wavelength. e exact solution of the above equation is as follows: 0.0. e Plank's radiation model described in (44) is simulated with several numerical algorithms while assuming two different initial guesses. e maximum number of iterations in each case is set to be N � 4. It can be observed in Table 1 that the accuracy is maximum for P12 at the cost of a slightly higher amount of CPU time, regardless of the initial estimate's location. e eleventh-order method N11 did not converge with the initial estimate taken to be x 0 � 0.49 while it converges second to P12 when the initial guess x 0 � 0.75 is taken. With this second initial guess, the CPU time consumption also slightly increases for N4 and N11.
Problem 6. Fraction Conversion of Nitrogen-Hydrogen to Ammonia [26]. is nonlinear scalar problem depicts the fractional conversion of nitrogen-hydrogen to ammonia and has appeared in several research works conducted in the past and recent literature. In this experiment, we set the pressure value to be 250 atmospheric pressure while the temperature is set to 500 degrees Celsius. In terms of a nonlinear function, the model mentioned above has the following polynomial form: It has been identified in the recent work [27] that one of its positive real roots lies in the open unit interval (0, 1) which is estimated to be 0.2777595428. For this nonlinear model, the numerical simulations are shown in Table 2 while the number of iterations for each method under consideration is set to be 7. Two different initial guesses, that is, x 0 � 0.5 (near to the exact solution) and x 0 � 0.95 (away from the exact solution), are chosen. It is seen that the fourth-order method N4 converges towards some other solution for the first initial guess while the method abbreviated as N11 failed after three iterations while the same is the case for the W8 method, but the method W8 produced the correct approximate solution till four iterations and failed after that.
e Newton method N2 has the comparatively most significant absolute error at the fourth iteration compared to other methods. Nonetheless, the proposed hybrid method, in addition to a few more methods, always converged towards the required solution. More so, the proposed method has achieved the minor error tolerance with a reasonable amount of time. Hence, it can be concluded that the initial location of the estimate does not matter much when it comes to the proposed hybrid method devised in this article. Figure 4: Polynomiographs for q 1 (t) using different methods.
Problem 7. Application of mechanical engineering [28]: ere are several fields wherein the use of thermodynamics is extensively required. e particular areas include mechanical, civil, mechatronics, electronic, chemical engineering, and many others. e fourth-order polynomial is used to show a relation between the zeropressure specific heat of dry air c p (kJ/kgK) to temperature x: where c p � 1.2 is used. As described in [28], the above nonlinear model given in terms of fourth-order polynomial has two real distinct roots given as: r 1 � 1126.009751 and r 2 � − 1289.950382. It can be observed in Table 3 that each method converges for the initial guesses chosen to determine the approximate solution of the above model. e eleventh-order method N11 performed better from an accuracy viewpoint, followed by the proposed hybrid method. Looking at the CPU values, it is clear that the methods N11 and P12 take the same amount of Figure 5: Polynomiographs for q 2 (t) using different methods.
Complexity time, thereby being equally time-efficient for this particular fourth-order polynomial. Moreover, the absolute functional values are identical for both methods, including others such as N9 and P10.
Problem 8. Neurophysiology Application: As a final experiment, we consider a six-dimensional nonlinear system first proposed in [29] and later was used by several researchers for the simulation purpose of their newly developed algorithms. See, for example, [30], and some cited references therein. e nonlinear model consists of the following six equations: e constants c i in the above model can be randomly chosen. In our experiment, we considered c i � 0, i � 1, . . . , 4. Figure 6: Polynomiograph for q 3 (t) using different methods. Figure 7: Polynomiograph for q 4 (t) using different methods.  e simulations for the above neurophysiology application system are shown in Table 4 wherein the two columns represent the absolute error at the last iteration and the CPU time consumption. Two more methods, including fifth-order Halley's (HM5) in [6] and a sixth-order Hameer-Mujtaba method (HM6) in [17], is used for the simulations of the above system. It is evident from Table 4 that the accuracy is much higher for the proposed approach in comparison to other competitive methods, while the CPU time is also reasonable.

Concluding Remarks and Future Directions
A new four-step nonlinear method for solving f(x) � 0 type models is introduced in this research work with twelfth-order convergence, and seven function evaluations are required per iteration. e theoretical order of convergence for the proposed hybrid method is proved under both scalar and vector cases, along with an asymptotic error constant. Comparison with various existing numerical methods discloses the better performance of the proposed approach when the absolute errors, the absolute functional value computed at the last iteration, and the time of machine in seconds are taken into consideration.
e proposed method brings out the slightest absolute error regardless of initial conditions chosen for the simulations of the nonlinear models that belong to real-world scenarios from science and engineering. e rapid convergence of the proposed hybrid method is confirmed with the aid of polynomiography when the method is applied to some complex-valued polynomials. Declarations

Data Availability
All the data required for this paper is included within the paper.

Ethical Approval
We with this affirm that the contents of this article are original. Furthermore, it has been neither published elsewhere fully or partially in any language nor submitted for publication (fully or partially) elsewhere simultaneously. It contains no matter that is scandalous, obscene, fraud, plagiarism, libelous, or otherwise contrary to law.

Conflicts of Interest
e authors do not have any conflicts of interest.

Authors' Contributions
All authors contributed equally in this paper.