A New High-Order Fractional Parallel Iterative Scheme for Solving Nonlinear Equations

Abstract

Solving fractional-order nonlinear equations is crucial in engineering, where precision and accuracy are essential. This study introduces a novel fractional parallel technique for solving nonlinear equations. To enhance convergence, we incorporate a simple root-finding method of order 3γ + 1 as a correction term in the parallel scheme. Theoretical analysis shows that the parallel scheme achieves a convergence order of 6γ + 3. Using a dynamical system approach, we identify optimal parameter values, and the symmetry in the dynamical planes for different fractional parameters demonstrates the method’s stability and consistency in handling nonlinear problems. These parameter values are applied to the parallel scheme, yielding highly consistent results. Several engineering problems are examined to assess the method’s efficiency, stability, and consistency compared to existing methods.

1. Introduction

In 1695, Leibniz and L’Hôpital introduced the concept of the semi-derivative [1,2,3], which laid the foundation for both classical and fractional calculus. Fractional calculus [4], an extension of classical calculus to non-integer orders, has gained attention for its effectiveness in modeling complex systems with non-local and memory-dependent behaviors. This makes it particularly valuable in engineering, where fractional-order nonlinear equations are used to simulate systems with memory effects and anomalous diffusion.

One form of the nonlinear problems we are particularly interested in is represented by the following equation:

$$\left\{\begin{array}{c}\left[{}_{\check{c}}⅁_{{\gamma }_1}^{\gamma }\right]y\left(x\right)+\hslash \left(y\left(x\right)\right)=g\left(x\right),\\ y\left(x_0\right)=y_0,y^{\prime }\left(x_0\right)=y_1,{\theta }_0^{[\ast ]}<\gamma \le {\theta }_1^{[\ast ]},\end{array}\right.$$

(1)

In Equation (1), $x_0$ is the initial point and fractional parameter $\gamma$ lies within the interval $\left({\theta }_0^{[\ast ]},{\theta }_1^{[\ast ]}\right]$. While obtaining exact solutions to these equations would be computationally attractive due to their precision and lack of approximation, the inherent complexities of fractional-order systems often make such solutions unattainable. In these cases, analytical methods, such as series expansions and integral transformations, play a critical role. These methods provide closed-form or semi-closed-form solutions, which not only offer valuable insights into the system’s behavior, but also serve as reference points for validating and guiding numerical solutions. By doing so, they help ensure that numerical approaches remain accurate and reliable when exact solutions are not feasible.

The Caputo derivative retains the benefits of classical calculus while enabling the development of more detailed models that can accurately capture the complex dynamics observed in various real-world applications, including finance [5], engineering [6,7], physics [8], and many more. Due to their inherent complexity, nonlinear equations remain one of the most challenging problems in science and engineering. Unlike the Caputo derivative, most other fractional derivatives do not satisfy the condition $\left[{⅁}^{\gamma }\right]\left(1\right)=0$ when $\gamma$ is not a natural number. Because of this, we will explore the fundamental concepts of fractional calculus and discuss the fractional iterative approach for solving nonlinear equations using Caputo-type derivatives.

Definition 1.

The Gamma function, also known as the generalized factorial function [9], is defined as:

$$\lceil \left(x\right)=\underset{0}{\overset{+\infty }{\int }}{u}^{x-1}{e}^{-u}du,$$

(2)

where $x>0$. The Gamma function satisfies the properties $\lceil \left(1\right)=1$ and $\lceil (n+1)=n!$ for any positive integer n.

Definition 2.

Let $\hslash :\mathbb{R}\to \mathbb{R},$ where

$$\hslash \in {\mathbb{C}}^{+\infty }\left(\left[{\gamma }_1,x\right]\right),-\infty <{\gamma }_1<x<+\infty ,{\gamma }_1\ge 0,m=\left[{\gamma }_1\right]+1.$$

(3)

The Caputo fractional derivative [10] of order ${\gamma }_1$ is defined as:

$$\left[{}_{\check{c}}⅁_{{\gamma }_1}^{\gamma }\right]\hslash \left(x\right)=\left\{\begin{array}{c}{\displaystyle \frac{1}{\lceil (m-\gamma )}}\underset{{\gamma }_1}{\overset{x}{\int }}\left[{\displaystyle \frac{d^m}{d{t}^m}}\hslash \left(t\right){\displaystyle \frac{1}{{\left(x-t\right)}^{\gamma -m+1}}}\right]dt,\gamma \notin N,\\ {\displaystyle \frac{d^{m-1}}{d{t}^{m-1}}}\hslash \left(x\right),\gamma =m-1\in N\cup \left\{0\right\},\end{array}\right.$$

(4)

where $\lceil \left(x\right)$ is a gamma function for $x>0.$

Theorem 1.

Suppose that $\left[{}_{\check{c}}⅁_{{\gamma }_1}^{m\gamma }\right]\hslash \left(x\right)\in C\left[{\gamma }_1,{\gamma }_2\right]$ for $m=1,\dots ,n+1$ where $\gamma \in \left(0,1\right]$. Then, the Generalized Taylor Formula [11] is given by:

$$\hslash \left(x\right)=\left[\sum _{i=0}^n\left[{}_{\check{c}}⅁_{{\gamma }_1}^{i\gamma }\right]\hslash \left({\gamma }_1\right){\displaystyle \frac{{\left(x-{\gamma }_1\right)}^{i\gamma }}{\lceil (i\gamma +1)}}+\left[{}_{\check{c}}⅁_{{\gamma }_1}^{\left(n+1\right)\gamma }\right]\hslash \left(\xi \right){\displaystyle \frac{{\left(x-{\gamma }_1\right)}^{\left(n+1\right)\gamma }}{\lceil (\left(n+1\right)\gamma +1)}}\right],$$

(5)

where

$${\gamma }_1\le \xi \le x,\forall x\in \left({\gamma }_1,{\gamma }_2\right],$$

(6)

and

$$\left[{}_{\check{c}}⅁_{{\gamma }_1}^{n\gamma }\right]=\left[\left[{}_{\check{c}}⅁_{{\gamma }_1}^{\gamma }\right].\left[{}_{\check{c}}⅁_{{\gamma }_1}^{\gamma }\right]\dots \left[{}_{\check{c}}⅁_{{\gamma }_1}^{\gamma }\right](n-times)\right].$$

(7)

Consider the Caputo-type Taylor expansion of $\hslash \left(x\right)$ around ${\gamma }_1=\xi$:

$$\hslash \left(x\right)=\left[{\displaystyle \frac{\left[{}_{\check{c}}⅁_{\xi }^{1\gamma }\right]\hslash \left(\xi \right)}{\lceil (+1)}}{\left(x-\xi \right)}^{\gamma }+{\displaystyle \frac{\left[{}_{\check{c}}⅁_{\xi }^{2\gamma }\right]\hslash \left(\xi \right)}{\lceil (2\gamma +1)}}{\left(x-\xi \right)}^{2\gamma }+O{\left(x-\xi \right)}^{3\gamma }\right].$$

(8)

By factoring out ${\displaystyle \frac{\left[{}_{\check{c}}⅁_{\xi }^{1\gamma }\right]\hslash \left(\xi \right)}{\lceil (\gamma +1)}}$, we obtain:

$$\hslash \left(x\right)=\left[{\displaystyle \frac{\left[{}_{\check{c}}⅁_{\xi }^{1\gamma }\right]\hslash \left(\xi \right)}{\lceil (\gamma +1)}}\left[\begin{array}{c}{\left(x-\xi \right)}^{\gamma }+\\ {\check{c}}_2{\left(x-\xi \right)}^{2\gamma }\end{array}\right]+O{\left(x-\xi \right)}^{3\gamma },\right]$$

(9)

where

$${\check{c}}_j={\displaystyle \frac{\lceil (\gamma +1)}{\lceil (m\gamma +1)}}{\displaystyle \frac{\left[{}_{\check{c}}⅁_{\xi }^{m\gamma }\right]\hslash \left(\xi \right)}{\left[{}_{\check{c}}⅁_{\xi }^{\gamma }\right]\hslash \left(\xi \right)}},m=2,3,\dots$$

(10)

The corresponding Caputo-type derivative of $\hslash \left(x\right)$ around ξ is:

$$\left[{}_{\check{c}}⅁_{\xi }^{\gamma }\right]\hslash \left(x\right)=\left[{\displaystyle \frac{\left[{}_{\check{c}}⅁_{\xi }^{1\gamma }\right]\hslash \left(\xi \right)}{\lceil (\gamma +1)}}\left[\begin{array}{c}\lceil (\gamma +1)+\\ {\displaystyle \frac{\lceil (2\gamma +1)}{\lceil (\gamma +1)}}{\check{c}}_2{\left(x-\xi \right)}^{\gamma }\end{array}\right]+O{\left(x-\xi \right)}^{2\gamma }\right].$$

(11)

These expansions are used in the convergence analysis of the proposed method.

The rest of the paper is structured as follows: Section 2 discusses the development, convergence, and complex dynamical analysis of fractional-order schemes for solving Equation (1). Section 3 evaluates the effectiveness and stability of the proposed method through numerical examples and comparisons with other techniques. Finally, the conclusions are presented in Section 5.

2. Construction and Dynamic Analysis of Fractional Iterative Scheme

The Newton method [12,13] refines an initial guess for solving a nonlinear problem by using a sequence of linear approximations. This technique, known for its rapid convergence—often quadratic—has been widely used since its development in the late 14th century to solve nonlinear equations. The fractional Newton method is an extension of the traditional Newton method, designed for nonlinear equations involving fractional derivatives. It incorporates fractional calculus operators, such as fractional integrals and derivatives, into the iterative process. This approach has demonstrated strong convergence properties and has been effective in solving complex nonlinear equations where other methods may fail. Several studies, including those by Torres-Hernandez et al. [14], Akgül et al. [15], Gajori et al. [16], Kumar et al. [17], and Candelario et al. [18] describe a fractional version of the Newton technique with various fractional derivatives.

Based on both stability and efficiency, the method presented in [19] is preferable to other well-known root-finding methods for nonlinear equations, such as those in [20,21,22,23] and the references cited therein. Nonlinear equations, particularly those involving fractional approaches, have gained significant attention due to their wide applications in addressing complex engineering problems. The numerical algorithm discussed exhibits local convergence in approximating simple solutions to fractional-order nonlinear problems. In addition to the basic iterative strategies for solving (1) with integer or fractional convergence orders, there are parallel schemes with global convergence that can find all solutions to a nonlinear equation simultaneously. These schemes are more stable, consistent, and are particularly useful in parallel computing (see, e.g., [24,25,26] and the references cited therein). However, a comprehensive literature review reveals that little work has focused on fractional-order numerical techniques for solving fractional nonlinear problems.

The primary goal of the current research is to enhance the convergence rate, stability, and accuracy of these parallel schemes. Challenges remain, particularly in ensuring efficient computation in higher dimensions. A comprehensive review of these trends, along with their limitations—such as handling fractional-order nonlinearities—provides context for the new numerical scheme. The proposed fractional family of parallel numerical schemes addresses several limitations of traditional iterative methods in the following ways:

• Parallel computation: Unlike simple root-finding methods, parallel approaches compute all solutions simultaneously.
• Global convergence: Parallel schemes exhibit global convergence, as all roots contribute to solving nonlinear problems.
• Robustness: Parallel schemes outperform other methods in terms of robustness and sensitivity to initial values for exact solutions.
• Efficiency in parallel computing: Fractional-order parallel methods are highly suitable for multicore and distributed systems, as updates for each solution can be computed independently. In contrast to simpler schemes, which often require shared access to derivatives and global convergence checks, parallel updates can be handled with minimal communication overhead.
• Accuracy and consistency: Compared to simple iterative methods for solving (1), fractional parallel approaches—relying on differences between approximated roots—tend to be more accurate and consistent, requiring fewer modifications to maintain global convergence.
• Memory efficiency: Fractional-order parallel methods are more memory-efficient in practical implementations compared to basic schemes.

On the other hand, when applied to nonlinear fractional-order differential equations with movable conditions, the method may exhibit slower convergence rates and higher computational demands. These challenges can introduce complex dynamics that might affect the overall effectiveness of the approach. Therefore, while the method is generally reliable, a careful analysis of the problem’s specific nature and conditions is essential to optimize its performance. In most cases, the method performs well, but in certain instances, additional strategies may be necessary to improve convergence and reduce computational costs.

Building on the method described in [19], we consider the following fractional iterative scheme:

$$v^{\left[l\right]}=x^{\left[l\right]}-{\displaystyle \frac{\hslash \left(x^{\left[l\right]}\right)}{{\hslash }^{\prime }\left(x^{\left[l\right]}\right)}}\left(\left({\displaystyle \frac{\hslash \left({\varsigma }^{\left[l\right]}\right)}{\hslash \left(x^{\left[l\right]}\right)}}\right)+2{\left({\displaystyle \frac{\hslash \left({\varsigma }^{\left[l\right]}\right)}{\hslash \left(x^{\left[l\right]}\right)}}\right)}^2+{\beta }^{[\ast ]}{\left({\displaystyle \frac{\hslash \left({\varsigma }^{\left[l\right]}\right)}{\hslash \left(x^{\left[l\right]}\right)}}\right)}^4\right),$$

(12)

where ${\varsigma }^{\left[l\right]}=x^{\left[l\right]}-{\displaystyle \frac{\hslash \left(x^{\left[l\right]}\right)}{{\hslash }^{\prime }\left(x^{\left[l\right]}\right)}}.$ This method has a convergence order of 4 (abbreviated as M^[γ]), which satisfies the following error equation:

$$e_i^{\left[\ast \right]}=\left({\beta }^{[\ast ]}{\check{c}}_2^3+5{\check{c}}_2^3-{\check{c}}_2{\check{c}}_3\right){\left(e^{\left[l\right]}\right)}^4+O\left({\left(e^{\left[l\right]}\right)}^5\right),$$

(13)

where $e_i^{\left[\ast \right]}=x^{\left[1\right]}-\xi$, $e^{\left[l\right]}=x^{\left[l\right]}-\xi$, and ${\check{c}}_j={\displaystyle \frac{1}{m!}}{\displaystyle \frac{\left[{}_{\check{c}}⅁_{\xi }^m\right]\hslash \left(\xi \right)}{\left[{}_{\check{c}}⅁_{\xi }^{\gamma }\right]\hslash \left(\xi \right)}},m\ge 2,3,\dots$ The fractional version of this method is given by:

$$v^{\left[l\right]}=x^{\left[l\right]}-{\left(\lceil (\gamma +1){\displaystyle \frac{\hslash \left(x^{\left[l\right]}\right)}{{\hslash }^{\prime }\left(x^{\left[l\right]}\right)}}\left(\begin{array}{c}\left({\displaystyle \frac{\hslash \left({\varsigma }^{\left[l\right]}\right)}{\hslash \left(x^{\left[l\right]}\right)}}\right)+2{\left({\displaystyle \frac{\hslash \left({\varsigma }^{\left[l\right]}\right)}{\hslash \left(x^{\left[l\right]}\right)}}\right)}^2\\ +{\beta }^{[\ast ]}{\left({\displaystyle \frac{\hslash \left({\varsigma }^{\left[l\right]}\right)}{\hslash \left(x^{\left[l\right]}\right)}}\right)}^4\end{array}\right)\right)}^{1/\gamma },$$

(14)

where ${\varsigma }^{\left[l\right]}=x^{\left[l\right]}-{\left(\lceil (\gamma +1){\displaystyle \frac{\hslash \left(x^{\left[l\right]}\right)}{\left[{}_{\check{c}}⅁_{{\gamma }_1}^{\gamma }\right]\hslash \left(x^{\left[l\right]}\right)}}\right)}^{1/\gamma }$. We abbreviate this method as MM^[γ].

2.1. Convergence Analysis

For the iterative scheme given by Equation (14), we establish the following theorem to determine its order of convergence.

Theorem 2.

Let

$$\hslash :⅁\subseteq \mathbb{R}\to \mathbb{R}$$

(15)

be a continuous function, with $\left[{}_{\check{c}}⅁_{{\gamma }_1}^{m\gamma }\right]\hslash \left(x\right)$ of order m for any $\gamma \ge 0$ and $\gamma \in \left(0,1\right]$ containing the exact root ξ of $\hslash \left(x\right)$. Furthermore, given a sufficiently close starting value $x_0$, the convergence order of the Caputo-type fractional iterative scheme

$$v^{\left[l\right]}=x^{\left[l\right]}-{\left(\lceil (\gamma +1){\displaystyle \frac{\hslash \left(x^{\left[l\right]}\right)}{{\hslash }^{\prime }\left(x^{\left[l\right]}\right)}}\left(\begin{array}{c}\left({\displaystyle \frac{\hslash \left({\varsigma }^{\left[l\right]}\right)}{\hslash \left(x^{\left[l\right]}\right)}}\right)\\ +2{\left({\displaystyle \frac{\hslash \left({\varsigma }^{\left[l\right]}\right)}{\hslash \left(x^{\left[l\right]}\right)}}\right)}^2\\ +{\beta }^{[\ast ]}{\left({\displaystyle \frac{\hslash \left({\varsigma }^{\left[l\right]}\right)}{\hslash \left(x^{\left[l\right]}\right)}}\right)}^4\end{array}\right)\right)}^{1/\gamma },$$

(16)

is least $2\gamma +1$, with the error equation given by:

$$e^{\left[\ast \right]}=e^{\left[\ast \right]}=\left(\begin{array}{c}{\pounds }_1^{[\ast ]}+{\pounds }_2^{[\ast ]}\\ +{\pounds }_3^{[\ast ]}+{\pounds }_4^{[\ast ]}\\ +{\pounds }_5^{[\ast ]}+{\pounds }_6^{[\ast ]}\end{array}\right){\left[e^{\left[l\right]}\right]}^{3\gamma +1}+O\left({\left[e^{\left[l\right]}\right]}^{4\gamma +1}\right),$$

(17)

where the constants ${\pounds }_i^{[\ast ]}$ are defined as follows:

${\pounds }_1^{[\ast ]}=-9{\displaystyle \frac{9{c_2}^3{\left(2^{\gamma }\right)}^4{\left(\lceil \left(\gamma +1/2\right)\right)}^2}{{\left(\lceil \left(\gamma \right)\right)}^2{\gamma }^2\pi }}$, ${\pounds }_2^{[\ast ]}=+{\displaystyle \frac{4{c_2}^3{\left(2^{\gamma }\right)}^6{\left(\lceil \left(\gamma +1/2\right)\right)}^3}{{\left(\lceil \left(\gamma \right)\right)}^3{\gamma }^3{\pi }^{3/2}}}$,

${\pounds }_3^{[\ast ]}=-{\displaystyle \frac{c_3{\left(3^{\gamma }\right)}^3\sqrt{3}\lceil \left(\gamma +1/3\right)\lceil \left(\gamma +2/3\right)c_2}{{\left(\lceil \left(\gamma \right)\right)}^2{\gamma }^2\pi }}$, ${\pounds }_4^{[\ast ]}={\displaystyle \frac{3}{2}}{\displaystyle \frac{c_3{\left(3^{\gamma }\right)}^3\sqrt{3}\lceil \left(\gamma +1/3\right)\lceil \left(\gamma +2/3\right)c_2}{\lceil \left(\gamma \right)\gamma \sqrt{\pi }{\left(2^{\gamma }\right)}^2\lceil \left(\gamma +1/2\right)}}$,

${\pounds }_5^{[\ast ]}=-4c_2{c}_3+{c_2}^3$, ${\pounds }_6^{[\ast ]}={\displaystyle \frac{3{\left(2^{\gamma }\right)}^2\lceil \left(\gamma +1/2\right)c_2{c}_3}{\lceil \left(\gamma \right)\gamma \sqrt{\pi }}}+{\displaystyle \frac{4{c_2}^3{\left(2^{\gamma }\right)}^2\lceil \left(\gamma +1/2\right)}{\lceil \left(\gamma \right)\gamma \sqrt{\pi }}}.$

Here, ${\check{c}}_j={\displaystyle \frac{\lceil (\gamma +1)}{\lceil (m\gamma +1)}}{\displaystyle \frac{\left[{}_{\check{c}}⅁_{\xi }^{m\gamma }\right]\hslash \left(\xi \right)}{\left[{}_{\check{c}}⅁_{\xi }^{\gamma }\right]\hslash \left(\xi \right)}},m\ge 2,3,\dots$ and ${\left[\lceil \left(.\right)\right]}^n{=\lceil }^n(.).$

Proof. 

Let $\xi$ be a root of ℏ and assume $x^{\left[l\right]}=\xi +e^{\left[l\right]}$. Using a Taylor series expansion of $\hslash \left(x^{\left[l\right]}\right)$ and $\left[{}_{\check{c}}⅁_{{\gamma }_1}^{\gamma }\right]\hslash \left(x^{\left[l\right]}\right)$ around $x=\xi ,$ and considering that $\hslash \left(\xi \right)=0,$ we obtain:

$$\hslash \left(x^{\left[l\right]}\right)={\displaystyle \frac{\left[{}_{\check{c}}⅁_{\xi }^{1\gamma }\right]\hslash \left(\xi \right)}{\lceil (\gamma +1)}}\left[{\left[e^{\left[l\right]}\right]}^{\gamma }+{\check{c}}_2{\left[e^{\left[l\right]}\right]}^{2\gamma }+{\check{c}}_3{\left[e^{\left[l\right]}\right]}^{3\gamma }\right]+O\left({\left[e^{\left[l\right]}\right]}^{4\gamma }\right),$$

(18)

$$\left[{}_{\check{c}}⅁_{{\gamma }_1}^{\gamma }\right]\hslash \left(x^{\left[l\right]}\right)={\displaystyle \frac{\left[{}_{\check{c}}⅁_{\xi }^{1\gamma }\right]\hslash \left(\xi \right)}{\lceil (\gamma +1)}}\left[\begin{array}{c}\lceil (\gamma +1)+{\displaystyle \frac{\lceil (2\gamma +1)}{\lceil (\gamma +1)}}{\check{c}}_2{\left[e^{\left[l\right]}\right]}^{\gamma }\\ +{\displaystyle \frac{\lceil (3\gamma +1)}{\lceil (2\gamma +1)}}{\check{c}}_3{\left[e^{\left[l\right]}\right]}^{2\gamma }\end{array}\right]+O\left({\left[e^{\left[l\right]}\right]}^{3\gamma }\right),$$

(19)

$${\displaystyle \frac{1}{\left[{}_{\check{c}}⅁_{{\gamma }_1}^{\gamma }\right]\hslash \left(x^{\left[l\right]}\right)}}={\displaystyle \frac{1}{\lceil \left(\gamma +1\right)}}-{\displaystyle \frac{\lceil \left(2\gamma +1\right)c_2{\left[e^{\left[l\right]}\right]}^{\gamma }}{{\left(\lceil \left(\gamma +1\right)\right)}^3}}+{\Theta }_1^{[\ast ]}+{\Theta }_2^{[\ast ]}+{\Theta }_3^{[\ast ]}+\dots$$

(20)

where ${\Theta }_1^{[\ast ]}={\displaystyle \frac{{\left[e^{\left[l\right]}\right]}^{2\gamma }}{\lceil \left(\gamma +1\right)}}\left(\begin{array}{c}-{\displaystyle \frac{\lceil \left(3\gamma +1\right)c_3}{\lceil \left(\gamma +1\right)\lceil \left(2\gamma +1\right)}}\\ +{\displaystyle \frac{{\left(\lceil \left(2\gamma +1\right)\right)}^2{c_2}^2}{{\left(\lceil \left(\gamma +1\right)\right)}^4}}\end{array}\right)$,

${\Theta }_2^{[\ast ]}={\displaystyle \frac{{\left[e^{\left[l\right]}\right]}^{3\gamma }}{\lceil \left(\gamma +1\right)}}\left({\displaystyle \frac{c_2\lceil \left(3\gamma +1\right)c_3}{{\left(\lceil \left(\gamma +1\right)\right)}^3}}-{\displaystyle \frac{\left(\begin{array}{c}{\left(\lceil \left(2\gamma +1\right)\right)}^3{c_2}^2-\\ \lceil \left(3\gamma +1\right)c_3{\left(\lceil \left(\gamma +1\right)\right)}^3\end{array}\right)c_2}{{\left(\lceil \left(\gamma +1\right)\right)}^6}}\right),$

${\Theta }_3^{[\ast ]}={\displaystyle \frac{{\left[e^{\left[l\right]}\right]}^{3\gamma }}{\lceil \left(\gamma +1\right)}}\left({\displaystyle \frac{c_2\lceil \left(3\gamma +1\right)c_3}{{\left(\lceil \left(\gamma +1\right)\right)}^3}}-{\displaystyle \frac{\left(\begin{array}{c}{\left(\lceil \left(2\gamma +1\right)\right)}^3{c_2}^2-\\ \lceil \left(3\gamma +1\right)c_3{\left(\lceil \left(\gamma +1\right)\right)}^3\end{array}\right)c_2}{{\left(\lceil \left(\gamma +1\right)\right)}^6}}\right).$

Thus, by multiplying Equations (18) and (20), we obtain

$${\displaystyle \frac{\hslash \left(x^{\left[l\right]}\right)}{\left[{}_{\check{c}}⅁_{{\gamma }_1}^{\gamma }\right]\hslash \left(x^{\left[l\right]}\right)}}=\left[\begin{array}{c}{\displaystyle \frac{{\left[e^{\left[l\right]}\right]}^{\gamma }}{\lceil \left(\gamma \right)\gamma }}+\left({\Theta }_4^{[\ast ]}\right){\left[e^{\left[l\right]}\right]}^{2\gamma }\\ +\left({\Theta }_5^{[\ast ]}\right){\left[e^{\left[l\right]}\right]}^{3\gamma }+\dots \end{array}\right],$$

(21)

where ${\Theta }_4^{[\ast ]}={\displaystyle \frac{c_2}{\lceil \left(\gamma \right)\gamma }}-{\displaystyle \frac{{\left(2^{\gamma }\right)}^2\lceil \left(\gamma +1/2\right)c_2}{{\gamma }^2\sqrt{\pi }{\left(\lceil \left(\gamma \right)\right)}^2}},$

${\Theta }_5^{[\ast ]}=\left(\begin{array}{c}{\displaystyle \frac{c_3}{\lceil \left(\gamma \right)\gamma }}-{\displaystyle \frac{{\left(2^{\gamma }\right)}^2\lceil \left(\gamma +1/2\right){c_2}^2}{{\gamma }^2\sqrt{\pi }{\left(\lceil \left(\gamma \right)\right)}^2}}-\\ {\displaystyle \frac{1}{2}}{\displaystyle \frac{c_3{\left(3^{\gamma }\right)}^3\sqrt{3}\lceil \left(\gamma +1/3\right)\lceil \left(\gamma +2/3\right)}{{\gamma }^2\sqrt{\pi }{\left(\lceil \left(\gamma \right)\right)}^2{\left(2^{\gamma }\right)}^2\lceil \left(\gamma +1/2\right)}}\\ +{\displaystyle \frac{{\left(2^{\gamma }\right)}^4{\left(\lceil \left(\gamma +1/2\right)\right)}^2{c_2}^2}{{\left(\lceil \left(\gamma \right)\right)}^3{\gamma }^3\pi }}\end{array}\right),$

$${\varsigma }^{\left[l\right]}=\left(-c_2+{\displaystyle \frac{{\left(2^{\gamma }\right)}^2\lceil \left(\gamma +1/2\right)c_2}{\lceil \left(\gamma \right)\gamma \sqrt{\pi }}}\right){\left[e^{\left[l\right]}\right]}^{\gamma +1}+{\Theta }_6^{[\ast ]}{\left[e^{\left[l\right]}\right]}^{2\gamma +1}+\dots$$

(22)

and ${\Theta }_6^{[\ast ]}=\left(\begin{array}{c}-c_3+{\displaystyle \frac{{\left(2^{\gamma }\right)}^2\lceil \left(\gamma +1/2\right){c_2}^2}{\lceil \left(\gamma \right)\gamma \sqrt{\pi }}}+\\ {\displaystyle \frac{1}{2}}{\displaystyle \frac{c_3{\left(3^{\gamma }\right)}^3\sqrt{3}\lceil \left(\gamma +1/3\right)\lceil \left(\gamma +2/3\right)}{\lceil \left(\gamma \right)\gamma \sqrt{\pi }{\left(2^{\gamma }\right)}^2\lceil \left(\gamma +1/2\right)}}\\ -{\displaystyle \frac{{\left(2^{\gamma }\right)}^4{\left(\lceil \left(\gamma +1/2\right)\right)}^2{c_2}^2}{{\left(\lceil \left(\gamma \right)\right)}^2{\gamma }^2\pi }}\end{array}\right).$ We have:

$$\hslash \left({\varsigma }^{\left[l\right]}\right)=\left[\begin{array}{c}{\Lambda }^{\left\{2\right\}}{\left[e^{\left[l\right]}\right]}^{4\gamma +1}+{\Lambda }^{\left\{3\right\}}{\left[e^{\left[l\right]}\right]}^{3\gamma +1}+-1/2{\Lambda }^{\left\{6\right\}}{\left[e^{\left[l\right]}\right]}^{2\gamma +1}\\ +{\Lambda }^{\left\{7\right\}}{\left[e^{\left[l\right]}\right]}^{\gamma +1}+\dots \end{array}\right],$$

(23)

where ${\Lambda }^{\left\{1\right\}}=\left(\begin{array}{c}2c_3{\left(\lceil \left(\gamma \right)\right)}^2{\gamma }^2{\pi }^{3/2}{\left(2^{\gamma }\right)}^2\lceil \left(\gamma +1/2\right)-\\ 2{\left(2^{\gamma }\right)}^4{\left(\lceil \left(\gamma +1/2\right)\right)}^2{c_2}^2\lceil \left(\gamma \right)\gamma \pi \\ -c_3{\left(3^{\gamma }\right)}^3\sqrt{3}\lceil \left(\gamma +1/3\right)\\ \lceil \left(\gamma +2/3\right)\lceil \left(\gamma \right)\gamma \pi +\\ 2{\left(2^{\gamma }\right)}^6{\left(\lceil \left(\gamma +1/2\right)\right)}^3{c_2}^2\sqrt{\pi }\end{array}\right),$

${\Lambda }^{\left\{2\right\}}={\displaystyle \frac{{c_2}^2\left(\lceil \left(\gamma \right)\gamma \sqrt{\pi }-{\left(2^{\gamma }\right)}^2\lceil \left(\gamma +1/2\right)\right)\left({\Lambda }^{\left\{1\right\}}\right)}{{\left(\lceil \left(\gamma \right)\right)}^3{\gamma }^3{\pi }^2{\left(2^{\gamma }\right)}^2\lceil \left(\gamma +1/2\right)}},$

${\Lambda }^{\left\{3\right\}}={\displaystyle \frac{{c_2}^3{\left(\lceil \left(\gamma \right)\gamma \sqrt{\pi }-{\left(2^{\gamma }\right)}^2\lceil \left(\gamma +1/2\right)\right)}^2}{{\left(\lceil \left(\gamma \right)\right)}^2{\gamma }^2\pi }},$

${\Lambda }^{\left\{4\right\}}=\left(\begin{array}{c}2c_3{\left(\lceil \left(\gamma \right)\right)}^2{\gamma }^2{\pi }^{3/2}{\left(2^{\gamma }\right)}^2\lceil \left(\gamma +1/2\right)\\ -2{\left(2^{\gamma }\right)}^4{\left(\lceil \left(\gamma +1/2\right)\right)}^2{c_2}^2\\ \lceil \left(\gamma \right)\gamma \pi -c_3{\left(3^{\gamma }\right)}^3\sqrt{3}\lceil \left(\gamma +1/3\right)\\ \lceil \left(\gamma +2/3\right)\lceil \left(\gamma \right)\gamma \pi +\\ 2{\left(2^{\gamma }\right)}^6{\left(\lceil \left(\gamma +1/2\right)\right)}^3{c_2}^2\sqrt{\pi },\end{array}\right),$

${\Lambda }^{\left\{5\right\}}={\left(\lceil \left(\gamma \right)\right)}^2{\gamma }^2{\pi }^{3/2}{\left(2^{\gamma }\right)}^2\lceil \left(\gamma +1/2\right),$

${\Lambda }^{\left\{6\right\}}={\displaystyle \frac{{\Lambda }^{\left\{4\right\}}}{{\Lambda }^{\left\{5\right\}}}},$ and

${\Lambda }^{\left\{7\right\}}=-{\displaystyle \frac{c_2\left(\lceil \left(\gamma \right)\gamma \sqrt{\pi }-{\left(2^{\gamma }\right)}^2\lceil \left(\gamma +1/2\right)\right)}{\lceil \left(\gamma \right)\gamma \sqrt{\pi }}}.$ Thus,

$${\displaystyle \frac{1}{\hslash \left(x^{\left[l\right]}\right)}}={e_n}^{-1}-c_2+\left({c_2}^2-c_3\right)e_n+\left(\begin{array}{c}{c}_2{c}_3\\ -c_4+\\ \left(\begin{array}{c}-{c_2}^2\\ +c_3\end{array}\right)c_2\end{array}\right){e_n}^2+\dots$$

(24)

$${\displaystyle \frac{\hslash \left({\varsigma }^{\left[l\right]}\right)}{\hslash \left(x^{\left[l\right]}\right)}}=\left[\begin{array}{c}{\Lambda }^{\left\{8\right\}}{\left[e^{\left[l\right]}\right]}^{\gamma }+{\Lambda }^{\left\{9\right\}}{\left[e^{\left[l\right]}\right]}^{\gamma +1}+{\Lambda }^{\left\{10\right\}}{\left[e^{\left[l\right]}\right]}^{2\gamma +1}\\ +{\Lambda }^{\left\{11\right\}}{\left[e^{\left[l\right]}\right]}^{3\gamma +1}+\dots \end{array}\right]$$

(25)

where ${\Lambda }^{\left\{8\right\}}=\left(-c_2+{\displaystyle \frac{{\left(2^{\gamma }\right)}^2\lceil \left(\gamma +1/2\right)c_2}{\lceil \left(\gamma \right)\gamma \sqrt{\pi }}}\right),$

${\Lambda }^{\left\{9\right\}}=\left(\begin{array}{c}-{\displaystyle \frac{{\left(2^{\gamma }\right)}^4{\left(\lceil \left(\gamma +1/2\right)\right)}^2{c_2}^2}{{\left(\lceil \left(\gamma \right)\right)}^2{\gamma }^2\pi }}+\\ {\displaystyle \frac{1}{2}}{\displaystyle \frac{c_3{\left(3^{\gamma }\right)}^3\sqrt{3}\lceil \left(\gamma +1/3\right)\lceil \left(\gamma +2/3\right)}{\lceil \left(\gamma \right)\gamma \sqrt{\pi }{\left(2^{\gamma }\right)}^2\lceil \left(\gamma +1/2\right)}}\\ -c_3+{c_2}^2\end{array}\right),$

${\Lambda }^{\left\{10\right\}}=\left(\begin{array}{c}-{\displaystyle \frac{2{c_2}^3{\left(2^{\gamma }\right)}^2\lceil \left(\gamma +1/2\right)}{\lceil \left(\gamma \right)\gamma \sqrt{\pi }}}+\\ {\displaystyle \frac{2{c_2}^3{\left(2^{\gamma }\right)}^4{\left(\lceil \left(\gamma +1/2\right)\right)}^2}{{\left(\lceil \left(\gamma \right)\right)}^2{\gamma }^2\pi }}-\\ {\displaystyle \frac{1}{2}}{\displaystyle \frac{c_3{\left(3^{\gamma }\right)}^3\sqrt{3}\lceil \left(\gamma +1/3\right)\lceil \left(\gamma +2/3\right)c_2}{\lceil \left(\gamma \right)\gamma \sqrt{\pi }{\left(2^{\gamma }\right)}^2\lceil \left(\gamma +1/2\right)}}\\ -{\displaystyle \frac{{\left(2^{\gamma }\right)}^2\lceil \left(\gamma +1/2\right)c_2{c}_3}{\lceil \left(\gamma \right)\gamma \sqrt{\pi }}}\\ +2c_2{c}_3\end{array}\right),$

${\Lambda }^{\left\{11\right\}}=\left(\begin{array}{c}{\displaystyle \frac{2{c_2}^4{\left(2^{\gamma }\right)}^4{\left(\lceil \left(\gamma +1/2\right)\right)}^2}{{\left(\lceil \left(\gamma \right)\right)}^2{\gamma }^2\pi }}-{c_2}^2{c}_3\\ -{\displaystyle \frac{2{c_2}^4{\left(2^{\gamma }\right)}^6{\left(\lceil \left(\gamma +1/2\right)\right)}^3}{{\left(\lceil \left(\gamma \right)\right)}^3{\gamma }^3{\pi }^{3/2}}}-\\ {\displaystyle \frac{{c_3}^2{\left(3^{\gamma }\right)}^3\sqrt{3}\lceil \left(\gamma +1/3\right)\lceil \left(\gamma +2/3\right)}{\lceil \left(\gamma \right)\gamma \sqrt{\pi }{\left(2^{\gamma }\right)}^2\lceil \left(\gamma +1/2\right)}}-\\ {\displaystyle \frac{1}{2}}{\displaystyle \frac{{c_2}^2{c}_3{\left(3^{\gamma }\right)}^3\sqrt{3}\lceil \left(\gamma +1/3\right)\lceil \left(\gamma +2/3\right)}{\lceil \left(\gamma \right)\gamma \sqrt{\pi }{\left(2^{\gamma }\right)}^2\lceil \left(\gamma +1/2\right)}}+\\ {\displaystyle \frac{{c_2}^2{c}_3{\left(3^{\gamma }\right)}^3\sqrt{3}\lceil \left(\gamma +1/3\right)\lceil \left(\gamma +2/3\right)}{{\left(\lceil \left(\gamma \right)\right)}^2{\gamma }^2\pi }}\\ -{\displaystyle \frac{{\left(2^{\gamma }\right)}^2\lceil \left(\gamma +1/2\right){c_2}^2{c}_3}{\lceil \left(\gamma \right)\gamma \sqrt{\pi }}}+{c_3}^2+\\ {\displaystyle \frac{{\left(2^{\gamma }\right)}^4{\left(\lceil \left(\gamma +1/2\right)\right)}^2{c_2}^2{c}_3}{{\left(\lceil \left(\gamma \right)\right)}^2{\gamma }^2\pi }}\\ -{\displaystyle \frac{{\left(2^{\gamma }\right)}^2\lceil \left(\gamma +1/2\right)c_2{c}_4}{\lceil \left(\gamma \right)\gamma \sqrt{\pi }}}+c_2{c}_4\end{array}\right),$

$$\left(\begin{array}{c}\left({\displaystyle \frac{\hslash \left({\varsigma }^{\left[l\right]}\right)}{\hslash \left(x^{\left[l\right]}\right)}}\right)+\\ 2{\left({\displaystyle \frac{\hslash \left({\varsigma }^{\left[l\right]}\right)}{\hslash \left(x^{\left[l\right]}\right)}}\right)}^2\\ +{\beta }^{[\ast ]}{\left({\displaystyle \frac{\hslash \left({\varsigma }^{\left[l\right]}\right)}{\hslash \left(x^{\left[l\right]}\right)}}\right)}^4\end{array}\right)={\Lambda }_1^{\{\ast \}}{\left[e^{\left[l\right]}\right]}^{\gamma }+{\Lambda }_2^{\{\ast \}}{\left[e^{\left[l\right]}\right]}^{\gamma +1}+{\Lambda }_3^{\{\ast \}}{\left[e^{\left[l\right]}\right]}^{2\gamma +1}+\dots$$

(26)

where ${\Lambda }_1^{\{\ast \}}=\left(-c_2+{\displaystyle \frac{{\left(2^{\gamma }\right)}^2\lceil \left(\gamma +1/2\right)c_2}{\lceil \left(\gamma \right)\gamma \sqrt{\pi }}}\right)$

${\Lambda }_2^{\{\ast \}}=\left(\begin{array}{c}{\displaystyle \frac{{\left(2^{\gamma }\right)}^4{\left(\lceil \left(\gamma +1/2\right)\right)}^2{c_2}^2}{{\left(\lceil \left(\gamma \right)\right)}^2{\gamma }^2\pi }}-4{\displaystyle \frac{4{\left(2^{\gamma }\right)}^2\lceil \left(\gamma +1/2\right){c_2}^2}{\lceil \left(\gamma \right)\gamma \sqrt{\pi }}}\\ +{\displaystyle \frac{1}{2}}{\displaystyle \frac{c_3{\left(3^{\gamma }\right)}^3\sqrt{3}\lceil \left(\gamma +1/3\right)\lceil \left(\gamma +2/3\right)}{\lceil \left(\gamma \right)\gamma \sqrt{\pi }{\left(2^{\gamma }\right)}^2\lceil \left(\gamma +1/2\right)}}-c_3+3{c_2}^2\end{array}\right)$

${\Lambda }_3^{\{\ast \}}=\left(\begin{array}{c}{\displaystyle \frac{2c_3{\left(3^{\gamma }\right)}^3\sqrt{3}\lceil \left(\gamma +1/3\right)\lceil \left(\gamma +2/3\right)c_2}{{\left(\lceil \left(\gamma \right)\right)}^2{\gamma }^2\pi }}-\\ 4{c_2}^3-5{\displaystyle \frac{{\left(2^{\gamma }\right)}^2\lceil \left(\gamma +1/2\right)c_2{c}_3}{\lceil \left(\gamma \right)\gamma \sqrt{\pi }}}\\ -{\displaystyle \frac{5}{2}}{\displaystyle \frac{c_3{\left(3^{\gamma }\right)}^3\sqrt{3}\lceil \left(\gamma +1/3\right)\lceil \left(\gamma +2/3\right)c_2}{\lceil \left(\gamma \right)\gamma \sqrt{\pi }{\left(2^{\gamma }\right)}^2\lceil \left(\gamma +1/2\right)}}-\\ {\displaystyle \frac{4{c_2}^3{\left(2^{\gamma }\right)}^6{\left(\lceil \left(\gamma +1/2\right)\right)}^3}{{\left(\lceil \left(\gamma \right)\right)}^3{\gamma }^3{\pi }^{3/2}}}\\ +2{\displaystyle \frac{{c_2}^3{\left(2^{\gamma }\right)}^2\lceil \left(\gamma +1/2\right)}{\lceil \left(\gamma \right)\gamma \sqrt{\pi }}}+6c_2{c}_3\\ +6{\displaystyle \frac{{c_2}^3{\left(2^{\gamma }\right)}^4{\left(\lceil \left(\gamma +1/2\right)\right)}^2}{{\left(\lceil \left(\gamma \right)\right)}^2{\gamma }^2\pi }}\end{array}\right).$ Therefore,

$$h_1^{[\ast ]}={\displaystyle \frac{\hslash \left(x^{\left[l\right]}\right)}{\left[{}_{\check{c}}⅁_{{\gamma }_1}^{\gamma }\right]\hslash \left(x^{\left[l\right]}\right)}}\left(\begin{array}{c}\left({\displaystyle \frac{\hslash \left({\varsigma }^{\left[l\right]}\right)}{\hslash \left(x^{\left[l\right]}\right)}}\right)+2{\left({\displaystyle \frac{\hslash \left({\varsigma }^{\left[l\right]}\right)}{\hslash \left(x^{\left[l\right]}\right)}}\right)}^2\\ +{\beta }^{[\ast ][\ast ]}{\left({\displaystyle \frac{\hslash \left({\varsigma }^{\left[l\right]}\right)}{\hslash \left(x^{\left[l\right]}\right)}}\right)}^4\end{array}\right),$$

(27)

$$h_1^{[\ast ]}={\Lambda }^{\left\{12\right\}}{\left[e^{\left[l\right]}\right]}^{\gamma +1}+{\Lambda }^{\left\{13\right\}}{\left[e^{\left[l\right]}\right]}^{2\gamma +1}+{\Lambda }^{\left\{14\right\}}{\left[e^{\left[l\right]}\right]}^{3\gamma +1}+\dots$$

(28)

where ${\Lambda }^{\left\{12\right\}}=\left({\displaystyle \frac{{\left(2^{\gamma }\right)}^2\lceil \left(\gamma +1/2\right)c_2}{{\gamma }^2\sqrt{\pi }{\left(\lceil \left(\gamma \right)\right)}^2}}-{\displaystyle \frac{c_2}{\lceil \left(\gamma \right)\gamma }}\right),$

${\Lambda }^{\left\{13\right\}}=\left(\begin{array}{c}-{\displaystyle \frac{2{\left(2^{\gamma }\right)}^2\lceil \left(\gamma +1/2\right){c_2}^2}{{\gamma }^2\sqrt{\pi }{\left(\lceil \left(\gamma \right)\right)}^2}}+{\displaystyle \frac{1}{2}}{\displaystyle \frac{c_3{\left(3^{\gamma }\right)}^3\sqrt{3}\lceil \left(\gamma +1/3\right)\lceil \left(\gamma +2/3\right)}{{\gamma }^2\sqrt{\pi }{\left(\lceil \left(\gamma \right)\right)}^2{\left(2^{\gamma }\right)}^2\lceil \left(\gamma +1/2\right)}}\\ -{\displaystyle \frac{c_3}{\lceil \left(\gamma \right)\gamma }}+{\displaystyle \frac{2{c_2}^2}{\lceil \left(\gamma \right)\gamma }}\end{array}\right),$

${\Lambda }^{\left\{14\right\}}=\left(\begin{array}{c}{\displaystyle \frac{9{c_2}^3{\left(2^{\gamma }\right)}^4{\left(\lceil \left(\gamma +1/2\right)\right)}^2}{{\left(\lceil \left(\gamma \right)\right)}^3{\gamma }^3\pi }}-\\ {\displaystyle \frac{4{c_2}^3{\left(2^{\gamma }\right)}^6{\left(\lceil \left(\gamma +1/2\right)\right)}^3}{{\left(\lceil \left(\gamma \right)\right)}^4{\gamma }^4{\pi }^{3/2}}}+\\ {\displaystyle \frac{c_3{\left(3^{\gamma }\right)}^3\sqrt{3}\lceil \left(\gamma +1/3\right)\lceil \left(\gamma +2/3\right)c_2}{{\left(\lceil \left(\gamma \right)\right)}^3{\gamma }^3\pi }}-\\ {\displaystyle \frac{3}{2}}{\displaystyle \frac{c_3{\left(3^{\gamma }\right)}^3\sqrt{3}\lceil \left(\gamma +1/3\right)\lceil \left(\gamma +2/3\right)c_2}{{\gamma }^2\sqrt{\pi }{\left(\lceil \left(\gamma \right)\right)}^2{\left(2^{\gamma }\right)}^2\lceil \left(\gamma +1/2\right)}}+\\ {\displaystyle \frac{4c_2{c}_3}{\lceil \left(\gamma \right)\gamma }}-{\displaystyle \frac{{c_2}^3}{\lceil \left(\gamma \right)\gamma }}-\\ {\displaystyle \frac{3{\left(2^{\gamma }\right)}^2\lceil \left(\gamma +1/2\right)c_2{c}_3}{{\gamma }^2\sqrt{\pi }{\left(\lceil \left(\gamma \right)\right)}^2}}-{\displaystyle \frac{4{c_2}^3{\left(2^{\gamma }\right)}^2\lceil \left(\gamma +1/2\right)}{{\gamma }^2\sqrt{\pi }{\left(\lceil \left(\gamma \right)\right)}^2}}\end{array}\right),$

$${\left(\lceil \left(\gamma +1\right)h_1^{[\ast ]}\right)}^{{\displaystyle \frac{1}{\gamma }}}={\Lambda }^{\left\{15\right\}}{\left[e^{\left[l\right]}\right]}^{2\gamma +1}+{\Lambda }^{\left\{16\right\}}{\left[e^{\left[l\right]}\right]}^{3\gamma +1}+\dots$$

(29)

where ${\Lambda }^{\left\{15\right\}}=\left({\displaystyle \frac{3{\left(2^{\gamma }\right)}^2\lceil \left(\gamma +1/2\right){c_2}^2}{\lceil \left(\gamma \right)\gamma \sqrt{\pi }}}-{\displaystyle \frac{{\left(2^{\gamma }\right)}^4{\left(\lceil \left(\gamma +1/2\right)\right)}^2{c_2}^2}{{\left(\lceil \left(\gamma \right)\right)}^2{\gamma }^2\pi }}-2{c_2}^2\right),$

${\Lambda }^{\left\{16\right\}}=\left(\begin{array}{c}-{\displaystyle \frac{9{c_2}^3{\left(2^{\gamma }\right)}^4{\left(\lceil \left(\gamma +1/2\right)\right)}^2}{{\left(\lceil \left(\gamma \right)\right)}^2{\gamma }^2\pi }}+\\ 4{\displaystyle \frac{4{c_2}^3{\left(2^{\gamma }\right)}^6{\left(\lceil \left(\gamma +1/2\right)\right)}^3}{{\left(\lceil \left(\gamma \right)\right)}^3{\gamma }^3{\pi }^{3/2}}}\\ -{\displaystyle \frac{c_3{\left(3^{\gamma }\right)}^3\sqrt{3}\lceil \left(\gamma +1/3\right)\lceil \left(\gamma +2/3\right)c_2}{{\left(\lceil \left(\gamma \right)\right)}^2{\gamma }^2\pi }}+\\ {\displaystyle \frac{3}{2}}{\displaystyle \frac{c_3{\left(3^{\gamma }\right)}^3\sqrt{3}\lceil \left(\gamma +1/3\right)\lceil \left(\gamma +2/3\right)c_2}{\lceil \left(\gamma \right)\gamma \sqrt{\pi }{\left(2^{\gamma }\right)}^2\lceil \left(\gamma +1/2\right)}}\\ -4c_2{c}_3+{c_2}^3+\\ {\displaystyle \frac{3{\left(2^{\gamma }\right)}^2\lceil \left(\gamma +1/2\right)c_2{c}_3}{\lceil \left(\gamma \right)\gamma \sqrt{\pi }}}+4{\displaystyle \frac{4{c_2}^3{\left(2^{\gamma }\right)}^2\lceil \left(\gamma +1/2\right)}{\lceil \left(\gamma \right)\gamma \sqrt{\pi }}}\end{array}\right),$

$$v^{\left[l\right]}=x^{\left[l\right]}-{\left(\lceil (\gamma +1){\displaystyle \frac{\hslash \left(x^{\left[l\right]}\right)}{{\hslash }^{\prime }\left(x^{\left[l\right]}\right)}}\left(\begin{array}{c}\left({\displaystyle \frac{\hslash \left({\varsigma }^{\left[l\right]}\right)}{\hslash \left(x^{\left[l\right]}\right)}}\right)+2{\left({\displaystyle \frac{\hslash \left({\varsigma }^{\left[l\right]}\right)}{\hslash \left(x^{\left[l\right]}\right)}}\right)}^2\\ +{\beta }^{[\ast ][\ast ]}{\left({\displaystyle \frac{\hslash \left({\varsigma }^{\left[l\right]}\right)}{\hslash \left(x^{\left[l\right]}\right)}}\right)}^4\end{array}\right)\right)}^{1/\gamma },$$

(30)

$$e^{\left[\ast \right]}=e^{\left[\ast \right]}=\left(\begin{array}{c}{\pounds }_1^{[\ast ]}+{\pounds }_2^{[\ast ]}\\ +{\pounds }_3^{[\ast ]}\\ +{\pounds }_4^{[\ast ]}\\ +{\pounds }_5^{[\ast ]}\\ +{\pounds }_6^{[\ast ]}\end{array}\right){\left[e^{\left[l\right]}\right]}^{3\gamma +1}+O\left({\left[e^{\left[l\right]}\right]}^{4\gamma +1}\right),$$

(31)

where ${\pounds }_1^{[\ast ]}=-9{\displaystyle \frac{9{c_2}^3{\left(2^{\gamma }\right)}^4{\left(\lceil \left(\gamma +1/2\right)\right)}^2}{{\left(\lceil \left(\gamma \right)\right)}^2{\gamma }^2\pi }}$, ${\pounds }_2^{[\ast ]}=+{\displaystyle \frac{4{c_2}^3{\left(2^{\gamma }\right)}^6{\left(\lceil \left(\gamma +1/2\right)\right)}^3}{{\left(\lceil \left(\gamma \right)\right)}^3{\gamma }^3{\pi }^{3/2}}}$, ${\pounds }_3^{[\ast ]}=-{\displaystyle \frac{c_3{\left(3^{\gamma }\right)}^3\sqrt{3}\lceil \left(\gamma +1/3\right)\lceil \left(\gamma +2/3\right)c_2}{{\left(\lceil \left(\gamma \right)\right)}^2{\gamma }^2\pi }}$, ${\pounds }_4^{[\ast ]}={\displaystyle \frac{3}{2}}{\displaystyle \frac{c_3{\left(3^{\gamma }\right)}^3\sqrt{3}\lceil \left(\gamma +1/3\right)\lceil \left(\gamma +2/3\right)c_2}{\lceil \left(\gamma \right)\gamma \sqrt{\pi }{\left(2^{\gamma }\right)}^2\lceil \left(\gamma +1/2\right)}}$,

${\pounds }_5^{[\ast ]}=-4c_2{c}_3+{c_2}^3$, ${\pounds }_6^{[\ast ]}={\displaystyle \frac{3{\left(2^{\gamma }\right)}^2\lceil \left(\gamma +1/2\right)c_2{c}_3}{\lceil \left(\gamma \right)\gamma \sqrt{\pi }}}+{\displaystyle \frac{4{c_2}^3{\left(2^{\gamma }\right)}^2\lceil \left(\gamma +1/2\right)}{\lceil \left(\gamma \right)\gamma \sqrt{\pi }}}.$

Thus, the theorem is proven.    □

2.2. Dynamical Analysis of the Fractional Scheme MM^γ

Ensuring the robustness and reliability of solving nonlinear equations requires a careful examination of the dynamical behavior of solution techniques [19]. Dynamical analysis in this context refers to a method’s ability to converge on the true root from an initial guess, even when there are small perturbations or errors in the calculations. Single root-finding algorithms generally have local convergence around the root, making them effective when the initial guess is close enough. However, their stability depends on the characteristics of the function and the initial approximation. If the function is poorly behaved or the initial guess is far from the root, these methods may diverge or converge to an incorrect solution.

The stability of a root-finding method can be analyzed using concepts from complex dynamical systems, which examine how sensitive the root is to variations in the input, as well as the convergence criterion, which evaluates how quickly the method approaches the root. The choice of parameters plays a crucial role in the method’s performance, particularly within the convergence zone (the red region) on the parametric plane. Stability and faster convergence are achieved when optimal parameters are selected from this region. However, choosing parameters near the boundary of this zone, or from the divergence zone, may lead to reduced stability or slower convergence rates. Selecting parameters from the divergence region can result in failure or divergence, as the method may become unstable. To mitigate the effects of computational errors and ensure reliable performance, stability is typically maintained by carefully balancing the method’s convergence properties with the appropriate selection of the initial guess, function parameters, and stopping criteria.

The following rational map $R\left(x\right)=\left(x-a\right)\left(x-b\right)$ is obtained as:

$$R\left(x\right)=y=\left[x-2^{-{\gamma }^{-1}}{\left({\displaystyle \frac{\lceil \left(\gamma +1\right)\left(x^2-1\right)}{x}}\right)}^{{\gamma }^{-1}}-2^{-{\gamma }^{-1}}{\wp }^{[\ast ]}\right],$$

(32)

where ${\wp }^{[\ast ]}=$${\left({\displaystyle \frac{\lceil \left(\gamma +1\right)}{{\left(x^2-1\right)}^{3}x}}\right)}^{{\gamma }^{-1}}\left(\begin{array}{c}{\wp }_{11}^{\left[i\right]}+{\wp }_{12}^{\left[i\right]}\\ +{\wp }_{13}^{\left[i\right]}\\ +{\wp }_{14}^{\left[i\right]}+{\wp }_{15}^{\left[i\right]}\\ {\wp }_{16}^{\left[i\right]}+{\wp }_{17}^{\left[i\right]}\\ +{\wp }_{18}^{\left[i\right]}\\ +{\wp }_{19}^{\left[i\right]}+{\wp }_{20}^{\left[i\right]}\\ {\wp }_{21}^{\left[i\right]}+{\wp }_{22}^{\left[i\right]}\\ +{\wp }_{23}^{\left[i\right]}\\ +{\wp }_{14}^{\left[i\right]}+{\wp }_{15}^{\left[i\right]}\end{array}\right),$

$${\wp }_{11}^{\left[i\right]}=\left(\begin{array}{c}3+{16}^{{\displaystyle \frac{\gamma -1}{\gamma }}}{\left({\displaystyle \frac{\lceil \left(\gamma +1\right)\left(x^2-1\right)}{x}}\right)}^{4{\gamma }^{-1}}{\beta }^{[\ast ][\ast ]}{x}^4-\\ {42}^{{\displaystyle \frac{\gamma -1}{\gamma }}}{\left({\displaystyle \frac{\lceil \left(\gamma +1\right)\left(x^2-1\right)}{x}}\right)}^{{\gamma }^{-1}}{\beta }^{[\ast ][\ast ]}{x}^7+\\ {122}^{{\displaystyle \frac{\gamma -1}{\gamma }}}{\left({\displaystyle \frac{\lceil \left(\gamma +1\right)\left(x^2-1\right)}{x}}\right)}^{{\gamma }^{-1}}{\beta }^{[\ast ][\ast ]}{x}^5-\\ {124}^{-{\gamma }^{-1}}{\left({\displaystyle \frac{\lceil \left(\gamma +1\right)\left(x^2-1\right)}{x}}\right)}^{2{\gamma }^{-1}}{\beta }^{[\ast ][\ast ]}{x}^4\end{array}\right),$$

$${\wp }_{12}^{\left[i\right]}=\left(\begin{array}{c}-{122}^{{\displaystyle \frac{\gamma -1}{\gamma }}}{\left({\displaystyle \frac{\lceil \left(\gamma +1\right)\left(x^2-1\right)}{x}}\right)}^{{\gamma }^{-1}}{\beta }^{[\ast ][\ast ]}{x}^3+\\ {124}^{-{\gamma }^{-1}}{\left({\displaystyle \frac{\lceil \left(\gamma +1\right)\left(x^2-1\right)}{x}}\right)}^{2{\gamma }^{-1}}{\beta }^{[\ast ]}{x}^2\\ +{42}^{{\displaystyle \frac{\gamma -1}{\gamma }}}{\left({\displaystyle \frac{\lceil \left(\gamma +1\right)\left(x^2-1\right)}{x}}\right)}^{{\gamma }^{-1}}{\beta }^{[\ast ]}x-\\ {48}^{{\displaystyle \frac{\gamma -1}{\gamma }}}{\left({\displaystyle \frac{\lceil \left(\gamma +1\right)\left(x^2-1\right)}{x}}\right)}^{3{\gamma }^{-1}}{\beta }^{[\ast ]}{x}^5\end{array}\right),$$

$$\left(\begin{array}{c}+{74}^{{\displaystyle \frac{\gamma -1}{\gamma }}}{\left({\displaystyle \frac{\lceil \left(\gamma +1\right)\left(x^2-1\right)}{x}}\right)}^{2{\gamma }^{-1}}{\beta }^{[\ast ][\ast ]}{x}^6+\\ {48}^{{\displaystyle \frac{\gamma -1}{\gamma }}}{\left({\displaystyle \frac{\lceil \left(\gamma +1\right)\left(x^2-1\right)}{x}}\right)}^{3{\gamma }^{-1}}{\beta }^{[\ast ][\ast ]}{x}^3+\\ {464}^{-{\gamma }^{-1}}{\left({\displaystyle \frac{\lceil \left(\gamma +1\right)\left(x^2-1\right)}{x}}\right)}^{6{\gamma }^{-1}}{\beta }^{[\ast ][\ast ]}{x}^2\end{array}\right),$$

$${\wp }_{13}^{\left[i\right]}=\left(\begin{array}{c}+{616}^{-{\gamma }^{-1}}{\left({\displaystyle \frac{\lceil \left(\gamma +1\right)\left(x^2-1\right)}{x}}\right)}^{4{\gamma }^{-1}}{\beta }^{[\ast ]}{x}^4-\\ {124}^{{\displaystyle \frac{\gamma -1}{\gamma }}}{\left({\displaystyle \frac{\lceil \left(\gamma +1\right)\left(x^2-1\right)}{x}}\right)}^{2{\gamma }^{-1}}{\beta }^{[\ast ]}{x}^4\\ -{1216}^{-{\gamma }^{-1}}{\left({\displaystyle \frac{\lceil \left(\gamma +1\right)\left(x^2-1\right)}{x}}\right)}^{4{\gamma }^{-1}}{\beta }^{[\ast ]}{x}^2\\ +{64}^{{\displaystyle \frac{\gamma -1}{\gamma }}}{\left({\displaystyle \frac{\lceil \left(\gamma +1\right)\left(x^2-1\right)}{x}}\right)}^{2{\gamma }^{-1}}{\beta }^{[\ast ]}{x}^2\end{array}\right),$$

$${\wp }_{14}^{\left[i\right]}=\left(\begin{array}{c}-4^{{\displaystyle \frac{\gamma -1}{\gamma }}}{\left({\displaystyle \frac{\lceil \left(\gamma +1\right)\left(x^2-1\right)}{x}}\right)}^{2{\gamma }^{-1}}-\\ 4^{{\displaystyle \frac{\gamma -1}{\gamma }}}{\left({\displaystyle \frac{\lceil \left(\gamma +1\right)\left(x^2-1\right)}{x}}\right)}^{3{\gamma }^{-1}}{2}^{{\displaystyle \frac{\gamma -1}{\gamma }}}{x}^5+\\ {84}^{-{\gamma }^{-1}}{\left({\displaystyle \frac{\lceil \left(\gamma +1\right)\left(x^2-1\right)}{x}}\right)}^{3{\gamma }^{-1}}{2}^{{\displaystyle \frac{\gamma -1}{\gamma }}}{x}^3\end{array}\right),$$

$${\wp }_{15}^{\left[i\right]}=\left(\begin{array}{c}-4^{{\displaystyle \frac{\gamma -1}{\gamma }}}{\left({\displaystyle \frac{\lceil \left(\gamma +1\right)\left(x^2-1\right)}{x}}\right)}^{3{\gamma }^{-1}}{2}^{{\displaystyle \frac{\gamma -1}{\gamma }}}x-\\ {124}^{{\displaystyle \frac{-2+\gamma }{\gamma }}}{\left({\displaystyle \frac{\lceil \left(\gamma +1\right)\left(x^2-1\right)}{x}}\right)}^{4{\gamma }^{-1}}{\beta }^{[\ast ]}{x}^2\\ +{124}^{{\displaystyle \frac{-2+\gamma }{\gamma }}}{\left({\displaystyle \frac{\lceil \left(\gamma +1\right)\left(x^2-1\right)}{x}}\right)}^{4{\gamma }^{-1}}{\beta }^{[\ast ]}{x}^4\\ -4^{{\displaystyle \frac{\gamma -1}{\gamma }}}{\left({\displaystyle \frac{\lceil \left(\gamma +1\right)\left(x^2-1\right)}{x}}\right)}^{2{\gamma }^{-1}}{\beta }^{[\ast ]}\end{array}\right),$$

$${\wp }_{16}^{\left[i\right]}=\left(\begin{array}{c}-{416}^{-{\gamma }^{-1}}{\left({\displaystyle \frac{\lceil \left(\gamma +1\right)\left(x^2-1\right)}{x}}\right)}^{4{\gamma }^{-1}}{x}^2+\\ {24}^{{\displaystyle \frac{\gamma -1}{\gamma }}}{\left({\displaystyle \frac{\lceil \left(\gamma +1\right)\left(x^2-1\right)}{x}}\right)}^{2{\gamma }^{-1}}{x}^2\\ +{616}^{-{\gamma }^{-1}}{\left({\displaystyle \frac{\lceil \left(\gamma +1\right)\left(x^2-1\right)}{x}}\right)}^{4{\gamma }^{-1}}{\beta }^{[\ast ]}-\\ 4^{{\displaystyle \frac{2\gamma -1}{\gamma }}}{\left({\displaystyle \frac{\lceil \left(\gamma +1\right)\left(x^2-1\right)}{x}}\right)}^{2{\gamma }^{-1}}{x}^4\end{array}\right),$$

$${\wp }_{17}^{\left[i\right]}=\left(\begin{array}{c}-{464}^{-{\gamma }^{-1}}{\left({\displaystyle \frac{\lceil \left(\gamma +1\right)\left(x^2-1\right)}{x}}\right)}^{6{\gamma }^{-1}}{\beta }^{[\ast ]}+\\ {216}^{-{\gamma }^{-1}}{\left({\displaystyle \frac{\lceil \left(\gamma +1\right)\left(x^2-1\right)}{x}}\right)}^{4{\gamma }^{-1}}{x}^4\\ +3{34}^{{\displaystyle \frac{\gamma -1}{\gamma }}}{\left({\displaystyle \frac{\lceil \left(\gamma +1\right)\left(x^2-1\right)}{x}}\right)}^{2{\gamma }^{-1}}{x}^6+\\ {256}^{-{\gamma }^{-1}}{\left({\displaystyle \frac{\lceil \left(\gamma +1\right)\left(x^2-1\right)}{x}}\right)}^{8{\gamma }^{-1}}{\beta }^{[\ast ]}\end{array}\right),$$

$${\wp }_{18}^{\left[i\right]}=\left(\begin{array}{c}+{216}^{-{\gamma }^{-1}}{\left({\displaystyle \frac{\lceil \left(\gamma +1\right)\left(x^2-1\right)}{x}}\right)}^{4{\gamma }^{-1}}+3x^8-\\ {52}^{{\displaystyle \frac{\gamma -1}{\gamma }}}{\left({\displaystyle \frac{\lceil \left(\gamma +1\right)\left(x^2-1\right)}{x}}\right)}^{{\gamma }^{-1}}{x}^7\\ +4^{-{\gamma }^{-1}}{\left({\displaystyle \frac{\lceil \left(\gamma +1\right)\left(x^2-1\right)}{x}}\right)}^{2{\gamma }^{-1}}{x}^6+\\ {122}^{{\displaystyle \frac{\gamma -1}{\gamma }}}{\left({\displaystyle \frac{\lceil \left(\gamma +1\right)\left(x^2-1\right)}{x}}\right)}^{{\gamma }^{-1}}{x}^5\end{array}\right),$$

$${\wp }_{19}^{\left[i\right]}=\left(\begin{array}{c}+{62}^{-{\gamma }^{-1}}{\left({\displaystyle \frac{\lceil \left(\gamma +1\right)\left(x^2-1\right)}{x}}\right)}^{{\gamma }^{-1}}{x}^5-\\ {154}^{-{\gamma }^{-1}}{\left({\displaystyle \frac{\lceil \left(\gamma +1\right)\left(x^2-1\right)}{x}}\right)}^{2{\gamma }^{-1}}{x}^4\\ -{122}^{{\displaystyle \frac{\gamma -1}{\gamma }}}{\left({\displaystyle \frac{\lceil \left(\gamma +1\right)\left(x^2-1\right)}{x}}\right)}^{{\gamma }^{-1}}{x}^3-\\ {62}^{-{\gamma }^{-1}}{\left({\displaystyle \frac{\lceil \left(\gamma +1\right)\left(x^2-1\right)}{x}}\right)}^{{\gamma }^{-1}}{x}^3\end{array}\right),$$

$${\wp }_{20}^{\left[i\right]}=\left(\begin{array}{c}+{154}^{-{\gamma }^{-1}}{\left({\displaystyle \frac{\lceil \left(\gamma +1\right)\left(x^2-1\right)}{x}}\right)}^{2{\gamma }^{-1}}{x}^2+\\ 5{52}^{{\displaystyle \frac{\gamma -1}{\gamma }}}{\left({\displaystyle \frac{\lceil \left(\gamma +1\right)\left(x^2-1\right)}{x}}\right)}^{{\gamma }^{-1}}x\\ +{\beta }^{[\ast ]}{x}^8-4{\beta }^{[\ast ]}{x}^6+\\ 6{\beta }^{[\ast ]}{x}^4-4{\beta }^{[\ast ]}{x}^2+18x^4\\ -12x^6-4^{-{\gamma }^{-1}}{\left({\displaystyle \frac{\lceil \left(\gamma +1\right)\left(x^2-1\right)}{x}}\right)}^{2{\gamma }^{-1}}\end{array}\right),$$

$${\wp }_{21}^{\left[i\right]}=\left(\begin{array}{c}-{1216}^{-{\gamma }^{-1}}{\left({\displaystyle \frac{\lceil \left(\gamma +1\right)\left(x^2-1\right)}{x}}\right)}^{5{\gamma }^{-1}}{2}^{{\displaystyle \frac{\gamma -1}{\gamma }}}{\beta }^{[\ast ]}{x}^3-\\ {124}^{-{\gamma }^{-1}}{\left({\displaystyle \frac{\lceil \left(\gamma +1\right)\left(x^2-1\right)}{x}}\right)}^{3{\gamma }^{-1}}{2}^{{\displaystyle \frac{\gamma -1}{\gamma }}}{\beta }^{[\ast ]}{x}^5\\ +{1216}^{-{\gamma }^{-1}}{\left({\displaystyle \frac{\lceil \left(\gamma +1\right)\left(x^2-1\right)}{x}}\right)}^{5{\gamma }^{-1}}{2}^{{\displaystyle \frac{\gamma -1}{\gamma }}}{\beta }^{[\ast ]}x\\ +{244}^{-{\gamma }^{-1}}{\left({\displaystyle \frac{\lceil \left(\gamma +1\right)\left(x^2-1\right)}{x}}\right)}^{3{\gamma }^{-1}}\end{array}\right),$$

$${\wp }_{22}^{\left[i\right]}=\left(\begin{array}{c}{2}^{{\displaystyle \frac{\gamma -1}{\gamma }}}{\beta }^{[\ast ]}{x}^3-\\ {124}^{-{\gamma }^{-1}}{\left({\displaystyle \frac{\lceil \left(\gamma +1\right)\left(x^2-1\right)}{x}}\right)}^{3{\gamma }^{-1}}{2}^{{\displaystyle \frac{\gamma -1}{\gamma }}}{\beta }^{[\ast ]}x-\\ 4^{{\displaystyle \frac{\gamma -1}{\gamma }}}{\left({\displaystyle \frac{\lceil \left(\gamma +1\right)\left(x^2-1\right)}{x}}\right)}^{5{\gamma }^{-1}}{8}^{{\displaystyle \frac{\gamma -1}{\gamma }}}{\beta }^{[\ast ]}{x}^3\\ +{616}^{-{\gamma }^{-1}}{\left({\displaystyle \frac{\lceil \left(\gamma +1\right)\left(x^2-1\right)}{x}}\right)}^{6{\gamma }^{-1}}{4}^{{\displaystyle \frac{\gamma -1}{\gamma }}}{\beta }^{[\ast ]}{x}^2\end{array}\right),$$

$$\begin{array}{ccc}\hfill {\wp }_{23}^{\left[i\right]}& =& -{464}^{-{\gamma }^{-1}}{\left({\displaystyle \frac{\lceil \left(\gamma +1\right)\left(x^2-1\right)}{x}}\right)}^{7{\gamma }^{-1}}{2}^{{\displaystyle \frac{\gamma -1}{\gamma }}}{\beta }^{[\ast ]}x\hfill \\ & & +{\beta }^{[\ast ]}-12x^2.\hfill \end{array}$$

For $\gamma \approx 1$, we have

$$R\left(x\right)=-{\displaystyle \frac{{\phi }_1^{[\ast ]}{x}^{10}+{\phi }_2^{[\ast ]}{x}^8+{\phi }_3^{[\ast ]}{x}^6+{\phi }_4^{[\ast ]}{x}^4+(5x^2-1){\beta }^{[\ast ]}}{512x^9}},$$

(33)

where ${\phi }_1^{[\ast ]}=\left({\beta }^{[\ast ]}-160\right)$, ${\phi }_1^{[\ast ]}=\left(-5{\beta }^{[\ast ]}-480\right)$, ${\phi }_1^{[\ast ]}=\left(10{\beta }^{[\ast ]}+160\right)$, ${\phi }_1^{[\ast ]}=\left(-10{\beta }^{[\ast ]}-32\right)$ and a, b belong to $\mathbb{C}.$ Thus, $R\left(x\right)$ depends on a, b, and x. Using the Möbius transformation $M\left(x\right)={\displaystyle \frac{x-a}{x-b}},$ we see that $R\left(x\right)$ is conjugate with the operator, and for $\gamma \approx 1,$ we have [19]

$$\mathrm{O}\left(x\right)={\displaystyle \frac{x^4\left(-x^6-8x^5-27x^4-48x^3+{\beta }^{[\ast ]}x-47x^2-24x-5\right)}{{\left({\beta }^{[\ast ]}{x}^5-5x^6-24x^5-47x^{4.}-48x^3-27x^2-8x-1\right)}^3}},$$

(34)

which is independent of a and b, as shown using Maple with 0 or ∞. Thus, $R\left(x\right)$ fits precisely with:

$$R^{[\ast ]}\left(x\right)={\displaystyle \frac{{\sum }_{i=0}^n\left(a_i{x}^i\right)}{{\sum }_{i=0}^n\left(a_{n-i}{x}^i\right)}};{\left\{a_i\right\}}_{i=0}^n\in R,$$

(35)

which exhibits interesting properties [27].

Proposition 1.

The fixed points of $\mathrm{O}\left(\tau \right)$ are as follows:

• ${\tau }_0=0;{\tau }_{\infty }=\infty$ are super-attracting points.
• ${\tau }_1=1$ is a repelling point, and
• The critical points 0 and 1 are super-attracting and repelling points, respectively, for $\gamma \approx 1$, and are given by:

  $$\begin{array}{ccc}\hfill C_1^{[\ast ]}& =& {\displaystyle \frac{-16.0+{\beta }^{[\ast ]}}{16.0}}-{\displaystyle \frac{\sqrt{{\beta }^{[\ast ]}\left(5{\beta }^{[\ast ]}-288\right)}}{16\sqrt{5.0}}}-{\displaystyle \frac{\sqrt{{\beta }_1^{[\ast ]}}}{8\sqrt{10.0}}},\hfill \\ \hfill C_2^{[\ast ]}& =& {\displaystyle \frac{-16.0+{\beta }^{[\ast ]}}{16.0}}-{\displaystyle \frac{\sqrt{{\beta }^{[\ast ]}\left(5{\beta }^{[\ast ]}-288\right)}}{16\sqrt{5.0}}}+{\displaystyle \frac{\sqrt{{\beta }_1^{[\ast ]}}}{8\sqrt{10.0}}},\hfill \\ \hfill C_3^{[\ast ]}& =& {\displaystyle \frac{-16.0+{\beta }^{[\ast ]}}{16.0}}-{\displaystyle \frac{\sqrt{{\beta }^{[\ast ]}\left(5{\beta }^{[\ast ]}-288\right)}}{16\sqrt{5.0}}}+{\displaystyle \frac{\sqrt{{\beta }_2^{[\ast ]}}}{2}},\hfill \\ \hfill C_4^{[\ast ]}& =& {\displaystyle \frac{-16.0+{\beta }^{[\ast ]}}{16.0}}-{\displaystyle \frac{\sqrt{{\beta }^{[\ast ]}\left(5{\beta }^{[\ast ]}-288\right)}}{16\sqrt{5.0}}}-{\displaystyle \frac{\sqrt{{\beta }_2^{[\ast ]}}}{2}},\hfill \end{array}$$

  where ${\beta }_1^{[\ast ]}=-224{\beta }^{[\ast ]}+5{\left({\beta }^{[\ast ]}\right)}^2-{\displaystyle \frac{4608\sqrt{5}{\beta }^{[\ast ]}}{\sqrt{5{\beta }^{[\ast ]}-288}}}+{\displaystyle \frac{268\sqrt{5}{\left({\beta }^{[\ast ]}\right)}^{2/3}}{\sqrt{5{\beta }^{[\ast ]}-288}}}-{\displaystyle \frac{5\sqrt{5}{\left({\beta }^{[\ast ]}\right)}^{5/2}}{\sqrt{5{\beta }^{[\ast ]}-288}}}$, ${\beta }_2^{[\ast ]}=-2+\left({\displaystyle \frac{16-{\left({\beta }^{[\ast ]}\right)}^2}{32}}\right)-{\displaystyle \frac{30+2{\beta }^{[\ast ]}}{5}}+{\displaystyle \frac{{\beta }_3^{[\ast ]}}{\sqrt{5{\beta }^{[\ast ]}-288}}}$, ${\beta }_3^{[\ast ]}=2\sqrt{5}\left(\left({\beta }^{[\ast ]}-32\right)-{\displaystyle \frac{{\left(64-{\beta }^{[\ast ]}\right)}^2}{64}}+{\displaystyle \frac{2\left(16-{\beta }^{[\ast ]}\right)\left(15+{\beta }^{[\ast ]}\right)}{5}}\right)$.

The stability function for $\gamma \approx 1$ in the single root-finding method is given by:

$$\left(1,{\beta }^{[\ast ]}\right)={\displaystyle \frac{-81,920+512{\beta }^{[\ast ]}}{{\left(-160+{\beta }^{[\ast ]}\right)}^2}}.$$

(36)

For other values of the fractional parameter $\gamma \in \left[0,1\right]$, the stability regions are examined, as shown in Figure 1a–c.

The stability region of the fractional iterative method MM^γ for finding the root of nonlinear equations illustrates the initial guesses or parameters for which the method converges to the correct root, as seen in Figure 1a–c. The behavior in which the fractions contribute to assessing the method’s resilience and effectiveness in solving nonlinear equations under a range of $\gamma$ values is seen in Figure 1a,c. The role of the fractional parameter $\gamma$ in determining the method’s stability and effectiveness is highlighted in these figures. The method is more stable when $\gamma =1$, with stability decreasing as $\gamma$ approaches zero.

By using critical points as starting values, parametric planes are generated (see Figure 2a–c), which depict the behavior of the method’s iterations, including patterns such as convergence, divergence, and chaotic behavior (Figure 3a–h and Figure 4a–c). These figures highlight stable and unstable behavior at various ${\beta }^{[\ast ]}$ values and $\gamma =1$. In Figure 2, Figure 3 and Figure 4, small circles represent fixed points, squares with stars indicate super-attracting fixed points, and squares denote critical points.

Parametric planes are essential for understanding the convergence behavior of iterative methods used to solve nonlinear equations. They provide insights into the stability and efficiency of iterations, helping to identify optimal parameters for fast and accurate convergence. These planes also reveal divergence zones, allowing for the avoidance of parameters that may cause instability or failure. By analyzing parametric planes, the optimal parameters for minimizing the number of iterations and improving solution reliability can be determined.

The figures also illustrate the method’s sensitivity to parameter changes, showing whether the approach is robust or requires careful tuning. Additionally, they display the nature of convergence—whether linear, super-linear, or quadratic—indicating how quickly the method reaches the solution. This information is crucial for selecting and optimizing parameters to ensure efficient and reliable problem-solving.

Understanding the dynamics of iterative methods is essential for solving nonlinear problems over time. Dynamic analysis helps evaluate how initial estimates affect outcomes, convergence rates, and overall stability. By analyzing these dynamics, one can ensure the reliability and efficiency of the approach, avoiding issues like slow convergence or divergence. In applied sciences and engineering, dynamic analysis plays a critical role in generating accurate solutions from iterative processes.

Dynamic planes offer insights into solution stability, convergence patterns, and behavior throughout the iterations. These planes reveal zones of attraction, where iterations converge, and regions of repulsion, where they diverge. They also highlight fixed points, cycles, and chaotic behavior, providing a deeper understanding of the methods’ sensitivity to initial conditions. This analysis is crucial for optimizing iterative methods and ensuring consistency in solving nonlinear equations. Similar to parametric planes, dynamic planes are generated to evaluate the stability of different schemes. Figure 3a–h shows the stable behavior of the schemes for various values of $\gamma$ and ${\beta }^{[\ast ]}.$ Choosing ${\beta }^{[\ast ]}$ values from the divergence zone of the parametric planes leads to unstable behavior, as illustrated in Figure 4a,c for different $\gamma$ values.

Using the dynamical and parametric plane notation, we find the best parameter value that accelerates the rate of convergence of the simple root-finding method. Then, using the newly developed stable fractional-order scheme MM^γ as a correction factor, we propose a novel parallel fractional scheme for analyzing (1) in the following section.

3. Development and Analysis of Fractional Parallel Scheme

Among derivative-free simultaneous methods, the Weierstrass–Durand–Kerner (also known as Weierstrass–Dochive) method [28] is one of the most effective approaches. It is defined as:

$$x_i^{\left[l+1\right]}=x_i^{\left[l\right]}-w\left(x_i^{\left[l\right]}\right),$$

(37)

where

$$w\left(x_i^{\left[l\right]}\right)={\displaystyle \frac{\hslash \left(x_i^{\left[l\right]}\right)}{\underset{\stackrel{j=1}{j\ne i}}{\stackrel{n}{\Pi }}(x_i^{\left[l\right]}-x_j^{\left[l\right]})}},(i,j=1,\dots ,n),$$

(38)

is known as Weierstrass’ correction. This method exhibits local quadratic convergence.

In 1977, Ehrlich [29] introduced the following third-order simultaneous method:

$$x_i^{\left[l+1\right]}=x_i^{\left[l\right]}-{\displaystyle \frac{1}{\frac{1}{N_i\left(x_i^{\left[l\right]}\right)}-{\displaystyle \sum _{\stackrel{j=1}{j\ne i}}^n}\left(\frac{1}{\left(x_i^{\left[l\right]}-x_j^{\left[l\right]}\right)}\right)}},$$

(39)

where $x_j^{\left[l\right]}=u_j^{\left[l\right]}$ is used as a correction in (15). Petkovic et al. [30] improved the convergence order of this method from three to six:

$$x_i^{\left[l+1\right]}=x_i^{\left[l\right]}-{\displaystyle \frac{1}{\frac{1}{N_i\left(x_i^{\left[l\right]}\right)}-{\displaystyle \sum _{\stackrel{j=1}{j\ne i}}^n}\left(\frac{1}{\left(x_i^{\left[l\right]}-x_j^{\left[l\right]}\right)}\right)}},$$

(40)

with the correction term: $u_j^{\left[l\right]}=x_i^{\left[l\right]}-{\displaystyle \frac{\hslash \left(s_j^{\left[l\right]}\right)-\hslash \left(x_i^{\left[l\right]}\right)}{2\ast \hslash \left(s_j^{\left[l\right]}\right)-\hslash \left(x_i^{\left[l\right]}\right)}}{\displaystyle \frac{\hslash \left(x_i^{\left[l\right]}\right)}{{\hslash }^{\prime }\left(x_i^{\left[l\right]}\right)}}$, where $s_j^{\left[l\right]}=x_j^{\left[l\right]}-{\displaystyle \frac{\hslash \left(x_j^{\left[l\right]}\right)}{{\hslash }^{\prime }\left(x_j^{\left[l\right]}\right)}}.$

Petkovic et al. [31] further accelerated the convergence order of this method from three to ten, giving the method (PM^[∗]) as:

$$x_i^{\left[l+1\right]}=x_i^{\left[l\right]}-{\displaystyle \frac{1}{\frac{1}{N_i\left(x_i^{\left[l\right]}\right)}-{\displaystyle \sum _{\stackrel{j=1}{j\ne i}}^n}\left(\frac{1}{\left(x_i^{\left[l\right]}-T_j^{\left[l\right]}\right)}\right)}},$$

(41)

where $T_j^{\left[l\right]}={\gamma }_j^{\left[l\right]}-{\displaystyle \frac{\left({\varsigma }_j^{\left[l\right]}-{\gamma }_j^{\left[l\right]})\hslash \left({\gamma }_j^{\left[l\right]}\right)\right)\left(\frac{\hslash \left(x_j^{\left[l\right]}\right)}{{\hslash }^{\prime }\left(x_j^{\left[l\right]}\right)}\right)}{{\left(\hslash \left(x_j^{\left[l\right]}\right)-\hslash \left({\gamma }_j^{\left[l\right]}\right)\right)}^2}}\left[\hslash \left({\varsigma }_j^{\left[l\right]}\right)-{\displaystyle \frac{\hslash \left(x_j^{\left[l\right]}\right)}{\hslash \left({\varsigma }_j^{\left[l\right]}\right)-\hslash \left({\gamma }_j^{\left[l\right]}\right)}}\right]$, ${\gamma }_j^{\left[l\right]}={\varsigma }_j^{\left[l\right]}-{\displaystyle \frac{\hslash \left(x_j^{\left[l\right]}\right)\hslash \left({\varsigma }_j^{\left[l\right]}\right)\left(\frac{\hslash \left(x_j^{\left[l\right]}\right)}{{\hslash }^{\prime }\left(x_j^{\left[l\right]}\right)}\right)}{{\left(\hslash \left(x_j^{\left[l\right]}\right)-\hslash \left({\varsigma }_j^{\left[l\right]}\right)\right)}^2}}$, ${\varsigma }_j^{\left[l\right]}=x_j^{\left[l\right]}-{\displaystyle \frac{\hslash \left(x_j^{\left[l\right]}\right)}{{\hslash }^{\prime }\left(x_j^{\left[l\right]}\right)}}.$

Using the method (38) in (16), we developed a new parallel scheme MM^[∗] for simultaneously finding all solutions of Equation (1). The updated iterative scheme is:

$$s_i^{\left[l\right]}={\varsigma }_i^{\left[l\right]}-{\displaystyle \frac{\hslash \left({\varsigma }_i^{\left[l\right]}\right)}{\underset{\stackrel{j=1}{j\ne i}}{\stackrel{n}{\Pi }}({\varsigma }_i^{\left[l\right]}-{\varsigma }_j^{\left[l\right]})}}\left(\begin{array}{c}\left({\displaystyle \frac{\hslash \left({\varsigma }_i^{\left[l\right]}\right)}{\hslash \left(x_i^{\left[l\right]}\right)}}\right)+2{\left({\displaystyle \frac{\hslash \left({\varsigma }_i^{\left[l\right]}\right)}{\hslash \left(x_i^{\left[l\right]}\right)}}\right)}^2\\ +{\beta }^{[\ast ]}{\left({\displaystyle \frac{\hslash \left({\varsigma }_i^{\left[l\right]}\right)}{\hslash \left(x_i^{\left[l\right]}\right)}}\right)}^4\end{array}\right),$$

(42)

where ${\varsigma }_i^{\left[k\right]}=x_i^{\left[l\right]}-{\displaystyle \frac{\hslash \left(x_i^{\left[l\right]}\right)}{\underset{\stackrel{j=1}{j\ne i}}{\stackrel{n}{\Pi }}(x_i^{\left[l\right]}-v_j^{\left[l\right]})}}$ and

$v_j^{\left[l\right]}=x_j^{\left[l\right]}-{\left(\lceil (\gamma +1){\displaystyle \frac{\hslash \left(x_j^{\left[l\right]}\right)}{{\hslash }^{\prime }\left(x_j^{\left[l\right]}\right)}}\left(\begin{array}{c}\left({\displaystyle \frac{\hslash \left({\varsigma }_j^{\left[l\right]}\right)}{\hslash \left(x_j^{\left[l\right]}\right)}}\right)\\ +2{\left({\displaystyle \frac{\hslash \left({\varsigma }_j^{\left[l\right]}\right)}{\hslash \left(x_j^{\left[l\right]}\right)}}\right)}^2\\ +{\beta }^{[\ast ]}{\left({\displaystyle \frac{\hslash \left({\varsigma }_j^{\left[l\right]}\right)}{\hslash \left(x_j^{\left[l\right]}\right)}}\right)}^4\end{array}\right)\right)}^{1/\gamma },$

${\varsigma }_j^{\left[l\right]}=x_j^{\left[l\right]}-{\left(\lceil (\gamma +1){\displaystyle \frac{\hslash \left(x_j^{\left[l\right]}\right)}{\left[{}_{\check{c}}⅁_{{\gamma }_1}^{\gamma }\right]\hslash \left(x_j^{\left[l\right]}\right)}}\right)}^{1/\gamma };{\beta }^{[\ast ]}\in R.$

The following theorem defines the local order of convergence for the inverse fractional scheme.

Theorem 3.

Let ${\zeta }_1,\dots ,{\zeta }_{\sigma }$ be simple roots of a nonlinear equation, and assume that the initial distinct estimates $x_1^{\left[0\right]},\dots ,x_n^{\left[0\right]}$ are sufficiently close to the true roots. Then, the MM^[∗] method achieves a convergence order of $6\gamma +3$.

Proof. 

Let $e_i=x_i^{\left[l\right]}-{\zeta }_i,e_{\varsigma }={\varsigma }_i^{\left[l\right]}-{\zeta }_i$ and $e_i^{[\ast ]}=s_i^{\left[l\right]}-{\zeta }_i$ represent the errors in $x_i^{\left[l\right]},$${\varsigma }_i^{\left[l\right]},$ and $s_i^{\left[l\right]}$, respectively. From the first-step of MM^[∗], we have:

$${\varsigma }_i^{\left[l\right]}-{\zeta }_i=x_i^{\left[l\right]}-{\zeta }_i-{\displaystyle \frac{\hslash \left(x_j^{\left[l\right]}\right)}{\underset{\stackrel{j=1}{j\ne i}}{\stackrel{n}{\Pi }}(x_i^{\left[l\right]}-v_j^{\left[l\right]})}},$$

(43)

$$e_i^{\prime }=e_i-e_i\underset{\stackrel{j=1}{j\ne i}}{\prod ^n}\left[{\displaystyle \frac{x_i^{\left[l\right]}-{\zeta }_j}{x_i^{\left[l\right]}-v_j^{\left[l\right]}}}\right],$$

(44)

$$e_{\varsigma }=e_i\left(1-\underset{\stackrel{j=1}{j\ne i}}{\prod ^n}\left[{\displaystyle \frac{x_i^{\left[l\right]}-{\zeta }_j}{x_i^{\left[l\right]}-v_j^{\left[l\right]}}}\right]\right),$$

(45)

where $\underset{\stackrel{j=1}{j\ne i}}{\prod ^n}\left[{\displaystyle \frac{x_i^{\left[l\right]}-{\zeta }_j}{x_i^{\left[l\right]}-u_j^{\left[l\right]}}}\right]=1+{\displaystyle \frac{v_j^{\left[l\right]}-{\zeta }_j}{x_i^{\left[l\right]}-v_j^{\left[l\right]}}}={\left(1+O\left({\left[e_i\right]}^{3\gamma +1}\right)\right)}^{l+1}$ and, from (31), $v_j^{\left[l\right]}-{\zeta }_j=O\left(\left|{\left[e_i\right]}^{3\gamma +1}\right|\right)$. Thus,

$$e_{\varsigma }=e_i\left(1+(l-1)O\left(\left|{\left[e_i\right]}^{3\gamma +1}\right|\right)\right),$$

(46)

Assuming $\left|e_i^{2\gamma +1}\right|=\left|e_j^{2\gamma +1}\right|,$ we have

$$e_{\varsigma }=O\left(\left|{\left[e_i\right]}^{3\gamma +2}\right|\right).$$

(47)

In the second step, we obtain:

$$s_i^{\left[l\right]}-{\zeta }_i=x_i^{\left[l\right]}-{\zeta }_i-{\displaystyle \frac{\hslash \left({\varsigma }_i^{\left[l\right]}\right)}{\underset{\stackrel{j=1}{j\ne i}}{\stackrel{n}{\Pi }}({\varsigma }_i^{\left[l\right]}-{\varsigma }_j^{\left[l\right]})}}\left(\begin{array}{c}\left({\displaystyle \frac{\hslash \left({\varsigma }_i^{\left[l\right]}\right)}{\hslash \left(x_i^{\left[l\right]}\right)}}\right)\\ +2{\left({\displaystyle \frac{\hslash \left({\varsigma }_i^{\left[l\right]}\right)}{\hslash \left(x_i^{\left[l\right]}\right)}}\right)}^2\\ +{\beta }^{[\ast ]}{\left({\displaystyle \frac{\hslash \left({\varsigma }_i^{\left[l\right]}\right)}{\hslash \left(x_i^{\left[l\right]}\right)}}\right)}^4\end{array}\right),$$

(48)

leading to:

$$e_i^{[\ast ]}=e_{\varsigma }-{\displaystyle \frac{\hslash \left({\varsigma }_i^{\left[l\right]}\right)}{\underset{\stackrel{j=1}{j\ne i}}{\stackrel{n}{\Pi }}({\varsigma }_i^{\left[l\right]}-{\varsigma }_j^{\left[l\right]})}}\left(\begin{array}{c}\left({\displaystyle \frac{\hslash \left({\varsigma }_i^{\left[l\right]}\right)}{\hslash \left(x_i^{\left[l\right]}\right)}}\right)+2{\left({\displaystyle \frac{\hslash \left({\varsigma }_i^{\left[l\right]}\right)}{\hslash \left(x_i^{\left[l\right]}\right)}}\right)}^2\\ +{\beta }^{[\ast ]}{\left({\displaystyle \frac{\hslash \left({\varsigma }_i^{\left[l\right]}\right)}{\hslash \left(x_i^{\left[l\right]}\right)}}\right)}^4\end{array}\right),$$

(49)

or

$$e_i^{[\ast ]}=e_{\varsigma }-e_{\varsigma }\underset{\stackrel{j=1}{j\ne i}}{\prod ^n}\left[{\displaystyle \frac{x_i^{\left[l\right]}-{\zeta }_j}{x_i^{\left[l\right]}-u_j^{\left[l\right]}}}\right]\left(\begin{array}{c}\left({\displaystyle \frac{e_{\varsigma }{\vartheta }^{[\ast ]}}{e_i}}\right)+\\ 2{\left({\displaystyle \frac{e_{\varsigma }{\vartheta }^{[\ast ]}}{e_i}}\right)}^2\\ +{\beta }^{[\ast ]}{\left({\displaystyle \frac{e_{\varsigma }{\vartheta }^{[\ast ]}}{e_i}}\right)}^4\end{array}\right),$$

(50)

where ${\displaystyle \frac{\hslash \left({\varsigma }_i^{\left[l\right]}\right)}{\hslash \left(x_i^{\left[l\right]}\right)}=\frac{e_{\varsigma }{\vartheta }^{[\ast ]}}{e_i}}$. Thus,

$$e_i^{[\ast ]}=e_{\varsigma }\left(1-{\displaystyle \frac{e_i}{e_{\varsigma }}}\underset{\stackrel{j=1}{j\ne i}}{\prod ^n}\left[{\displaystyle \frac{x_i^{\left[l\right]}-{\zeta }_j}{x_i^{\left[l\right]}-u_j^{\left[l\right]}}}\right]\left(\begin{array}{c}\left({\displaystyle \frac{{\vartheta }^{[\ast ]}}{e_i}}\right)+\\ 2e_{\varsigma }{\left({\displaystyle \frac{{\vartheta }^{[\ast ]}}{e_i}}\right)}^2+\\ {\beta }^{[\ast ]{\left(e_{\varsigma }\right)}^3}{\left({\displaystyle \frac{{\vartheta }^{[\ast ]}}{e_i}}\right)}^4\end{array}\right)\right).$$

(51)

leading to:

$$e_i^{[\ast ]}=e_{\varsigma }\left(1-{\vartheta }^{[\ast ]}\left(\begin{array}{c}\left({\displaystyle \frac{{\vartheta }^{[\ast ]}}{e_i}}\right)+2e_{\varsigma }{\left({\displaystyle \frac{{\vartheta }^{[\ast ]}}{e_i}}\right)}^2\\ +{\beta }^{[\ast ]{\left(e_{\varsigma }\right)}^3}\left({\displaystyle \frac{{\left({\vartheta }^{[\ast ]}\right)}^4}{{\left(e_i\right)}^3}}\right)\end{array}\right)\right),$$

(52)

$$e_i^{[\ast ]}=e_{\varsigma }\left(1-O\left(\left|{\left[e_i\right]}^{3\gamma +1}\right|\right)\left({\displaystyle \frac{{\Theta }^{[\ast \ast ]}}{{\left(e_i\right)}^3}}\right)\right),$$

(53)

where ${\Theta }^{[\ast \ast ]}=\left({\left(e_i\right)}^2{\vartheta }^{[\ast ]}+2e_i{e}_{\varsigma }{\vartheta }^{[\ast ]}+{\beta }^{[\ast ]}{\left(e_{\varsigma }\right)}^3{\left({\vartheta }^{[\ast ]}\right)}^4\right).$ We conclude that:

$$e_i^{[\ast ]}=e_{\varsigma }\left(O\left(\left|{\left[e_i\right]}^{3\gamma +1}\right|\right)\right),$$

(54)

$$e_i^{[\ast ]}=O\left(\left|{\left[e_i\right]}^{6\gamma +3}\right|\right).$$

(55)

This completes the proof.    □

4. Numerical Results

In this section, we demonstrate the effectiveness and stability of the method by examining several engineering applications using the following termination criteria implemented in Maple 18:

$$\left(i\right)e_i=∥x_i^{\left[1+1\right]}-x_i^{\left[1\right]}∥<\in ={10}^{-18}$$

(56)

where $e_i$ represents the residual error norm-2. We compare our method with the PM^[∗] method and the approach presented in [32], referred to as ZM^[∗]. Using Algorithm 1, we find all solutions to the nonlinear problem (1) simultaneously.

Algorithm 1 Parallel fractional numerical scheme MM^[∗] implementation on computer for solving nonlinear problems.

Algorithm 1 Parallel fractional numerical scheme MM[∗].

4.1. Example 1

Fractional differential equations are particularly useful in modeling complex engineering systems, such as the quarter car suspension model, where memory effects and damping properties play a critical role. Unlike traditional integer-order equations, fractional equations provide more accurate representations of systems with hereditary behaviors, such as viscoelasticity and energy dissipation in suspension systems. In this context, fractional differential equations are employed to design suspension models with improved performance and stability, making them powerful tools for addressing the challenges in automotive engineering and control system design. This leads to the following fractional differential equation:

$$\left\{\begin{array}{c}\left[{}_{\check{c}}⅁_{{\gamma }_1}^{\gamma }\right]y\left(x\right)+\alpha {\left(U^{[\ast ]}-y\left(x\right)\right)}^4=0,\\ y\left(0\right)=0,0<\gamma \le 1,x\ge 0\end{array}\right.$$

(57)

Using the method from [32], we simulate Equation (57) using the following nonlinear approximation:

$$y\left(x\right)\approx \left[\begin{array}{c}{\displaystyle \frac{\alpha {\left(U^{[\ast ]}\right)}^4}{\lceil \left(\gamma +1\right)}}{x}^{\gamma }-{\displaystyle \frac{4{\alpha }^2{\left(U^{[\ast ]}\right)}^7}{\lceil \left(3\gamma +1\right)}}{x}^{2\gamma }\\ +\left({\displaystyle \frac{16{\alpha }^3{\left(U^{[\ast ]}\right)}^{10}\lceil \left(2\gamma +1\right)}{{\lceil }^2\left(\gamma +1\right)\lceil \left(3\gamma +1\right)}}\right)x^{3\gamma }+\dots \end{array}\right].$$

(58)

If we set $U^{[\ast ]}=1$, $\alpha =1$, and $\gamma =1$, Equation (58) becomes:

$$y_5\left(x\right)\approx x-2x^2+{\displaystyle \frac{14x^3}{3}}-{\displaystyle \frac{35x^4}{3}}+{\displaystyle \frac{91}{3}}{x}^5.$$

(59)

Equation (59) has five solutions:

$$\begin{array}{ccc}\hfill {\zeta }_1& =& 0.0,\hfill \\ \hfill {\zeta }_{2,3}& =& -0.152018637\pm 03931364228i,\hfill \\ \hfill {\zeta }_{4,5}& =& 0.314432633\pm 0.25883527180i.\hfill \end{array}$$

The Caputo-type derivative of Equation (59) is:

$$\left[{}_{\check{c}}⅁_{{\gamma }_1}^{\gamma }\right]y_5\left(x\right)=\left[\begin{array}{c}{\displaystyle \frac{91}{3}}{\displaystyle \frac{\lceil \left(6\right)}{\lceil (6-\gamma )}}{x}^{5-\gamma }-{\displaystyle \frac{35}{3}}{\displaystyle \frac{\lceil \left(5\right)}{\lceil (5-\gamma )}}{x}^{4-\gamma }\\ +{\displaystyle \frac{14}{3}}{\displaystyle \frac{\lceil \left(4\right)}{\lceil (4-\gamma )}}{x}^{3-\gamma }\\ -2{\displaystyle \frac{\lceil \left(3\right)}{\lceil (3-\gamma )}}{x}^{2-\gamma }+{\displaystyle \frac{\lceil \left(2\right)}{\lceil (2-\gamma )}}{x}^{1-\gamma }\end{array}\right].$$

(60)

For the dynamical analysis, we set $\gamma ={\displaystyle \frac{1}{2}}$ in Equation (58) and obtain the following nonlinear equation:

$$y_5\left(x\right)\approx {\displaystyle \frac{2\sqrt{x}}{\sqrt{\pi }}}-4x+{\displaystyle \frac{32\left(3+2\pi \right)x^{\frac{3}{2}}}{3{\pi }^{\frac{3}{2}}}}+\left(68-{\displaystyle \frac{60}{\pi }}\right)x^2,$$

(61)

The exact root of Equation (61) is zero. The dynamical analysis of the simple root-finding method for this equation includes measurements of elapsed time, percentage convergence, total points, and the number of iterations for different ${\beta }^{[\ast ]}$ values, as shown in

Table 1. Dynamical outcomes of MM^[γ] for solving (61).

${\mathit{\beta }}^{\mathbf{[}\mathbf{\ast }\mathbf{]}}$	Iterations	Total Points	Converging Points	Elapsed Time
0.9	21	640,000	54%	2.14584165
0.7	22	640,000	59%	3.65483535
0.5	27	640,000	58%	5.36528456

and Figure 5a–c.

Using the information from Table 1 and the optimal parameter values from the dynamical analysis, these values are applied in parallel schemes to find all solutions to the nonlinear problems. Using initial guesses close to the exact solution $\in ={10}^{-2}$, the convergence rate of the simultaneous methods increases, leading to faster convergence to the exact roots with fewer iterations.

In

Table 2. Numerical results for Example 1.

Err	PM[∗]	ZM[∗]	MM[∗]
$e_1^{\left[5\right]}$	$9.31\times {10}^{-45}$	$1.11\times {10}^{-99}$	$1.1\times {10}^{-104}$
$e_2^{\left[5\right]}$	$2.09\times {10}^{-55}$	$3.02\times {10}^{-76}$	$0.0$
$e_3^{\left[5\right]}$	$0.17\times {10}^{-76}$	$0.12\times {10}^{-85}$	$0.0$
$e_4^{\left[5\right]}$	$0.07\times {10}^{-39}$	$6.27\times {10}^{-101}$	$4.9\times {10}^{-165}$
$e_4^{\left[5\right]}$	$1.07\times {10}^{-39}$	$6.27\times {10}^{-101}$	$4.9\times {10}^{-165}$
${\sigma }_i^{\left[4\right]}$	$8.127654$	$9.001923$	$9.125432$

, ${\sigma }_i^{\left[4\right]}$ represents the local computational order of convergence [33]. The table clearly demonstrates that when $\gamma \approx 1$, the residual error for each root, along with the order and rate of convergence, is significantly improved compared to existing methods. This indicates higher efficiency and stability of the approach.

Table 2 and

Table 3. Outcomes of the fractional scheme MM^[∗] for different γ values for solving engineering application 1.

$\mathit{\gamma }$	Iteration	$∥{\mathit{x}}_{\mathit{i}}^{\left[1+1\right]}-{\mathit{x}}_{\mathit{i}}^{\left[1\right]}∥$	${\mathit{\sigma }}_{\mathit{i}}^{[\mathit{n}-1]}$	CPU Time
$0.1$	36	$0.2\times {10}^{-5}$	$3.2132$	$0.097$
$0.2$	24	$1.4\times {10}^{-6}$	$4.2654$	$0.096$
$0.3$	22	$0.7\times {10}^{-6}$	$4.8236$	$0.076$
$0.4$	12	$0.5\times {10}^{-4}$	$5.4455$	$0.066$
$0.5$	10	$5.4\times {10}^{-6}$	$6.0432$	$0.035$
$0.6$	08	$1.3\times {10}^{-7}$	$6.6653$	$0.038$
$0.7$	08	$0.6\times {10}^{-9}$	$7.2124$	$0.015$
$0.8$	07	$0.3\times {10}^{-11}$	$7.5434$	$0.004$
$0.9$	06	$2.6\times {10}^{-16}$	$8.1325$	$0.001$

clearly show that across a range of $\gamma$, MM^[∗] outperforms ZM^[∗] and PM^[∗] in terms of efficiency, computation time (CPU time), number of iterations (n), and residual error $\left(∥x_i^{\left[1+1\right]}-x_i^{\left[1\right]}∥\right)$. The convergence order of a fractional parallel iterative scheme is influenced by the fractional parameter, which captures memory effects. The convergence order of parallel systems is determined by the order of the fractional derivative, with faster convergence observed as $\gamma$ increases. Our method, MM^[∗], outperforms the existing methods ZM^[∗] and PM^[∗] in the literature, achieving a convergence order of 9.002313 when $\gamma \approx 0.99$.

4.2. Example 2

In order to extend classical differential equations and enable more accurate descriptions of physical systems, such as electrical circuits, fluid flow, and viscoelastic materials, fractional calculus is added. This method also enables a more flexible representation of anomalous diffusion and wave propagation, which improves predictions in fields such as biology and materials science and leads to the following fractional initial value problem [34]:

$$\left\{\begin{array}{c}\left[{}_{\check{c}}⅁_{{\gamma }_1}^{\gamma }\right]y\left(x\right)+e^{y\left(x\right)}=0,\\ y\left(0\right)=0,y^{\prime }\left(0\right)=0,1<\gamma \le 2,x\ge 0\end{array}\right.$$

(62)

Using the method from [32], we simulate Equation (62) with the following nonlinear approximation:

$$y_n\left(x\right)\approx \left[\begin{array}{c}-{\displaystyle \frac{1}{\lceil \left(\gamma +1\right)}}{x}^{\gamma }+{\displaystyle \frac{1}{\lceil \left(3\gamma +1\right)}}{x}^{2\gamma }\\ +\left({\displaystyle \frac{2{\left(\lceil \left(\gamma +1\right)\right)}^2\lceil \left(2\gamma +1\right)}{{2\lceil }^2\left(\gamma +1\right)\lceil \left(3\gamma +1\right)}}\right)x^{3\gamma }+\dots \end{array}\right].$$

(63)

For $\gamma =1$, Equation (63) simplifies to:

$$y_3\left(x\right)\approx -x+{\displaystyle \frac{1}{2}}{x}^2-2x^3,$$

(64)

Equation (64) has the following solutions:

$$\begin{array}{ccc}\hfill {\zeta }_1& =& 0.0,\hfill \\ \hfill {\zeta }_{2,3}& =& 0.1250000\pm 0.695970545i.\hfill \end{array}$$

The Caputo-type derivative of Equation (63) is:

$$\left[{}_{\check{c}}⅁_{{\gamma }_1}^{\gamma }\right]y_5\left(x\right)=\left[\begin{array}{c}-{\displaystyle \frac{2\lceil \left(4\right)}{\lceil (4-\gamma )}}{x}^{3-\gamma }+{\displaystyle \frac{1}{2}}{\displaystyle \frac{\lceil \left(3\right)}{\lceil (3-\gamma )}}{x}^{2-\gamma }\\ -{\displaystyle \frac{\lceil \left(2\right)}{\lceil (2-\gamma )}}{x}^{1-\gamma }\end{array}\right].$$

(65)

For the dynamical analysis, we set $\gamma ={\displaystyle \frac{1}{2}}$ in Equation (63) and obtain the following nonlinear equation:

$$y_5\left(x\right)\approx x-{\displaystyle \frac{2\sqrt{x}}{\sqrt{\pi }}}-{\displaystyle \frac{2\left(\frac{1}{2}\pi +1\right)x^{\frac{3}{2}}}{\pi \left(\frac{3}{4}\sqrt{\pi }\right)}},$$

(66)

The exact root of Equation (66) is zero. The dynamical analysis of the simple root-finding method for this equation includes elapsed time, percentage convergence, total points, and the number of iterations for different ${\beta }^{[\ast ]}$ values, as shown in

Table 4. Dynamical outcomes of MM^[γ] for solving (66).

${\mathit{\beta }}^{\mathbf{[}\mathbf{\ast }\mathbf{]}}$	Iterations	Total Points	Converging Points	Elapsed Time
0.9	21	640,000	54%	1.33464564
0.7	18	640,000	59%	2.67454366
0.5	20	640,000	58%	5.233554675

.

Using initial guesses close to the exact solution (up to $\in ={10}^{-2}$), the convergence rate of the simultaneous methods improves, reaching the exact roots with fewer iterations. Table 4 clearly shows that when $\gamma \approx 1$, the residual error for each root, as well as the convergence order and rate of convergence, are superior to existing methods, indicating better efficiency and stability. The numerical outcomes for random initial guesses are presented in

Table 5. Numerical results for Example 2.

Err	PM[∗]	ZM[∗]	MM[∗]
$e_1^{\left[3\right]}$	$0.2\times {10}^{-53}$	$1.1\times {10}^{-51}$	$0.0$
$e_2^{\left[3\right]}$	$1.4\times {10}^{-44}$	$0.0$	$0.0$
$e_3^{\left[3\right]}$	$0.6\times {10}^{-45}$	$5.0\times {10}^{-57}$	$5.0\times {10}^{-65}$
$e_4^{\left[3\right]}$	$0.6\times {10}^{-45}$	$5.0\times {10}^{-57}$	$5.0\times {10}^{-65}$
$e_5^{\left[3\right]}$	$0.6\times {10}^{-45}$	$5.0\times {10}^{-57}$	$5.0\times {10}^{-65}$
${\sigma }_i^{\left[2\right]}$	$7.993434$	$8.991923$	$9.001232$

and Figure 6a–c.

In Table 5, ${\sigma }_i^{[n-1]}$ represents the local computational order of convergence. The table clearly demonstrates that when $\gamma \approx 1$, the residual error for each root, as well as the convergence order and rate of convergence, are significantly improved compared to existing methods. This highlights the method’s enhanced efficiency and stability.

The numerical results in

Table 6. Outcomes of the fractional scheme MM^[∗] for different γ values for solving engineering application 2.

$\mathit{\gamma }$	n	$∥{\mathit{x}}_{\mathit{i}}^{\left[1+1\right]}-{\mathit{x}}_{\mathit{i}}^{\left[1\right]}∥$	${\mathit{\sigma }}_{\mathit{i}}^{[\mathit{n}-1]}$	CPU Time
$0.1$	05	$0.3\times {10}^{-4}$	$3.0976$	$0.097$
$0.2$	05	$1.5\times {10}^{-8}$	$4.2004$	$0.086$
$0.3$	06	$0.7\times {10}^{-8}$	$4.8546$	$0.088$
$0.4$	06	$7.4\times {10}^{-9}$	$5.4132$	$0.076$
$0.5$	05	$0.5\times {10}^{-11}$	$5.9962$	$0.063$
$0.6$	06	$0.4\times {10}^{-15}$	$6.2353$	$0.058$
$0.7$	06	$0.1\times {10}^{-19}$	$6.9864$	$0.045$
$0.8$	05	$0.4\times {10}^{-20}$	$7.1243$	$0.038$
$0.9$	05	$3.7\times {10}^{-27}$	$7.9876$	$0.001$

show that across a range of $\gamma$ values, MM^[∗] outperforms both ZM^[∗∗] and PM^[∗] in terms of efficiency, CPU time (in seconds), number of iterations, and residual error. The convergence order of a fractional parallel iterative scheme is influenced by the fractional parameter, which captures memory effects. The convergence order of parallel systems is determined by the order of the fractional derivative, with faster convergence observed as $\gamma$ increases. Our method, MM^[∗], outperforms the existing methods ZM^[∗] and PM^[∗] from the literature, achieving a convergence order of 9.002313 when $\gamma \approx 0.99$.

5. Conclusions

We developed a new Caputo-type fractional scheme with a convergence order of $3\gamma +1$ to find simple solutions to Equation (1) and extended it into a fractional parallel numerical scheme for solving all nonlinear fractional problems. The parallel schemes were shown to have a convergence order of $6\gamma +3$ based on our convergence analysis. To improve the convergence rate of these fractional schemes, dynamical planes were used to select optimal initial guesses, ensuring convergence to exact solutions. Several nonlinear problems were tested to assess the stability and consistency of MM^[∗] in comparison to ZM^[∗∗], PM^[∗]. The numerical results demonstrate that the MM^[∗] method is more stable and consistent in terms of residual error, CPU time, and error graphs for different $\gamma$ values, outperforming both ZM^[∗∗] and PM^[∗]. Future research will focus on applying our fractional-order parallel technique to advanced control systems, biological processes, financial models, complex networks, and chaotic dynamics. Additionally, integrating fractional-order parallel schemes with machine learning [35,36] and enhancing computational efficiency could further expand its applicability and impact.