Robust Fractional Quantum Two-Step Schemes with Enhanced Stability for Nonlinear Equations

1. Introduction

Nonlinear equations constitute a fundamental mathematical framework for modeling equilibrium states, steady configurations, and implicit relationships arising in science and engineering. Given a nonlinear operator

$$f\left(x\right)=0,f:\Omega \subset \mathbb{R}\to \mathbb{R},$$

(1)

the objective is to determine $x^{*}\in \Omega$ such that $f\left(x^{*}\right)=0$. This abstract formulation arises in numerous contexts, including steady states of differential systems, discretized fractional partial differential equations, nonlinear inverse problems, and parameter identification models. In practice, the structural properties of f, such as smoothness, nonlocality, and memory dependence, strongly influence the performance and stability of numerical solvers designed to approximate its roots [1,2,3].

Classical calculus assumes locality: the derivative at a point depends solely on infinitesimal neighborhood information. However, many real-world systems exhibit hereditary effects, long-range interactions, and scaling symmetries that cannot be adequately represented within this local framework. Fractional calculus [4] generalizes differentiation to non-integer orders, thereby incorporating memory effects through integral operators. For example, the Caputo derivative of order $\alpha \in (0,1)$ introduces nonlocal dependence through a power-law kernel, allowing past states of a system to influence its present behavior. Such operators arise naturally in models of viscoelasticity, anomalous diffusion, bioheat transfer, and population dynamics with memory.

Parallel to fractional calculus, quantum calculus (or q-calculus) eliminates the classical limit process and replaces differentiation by finite scaling differences [5]. The Jackson q-derivative is defined as

$$D_{q}f\left(x\right)={\displaystyle \frac{f\left(qx\right)-f\left(x\right)}{(q-1)x}},q\ne 1,$$

(2)

and captures discrete self-similarity and multiplicative scaling effects. Fractional quantum operators combine these two extensions, enabling the simultaneous representation of memory and scaling phenomena [6,7]. From a theoretical perspective, these operators enlarge the functional space in which nonlinear operators act and provide additional structural flexibility for the construction of iterative numerical schemes.

Root-finding methods are typically formulated as fixed-point iterations,

$$x_{k+1}=G\left(x_k\right),$$

(3)

where $G:\Omega \to \Omega$ satisfies $G\left(x^{*}\right)=x^{*}$. The convergence of (3) is governed by the spectral radius of the derivative,

$$\rho \left(G^{\prime }\left(x^{*}\right)\right)<1.$$

Among classical approaches, Newton’s method

$$x_{k+1}=x_k-{\displaystyle \frac{f\left(x_k\right)}{f^{\prime }\left(x_k\right)}},$$

(4)

is widely used due to its quadratic convergence under standard regularity assumptions. However, its theoretical guarantees are essentially local, and the size of the basin of attraction depends strongly on the nonlinear structure of f. For highly nonlinear or stiff operators, the associated iterative map may exhibit multiple attracting cycles, chaotic regions, or narrow attraction domains [8].

From the perspective of complex dynamics [9], iterative schemes generate intricate fractal basin structures in the complex plane. Such analyses provide valuable insight into the global stability of iterative schemes beyond classical convergence order. Consequently, the convergence behavior of a method depends not only on its local order but also on the global dynamics of the iteration map G. This observation motivates the development of iterative frameworks with additional structural parameters capable of reshaping attraction domains and improving global convergence behavior. In this context, fractional and quantum calculus provide promising tools for enhancing both the flexibility and stability of iterative numerical schemes.

Building on these observations, several studies have extended classical iterative methods by incorporating fractional and q-calculus operators in order to improve convergence, stability, and global dynamical behavior. For instance, Ali et al. investigate iterative schemes based on fractional derivatives and report improved efficiency compared with classical root-finding methods [10]. Similarly, Sana et al. develop q-iterative techniques derived from q-Taylor expansions and coupled system formulations, demonstrating enhanced convergence properties [11]. Other q-analogue approaches have also been proposed, achieving higher convergence rates under suitable choices of the parameter q [12,13]. Moreover, multi-step and parallel q-fractional algorithms have been introduced to approximate multiple roots simultaneously while providing insights into the global convergence behavior of the associated iterative maps [14].

Despite these advances, a unified framework that systematically integrates fractional and quantum operators in order to control the geometry of attraction domains and the dynamical properties of iterative schemes remains largely unexplored. More specifically, while fractional Newton-type methods [15,16] and q-deformed schemes [17] have been investigated separately in the literature, several theoretical and computational questions remain open. In particular, the unified integration of fractional and quantum operators within a single iterative framework—and the resulting impact on convergence and stability properties—has not yet been systematically studied.

• The absence of a unified framework combining Caputo-type memory operators with quantum scaling derivatives within multi-step iterative structures.
• Limited theoretical understanding of how fractional and quantum parameters influence the dynamical system generated by the iteration map.
• Insufficient exploration of parameterized iterative families as mechanisms for controlling stability regions and convergence behavior.
• A lack of systematic investigation of the geometry of attraction basins under fractional–quantum perturbations.

Existing approaches [18,19,20] and references therein typically treat the fractional order $\alpha$ and the quantum parameter q as fixed constants, rather than as intrinsic regulators of the convergence dynamics. In contrast, we introduce a parameterized family of two-step fractional quantum iterative schemes in which these quantities, together with additional structural parameters, play an active role in shaping the iterative process. The proposed framework combines fractional and quantum operators within a flexible two-step structure that enables improved control of accuracy, stability, and convergence behavior.

The resulting formulation defines a family of iterative maps whose fixed points correspond to solutions of the nonlinear equation. This perspective provides a natural and rigorous basis for both convergence analysis and the investigation of fractal basin structures associated with the dynamical behavior of the proposed schemes. The proposed framework therefore provides a new perspective for the design and analysis of fractional–quantum iterative methods. The principal theoretical contributions of this work can be summarized as follows:

1.

Formulation of a unified fractional–quantum operator suitable for the construction of iterative root-finding schemes.

2.

Development of a parameterized two-step iterative framework that provides additional structural flexibility in controlling convergence behavior.

3.

Derivation of the local convergence properties through a perturbation analysis of the classical Taylor expansion in the presence of fractional operators.

4.

Dynamical systems analysis of the proposed schemes, including the study of stability and basin geometry in the complex plane.

5.

Identification of fractional and quantum parameters as geometric regulators capable of modifying the structure of attraction domains.

The novelty of this study lies in interpreting fractional and quantum operators not merely as modeling tools, but as intrinsic components of iterative dynamics. By embedding Caputo-type memory effects and q-scaling operators directly into the iterative structure, we construct a dynamical framework whose convergence behavior can be continuously tuned through the parameters $(\alpha ,q,{\beta }_1,{\beta }_2)$, where ${\beta }_1$ and ${\beta }_2$ are additional structural parameters utilized by the proposed schemes.

To the best of our knowledge, the systematic combination of the following elements has not yet been explored within a unified framework for solving nonlinear equations:

• Caputo-type fractional derivatives,
• quantum (q-) difference operators,
• parameterized two-step iterative schemes,
• and complex fractal basin analysis.

The paper is organized as follows. Section 2 introduces the fractional quantum operators and establishes their main analytical properties. Section 3 develops the proposed family of two-step iterative schemes and presents the corresponding convergence analysis. Section 4 investigates the stability and dynamical behavior of the methods, including the analysis of fractal basin structures. Section 5 provides numerical experiments on nonlinear problems arising in engineering and biomedical applications. Finally, Section 6 concludes the paper and outlines directions for future research.

2. Preliminaries

In this section, we recall the main concepts of $\mathfrak{q}$-calculus and fractional $\mathfrak{q}$-calculus that form the analytical foundation for the fractional $\mathfrak{q}$-iterative schemes developed later in the paper. Throughout the paper, we assume $\mathfrak{q}\in (0,1)$ and consider $0<a<x<b$. We denote by $\mathbb{N}=\{1,2,\dots \}$ the set of positive integers and by ${\mathbb{N}}_0=\mathbb{N}\cup \left\{0\right\}$ the set of non-negative integers.

Basic Elements of $\mathfrak{q}$-Calculus

The concept of $\mathfrak{q}$-analogues extends classical algebraic quantities by introducing a deformation parameter $\mathfrak{q}$. These generalized quantities reduce to their classical counterparts in the limit $\mathfrak{q}\to 1^{-}$.

Definition 1 

($\mathfrak{q}$-Number and $\mathfrak{q}$-Factorial [21]). The $\mathfrak{q}$-number associated with $r\in \mathbb{R}$ is defined by

$${\left[r\right]}_{\mathfrak{q}}={\displaystyle \frac{1-{\mathfrak{q}}^r}{1-\mathfrak{q}}}.$$

(5)

For $n\in \mathbb{N}$, the $\mathfrak{q}$-factorial is defined recursively as

$${\left[n\right]}_{\mathfrak{q}}!=\prod _{k=1}^n{\left[k\right]}_{\mathfrak{q}},{\left[0\right]}_{\mathfrak{q}}!:=1.$$

(6)

The $\mathfrak{q}$-shifted power plays a fundamental role in the construction of fractional kernels and related operators.

Definition 2 

($\mathfrak{q}$-Shifted Power [22]). For $a,b\in \mathbb{R}$ and $n\in \mathbb{N}$,

$${(a-b)}_{\mathfrak{q}}^n=\prod _{j=0}^{n-1}(a-b{\mathfrak{q}}^j).$$

(7)

More generally, for $\theta \in \mathbb{R}$,

$${(a-b)}_{\mathfrak{q}}^{\theta }=a^{\theta }\prod _{j=0}^{\infty }{\displaystyle \frac{a-b{\mathfrak{q}}^j}{a-b{\mathfrak{q}}^{j+\theta }}},$$

(8)

whenever the infinite product converges.

The central operator in quantum calculus is the Jackson $\mathfrak{q}$-derivative, which replaces the classical limit process by multiplicative scaling.

Definition 3 

(Jackson $\mathfrak{q}$-Derivative [23]). Let $f:[0,b]\to \mathbb{R}$. The Jackson $\mathfrak{q}$-derivative is defined by

$$\left(D_{\mathfrak{q}}f\right)\left(x\right)={\displaystyle \frac{f\left(\mathfrak{q}x\right)-f\left(x\right)}{(\mathfrak{q}-1)x}},x\ne 0,$$

(9)

with

$$\left(D_{\mathfrak{q}}f\right)\left(0\right)=\underset{x\to 0}{lim}\left(D_{\mathfrak{q}}f\right)\left(x\right),$$

(10)

provided the limit exists. Higher-order derivatives are defined recursively as

$$D_{\mathfrak{q}}^{n}f:=D_{\mathfrak{q}}\left(D_{\mathfrak{q}}^{n-1}f\right),n\in \mathbb{N}.$$

(11)

Jackson Integral and $\mathfrak{q}$-Gamma Function

The integral counterpart of the Jackson $\mathfrak{q}$-derivative is the Jackson $\mathfrak{q}$-integral, which restores an inversion relationship analogous to the fundamental theorem of classical calculus.

Definition 4 

(Jackson $\mathfrak{q}$-Integral [24]). The Jackson $\mathfrak{q}$-integral over $[0,x]$ is defined by

$${\int }_0^{x}f\left(t\right)d_{\mathfrak{q}}t=(1-\mathfrak{q})x\sum _{k=0}^{\infty }{\mathfrak{q}}^{k}f\left(x{\mathfrak{q}}^k\right),$$

(12)

whenever the series converges.

The following identity expresses the inverse relation between the $\mathfrak{q}$-derivative and the Jackson integral.

2.0.0.1. Proposition 2.1 (Fundamental Theorem of $\mathfrak{q}$-Calculus [25]).

If f is continuous at 0, then

$$D_{\mathfrak{q}}\left({\int }_0^{x}f\left(t\right)d_{\mathfrak{q}}t\right)=f\left(x\right).$$

(13)

To extend $\mathfrak{q}$-calculus operators to fractional orders, the $\mathfrak{q}$-Gamma function is introduced as a normalization factor.

Definition 5 

($\mathfrak{q}$-Gamma Function [26]). For $\alpha >0$, the $\mathfrak{q}$-Gamma function is defined by

$${\Gamma }_{\mathfrak{q}}\left(\alpha \right)={(1-\mathfrak{q})}^{1-\alpha }{\displaystyle \frac{{(\mathfrak{q};\mathfrak{q})}_{\infty }}{{({\mathfrak{q}}^{\alpha };\mathfrak{q})}_{\infty }}},$$

(14)

where

$${(a;\mathfrak{q})}_{\infty }=\prod _{k=0}^{\infty }(1-a{\mathfrak{q}}^k).$$

(15)

It satisfies the functional relation

$${\Gamma }_{\mathfrak{q}}(\alpha +1)={\left[\alpha \right]}_{\mathfrak{q}}{\Gamma }_{\mathfrak{q}}\left(\alpha \right).$$

(16)

Fractional $\mathfrak{q}$-Operators

Fractional $\mathfrak{q}$-operators extend the concepts of $\mathfrak{q}$-calculus to non-integer orders and provide the analytical framework required for the construction of fractional $\mathfrak{q}$-iterative schemes.

Definition 6 

(Riemann–Liouville Fractional $\mathfrak{q}$-Integral [27]). For $\alpha >0$, the Riemann–Liouville fractional $\mathfrak{q}$-integral of order α is defined as

$$\left(I_{\mathfrak{q},a}^{\left[\alpha \right]}f\right)\left(x\right)={\displaystyle \frac{1}{{\Gamma }_{\mathfrak{q}}\left(\alpha \right)}}{\int }_a^x{(x-t)}_{\mathfrak{q}}^{\alpha -1}f\left(t\right)d_{\mathfrak{q}}t.$$

(17)

Based on this operator, the corresponding Riemann–Liouville fractional $\mathfrak{q}$-derivative is defined as follows.

Definition 7 

(Riemann–Liouville Fractional $\mathfrak{q}$-Derivative [28]). Let $\alpha >0$ and $n=\lceil \alpha \rceil$. Then

$$\left(D_{\mathfrak{q},a}^{\left[\alpha \right]}f\right)\left(x\right)=D_{\mathfrak{q}}^n\left(I_{\mathfrak{q},a}^{n-\alpha }f\right)\left(x\right).$$

(18)

An alternative formulation frequently used in applications is the Caputo-type fractional $\mathfrak{q}$-derivative, which differs from the Riemann–Liouville derivative in the order of differentiation and integration.

Definition 8 

(Caputo Fractional $\mathfrak{q}$-Derivative [29]). Let $\alpha >0$ and $n=\lceil \alpha \rceil$. Then

$$\left({}^{C}D_{\mathfrak{q},a}^{\left[\alpha \right]}f\right)\left(x\right)=I_{\mathfrak{q},a}^{n-\alpha }\left(D_{\mathfrak{q}}^{n}f\right)\left(x\right).$$

(19)

Key Properties and Fractional $\mathfrak{q}$-Taylor Formula

The following identities describe the interaction between classical $\mathfrak{q}$-derivatives and fractional $\mathfrak{q}$-operators introduced in the previous subsection.

2.0.0.2. Proposition 2.2.

For $\alpha \in {\mathbb{R}}^{+}\setminus \mathbb{N}$ and $0<a<x<b$,

$${}^{C}D_{\mathfrak{q},a}^{\left[\alpha \right]}\left(D_{\mathfrak{q}}f\right)\left(x\right)={}^{C}D_{\mathfrak{q},a}^{[\alpha +1]}f\left(x\right),$$

(20)

and

$$D_{\mathfrak{q}}\left({}^{C}D_{\mathfrak{q},a}^{\left[\alpha \right]}f\right)\left(x\right)={}^{C}D_{\mathfrak{q},a}^{[\alpha +1]}f\left(x\right)+{\displaystyle \frac{D_{\mathfrak{q}}^{\lceil \alpha \rceil }f\left(a\right)}{{\Gamma }_{\mathfrak{q}}(\lceil \alpha \rceil -\alpha )}}{(x-a)}_{\mathfrak{q}}^{\lceil \alpha \rceil -\alpha -1}.$$

(21)

The fractional $\mathfrak{q}$-Taylor expansion provides the fundamental analytical tool for deriving local error equations in fractional $\mathfrak{q}$-iterative methods.

Theorem 2.1 

(Fractional $\mathfrak{q}$-Taylor Expansion [30]). Let f be sufficiently smooth and $\alpha >0$. Then

$$f\left(x\right)=\sum _{k=0}^{n-1}{\displaystyle \frac{{(x-a)}_{\mathfrak{q}}^k}{{\left[k\right]}_{\mathfrak{q}}!}}{D}_{\mathfrak{q}}^{k}f\left(a\right)+{\displaystyle \frac{1}{{\Gamma }_{\mathfrak{q}}\left(\alpha \right)}}{\int }_a^x{(x-t)}_{\mathfrak{q}}^{\alpha -1}\left({}^{C}D_{\mathfrak{q},a}^{\alpha }f\right)\left(t\right)d_{\mathfrak{q}}t,$$

(22)

where $n=\lceil \alpha \rceil$.

3. Development and Analysis of the $\mathfrak{q}$-Fractional Scheme

3.1. Proposed Fractional Weighted Iterative Method

Let $\xi$ be a simple root of $f:\mathbb{R}\to \mathbb{R}$, and define the error at the k-th iteration by

$$e_k:=x_k-\xi .$$

We now introduce the proposed two-step weighted iterative method, abbreviated as ${\mathrm{MSB}}_{\mathfrak{q}:\alpha }$, for the solution of nonlinear equations:

$$\left\{\begin{array}{cc}\hfill y_k& =x_k-{\left({\Gamma }_{\mathfrak{q}}(\alpha +1){\displaystyle \frac{f\left(x_k\right)}{\left[{}^{C}D_{\mathfrak{q}}^{\left[\alpha \right]}f\right]\left(x_k\right)}}\left({\displaystyle \frac{1}{1-{\displaystyle \frac{{\beta }_{1}f\left(x_k\right)}{1+{\beta }_{2}f\left(x_k\right)}}}}-f\left(x_k\right)\right)\right)}^{{\displaystyle \frac{1}{\alpha }}},\hfill \\ \hfill x_{k+1}& =y_k-{\left({\Gamma }_{\mathfrak{q}}(\alpha +1){\displaystyle \frac{f\left(x_k\right)}{\left[{}^{C}D_{\mathfrak{q}}^{\left[\alpha \right]}f\right]\left(x_k\right)}}\left({\displaystyle \frac{f\left(y_k\right)}{f\left(x_k\right)}}+2{\left({\displaystyle \frac{f\left(y_k\right)}{f\left(x_k\right)}}\right)}^2\right)\right)}^{{\displaystyle \frac{1}{\alpha }}}.\hfill \end{array}\right.$$

(23)

Here, ${\beta }_1,{\beta }_2\in \mathbb{R}$ are free parameters that regulate the convergence behavior of the scheme.

3.2. Convergence Theorem

The following theorem establishes the convergence order of the proposed iterative scheme.

Theorem 3.1. 

Let f be sufficiently smooth in a neighborhood of the simple root ξ, and let the initial approximation $x_0$ be sufficiently close to ξ. Then the sequence $\left\{x_k\right\}$ generated by the iterative scheme (23) satisfies

$$e_{k+1}=({\mathcal{D}}_1+{\mathcal{D}}_2+{\mathcal{D}}_3+{\mathcal{D}}_4+{\mathcal{D}}_5)\left(e_k^{3\alpha +1};\mathfrak{q}\right)+O\left(e_k^{4\alpha +1};\mathfrak{q}\right),$$

(24)

where the coefficients ${\mathcal{D}}_1,\dots ,{\mathcal{D}}_5$ depend on the Taylor coefficients of f at ξ and on the parameters α, ${\beta }_1$, and ${\beta }_2$. Consequently, the proposed method has convergence order $3\alpha +1$.

Proof. 

The proof is based on a series expansion of f around the simple root $\xi$ in terms of the error $e_k$. We write the expansion in a fully normalized form by factoring out the leading term:

$$f\left(x_k\right)={\displaystyle \frac{{}^{C}D_{\mathfrak{q}}f\left(\xi \right)}{{\Gamma }_{\mathfrak{q}}(\alpha +1)}}\left[\left(e_k^{\alpha };\mathfrak{q}\right)+c_2\left(e_k^{2\alpha };\mathfrak{q}\right)+c_3\left(e_k^{3\alpha };\mathfrak{q}\right)+c_4\left(e_k^{4\alpha };\mathfrak{q}\right)+O\left(e_k^{5\alpha };\mathfrak{q}\right)\right],$$

(25)

where

$$c_j={\displaystyle \frac{{\Gamma }_{\mathfrak{q}}(\alpha +1)}{{\Gamma }_{\mathfrak{q}}(j\alpha +1)}}{\displaystyle \frac{\left[{}^{C}D_{\mathfrak{q}}^{\left[j\right]}f\right]\left(\xi \right)}{j!\left[{}^{C}D_{\mathfrak{q}}^{\left[1\right]}f\right]\left(\xi \right)}},j\ge 2.$$

Next, we expand the Caputo $\mathfrak{q}$-derivative at the iterate $x_k$.

$$\begin{array}{cc}\hfill \left[{}^{C}D_{\mathfrak{q}}f\right]\left(x_k\right)& ={\displaystyle \frac{\left[{}^{C}D_{\mathfrak{q}}^{\left[1\right]}f\right]\left(\xi \right)}{{\Gamma }_{\mathfrak{q}}(\alpha +1)}}[{\Gamma }_{\mathfrak{q}}(\alpha +1)+{\displaystyle \frac{{\Gamma }_{\mathfrak{q}}(2\alpha +1)c_2\left(e_k^{\alpha };\mathfrak{q}\right)}{{\Gamma }_{\mathfrak{q}}(\alpha +1)}}\hfill \\ \hfill & +{\displaystyle \frac{{\Gamma }_{\mathfrak{q}}(3\alpha +1)c_3\left(e_k^{2\alpha };\mathfrak{q}\right)}{{\Gamma }_{\mathfrak{q}}(2\alpha +1)}}+O\left(e_k^{3\alpha };\mathfrak{q}\right)].\hfill \end{array}$$

(26)

Next, the reciprocal of the derivative is expanded as a power series in terms of the error $e_k$:

$$\begin{array}{cc}\hfill {\displaystyle \frac{1}{\left({}^{C}D_{\mathfrak{q}}f\right)\left(x_k\right)}}& ={\displaystyle \frac{1}{{\Gamma }_{\mathfrak{q}}(\alpha +1)}}-{\displaystyle \frac{{\Gamma }_{\mathfrak{q}}(2\alpha +1)c_2\left(e_k^{\alpha };\mathfrak{q}\right)}{{\Gamma }_{\mathfrak{q}}{(\alpha +1)}^3}}\hfill \\ \hfill & +{\displaystyle \frac{\left(e_k^{2\alpha };\mathfrak{q}\right)}{{\Gamma }_{\mathfrak{q}}(\alpha +1)}}\left(-{\displaystyle \frac{{\Gamma }_{\mathfrak{q}}(3\alpha +1)c_3}{{\Gamma }_{\mathfrak{q}}(\alpha +1){\Gamma }_{\mathfrak{q}}(2\alpha +1)}}+{\displaystyle \frac{{\Gamma }_{\mathfrak{q}}{(2\alpha +1)}^2{c}_2^2}{{\Gamma }_{\mathfrak{q}}{(\alpha +1)}^4}}\right)\hfill \\ \hfill & +{\displaystyle \frac{\left(e_k^{3\alpha };\mathfrak{q}\right)}{{\Gamma }_{\mathfrak{q}}(\alpha +1)}}\left({\displaystyle \frac{c_2{\Gamma }_{\mathfrak{q}}(3\alpha +1)c_3}{{\Gamma }_{\mathfrak{q}}{(\alpha +1)}^3}}-{\displaystyle \frac{{\mathcal{C}}_1}{{\Gamma }_{\mathfrak{q}}{(\alpha +1)}^6}}\right)+O\left(e_k^{4\alpha };\mathfrak{q}\right).\hfill \end{array}$$

(27)

where ${\mathcal{C}}_1=\left({\left({\Gamma }_{\mathfrak{q}}\left(2\alpha +1\right)\right)}^3{c_2}^2-{\Gamma }_{\mathfrak{q}}\left(3\alpha +1\right)c_3{\left({\Gamma }_{\mathfrak{q}}\left(\alpha +1\right)\right)}^3\right)c_2$. Using the previous expansion, the Newton-type correction term can be written as

$$\begin{array}{cc}\hfill {\displaystyle \frac{f\left(x_k\right)}{\left({}^{C}D_{\mathfrak{q}}f\right)\left(x_k\right)}}& ={\displaystyle \frac{\left(e_k^{\alpha };\mathfrak{q}\right)}{{\Gamma }_{\mathfrak{q}}(\alpha +1)}}\hfill \\ \hfill & +\left({\displaystyle \frac{c_2}{{\Gamma }_{\mathfrak{q}}(\alpha +1)}}-{\displaystyle \frac{2^{2\alpha }{c}_2}{\sqrt{\pi }{\Gamma }_{\mathfrak{q}}{(\alpha +1)}^2}}{\Gamma }_{\mathfrak{q}}\left(\alpha +\frac{1}{2}\right)\right)\left(e_k^{2\alpha };\mathfrak{q}\right)\hfill \\ \hfill & +({\displaystyle \frac{c_3}{{\Gamma }_{\mathfrak{q}}(\alpha +1)}}-{\displaystyle \frac{2^{2\alpha }{c}_2^2}{\sqrt{\pi }{\Gamma }_{\mathfrak{q}}{(\alpha +1)}^2}}{\Gamma }_{\mathfrak{q}}\left(\alpha +\frac{1}{2}\right)\hfill \\ \hfill & +{\displaystyle \frac{2^{4\alpha }{c}_2^2}{\pi {\Gamma }_{\mathfrak{q}}{(\alpha +1)}^3}}{\Gamma }_{\mathfrak{q}}{\left(\alpha +\frac{1}{2}\right)}^2\hfill \\ \hfill & -{\displaystyle \frac{3^{3\alpha }\sqrt{3}}{2\sqrt{\pi }{2}^{2\alpha }{\Gamma }_{\mathfrak{q}}{(\alpha +1)}^2}}{\Gamma }_{\mathfrak{q}}\left(\alpha +\frac{1}{3}\right){\Gamma }_{\mathfrak{q}}\left(\alpha +\frac{2}{3}\right))\left(e_k^{3\alpha };\mathfrak{q}\right)\hfill \\ \hfill & +O\left(e_k^{4\alpha };\mathfrak{q}\right).\hfill \end{array}$$

(28)

Intermediate iteration $y_k$. We now expand the first step of the proposed method.

$$\begin{array}{cc}\hfill {\left(1-{\displaystyle \frac{{\beta }_{1}f\left(x_k\right)}{1+{\beta }_{2}f\left(x_k\right)}}\right)}^{-1}& =1+{\beta }_1\left(e_k^{\alpha };\mathfrak{q}\right)\hfill \\ \hfill & +({\beta }_1^2-{\beta }_1{\beta }_2+{\beta }_1{c}_2)\left(e_k^{2\alpha };\mathfrak{q}\right)\hfill \\ \hfill & +\left({\beta }_1{\beta }_2^2-2{\beta }_1{\beta }_2{c}_2+{\beta }_1{c}_3-({\beta }_1^2-{\beta }_1{\beta }_2+{\beta }_1{c}_2){\beta }_1\right)\left(e_k^{3\alpha };\mathfrak{q}\right)\hfill \\ \hfill & +O\left(e_k^{4\alpha };\mathfrak{q}\right).\hfill \end{array}$$

(29)

From this expansion we obtain

$$\begin{array}{cc}\hfill G_k& =1+({\beta }_1-1)\left(e_k^{\alpha };\mathfrak{q}\right)+({\beta }_1^2-{\beta }_1{\beta }_2+{\beta }_1{c}_2-c_2)\left(e_k^{2\alpha };\mathfrak{q}\right)\hfill \\ \hfill & +({\beta }_1^3-2{\beta }_1^2{\beta }_2+2{\beta }_1^2{c}_2+{\beta }_1{\beta }_2^2-2{\beta }_1{\beta }_2{c}_2+{\beta }_1{c}_3-c_3)\left(e_k^{3\alpha };\mathfrak{q}\right)\hfill \\ \hfill & +O\left(e_k^{4\alpha };\mathfrak{q}\right),\hfill \end{array}$$

(30)

where

$$G_k=\left({\left(1-{\displaystyle \frac{{\beta }_{1}f\left(x_k\right)}{1+{\beta }_{2}f\left(x_k\right)}}\right)}^{-1}-f\left(x_k\right)\right).$$

For convenience, we introduce the notation

$$\begin{array}{cc}\hfill F_k& :={\displaystyle \frac{f\left(x_k\right)}{\left({}^{C}D_{\mathfrak{q}}f\right)\left(x_k\right)}}\left({\left(1-{\displaystyle \frac{{\beta }_{1}f\left(x_k\right)}{1+{\beta }_{2}f\left(x_k\right)}}\right)}^{-1}-f\left(x_k\right)\right),\hfill \end{array}$$

(31)

For notational convenience, we define

$$G_x:={\Gamma }_{\mathfrak{q}}\left(x\right),$$

and this convention is used consistently throughout the manuscript. In particular,

$$\begin{array}{cc}\hfill G_{\alpha +1}& ={\Gamma }_{\mathfrak{q}}(\alpha +1),G_{\alpha +\frac{1}{2}}={\Gamma }_{\mathfrak{q}}\left(\alpha +\frac{1}{2}\right),\hfill \\ \hfill G_{\alpha +\frac{1}{3}}& ={\Gamma }_{\mathfrak{q}}\left(\alpha +\frac{1}{3}\right),G_{\alpha +\frac{2}{3}}={\Gamma }_{\mathfrak{q}}\left(\alpha +\frac{2}{3}\right).\hfill \end{array}$$

(32)

Substituting the previous expansions yields

$$\begin{array}{cc}\hfill F_k& ={\displaystyle \frac{\left(e_k^{\alpha };\mathfrak{q}\right)}{G_{\alpha +1}}}+\left({\displaystyle \frac{{\beta }_1+c_2-1}{G_{\alpha +1}}}-{\displaystyle \frac{2^{2\alpha }{c}_2}{\sqrt{\pi }{G}_{\alpha +1}^2}}{G}_{\alpha +1/2}\right)\left(e_k^{2\alpha };\mathfrak{q}\right)\hfill \\ \hfill & +(-{\displaystyle \frac{2^{2\alpha }{c}_2{\beta }_1}{\sqrt{\pi }{G}_{\alpha +1}^2}}{G}_{\alpha +1/2}-{\displaystyle \frac{3^{3\alpha }\sqrt{3}{c}_3}{2\sqrt{\pi }{2}^{2\alpha }{G}_{\alpha +1}^2}}{\displaystyle \frac{G_{\alpha +1/3}{G}_{\alpha +2/3}}{G_{\alpha +1/2}}}\hfill \\ \hfill & -{\displaystyle \frac{2c_2}{G_{\alpha +1}}}+{\displaystyle \frac{{\beta }_1^2+c_3}{G_{\alpha +1}}}-{\displaystyle \frac{2^{2\alpha }{c}_2^2}{{\alpha }^2\sqrt{\pi }{G}_{\alpha +1}^2}}{G}_{\alpha +1/2}\hfill \\ \hfill & +{\displaystyle \frac{2^{4\alpha }{c}_2^2}{{\alpha }^3\pi G_{\alpha +1}^3}}{G}_{\alpha +1/2}^2+{\displaystyle \frac{2^{2\alpha }{c}_2}{\sqrt{\pi }{G}_{\alpha +1}^2}}{G}_{\alpha +1/2}+{\displaystyle \frac{2{\beta }_1{c}_2-{\beta }_1{\beta }_2}{G_{\alpha +1}}})\left(e_k^{3\alpha };\mathfrak{q}\right)\hfill \\ \hfill & +O\left(e_k^{4\alpha };\mathfrak{q}\right).\hfill \end{array}$$

(33)

From the previous expansion, we obtain

$$\begin{array}{cc}\hfill y_k-\xi & =x_k-\xi -{\left(G_{\alpha +1}{F}_k\right)}^{{\displaystyle \frac{1}{\alpha }}}=\left(-{\beta }_1-c_2+1+{\displaystyle \frac{2^{2\alpha }{c}_2}{G_{\alpha +1}\sqrt{\pi }}}{G}_{\alpha +1/2}\right)\left(e_k^{\alpha +1};\mathfrak{q}\right)\hfill \\ \hfill & +({\displaystyle \frac{2^{2\alpha }{c}_2{\beta }_1}{G_{\alpha +1}\sqrt{\pi }}}{G}_{\alpha +1/2}-2{\beta }_1{c}_2+{\displaystyle \frac{3^{3\alpha }\sqrt{3}{c}_3}{2G_{\alpha +1}\sqrt{\pi }{2}^{2\alpha }}}{\displaystyle \frac{G_{\alpha +1/3}{G}_{\alpha +2/3}}{G_{\alpha +1/2}}}\hfill \\ \hfill & +2c_2-{\beta }_1^2-c_3+{\beta }_1{\beta }_2+{\displaystyle \frac{2^{2\alpha }{c}_2^2}{G_{\alpha +1}\sqrt{\pi }}}{G}_{\alpha +1/2}-{\displaystyle \frac{2^{4\alpha }{c}_2^2}{G_{\alpha +1}^2\pi }}{G}_{\alpha +1/2}^2\hfill \\ \hfill & -{\displaystyle \frac{2^{2\alpha }{c}_2}{G_{\alpha +1}\sqrt{\pi }}}{G}_{\alpha +1/2})\left(e_k^{2\alpha +1};\mathfrak{q}\right)+O\left(e_k^{3\alpha +1};\mathfrak{q}\right).\hfill \end{array}$$

(34)

Hence,

$$\begin{array}{cc}\hfill y_k-\xi & =\left(-{\beta }_1-c_2+1+{\displaystyle \frac{2^{2\alpha }{c}_2}{G_{\alpha +1}\sqrt{\pi }}}{G}_{\alpha +1/2}\right)\left(e_k^{\alpha +1};\mathfrak{q}\right)\hfill \\ \hfill & +({\displaystyle \frac{2^{2\alpha }{c}_2{\beta }_1}{G_{\alpha +1}\sqrt{\pi }}}{G}_{\alpha +1/2}-2{\beta }_1{c}_2+{\displaystyle \frac{3^{3\alpha }\sqrt{3}{c}_3}{2G_{\alpha +1}\sqrt{\pi }{2}^{2\alpha }}}{\displaystyle \frac{G_{\alpha +1/3}{G}_{\alpha +2/3}}{G_{\alpha +1/2}}}\hfill \\ \hfill & +2c_2-{\beta }_1^2-c_3+{\beta }_1{\beta }_2+{\displaystyle \frac{2^{2\alpha }{c}_2^2}{G_{\alpha +1}\sqrt{\pi }}}{G}_{\alpha +1/2}-{\displaystyle \frac{2^{4\alpha }{c}_2^2}{G_{\alpha +1}^2\pi }}{G}_{\alpha +1/2}^2\hfill \\ \hfill & -{\displaystyle \frac{2^{2\alpha }{c}_2}{G_{\alpha +1}\sqrt{\pi }}}{G}_{\alpha +1/2})\left(e_k^{2\alpha +1};\mathfrak{q}\right)+O\left(e_k^{3\alpha +1};\mathfrak{q}\right).\hfill \end{array}$$

(35)

Substituting this expression into $f\left(y_k\right)$ yields

$$\begin{array}{cc}\hfill f\left(y_k\right)& =-{\displaystyle \frac{\left(e_k^{\alpha +1};\mathfrak{q}\right)}{G_{\alpha +1}\sqrt{\pi }}}\left(({\beta }_1+c_2-1)G_{\alpha +1}\sqrt{\pi }-2^{2\alpha }{G}_{\alpha +1/2}{c}_2\right)\hfill \\ \hfill & +{\displaystyle \frac{\left(e_k^{2\alpha +1};\mathfrak{q}\right)}{2G_{\alpha +1}^2{\pi }^{3/2}{2}^{2\alpha }}}[-2·2^{6\alpha }{G}_{\alpha +1/2}^3{c}_2^2\sqrt{\pi }+2·2^{4\alpha }{G}_{\alpha +1/2}^2{G}_{\alpha +1}\pi (c_2{\beta }_1+c_2^2-c_2)\hfill \\ \hfill & +c_3{3}^{3\alpha }\sqrt{3}{G}_{\alpha +1/3}{G}_{\alpha +2/3}{G}_{\alpha +1}\pi +2G_{\alpha +1/2}·2^{2\alpha }{G}_{\alpha +1}^2{\pi }^{3/2}\hfill \\ \hfill & ({\beta }_1{\beta }_2-{\beta }_1^2-2{\beta }_1{c}_2+2c_2-c_3)]G_{\alpha +1/2}^{-1}+O\left(e_k^{3\alpha +1};\mathfrak{q}\right).\hfill \end{array}$$

(36)

Consequently,

$$\begin{array}{cc}\hfill {\displaystyle \frac{f\left(y_k\right)}{f\left(x_k\right)}}& ={\mathcal{A}}_1\left(e_k^{2\alpha +1};\mathfrak{q}\right)+{\mathcal{A}}_2\left(e_k^{3\alpha +1};\mathfrak{q}\right)+O\left(e_k^{4\alpha +1};\mathfrak{q}\right).\hfill \end{array}$$

(37)

The coefficients ${\mathcal{A}}_i$ and ${\mathcal{D}}_i$, for $i=1,\dots ,5$, were derived and verified using symbolic computation in Maple 2020, ensuring their correctness. In particular, the coefficients ${\mathcal{A}}_1$ and ${\mathcal{A}}_2$ are given by

$$\begin{array}{cc}\hfill {\mathcal{A}}_1& =2{\beta }_1^2{\beta }_2-{\beta }_1{\beta }_2^2-2{\beta }_1{c}_3+{\displaystyle \frac{3^{3\alpha }\sqrt{3}{\beta }_1{c}_3}{2G_{\alpha +1}\sqrt{\pi }{2}^{2\alpha }}}{\displaystyle \frac{G_{\alpha +1/3}{G}_{\alpha +2/3}}{G_{\alpha +1/2}}}+{\displaystyle \frac{3^{3\alpha }\sqrt{3}{c}_2{c}_3}{2G_{\alpha +1}\sqrt{\pi }{2}^{2\alpha }}}{\displaystyle \frac{G_{\alpha +1/3}{G}_{\alpha +2/3}}{G_{\alpha +1/2}}}\hfill \\ \hfill & +{\displaystyle \frac{2^{2\alpha }{c}_2{c}_3}{G_{\alpha +1}\sqrt{\pi }}}{G}_{\alpha +1/2}+{\displaystyle \frac{2^{2\alpha }{c}_2{\beta }_1^2}{G_{\alpha +1}\sqrt{\pi }}}{G}_{\alpha +1/2}-{\displaystyle \frac{2^{4\alpha }{c}_2^2{\beta }_1}{G_{\alpha +1}^{2}2\pi }}{G}_{\alpha +1/2}^2-2{\beta }_1^2{c}_2+c_2^2{\beta }_1-{\displaystyle \frac{2^{2\alpha }{c}_2{\beta }_1{\beta }_2}{G_{\alpha +1}\sqrt{\pi }}}{G}_{\alpha +1/2}\hfill \\ \hfill & -{\displaystyle \frac{3^{3\alpha }\sqrt{3}{c}_3}{2G_{\alpha +1}\sqrt{\pi }{2}^{2\alpha }}}{\displaystyle \frac{G_{\alpha +1/3}{G}_{\alpha +2/3}}{G_{\alpha +1/2}}}+3{\beta }_1{\beta }_2{c}_2-2{\beta }_1{c}_2-{\displaystyle \frac{3^{3\alpha }\sqrt{3}{c}_2{c}_3}{G_{\alpha +1}^2\pi }}{G}_{\alpha +1/3}{G}_{\alpha +2/3}-c_2^2-{\beta }_1^3\hfill \\ \hfill & +{\displaystyle \frac{2^{4\alpha }{c}_2^2}{G_{\alpha +1}^2\pi }}{G}_{\alpha +1/2}^2-{\displaystyle \frac{2^{2\alpha +1}{c}_2^3}{G_{\alpha +1}\sqrt{\pi }}}{G}_{\alpha +1/2}+{\displaystyle \frac{2^{6\alpha }{c}_2^3}{G_{\alpha +1}^3{\pi }^{3/2}}}{G}_{\alpha +1/2}^3.\hfill \end{array}$$

$$\begin{array}{cc}\hfill {\mathcal{A}}_2& =-{\displaystyle \frac{3^{3\alpha }\sqrt{3}{c}_2^2{c}_3}{2G_{\alpha +1}\sqrt{\pi }{2}^{2\alpha }}}{\displaystyle \frac{G_{\alpha +1/3}{G}_{\alpha +2/3}}{G_{\alpha +1/2}}}+{\displaystyle \frac{3^{3\alpha }\sqrt{3}{c}_2{c}_3}{2G_{\alpha +1}\sqrt{\pi }{2}^{2\alpha }}}{\displaystyle \frac{G_{\alpha +1/3}{G}_{\alpha +2/3}}{G_{\alpha +1/2}}}-{\displaystyle \frac{3^{3\alpha }\sqrt{3}{\beta }_1{c}_2{c}_3}{2G_{\alpha +1}\sqrt{\pi }{2}^{2\alpha }}}{\displaystyle \frac{G_{\alpha +1/3}{G}_{\alpha +2/3}}{G_{\alpha +1/2}}}\hfill \\ \hfill & -{\displaystyle \frac{2^{2\alpha }{c}_2^2{\beta }_1^2}{G_{\alpha +1}\sqrt{\pi }}}{G}_{\alpha +1/2}-2c_2{c}_3+{\displaystyle \frac{2^{4\alpha }{c}_2^3{\beta }_1}{G_{\alpha +1}^2\pi }}{G}_{\alpha +1/2}^2-{\displaystyle \frac{2^{2\alpha }{c}_2^2{c}_3}{G_{\alpha +1}\sqrt{\pi }}}{G}_{\alpha +1/2}\hfill \\ \hfill & +{\beta }_1^3{c}_2+c_2{c}_4+2{\beta }_1^2{c}_2^2-c_2^3{\beta }_1+2{\beta }_1{c}_2^2+{\displaystyle \frac{2^{2\alpha }{c}_2^2{\beta }_1{\beta }_2}{G_{\alpha +1}\sqrt{\pi }}}{G}_{\alpha +1/2}\hfill \\ \hfill & +2{\beta }_1{c}_2{c}_3-2{\beta }_1^2{\beta }_2{c}_2-3{\beta }_1{\beta }_2{c}_2^2+{\beta }_1{\beta }_2^2{c}_2+{\displaystyle \frac{3^{3\alpha }\sqrt{3}{c}_2^2{c}_3}{G_{\alpha +1}^2\pi }}{G}_{\alpha +1/3}{G}_{\alpha +2/3}\hfill \\ \hfill & -c_2^4-c_2^2+c_2^3-{\displaystyle \frac{2^{6\alpha }{c}_2^4}{G_{\alpha +1}^3{\pi }^{3/2}}}{G}_{\alpha +1/2}^3-{\displaystyle \frac{2^{4\alpha }{c}_2^3}{G_{\alpha +1}^2\pi }}{G}_{\alpha +1/2}^2+{\displaystyle \frac{2^{2\alpha +1}{c}_2^4}{G_{\alpha +1}\sqrt{\pi }}}{G}_{\alpha +1/2}.\hfill \end{array}$$

Next, substituting the expansion of ${\displaystyle \frac{f\left(y_k\right)}{f\left(x_k\right)}}$ into the second correction term gives

$$\begin{array}{cc}\hfill \left({\displaystyle \frac{f\left(y_k\right)}{f\left(x_k\right)}}+2{\left({\displaystyle \frac{f\left(y_k\right)}{f\left(x_k\right)}}\right)}^2\right)& ={\mathcal{A}}_1\left(e_k^{2\alpha +1};\mathfrak{q}\right)+{\mathcal{A}}_2\left(e_k^{3\alpha +1};\mathfrak{q}\right)+\dots .\hfill \end{array}$$

(38)

Accordingly, we define

$$F_k^{*}={\displaystyle \frac{f\left(x_k\right)}{\left({}^{C}D_{\mathfrak{q}}f\right)\left(x_k\right)}}\left({\displaystyle \frac{f\left(y_k\right)}{f\left(x_k\right)}}+2{\left({\displaystyle \frac{f\left(y_k\right)}{f\left(x_k\right)}}\right)}^2\right),$$

and obtain the expansion

$$\begin{array}{cc}\hfill F_k^{*}& =2{\mathcal{A}}_3\left(e_k^{\alpha +1};\mathfrak{q}\right)+2{\mathcal{A}}_4\left(e_k^{2\alpha +1};\mathfrak{q}\right)+{\mathcal{A}}_5\left(e_k^{3\alpha +1};\mathfrak{q}\right)+O\left(e_k^{4\alpha +1};\mathfrak{q}\right).\hfill \end{array}$$

(39)

The coefficients ${\mathcal{A}}_3$, ${\mathcal{A}}_4$, and ${\mathcal{A}}_5$ are given by

$$\begin{array}{cc}\hfill {\mathcal{A}}_3& =-{\beta }_1-c_2+1+{\displaystyle \frac{2^{2\alpha }{G}_{\alpha +1/2}{c}_2}{G_{\alpha +1}\sqrt{\pi }}},\hfill \end{array}$$

$$\begin{array}{cc}\hfill {\mathcal{A}}_4& ={\displaystyle \frac{2^{2\alpha }{G}_{\alpha +1/2}{c}_2{\beta }_1}{G_{\alpha +1}\sqrt{\pi }}}-2{\beta }_1{c}_2+{\displaystyle \frac{1}{2}}{\displaystyle \frac{c_3{3}^{3\alpha }\sqrt{3}{G}_{\alpha +1/3}{G}_{\alpha +2/3}}{G_{\alpha +1}\sqrt{\pi }{2}^{2\alpha }{G}_{\alpha +1/2}}}\hfill \\ \hfill & +2c_2-{\beta }_1^2-c_3+{\beta }_1{\beta }_2+{\displaystyle \frac{2^{2\alpha }{G}_{\alpha +1/2}{c}_2^2}{G_{\alpha +1}\sqrt{\pi }}}-{\displaystyle \frac{2^{4\alpha }{G}_{\alpha +1/2}^2{c}_2^2}{G_{\alpha +1}^2\pi }}-{\displaystyle \frac{2^{2\alpha }{G}_{\alpha +1/2}{c}_2}{G_{\alpha +1}\sqrt{\pi }}}.\hfill \end{array}$$

$$\begin{array}{cc}\hfill {\mathcal{A}}_5& ={\displaystyle \frac{2^{2\alpha }{c}_2{c}_3}{\sqrt{\pi }{G}_{\alpha +1}^2}}{G}_{\alpha +1/2}-{\displaystyle \frac{2^{4\alpha }{c}_2^2{\beta }_1}{G_{\alpha +1}^3\pi }}{G}_{\alpha +1/2}^2-{\displaystyle \frac{{\beta }_1{\beta }_2^2}{G_{\alpha +1}}}+2{\displaystyle \frac{{\beta }_1^2{\beta }_2}{G_{\alpha +1}}}-2{\displaystyle \frac{{\beta }_1^2{c}_2}{G_{\alpha +1}}}+{\displaystyle \frac{c_2^2{\beta }_1}{G_{\alpha +1}}}\hfill \\ \hfill & +{\displaystyle \frac{2^{2\alpha }{c}_2{\beta }_1^2}{\sqrt{\pi }{G}_{\alpha +1}^2}}{G}_{\alpha +1/2}-2{\displaystyle \frac{{\beta }_1{c}_2}{G_{\alpha +1}}}-2{\displaystyle \frac{{\beta }_1{c}_3}{G_{\alpha +1}}}-{\displaystyle \frac{3^{3\alpha }\sqrt{3}{c}_3}{2\sqrt{\pi }{G}_{\alpha +1}^2{2}^{2\alpha }}}{\displaystyle \frac{G_{\alpha +1/3}{G}_{\alpha +2/3}}{G_{\alpha +1/2}}}\hfill \\ \hfill & +{\displaystyle \frac{3^{3\alpha }\sqrt{3}{c}_2{c}_3}{2\sqrt{\pi }{G}_{\alpha +1}^2{2}^{2\alpha }}}{\displaystyle \frac{G_{\alpha +1/3}{G}_{\alpha +2/3}}{G_{\alpha +1/2}}}+{\displaystyle \frac{3^{3\alpha }\sqrt{3}{\beta }_1{c}_3}{2\sqrt{\pi }{G}_{\alpha +1}^2{2}^{2\alpha }}}{\displaystyle \frac{G_{\alpha +1/3}{G}_{\alpha +2/3}}{G_{\alpha +1/2}}}+3{\displaystyle \frac{{\beta }_1{\beta }_2{c}_2}{G_{\alpha +1}}}\hfill \\ \hfill & -{\displaystyle \frac{2^{2\alpha }{c}_2{\beta }_1{\beta }_2}{\sqrt{\pi }{G}_{\alpha +1}^2}}{G}_{\alpha +1/2}-{\displaystyle \frac{c_4}{G_{\alpha +1}}}-{\displaystyle \frac{{\beta }_1^3}{G_{\alpha +1}}}+{\displaystyle \frac{c_2^3}{G_{\alpha +1}}}-2{\displaystyle \frac{2^{2\alpha }{G}_{\alpha +1/2}{c}_2^3}{{\alpha }^2\sqrt{\pi }{G}_{\alpha +1}^2}}+{\displaystyle \frac{2^{4\alpha }{c}_2^2}{G_{\alpha +1}^3\pi }}{G}_{\alpha +1/2}^2\hfill \\ \hfill & -{\displaystyle \frac{c_2^2}{G_{\alpha +1}}}+{\displaystyle \frac{c_2}{G_{\alpha +1}}}+2{\displaystyle \frac{c_3}{G_{\alpha +1}}}+{\displaystyle \frac{2^{6\alpha }{c}_2^3}{G_{\alpha +1}^4{\pi }^{3/2}}}{G}_{\alpha +1/2}^3-{\displaystyle \frac{3^{3\alpha }\sqrt{3}{c}_2{c}_3}{G_{\alpha +1}^3\pi }}{G}_{\alpha +1/3}{G}_{\alpha +2/3}.\hfill \end{array}$$

Finally, substituting these expansions into the second step of the iteration yields

$$\begin{array}{cc}\hfill x_{k+1}-\xi & =e_{k+1}=({\mathcal{D}}_1+{\mathcal{D}}_2+{\mathcal{D}}_3+{\mathcal{D}}_4+{\mathcal{D}}_5)\left(e_k^{3\alpha +1};\mathfrak{q}\right)+O\left(e_k^{4\alpha +1};\mathfrak{q}\right),\hfill \end{array}$$

(40)

where   

$$\begin{array}{cc}\hfill {\mathcal{D}}_1& :={\displaystyle \frac{2^{4\alpha }{c}_2^2{\beta }_1}{G_{\alpha +1}^3\pi }}{G}_{\alpha +1/2}^2-{\displaystyle \frac{2^{4\alpha }{c}_2^2{\beta }_1}{G_{\alpha +1}^2\pi }}{G}_{\alpha +1/2}^2+{\displaystyle \frac{3^{3\alpha }\sqrt{3}{c}_2{c}_3}{G_{\alpha +1}^3\pi }}{G}_{\alpha +1/3}{G}_{\alpha +2/3}-{\displaystyle \frac{3^{3\alpha }\sqrt{3}{c}_2{c}_3}{G_{\alpha +1}^2\pi }}\times \hfill \\ \hfill & G_{\alpha +1/3}{G}_{\alpha +2/3},\hfill \end{array}$$

$$\begin{array}{cc}\hfill {\mathcal{D}}_2& :={\displaystyle \frac{c_4+{\beta }_1^3-c_2}{G_{\alpha +1}}}+3{\beta }_1{\beta }_2{c}_2+2{\displaystyle \frac{2^{2\alpha }{G}_{\alpha +1/2}{c}_2^3}{G_{\alpha +1}^2\sqrt{\pi }}}-2{\displaystyle \frac{2^{2\alpha }{G}_{\alpha +1/2}{c}_2^2}{G_{\alpha +1}\sqrt{\pi }}}-{\displaystyle \frac{2^{6\alpha }{c}_2^3}{G_{\alpha +1}^4{\pi }^{3/2}}}{G}_{\alpha +1/2}^3\hfill \\ \hfill & -{\displaystyle \frac{2^{4\alpha }{c}_2^2}{G_{\alpha +1}^3\pi }}{G}_{\alpha +1/2}^2+{\displaystyle \frac{2^{6\alpha }{c}_2^3}{G_{\alpha +1}^3{\pi }^{3/2}}}{G}_{\alpha +1/2}^3-{\displaystyle \frac{2^{4\alpha }{c}_2^3}{G_{\alpha +1}^2\pi }}{G}_{\alpha +1/2}^2-{\displaystyle \frac{2^{2\alpha }{c}_2{\beta }_1^2}{G_{\alpha +1}^2\sqrt{\pi }}}{G}_{\alpha +1/2},\hfill \end{array}$$

$$\begin{array}{cc}\hfill {\mathcal{D}}_3& :=2{\displaystyle \frac{2^{2\alpha }{G}_{\alpha +1/2}{c}_2^2{\beta }_1}{G_{\alpha +1}\sqrt{\pi }}}+{\displaystyle \frac{2^{2\alpha }{c}_2{c}_3}{G_{\alpha +1}\sqrt{\pi }}}{G}_{\alpha +1/2}-{\displaystyle \frac{2^{2\alpha }{c}_2{c}_3}{G_{\alpha +1}^2\sqrt{\pi }}}{G}_{\alpha +1/2}-{\beta }_1^3-c_4-{\displaystyle \frac{c_2^3}{G_{\alpha +1}}}-2{\displaystyle \frac{c_3}{G_{\alpha +1}}},\hfill \end{array}$$

$$\begin{array}{cc}\hfill {\mathcal{D}}_4& :={\displaystyle \frac{c_2^2}{G_{\alpha +1}}}+c_2^2-3{\displaystyle \frac{{\beta }_1{\beta }_2{c}_2}{G_{\alpha +1}}}+2{\beta }_1^2{\beta }_2-2c_3{\beta }_1-3c_2{\beta }_1^2-{\beta }_1{\beta }_2^2-c_2^2{\beta }_1-2{\displaystyle \frac{{\beta }_1^2{\beta }_2}{G_{\alpha +1}}}+{\displaystyle \frac{{\beta }_1{\beta }_2^2}{G_{\alpha +1}}},\hfill \\ \hfill {\mathcal{D}}_5& :=2{\displaystyle \frac{c_3{\beta }_1}{G_{\alpha +1}}}+2{\displaystyle \frac{c_2{\beta }_1^2}{G_{\alpha +1}}}-{\displaystyle \frac{c_2^2{\beta }_1}{G_{\alpha +1}}}+2{\displaystyle \frac{c_2{\beta }_1}{G_{\alpha +1}}}+2c_3.\hfill \end{array}$$

This proves the stated error equation and therefore establishes convergence of order $3\alpha +1$.    □

4. Complex Dynamical Stability of the Proposed Scheme

The global behavior of iterative root-finding schemes can be analyzed through the study of their associated iteration functions in the complex plane [31,32]. In particular, the analysis of fixed points, strange fixed points, and critical points provides valuable insight into the stability and convergence properties of the method.

To analyze these properties for the proposed method, we consider the following two-step scheme:

$$\left\{\begin{array}{cc}\hfill y& =x-{\left({\Gamma }_{\mathfrak{q}}(\alpha +1){\displaystyle \frac{f\left(x\right)}{\left({}^{C}D_{\mathfrak{q}}f\right)\left(x\right)}}\left({\displaystyle \frac{1}{1-\frac{{\beta }_{1}f\left(x\right)}{1+{\beta }_{2}f\left(x\right)}}}-f\left(x\right)\right)\right)}^{{\displaystyle \frac{1}{\alpha }}},\hfill \\ \hfill R\left(x\right)& =y-{\left({\Gamma }_{\mathfrak{q}}(\alpha +1){\displaystyle \frac{f\left(x\right)}{\left({}^{C}D_{\mathfrak{q}}f\right)\left(x\right)}}\left({\displaystyle \frac{f\left(y\right)}{f\left(x\right)}}+2{\left({\displaystyle \frac{f\left(y\right)}{f\left(x\right)}}\right)}^2\right)\right)}^{{\displaystyle \frac{1}{\alpha }}}.\hfill \end{array}\right.$$

(41)

To study the associated complex dynamics, we consider the quadratic polynomial

$$f\left(z\right)=(z-a)(z-b),a,b\in \mathbb{C}.$$

Following the standard approach in complex dynamical analysis, we apply the Möbius transformation

$$x={\displaystyle \frac{z-a}{z-b}},$$

(42)

which maps the roots of the polynomial as follows:

$$a\mapsto 0,b\mapsto \infty .$$

Under this transformation, the proposed scheme induces a rational iteration function

$$x_{h+1}=R\left(x_h\right),$$

(43)

which determines the complex dynamical behavior of the scheme.

4.1. Rational Iteration Map

After simplification, the iteration function associated with the proposed scheme can be written as the rational map

$$R(x;{\beta }_1,{\beta }_2,\alpha ,\mathfrak{q})=-{\displaystyle \frac{x^{3}A(x;{\beta }_1,{\beta }_2,\alpha ,\mathfrak{q})B(x;{\beta }_1,{\beta }_2,\alpha ,\mathfrak{q})}{32C(x;{\beta }_1,{\beta }_2,\alpha ,\mathfrak{q})D(x;{\beta }_1,{\beta }_2,\alpha ,\mathfrak{q})}},$$

(44)

where $A(x;{\beta }_1,{\beta }_2,\alpha ,\mathfrak{q})$, $B(x;{\beta }_1,{\beta }_2,\alpha ,\mathfrak{q})$, $C(x;{\beta }_1,{\beta }_2,\alpha ,\mathfrak{q})$, and $D(x;{\beta }_1,{\beta }_2,\alpha ,\mathfrak{q})$ are polynomials in x depending on the parameters. The polynomial components appearing in the rational map (44) are provided explicitly in Appendix A.

For the purpose of dynamical analysis, and to obtain a more tractable form of the iteration map, we consider the classical case $\alpha =1$ and $\mathfrak{q}=1$. Under these assumptions, the iteration function reduces to a rational function, whose fully expanded form is reported in Appendix A.

4.2. Fixed Points

Fixed points of the iteration map are the solutions of

$$R\left(x\right)=x.$$

(45)

Substituting the rational map (44) into (45) yields the equation

$$x+{\displaystyle \frac{x^{3}A(x;{\beta }_1,{\beta }_2,1,1)B(x;{\beta }_1,{\beta }_2,1,1)}{32C(x;{\beta }_1,{\beta }_2,1,1)D(x;{\beta }_1,{\beta }_2,1,1)}}=0.$$

(46)

From the Möbius transformation (42), the roots of the original polynomial are mapped to

$$x=0,x=\infty .$$

These points are fixed points of the iteration function and correspond to the roots of the original polynomial. They are therefore referred to as root fixed points.

All other solutions of equation (45) are referred to as strange fixed points, since they do not correspond to roots of the underlying nonlinear equation but instead arise from the intrinsic dynamical structure of the iteration map.

4.3. Stability of Fixed Points

The stability of a fixed point $x^{*}$ is determined by the multiplier

$$\lambda =R^{\prime }\left(x^{*}\right).$$

(47)

According to the modulus of the multiplier, fixed points are classified as follows:

$$\left\{\begin{array}{cc}\left|\lambda \right|<1,\hfill & \mathrm{attracting},\hfill \\ \left|\lambda \right|>1,\hfill & \mathrm{repelling},\hfill \\ \left|\lambda \right|=1,\hfill & \mathrm{neutral}.\hfill \end{array}\right.$$

Differentiating the rational map (44) using the quotient rule yields

$$R^{\prime }\left(x\right)=-{\displaystyle \frac{1}{32}}{\displaystyle \frac{{\mathcal{N}}_1-{\mathcal{N}}_2}{C{(x;{\beta }_1,{\beta }_2,1,1)}^{2}D{(x;{\beta }_1,{\beta }_2,1,1)}^2}}.$$

(48)

where

$${\mathcal{N}}_1={\left(x^{3}A(x;{\beta }_1,{\beta }_2,1,1)B(x;{\beta }_1,{\beta }_2,1,1)\right)}^{\prime }C(x;{\beta }_1,{\beta }_2,1,1)D(x;{\beta }_1,{\beta }_2,1,1),$$

and

$${\mathcal{N}}_2=x^{3}A(x;{\beta }_1,{\beta }_2,1,1)B(x;{\beta }_1,{\beta }_2,1,1){\left(C(x;{\beta }_1,{\beta }_2,1,1)D(x;{\beta }_1,{\beta }_2,1,1)\right)}^{\prime }.$$

Evaluating the derivative at the root fixed point $x=0$ gives

$$R^{\prime }\left(0\right)=0,$$

which shows that the root is a superattracting fixed point. This confirms the strong local convergence of the method toward the root.

4.4. Critical Points

Critical points play a central role in determining the global dynamics of the iteration map. They are obtained as the solutions of

$$R^{\prime }(x;{\beta }_1,{\beta }_2,1,1)=0.$$

(49)

Using the derivative expression (48), the critical points satisfy

$$\begin{array}{cc}\hfill & {\left(x^{3}A(x;{\beta }_1,{\beta }_2,1,1)B(x;{\beta }_1,{\beta }_2,1,1)\right)}^{\prime }C(x;{\beta }_1,{\beta }_2,1,1)D(x;{\beta }_1,{\beta }_2,1,1)\hfill \\ \hfill & -x^{3}A(x;{\beta }_1,{\beta }_2,1,1)B(x;{\beta }_1,{\beta }_2,1,1){\left(C(x;{\beta }_1,{\beta }_2,1,1)D(x;{\beta }_1,{\beta }_2,1,1)\right)}^{\prime }=0.\hfill \end{array}$$

(50)

Among these solutions, the points

$$x=0,x=\infty$$

are called fixed critical points, since they coincide with the root fixed points of the iteration map. All other solutions are referred to as free critical points. The orbits of the free critical points play a crucial role in determining the global convergence behavior and the structure of the basins of attraction.

For illustration, we consider the parameter choice ${\beta }_1=1$. Under this assumption, the critical points of the iteration map can be explicitly characterized as

$$\begin{array}{cc}\hfill {\mathrm{cr}}_1:x& =0(\mathrm{multiplicity}3),\hfill \\ \hfill {\mathrm{cr}}_2:x& =-1(\mathrm{multiplicity}4),\hfill \\ \hfill {\mathrm{cr}}_3:x& =1(\mathrm{multiplicity}8),\hfill \\ \hfill {\mathrm{cr}}_4:x& =RootOf\left(P(Z;{\beta }_1,{\beta }_2)\right),\hfill \\ \hfill {\mathrm{cr}}_5:x& =-2{\beta }_2+3\pm 2\sqrt{{\beta }_2^2-3{\beta }_2+2}(\mathrm{each}\mathrm{of}\mathrm{multiplicity}3),\hfill \end{array}$$

where $P(Z;{\beta }_1,{\beta }_2)$ is a polynomial of degree 18 given in Appendix A.

This connection will be illustrated in the parameter planes shown in Figure 1.

4.5. Parameter Planes

While dynamical planes describe the behavior of the method for fixed parameter values, parameter planes analyze the stability of the method with respect to variations in the parameters.

Let

$$({\beta }_1,{\beta }_2)\in {\mathbb{C}}^2.$$

For each parameter pair, the orbit of a selected free critical point c is computed as

$$c,R\left(c\right),R^2\left(c\right),\dots .$$

The parameter plane is defined as

$$\mathcal{P}=\left\{({\beta }_1,{\beta }_2):\underset{h\to \infty }{lim}{R}^{\left(h\right)}\left(c\right)={\xi }_i,i=1,2\right\},$$

(51)

where ${\xi }_1$ and ${\xi }_2$ denote the two roots of the transformed quadratic polynomial.

Stable parameter values produce dynamical planes with well-defined basins of attraction, whereas unstable parameter choices lead to chaotic or fractal structures in the complex plane.

Table 1 summarizes the main dynamical characteristics observed in the parameter planes associated with the critical points ${\mathrm{cr}}_1$ and ${\mathrm{cr}}_2$. The results indicate a tendency toward larger basin coverage and improved convergence behavior for the parameter plane corresponding to ${\mathrm{cr}}_2$, suggesting that the orbit of this critical point leads to a wider range of stable parameter values.

Furthermore, the parameter planes shown in Figure 1a,b illustrate how the stability of the method depends on the choice of the parameters. In particular, the regions associated with the critical point ${\mathrm{cr}}_2$ exhibit slightly larger areas of convergence, suggesting a wider range of stable parameter values compared to ${\mathrm{cr}}_1$.

4.6. Dynamical Planes

For fixed parameter values $({\beta }_1,{\beta }_2)$, the iteration function

$$x_{k+1}=R\left(x_k\right)$$

defines a discrete dynamical system on the complex plane.

Let $\xi$ denote a root of the nonlinear equation. The basin of attraction of $\xi$ is defined as

$$\mathcal{A}\left(\xi \right)=\left\{x_0\in \mathbb{C}:\underset{h\to \infty }{lim}{R}^{\left(h\right)}\left(x_0\right)=\xi \right\}.$$

(52)

The dynamical plane is obtained by coloring each initial point $x_0$ in the complex plane according to the root to which the iteration converges. Points converging to different roots are assigned different colors, while points that do not converge within a prescribed number of iterations are represented by a separate color.

In this study, the parameter pair $({\beta }_1,{\beta }_2)$ is selected from the stable regions identified in the parameter-plane analysis. Under these choices, the iteration map satisfies the stability condition derived from the derivative criterion (48). Consequently, the generated dynamical planes clearly illustrate the basins of attraction associated with the roots of the nonlinear equation.

Figure 2a–e presents the dynamical planes of the proposed scheme for several values of the fractional parameter $\alpha$. It can be observed that, as $\alpha$ increases, the basins of attraction associated with the roots become larger and more connected. At the same time, the fractal boundaries separating different basins gradually decrease, indicating a reduction in chaotic behavior and an improvement in the stability of the iterative process.

This behavior suggests that larger values of $\alpha$ are associated with more robust convergence properties of the proposed fractional $\mathfrak{q}$-scheme. In particular, the attraction regions become smoother and occupy larger portions of the complex plane, implying that a broader range of initial guesses converges successfully to the desired roots.

Table 2 provides quantitative measures supporting the graphical observations obtained from the dynamical planes. As the fractional parameter $\alpha$ increases, the percentage of convergent initial guesses grows significantly, indicating that larger portions of the complex plane belong to stable basins of attraction. At the same time, the average number of iterations decreases, demonstrating faster convergence of the proposed iterative method.

The computational time also decreases slightly because fewer iterations are required to reach the prescribed tolerance, while the memory consumption remains nearly constant. This shows that the improvement in convergence efficiency does not introduce additional computational overhead.

These results highlight a strong correspondence between the parameter-plane and dynamical-plane analyses. The stable regions identified in the parameter planes guide the selection of suitable parameter values, which in turn produce dynamical planes with well-structured basins of attraction and reduced chaotic behavior. These graphical representations provide a clear description of the stability and convergence properties of the proposed iterative scheme. For completeness and reproducibility, the computational procedure used to construct both the dynamical and parameter planes is summarized in Algorithm 1.

Algorithm 1:Construction of Dynamical and Parameter Planes for the ${\mathrm{MSB}}_{\mathfrak{q}:\alpha }$ Scheme

Require:  Nonlinear function f, Caputo $\mathfrak{q}$-derivative $\left({}^{C}D_{\mathfrak{q}}f\right)$, fractional parameter $\alpha$, quantum parameter $\mathfrak{q}$, grid $\Omega \subset \mathbb{C}$, parameter pair $({\beta }_1,{\beta }_2)$, maximum number of iterations $N_{max}$, tolerance $\epsilon$ Ensure:  Dynamical-plane and parameter-plane visualizations 1: Stage I: Dynamical Plane Generation 2: for all initial points $x_0\in \Omega$ do 3:     Set $k=0$ 4:     while $k<N_{max}$ do 5:         Compute first update $$y_k=x_k-{\left({\Gamma }_{\mathfrak{q}}(\alpha +1){\displaystyle \frac{f\left(x_k\right)}{\left({}^{C}D_{\mathfrak{q}}f\right)\left(x_k\right)}}\left({\displaystyle \frac{1}{1-{\displaystyle \frac{{\beta }_{1}f\left(x_k\right)}{1+{\beta }_{2}f\left(x_k\right)}}}}-f\left(x_k\right)\right)\right)}^{{\displaystyle \frac{1}{\alpha }}}$$ 6:         Compute second update $$x_{k+1}=y_k-{\left({\Gamma }_{\mathfrak{q}}(\alpha +1){\displaystyle \frac{f\left(x_k\right)}{\left({}^{C}D_{\mathfrak{q}}f\right)\left(x_k\right)}}\left({\displaystyle \frac{f\left(y_k\right)}{f\left(x_k\right)}}+2{\left({\displaystyle \frac{f\left(y_k\right)}{f\left(x_k\right)}}\right)}^2\right)\right)}^{{\displaystyle \frac{1}{\alpha }}}$$ 7:         if $|x_{k+1}-x_k|<\epsilon$or$\left|f\right(x_{k+1}\left)\right|<\epsilon$ then 8:            Identify the root $r_k$ reached by the orbit (e.g., by comparison with known roots) 9:            Assign to $x_0$ the color associated with root $r_k$ 10:            break 11:         end if 12:         $k\leftarrow k+1$ 13:     end while 14:     if $k=N_{max}$ then 15:         Mark $x_0$ as non-convergent 16:     end if 17: end for 18: Stage II: Parameter Plane Generation 19: Fix the reduced case $\alpha =1$ and $\mathfrak{q}=1$ 20: for all parameter pairs $({\beta }_1,{\beta }_2)\in \mathcal{P}$ do 21:     Select one free critical point c (e.g., ${\mathrm{cr}}_1$ or ${\mathrm{cr}}_2$) 22:     Set $x_0=c$, $k=0$ 23:     while $k<N_{max}$ do 24:         Compute $x_{k+1}=R\left(x_k\right)$ using the reduced map (44) 25:         if $|x_{k+1}-x_k|<\epsilon$or$\left|f\right(x_{k+1}\left)\right|<\epsilon$ then 26:            Identify the attracting root 27:            Assign to $({\beta }_1,{\beta }_2)$ the corresponding color 28:            break 29:         end if 30:         $k\leftarrow k+1$ 31:     end while 32:     if $k=N_{max}$ then 33:         Mark $({\beta }_1,{\beta }_2)$ as divergent 34:     end if 35: end for

5. Numerical Experiments and Performance Analysis

In this section, we present a comprehensive numerical study of the proposed Caputo $\mathfrak{q}$-fractional iterative scheme (41). The numerical experiments are designed with two main objectives: first, to validate the theoretical convergence properties established in the previous sections, and second, to compare the computational performance of the proposed method with several existing $\mathfrak{q}$-type iterative schemes. All experiments are carried out using high-precision arithmetic in order to accurately capture the asymptotic error behavior and minimize round-off effects.

5.1. Computational Environment

All numerical simulations were performed in MATLAB R2024a using the Symbolic Math Toolbox with variable-precision arithmetic (VPA). The working precision was set to

$$\mathtt{digits}\left(350\right),$$

which provides numerical accuracy of approximately ${10}^{-300}$.

The computations were carried out on a workstation with the following specifications:

• Processor: Intel Core i7 (12th Generation) @ 2.40 GHz
• RAM: 16 GB DDR4
• Operating System: Windows 11 (64-bit)
• Software: MATLAB R2024a with Symbolic Math Toolbox

The use of high-precision arithmetic ensures that the numerical errors observed during the iterative process are primarily due to the convergence behavior of the methods rather than floating-point limitations.

For reproducibility, all numerical codes were implemented using symbolic evaluation of the nonlinear functions and their derivatives under VPA arithmetic, thereby ensuring uniform precision throughout the computations. The overall computational procedure follows Algorithm 2 and the flowchart shown in Figure 3.

Algorithm 2:High-Precision Root Finding using the Proposed ${\mathrm{MSB}}_{\mathfrak{q}:\alpha }$ Method

Require:  Nonlinear function $f\left(x\right)$, Caputo $\mathfrak{q}$-derivative $\left({}^{C}D_{\mathfrak{q}}^{\left[\alpha \right]}f\right)\left(x\right)$, initial guesses ${\left\{x_0^{\left(i\right)}\right\}}_{i=1}^5$, parameters $({\beta }_1,{\beta }_2)$, maximum number of iterations $N_{max}$, tolerance $\epsilon ={10}^{-300}$ Ensure:  Error convergence history for each initial guess 1: Compute all roots ${\xi }_j$ of $f\left(x\right)$ using a high-precision solver 2: for$k=0$ to $N_{max}-1$ do 3:     for $i=1$ to 5 do 4:         Evaluate $$f\left(x_k^{\left(i\right)}\right),\left({}^{C}D_{\mathfrak{q}}^{\left[\alpha \right]}f\right)\left(x_k^{\left(i\right)}\right)$$ 5:         Compute the intermediate iterate $$y_k^{\left(i\right)}=x_k^{\left(i\right)}-{\left({\Gamma }_{\mathfrak{q}}(\alpha +1){\displaystyle \frac{f\left(x_k^{\left(i\right)}\right)}{\left({}^{C}D_{\mathfrak{q}}^{\left[\alpha \right]}f\right)\left(x_k^{\left(i\right)}\right)}}\left({\left(1-{\displaystyle \frac{{\beta }_{1}f\left(x_k^{\left(i\right)}\right)}{1+{\beta }_{2}f\left(x_k^{\left(i\right)}\right)}}\right)}^{-1}-f\left(x_k^{\left(i\right)}\right)\right)\right)}^{{\displaystyle \frac{1}{\alpha }}}$$ 6:         Compute the next iterate $$x_{k+1}^{\left(i\right)}=y_k^{\left(i\right)}-{\left({\Gamma }_{\mathfrak{q}}(\alpha +1){\displaystyle \frac{f\left(x_k^{\left(i\right)}\right)}{\left({}^{C}D_{\mathfrak{q}}^{\left[\alpha \right]}f\right)\left(x_k^{\left(i\right)}\right)}}\left({\displaystyle \frac{f\left(y_k^{\left(i\right)}\right)}{f\left(x_k^{\left(i\right)}\right)}}+2{\left({\displaystyle \frac{f\left(y_k^{\left(i\right)}\right)}{f\left(x_k^{\left(i\right)}\right)}}\right)}^2\right)\right)}^{{\displaystyle \frac{1}{\alpha }}}$$ 7:         Compute the error with respect to the nearest root $$e_{k+1}^{\left(i\right)}=\underset{j}{min}\left|x_{k+1}^{\left(i\right)}-{\xi }_j\right|$$ 8:         Store $e_{k+1}^{\left(i\right)}$ 9:     end for 10:     if $e_{k+1}^{\left(i\right)}<\epsilon$ for all $i=1,\dots ,5$ then 11:         Stop iteration 12:     end if 13: end for 14: Plot the error curves using a semilogarithmic scale

5.2. Performance Metrics

To provide a comprehensive comparison between the proposed scheme and existing methods, several performance indicators are considered.

5.2.1. Absolute Error

The absolute error at iteration k is defined as

$$e_k=|x_k-\xi |,$$

(53)

where $\xi$ denotes the exact root of the nonlinear equation.

5.2.2. Computational Order of Convergence (COC)

The computational order of convergence is estimated using

$$\mathrm{COC}\approx {\displaystyle \frac{ln\left({\displaystyle \frac{e_{k+1}}{e_k}}\right)}{ln\left({\displaystyle \frac{e_k}{e_{k-1}}}\right)}}.$$

(54)

This quantity provides a numerical estimate of the theoretical order of convergence.

5.2.3. Residual Error

The residual error is defined as

$$R_k=\left|f\left(x_k\right)\right|,$$

(55)

which measures how accurately the approximate solution satisfies the nonlinear equation $f\left(x\right)=0$.

5.2.4. Number of Iterations

The total number of iterations required to achieve the prescribed tolerance

$$e_k<{10}^{-300}$$

is recorded for each method.

5.3. Methods Used for Comparison

To evaluate the efficiency of the proposed scheme, we compare it with several existing $\mathfrak{q}$–type iterative methods available in the literature, including

• Sana et al. [11] — third- and fourth-order methods, abbreviated as ${\mathrm{SAN}}_{\mathfrak{q}:3}$ and ${\mathrm{SAN}}_{\mathfrak{q}:4}$, respectively;
• Shams et al. [14] — a fourth-order method, abbreviated as ${\mathrm{QSB}}_{\mathfrak{q}:4}$.

Each method is implemented using the same computational precision and stopping criteria in order to ensure a fair and consistent comparison.

5.4. Stopping Criteria

The iterative process is terminated when one of the following conditions is satisfied:

1.

The difference between two successive iterates satisfies

$$|x_{k+1}-x_k|<{10}^{-300}.$$

2.

The maximum number of iterations $N_{max}=50$ is reached.

To illustrate the convergence behavior of the considered methods, the absolute error is plotted against the iteration number on a semi–logarithmic scale.

In particular, for each initial guess $x_0$ the sequence of errors $e_k$ is recorded and represented using the semilogy function. This representation is particularly useful for observing the asymptotic convergence phase.

5.5. Application: Nonlinear Continuous Stirred Tank Reactor [33]

Highly nonlinear algebraic equations frequently arise in chemical engineering, particularly in the steady–state analysis of chemical reactors. One of the most widely studied systems is the continuous stirred tank reactor (CSTR), where nonlinearities emerge due to the Arrhenius temperature dependence of the reaction rate and the coupling between heat generation and heat removal mechanisms. Determining the steady–state reactor temperature is crucial for safe operation, since multiple equilibrium states may exist depending on the operating conditions.

Consider an irreversible first–order exothermic reaction occurring in a CSTR. The steady–state energy balance leads to the nonlinear equation

$$f\left(T\right)={\displaystyle \frac{F\rho C_p(T-T_f)}{V}}+{\displaystyle \frac{UA(T-T_c)}{V}}-(-\Delta H)k_0{e}^{-E/\left(RT\right)}{C}_A=0,$$

(56)

where T denotes the reactor temperature. The parameters represent physical properties of the reactor and reaction system. In particular:

• F : volumetric flow rate
• $\rho$ : fluid density
• $C_p$ : heat capacity
• $T_f$ : feed temperature
• $T_c$ : coolant temperature
• U : overall heat transfer coefficient
• A : heat transfer area
• V : reactor volume
• $\Delta H$ : heat of reaction
• $k_0$ : pre–exponential reaction constant
• E : activation energy
• R : universal gas constant
• $C_A$ : reactant concentration

For the numerical experiments, the following representative parameter values are used:

$$F=1{\mathrm{m}}^3/\mathrm{s},\rho =1000{\mathrm{kg}/\mathrm{m}}^3,C_p=4.2\mathrm{kJ}/\mathrm{kgK},$$

$$T_f=350\mathrm{K},T_c=300\mathrm{K},U=500{\mathrm{W}/\mathrm{m}}^2\mathrm{K},$$

$$A=5{\mathrm{m}}^2,V=1{\mathrm{m}}^3,k_0=7.2\times {10}^{10},$$

$$E=8.314\times {10}^4\mathrm{J}/\mathrm{mol},C_A=1\mathrm{mol}/\mathrm{L},\Delta H=-5\times {10}^4\mathrm{J}/\mathrm{mol}.$$

The steady–state temperature of the reactor is obtained by solving the nonlinear equation $f\left(T\right)=0$. Due to the strong exponential term appearing in the Arrhenius reaction rate, the equation may admit multiple solutions corresponding to different thermal equilibrium states.

Substituting the parameter values into the reactor energy balance equation yields

$$f\left(T\right)=4200(T-350)+2500(T-300)-3.6\times {10}^{15}{e}^{-10000/T}=0.$$

(57)

Combining the linear terms leads to the simplified nonlinear equation

$$f\left(T\right)=6700T-2.22\times {10}^6-3.6\times {10}^{15}{e}^{-10000/T}=0.$$

(58)

Because of the strong exponential nonlinearity arising from the Arrhenius reaction term, the equation does not admit a closed–form analytical solution. Therefore, numerical root–finding methods must be employed to determine the steady–state reactor temperature.

Using high–precision numerical computation, the solution of the nonlinear equation is obtained as

$$T^{*}\approx 5.373134322\times {10}^{11}\mathrm{K}.$$

This value corresponds to the steady–state solution of the nonlinear equation under the specified parameter configuration and is used as the reference root in the subsequent numerical experiments.

5.5.1. Efficiency Analysis

The efficiency of the considered schemes is evaluated using the absolute error obtained across successive iterations, as reported in Table 3. The proposed method ${\mathrm{MSB}}_{\mathfrak{q}:\alpha }$ demonstrates superior accuracy compared with the existing schemes ${\mathrm{SAN}}_{\mathfrak{q}:3}$, ${\mathrm{SAN}}_{\mathfrak{q}:4}$, and ${\mathrm{QSB}}_{\mathfrak{q}:4}$. In particular, the proposed method achieves extremely small residual errors within a limited number of iterations, indicating a faster asymptotic convergence behavior relative to the competing methods.

5.5.2. Stability and Computational Cost

The stability and computational efficiency of the iterative schemes are evaluated using several performance indicators, including the maximum error, CPU time, memory consumption, and the number of arithmetic operations. The results reported in Table 4 indicate that the proposed method exhibits improved numerical stability while maintaining competitive computational cost compared with the existing schemes.

5.5.3. Consistency for Fractional and Quantum Parameters

To examine the robustness of the proposed method, different values of the fractional parameter $\alpha$ and the quantum parameter $\mathfrak{q}$ are considered. To verify the consistency of the proposed scheme, we consider the initial vector

$${\mathbf{x}}_0=\mathrm{vpa}\left([15.37\times {10}^5,11.37\times {10}^3,41.37\times {10}^8,6.3732\times {10}^9,7.37\times {10}^9]\right).$$

The outcomes are presented in Table 5 and Figure 4. As illustrated in Figure 4, ${\mathrm{MSB}}_{q:\alpha }$ provides improved residual accuracy and overall performance relative to ${\mathrm{SAN}}_{\mathfrak{q}:3}$–${\mathrm{SAN}}_{\mathfrak{q}:4}$ and ${\mathrm{QSB}}_{\mathfrak{q}:4}$.

5.5.4. Dynamical and Computational Performance Comparison

The results presented in Table 6 and Figure 5 summarize the dynamical and computational performance of the considered iterative schemes. The comparison is based on the resolution of the fractal basin of attraction, the percentage of convergent and divergent points, the elapsed computational time, the number of arithmetic operations, CPU time, and memory consumption. The fractal basin plots illustrate the global convergence behavior of each method in the complex plane.

A larger convergence region and fewer divergence points indicate greater numerical stability and robustness of the iterative scheme. The column F-Resolution qualitatively describes the visual resolution of the fractal basin of attraction obtained from the dynamical plane analysis. From Table 6, it can be observed that the proposed Caputo fractional quantum scheme ${\mathrm{MSB}}_{\mathfrak{q}:\alpha }$ achieves the highest convergence ratio (96%) and the smallest percentage of divergence points (4%). In addition, the proposed scheme requires fewer arithmetic operations and exhibits lower CPU time and memory consumption compared with the competing methods. Consequently, the proposed ${\mathrm{MSB}}_{\mathfrak{q}:\alpha }$ scheme demonstrates superior overall performance in terms of convergence behavior and computational efficiency when compared with ${\mathrm{SAN}}_{\mathfrak{q}:3}$, ${\mathrm{SAN}}_{\mathfrak{q}:4}$, and ${\mathrm{QSB}}_{\mathfrak{q}:4}$. These characteristics make the proposed method an efficient and reliable approach for solving nonlinear equations arising in applied dynamical systems and scientific computing.

5.6. Application: Biomedical Engineering Oxygen Transport Model [34]

Nonlinear equations also play a fundamental role in biomedical engineering, particularly in physiological modeling, pharmacokinetics, and biomedical signal processing. Many biological processes exhibit nonlinear behavior due to saturation effects, cooperative binding, and feedback mechanisms. As a result, determining equilibrium states in biomedical models frequently requires solving nonlinear algebraic equations.

A classical example arises in modeling oxygen transport in blood through the Hill equation, which describes the nonlinear relationship between oxygen partial pressure and hemoglobin saturation. The steady–state condition can be formulated as

$$f\left(x\right)={\displaystyle \frac{x^n}{P_{50}^n+x^n}}-S=0,$$

(59)

where x represents the oxygen partial pressure, $P_{50}$ denotes the partial pressure at which hemoglobin is 50% saturated, n is the Hill coefficient characterizing cooperative binding, and S is the measured oxygen saturation level.

For the numerical simulations, typical physiological parameters are selected as

$$n=2.7,P_{50}=26,S=0.75.$$

Solving the nonlinear equation $f\left(x\right)=0$ determines the oxygen partial pressure corresponding to the observed saturation level in blood. The proposed iterative scheme ${\mathrm{MSB}}_{\mathfrak{q}:\alpha }$ is applied to compute this root and its performance is compared with the existing methods ${\mathrm{SAN}}_{\mathfrak{q}:3}$, ${\mathrm{SAN}}_{\mathfrak{q}:4}$, and ${\mathrm{QSB}}_{\mathfrak{q}:4}$.

For the chosen parameter values, the nonlinear equation has the solution

$$x^{*}\approx 39.05574476,$$

which represents the oxygen partial pressure corresponding to the prescribed saturation level, computed to eight decimal places. The numerical comparison illustrates the robustness and rapid convergence behavior of the proposed method when applied to nonlinear biomedical models.

5.6.1. Efficiency Analysis

The efficiency of the considered numerical schemes is evaluated using the absolute error obtained at successive iterations. The comparison results are reported in Table 7, where ${\epsilon }_k=|x_k-\xi |$ denotes the absolute error at the kth iteration.

From the data presented in Table 7, it can be observed that the proposed scheme ${\mathrm{MSB}}_{\mathfrak{q}:\alpha }$ achieves significantly smaller absolute errors within a limited number of iterations when compared with the competing methods ${\mathrm{SAN}}_{\mathfrak{q}:3}$, ${\mathrm{SAN}}_{\mathfrak{q}:4}$, and ${\mathrm{QSB}}_{\mathfrak{q}:4}$. In particular, the proposed method rapidly reduces the approximation error to extremely small magnitudes, demonstrating improved convergence efficiency and strong numerical robustness for the considered nonlinear chemical reactor model. Values reported as $0.0$ correspond to errors below the working precision.

5.6.2. Stability and Computational Cost

The stability and computational efficiency of the iterative schemes are further analyzed using several performance indicators, including the maximum absolute error, CPU time, memory usage, and the total number of arithmetic operations required during the iterative process.

The results summarized in Table 8 show that the proposed ${\mathrm{MSB}}_{\mathfrak{q}:\alpha }$ method achieves an extremely small maximum error while maintaining competitive computational cost. In particular, the method requires a moderate number of arithmetic operations and reasonable memory usage while preserving strong numerical stability. Figure 6 illustrates the decay of the absolute error for different combinations of the parameters $(\alpha ,q)$ and several initial guesses when solving the nonlinear equation (59).

These results indicate that the ${\mathrm{MSB}}_{\mathfrak{q}:\alpha }$ scheme provides an efficient balance between accuracy and computational cost when applied to nonlinear problems arising in chemical reaction engineering.

5.6.3. Consistency for Fractional and Quantum Parameters

To investigate the robustness of the proposed algorithm, the numerical performance of the scheme is evaluated for different combinations of the fractional parameter $\alpha$ and the quantum parameter $\mathfrak{q}$.

For the consistency analysis, the initial vector is chosen as

$${\mathbf{x}}_0=\mathrm{vpa}\left([29.3,29.7,28.8,25.9,34.1]\right).$$

The obtained results are reported in Table 9 and illustrated in Figure 6. The numerical experiments indicate that the ${\mathrm{MSB}}_{\mathfrak{q}:\alpha }$ scheme maintains strong convergence characteristics across a wide range of parameter combinations.

In particular, the absolute error curves displayed in Figure 6 show that the proposed method consistently achieves rapid error decay, indicating stable convergence behavior for different values of the fractional and quantum parameters.

5.6.4. Fractal Basin Analysis and Dynamical Performance Comparison

To further examine the global dynamical behavior of the considered iterative schemes, fractal basin structures and convergence regions are analyzed. This analysis provides additional insight into the stability and reliability of the numerical methods when applied to nonlinear equations.

The comparative dynamical performance of the schemes is summarized in Table 10, while the corresponding fractal basin structures are illustrated in Figure 7. The comparison considers several indicators including fractal basin resolution, the percentage of convergence and divergence points, elapsed computational time, arithmetic operation counts, CPU time, and memory consumption.

The fractal basin plots indicate that the proposed Caputo fractional quantum scheme ${\mathrm{MSB}}_{\mathfrak{q}:\alpha }$ exhibits larger convergence regions and smaller divergence areas compared with the competing schemes. Furthermore, the proposed method requires fewer arithmetic operations and reduced computational cost while maintaining strong numerical stability.

Consequently, the ${\mathrm{MSB}}_{\mathfrak{q}:\alpha }$ method provides a robust and computationally efficient numerical tool for solving nonlinear equations arising in chemical reaction engineering and related applied dynamical systems, outperforming ${\mathrm{SAN}}_{\mathfrak{q}:3}$, ${\mathrm{SAN}}_{\mathfrak{q}:4}$, and ${\mathrm{QSB}}_{\mathfrak{q}:4}$.

5.7. Application in Diode Circuit Operating Point Analysis [35]

Nonlinear equations frequently arise in electrical engineering applications, particularly in the analysis of semiconductor circuits, nonlinear oscillators, and diode–based electronic systems. Determining the operating point of such circuits typically requires solving nonlinear transcendental equations that cannot be handled analytically. Efficient numerical root–finding techniques therefore play a crucial role in circuit simulation and electrical system design.

To demonstrate the applicability of the proposed Caputo fractional quantum iterative scheme ${\mathrm{MSB}}_{\mathfrak{q}:\alpha }$, we consider the nonlinear diode circuit equation

$$f\left(x\right)=I_s\left(e^{x/V_T}-1\right)-{\displaystyle \frac{V-x}{R}},$$

(60)

where $I_s$ denotes the saturation current of the diode, $V_T$ represents the thermal voltage, V is the supply voltage, and R is the circuit resistance. Solving the nonlinear equation $f\left(x\right)=0$ yields the steady–state operating voltage of the nonlinear circuit.

For the numerical experiments, representative circuit parameters are selected as

$$I_s={10}^{-12},V_T=0.0259,R=1000\Omega ,V=5\mathrm{V}.$$

The proposed ${\mathrm{MSB}}_{\mathfrak{q}:\alpha }$ scheme is compared with the existing iterative methods ${\mathrm{SAN}}_{\mathfrak{q}:3}$, ${\mathrm{SAN}}_{\mathfrak{q}:4}$, and ${\mathrm{QSB}}_{\mathfrak{q}:4}$. The comparison is carried out in terms of convergence efficiency, numerical stability, and dynamical behavior.

5.7.1. Efficiency Analysis

The efficiency of the considered iterative schemes is evaluated using the absolute error at successive iterations. The numerical results are summarized in Table 11, where ${\epsilon }_k=|x_k-\xi |$ denotes the absolute error at the kth iteration.

From the results presented in Table 11, it can be observed that the proposed scheme ${\mathrm{MSB}}_{\mathfrak{q}:\alpha }$ achieves substantially smaller absolute errors within a limited number of iterations when compared with the competing methods ${\mathrm{SAN}}_{\mathfrak{q}:3}$, ${\mathrm{SAN}}_{\mathfrak{q}:4}$, and ${\mathrm{QSB}}_{\mathfrak{q}:4}$.

Values reported as $0.0$ correspond to errors below the working precision.

In particular, the ${\mathrm{MSB}}_{\mathfrak{q}:\alpha }$ scheme rapidly reduces the approximation error to extremely small magnitudes under high-precision arithmetic, demonstrating improved convergence efficiency and strong numerical robustness for the nonlinear diode circuit equation.

5.7.2. Stability and Computational Cost

The stability and computational efficiency of the considered iterative schemes are evaluated using several performance indicators, including maximum absolute error, CPU time, memory usage, and the total number of arithmetic operations required during the iterative process. The corresponding results are summarized in Table 12, where the column “Max-Absolute Error” represents the maximum value of $|x_k-\xi |$ obtained during the iteration process.

From the data reported in Table 12, it can be observed that the proposed ${\mathrm{MSB}}_{\mathfrak{q}:\alpha }$ scheme maintains extremely small absolute errors while requiring competitive computational resources. In particular, the proposed method exhibits moderate CPU time and memory consumption while preserving strong numerical stability.

These results indicate that the ${\mathrm{MSB}}_{\mathfrak{q}:\alpha }$ method provides an efficient balance between numerical accuracy and computational cost when applied to nonlinear equations arising in diode circuit analysis.

5.7.3. Consistency for Fractional and Quantum Parameters

The robustness of the ${\mathrm{MSB}}_{\mathfrak{q}:\alpha }$ scheme is further investigated for different combinations of the fractional parameter $\alpha$ and the quantum parameter $\mathfrak{q}$. For the consistency analysis, the initial vector is chosen as

$${\mathbf{x}}_0=\mathrm{vpa}\left([0.7,1.0,1.5,2.0,0.6]\right).$$

The numerical outcomes are reported in Table 13 and illustrated in Figure 8. The numerical outcomes are reported in Table 13 and illustrated in Figure 8. The error curves displayed in Figure 8 indicate that the ${\mathrm{MSB}}_{\mathfrak{q}:\alpha }$ scheme exhibits rapid error decay and stable convergence behavior for different values of the fractional and quantum parameters.

The error curves displayed in Figure 8 show that the ${\mathrm{MSB}}_{\mathfrak{q}:\alpha }$ scheme exhibits rapid error decay and stable convergence behavior for different values of the fractional and quantum parameters. These observations confirm that the method maintains stable and consistent convergence across a wide range of parameter settings.

5.7.4. Dynamical Behavior and Computational Performance Comparison

The global dynamical behavior of the considered iterative schemes is analyzed through fractal basin representations, convergence–divergence distributions, and computational performance indicators. These graphical and quantitative analyses provide valuable insights into the stability, robustness, and efficiency of the iterative processes when applied to nonlinear equations.

Figure 9 illustrates the fractal basins of attraction associated with the methods ${\mathrm{SAN}}_{\mathfrak{q}:3}$, ${\mathrm{SAN}}_{\mathfrak{q}:4}$, ${\mathrm{QSB}}_{\mathfrak{q}:4}$, and the proposed ${\mathrm{MSB}}_{\mathfrak{q}:\alpha }$ scheme. In these plots, each point in the complex plane represents an initial guess, while the corresponding color indicates the root toward which the iterative process converges. Regions with smooth color transitions correspond to stable convergence toward the root, whereas irregular or fragmented boundaries indicate sensitive dynamical behavior, where small perturbations in the initial guess may lead to different convergence outcomes or divergence.

A visual comparison of the fractal patterns reveals that the ${\mathrm{MSB}}_{\mathfrak{q}:\alpha }$ scheme exhibits significantly larger and smoother convergence regions than the competing methods. The corresponding basin boundaries are also less irregular, indicating reduced sensitivity to the choice of initial guesses and improved numerical stability of the iterative process.

Table 14 summarizes several dynamical and computational indicators used to compare the considered iterative schemes. In particular, the percentage of convergence points (C-Points) and divergence points (D-Points) measures the stability of each method over a large set of initial guesses in the complex plane. Additional performance metrics, including execution time, arithmetic operation counts, CPU time, and memory usage, provide further information on the computational efficiency of the corresponding algorithms.

The results show that the ${\mathrm{MSB}}_{\mathfrak{q}:\alpha }$ scheme achieves the highest convergence percentage (97%) and the lowest divergence rate (3%) among the tested methods. At the same time, it requires fewer arithmetic operations and exhibits lower execution time and memory usage. These characteristics confirm the improved efficiency of the proposed scheme compared with ${\mathrm{SAN}}_{\mathfrak{q}:3}$, ${\mathrm{SAN}}_{\mathfrak{q}:4}$, and ${\mathrm{QSB}}_{\mathfrak{q}:4}$.

Overall, the obtained results demonstrate that the Caputo fractional quantum iterative scheme ${\mathrm{MSB}}_{\mathfrak{q}:\alpha }$ provides enhanced dynamical stability, larger convergence domains, and improved computational efficiency. Consequently, the proposed method represents a reliable and efficient numerical tool for solving nonlinear equations arising in electrical circuit analysis and other nonlinear dynamical systems.

6. Conclusion

This study introduced a novel Caputo fractional quantum iterative scheme, denoted by ${\mathrm{MSB}}_{\mathfrak{q}:\alpha }$, for solving nonlinear equations. The proposed method combines the Caputo fractional operator with the quantum derivative framework to construct an efficient root–finding algorithm. This hybrid formulation introduces additional flexibility through the fractional parameter $\alpha$ and the quantum parameter $\mathfrak{q}$, enabling improved convergence behavior and enhanced numerical stability compared with existing iterative schemes such as ${\mathrm{SAN}}_{\mathfrak{q}:3}$, ${\mathrm{SAN}}_{\mathfrak{q}:4}$, and ${\mathrm{QSB}}_{\mathfrak{q}:4}$.

Extensive numerical experiments (Table 3, Table 4, Table 5, Table 6, Table 7, Table 8, Table 9, Table 10, Table 11, Table 12, Table 13 and Table 14 and Figure 3, Figure 4, Figure 5, Figure 6, Figure 7, Figure 8 and Figure 9) were conducted to validate the theoretical properties of the proposed method. The error analyses demonstrate that the ${\mathrm{MSB}}_{\mathfrak{q}:\alpha }$ scheme exhibits faster error reduction and achieves high numerical accuracy within fewer iterations than the competing methods. Furthermore, the computational performance results indicate that the proposed scheme requires reduced CPU time, fewer arithmetic operations, and lower memory usage while maintaining high numerical precision (see e.g., Table 4, Table 8 and Table 12, respectively).

The robustness of the method was further investigated through parametric and dynamic studies (Table 1 and Figure 1 and Figure 2) involving different values of the fractional parameter $\alpha$ and the quantum parameter $\mathfrak{q}$. The results show that the ${\mathrm{MSB}}_{\mathfrak{q}:\alpha }$ scheme maintains stable and consistent convergence across a wide range of parameter settings. Moreover, the dynamical and fractal analyses from Table 6, Table 10, Table 14 and Figure 5, Figure 7 and Figure 9 reveal that the proposed method possesses larger convergence basins, reduced divergence regions, and improved global dynamical stability compared with the existing iterative schemes.

To illustrate its practical applicability, the proposed method was applied to nonlinear models arising in chemical reaction engineering, biomedical systems, and electrical circuit analysis. In all cases, the ${\mathrm{MSB}}_{\mathfrak{q}:\alpha }$ scheme demonstrated superior convergence speed, improved stability, and higher computational efficiency compared with the considered numerical methods as presented in Table 3, Table 4, Table 7, Table 8, Table 11, Table 12, Figure 4, Figure 6 and Figure 8 respectively.

Despite these promising results, several aspects deserve further investigation. The performance of the proposed method may depend on the appropriate selection of the fractional and quantum parameters, and the current formulation has been primarily developed for scalar nonlinear equations. Future work may therefore focus on extending the ${\mathrm{MSB}}_{\mathfrak{q}:\alpha }$ framework to systems of nonlinear equations and large–scale multidimensional problems. In addition, adaptive strategies for selecting optimal fractional and quantum parameters could further improve the efficiency of the algorithm. Further studies may also explore deeper dynamical properties of the scheme, including high–resolution fractal basin structures and global convergence characteristics.

Overall, the proposed Caputo fractional quantum iterative scheme provides an efficient, stable, and computationally attractive approach for solving nonlinear equations. The integration of fractional and quantum operators opens promising directions for the development of advanced numerical algorithms in applied mathematics, computational science, and engineering applications.

All authors declare that this work complies with ethical guidelines set by the Committee on Publication Ethics (COPE). The authors declare that they have no competing financial interests and personal relationships that could have appeared to influence the research in this paper.