A Novel Three-Step Numerical Solver for Physical Models under Fractal Behavior

Abstract

In this paper, we suggest an iterative method for solving nonlinear equations that can be used in the physical sciences. This response is broken down into three parts. Our methodology is inspired by both the standard Taylor’s method and an earlier Halley’s method. Three evaluations of the given function and two evaluations of its first derivative are all that are needed for each iteration with this method. Because of this, the unique methodology can complete its goal far more quickly than many of the other methods currently in use. We looked at several additional practical research models, including population growth, blood rheology, and neurophysiology. Polynomiographs can be used to show the convergence zones of certain polynomials with complex values. Polynomiographs are produced as a byproduct, and these end up having an appealing look and being artistically engaging. The twisting of polynomiographs is symmetric when the parameters are all real and asymmetric when some of the parameters are imaginary.

1. Introduction

Nonlinear equations of the form

$$\Phi \left(x\right)=0,x\in \mathbb{R},$$

(1)

are commonplace in many subfields of science and engineering. The current era of computer research, in which one must acquire maximal results in the shortest amount of time [1], necessitates the development of innovative and effective iterative approaches for the resolution of nonlinear equations and systems. The new approaches are anticipated to be higher-order convergent; however, the efficiency of time consumption and the reduction in the amount of computing information required are the primary concerns. In the field of computational sciences, one of the most crucial and topical challenges is to find a solution to the nonlinear equations that regulates the natural events that arise in real-world problems. There are a great number of non-linear equations for which it is extremely difficult to obtain exact solutions. When dealing with these kinds of problems, it is important to come up with plans for getting close to a solution.

The investigation of productive approaches to the solution of nonlinear equations and systems appears to be one of the most significant challenges facing the fields of mathematics and engineering sciences. The solution of nonlinear equations and systems is not only a distinct subfield of mathematical research, but it also has an impact on the progression of mathematical knowledge in other subfields. For instance, the standard discretization processes for nonlinear differential and integral equations [2,3] result in systems of nonlinear equations that need to be solved by iterative nonlinear solvers. These systems of equations can then be analyzed and interpreted. Numerous studies have focused their attention on nonlinear equation systems that demonstrate real-world applications. For instance, there are issues with blood flow and rheology in [4,5]. In [6,7], a multitude of difficult and well-known application issues of the system of nonlinear equations have been addressed. The systems of nonlinear ordinary and partial differential equations [8,9,10,11,12], when taken in steady-state, are converted to nonlinear algebraic equations and thereby require numerical solvers such as the one proposed herein. These problems involve chemical equilibrium, neurophysiology, combustion, economic modeling, and kinematics. There are some nonlinear equations and systems that do not have exact solutions. These equations and systems arise from the mathematical modeling of physical systems.

For this reason, mathematicians have concentrated their efforts on developing nonlinear methods that are robust. Researchers have focused their attention on developing strategies that involve an improvement in the order of convergence, the utilization of a smaller number of operations, and a reduction in the number of function and higher-order derivative evaluations per iteration to create reliable modifications of conventional methods. The Taylor’s expansion, which was utilized for the Newton method, the Adomian decomposition, the Homotopy perturbation approach, the quadrature rules, and the variational iteration method are all examples of the many methods that can be employed to create novel nonlinear solvers [13,14,15,16,17,18]. The quadratic convergence of the second-order classical Newton method has led to it being given considerable weight. This is an old approach that can be used to solve nonlinear equations. Each repetition of the open-type Newton technique requires two evaluations: a functional evaluation and an evaluation of the first-order derivative to approximate the root. According to the hypothesis put forward by Kung and Traub [19], the order of convergence for the best iterative approach is $2^{\alpha -1}$, where $\alpha$ refers to the number of evaluations performed during each iteration. In this regard, the Newton technique might also be considered an ideal approach. Although the objective of current research studies has been to suggest improvements to the Newton method and its existing versions that are also optimal, it should be noted that approaches with a higher degree of convergence are also of significant practical value. Since the Kung and Traub argument is only a hypothesis, many studies have focused on improving the order in which techniques converge rather than on finding the optimal order. When ranking the different approaches, the efficiency index $(\rho )$ is another factor that must be considered. When figuring out a technique’s efficiency index, the number of evaluations that happen during each iteration ($\alpha$) and the order ($\omega$) in which the results converge are both taken into account. It may be noted that the efficiency index is computed as $\rho ={\omega }^{1/\alpha }$.

The following outline constitutes the framework of this study. In the next Section 2, you can find information on several existing schemes. After that comes the formulation of the suggested method (Section 3), followed by the determination of its order of convergence and the asymptotic error terms for both scalar and vector equations in Section 4. Some theoretical aspects are discussed, along with the suggested scheme’s stability analysis and basins of attraction in Section 5. In Section 6, there is a presentation of some numerical experiments that involve applied problems. Section 7, which talks about the most important findings and where future research should go, is where the work comes to a close.

2. Existing Iterative Methods

In this section, we will conduct a quick overview of a few well-known numerical schemes that are widely used to identify approximate solutions to nonlinear models. These schemes can be found in a variety of mathematical fields. There is absolutely no analytical method that can produce correct findings when it comes to the solution of such nonlinear models. As a result, the utilization of numerical strategies is required of us.

When someone thinks of solving a nonlinear equation with a numerical scheme, the first scheme that comes to mind is the traditional Newton–Raphson, which has optimal second-order convergence. This is because the Newton–Raphson scheme has been around for a long time. In a nutshell, we denoted this method with the symbol CNM2 since it employs two function evaluations every iteration: one function evaluation and one first-order derivative evaluation. The outline of the plan is as follows:

$$x_{n+1}=x_n-{\displaystyle \frac{\Phi \left(x_n\right)}{{\Phi }^{\prime }\left(x_n\right)}},n=0,1,2,\cdots ,$$

(2)

where ${\Phi }^{\prime }\left(x_n\right)\ne 0$. In 2021, the authors of [20] proposed a new eight-order method with five function evaluations for solving non-linear equations. The method is shown below and abbreviated as VNM8:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill & y_n=x_n-{\displaystyle \frac{\Phi \left(x_n\right)}{{\Phi }^{\prime }\left(x_n\right)}},\hfill \\ \hfill & z_n=y_n-{\displaystyle \frac{\Phi {\left(x_n\right)}^2\Phi \left(y_n\right)}{\Phi {\left(x_n\right)}^2{\Phi }^{\prime }\left(x_n\right)-2\Phi \left(x_n\right){\Phi }^{\prime }\left(x_n\right)\Phi \left(y_n\right)+{\Phi }^{\prime }\left(x_n\right)\Phi {\left(y_n\right)}^2}},\hfill \\ \hfill & x_{n+1}=z_n-{\displaystyle \frac{\Phi \left(z_n\right)}{{\Phi }^{\prime }\left(z_n\right)}},n=0,1,\cdots \hfill \end{array}\end{array}$$

(3)

In 2010, Xia Wang and Liping Liu [21] presented a new eight-order iterative method that required four function evaluations per iteration for solving nonlinear equations. The method is shown below and abbreviated as PNM8:

$$\begin{array}{cc}\hfill y_n& =x_n-{\displaystyle \frac{\Phi \left(x_n\right)}{{\Phi }^{\prime }\left(x_n\right)}},\hfill \\ \hfill z_n& =x_n-{\displaystyle \frac{\Phi \left(x_n\right)}{{\Phi }^{\prime }\left(x_n\right)}}{\displaystyle \frac{4\Phi {\left(x_n\right)}^2-5\Phi \left(x_n\right)\Phi \left(y_n\right)-\Phi {\left(y_n\right)}^2}{4\Phi {\left(x_n\right)}^2-9\Phi \left(x_n\right)\Phi \left(y_n\right)}},\hfill \\ \hfill x_{n+1}& =z_n-{\displaystyle \frac{\Phi \left(z_n\right)}{{\Phi }^{\prime }\left(x_n\right)}}\left[1+4{\displaystyle \frac{\Phi \left(z_n\right)}{\Phi \left(x_n\right)}}\right]\left[1+{\displaystyle \frac{8\Phi \left(y_n\right)}{4\Phi \left(x_n\right)-11\Phi \left(y_n\right)}}+{\displaystyle \frac{\Phi \left(z_n\right)}{\Phi \left(y_n\right)}}\right].\hfill \end{array}\begin{array}{c}n=0,1,\dots \hfill \end{array}$$

(4)

In 2020, Naseem, A. and et al. [22] presented a ninth-order iterative method that required six function evaluations per iteration for solving nonlinear equations. The method is shown below and abbreviated as PNM9:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill & y_n=x_n-{\displaystyle \frac{\Phi \left(x_n\right)}{{\Phi }^{\prime }\left(x_n\right)}}-{\displaystyle \frac{{\Phi }^2\left(x_n\right){\Phi }^{″}\left(x_n\right)}{2{\Phi }^{\prime 3}\left(x_n\right)}},\hfill \\ & z_n=y_n-{\displaystyle \frac{\Phi \left(y_n\right)}{{\Phi }^{\prime }\left(y_n\right)}},\hfill \\ \hfill & x_{n+1}=z_n-{\displaystyle \frac{\Phi \left(z_n\right)}{{\Phi }^{\prime }\left(y_n\right)-\beta \Phi \left(y_n\right)}},n=0,1,\cdots ,\hfill \end{array}\end{array}$$

(5)

where $\beta =1.$

In 2011, Hu et al. [23] presented a three-step iterative method of order nine that required five function evaluations for solving nonlinear equations. The method is shown below and abbreviated as PNR9:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill & y_n=x_n-{\displaystyle \frac{\Phi \left(x_n\right)}{{\Phi }^{\prime }\left(x_n\right)}},\hfill \\ & z_n=y_n-\left(1+{\left({\displaystyle \frac{\Phi \left(y_n\right)}{\Phi \left(x_n\right)}}\right)}^2\right)\left({\displaystyle \frac{\Phi \left(y_n\right)}{{\Phi }^{\prime }\left(y_n\right)}}\right),\hfill \\ \hfill & x_{n+1}=z_n-\left(1+2{\left({\displaystyle \frac{\Phi \left(y_n\right)}{\Phi \left(x_n\right)}}\right)}^2+2\left({\displaystyle \frac{\Phi \left(z_n\right)}{\Phi \left(y_n\right)}}\right)\right)\left({\displaystyle \frac{\Phi \left(z_n\right)}{{\Phi }^{\prime }\left(y_n\right)}}\right),n=0,1,\cdots \hfill \end{array}\end{array}$$

(6)

In 2020, Cordero, A. et al. [17] presented some one-step variants of Halley’s method with memory. One of such variants of order four that requires four function evaluations per iteration for solving nonlinear equations is shown below:

$$\begin{array}{cc}\hfill x_{n+1}& =x_n-{\displaystyle \frac{6\Phi \left(x_n\right){\Phi }^{\prime }{\left(x_n\right)}^2-3\Phi {\left(x_n\right)}^2{\Phi }^{″}\left(x_n\right)}{6{\Phi }^{\prime }{\left(x_n\right)}^3-6\Phi \left(x_n\right){\Phi }^{\prime }\left(x_n\right){\Phi }^{″}\left(x_n\right)+\Phi {\left(x_n\right)}^2{\Phi }^{‴}\left(x_n\right)}},\hfill \end{array}\begin{array}{c}n=0,1,\dots \hfill \end{array}$$

(7)

3. Proposed Iterative Technique

In this section, we have attempted to derive a nonlinear iterative method to solve nonlinear mathematical models of the type $\Phi \left(x\right)=0$, which is supposed to be a differentiable real-valued function. Let us further assume that $\gamma \in D\subset \mathbb{R}$ is a simple zero of $\Phi \left(x\right)=0$. Moreover, we assume that $x_n$ be an initial estimate which is close enough to the exact root $\gamma$. After being motivated by a recent research study carried out in [24], the cubic approximation of $\Phi \left(x\right)$ via Taylor’s expansion about the point $x_n$ can be written as:

$$\Phi \left(x\right)=\Phi \left(x_n\right)+(x-x_n){\Phi }^{\prime }\left(x_n\right)+{\displaystyle \frac{{(x-x_n)}^2}{2!}}{\Phi }^{″}\left(x_n\right)+{\displaystyle \frac{{(x-x_n)}^3}{3!}}{\Phi }^{‴}\left(x_n\right).$$

(8)

Neglecting the quadratic term in (8) and assuming $\Phi \left(x_{n+1}\right)=0$, we obtain the following relation:

$$\Phi \left(x_n\right)+(x_{n+1}-x_n){\Phi }^{\prime }\left(x_n\right)+{\displaystyle \frac{{(x_{n+1}-x_n)}^3}{6}}{\Phi }^{‴}\left(x_n\right)=0.$$

(9)

Solving the above equation for $x_{n+1}$ only in the linear part, we get

$$x_{n+1}=x_n-{\displaystyle \frac{\Phi \left(x_n\right)}{{\Phi }^{\prime }\left(x_n\right)}}-{\displaystyle \frac{{\Phi }^{‴}\left(x_n\right)}{6{\Phi }^{\prime }\left(x_n\right)}}{(x_{n+1}-x_n)}^3.$$

(10)

Substituting the variant of the Halley’s method given in (7) for $x_{n+1}$ only on the right hand side of (10), we obtain

$$\begin{array}{cc}\hfill & x_{n+1}=x_n-{\displaystyle \frac{\Phi \left(x_n\right)}{{\Phi }^{\prime }\left(x_n\right)}}-{\displaystyle \frac{1}{48}}\left[{\displaystyle \frac{{\Phi }^{‴}\left(x_n\right)\Phi {\left(x_n\right)}^3{\left(2{\Phi }^{\prime }{\left(x_n\right)}^2-\Phi \left(x_n\right){\Phi }^{″}\left(x_n\right)\right)}^3}{{\Phi }^{\prime }\left(x_n\right){\left({\Phi }^{\prime }{\left(x_n\right)}^3-\Phi \left(x_n\right){\Phi }^{\prime }\left(x_n\right){\Phi }^{″}\left(x_n\right)+\frac{1}{6}\Phi {\left(x_n\right)}^2{\Phi }^{‴}\left(x_n\right)\right)}^3}}\right].\hfill \end{array}$$

(11)

To get rid of the third-order derivative in (11), we use the following identity taken from [24] as:

$${\Phi }^{‴}\left(x_n\right)={\displaystyle \frac{6\Phi \left(x_n\right){\Phi }^{\prime }{\left(x_n\right)}^3-3\Phi {\left(x_n\right)}^2{\Phi }^{\prime }\left(x_n\right){\Phi }^{″}\left(x_n\right)}{2\Phi {\left(x_n\right)}^3}}.$$

(12)

Substituting (12) into (11), we get

$$\begin{array}{cc}\hfill & x_{n+1}=x_n-{\displaystyle \frac{\Phi \left(x_n\right)}{{\Phi }^{\prime }\left(x_n\right)}}-\hfill \\ \hfill & \left[{\displaystyle \frac{4\Phi {\left(x_n\right)}^3{\left(\Phi \left(x_n\right){\Phi }^{″}\left(x_n\right)-2\Phi {\left(x_n\right)}^2\right)}^3\left(-\frac{1}{2}\Phi {\left(x_n\right)}^2{\Phi }^{\prime }\left(x_n\right){\Phi }^{″}\left(x_n\right)+\Phi \left(x_n\right){\Phi }^{\prime }{\left(x_n\right)}^3\right)}{{\Phi }^{\prime }\left(x_n\right){\left(-5\Phi {\left(x_n\right)}^2{\Phi }^{\prime }\left(x_n\right){\Phi }^{″}\left(x_n\right)+4{\Phi }^{\prime }{\left(x_n\right)}^3\Phi \left(x_n\right)+2\Phi \left(x_n\right)+2\Phi \left(x_n\right){\Phi }^{\prime }{\left(x_n\right)}^3\right)}^3}}\right].\hfill \end{array}$$

(13)

Using the classical Newton’s method $z_n=x_n-{\displaystyle \frac{\Phi \left(x_n\right)}{{\Phi }^{\prime }\left(x_n\right)}}$, the iterative method given in [25], and the proposed one-step iterative method given in (13), we obtain the final form of the proposed three-step method as follows:

$$\begin{array}{cc}\hfill & z_n=x_n-{\displaystyle \frac{\Phi \left(x_n\right)}{{\Phi }^{\prime }\left(x_n\right)}},\hfill \\ \hfill & y_n=z_n-\left(1+{\displaystyle \frac{2\Phi \left(z_n\right)}{\Phi \left(x_n\right)-2\Phi \left(z_n\right)}}\right)\left({\displaystyle \frac{\Phi \left(z_n\right)}{{\Phi }^{\prime }\left(x_n\right)}}\right),\hfill \\ \hfill & x_{n+1}=y_n-{\displaystyle \frac{\Phi \left(y_n\right)}{{\Phi }^{\prime }\left(y_n\right)}}-\left[{\displaystyle \frac{2\Phi {\left(y_n\right)}^3{\left(\Phi \left(x_n\right)S-2\Phi {\left(x_n\right)}^2\right)}^3\left(\Phi {\left(y_n\right)}^2{\Phi }^{\prime }\left(y_n\right)S-2\Phi \left(y_n\right){\Phi }^{\prime }{\left(y_n\right)}^3\right)}{{\Phi }^{\prime }\left(y_n\right){\left(5\Phi {\left(x_n\right)}^2{\Phi }^{\prime }\left(x_n\right)S+4{\Phi }^{\prime }{\left(x_n\right)}^3\Phi \left(x_n\right)-2\Phi \left(x_n\right){\Phi }^{\prime }{\left(x_n\right)}^3\right)}^3}}\right].\hfill \\ \hfill n& =0,1,\cdots \hfill \end{array}$$

(14)

It may also be noted that the letter S in the third step of the above equations stands for the finite-difference approximation as ${\displaystyle \frac{{\Phi }^{\prime }\left(y_n\right)-{\Phi }^{\prime }\left(x_n\right)}{y_n-x_n}}$, which has been utilized for removing the second-order derivative ${\Phi }^{″}\left(x_n\right)$. In addition, the method in (14) is abbreviated as PTNM, which stands for Proposed Taylor’s expansion-based Numerical Method. The flowchart of the above eighth-order proposed three-step iterative method is shown in Figure 1. Moreover, the convergence order, the number of function evaluations per iteration, and the efficiency index of the proposed method are computed as 8, 5, and  1.5157, respectively. A comparative analysis of some methods under consideration based on these parameters is depicted in Figure 2.

4. Convergence Analysis

When a convergence theorem for an iterative method claims convergence while presuming the existence of a solution $x^{*}$ and a suitably chosen initial approximation close enough to $x^{*}$, it is referred to as a local convergence theorem. In contrast, a convergence theorem [26] pp. 575–600, does not presuppose the existence of any solution a prior but instead presupposes that specific criteria are true at the beginning point. A semilocal convergence theorem is what it is termed.

This section has been devoted to proving the convergence analysis based on the assumption of the existence of at least eighth-order derivatives of $\Phi \left(x\right)=0$. The necessary order of convergence has been obtained through the application of Taylor’s series expansion with a single variable. It is important to note that the convergence analysis is approached in a manner that is comparable to that of a large number of other previously published publications and that the primary interest in the creation of higher-order methods is of an intellectual nature. Even if higher-order algorithms are more difficult to understand, their efficiency can be assessed; this is why we have included the CPU time, which can be found in Section 6 of the article.

Theorem 1.

Suppose that $\gamma \in P$ is the exact root of a differentiable function $\Phi :P\subset \mathbb{R}\to \mathbb{R}$ for an open interval $P.$ Then, the three-step proposed method given in (14) has eighth-order convergence and the asymptotic error term is determined to be:

$$e_{n+1}=c_2^3(c_2^2-c_3)e_n^8+\mathcal{O}\left(e_n^9\right),$$

(15)

where $e_n=x_n-\gamma$, and $c_r={\displaystyle \frac{{\Phi }^{\left(r\right)}(\gamma )}{r!}}$, $r=1,2,3,\cdots .$

Expanding $\Phi \left(x_n\right)$ via series of Taylor about $\gamma$, we obtain

$$\Phi \left(x_n\right)=c_1{e}_n+c_2{e}_n^2+c_3{e}_n^3+\mathcal{O}\left(e_n^4\right).$$

(16)

Expanding ${\Phi }^{\prime }\left(x_n\right)$ via series of Taylor about $\gamma$, we obtain

$${\Phi }^{\prime }\left(x_n\right)=c_1+2c_2{e}_n+3c_3{e}_n^2+\mathcal{O}\left(e_n^4\right).$$

(17)

Applying Taylor series for function ${\displaystyle \frac{1}{{\Phi }^{\prime }\left(x_n\right)}}$ about $\gamma$, we have

$${\displaystyle \frac{1}{{\Phi }^{\prime }\left(x_n\right)}}=\frac{-3c_1{c}_3{e}_n^2+4c_2^2{e}_n^2-2c_1{c}_2{e}_n+c_1^2}{c_1^3}+\mathcal{O}\left(e_n^3\right).$$

(18)

Multiplying (16) and (18), we obtain:

$${\displaystyle \frac{\Phi \left(x_n\right)}{{\Phi }^{\prime }\left(x_n\right)}}=\frac{e_n\left(c_2{e}_n+c_1\right)\left(-3c_1{c}_3{e}_n^2+4c_2^2{e}_n^2-2c_1{c}_2{e}_n+c_1^2\right)}{c_1^3}.$$

(19)

Substituting (19) in the first-step of (14), we obtain

$${\widehat{e}}_n=\frac{e_n^2\left(3c_1{c}_2{c}_3{e}_n^2-4c_2^3{e}_n^2+3c_1^2{c}_3{e}_n-2c_1{c}_2^2{e}_n+c_2{c}_1^2\right)}{c_1^3},$$

(20)

where ${\widehat{e}}_n=y_n-\gamma$. Using the Taylor’s series for $\Phi \left(y_n\right)$ around $\gamma$, we obtain:

$$\Phi \left(y_n\right)=c_1{\widehat{e}}_n+c_2{\widehat{e}}_n^2+\mathcal{O}\left({\widehat{e}}_n^3\right).$$

(21)

Expanding ${\Phi }^{\prime }\left(y_n\right)$ via series of Taylor about $\gamma$, we obtain

$${\Phi }^{\prime }\left(y_n\right)=c_1+2c_2{\widehat{e}}_n+3c_3^2{\widehat{e}}_n+\mathcal{O}\left({\widehat{e}}_n^4\right).$$

(22)

Applying Taylor series for function ${\displaystyle \frac{1}{{\Phi }^{\prime }\left(y_n\right)}}$ about $\gamma$, we have

$${\displaystyle \frac{1}{{\Phi }^{\prime }\left(y_n\right)}}=\frac{-3c_1{c}_3{\widehat{e}}_n^2+4c_2^2{\widehat{e}}_n^2-2c_1{c}_2{\widehat{e}}_n+c_1^2}{c_1^3}+\mathcal{O}\left({\widehat{e}}_n^3\right).$$

(23)

Multiplying (21) and (23), we obtain:

$${\displaystyle \frac{\Phi \left(y_n\right)}{{\Phi }^{\prime }\left(y_n\right)}}=\frac{{\widehat{e}}_n\left(c_2{\widehat{e}}_n+c_1\right)\left(-3c_1{c}_3{\widehat{e}}_n^2+4c_2^2{\widehat{e}}_n^2-2c_1{c}_2{\widehat{e}}_n+c_1^2\right)}{c_1^3},$$

(24)

where ${\tilde{e}}_n=z_n-\gamma$. Using the Taylor’s series for $\Phi \left(z_n\right)$ around $\gamma$, we obtain

$$\Phi \left(z_n\right)=c_1{\tilde{e}}_n+c_2{\tilde{e}}_n^2+\mathcal{O}\left({\tilde{e}}_n^3\right).$$

(25)

Expanding the Taylor series using Equations (17) and (22), we obtain the expansion for the finite-difference quotient S as follows:

$$\begin{array}{cc}\hfill & S=2c_2+3c_3{e}_n+(-2c_2(c_2^3+c_2{c}_3)-2c_2(-c_2^3-c_2(4c_2^2-3c_3)-2c_2(-4c_2^2+2c_3)\hfill \\ \hfill & +2c_2(2c_2^2+2c_3)-2c_2(-c_2^2+2c_3)e_n^3+\mathcal{O}\left(e_n^4\right).\hfill \end{array}$$

(26)

Using the Equations (16), (18) and (25) in the second step of the proposed method (14), we obtain

$$\widehat{e}=(-c_2^3-c_2(4c_2^2-c_3)-2c_2(-4c_2^2+2c_3)+2c_2(-2c_2^2+2c_3)-2c_2(-c_2^2+2c_3)e_n^4+\mathcal{O}\left(e_n^5\right).$$

(27)

Using the Equations (16)–(18), (21)–(23) and (26) in the third step of the proposed method (14), we obtain

$$e_{n+1}=c_2^3(c_2^2-c_3)e_n^8+\mathcal{O}\left(e_n^9\right).$$

(28)

The aforementioned equation demonstrates that the suggested approach (14) converges to the eighth order.

5. Polynomiography: Fractal and Non-Fractal Images

Creative design and pattern development are two of the trickiest challenges in CAD. The process of creating patterns includes many different steps, such as analysis, creativity, and development. A designer must address all of these factors to create an engaging pattern that can be used in jewelry design, carpet design, tapestry design, etc. That is why it is so inspiring and practical to find novel ways to generate a wide range of fascinating patterns. Mathematical theory is one potential source for such procedures. As a type of mathematical object, polynomials have the potential to produce a wide range of interesting and lovely patterns. Polynomiography is commonly used to create patterns using polynomials.

Polynomiography is a method that combines mathematics with art to develop a new kind of visual art. The visuals that are produced are the outcome of an algorithmic visualization of repeated attempts to solve a polynomial problem. At the turn of the 21st century, Dr. Bahman Kalantari was the one who first articulated the meaning of this phrase [27]. Dr. Bahman Kalantari’s research on polynomial root-finding served as the impetus for the development of the concepts of polynomiography. Polynomial root-finding is an age-old and time-honored field of study, but it never ceases to produce new insights with the advent of new generations of mathematicians and scientists. The term “polynomiography”, which is a combination of the words “polynomial” and the suffix “graphy”, was initially developed by Dr. Kalantari. A picture that was generated independently as a consequence of polynomiography is referred to as a “polynomiograph”. “An iterative process for constructing two-dimensional colored drawings (polynomiographs) that represent polynomials”, is how its definition reads.

Fractals are connected to polynomials via an area of mathematics known as polynomiography. Polynomiography is described as “the art and science of visualizing the zeros of complex polynomials utilizing fractal and non-fractal pictures formed using the mathematical convergence qualities of iteration functions. Instead of being a fractal image, a polynomiograph is a visual depiction of a polynomial’s roots. Finding the roots of small to intermediate degree polynomials has received a lot of research, while big polynomials have received less attention since they are less frequently encountered in mathematical studies. The fractal patterns that were created were different and could be used in many different ways, such as in the textile and ceramics industries.

To generate polynomiographs over the complex plane $\mathbb{C}$ utilizing the computer program, we made use of a rectangle R along with the dimensions $[-2,2] \times [-2,2]$, and the maximum number of iterations of $N=50$ with accuracy $ϵ=0.001$. This was done by taking some complex polynomials as shown in Equation (29). The areas of the image where the algorithm did not successfully converge have been assigned the color black. The pixel density of the produced visual representations is defined by the dividing of R; for example, if we partition the rectangle R into a grid of 2000, the plotting polynomiographs would have a higher resolution as a result. We used the following four complex polynomials including their roots (both real and complex) to show graphic elements in the complex plane:

$$\begin{array}{cc}\hfill {\Phi }_1\left(z\right)& =z^3-1,z_{1,2,3}=1,-0.5\pm 0.8660i,\hfill \\ \hfill {\Phi }_2\left(z\right)& =z^4-1,z_{1,\dots ,4}=\pm 1,\pm i,\hfill \\ \hfill {\Phi }_3\left(z\right)& =z^5+z,z_{1,\dots ,5}=0,0.7071\pm 0.7071i,-0.7071\pm 0.7071i,\hfill \\ \hfill {\Phi }_4\left(z\right)& =45z^{15}+45z^8-1,\hfill \\ \hfill z_{1,\dots ,15}& =0.6187,-0.6243,-0.9967,0.4367\pm 0.4393i,-0.4423\pm 0.4393i,-0.0028\pm 0.6212i,\hfill \\ & 0.2257\pm 0.975i,-0.6203\pm 0.7819i,0.9041\pm 0.4339i.\hfill \end{array}$$

(29)

The images obtained are shown in Figure 3, Figure 4, Figure 5 and Figure 6 for the complex-valued polynomials given in Equation (29), respectively while considering the iterative techniques as follows:

(a)

Classical Newton Method of second order convergence in (2).

(b)

A method with eighth-order convergence in (3).

(c)

A method with eighth-order convergence in (4).

(d)

A method with ninth-order convergence in (5).

(e)

A method with ninth-order convergence in (6).

(f)

The proposed method with eighth-order convergence in (14).

It is worth to be noted that in each Figure, plot (a) shows the fractal behaviour of the classical Newton Method of second-order convergence, plot (b) shows the fractal behaviour of the eighth-order method given in (3), plot (c) shows the fractal behaviour of the eighth-order method given in (4), plot (d) shows the fractal behaviour of the ninth-order method given in (5), plot (e) shows the fractal behaviour of the ninth-order method given in (6), and plot (f) shows the fractal behaviour of the proposed eighth-order method given in (14). Darker regions in each plot indicate divergence regions that do not appear in fractal pictures generated using the suggested eighth-order approach. It demonstrates how, with a good first approximation, the approach can converge more quickly.

It may also be noted that the pseudo-code for the polynomiography via the proposed three-step method (14) is given in Algorithm 1. It has been explained how the different coloring is assigned in some regions and how the colormap works in case of different colors and shadings within the phase planes of several complex polynomials under consideration. Pseudo-code offers a way to find polynomiographs. The proposed three-step procedure is iterated in the algorithm for each point in the region $R\subset \mathbb{C}$ under consideration. We assume the resulting sequence converges to a root of p and stops iterating if the modulus of the difference between two successive points in the iteration process is less than the provided precision $(ϵ)$. If the maximum number of cycles (N) is reached, then the generated sequence is assumed to not converge to a root of p. At the very end, we assigned a color to the location that was being considered by utilizing the colormap that was presented to us together with the iteration number when we exited the while loop.    

6. Physical Models for Numerical Simulations

In this section, we will conduct some numerical simulations to highlight the performance of the suggested iterative technique with eighth-order convergence, which is given in the equation. There are many distinct kinds of nonlinear functions that are taken into account. Only applicable nonlinear models, such as those used in medical science, healthcare, population dynamics, and blood flow, are discussed in this article in order to provide the findings with more support and credibility. We utilize the stopping criterion $\left|\Phi \right(x_N\left)\right|\le \epsilon$ to determine when to end the number of iterations, where the tolerance is set to $\epsilon ={10}^{-300}$, and the needed accuracy is specified to be as high as 4000.

The following tabular representations of numerical results obtained from nonlinear practical models have been created, with each table containing an initial guess, a total number of function evaluations to achieve the tolerance, absolute error, the absolute value of the underlying function at the final approximated solution, and the amount of CPU time consumed by each approach, measured in seconds. The required numerical results were accomplished by using the program Mathematica 12.1 while operating under the Windows 10 operating system. The computer has an Intel(R) Core(TM) i7-1065G7 CPU running at a speed of 1.30 GHz. The installed memory was 24.0 GB.

Numerical Results and Discussions

Table 1. Comparison of several methods with the proposed method for the initial guess $x_0=5.5$ for the blood rheology model given in (30) having an approximate solution to be $x=3.8200$.

Method	FE	$\mathit{\epsilon }=|{\mathit{x}}_{\mathit{N}}-{\mathit{x}}_{\mathit{e}}|$	$|{\mathit{\Phi }}_1\left({\mathit{x}}_{\mathit{N}}\right)|$	Time
CNM2	28	1.48 $\times {10}^{-443}$	2.70 $\times {10}^{-884}$	1.4100 $\times {10}^{-1}$
PNM8	20	1.49 $\times {10}^{-492}$	1.93 $\times {10}^{-3933}$	2.8200 $\times {10}^{-1}$
VNM8	30	1.46 $\times {10}^{-1527}$	1.00 $\times {10}^{-3998}$	1.8700 $\times {10}^{-1}$
PNM9	30	1.16 $\times {10}^{-1895}$	1.00 $\times {10}^{-3998}$	2.0400 $\times {10}^{-1}$
PNR9	30	4.34 $\times {10}^{-316}$	9.26 $\times {10}^{-2835}$	1.2500 $\times {10}^{-1}$
PTNM	30	3.00 $\times {10}^{-2328}$	1.00 $\times {10}^{-3998}$	3.1200 $\times {10}^{-1}$

,

Table 2. Comparison of several methods with the proposed method with the initial guess $x_0=2$ for the blood rheology model given in (33) having an approximate solution to be $x=1$.

Method	FE	$\mathit{\epsilon }=|{\mathit{x}}_{\mathit{N}}-{\mathit{x}}_{\mathit{e}}|$	$|{\mathit{\Phi }}_2\left({\mathit{x}}_{\mathit{N}}\right)|$	Time
CNM2	46	5.56 $\times {10}^{-382}$	2.88 $\times {10}^{-762}$	1.0900 $\times {10}^{-1}$
PNM8	36	6.41 $\times {10}^{-882}$	3.00 $\times {10}^{-3999}$	1.2500 $\times {10}^{-1}$
VNM8	45	7.42 $\times {10}^{-857}$	7.00 $\times {10}^{-3999}$	7.8000 $\times {10}^{-2}$
PNM9	45	4.04 $\times {10}^{-845}$	2.00 $\times {10}^{-3999}$	7.8000 $\times {10}^{-2}$
PNR9	54	2.39 $\times {10}^{-1231}$	6.00 $\times {10}^{-3999}$	6.2000 $\times {10}^{-2}$
PTNM	45	1.09 $\times {10}^{-1514}$	7.00 $\times {10}^{-3999}$	1.4100 $\times {10}^{-1}$

and

Table 3. Comparison of several methods with the proposed method with the initial guess $x_0=5.5$ for the population growth model given in (36) having an approximate solution to be $x=1.01\times {10}^{-1}$.

Method	FE	$\mathit{\epsilon }=|{\mathit{x}}_{\mathit{N}}-{\mathit{x}}_{\mathit{e}}|$	$|{\mathit{\Phi }}_3\left({\mathit{x}}_{\mathit{N}}\right)|$	Time
CNM2	32	4.29 $\times {10}^{551}$	1.16 $\times {10}^{1098}$	2.0300 $\times {10}^1$
PNM8	diverge	-	-	-
VNM8	30	5.37 $\times {10}^{480}$	8.23 $\times {10}^{3834}$	2.5000 $\times {10}^1$
PNM9	30	1.90 $\times {10}^{508}$	3.10 $\times {10}^{3996}$	1.5600 $\times {10}^1$
PNR9	diverge	-	-	-
PTNM	30	2.82 $\times {10}^{995}$	6.00 $\times {10}^{3997}$	2.0300 $\times {10}^1$

show that the absolute error at the final iteration after simulating the chosen real models with the proposed method is the smallest of all the errors produced by other existing methods, as is the absolute functional value while consuming a reasonable amount of machine time measured in seconds. The number of function evaluations required by PTNM is appropriate for the first assumptions and, in some cases, it outperforms higher-order approaches. In addition, it is important to highlight that the suggested technique does not contain any aspects of divergence, in contrast to certain current methods, which do have such elements, as shown in Table 3. In conclusion, it has been demonstrated that the procedure described in (14) is applicable in a variety of real-life and physical circumstances where one must interact with a nonlinear phenomenon.

Problem 1.

We have taken a blood rheology nonlinear model from [4] as shown below:

$${\Phi }_1\left(x\right)={\displaystyle \frac{x^8}{441}}-{\displaystyle \frac{8x^5}{63}}-0.05714285714x^4+{\displaystyle \frac{16x^2}{9}}-3.624489796x+0.36,$$

(30)

where x shows the plug flow of Caisson fluid flow. Caisson fluid is used to represent blood despite the fact that it is a non-Newtonian fluid. The Caisson fluid model hypothesizes that a simple fluid like water or blood will move through a tube in such a way that its central core will move as a plug with very little distortion and that a velocity gradient will develop along the tube’s walls. It may be noted that the blood comes under the category of Caisson fluid. The numerical simulations for the model (30) are presented in Table 1.

Problem 2.

Law of Blood Flow [5]:

$${\Phi }_2\left(x\right)={\displaystyle \frac{P}{\eta l}}\left(R^2-x^2\right),$$

(31)

where $P=4000,\eta =0.027,R=0.008,$ and $l=2$ are chosen for the simulations with $x\in [0,R]$ being the distance to be determined. With these parameters, the above model turns out to be

$$R_f{x}^3-20p{(1-x)}^2=0,$$

(32)

where $R_f$ stands for the radius of the fiber, p shows the specific hydraulic permeability, and $x\in [0,1]$ is the porosity of the medium. If we assume $R_f=100 \times {10}^{-9}$ and $p=0.4655$, we obtain the following third-degree polynomial:

$${\Phi }_2\left(x\right)=-100 \times {10}^{-9}{x}^3+9.3100x^2-18.6200x+9.3100.$$

(33)

The numerical simulations for the model (33) are presented in Table 2.

Problem 3.

In the study of population dynamics, one of the models that is utilized most commonly is the mathematical representation of population growth, referred to as the population growth (PG) model. Using a first-order ordinary differential equation, the following expression can be used to describe the model:

$$P{G}^{\prime }\left(t\right)=\kappa PG\left(t\right)+\nu ,$$

(34)

where $PG\left(t\right)$ stands for the population growth at time t, κ shows the constant birth rate while ν is for constant immigration rate. The above linear ordinary differential equation has the following closed-form solution:

$$PG\left(t\right)=P{G}_{0}exp(\kappa t)+{\displaystyle \frac{\nu }{\kappa }}\left(exp(\kappa t)-1\right),$$

(35)

where $P{G}_0$ is the population at $t=0$ (initial population). In one particular piece of research, the problem is stated as follows: Let us say that a particular community had a starting population of 1,000,000 people, that 435,000 people immigrated into the community in the first year, and that the town had a total population of 1,564,000 people by the conclusion of the first year. Determine the rate of births $(\kappa )$ in this population. Using the initial condition and values of the other parameters given as stated, the birth rate κ can be determined with the help of the following nonlinear equation which is the closed-form solution of the model (34) after applying the stated situation:

$${\Phi }_3\left(x\right)=1564-{10}^{3}exp\left(x\right)-{\displaystyle \frac{435}{x}}\left(exp\left(x\right)-1\right)=0,$$

(36)

where $x=\kappa$ is the required birth rate. The numerical simulations for the model (36) are presented in Table 3.

7. Concluding Remarks and Future Directions

The success of this investigation hinges on the development of a brand-new iterative strategy that converges at the eighth-order level after three iterations. The projected computational order of convergence is consistent with the theoretical demonstration of convergence provided by Taylor’s series expansion for nonlinear equations. Thus, the PTNM method can be applied to problems with either a single nonlinear equation or a system of such equations. In addition, dynamical aspects of PTNM can be investigated using basins of attraction, which, when applied to complex-valued functions, yield aesthetically pleasing phase plane diagrams. This demonstrates the method’s validity for making first estimates that are statistically very close to the underlying nonlinear model. We also get polynomiographs that are both aesthetically pleasing and intellectually stimulating. The real values of parameters have a symmetrical effect, but the imaginary values of parameters distort polynomials. In the end, some medically-relevant nonlinear equations were applied to PTNM and other well-known optimum and non-optimal techniques in the sense of King–Traub. It has been found that PTNM yields better results than other methods in the vast majority of cases. This is notably the case concerning the N-th iteration threshold for achieving the target precision, error, and functional value. It is also worth noting that the proposed method will always converge, regardless of how close to or far from the starting guess the approximate solution is. Remembering this is essential. The PTNM algorithm converges to the eighth order and is an iterative method that can compete with other algorithms. It can be used to find solutions to nonlinear systems and equations.

As far as we can tell, high-order approaches are solely of interest in the academic community since very precise approximations to answers are not needed in most real-world contexts. However, these methods are rather complex and do not significantly increase productivity relative to simpler alternatives. The method proposed here is also one of the memory-free techniques discussed in the article. In the case of nonlinear systems, it also requires the evaluation of additional Jacobian matrices, which increases the computing burden. To reduce the computational complexity of the current method, we plan to make a change in future research by proposing a finite-difference approximation for the first-order derivative. We will also look more closely at the proposed method to see if it has semi-local convergence and better visualization with a variety of convergence tests including distance and non-distance conditions over different types of metrics. In addition, the three-step numerical solver that was proposed can be extended to the concept of fractional calculus [28], in which the traditional derivative is replaced by the fractional Riemann–Liouville or Caputo derivatives. This allows the solver to be used to solve problems involving fractional calculus.