Algorithms of predictor-corrector type with convergence and stability analysis for solving nonlinear systems

1.   Introduction

In the field of applied mathematics, solving algebraic and non-algebraic nonlinear equations is an open challenge for mathematicians and scientists. The exact solution of many nonlinear equations is nearly impossible. Therefore, it is important to develop new algorithms for dealing with such equations. Hence, researchers always keep searching for new numerical algorithms for getting the approximate solution of the nonlinear equations of the following type (if taken in univariate case):

where $\psi$ is a Fréchet-differentiable operator defined on a nonempty, open convex subset $\mathcal{S}$ of a Banach space $\mathcal{B}$ with values in a Banach space $\mathcal{B'}$. Similarly, if taken into the multivariate case (a system of $n$ nonlinear equations in $n$ unknowns) then one will have the following structure in its vector form:

where ${\boldsymbol{\psi}}:\mathbb{R}^n \rightarrow \mathbb{R}^n$.

The Newton-Raphson (NR) algorithm is a common iterative algorithm for finding approximate solutions to nonlinear equations with variables that have real values ^[1]. The fundamental notion of this concept is centered on the use of linear approximations. When presented with an initial estimate for the root of a function, denoted as $\psi(x) = 0$, it is possible to create a linear approximation of the function in the vicinity of this estimate by employing the Taylor series. This approximation facilitates the formulation of the subsequent estimate for the root, which is then iteratively revised to approximate the solution. The algorithm commences by establishing an initial approximation and thereafter iteratively enhances it through the evaluation of the function and its derivative. The aforementioned approach is notable for its expeditious convergence in close proximity to the true root, its simplicity, and its wide-ranging applicability across functions.

Nevertheless, the NR approach does have its limitations. It is worth noting that the convergence of the algorithm is not always assured and is strongly dependent on the initial estimate. The computational demand or complexity of computing the derivative at each iteration can be significant for certain functions. Moreover, it is worth noting that the procedure may encounter difficulties or deviate from its intended path when launched in proximity to an inflection point or when the derivative approaches zero. When confronted with a function that has numerous roots, the outcome of the approach can differ, as it may converge to any of these roots depending on the initial approximation. To enhance its effectiveness, it is crucial to carefully choose the initial estimate and be aware of the potential constraints of the procedure. Adjustments, like introducing a damping factor, can also be incorporated to enhance its robustness. The classical Newton's technique to find the approximate solution of the non-linear equation (1.1) is given as:

This second-order one-step approach has two function evaluations per iteration. The theoretical limit of the one-point techniques is beaten by the multi-step iterative approaches of both higher convergence order and better efficiency index. The formula for the efficiency index is $\kappa = r^{1/p}$, where $r$ is the order of the algorithm and $p$ is the number of function evaluations per iteration. The authors proposed solving a system of two equations to solve a scalar problem in ^[2]. So, the linked system's solutions on the identity line solve the problem. In problematic circumstances where the scalar Newton approach fails, this strategy succeeds. The authors in ^[3], proposed an alternate algorithm for solving nonlinear equation systems when Newton's algorithm fails. Based on the paper ^[2], the proposed technique extends systems. The theory behind Newton's algorithm is used to solve an associated system to approximate a system of equations. Many such iterative algorithms have been published.

Researchers have focused on improving the convergence order and reducing function evaluations, leading to the development of multi-step iterative algorithms. The authors introduced a three-step iterative nonlinear approach for nonlinear equations and systems in the paper ^[4]. Iterations of the proposed technique need three function evaluations and two first-order derivative checks. The theory proves the proposed strategy is sixth-order convergent. The suggested algorithm is compared to others based on error distributions, computational efficiency, and CPU times. Qureshi et al. in ^[5] combined existing algorithms to create an optimal fourth-order algorithm for solving scalar and vector forms of problems. A new root-finding algorithm ^[6] that combines forward and finite-difference methods, is efficient, derivative-free, and cheaper per iteration. Quantitative and graphical studies reveal that the approach is quintic-order convergent and outperforms previous algorithms. An eighth-order three-step algorithm after merging second-order NR with an existing fourth-order algorithm has been devised in ^[7] that is later proven to be efficient while solving single-variable nonlinear physical models. In a research study ^[8], a twelfth-order numerical algorithm was proposed to solve the nonlinear equations of the type (1.1), however, the algorithm was four-step and thus required seven function evaluations per iteration. Algorithms for solving nonlinear algebraic equations are also needed while discretizing ordinary and partial differential equations ^[9,10,11,12].

The following are some of the well-known iterative algorithms having fourth- and eighth-order convergence selected for comparison with algorithms developed in this research work. For example, Hueso et al. in ^[13] proposed an iterative fourth-order algorithm, i.e., Jarratt's method. Newton's step is used as a predictor of the following form and abbreviated as PJNM:

where $j = 1, 2, \ldots$.

Behl in 2015 ^[14] constructed an optimal algorithm of fourth-order convergence of the following form and abbreviated as KTNM:

where $j = 1, 2, \ldots$.

Similarly, Özban and Kaya in ^[15] proposed two iterative fourth-order algorithms of the following form and abbreviated as CNM1 and CNM2:

and

where $j = 1, 2, \ldots$.

In 2021, Kong-ied proposed a new eight-order algorithm ^[16], abbreviated as ONM:

where $j = 1, 2, \ldots$.

Motivated by the current studies in this direction, we attempt to propose two algorithms having eighth- and fourth-order convergence for solving nonlinear models of the type (1.1). The development of the algorithms was extensively supported by alternating Taylor's series expansion and the classical Halley algorithm. Moreover, the developed algorithms are equally applicable to both univariate and multivariate cases.

This article is structured as follows: In Section 2, we present the basic construction of both algorithms including their local convergence analysis based on Taylor's expansion, while Section 3 presents the local and semilocal convergence analysis for the proposed algorithms under consideration. A detailed visual analysis via the polynomiography of the algorithms is presented in Section 4. To prove the better performance of the proposed algorithms, some numerical experiments (single-variable and system), including both academic and real-life situations, are conducted in Section 5, whereas the concluding remarks with some future directions are described in Section 6.

2.   Derivation of the proposed algorithms

In this section, we present the construction of the proposed two algorithms. Moreover, we perform the convergence analysis of the algorithms via Taylor's expansion.

2.1.   Construction of the eighth-order algorithm

Suppose that $\psi$ is a real function defined on the interval $W \subset \mathbb{R}$. Let $\xi \in W$ be an exact simple zero for the non-linear equation $\psi(x) = 0$ and further let $x_j$ be an initial guess close to the exact root of the nonlinear equation. By the Taylor series expansion, the function $\psi$ is expanded to three terms around the point $x_j$ as follows:

For getting the next approximation $x_{j+1}$ for the root of $\psi(x) = 0$ in Eq (2.1), we assume that $\psi(x_{j+1}) = 0$. Hence, we get

After simplification, we get

From ^[17], we have

Substituting Eq (2.4) into Eq (2.3), we get

After simplification, we get

Using the second-order Newton-Raphson algorithm as a predictor and the Eq (2.6) as a corrector, we develop a predictor-corrector (PC) type algorithm as given below:

where $j = 1, 2, \ldots$. The new PC numerical algorithm given by (2.7) has an eighth-order of convergence and is denoted as PCNM8. The proposed PCNM8 algorithm has five function evaluations (two evaluations of functions, two first-order derivatives, and one second-order derivative), and the efficiency index of this algorithm is $r^{1/5} \approx 1.5157$, where $r$ stands for the order of the algorithm. It may be noted that the efficiency index of the proposed eighth-order algorithm is higher than the classical NR algorithm.

2.2.   Construction of the fourth-order algorithm

In several research studies, researchers have employed different strategies, including finite differences, weight functions, regularization, and the secant algorithm, among others; to get rid of higher-order derivatives. Following the same approach, we attempt to use several possible approaches to remove the second-order derivative in the above-developed eighth-order algorithm in Eq (2.7). For this purpose, the following structure:

for the second derivative taken from ^[18] is chosen due to its simplicity and better simulation results.

Substituting Eq (2.8) into Eq (2.7), we get

Hence, the algorithm (2.9) is a modified algorithm that is now free from the second derivative and abbreviated as PCNM4.

2.3.   Convergence analysis of proposed algorithms via Taylor expansion

Theorem 1. Let $\xi \in W$ be a simple root of a differential function $\psi: W \subset \mathbb{R} \to \mathbb{R}$, where $W$ is an open interval. If $x_0$ is an initial guess considerably near $\xi$, then the modified algorithm defined by (2.7) has eighth-order convergence and satisfies the following error equation:

where $e_j = x_j - \xi$.

Proof. Since $\xi$ is a root of $\psi$ and $e_j = x_j- \xi$ is the error at $j^{th}$ iteration, we can expand $\psi(x_j)$ in power of $e_j$ by Taylor's series expansion as follows:

By Taylor's series for $\frac{1}{\psi'(x_j)}$ about $\xi$, we obtain

Multiplying (2.11) and (2.12), we get

Using (2.13) in the first step of (2.7), we get

where $\sigma_j = y_j-\xi$. By Taylor's series for $\psi(y_j)$ about $\xi$, we obtain

By Taylor's series for $\frac{1}{\psi'(y_j)}$ about $\xi$, we obtain

Multiplying (2.15) and (2.16), we get

By Taylor's series for $\psi'(y_j)$ about $\xi$, we obtain

By Taylor's series for $\psi''(y_j)$ about $\xi$, we obtain

Squaring (2.15), we have

Squaring (2.18), we have

Using the values of $\psi^2(y_j)$, $\psi'(y_j)$ and $\psi''(y_j)$ in numerator, i.e., $2 \psi^2(y_j)$$\psi'(y_j)$$\psi''(y_j)$ of the second step of (2.7), we get

Using the values of $\psi(y_j)$, $\psi'^2(y_j)$ and $\psi''(y_j)$ in denominator, i.e., $[2 \psi'^2(y_j)-\psi(y_j)\psi''(y_j)]^{2}$ of the second step of (2.7), we get

Finally, by using (2.18), (2.22) and (2.23) in the second step of Eq (2.7), we get

By substituting $\sigma_j = \frac{\psi''(\xi) e_j^2}{2\psi'(\xi)} + \frac{\psi''(\xi)^2 e_j^3}{2\psi'(\xi)^2} + \mathcal{O}(e_j)^4$ in (2.24), we get

Hence, Eq (2.25) shows that the proposed algorithm (2.7) has eighth-order of convergence. □

As the convergence analysis carried out for the proposed eighth-order algorithm (2.7), the convergence analysis for the proposed fourth-order algorithm is addressed in the same way.

Theorem 2. Let $\xi \in W$ be a simple root of a differential function $\psi: W \subset \mathbb{R} \to \mathbb{R}$, where $W$ is an open interval. Let $x_0$ be an initial guess considerably near to the exact root $\xi$. Then, the modified algorithm defined by (2.9) has fourth-order convergence and satisfies the following error equation:

Proof. Since $\xi$ is a root of $\psi$ and $e_j = x_j-\xi$ is the error at $j^{th}$ iteration, we can expand $\psi(x_j)$ in power of $e_j$ by Taylor's series expansion as follows:

By Taylor's series for $\frac{1}{\psi'(x_j)}$ about $\xi$, we obtain

Multiplying (2.27) and (2.28), we get

Using (2.29) in the first step of (2.9), we get

By Taylor's series for $\psi(y_j)$ about $\xi$, we obtain

By finding the derivative of $\psi(x_j)$ and $\psi(y_j)$ respectively, we have

By Taylor's series for $\frac{1}{\psi'(y_j)}$ about $\xi$, we obtain

Multiplying (2.31) and (2.34), we get

By putting all the above equations into second step of (2.9), we get

The error equation (2.36) shows that the proposed algorithm (2.9) has fourth-order convergence. □

The flowchart of the proposed two-step iterative algorithms is shown in Figure 1.

3.   Local and semi-local convergence

There are certain limitations with the local convergence analysis performed in Section 2 for the algorithm (2.7).

$(L_1)$ The algorithm (2.7) can be also used to solve equations on the real line.

$(L_2)$ The existence of at least fourth derivative is required to show convergence, and the solution has to be simple. Consider the example, for $W = [-1.5, 1.5],$ and function $\psi$ defined on $W$ by $\psi(x) = x^4\log(x)+x^5-x^4$ for $x\neq 0,$ and $\psi(0) = 0$. Then, clearly $\zeta = 1$ solves the equation $\psi(x) = 0,$ but $\psi^{(4)}$ does not exist at $x = 0$. Consequently, the results of Section 2 cannot guarantee convergence to $\zeta$ although algorithm converges for say $x_0 = 1\in W$.

$(L_3)$ There are no commutable error bounds on the distances $||x_j-\zeta||$. That is, we do not know a priori how many iterates are needed to reach a certain error tolerance.

$(L_4)$ There are no uniqueness of the solution results.

$(L_5)$ The more interesting semi-local convergence is not considered in Section 2.

$(L_6)$ There is no radius of convergence. Thus, we do not know how to pick the initial point $x_0$.

We positively address the limitations as follows:

$(L_1)'$ The convergence is given for Banach space valued equations.

$(L_2)'$ The convergence conditions involve only the operators on the algorithm (3.1), i.e., $F$, $F'$ and $F''$.

$(L_3)'$ Computable error bounds on $||x_j-\zeta||$ are provided. Hence, we know in advance the number of iterations to be executed to obtain a certain accuracy.

$(L_4)'$ A certain region is specified containing only one solution of the equation $Q(x) = 0.$

$(L_5)'$ The semi local convergence analysis of algorithm (2.7) is developed using majorizing sequences.

$(L_6)'$ The radius of convergence is specified. This is how we extend the applicability of the algorithm (2.7).

The algorithm (2.7) in a Banach space setting is given for starter $x_0\in D$, and each $j = 0, 1, 2, \ldots$ by

where the operators $Q'$ and $Q''$ denote the first and the second Fréchet derivative of $Q$, respectively, to solve the equation $Q(x) = 0$, where $Q : D \subset W_1 \rightarrow W_2$ is a twice Fréchet differentiable operator, $D$ is open and convex set and $W_1, W_2$ are Banach spaces.

Notice that $Q'$ is linear and $Q''$ is a bilinear operator ^[19].

It is worth noting that (2.9) cannot be written in Banach space, since the long denominator is not a linear operator.

3.1.   Local analysis

The local convergence analysis is based on some conditions.

Suppose:

$(C_1)$ There exists a function $\kappa_0:[0, +\infty)\rightarrow \mathbb{R}$ which is continuous and also non-decreasing (CFND) such that the equation $\kappa_0(t)-1 = 0$ has a smallest solution which is positive (SSP). Denote such solution by $R_0$. Set $T_0 = [0, R_0)$.

$(C_2)$ There exists CFND $\kappa:T_0 \rightarrow \mathbb{R}$ such that for $h_1:T_0 \rightarrow \mathbb{R}$ defined by

the equation $h_1(t)-1 = 0$ has a SSP. Denote such solution by $t_1.$

$(C_3)$ The equation $\kappa_0(h_1(t))-1 = 0$ has a SSP in the interval $T_0.$ Denote such solution by $R_1.$ Set $T_1 = [0, R_1)$.

$(C_4)$ For $h_2:T_1\rightarrow \mathbb{R}$ defined by

the equation $h_2(t)-1 = 0$ has a SSP. Denote such solution by $t_2$.

$(C_5)$ Let $\alpha, \beta > 0$ be given constants, and $\kappa_2:T_1\rightarrow \mathbb{R}$ be an CFND. For $p:T_1 \rightarrow \mathbb{R}$ defined by

the equation $p(t)-1 = 0$ has SSP. Denote such solution by $R_2$.

Set $T_2 = [0, R_2)$. Consider function $q:T_2\rightarrow \mathbb{R}$ defined by

$(C_6)$ The equation $\kappa_0(t h_2(t))-1 = 0$ has an SSP. Denote such solution by $R_2$. Set $T_3 = [0, R_3)$.

$(C_7)$ Let $\kappa_1 :T_3 \rightarrow \mathbb{R}$ and define functions $\bar{\kappa}:T_3 \rightarrow \mathbb{R}$ and $h_3 :T_3 \rightarrow \mathbb{R}$

and

The equation $h_3(t)-1 = 0$ has an SSP. Denote such solution by $t_3$.

Set

and $T = [0, t_0)$. Then it follows by these conditions that for each $t\in T$

and

The constant $t_0$ is proven to be a radius of convergence for the algorithm (3.1) in Theorem 1. First, we relate the developed functions $\kappa, \kappa_0, \kappa_1$ and $\kappa_2$ to the operators on the algorithm (3.1).

$(C_8)$ There exists a solution $\zeta \in D$ of the Eq (3.1), an invertible linear operator $M$ such that $||M|| \leq \alpha$ and $||M^{-1}||\leq \beta$.

$(C_9)$ $||Q'(x)-M||\leq \frac{1}{\beta} \kappa_0(||x-\zeta||)$ for each $x\in D.$ Define the region $D_0 = D \cap S(\zeta, R_0)$, where $S(\zeta, R_0)$ stands for a ball that is open having center $\zeta$ and radius $R_0 > 0$. It follows by this condition and $(C_1)$ that for $x = \zeta$

Thus, from the standard lemma by Banach on linear, and invertible operators ^[20,21] $Q'(x)$ is invertible, and

$(C_{10})$ $||Q'(x)-Q'(y)|| \leq \frac{1}{\beta} \kappa(||x-y||)$ for each $x, y \in D_0.$

Set $D_1 = D \cap S(\zeta, R_3)$.

$(C_{11})$ $||Q'(x)|| \leq \frac{1}{\beta}\kappa_1(||x-\zeta||)$ and $||Q''(x)|| \leq \frac{1}{\beta}\kappa_2(||x-\zeta||)$ for each $x\in D_1$.

$(C_{12})$ $S[\zeta, t_0]\subset D$, where $S[\zeta, t_0]$ stands for the closure of $S(\zeta, t_0).$

Notice that the first condition in $(C_{11})$ can be dropped, since

Hence, the function $\kappa_1$ can be defined by $\kappa_1(t) = \alpha + \frac{1}{\beta}\kappa_0(t).$ Moreover, the function $\bar{\kappa}$ given in the condition $(C_7)$ is chosen in particular to be the smallest of the functions $\kappa_1$ and $\alpha + \frac{1}{\beta}\kappa_0(t).$ The developed notation and the conditions $(C_1)$–$(C_{12})$ shall be used in the local analysis of convergence for the algorithm (3.1).

Theorem 3. Suppose that the conditions $(C_1)$–$(C_{12})$ are satisfied, and $x_0 \in S_0 = S(\zeta, t_0) \setminus \{\zeta\}$. Then, the following assertions hold for each $j = 0, 1, 2, \ldots$

and the sequence $\{x_j\}$ is convergent to $\zeta.$

Proof. By assuming $x_0\in S_0,$ the assertion (3.9) is satisfied for $j = 0$. Let $x\in S(\zeta, t_0)$. Then, by the condition $(C_1)$, $Q'(x)$ is invertible, and the estimate (3.8) holds. In particular, for $x = x_0$ the linear operator $Q'(x_0)$ is invertible, since $x_0\in S_0$. Thus, the iterate $y_0$ is well defined. Then, we can write by the first substep of the algorithm (3.1) for $j = 0$:

By (3.8), the condition (3.9), the definition of the function $h_1$, (3.7), and (3.13), we have in turn

so the assertion (3.10) is satisfied, the iterate $y_0\in S(\zeta, t_0)$, and consequently (3.8) holds if $x = y_0$. Thus, the iterate $z_0$ is well defined. By the second substep of the algorithm (3.1) for $j = 0$:

leading as in (3.13) and (3.14) to

showing the assertion (3.11), the iterate $z_0 \in S(\zeta, t_0)$, and the estimate (3.8) holds for $x = z_0$. We shall show the invertability of the operator $A_0$, which implies the existence of the iterate $x_1$. Consider the identity

In view of $(C_7)$, the second condition in $(C_{11})$, $(C_{5})$, (3.5)–(3.7), and (3.17), we get

where we also use $S_0||A_0^{-1}||\leq q_0$, (3.8)

Thus,

Then, we can write by the third substep of algorithm (3.1)

It follows by $(C_8)$–$(C_{11})$, (3.7) (for $j = 3$), (3.14), (3.15), (3.19), and (3.20) that

showing the assertion (3.12) and $x_1 \in S(\zeta, t_0)$. Simply switch $x_0, y_0, z_0, x_1$ by $x_k, y_k, z_k, x_{k+1}$ in the preceding calculations to complete the induction for the assertions (3.9)–(3.12). Then, from the estimation

we deduce that $x_{k+1}\in S(\zeta, t_0)$ and

□

Next, a neighborhood of $\zeta$ is determined that contains no other solution.

Proposition 1. Suppose there exists a solution $\zeta\in S(\zeta, t_4)\subset D$ of the equation $Q(x) = 0,$ such that the condition $(C_9)$ holds on $S(\zeta, t_4)$, and there exists $t_5\geq t_4$ so that

Define the region $D_2 = D \cap S[\zeta, t_5].$ Then, the only solution of the Eq (1.1) in the region $D_2$ is $\zeta$.

Proof. Let $E = \int_0^1Q'(\zeta+\theta(u-\zeta))d\theta$ for some $u\in D_2$ with $Q(u) = 0.$ By applying the condition $(C_9)$, and the assumption (3.22), we obtain in turn

so $E$ is invertible. Then, from the identity

we deduce that $u = \zeta.$ □

3.2.   Semi-local analysis

This analysis is based on majorizing sequences. The role of $x_0$ and the $\kappa$ function is replaced by $\zeta$ and the $\vartheta$.

Suppose:

$(H_1)$ There exists CFND $\vartheta_0:[0, +\infty)\rightarrow \mathbb{R}$ such that equ ation $\vartheta_0(t)-1 = 0$ has an SSP. Denote such solution by $\delta_0$.

$(H_2)$ There exist CFND $\vartheta, \vartheta_1, \vartheta_2:[0, \delta_0)\rightarrow \mathbb{R}.$ Define the sequences $\{a_j\}, \{b_j\}, \{c_j\}$ for $a_0 = 0,$ some $b_0\geq 0$, and each $j = 0, 1, \ldots$, by

These scalar sequences are shown to be majorizing for the algorithm (3.1) in Theorem 2. However, a convergence condition for these sequences is needed.

$(H_3)$ There exists $\delta\in [0, \delta_0)$ such that for each $j = 0, 1, 2, \ldots$

It follows by this definition and (3.23) that $0\leq a_j \leq b_j \leq c_j \leq a_{j+1}\leq \delta,$ and there exists $a^*\in [0, \delta]$ such that

The functions $\vartheta$ relate to the operator on the algorithm.

$(H_4)$ There exist constants $\alpha, \beta > 0$, and an invertiable operator $M$.

$(H_5)$ For some $x_0\in D,$ and each $x\in D$

Set $D_3 = D \cap S(x_0, \delta_0)$.

Notice that if $x = x_0,$ we get $||M^{-1}||\, ||Q'(x_0)-M|| \leq \vartheta_0(0) < 1,$ by the definition of $\delta_0$. Thus, the linear operator $Q'(x_0)$ is invertible. Set $||Q'(x_0)^{-1}Q(x_0)||\leq b_0$.

$(H_6)$ For each $x, y\in D_3$

$(H_7)$ $S[x_0, a^*]\subset D$.

Next, the semilocal analysis of convergence for the algorithm (3.1) is developed under the conditions $(H_1)$–$(H_7)$.

Theorem 4. Suppose that the conditions $(H_1)$–$(H_7)$ hold. Then, the sequence $\{x_j\}$ generated by the algorithm (3.1), exists in $S(x_0, a^*)$, stays in $S(x_0, a^*)$ for each $j = 0, 1, 2, \ldots$, and is convergent to a solution $\zeta \in S[x_0, a^*]$ of the equation $Q(x) = 0,$ such that for each $j = 0, 1, 2, \ldots$

Proof. Mathematical induction is utilized to show the assertions

The assertion (3.26) holds if $j = 0,$ since by the first substep of algorithm (3.1), and the definition of $b_0$:

Thus, the iterate $y_0\in S(x_0, a^*)$. We can write by the first substep of (3.1) that

Using the first condition in $(H_6)$

By the condition $(H_4)$ for $x = y_0$

Thus, $Q(y_0)$ is invertible, and

and the iterate $z_0$ is well defined by the second substep of (3.1). Then, we can also write

In view of (3.23) and (3.29)–(3.31), we get

and

Hence, the assertion (3.29) holds if $j = 0$, and $z_0\in S(x_0, a^*).$

As in the local case:

so, $A_0$ is invertible, and

Notice, also that $Q'(z_0)$ is invertible, so the iterate $x_1$ is well defined by the third substep of (3.1). Then, we have by (3.1)

leading to

and

Thus, the assertion (3.30) holds for $j = 0,$ and the iterate $x_1\in S(z_0, a^*)$. We can write

leading to

so

and

The induction for the assertions (3.25)–(3.27) is completed if $x_0$, $y_0$, $z_0$, $x_1$ are replaced by $x_m$, $y_m$, $z_m$, $x_{m+1}$, respectively. Then, we have

Consequently, the sequence $\{x_m\}$ is Cauchy in a Banach space and such it is convergent to some $\zeta\in S[x_0, a^*]$. By letting $m\rightarrow +\infty$ in (3.33) we get $Q(\zeta) = 0.$ Moreover, by the estimate

we obtain (3.24) by letting $k\rightarrow +\infty$ in (3.34). □

As in the local case a neighborhood of $x_0$ is special containing only one solution in the next result.

Proposition 2. Suppose there exists a solution $y^*\in S(x_0, t_7)$ of the Eq (3.1).

The condition $(H_5)$ holds in ball $S(x_0, t_7),$ and there exists $t_8\geq t_7$ such that

Define the region $D_4 = D \cap S[x_0, t_8].$

Then, the only solution of the equation $Q(x) = 0$ in the region $D_4$ is $y^*$.

Proof. Let $E_1 = \int_0^1 Q'(y^*+\theta (z^*-y^*)d\theta$ for some $z^*\in D_4$ with $Q(z^*) = 0.$ It follows by the condition $(H_5)$ and (3.36)

Therefore, the operator $E_1$ is invertible. Then, from the identity

we deduce that $z^* = y^*.$ □

Remark 1. (i) Popular choices for $M = Q'(\zeta)$ (local case) and $M = Q'(x_0)$ (semilocal case). However, these choices are not the necessary the most flexible. Other choices exist ^[22,23].

(ii) Note that not all the conditions $(H_1)$–$(H_7)$ are assumed in Proposition 1. However, if they are, we can take $y^* = \zeta,$ and $t_7 = a^*$.

(iii) The constant $\delta_0$ can replace the limit point $a^*$ is the condition $(H_7)$.

4.   Visual analysis via polynomiography

The analysis of stability and speed of convergence is an essential aspect of analyzing every root-finding algorithm. The standard way of performing such an analysis is the use of visual approach via the so-called polynomiography ^[24,25]. In polynomiography, we generate images called polynomiographs using a general algorithm presented in Algorithm 1. Depending on the coloring algorithm in the last step of the algorithm, we can visualize various aspects of the root-finding algorithm.

In this section, we use a coloring algorithm that shows basins of attraction and speed of convergence in the same polynomiograph ^[26]. In this algorithm, each root of the considered polynomial gets a distinct color. We also need an additional color to depict the non-convergent points (we use the black color for this aim). Now, to color the given starting point $z_0$, we first check whether the root-finding algorithm has converged to any of the roots, i.e., if $n < N$, where $n$ is the number of the last iterate. If the algorithm has converged, then, for the found root, we find the closest root of the polynomial and use its color to color $z_0$, else we use the additional color. In this way, we show basins of attraction. To add the information on the speed of convergence, we change the shade of the color based on the number of performed iteration. In this way, a light color will show a small number of iterations and a dark color a large number of iterations.

Based on polynomiographs, we can compute some numerical measures that will allow the analysis to be broadened. The three most widely used measures are: The average number of iterations (ANI) ^[27], convergence area index (CAI) ^[27], and the computations time ^[28]. To compute the average number of iterations, we take the information on the speed of convergence in the polynomiograph because it contains the number of iterations needed to find the root for all points in the considered area $A$. CAI is the ratio of the number of starting points that converged to any root to the number of all points in the considered area $A$. The CAI value is a real number between $0$ and $1$, and it gives us information on the percentage of points in $A$ that have converged to roots.

The polynomiographs used for the analysis presented in this section were generated for the following common parameters: $N = 30$, $\varepsilon = 0.001$, and polynomiograph resolution $1000 \times 1000$ pixels. In our analysis, we use three polynomials:

● $p_3(z) = z^3 - 1$, roots: $1$, $-0.5 + \frac{\sqrt{3}}{2}i$, $-0.5 - \frac{\sqrt{3}}{2}i$, the area $A = [-2, 2]^2$,

● $p_4(z) = z^4 -10 z^2 + 9$, roots: $-3$, $-1$, $1$, $3$, the area $A = [-4, 4]^2$,

● $p_5(z) = z^5 - z$, roots: $0$, $-1$, $1$, $-i$, $i$, the area $A = [-2, 2]^2$.

The program for generating polynomiographs was implemented in Mathematica 13.2 using the parallelization option of the Compile command. Moreover, the polynomiographs were generated on the computer with the specifications: Intel i5-9600K (@3.70 GHz), 32 GB DDR4 RAM, and Windows 10 (64 bit). Apart from the two proposed algorithms, PCNM8 and PCNM4, we generated polynomiographs for the PJNM, KTNM, CNM1, CNM2, and ONM algorithms for comparison purposes.

The obtained results for the $p_3$ polynomial are presented in Figure 2 (the polynomiographs) and in Table 1 (the numerical measures). From the polynomiographs, we see that each algorithm has different basins of attraction. However, in each case, we see three characteristic braids. The degree of interweaving of the basins tells us about the stability of the algorithm (the less the interweaving, the better the stability of the algorithm). Based on the polynomiographs, we see that the best stability is obtained by the PJNM algorithm, followed by the PCNM8 and ONM algorithms. For the CNM1 algorithm, we also see a low degree of interweaving, but there are also non-convergent points (the black areas) which are not visible in the three algorithms mentioned. This is confirmed by the CAI measure presented in Table 1. The PJNM, PCNM8 and ONM algorithms all obtained values equal to $1.0$, whereas for the CNM1 algorithm, the CAI value is less than $1.0$, namely it is equal to $0.967$. The lowest value of CAI was obtained by the CNM2 algorithm ($0.946$). When it comes to the speed of convergence, we see that the lightest colors are visible in the polynomiographs for the PCNM8 and ONM algorithms, followed by the PCNM4 and CNM1 algorithms. Again, we can confirm this observation by looking at the values of ANI gathered in Table 1. The lowest value of ANI was obtained by the PCNM8 algorithm ($1.496$), and the second lowest value equal to $1.624$ was obtained by the ONM algorithm. Now, if we look at the generation time of polynomiographs, we see that the fastest algorithm is the CNM1 algorithm ($0.762$ s). The two best algorithms in terms of ANI, i.e., the PCNM8 and ONM, obtained slightly slower times equal to $0.778$ and $0.775$ s, respectively. The slowest algorithm is the KTNM algorithm ($0.914$ s).

In the next example, we present and analyze the results obtained for the $p_4$ polynomial. In Figure 3, we see the polynomiographs for this polynomial, whereas in Table 2, we gather the values of numerical measures. In the polynomiographs, we see four basins of attraction. The most significant difference in those basins is visible at the boundaries, where the interweaving of the basins appears. The lowest degree of interweaving can be observed for the PCNM4, PJNM, and CNM1 algorithms. Thus, these three algorithms obtained the best stability among the analyzed algorithms. However, for the CNM1 algorithm, we can also observe some non-convergent points in the polynomiograph, which causes the CAI value to be equal to $0.999$. The other two algorithms obtained a CAI value of $1.0$. There are also two other algorithms with full convergence (CAI value equal to $1.0$), namely the PCNM8 and ONM algorithms. In terms of speed of convergence, the best results are obtained by the ONM and PCNM8 algorithms, for which the ANI value is equal to $1.566$ and $1.660$, respectively. We can observe this in the polynomiographs because these two algorithms have the biggest light areas. On the other extreme, the worst speed of convergence can be observed for the PJNM algorithm (the darkest colors and the value of ANI equal to $7.561$). If we look at the generation times, which reflect the real time of computations, we see that the fastest algorithm was the PCNM8 algorithm, followed by the CNM1 algorithm. These two algorithms obtained the times equal to $0.781$ and $0.840$ s, respectively. The slowest algorithm among the analyzed algorithms was the PJNM algorithm, which obtained a time equal to $1.085$ s.

In the last example, we visually analyze the results obtained for the $p_5$ polynomial. The polynomiographs generated for this polynomial are presented in Figure 4, and the calculated numerical measures are in Table 3. The best stability is observed for the PJNM algorithm, followed by the CNM1 algorithm. For these two algorithms, the interweaving of the basins is the smallest. However, for the CNM1 algorithm, we can observe areas of non-convergent points, whereas, for the PJNM algorithm, we do not see any such points. The values of CAI gathered in Table 3 confirm this observation, where the PJNM and CNM1 algorithms obtained values of $1.0$ and $0.997$, respectively. The highest value of CAI ($1.0$) is also obtained by three other analyzed algorithms, namely by the PCNM8, PCNM4, and ONM algorithms. The lowest value of $0.982$ is obtained by the KTNM algorithm, which is clearly visible in the polynomiograph, where we see many black areas. Now, if we analyze the speed of convergence, then we notice that the fastest algorithms are the ONM and PCNM8 algorithms, for which we see the lightest shades of the basins' colors. Consequently, these two algorithms obtained the lowest values of ANI, equal to $1.556$ and $2.013$, respectively. The worst speed of convergence is visible for the PJNM algorithm. Moreover, this algorithm also obtained the highest ANI value, which was equal to $5.678$. When it comes to the generation times, like for the $p_3$ polynomial, the shortest time of $0.834$ s was obtained by the CNM1 algorithm. The second best algorithm with a time equal to $0.838$ s is the PCNM8 algorithm.

5.   Numerical results

In this section, we investigate the effectiveness of newly proposed fourth- and eighth-order algorithms for solving several nonlinear real-world application problems including single and multi-variables. Several parameters have been specified to compare the proposed algorithms with some of the existing ones having similar order of convergence. The parameters include the number of iterations (NI), the number of function evaluations (NFE), the absolute error at the final iteration ($|\eta| = |x_N - x_{N-1}|$), the absolute value of the nonlinear function ($|\psi(x_N)|$), and the CPU time (in seconds). As the stopping criterion, we use the following tolerance:

Problem 1. For the first nonlinear model, we consider a path traversed by an electron in the air gap between two parallel plates considering the multi-factor effect, which is given by ^[29]

where $\mu_0$ and $v_0$ are the position and velocity of the electron at time $x$, $m$ and $c_0$ are the mass and the charge of the electron at rest, respectively, and $E_0\sin(\omega x + \sigma)$ is the radio frequency electric field between the plates.

If particular parameters are chosen, Eq (5.1) can be simplified as

The results of numerical simulations of an engineering application Problem 1 carried out for a fixed number of iterations, that is, when NI $= 7$ and $x_0 = 10.5$ are presented in Table 4. It is observed that the proposed fourth-order PCNM4 algorithm produces, at the final iteration, the smallest absolute error and the smallest absolute value of the function itself when compared to the other existing algorithms. This proposed fourth-order algorithm has the advantage of less CPU time consumption in comparison to other fourth-order existing algorithms. Similarly, in Table 5, it is observed that the existing algorithm of eight-order convergence (ONM) takes the fewest number of iterations (NI) in comparison to the other algorithms, while the proposed eight-order PCNM8 algorithm produces the smallest absolute error when compared to the other algorithms. In Table 5, we observe that the existing algorithms of order four PJNM, KTNM, CNM1, and CNM2 produce the smallest absolute errors when compared to the proposed fourth-order PCNM4, while the PCNM4 algorithm takes a fewer number of iterations in comparison to the other fourth-order existing algorithms.

Problem 2. NASA's "Wind" Satellite's Orbit. NASA intends to launch a spacecraft named Wind that will remain in a fixed position along a line from the earth to the sun so that the solar wind will travel by the satellite on its approach to the earth ^[30]. We devise the following equation based on the relevant physical laws:

where, $G = 6.67 \times 10^{-11}$ [Nm$^2$/kg$^2$], $M_s = 1.98 \times 10^{30}$[kg], $M_e = 5.98 \times 10^{24}$[kg], $m =$ the mass of satellite [kg], $R = 1.49 \times 10^{11}$[m], $T = 3.15576 \times 10^{7}$[s] and $\omega = \frac{2\pi}{T}$. The required solution to the non-linear real-world problem $\psi_2(r) = 0$ is $r^* = 1.4761775009615046593 \times 10^{11}$.

The numerical simulations of the real-life application problem carried out in Table 6 take a fixed number of iterations, meaning NI $= 11$ and $x_0 = 7.6$. It is observed that the existing fourth-order PJNM algorithm produces, at the final iteration, the smallest absolute error and the smallest absolute value of the function itself when compared to the other algorithms. However, the proposed fourth-order algorithm PCNM4 has the advantage of less CPU time consumption in comparison to other fourth-order existing algorithms. Similarly, by fixing an initial guess of $x_0 = 5.5$, the numerical results are presented in Table 7. It is observed that the proposed PCNM8 and existing algorithm ONM (both having eighth-order convergence) take the fewest number of iterations in comparison to the other algorithms, whereas the existing eight-order ONM algorithm produces the smallest absolute error when compared to the other algorithms. We further observe that the PCNM4 and PJNM fourth-order algorithms takes the fewest number iteration (NI) in comparison to the other fourth-order algorithms KJNM, CNM1, and CNM2.

Problem 3. The system of nonlinear equations from ^[31] is:

The numerical simulations for solving the system numerical result carried out in Table 8 take a fixed number of iterations, that is $NI = 9$, and initial guess for two by two non linear system are $x_1 = 0.1$ and $x_2 = 0.2$. It is observed that the proposed algorithm PCNM4 produces the smallest absolute error when compared to the other algorithms. The proposed fourth-order algorithm PCNM4 has the advantage of less CPU time consumption in comparison to other fourth-order existing algorithms.

Problem 4. The system of three nonlinear equations is ^[32]:

The numerical simulations for solving the system (5.5) carried out in Table 9 take a fixed number of iterations, that is $NI = 6$, and the initial guesses for the three-by-three nonlinear system are chosen to be $x_1 = 0.8$, $x_2 = 1.0$, $x_3 = 0.8$. Although all the algorithms converge, we observe that the proposed algorithm PCNM4 produces the smallest absolute error when compared to the other algorithms. The proposed fourth-order algorithm PCNM4 has the advantage of less CPU time consumption in comparison to other fourth-order existing algorithms.

Problem 5. Neurophysiology application ^[33]:

The nonlinear model consists of the following six equations:

The constants $c_i$ in (5.6) can be randomly chosen. In our experiment, we considered $c_i = 0$ for $i = 1, \dots, 4$.

The numerical simulations for solving the real-life application system used in medical science (Neurophysiology) are carried out in Table 10. Neurophysiology is the branch of physiology that focuses on the study of the nervous system, including the brain, spinal cord, and peripheral nerves. The numerical simulations carried out in Table 10 take a fixed number of iterations, that is $NI = 7$, and the initial guesses for the system given in (5.6) are chosen to be $x_1 = 0.9$, $x_2 = 0.8$, $x_3 = 0.7$, $x_4 = 0.5$, $x_5 = 0.3$, $x_6 = 0.1$. It is observed that the proposed algorithm, PCNM4, produces the smallest absolute error when compared to the other algorithms. The existing fourth-order algorithms taken for comparison, namely; PJNM, CNM1, and CNM2 diverge from the exact solution, while the algorithm KTNM converges but is comparatively slower than the proposed one. Hence, the performance of the proposed fourth-order algorithm for the nonlinear model (5.6) is far better than some of the existing approaches.

6.   Conclusions and future remarks

We introduce two predictor-corrector algorithms, amalgamating Taylor's series and the composition approach. One of these algorithms demonstrates eighth-order convergence alongside a high-efficiency index of approximately 1.5157, surpassing certain established techniques. The other algorithm achieves fourth-order convergence. Our comprehensive convergence analyses encompass both local and semilocal aspects. We use a number of complex polynomials and polynomiographs that show how much faster the modified algorithms converge, which is especially clear when we compare basins of attraction to other algorithms. Real-world applications in chemistry, astronomy, and neurology validate our simulations, where the proposed algorithms consistently outperform their numerical counterparts. It is intended to return to the proposed approach in the future to make any necessary adjustments to attain the best potential algorithm.

Author contributions

D. K. Almutairi: Validation, writing – review and editing; I. K. Argyros: Methodology, formal analysis, writing – review and editing; K. Gdawiec: Formal analysis, software, visualization, writing – original draft, writing – review and editing; S. Qureshi: Conceptualization, software, writing – original draft; A. Soomro: Software, investigation, writing – original draft; K. H. Jamali: Supervision, validation; M. Alquran: Supervision, resources, data curation; A. Tassaddiq: Conceptualization, methodology, validation, formal analysis, investigation, resources, writing – review and editing, project administration, funding acquisition. All authors have read and approved the final version of the manuscript for publication.

Acknowledgments

The author extends the appreciation to the Deanship of Postgraduate Studies and Scientific Research at Majmaah University for funding this research work through the project number (ER-2024-1216).

Conflict of interest

The authors declare that they have no conflict of interest.