A Robust and Versatile Numerical Framework for Modeling Complex Fractional Phenomena: Applications to Riccati and Lorenz Systems

Abstract

The fractional differential quadrature method (FDQM) with generalized Caputo derivatives is used in this paper to show a new numerical way to solve fractional Riccati equations and fractional Lorenz systems. Unlike previous FDQM applications that have primarily focused on linear problems, our work pioneers the use of this method for nonlinear fractional initial value problems. By combining Lagrange interpolation polynomials and discrete singular convolution (DSC) shape functions with the generalized Caputo operator, we effectively transform nonlinear fractional equations into algebraic systems. An iterative method is then utilized to address the nonlinearity. Our numerical results, obtained using MATLAB, demonstrate the exceptional accuracy and efficiency of this approach, with convergence rates reaching 10⁻⁸. Comparative analysis with existing methods highlights the superior performance of the DSC shape function in terms of accuracy, convergence speed, and reliability. Our results highlight the versatility of our approach in tackling a wider variety of intricate nonlinear fractional differential equations.

1. Introduction

Many phenomena in chemistry, biology, acoustics, psychology, control theory, rheology, damping laws, diffusion processes, and other fields of science have been successfully modeled using fractional-order derivatives in recent years. This is because fractional calculus can be used to successfully model a physical phenomenon that is dependent not only on the time instant but also on the prior time history [1,2,3,4,5,6]. Hence, numerous physical problems are defined by fractional differential equations (FDEs), and solving these equations has been the focus of several studies in recent years. Several techniques have recently been developed to solve FDEs, including numerical and analytical techniques. Various methods, including homotopy perturbation [7,8,9], homotopy analysis [10,11], Taylor matrix [12], Adomian decomposition [11], and Haar wavelet [13] methods, have been employed to solve the fractional-order Riccati differential equation. Unfortunately, the convergence region for the corresponding outcomes is relatively limited. The fractional derivative operator’s unique properties can make the numerical solution of fractional equations challenging, particularly in high-dimensional spaces. To address this challenge, numerical methods such as the Finite Difference Method (FDM) [14,15,16,17], Galerkin [18,19,20,21], Collocation [22,23,24,25], and finite volume element methods [26,27,28] have been utilized to tackle such fractional equations.

Liu et al. [29] introduced a radial basis function finite difference approach for studying the time fractional convection equation. Saadeh [30] employed the finite-difference and space finite-volume techniques to solve fractional diffusion equations. To address a two-dimensional space fractional diffusion equation, Tuan et al. [31] employed finite difference discretization with Caputo derivatives. Devshali and Arora [32] proposed differential transform and differential quadrature methods for solving the fractional diffusion equation. Odibat and Momani [9] utilized a modified homotopy perturbation technique to solve fractional Riccati differential equations (FRDEs). Khader [33] applied the fractional Chebyshev FDM to solve FRDEs. Li et al. [34] solved FRDEs via the quasi-linearization method. Sakar et al. [34] explained an iterative reproducing kernel Hilbert space technique to obtain the solutions of FRDEs. Agheli [35] explained an iterative reproducing kernel Hilbert space technique to obtain the solutions of FRDEs. Agheli [36] presented numerical solutions for solving FRDEs using trigonometric basic functions. Liu et al. [37] offered the Laplace transform and quadrature rule with Caputo sense to solve FRDEs.

In the last ten years, chaos has emerged as a popular topic in fractional calculus [38]. The chaotic behavior becomes more complicated because the equation contains fractional orders. Numerical methods were developed to analyze nonlinear dynamics to better understand physical phenomena. To find numerical solutions to various nonlinear fractional differential equations, the predictor–corrector method (P-C) was developed [39,40,41]. Fuzzy fractional differential and fractional delay equations [42,43,44,45] demonstrate the emergence of new trends in fractional differential equations.

Despite yielding fruitful results, finding more general chaotic differential equations remains an intriguing task. A generalized fractional derivative was recently proposed in [46]. Fractional derivatives have demonstrated superior performance compared to regular derivatives in several respects and may have even more real-world applications. One such application is in image encryption, where fractional differential equations have been suggested as a means of introducing chaos [47,48]. This means that image encryption results can be made more secure by using fractional chaotic equations with two parameters. This derivative was recently proposed in quantum mechanics [49]. Furthermore, two-parameter models in control theory and diffusion issues can have degrees of freedom in control and fitting. This derivative and its applications are depicted as a new direction in fractional calculus. Li and Chen [50] demonstrated the chaotic behaviors in the fractional order Chen system. Alomari [51] used the step homotopy analysis technique to solve the fractional chaotic Chen system. Luo and Wang [44] solved chaos in the fractional-order complex Lorenz system and its synchronization. Petráš [2] introduced a new classification of the fractional-order Lorenz-type systems. Erturka and Kumar [52] presented a solution for a COVID-19 model using new generalized Caputo-type fractional derivatives. Xu et al. [53] studied numerical and analytical solutions of a new generalized fractional diffusion equation. Kumar et al. [54] proposed a new technique to solve generalized Caputo-type FDEs with the example of a computer virus model.

The primary goal of this paper is to apply the novel fractional Differential Quadrature Method (FDQM) with generalized Caputo definition fractional to solve nonlinear initial value fractional problems. Two different shape functions, the Lagrange interpolation [55,56] and the regularized Shannon kernel [57,58,59,60], have been successfully employed to address initial value problems involving fractional derivatives. To exhibit the efficacy, efficiency, and capabilities of the proposed algorithm, two test problems were investigated, including FRDEs and the fractional Lorenz system. Then, by the proposed methods, the given problems are reduced to a system of nonlinear algebraic equations, and by solving this system via the iterative method, we obtain the solution of FRDEs and the fractional Lorenz system. Furthermore, we create a MATLAB code for each approach to obtain a numerical solution for the two problems under consideration. A comparison between the computed results and previous analytical and numerical [61,62,63,64,65] methodologies is included to demonstrate the validity and applicability of the proposed methods. Furthermore, we conducted some parametric investigations to showcase the reliability of our techniques in the presence of fractional order derivatives.

This paper introduces a novel fractional Differential Quadrature Method (FDQM) to solve nonlinear initial value fractional problems. This method employs the generalized Caputo fractional derivative and utilizes Lagrange interpolation and the Regularized Shannon kernel as shape functions. Numerical simulations demonstrate the method’s superior accuracy, efficiency, and versatility in handling fractional Riccati and Lorenz systems. The FDQM’s potential applications extend to various fields where it can be used to model complex systems with memory effects and nonlinearities [66].

2. Formulation of the Problem

The following two nonlinear fractional differential equations serve as examples to illustrate the capabilities of our proposed methods:

2.1. Our 1st Example Is the Fractional Riccati Equation

$${\displaystyle \frac{{\mathrm{d}}^{\alpha ,\rho }\upsilon \left(\mathrm{t}\right)}{\mathrm{d}\mathrm{t}}}=2\upsilon \left(\mathrm{t}\right)-{\upsilon }^2\left(\mathrm{t}\right)+1, \mathrm{when} \left(0< \mathrm{t} \mathrm{and} 0<\alpha \le 1\right)$$

(1)

where ${\displaystyle \frac{{\mathrm{d}}^{\alpha ,\rho }}{\mathrm{d}\mathrm{t}}}$ is the operator of the generalized Caputo-type fractional derivative [61].

Also, the initial condition for FRDE is:

$$\upsilon \left(0\right)=0$$

(2)

In addition, the exact solution of FRDE at $\alpha =1, \mathrm{and} \rho =1$ is given by [67]:

$$\upsilon \left(\mathrm{x},\mathrm{t}\right)=1+\sqrt{2} \tanh \left[\sqrt{2}\mathrm{t}+\log \sqrt{{\displaystyle \frac{\sqrt{2}-1}{\sqrt{2}+1}}}\right]$$

(3)

2.2. Our 2nd Example Is the Fractional Lorenz System

$${\displaystyle \frac{{\mathrm{d}}^{\mathsf{\alpha },\mathsf{\rho }}\mathrm{X}\left(\mathrm{t}\right)}{\mathrm{d}\mathrm{t}}}=\mathsf{\lambda }\left(\mathrm{Y}\left(\mathrm{t}\right)-\mathrm{X}\left(\mathrm{t}\right)\right)$$

(4)

$${\displaystyle \frac{{\mathrm{d}}^{\mathsf{\alpha },\mathsf{\rho }}\mathrm{Y}\left(\mathrm{t}\right)}{\mathrm{d}\mathrm{t}}}=\left(\mathsf{\varphi }-\mathsf{\lambda }\right)\mathrm{X}\left(\mathrm{t}\right)-\mathrm{X}\left(\mathrm{t}\right)\mathrm{Z}\left(\mathrm{t}\right)+\mathsf{\varphi }\mathrm{Y}\left(\mathrm{t}\right)$$

(5)

$${\displaystyle \frac{{\mathrm{d}}^{\mathsf{\alpha },\mathsf{\rho }}\mathrm{Z}\left(\mathrm{t}\right)}{\mathrm{d}\mathrm{t}}}=\mathrm{X}\left(\mathrm{t}\right)\mathrm{Y}\left(\mathrm{t}\right)-\mathsf{\beta } \mathrm{Z}\left(\mathrm{t}\right)$$

(6)

where $\mathsf{\lambda }, \mathsf{\varphi }, \mathrm{and} \mathsf{\beta }$ are constant parameters $\in \mathbb{R}$ that affect chaotic behavior.

Consequently, the fractional Lorenz system is subject to the following initial condition:

$$\mathrm{X}\left(0\right)={\mathrm{x}}_0, \mathrm{Y}\left(0\right)={\mathrm{y}}_0, \mathrm{Z}\left(0\right)={\mathrm{z}}_0$$

(7)

3. Method of Solution

This paper introduces a novel application of the FDQM to solve initial value fractional problems. We employ two distinct shape functions, Lagrange interpolation and the regularized Shannon kernel, in conjunction with the generalized Caputo fractional derivative to transform fractional problems into nonlinear algebraic systems.

We begin by defining a fractional derivative, of which several definitions exist. In this work, we utilize the recently proposed generalized Caputo definition.

3.1. Generalized Caputo-Kind Fractional Derivative [63]

The fractional derivative has good memory effects compared to ordinary calculus. FDEs are realized in model problems in fluid flow, viscoelasticity, finance, engineering, and other areas of applications.

Caputo’s Fractional Derivative

A concise overview of Caputo’s fractional derivative is presented in this section. This definition, which is derived from the Riemann–Liouville Fractional Derivative [68], is explained in greater detail in our prior publication [69].

Suppose $\alpha \in R^{+}$, If $\mathbb{N}$ is a positive integer, and $\mathbb{N}-1<\mathsf{\alpha }\le \mathbb{N}$. According to Riemann–Liouville fractional, which is one of the most researched definitions, the fraction derivative of a function $\upsilon \left(t\right)$ of order $\mathsf{\alpha }$ is defined as follows:

$$D_c^{\alpha }\upsilon (t)={\displaystyle \frac{1}{\Gamma (\mathbb{N}-\alpha )}} {\displaystyle \frac{d^{\mathbb{N}}}{d{t}^{\mathbb{N}}}} {\displaystyle \underset{c}{\overset{t}{\int }}{(t-x)}^{\mathbb{N}-\alpha -1}\hspace{0.33em}{\upsilon }^{\mathbb{N}}(x)\hspace{0.33em}dx}$$

(8)

Generalized Caputo’s Fractional Derivative of operator $D_{c_{+}}^{\alpha ,\rho }$ and order $\alpha$ is defined as [52,53,54]:

$$\begin{array}{l}{D}_{c_{+}}^{\alpha ,\rho }\upsilon (t)={\displaystyle \frac{{\rho }^{\alpha -\mathbb{N}+1}}{\Gamma (\mathbb{N}-\alpha )}}{\displaystyle \underset{c}{\overset{t}{\int }}{\mathrm{x}}^{\rho -1}{(t^{\rho }-x^{\rho })}^{\mathbb{N}-\alpha -1}\hspace{0.33em}{\left(x^{1-\rho }\frac{d}{dx}\right)}^{\mathbb{N}}\upsilon (x)\hspace{0.33em}d\mathrm{x},\hspace{0.33em}\hspace{0.33em}}\\ \hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\mathbb{N}-1<\alpha <\mathbb{N},\hspace{0.33em}\hspace{0.33em}\rho >0,\hspace{0.33em}\hspace{0.33em}c\ge 0\end{array}$$

(9)

where $\mathrm{c}$ is the lower limit of integration.

Consequently, the solution to Equation (9) can be written as [52,53,54]:

$$D_{c_{+}}^{\alpha ,\rho }{(t^{\rho }-x^{\rho })}^{\gamma }={\rho }^{\alpha }{\displaystyle \frac{\Gamma \left(\gamma +1\right)}{\Gamma \left(\gamma -\alpha +1\right)}}{\left(t^{\rho }-c^{\rho }\right)}^{\gamma -\alpha }$$

(10)

Next, we will discuss the differential quadrature technique using the specified shape functions.

3.2. Our First Shape Function Is Lagrange Interpolation Polynomial Based DQM (PDQM)

Within this partition, the functional values of any unknown $\mathrm{P}$ at a specific set of $\mathrm{N}$ grid points can be represented using this shape function, as described in [55,56]:

$$\mathrm{P}({\mathrm{t}}_{\mathrm{i}})={\displaystyle \sum _{\mathrm{j}=1}^{\mathrm{N}}\frac{1}{{\mathrm{t}}_{\mathrm{i}}-{\mathrm{t}}_{\mathrm{j}}}\times \frac{{\displaystyle \prod _{\mathrm{k}=1}^{\mathrm{N}}({\mathrm{t}}_{\mathrm{i}}-{\mathrm{t}}_{\mathrm{k}})}}{{\displaystyle \prod _{\mathrm{j}=1,\mathrm{j}\ne \mathrm{k}}^{\mathrm{N}}({\mathrm{t}}_{\mathrm{j}}-{\mathrm{t}}_{\mathrm{k}})}}\mathrm{P}(t_j)}\hspace{0.17em},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}(\mathrm{i}=1:\mathrm{N})$$

(11)

As a result, the following are the different derivatives of this unknown $\mathrm{P}$:

$${\displaystyle \frac{{\partial }^{\mathrm{r}}\mathrm{P}}{\partial {\mathrm{t}}^{\mathrm{r}}}}\left|\begin{array}{c} \\ \mathrm{t}={\mathrm{t}}_{\mathrm{i}}\end{array}\right.={\displaystyle \sum _{\mathrm{j}=\mathrm{1}}^{\mathrm{N}}{\mathrm{a}}_{\mathrm{i}\mathrm{j}}^{(\mathrm{r})}}\mathrm{P}({\mathrm{t}}_{\mathrm{j}})\hspace{0.17em},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}(\mathrm{i}=1:\mathrm{N})$$

(12)

where ${\mathrm{a}}_{\mathrm{i}\mathrm{j}}^{(\mathrm{r})}$ is the rth derivative weighting coefficient. However, determining the weighting coefficients is critical to DQM accuracy. As a result, they differ based on the shape function.

Differentiating Equation (12) results in the calculation of the weighting coefficients ${\mathrm{a}}_{\mathrm{i}\mathrm{j}}^{(1)}$ and ${\mathrm{a}}_{\mathrm{i}\mathrm{j}}^{(2)}$, representing the first and second derivatives.

$${\mathrm{a}}_{\mathrm{i}\mathrm{j}}^{(1)}=\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\left\{\begin{array}{c}{\displaystyle \frac{1}{({\mathrm{t}}_{\mathrm{i}}-{\mathrm{t}}_{\mathrm{j}})}}{\displaystyle \prod _{\begin{array}{l}\mathrm{k}=1,\\ \mathrm{k}\ne \mathrm{i},\mathrm{j}\end{array}}^{\mathrm{N}}\frac{({\mathrm{t}}_{\mathrm{i}}-{\mathrm{t}}_{\mathrm{k}})}{({\mathrm{t}}_{\mathrm{j}}-{\mathrm{t}}_{\mathrm{k}})}}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{i}\ne \mathrm{j}\\ \\ -{\displaystyle \sum _{\begin{array}{l}\mathrm{j}=1,\\ \mathrm{j}\ne \mathrm{i}\end{array}}^{\mathrm{N}}{\mathrm{a}}_{\mathrm{i}\mathrm{j}}^{(1)}}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\mathrm{i}=\mathrm{j}\end{array}\right.\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em},{ \mathrm{a}}_{\mathrm{i}\mathrm{j}}^{(2)}=\hspace{0.17em}\left[{\mathrm{a}}_{\mathrm{i}\mathrm{j}}^{(1)}\right]\left[{\mathrm{a}}_{\mathrm{i}\mathrm{j}}^{(1)}\right]$$

(13)

The distribution of grid points N, whether uniform or non-uniform, significantly influences the accuracy of the PDQM results. The non-uniform distribution is defined by the following equation, based on Chebyshev’s distribution:

$${\mathrm{x}}_{\mathrm{i}}={\displaystyle \frac{1}{2}}{\mathrm{L}}_{\mathrm{x}}\left(1-\cos \hspace{0.33em}\left({\displaystyle \frac{\mathsf{\pi }\left(\mathrm{i}-1\right)}{\mathrm{N}-\mathrm{1}}}\right)\right),\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}(\mathrm{i}=1:\mathrm{N})$$

(14)

3.3. Discrete Singular Convolution-Based DQM (DSCDQM)

In this part, according to previous research [59,70,71], singular convolution can be expressed as follows:

$$\mathrm{Y}(\mathrm{t})=(\mathrm{F}\ast \mathrm{H})(\mathrm{t})={\displaystyle \underset{-\infty }{\overset{\infty }{\int }}\mathrm{F}(\mathrm{t}-\mathrm{s}) \mathrm{H}(\mathrm{s}) \mathrm{d}\mathrm{s}}$$

(15)

$\mathrm{F}(\mathrm{t}-\mathrm{s})$ and $\mathrm{H}(\mathrm{t})$ denote a singular kernel and a test function space element, respectively.

The choice of kernel type determines the shape function used in this technique. Given the variety of available kernels, we have selected the kernel demonstrated to have the highest accuracy [69] to represent the functional values of the unknown $\mathrm{P}$ and its derivatives at a specified number of grid points $\mathrm{N}$:

•

Our second shape function is the Regularized Shannon kernel (DSCDQM–RSK)

$$\begin{array}{l}\mathrm{P}({\mathrm{t}}_{\mathrm{i}})={\displaystyle \sum _{\mathrm{j}=-\mathrm{M}}^{\mathrm{M}}\langle \frac{\sin \left[\frac{\mathsf{\pi }({\mathrm{t}}_{\mathrm{i}}-{\mathrm{t}}_{\mathrm{j}})}{\Delta }\right]}{\frac{\mathsf{\pi }({\mathrm{t}}_{\mathrm{i}}-{\mathrm{t}}_{\mathrm{j}})}{\Delta }}\hspace{0.17em}\hspace{0.17em}\exp (\frac{-{({\mathrm{t}}_{\mathrm{i}}-{\mathrm{t}}_{\mathrm{j}})}^2}{2{\mathsf{\sigma }}^2})\rangle \hspace{0.17em}\hspace{0.17em}\mathrm{P}({\mathrm{x}}_{\mathrm{j}})},\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}(\mathrm{i}=-\mathrm{N}:\mathrm{N}),\mathsf{\sigma }=(\mathrm{h}\times \Delta \hspace{0.17em})\hspace{0.17em}>0\end{array}$$

(16)

The parameters $\mathsf{\sigma }, \mathrm{h}, \mathrm{and} \Delta$ represent the Regularized Shannon factor, the computational parameter, and the mesh size, respectively.

$$\begin{array}{l}{\displaystyle \frac{{\partial }^{\mathrm{r}}\mathrm{P}}{\partial {\mathrm{t}}^{\mathrm{r}}}}\left|\mathrm{t}={\mathrm{t}}_{\mathrm{i}}\right.={\displaystyle \sum _{\mathrm{j}=1}^{\mathrm{N}}{\mathrm{a}}_{\mathrm{i}\mathrm{j}}^{(\mathrm{r})}}\mathrm{P}({\mathrm{t}}_{\mathrm{j}})\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.33em}\hspace{0.33em}\\ \left(\mathrm{i}=-\mathrm{N}:\mathrm{N}\right)\hfill \end{array}$$

(17)

Differentiating Equation (15) allows us to determine the coefficients ${\mathrm{a}}_{\mathrm{i}\mathrm{j}}^{(1)}\hspace{0.33em}\hspace{0.33em}\mathrm{and}{ \mathrm{a}}_{\mathrm{i}\mathrm{j}}^{(2)}$, as described in [72]:

$${\mathrm{a}}_{\mathrm{i}\mathrm{j}}^{(1)}=\left\{\begin{array}{c}{\displaystyle \frac{{(-1)}^{\mathrm{i}-\mathrm{j}}}{\Delta (\mathrm{i}-\mathrm{j})}}\exp (-{\Delta }^2({\displaystyle \frac{{(\mathrm{i}-\mathrm{j})}^2}{2{\mathsf{\sigma }}^2}})),\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{i}\ne \mathrm{j}\\ \\ 0\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.17em}\hspace{0.17em}\mathrm{i}=\mathrm{j}\end{array}\right.\hspace{0.17em}\hspace{0.17em},{\mathrm{a}}_{\mathrm{i}\mathrm{j}}^{(2)}=\left\{\begin{array}{c}\left[2{\displaystyle \frac{{(-1)}^{\mathrm{i}-\mathrm{j}+1}}{{\mathrm{a}}^2{(\mathrm{i}-\mathrm{j})}^2}}+{\displaystyle \frac{1}{{\mathsf{\sigma }}^2}}\right]\exp \left(-{\Delta }^2({\displaystyle \frac{{(\mathrm{i}-\mathrm{j})}^2}{2{\mathsf{\sigma }}^2}}\right),\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{i}\ne \mathrm{j}\\ \\ -{\displaystyle \frac{1}{{\mathsf{\sigma }}^2}}-{\displaystyle \frac{{\mathsf{\pi }}^2}{3{\Delta }^2}}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{i}=\mathrm{j}\end{array}\right.$$

(18)

The kernel type, grid points (N), and bandwidth (2K + 1) are all important parameters that affect the convergence and accuracy of the solutions, as our analysis shows.

Now, after mentioning DQM based on two shape functions, we will demonstrate the effect of the generalized Caputo’s fractional derivative, which is shown in Equation (9) on the FDQM in Equation (12), to determine the weighting coefficients $a_{ij}^{\alpha }$ for $\alpha \in \left(0,1\right]$ and $\mathsf{\rho }>0$, as follows [69]:

$$D_{c_{+}}^{\alpha ,\rho }\mathrm{P}(t)=\left\{\begin{array}{l}{\displaystyle \sum _{j=1}^N{a}_{ij}^{\alpha ,\rho } \mathrm{P}(t_j},x),\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}0<\alpha \le 1\hspace{0.33em}\hspace{0.33em},\hspace{0.33em}\hspace{0.17em}\rho >0\hspace{0.33em}\hspace{0.33em}\\ {\displaystyle \sum _{j=1}^N{a}_{ij}^{(1)} \mathrm{P}(t_j},x)\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\alpha =\rho =1\end{array}\right.\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}$$

(19)

Hence, the weighting coefficient ${\mathrm{a}}_{\mathrm{i}\mathrm{j}}^{\mathsf{\alpha },\mathsf{\rho }}$ is calculated as:

$${\mathrm{a}}_{\mathrm{i}\mathrm{j}}^{\mathsf{\alpha },\mathsf{\rho }}={\mathrm{A}}^{1-\mathsf{\alpha }}{\mathsf{\rho }}^{\mathsf{\alpha }}{ \mathrm{a}}_{\mathrm{i}\mathrm{j}}^{(1)}+{\displaystyle \frac{{\mathsf{\rho }}^{\mathsf{\alpha }}{\mathrm{a}}_{1,\mathrm{j}}^{(1)}}{\Gamma (2-\mathsf{\alpha })}}\hspace{0.33em}{({\mathrm{t}}^{\mathsf{\rho }}-{\mathrm{c}}^{\mathsf{\rho }})}^{1-\mathsf{\alpha }},\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}{\mathrm{A}}_{\mathrm{i}\mathrm{j}}={\mathrm{a}}_{\mathrm{i}\mathrm{j}}^{(1)}-{\mathrm{a}}_{1\mathrm{j}}^{(1)}$$

(20)

3.4. Algorithm: Fractional Differential Quadrature Method (FDQM) for Nonlinear Initial Value Problems

This pseudo-code (Algorithm 1) outlines the key steps for implementing the proposed numerical framework for solving nonlinear fractional differential equations using the FDQM approach.

Algorithm 1: Fractional Differential Quadrature Method (FDQM) for Nonlinear Initial Value Problems

Input: - Fractional order, constants in the fractional differential equation - Nonlinear fractional differential equation. - Initial conditions. - Grid points (N) - Shape functions (Lagrange, Regularized Shannon) - Tolerance for convergence (ε) Output: - The solution of the fractional differential equation numerically. Define the generalized Caputo fractional derivative operator. Initialize grid points x1, x2, …, xN based on Chebyshev distribution or uniform distribution. For each time step “t: ▪ Construct the shape functions using Lagrange interpolation and Regularized Shannon kernel. ▪ Formulate the algebraic system from the fractional differential equation using FDQM: Apply shape functions to approximate the unknown function and its derivatives. Substitute into the original FDE to derive a system of nonlinear algebraic equations. Initialize solution guess for the unknown function. While not converged (|| fnew − fold|| > ε) Solve the algebraic system iteratively: Use a numerical method (iterative differential quadrature method) to update the solution. Update fold with fnew End while Return the numerical solution for the fractional differential equation at the specified time steps. End Algorithm

4. Numerical Results

Now that it is easier to understand FDQM with different shape functions such as PDQM [55,56], and DSCDQM–RSK [59,70,71] based on the generalized Caputo definition fractional derivative, two examples will be given here and then will be discussed. In all these examples, MATLAB software(R2022b) is used for computations and graphs. The primary goal of this article is to learn about the performance, validity, efficiency, and accuracy of developed techniques by comparing the computed results to previous numerical and analytical solutions.

4.1. Problem 4.1

We introduce the first example fractional Riccati equation after substituting the Equations (19) and (20) for the proposed methods in Equation (1) as follows:

$$\sum _{\mathrm{j}=1}^{\mathrm{L}}{\mathrm{a}}_{\mathrm{i}\mathrm{j}}^{\mathsf{\alpha },\mathsf{\rho }} \mathsf{\upsilon }\left({\mathrm{t}}_{\mathrm{j}}\right)=2\sum _{\mathrm{j}=1}^{\mathrm{N}}{\mathsf{\delta }}_{\mathrm{ij}} \mathsf{\upsilon }\left({\mathrm{t}}_{\mathrm{j}}\right)-{\left(\sum _{\mathrm{j}=1}^{\mathrm{N}}{\mathsf{\delta }}_{\mathrm{ij}} \mathsf{\upsilon }\left({\mathrm{t}}_{\mathrm{j}}\right)\right)}^2+1$$

(21)

The governing Equation (21) is also used to deal with the initial condition (2). To solve the nonlinear problem, the iterative method is applied [55,72,73]. As a first step, the governing equation is solved as a linear system. Then we solve them iteratively as a nonlinear system until we reach the requisite convergence, which is as follows:

$$\left|\frac{{\mathsf{\upsilon }}_{\mathrm{m}+1}}{{\mathsf{\upsilon }}_{\mathrm{m}}}\right|<1, \mathrm{where} \mathrm{m}=0,1,2,\dots$$

(22)

Also, to assess the convergence and accuracy of the developed methods, we use the error computation method:

$${\mathrm{L}}_{\infty } \mathrm{Error}= \underset{1\le \mathrm{i}\le \mathrm{N}}{\max } \left|{\upsilon }_{\mathrm{numerical}}-{\upsilon }_{\mathrm{exact}}\right|$$

(23)

Now, the obtained results will be demonstrated as follows:

The effect of applying PDQM with uniform and non-uniform grid distributions on the computation of the fractional Riccati equation with fractions $\left(\mathsf{\alpha }=1, \mathsf{\rho }=1\right)$ at different grid points (N) and times (T) is shown in

Table 1. Computation of $\mathsf{\upsilon }\left(t\right)$ via uniform and non-uniform PDQM for fractional Riccati equation with fraction $\left(\mathsf{\alpha }=1, \mathsf{\rho }=1\right)$ at different grid points (N), and times (T).

T	PDQ Solutions							Previous Solutions
	N	Uniform	CPU (s)	Error	Non-Uniform	CPU (s)	Error	N	Earlier Numerical [64]	Exact [61]
1	4	1.59030488	0.016	0.0971	1.64023865	0.013	0.04721	10	1.68745117	1.68949839
	6	1.68745374	0.018	0.0015	1.67122383	0.017	0.01774	20	1.68896723
	8	1.68821427	0.018	0.0011	1.68941502	0.017	5.16 × 10−5	40	1.68936339
	9	1.68921673	0.019	0.0002	1.68948616	0.018	2.18 × 10−5	80	1.68946438
	11	1.68941775	0.020	7.2 × 10−5	1.68949815	0.018	8.29 × 10−6	160	1.68948986
	12	1.68948043	0.021	1.6 × 10−5	1.68949820	0.020	1.95 × 10−6	320	1.68949625
	13	1.68949755	0.021	3.1 × 10−7	1.68949839	0.020	1.2 × 10−8	640	1.68949786
2	8	2.34267643	0.029	0.001168	2.35647559	0.019	0.012631	20	2.35530727	2.35777165
	12	2.35721628	0.031	0.00056	2.35777266	0.020	3.73 × 10−6	40	2.35721255
	16	2.35757661	0.032	0.00014	2.35777175	0.021	6.14 × 10−5	80	2.35763805
	18	2.35773329	0.033	3.27 × 10−5	2.35777169	0.021	5.68 × 10−6	160	2.35773897
	22	2.35777129	0.034	8.08 × 10−6	2.35777165	0.022	7.72 × 10−6	320	2.35776357
	24	2.35777151	0.035	2.01 × 10−6	2.35777165	0.023	1.87 × 10−6	640	2.35776964
	26	2.35777167	0.036	5 × 10−7	2.35777165	0.024	2 × 10−8	1280	2.35777115
5	20	2.41421578	0.033	1.743 × 10−5	2.41420169	0.021	3.34 × 10−6	50	2.41419835	2.41420167
	21	2.41420238	0.033	1.37 × 10−6	2.41420169	0.021	6.8 × 10−7	100	2.41420101
	22	2.41420214	0.034	6.2 × 10−7	2.41420168	0.022	1.6 × 10−7	200	2.41420152
	23	2.41420177	0.034	1.4 × 10−7	2.41420167	0.022	4 × 10−8	400	2.41420163
	24	2.41420175	0.035	9 × 10−8	2.41420167	0.023	1 × 10−8	800	2.41420166
	25	2.41420171	0.035	4 × 10−8	2.41420167	0.023	1 × 10−8	1600	2.41420167
	26	2.41420169	0.036	2 × 10−8	2.41420167	0.024	1 × 10−8	3200	2.41420167

. Hence, it is found that non-uniform grid results are higher and more consistent with earlier solutions than uniform ones with an error ≤${10}^{-8}$, and execution time of about (0.024 s). Also, when the grid points increase with time, the accuracy increases; for example, at time (t = 1), we make (N = 13), and at time (t = 2), we make (N = 26) Furthermore, the maximum number of grids we use is significantly less than in previous studies (N = 3200).

Table 2. Computation of $\mathsf{\upsilon }\left(\mathrm{t}\right)$ via non-uniform PDQM and DSCDQM–RSK for fractional Riccati equation with time (T = 1) and fraction $\left(\mathsf{\alpha }=1, \mathsf{\rho }=1\right)$ at various grid points (N), regularized Shannon factor $\left(\mathsf{\sigma }=\mathrm{h}\times \Delta \right)$, and bandwidth (2K + 1).

N	2K + 1	Non-Uniform PDQM	DSCDQM–RSK
			$\mathsf{\sigma }=0.2\times \Delta$	$\mathsf{\sigma }=0.4\times \Delta$	$\mathsf{\sigma }=0.45\times \Delta$	$\mathsf{\sigma }=0.5\times \Delta$	CPU (s)
9	3	1.68948616	1.69745751	1.689880214	1.68950741	1.64023989	0.008
	5	1.68948616	1.69745647	1.689876547	1.68949956	1.64023942	0.01
	7	1.68948616	1.69745559	1.689875120	1.68949839	1.64023865	0.012
	9	1.68948616	1.69745559	1.689875120	1.68949839	1.64023865	0.014
	11	1.68948616	1.69745559	1.689875120	1.68949839	1.64023865	0.016
11	3	1.68949815	1.69745666	1.689877415	1.68949951	1.64023937	0.009
	5	1.68949815	1.69745559	1.689875120	1.68949839	1.64023865	0.01
	7	1.68949815	1.69745559	1.689875120	1.68949839	1.64023865	0.011
	9	1.68949815	1.69745559	1.689875120	1.68949839	1.64023865	0.012
	11	1.68949815	1.69745559	1.689875120	1.68949839	1.64023865	0.014
13	3	1.68949839	1.69745578	1.689876014	1.6894990	1.64023900	0.01
	5	1.68949839	1.69745559	1.689875120	1.68949839	1.64023865	0.012
	7	1.68949839	1.69745559	1.689875120	1.68949839	1.64023865	0.014
	9	1.68949839	1.69745559	1.689875120	1.68949839	1.64023865	0.016
	11	1.68949839	1.69745559	1.689875120	1.68949839	1.64023865	0.018
Earlier numerical solutions [64]			1.68949786           at (N = 640)
Exact [61]			1.68949839

compares non-uniform PDQM and DSCDQM–RSK for the fractional Riccati equation under various conditions of time (T = 1), fraction $\left(\mathsf{\alpha }=1, \mathsf{\rho }=1\right)$, different grid points (N), regularized Shannon factor $\left(\mathsf{\sigma }=\mathrm{h}\times \Delta \right)$, and bandwidth (2K + 1). To begin, Table 2 ensures that the best value of the regularized Shannon factor is $\mathsf{\sigma }=0.45\times \Delta$, with results matching previous studies and PDQM shown in Table 1 at the fewest grid points (N = 9), bandwidth (2K + 1 = 7) and performance time of about (0.018 s). DSCDQM–RSK is the best method overall, based on low grid points and performance time when compared to PDQM (N = 16) and previous studies (N = 640).

The efficiency, validity, and accuracy of the created methodologies are presently being explored by comparing the calculated results to earlier numerical and analytical solutions at various powers of fraction $\left(\mathsf{\alpha },\mathsf{\rho }\right)$, as shown in

Table 3. Computation of $\mathsf{\upsilon }\left(t\right)$ via non-uniform PDQM and DSCDQM–RSK for fractional Riccati equation at various times (T), fractions $\left(\mathsf{\alpha }\right)$, and $\mathsf{\rho }=1$.

T	Non-Uniform PDQM		DSCDQM–RSK		Previous Results [61,62,63,64,65]
	α = 0.5	α = 0.75	α = 0.5	α = 0.75	α = 0.5	α = 0.75
0.1	0.59149371	0.24512554	0.59149373	0.24512556	0.59149373	0.24512556
0.2	0.93141486	0.47450194	0.93141488	0.47450196	0.93141488	0.47450196
0.3	1.171926469	0.709154008	1.171926471	0.709154010	1.171926471	0.709154010
0.4	1.344407759	0.937441050	1.344407761	0.937441052	1.344407761	0.937441052
0.5	1.471501155	1.147807349	1.471501157	1.147807351	1.471501157	1.147807351
0.6	1.568070430	1.332985034	1.568070432	1.332985036	1.568070432	1.332985036
0.7	1.643596282	1.490535187	1.643596284	1.490535189	1.643596284	1.490535189
0.8	1.704182955	1.621592245	1.704182957	1.621592247	1.704182957	1.621592247
0.9	1.753855962	1.729220378	1.753855964	1.729220380	1.753855964	1.729220380
1	1.817133594	1.795344168	1.817133596	1.795344170	1.817133596	1.795344170

,

Table 4. Computation of $\mathsf{\upsilon }\left(t\right)$ via non-uniform PDQM and DSCDQM–RSK for fractional Riccati equation with fraction $\left(\mathsf{\alpha }=1\right)$ at various times (T), and fractions $\left(\mathsf{\rho }\right)$.

T	Non-Uniform PDQM		DSCDQM–RSK		Previous Results [61,62,63,64,65]
	ρ = 0.8	ρ = 1.2	ρ = 0.8	ρ = 1.2	ρ = 0.8	ρ = 1.2
0.1	0.14117992	0.09045268	0.14117994	0.09045271	0.14117994	0.09045271
0.2	0.31592641	0.195667845	0.31592645	0.19566787	0.31592645	0.19566787
0.3	0.52298485	0.315926409	0.52298488	0.315926411	0.52298488	0.315926411
0.4	0.75601439	0.450653813	0.75601442	0.450653816	0.75601442	0.450653816
0.5	1.00354951	0.59824597	1.00354953	0.59824599	1.00354953	0.59824599
0.6	1.25086733	0.75601439	1.25086736	0.75601442	1.25086736	0.75601442
0.7	1.48329584	0.92030072	1.48329586	0.92030075	1.48329586	0.92030075
0.8	1.68949839	1.08677371	1.68949842	1.08677374	1.68949842	1.08677374
0.9	1.86328744	1.25086733	1.86328746	1.25086736	1.86328746	1.25086736
1	2.00353694	1.40827080	2.00353696	1.40827081	2.00353696	1.40827081

,

Table 5. Computation of $\mathsf{\upsilon }\left(t\right)$ via non-uniform PDQM for fractional Riccati equation at different grid points (N), and fractions $\left(\mathsf{\alpha }, \mathsf{\rho }\right)$.

T	Non-Uniform				Earlier Numerical Solutions [62]
	N	α = 1, ρ = 0.9	α = 0.95, ρ = 0.75	α = 0.9, ρ = 1.2	N	α = 1, ρ = 0.9	α = 0.95, ρ = 0.75	α = 0.9, ρ = 1.2
1	4	1.80602802	1.96263496	1.39368786	10	1.84281224	2.06729863	1.52944766
	6	1.84319793	2.04896144	1.49050147	20	1.84491385	2.07202706	1.53119172
	8	1.84556010	2.06510470	1.51019655	40	1.84546411	2.07322261	1.53167452
	9	1.84561319	2.06904164	1.52085146	80	1.84560424	2.07352741	1.53180584
	11	1.84565025	2.07362649	1.52696219	160	1.84563955	2.07360571	1.53184129
	12	1.84565137	2.07363589	1.52945738	320	1.84564841	2.07362592	1.53185082
	13	1.84565137	2.07363256	1.53185408	640	1.84565063	2.07363115	1.53185339
	14	1.84565137	2.07363256	1.53185408	1280	1.84565119	2.07363250	1.53185407
2	8	2.28430436	2.32839830	2.20693947	20	2.36576348	2.34646084	2.26631061
	12	2.32381757	2.33437164	2.21000463	40	2.36763874	2.34834846	2.26840179
	16	2.34381756	2.34307946	2.22083045	80	2.36805246	2.34876916	2.26890814
	18	2.36382620	2.34836151	2.23103727	160	2.36815011	2.34887135	2.26903810
	22	2.36832617	2.34863032	2.24089052	320	2.36817385	2.34889710	2.26907235
	24	2.36818255	2.34889017	2.26879047	640	2.36817971	2.34890369	2.26908148
	26	2.36818153	2.34890584	2.26908459	1280	2.36818116	2.34890540	2.26908393
	27	2.36818153	2.34890584	2.26908459	2560	2.36818153	2.34890584	2.26908459

and

Table 6. Computation of $\mathsf{\upsilon }\left(t\right)$ via DSCDQM–RSK for fractional Riccati equation at different values of $\left(\mathsf{\alpha }, \mathsf{\rho }\right)$ and times.

α	ρ	T = 0.5	T = 1	T = 2	T = 2.5	T = 3
0.4	1	1.58967600	1.77525996	2.09537035	2.21247852	2.22642314
	1.1	1.54739283	1.74403298	2.07553739	2.19214756	2.20554878
	1.2	1.50756398	1.71405017	2.05317313	2.17021456	2.18875143
	1.4	1.43438179	1.65736347	2.02666948	2.14654879	2.16214787
	1.9	1.28182554	1.53134020	2.01520768	2.12958092	2.13478462
0.5	1	1.68300421	1.79935747	2.12065899	2.25983372	2.27664509
	1.1	1.66554391	1.75798866	2.10285480	2.23984120	2.25471201
	1.2	1.65154233	1.71762640	2.09852919	2.21874621	2.23789123
	1.4	1.62320069	1.66991063	2.07719592	2.20997411	2.21987423
	1.9	1.58526142	1.56342314	2.05661402	2.20278414	2.17645789
0.7	1	1.71613371	1.82286926	2.14596688	2.29139917	2.31617003
	1.1	1.69532577	1.79031758	2.12529159	2.27075056	2.29157030
	1.2	1.67462609	1.77811144	2.10299487	2.25941935	2.26895529
	1.4	0.65934200	1.73806497	2.08048316	2.22364828	2.24318890
	1.9	0.60930946	1.65504250	2.05990416	2.20861195	2.20862313
0.85	1	0.73587574	1.84052325	2.24870001	2.30150738	2.33501568
	1.1	0.71892188	1.82175475	2.21842773	2.28982006	2.31048992
	1.2	0.69728893	1.80614145	2.18322274	2.27519396	2.29534657
	1.4	0.67897490	1.77168741	2.16870824	2.23587987	2.27096577
	1.9	0.62629628	1.72273733	2.13184852	2.20356265	2.22777374
0.95	1	0.82414933	1.87127626	2.31830946	2.36707298	2.37441750
	1.1	0.73300489	1.84609224	2.28379937	2.35373760	2.35930723
	1.2	0.71834905	1.82980945	2.24083446	2.33559233	2.33776914
	1.4	0.68447266	1.80982469	2.21190242	2.28289347	2.31913452
	1.9	0.63618443	1.77131197	2.19457437	2.23111789	2.29258450

. Table 3, Table 4, Table 5 and Table 6 show that increasing the fraction power $\left(\mathsf{\alpha } \mathrm{or} \mathsf{\rho }\right)$ decreases the value of $\mathsf{\upsilon }\left(t\right)$, but increases with time for the fractional Riccati equation. Furthermore, the results show that DSCDQM–RSK outperforms non-uniform PDQM in terms of efficiency, validity, and accuracy.

In addition, the dynamic behaviors of the fractional Riccati equation with respect to the parameters $\left(\mathsf{\alpha } \mathrm{and} \mathsf{\rho }\right)$ and against the time variable t are depicted in Figure 1 and Figure 2.

4.2. Problem 4.2

We deal with the fractional Lorenz system after substituting Equations (4)–(6) with Equations (19) and (20) of the proposed methods:

$$\sum _{j=1}^L{a}_{ij}^{\mathsf{\alpha },\mathsf{\rho }} \mathrm{X}\left(t_j\right)=\mathsf{\lambda }\left[\sum _{j=1}^N{\mathsf{\delta }}_{\mathrm{ij}} \mathrm{Y}\left(t_j\right)-\sum _{j=1}^N{\mathsf{\delta }}_{\mathrm{ij}} \mathrm{X}\left(t_j\right)\right]$$

(24)

$$\begin{array}{l}\sum _{j=1}^L{a}_{ij}^{\mathsf{\alpha },\mathsf{\rho }} \mathrm{Y}\left(t_j\right)=\left(\mathsf{\varphi }-\mathsf{\lambda }\right)\sum _{j=1}^N{\mathsf{\delta }}_{\mathrm{ij}} \mathrm{X}\left(t_j\right)-\\ \sum _{j=1}^N{\mathsf{\delta }}_{\mathrm{ij}} \mathrm{X}\left(t_j\right)\sum _{j=1}^N{\mathsf{\delta }}_{\mathrm{ij}} \mathrm{Z}\left(t_j\right)+\mathsf{\varphi }\sum _{j=1}^N{\mathsf{\delta }}_{\mathrm{ij}} \mathrm{Y}\left(t_j\right)\end{array}$$

(25)

$$\sum _{j=1}^L{a}_{ij}^{\mathsf{\alpha },\mathsf{\rho }} \mathrm{Z}\left(t_j\right)=\sum _{j=1}^N{\mathsf{\delta }}_{\mathrm{ij}} \mathrm{X}\left(t_j\right)\sum _{j=1}^N{\mathsf{\delta }}_{\mathrm{ij}} \mathrm{Y}\left(t_j\right)-\mathsf{\beta }\sum _{j=1}^N{\mathsf{\delta }}_{\mathrm{ij}} \mathrm{Z}\left(t_j\right)$$

(26)

Dealing with the initial condition (7) is also done by substituting in the governing Equations (24)–(26). After that, we use Equation (22) to solve this system.

We will now begin to demonstrate the obtained results to explain the stability, reliability, convergence, and performance of FDQM using two types of shape functions with generalized Caputo sense, as follows:

Table 7. Computation of numerical solutions $X,Y,Z$ via uniform and non-uniform PDQM for fractional Chen system (Lorenz system) with time (T = 2) and fraction $\left(\mathsf{\alpha }=1, \mathsf{\rho }=1\right)$ at various grid points (N). $\left({\mathrm{x}}_0=2.5,{\mathrm{y}}_0=1,{\mathrm{z}}_0=0.5\right)$ and $\left(\mathsf{\lambda }=1,\mathsf{\beta }=2,\mathsf{\varphi }=1\right)$.

N	Uniform PDQM			CPU (s)	Non-uniform PDQM			CPU (s)
	X	Y	Z		X	Y	Z
4	0.837919933	0.623902547	0.302407988	0.02	0.837919846	0.623902388	0.302407749	0.02
5	0.771720200	0.527485798	0.248194100	0.021	0.771720184	0.527485752	0.248193998	0.021
6	0.761378603	0.49832669	0.2415685802	0.022	0.763378566	0.522326584	0.2475685753	0.022
7	0.760223666	0.497601296	0.2409074397	0.023	0.762223591	0.515601283	0.2469074378	0.023
9	0.762203989	0.500122597	0.242400811	0.025	0.762203974	0.500122578	0.2424007951	0.024
10	0.762216222	0.500169299	0.2424081878	0.026	0.762216161	0.500169278	0.2424081821	0.025
11	0.76221575	0.500167396	0.242407838	0.027	0.76221572	0.500167392	0.242407832	0.025
12	0.76221573	0.500167394	0.2424078434	0.028	0.76221572	0.500167392	0.2424078432	0.026
13	0.76221572	0.500167392	0.2424078432	0.029	0.76221572	0.500167392	0.2424078432	0.027
	Earlier numerical solutions [62,63]
	X		Y		Z
1280	0.76221649		0.50016919		0.24240833
	RK4 [64]
	0.76221572		0.50016739		0.24240783

explains the impact of grid points (N) on the obtained results X, Y, and Z via PDQM with uniform and non-uniform grid distributions of the fractional Lorenz system with fraction $\left(\mathsf{\alpha }=1, \mathsf{\rho }=1\right)$. It is remarkable that non-uniform grid results are higher and more consistent with earlier solutions [62,63] and RK4 [64] at N = 11 than uniform ones at N = 13 with error ≤${10}^{-8}$, and execution time of about (0.027 s).

Table 8. Computation of numerical solutions $X,Y,Z$ via DSCDQM–RSK for fractional Chen system (Lorenz system) with time (T = 2) and fraction $\left(\mathsf{\alpha }=1, \mathsf{\rho }=1\right)$ at various grid points (N), regularized Shannon factor $\left(\mathsf{\sigma }=\mathrm{h}\times \Delta \right)$, and bandwidth (2K + 1). $\left({\mathrm{x}}_0=2.5,{\mathrm{y}}_0=1,{\mathrm{z}}_0=0.5\right)$ and $\left(\mathsf{\lambda }=1,\mathsf{\beta }=2,\mathsf{\varphi }=1\right)$.

N	2K + 1	DSCDQM–RSK									CPU (s)
		σ = 0.4 × Δ			σ = 0.45 × Δ			σ = 0.47 × Δ
		X	Y	Z	X	Y	z	X	Y	Z
9	3	0.7738	0.5133	0.2549	0.7625	0.5008	0.2430	0.76225	0.5004	0.2425	0.01
	5	0.7735	0.5127	0.2540	0.7623	0.5006	0.2427	0.76223	0.5003	0.2424	0.012
	7	0.7732	0.5122	0.2536	0.76229	0.5005	0.2426	0.76221	0.5002	0.2424	0.014
	9	0.7732	0.5122	0.2536	0.76229	0.5005	0.2426	0.76221	0.5002	0.2424	0.016
	11	0.7732	0.5122	0.2536	0.76229	0.5005	0.2426	0.76221	0.5002	0.2424	0.018
11	3	0.7347	0.5127	0.2540	0.7623	0.5006	0.2427	0.76222	0.50023	0.2425	0.012
	5	0.734	0.5124	0.2537	0.76229	0.5005	0.2426	0.76221	0.5002	0.2424	0.014
	7	0.7732	0.5122	0.2536	0.76229	0.5005	0.2426	0.76221	0.5002	0.2424	0.016
	9	0.7732	0.5122	0.2536	0.76229	0.5005	0.2426	0.76221	0.5002	0.2424	0.018
	11	0.7732	0.5122	0.2536	0.76229	0.5005	0.2426	0.76221	0.5002	0.2424	0.020
13	3	0.7733	0.5122	0.2537	0.76231	0.5006	0.2427	0.762215	0.50026	0.24247	0.014
	5	0.7732	0.5122	0.2536	0.76229	0.5005	0.2426	0.76221	0.5002	0.2424	0.016
	7	0.7732	0.5122	0.2536	0.76229	0.5005	0.2426	0.76221	0.5002	0.2424	0.018
	9	0.7732	0.5122	0.2536	0.76229	0.5005	0.2426	0.76221	0.5002	0.2424	0.020
	11	0.7732	0.5122	0.2536	0.76229	0.5005	0.2426	0.76221	0.5002	0.2424	0.022
	Earlier numerical solutions [62,63]
	X			Y				Z
1280	0.76221649			0.50016919				0.24240833
	RK4 [64]
	0.76221572			0.50016739				0.24240783

demonstrates the effect of control values like grid points (N), regularized Shannon factor $\left(\mathsf{\sigma }=h\times \Delta \right)$, and bandwidth (2K + 1) on the obtained results by DSCDQM–RSK at time (T = 2), fraction $\left(\mathsf{\alpha }=1, \mathsf{\rho }=1\right)$. We found the best value of the regularized Shannon factor is $\mathsf{\sigma }=0.47\times \Delta$, with results matching previous studies and PDQM shown in Table 7 at the fewest grid points (N = 9), bandwidth (2K + 1 = 3), and performance time of about (0.022 s). Furthermore, from Table 7 and Table 8, it is noted that the maximum number of grids we use is significantly less than in previous studies (N = 1280).

Figure 3 and Figure 4 present the influence of time and fraction $\mathsf{\alpha }$ on the numerical results X, Y, and Z via non-uniform PDQM and DSCDQM–RSK at $\left({\mathrm{x}}_0=2.5,{\mathrm{y}}_0=1,{\mathrm{z}}_0=0.5\right)$ and $\left(\mathsf{\lambda }=1,\mathsf{\beta }=2,\mathsf{\varphi }=1\right)$. Thus, it is found that the dynamic behaviors of X, Y, and Z differ when the fraction $\mathsf{\alpha }$ change. This means that when the value $\mathsf{\alpha }$ decreases, the values of X, Y, and Z increase.

Also, Figure 3 and Figure 4 show that the Lorenz system, a classic model for chaotic dynamics, can be extended to fractional-order systems. This extension introduces significant differences in the observed chaotic behavior [74,75]:

Chaotic Regimes:

▪

The integer-order Lorenz system exhibits a well-defined chaotic regime within a specific parameter range.

▪

Fractional-order Lorenz systems often exhibit chaotic behavior over a wider range of parameters and fractional orders. This can lead to more complex and diverse chaotic dynamics.

Attractor Structure:

▪

The integer-order Lorenz system typically has a single strange attractor.

▪

Fractional-order Lorenz systems can exhibit multiple strange attractors or even the coexistence of different attractors, depending on the fractional order and system parameters.

Fractal Dimension:

▪

The fractal dimension of the strange attractor in the integer-order Lorenz system is generally between 2 and 3.

▪

The fractal dimension of the strange attractors in fractional-order Lorenz systems can vary more widely, often exceeding 3. This indicates a more complex and convoluted structure.

Sensitivity to Initial Conditions:

▪

The integer-order Lorenz system is highly sensitive to initial conditions, leading to the butterfly effect.

▪

Fractional-order Lorenz systems can exhibit even greater sensitivity to initial conditions, making long-term predictions even more challenging.

Memory Effects:

▪

The integer-order Lorenz system does not have memory effects.

▪

The fractional-order Lorenz system incorporates memory effects, which can influence the system’s dynamics and make it more resilient to perturbations.

So, fractional-order Lorenz systems can exhibit more complex and diverse chaotic behaviors compared to their integer-order counterparts. The increased sensitivity to initial conditions in fractional-order systems makes long-term predictions even more challenging. Fractional-order Lorenz systems can be used to model real-world phenomena with memory effects or nonlinearities that are not adequately captured by integer-order models.

Figure 4, Figure 5, Figure 6 and Figure 7 show the fractional Lorenz system’s dynamic behaviors as the values of fractions $\left(\mathsf{\alpha }, \mathsf{\rho }\right)$ at $\left({\mathrm{x}}_0=2.5,{\mathrm{y}}_0=1,{\mathrm{z}}_0=0.5\right)$ and $\left(\mathsf{\lambda }=1,\mathsf{\beta }=2,\mathsf{\varphi }=1\right)$ changes. Also, in these figures, we show projections of the fractional Lorenz system attractors calculated via DSCDQM–RSK when T = 100 for some values of the fractions $\left(\mathsf{\alpha }, \mathsf{\rho }\right)$. It is observed that when the fraction $\mathsf{\rho }$ increases, the chaotic behavior increases more than the fraction $\mathsf{\alpha }$ changes. Consequently, it is noted that the fractional Lorenz system may exhibit chaotic attractors similar to those of its integer-order counterpart when $\left(\mathsf{\alpha }=0.9, \mathsf{\rho }=0.8\right)$ and $\left(\mathsf{\alpha }=0.8, \mathsf{\rho }=1.2\right)$. Also, for smaller values of the fractions $\left(\mathsf{\alpha }, \mathsf{\rho }\right)$ the system loses its chaotic character.

To analyze the computational complexity of the provided code for solving the Riccati equation and the Lorenz system using the Fractional Differential Quadrature Method (FDQM), we will focus on memory space and simulation time complexities in detail. The configuration of the computer used to perform the simulation results is HP Probook 450 G8 Laptop—11th Intel Core i5-1135G7, 8 GB RAM, 512 GB PCIe NVMe SSD, 15.6” FHD (1920 × 1080), and Intel Iris X Graphics.

1. Memory Space Complexity

Variables:

▪

Grid Points: The function Chebyshev grid (N) generates N Chebyshev nodes, requiring O(N) space.

The DSCDQM–RSK method demonstrated optimal performance with N = 9 grid points, a bandwidth of 2k + 1 = 7, and a regularized Shannon factor of $\mathsf{\sigma }=0.45\times \Delta$ for the fractional Riccati equation. This configuration yielded a CPU time of 0.018 s. Similarly, for the fractional Lorenz system, the best results were obtained with N = 9 grid points, a bandwidth of 2k + 1 = 3, and a regularized Shannon factor of $\mathsf{\sigma }=0.47\times \Delta$, achieving a CPU time of 0.02 s.

In contrast, earlier numerical methods typically required significantly more grid points, often reaching N = 640, 1280, or even 3200.

▪

Weighting Coefficients: The function PDQM weights (N, t) creates a matrix A of size N × N. Therefore, it requires O(N²) space.

▪

Solution Vectors:

-

The solution υ, X, Y, Z vectorsfor the Riccati equation and Lorenz system are initialized as zero vectors of size N × 1, each requiring O(N) space.

-

Total for all solution vectors combined: 4 × O(N) = O(N).

Total Memory Space Complexity:

▪

The dominant term is O(N²) from the weighting coefficients matrix. Thus, the total memory space complexity is:

O(N²)

2. Simulation Time Complexity

For Solving the Riccati Equation:

Grid Point Generation:

▪

The grid points are generated in O(N) time.

Weighting Coefficients Calculation:

▪

The PDQM weights(N, t, 1) function computes the coefficients with a nested loop over N, resulting in O(N²) time complexity:

▪

Each entry in matrix A involves calculations that depend on N, leading to O(N²) complexity for the entire matrix.

Time Integration:

▪

The time integration loop runs for N-1 iterations, performing a constant time calculation for each iteration:

O(N)

Total Time Complexity for Riccati Equation:

▪

Combining these, we get:

O(N²) (weighting coefficients) + O(N) (grid points) + O(N) (integration) = O(N²)

For Solving the Lorenz System:

Grid Point Generation:

▪

Again, this takes O(N) time.

Weighting Coefficients Calculation:

▪

The calculations for Ax, Ay, and Az each take O(N²):

3 × O(N²) = O(N²)

Time Integration:

▪

Similar to the Riccati equation, the integration loop runs for N-1 iterations:

O(N)

Total Time Complexity for Lorenz System:

▪

Again combining these:

O(N²) (coefficients) + O(N) (grid points) + O(N) (integration) = O(N²)

3. Overall Complexity Summary

▪

Memory Complexity: O(N²)

▪

Time Complexity for Riccati Equation: O(N²)

▪

Time Complexity for Lorenz System: O(N²)

4. Real Numbers Example

For practical evaluation, consider the following:

▪

For N = 9:

Memory for weighting coefficients: 81 entries.

Memory for solution vectors: 36 entries.

Assuming each entry takes 8 bytes (for double precision), the memory usage would be:

▪

Weighting coefficients: 81 × 8 = 648 bytes (approximately 0.64 KB).

▪

Solution vectors: 36 × 8 = 288 bytes (approximately 0.28 KB).

Total memory usage for N = 9 would be approximately 1 KB.

Execution Time:

▪

If the operations in the loops take, say, 0.001 s per iteration:

For N = 9, the time for solving the Riccati equation and the Lorenz system would be dominated by the O(N²) term, leading to an estimated execution time of about 0.1 s (for illustration).

This detailed analysis provides insights into the computational complexity of the code, which is crucial for assessing performance in practical scenarios. The findings underscore the importance of optimizing the weighting coefficients and the iterative solvers for larger values of N.

The choice of time step in numerical methods significantly affects the accuracy of the solutions for differential equations, including those solved using the Differential Quadrature Method (DQM). Here’s how:

5. Stability and Convergence:

Our numerical methods are stable because small perturbations in the initial conditions lead to small perturbations in the numerical solution.

Consider the fractional differential equation (FDE) represented in the form:

$$D_c^{\alpha }U\left(t\right)=f(t,U(t\left)\right),t\in [0,T]$$

(27)

where $D_c^{\alpha }$ is the generalized Caputo fractional derivative.

Let $U\left(t\right)$ be the exact solution and $U_n\left(t\right)$ be the numerical solution obtained through the proposed methods. Introduce a perturbation ${ϵ}_n$ such that:

$$U_n\left(t\right)=U\left(t\right)+{ϵ}_n$$

(28)

By analyzing how the perturbation evolves over time:

$$D_c^{\alpha }\left(U_n\left(t\right)+{ϵ}_n\right)=f\left(t,U_n\left(t\right)+{ϵ}_n\right)$$

(29)

This leads to the error equation:

$$D_c^{\alpha }{ϵ}_n=f\left(t,U_n\left(t\right)+{ϵ}_n\right)-f(t,U(t\left)\right)$$

(30)

Assume $f$ satisfies a Lipschitz condition:

$$\left|f\left(t,U_1\right)-f\left(t,U_2\right)\right|\le L\left|U_1-U_2\right|$$

(31)

where $L$ is a constant.

By applying Gronwall’s inequality, show that:

$$\left|{ϵ}_n\right|\le C\cdot \Delta t^p$$

where $C$ is a constant and $p$ is the order of the method, ensuring that the solution remains bounded as $n\to \mathrm{\infty }$.

Our numerical method converges because the numerical solution approaches the exact solution as the grid refinement increases.

Let:

$$\underset{N\to \mathrm{\infty }}{\mathrm{l}\mathrm{i}\mathrm{m}}∥U_n\left(t\right)-U\left(t\right)∥=0$$

(32)

By conducting an error analysis between the numerical solution and the exact

$$E_n=U_n\left(t\right)-U\left(t\right)$$

(33)

By using Taylor expansion around $t_n$ to express $U\left(t\right)$:

$$U\left(t_n+\Delta t\right)=U\left(t_n\right)+\Delta t{U}^{\prime }\left(t_n\right)+O\left(\Delta t^2\right)$$

(34)

Relate this to the discretized version derived from the method:

$$U_n\left(t_n+\Delta t\right)=U_n\left(t_n\right)+O\left(\Delta t^p\right)$$

(35)

The method is stable and the achieved convergence rates reached 10⁻⁸, indicating a high level of precision in solving the nonlinear fractional initial value problems. The truncation error in our numerical methods is very small due to the Gaussian regularizer $\mathsf{\sigma }=0.45\times \Delta$ which depends on a small computational domain:

6. Accuracy of the Solution:

▪

Each step introduces local error, which accumulates over time. Smaller time steps help minimize this accumulation, resulting in a more accurate final solution (error ≤10⁻⁸)

▪

The global error, which is the total error over the entire integration period, also tends to decrease with smaller time steps, leading to better overall accuracy.

4.3. Stability Analysis

After applying our discretization schemes to Equations (1)–(6), we obtained an equivalent set of ordinary differential equations in the time domain:

$$\frac{\mathrm{d}\left[U\right]}{\mathrm{d}\mathrm{t}}=\mathrm{R}\left[U\right]+\left[\mathrm{K}\right]$$

(36)

where

• The vector {U} represents the unknown variables at the internal grid points, where υ, X, Y, and Z are the individual components;

  •

  The initial conditions are stored in the vector [K];

• R[U] is the right-hand side of Equations (1)∓(6); and
• $a_{ij}^{\left(1\right)}$ is the weighting coefficient matrix of the first derivative:

$$a_{ij}^{\left(1\right)}={\left[\begin{array}{ccc}\begin{array}{cc}{a}_{22}^{\left(1\right)}& a_{23}^{\left(1\right)}\\ a_{32}^{\left(1\right)}& a_{33}^{\left(1\right)}\end{array}& \cdots & \begin{array}{c}{a}_{2(n-1)}^{\left(1\right)}\\ a_{3(n-1)}^{\left(1\right)}\end{array}\\ ⋮& \ddots & ⋮\\ \begin{array}{cc}{a}_{(n-1)2}^{\left(1\right)}& a_{(n-1)3}^{\left(1\right)}\end{array}& \cdots & a_{(n-1)(n-1)}^{\left(1\right)}\end{array}\right]}_{\left(N-2\right)\times (N-2)}$$

(37)

The stability of our technique was evaluated by examining system (35). We employed eigenvalue analysis of the coefficient matrices $\left(a\right)$ to determine stability.

For the fractional differential equation represented in the form:

$$D_c^{\alpha }U\left(t\right)=f(t,U(t\left)\right)$$

discretizing this using the FDQM leads to a system of equations expressed as:

$$\left[a\right]\left[\mathrm{U}\right]=\mathrm{R}\left[U\right]$$

To analyze stability, we need to compute the eigenvalues $\lambda$ of the matrix $\left[a\right]$. The eigenvalue problem is given by:

$$\left[a\right]\left[\mathrm{v}\right]=\lambda \left[\mathrm{v}\right]$$

where $\mathrm{v}$ is the eigenvector associated with the eigenvalue $\lambda$.

The characteristic polynomial is obtained from:

$$\mathrm{d}\mathrm{e}\mathrm{t}(a-\mathsf{\lambda }I)=0$$

where $I$ is the identity matrix.

The stability of the numerical method hinges on the eigenvalues of $\left[a\right]$:

• If all eigenvalues satisfy $\left|\lambda \right|<1$, the method is stable.
• If any eigenvalue has $\left|\lambda \right|\ge 1$, the method may be unstable.

Figure 8 displays the stability regions for a numerical method at various fractional orders (α). Each subplot represents a different α value: 0.5, 0.7, 0.85, and 1. The plots show the eigenvalues of the coefficient matrix ($a$) in the complex plane (Real vs. Imaginary). At α = 0.5, the eigenvalues are all located on the real axis near −0.5. Since the magnitude of these eigenvalues is less than 1, this indicates stability for this fractional order. At α = 0.7, the eigenvalues form a V-shape, extending into both the positive and negative real axis. A portion of the eigenvalues have magnitudes greater than 1 (outside the unit circle), indicating instability in this region. The region of stability is limited to the portion of the V-shape within the unit circle. At α = 0.85, the eigenvalues form an inverted V-shape. Similar to (b), portions of the eigenvalues are outside the unit circle, indicating instability. The region of stability is again limited to the portion within the unit circle. At α = 1, the eigenvalues lie entirely on the negative real axis, forming a vertical line. All eigenvalues appear to be within the unit circle, suggesting stability for this fractional order (which corresponds to the standard integer-order case). Furthermore, Figure 9 depicts the error propagation in relation to time and fractional order.

5. Conclusions

In this present work, we have successfully investigated new numerical methods for solving the fractional Riccati equation and fractional Lorenz system. The novel numerical method is FDQM, which is based on two base functions: Lagrange interpolation polynomial and discrete singular convolution-Regularized Shannon kernel with a new generalized Caputo kind. These methods are used to transform the proposed problems into a nonlinear algebraic system. Then, the iterative method is employed to deal with the problem of nonlinearity. All numerical results were obtained using MATLAB. By comparing our results with those of existing methods, we demonstrated the superior accuracy, efficiency, and overall performance of our proposed techniques. The achieved convergence rates reached 10⁻⁸, indicating a high level of precision in solving the nonlinear fractional initial value problems. Error analysis showed that non-uniform grid distributions consistently outperformed uniform distributions, with maximum errors diminishing significantly as grid points increased. Our numerical results demonstrate that the DSC-RSK method achieved significantly faster convergence rates compared to other techniques. The best results of the DSC-RSK method are achieved when grid points are N = 9, bandwidth is 2k + 1 = 7, and the regularized Shannon factor is $\mathsf{\sigma }=0.45\times \Delta$ at CPU time = 0.018 s for the fractional Riccati equation. But for the fractional Lorenz system, the best results are when grid points are N = 9, bandwidth is 2k + 1 = 3, and the regularized Shannon factor is $\mathsf{\sigma }=0.47\times \Delta$ at CPU time = 0.02 s. Also, the proposed techniques have been successfully employed to explain the fractional systems’ dynamic behaviors. The numerical results demonstrate a strong dependence of the solution on the fractional derivative. The fractional parameters, $\mathsf{\alpha } \mathrm{and} \mathsf{\rho }$, offer significant advantages in studying the proposed problems with greater accuracy compared to traditional approaches. These techniques hold promise for solving more complex nonlinear equations and other differential applications involving fractional derivatives. The versatility of the FDQM was highlighted through its successful application to both the fractional Riccati equations and the fractional Lorenz system, demonstrating its potential for broader applications in fields requiring the modeling of complex dynamics. The findings suggest that further exploration of fractional orders could uncover even richer dynamics, with potential applications extending to areas such as control theory, physics, and finance.