High-Dimensional Chaotic Lorenz System: Numerical Treatment Using Changhee Polynomials of the Appell Type

Abstract

Presenting and simulating the numerical treatment of the nine-dimensional fractional chaotic Lorenz system is the goal of this work. The spectral collocation method (SCM), which makes use of Changhee polynomials of the Appell type, is the suggested approximation technique to achieve this goal. A rough formula for the Caputo fractional derivative is first derived, and it is used to build the numerical strategy for the suggested model’s solution. This procedure creates a system of algebraic equations from the model that was provided. We validate the effectiveness and precision of the provided approach by evaluating the residual error function (REF). We compare the results obtained with the fourth-order Runge–Kutta technique and other existing published work. The outcomes demonstrate that the technique used is a simple and effective tool for simulating such models.

1. Introduction

In 1963, chaotic behavior was discovered in differential equations representing meteorological phenomena [1], and chaos as a phenomenon was born. Since then, numerous applications in science and engineering have given rise to nonlinear chaotic and hyperchaotic systems. The positive Lyapunov exponents of the system are inversely correlated with the complexity of the chaotic behavior. A chaotic system essentially has one positive Lyapunov exponent; however, hyperchaotic systems are believed to have more than one positive Lyapunov exponent. When Rossler [2] discovered hyperchaotic behavior in differential equations for simulating chemical reactions, he proposed the idea of hyperchaos. Numerous chaotic and hyperchaotic systems developing in other disciplines have been researched since these two findings [3]. In comparison to regular chaotic systems, hyperchaotic systems typically exhibit more complex dynamical behaviors. Strong sensitivity to initial conditions and rapidly changing outcomes are traits of chaotic systems. Chaotic dynamical systems are challenging to quantitatively investigate due to these properties. Systems of ordinary differential equations (ODEs) displaying chaos have been solved using a variety of direct numerical approaches. These techniques include the differential transform method [4], power series method [5], barycentric Lagrange interpolation collocation method [6], and differential quadrature method [7]. Unfortunately, the majority of direct numerical approaches now in use converge slowly for these issues, which results in erroneous approximations.

The applications that can be modeled using the excellent tool known as fractional analysis cannot be listed [8,9,10]. In the study of so-called biology, where the scaling power law of fractional order consistently manifests as an empirical explanation of such complex processes, fractional analysis has attracted growing attention [11,12,13]. The chaotic Lorenz system is just one of these mathematical equations. The so-called fractional Lorenz model was created by adding a modification to the Lorenz model. There are a good number of research papers that studied such biological models; see, for example, [14]. With the use of Appell-type Changhee polynomials, we construct an approximate formula for the fractional derivative for the first time in this work. We then use this formula to solve the fractional Lorenz model using the SCM. The SCM involves an approximation of the solution by summing up basis functions and determining their coefficients by enforcing the differential equation at a limited number of collocation points [15]. Among the variety of base functions available for use in the SCM, the Appell-type Changhee polynomials (ACPs) are a set of orthogonal polynomials that can serve as effective basis functions. Utilizing the SCM with ACPs comes with the benefit of their remarkable convergence properties. The accuracy of the solution is rapidly improved by increasing the number of collocation points. Moreover, ACPs exhibit good stability properties, making them suitable for solving differential equations that are stiff or have rapidly varying solutions. ACPs have the added advantage of having a closed-form expression, which simplifies their computation and manipulation. This feature can result in substantial computational savings, especially when compared to other techniques that necessitate solving extensive systems of ODEs [15].

The main objectives of the present study are addressed by the following points:

• We describe the proposed 9Dim chaotic Lorenz system using a theoretical and numerical simulation.
• Particular attention is given to using ACPs to provide a suitable formula for the Caputo fractional (CF) derivative.
• The recommended method and this approximation are used to convert the model into a set of algebraic equations. Through the use of the Newton iteration method, this system is numerically solved.
• With alternative values for the parameter r and varied values for the fractional order $\alpha$, we provide a numerical simulation of the model under consideration using the suggested methodology.
• The REF is then introduced to calculate the solution’s error. Additionally, we compared the solution produced by the suggested method with that produced by (RK4) and other previously published research.

2. Preliminaries

2.1. Definitions in Fractional Calculus

Although there are several definitions for fractional-order integration and differentiation, the Riemann–Liouville and Caputo fractional derivative definitions, which are given here, are the most important ones. The theory of fractional calculus was developed by using these definitions.

Definition 1.

The function $\psi \left(t\right)$’s fractional integral of order $\nu \in {\mathbb{R}}^{+}$ in the Riemann–Liouville sense is defined as follows [16]:

$$I^{\nu }\psi \left(t\right)=\frac{1}{\Gamma \left(\nu \right)}{\int }_0^t{(t-\tau )}^{\nu -1}\psi \left(\tau \right)d\tau ,t>0,$$

where $\Gamma (.)$ is the gamma function.

Its properties include:

$$I^{\nu }{I}^{\gamma }\psi \left(t\right)=I^{\nu +\gamma }\psi \left(t\right),\nu ,\gamma >0,$$

$$I^{\nu }{t}^m=\frac{\Gamma (m+1)}{\Gamma (m+\nu +1)}{t}^{m+\nu },$$

and the linearity property:

$$I^{\nu }[c_1{\psi }_1\left(t\right)+c_2{\psi }_2\left(t\right)]=c_1{I}^{\nu }{\psi }_1\left(t\right)+c_2{I}^{\nu }{\psi }_2\left(t\right),$$

for some constants $c_1$ and $c_2$.

Definition 2.

The function $\psi \left(t\right)$’s fractional derivative of order ν in the Riemann–Liouville sense is given as follows [16]:

$$D^{\nu }\psi \left(t\right)=\frac{d^m}{d{t}^m}\left(I^{m-\nu }\psi \left(t\right)\right),m-1<\nu \le m,m\in \mathbb{N}.$$

The Riemann–Liouville definition suffers from various limitations when simulating certain real-world problems [17]. However, Caputo’s definition, which we use, as outlined in the following definition, was created to deal with such concerns.

Definition 3.

The Caputo fractional derivative of order ν of a function $\psi \left(t\right)$ is given as follows [17]:

$$D^{\nu }\psi \left(t\right)=\frac{1}{\Gamma (n-\nu )}{\int }_0^t\frac{{\psi }^{\left(n\right)}\left(\tau \right)}{{(t-\tau )}^{\nu -n+1}}d\tau ,n-1<\nu <n,n\in \mathbb{N}$$

such that:

$$D^{\nu }C=0,whereCisaconstant,$$

$$D^{\nu }{t}^{\theta }=\frac{\Gamma (\theta +1)}{\Gamma (\theta +1-\nu )}{t}^{\theta -\nu },if\theta \in \mathbb{N}\cup \left\{0\right\},and\theta \ge \lceil \nu \rceil ,$$

(1)

where $\lceil \nu \rceil$ denotes the ceil function. Additionally, this operator $D^{\nu }$ is linear:

$$D^{\nu }[c_1{\psi }_1\left(t\right)+c_2{\psi }_2\left(t\right)]=c_1{D}^{\nu }{\psi }_1\left(t\right)+c_2{D}^{\nu }{\psi }_2\left(t\right),$$

for some constants $c_1$ and $c_2$.

2.2. Some Concepts of Changhee Polynomials

If $C{h}_m\left(t\right)$ denotes Changhee polynomials and Ch${}_m=C{h}_m\left(0\right)$ for Changhee numbers, as we know, $C{h}_m\left(t\right)$ and Ch${}_m$ are defined using the generating function [18,19]:

$$\frac{2}{z+2}{(1-z)}^t=\sum _{m=0}^{\infty }C{h}_m\left(t\right)\frac{z^m}{m!}.$$

where

$$C{h}_m\left(t\right)=\sum _{\ell =0}^{\infty }{S}_1(m,\ell )E_{\ell }\left(t\right),$$

such that $S_1(m,\ell )$ and $E_{\ell }\left(t\right)$ are first-kind Sterling numbers and Euler polynomials, respectively. Let ${\mathrm{Ch}}_m^{*}\left(t\right)$ be the Appell-type Changhee polynomials that are given by [20]:

$$\frac{2}{z+2}{e}^{tz}=\sum _{m=0}^{\infty }{\mathrm{Ch}}_m^{*}\left(t\right)\frac{z^m}{m!}.$$

The Appell-type Changhee polynomial of degree m is introduced by:

$${\mathrm{Ch}}_m^{*}\left(t\right)=\sum _{j=0}^m\left(\begin{array}{c}m\hfill \\ j\hfill \end{array}\right){\mathrm{Ch}}_{m-j}^{*}{t}^j.$$

(2)

From Formula (2), one can easily obtain:

$$\frac{d}{dt}{\mathrm{Ch}}_m^{*}\left(t\right)=\frac{d}{dt}m{\mathrm{Ch}}_{m-1}^{*}\left(t\right),$$

(3)

and from (3), we have:

$${\mathrm{Ch}}_m^{*}\left(t\right)={\int }_0^{t}m{\mathrm{Ch}}_{m-1}^{*}\left(y\right)dy+{\mathrm{Ch}}_m^{*}.$$

Clearly, ${\mathrm{Ch}}_0^{*}=1$ and $2{\mathrm{Ch}}_m^{*}+m{\mathrm{Ch}}_{m-1}^{*}=0,\forall m\ge 1$.

Additionally, it can proven that:

$${\int }_0^1{\mathrm{Ch}}_n^{*}\left(t\right){\mathrm{Ch}}_m^{*}\left(t\right)dt=\sum _{i=0}^m\sum _{k=0}^{m-i}(\begin{array}{c}m\\ i\end{array})\frac{{(-1)}^{m-i-1}(m-i)\left(\begin{array}{c}m-i\\ k\end{array}\right){\mathrm{Ch}}_k^{*}\left(1\right){\mathrm{Ch}}_i^{*}}{(2(m-i)-k+1)\left(\begin{array}{c}2(m-i)-k\\ m-i\end{array}\right)}.$$

(4)

Consider the following set of ACPs ${\left\{{\mathrm{Ch}}_i^{*}\left(t\right)\right\}}_{i=1}^m\subset {\mathrm{L}}^2[0,1]$, and let $\Omega =Span{\left\{{\mathrm{Ch}}_i^{*}\left(t\right)\right\}}_{i=1}^m$ be an ${\mathrm{L}}^2[0,1]$ subspace with finite dimensions [20]. The function $u\left(t\right)$ of ${\mathrm{L}}^2[0,1]$ provides a good unique approximation in $\Omega$. If $u^{*}\left(t\right)$ is the unique approximation of $u\left(t\right)$, we can use the following error estimation:

$${∥u\left(t\right)-u^{*}\left(t\right)∥}_2\le {\parallel u\left(t\right)-v\left(t\right)\parallel }_2,\forall v\left(t\right)\in \Omega .$$

Since $\Omega \subset L^2[0,1]$ and is closed, then according to [17], we can find that $L^2[0,1]=\Omega \oplus$${\Omega }^{\perp }$, where ${\Omega }^{\perp }$ denotes the orthogonal complement of $\Omega$, so we have $u\left(t\right)=v\left(t\right)+g\left(t\right)$ and thus $g\left(t\right)=u\left(t\right)-v\left(t\right)$, which also means that $u\left(t\right)-u^{*}\left(t\right)\in {\Omega }^{\perp }$. Then:

$$〈u\left(t\right)-u^{*}\left(t\right),v\left(t\right)〉=0,\forall v\left(t\right)\in \Omega ,$$

(5)

where $\langle .,.\rangle$ denotes the inner product.

Since $u^{*}\left(t\right)\in \Omega$, then it can be expressed as follows:

$$u\left(t\right)\approx u^{*}\left(t\right)=\sum _{i=1}^m{c}_i{\mathrm{Ch}}_i^{*}\left(t\right)=C^T{\mathrm{Ch}}^{*}\left(t\right),$$

(6)

where

$$\begin{array}{c}C={\left[c_1,c_2,\cdots ,c_m\right]}^T,{\mathrm{Ch}}^{*}\left(t\right)={\left[{\mathrm{Ch}}_1^{*}\left(t\right),{\mathrm{Ch}}_2^{*}\left(t\right),\dots ,{\mathrm{Ch}}_m^{*}\left(t\right)\right]}^T.\end{array}$$

Let $v\left(t\right)={\mathrm{Ch}}_i^{*}\left(t\right)$; from (6) in (5), we have:

$$〈u\left(t\right)-C^T{\mathrm{Ch}}^{*}\left(t\right),{\mathrm{Ch}}_i^{*}\left(t\right)〉=0.$$

(7)

Additionally, from (6), we can see that:

$$〈u\left(t\right),{\mathrm{Ch}}^{*}\left(t\right)〉=C^T〈{\mathrm{Ch}}^{*}\left(t\right),{\mathrm{Ch}}^{*}\left(t\right)〉=C^T\mathbb{A},$$

(8)

where $\mathbb{A}=〈{\mathrm{Ch}}^{*}\left(t\right),{\mathrm{Ch}}^{*}\left(t\right)〉$ is an $m\times m$ matrix and defined by:

$$\mathbb{A}={\int }_0^t{\mathrm{Ch}}^{*}\left(\tau \right){\mathrm{Ch}}^{*T}\left(\tau \right)d\tau ,$$

where $\mathbb{A}$ can be obtained by using (4). Then, the solution for (8) can be given as:

$$C={\mathbb{A}}^{-1}〈u\left(t\right),{\mathrm{Ch}}^{*}\left(t\right)〉.$$

3. Approximation of the Derivative D^ν via ACPs

In this part, we present an approximate formula for $D^{\nu }\left(u^{*}\left(t\right)\right)$ given in (6) through the following theorem.

Theorem 1.

The Caputo fractional derivative of order $\nu >0$ for $u^{*}\left(t\right)$ given in (6) can be estimated by:

$$D^{\nu }{u}^{*}\left(t\right)=\sum _{i=\lceil \nu \rceil }^m\sum _{j=\lceil \nu \rceil }^i{c}_i{\kappa }_{i,j,\nu }{t}^{j-\nu },$$

(9)

such that

$${\kappa }_{i,j,\nu }=\frac{(i!){\mathrm{Ch}}_{i-j}^{*}}{(i-j)!\Gamma (j+1-\nu )}.$$

where ${\mathrm{Ch}}_{i-j}^{*}$ is the Changhee number.

Proof. 

Let ${\mathrm{Ch}}_{\mathrm{i}}^{*}\left(\mathrm{t}\right)$ be an Appell-type Changhee polynomial of degree i, with $i=0,1,\dots ,m$. By using (1) and (2), we can obtain:

$$\begin{array}{cc}\hfill D^{\nu }{u}^{*}\left(t\right)& =\sum _{i=0}^m{c}_i{D}^{\nu }{\mathrm{Ch}}_i^{*}\left(t\right)=\sum _{i=\lceil \nu \rceil }^m\sum _{j=\lceil \nu \rceil }^i{c}_i\frac{(i!){\mathrm{Ch}}_{i-j}^{*}}{(j!)(i-j)!}{D}^{\nu }{t}^j\hfill \\ & =\sum _{i=\lceil \nu \rceil }^m\sum _{j=\lceil \nu \rceil }^i{c}_i\frac{(i!){\mathrm{Ch}}_{i-j}^{*}}{(i-j)!\Gamma (j+1-\nu )}{t}^{j-\nu }\hfill \\ & =\sum _{i=\lceil \nu \rceil }^m\sum _{j=\lceil \nu \rceil }^i{c}_i{\kappa }_{i,j,\nu }{t}^{j-\nu },\hfill \end{array}$$

where ${\kappa }_{i,j,\nu }$ is given in (9), and this completes the proof. □

4. Numerical Implementation

This section provides an overview of how the suggested approach for solving the 9D Lorenz system will be put into practice. By utilizing a method similar to the well-known 3D Lorenz system, the 3D spatial Boussinesq–Oberbeck equations that control thermal convection were expanded three times to obtain the 9D Lorenz system [6]. The system comes from:

$$D^{\alpha }{\varphi }_1\left(t\right)=-\sigma b_1{\varphi }_1-\sigma b_2{\varphi }_7-{\varphi }_2{\varphi }_4+b_4{\varphi }_4^2+b_3{\varphi }_3{\varphi }_5,$$

(10)

$$D^{\alpha }{\varphi }_2\left(t\right)=-\sigma {\varphi }_2-0.5\sigma {\varphi }_9+{\varphi }_1{\varphi }_4-{\varphi }_2{\varphi }_5+{\varphi }_4{\varphi }_5,$$

(11)

$$D^{\alpha }{\varphi }_3\left(t\right)=-\sigma b_1{\varphi }_3+\sigma b_2{\varphi }_8+{\varphi }_2{\varphi }_4-b_4{\varphi }_2^2-b_3{\varphi }_1{\varphi }_5,$$

(12)

$$D^{\alpha }{\varphi }_4\left(t\right)=-\sigma {\varphi }_4+0.5\sigma {\varphi }_9-{\varphi }_2{\varphi }_3-{\varphi }_2{\varphi }_5+{\varphi }_4{\varphi }_5,$$

(13)

$$D^{\alpha }{\varphi }_5\left(t\right)=-\sigma b_5{\varphi }_5+0.5{\varphi }_2^2-0.5{\varphi }_4^2,$$

(14)

$$D^{\alpha }{\varphi }_6\left(t\right)=-b_6{\varphi }_6+{\varphi }_2{\varphi }_9-{\varphi }_4{\varphi }_9,$$

(15)

$$D^{\alpha }{\varphi }_7\left(t\right)=-r{\varphi }_1-b_1{\varphi }_7+2{\varphi }_5{\varphi }_8-{\varphi }_4{\varphi }_9,$$

(16)

$$D^{\alpha }{\varphi }_8\left(t\right)=r{\varphi }_3-b_1{\varphi }_8-2{\varphi }_5{\varphi }_7+{\varphi }_2{\varphi }_9,$$

(17)

$$D^{\alpha }{\varphi }_9\left(t\right)=-r{\varphi }_2+r{\varphi }_4-{\varphi }_9-2{\varphi }_2{\varphi }_6+2{\varphi }_4{\varphi }_6+{\varphi }_4{\varphi }_7-{\varphi }_2,$$

(18)

where r is a parameter value (the reduced Rayleigh number) that classifies the system as chaotic or hyperchaotic, $\sigma$ is a constant, and the constant parameters $b_i, i=1,2,\cdots ,6,$ (measures of the geometry of the square cell) are defined by [21]:

$$b_1=4\frac{1+a^2}{1+2a^2},b_2=\frac{1+2a^2}{2\left(1+a^2\right)},b_3=2\frac{1-a^2}{1+a^2}{b}_4=\frac{a^2}{1+a^2},b_5=\frac{8a^2}{1+2a^2},b_6=\frac{4}{1+2a^2},$$

(19)

where a denotes the wavenumber in the horizontal direction. We consider the following initial conditions:

$${\varphi }_q\left(0\right)={\varphi }_q^0,q=1,2,\dots ,9.$$

(20)

Based on the memory effect of fractional derivatives, we may more precisely characterize the influence of the spread of the solution in the future and history with the model (10)–(20) in its fractional form. As is also known, fractional differential equations (FDEs) are necessary to convert several chaotic models in order to effectively explain the dynamics of chaotic systems. Although mathematical models with integer derivatives play an important role and have significance in understanding the dynamics of chaotic systems, they have some limitations, such as the fact that these systems lack memory or non-local effects. FDEs are commonly utilized in the theory of complex systems, the investigation of anomalous phenomena in nature, and in general, when considering the features of a curve across a wide area. Last but not least, this type of model explains temporal delays, fractal features, etc. See [6,21] for more information about this model.

Now, we implement the proposed method for solving the model (10)–(18) through the following steps:

Let ${\varphi }_{k,m}\left(t\right),k=1,2,\cdots ,9$, be the approximation of ${\varphi }_k\left(t\right)$ in terms of ACPs as follows:

$${\varphi }_{k,m}\left(t\right)=\sum _{i=0}^m{a}_{k,i}{\mathrm{Ch}}_i^{*}\left(t\right),k=1,2,\cdots ,9.$$

(21)

Using (9) and (21) in (10)–(18),

$$\begin{array}{cc}\hfill \sum _{i=\lceil \nu \rceil }^m& \sum _{j=\lceil \nu \rceil }^i{a}_{1,i}{\kappa }_{i,j,\nu }{t}^{j-\nu }=-\left(\sum _{i=0}^N(\sigma b_1{a}_{1,i}+\sigma b_2{a}_{7,i}){\mathrm{Ch}}_i^{*}\left(t\right)\right)-\left(\sum _{i=0}^N{a}_{2,i}{\mathrm{Ch}}_i^{*}\left(t\right)\right).\hfill \\ & \left(\sum _{i=0}^N{a}_{4,i}{\mathrm{Ch}}_i^{*}\left(t\right)\right)+b_4{\left(\sum _{i=0}^N{a}_{4,i}{\mathrm{Ch}}_i^{*}\left(t\right)\right)}^2+b_3\left(\sum _{i=0}^N{a}_{3,i}{\mathrm{Ch}}_i^{*}\left(t\right)\right)\left(\sum _{i=0}^N{a}_{4,i}{\mathrm{Ch}}_i^{*}\left(t\right)\right),\hfill \end{array}$$

(22)

$$\begin{array}{cc}\hfill \sum _{i=\lceil \nu \rceil }^m& \sum _{j=\lceil \nu \rceil }^i{a}_{2,i}{\kappa }_{i,j,\nu }{t}^{j-\nu }=-\sigma \left(\sum _{i=0}^N{a}_{2,i}{\mathrm{Ch}}_i^{*}\left(t\right)\right)-0.5\sigma \left(\sum _{i=0}^N{a}_{9,i}{\mathrm{Ch}}_i^{*}\left(t\right)\right)\hfill \\ & +\left(\sum _{i=0}^N{a}_{1,i}{\mathrm{Ch}}_i^{*}\left(t\right)\right)\left(\sum _{i=0}^N{a}_{4,i}{\mathrm{Ch}}_i^{*}\left(t\right)\right)-\left(\sum _{i=0}^N(a_{2,i}-a_{4,i}){\mathrm{Ch}}_i^{*}\left(t\right)\right)\left(\sum _{i=0}^N{a}_{5,i}{\mathrm{Ch}}_i^{*}\left(t\right)\right),\hfill \end{array}$$

(23)

$$\begin{array}{cc}\hfill \sum _{i=\lceil \nu \rceil }^m& \sum _{j=\lceil \nu \rceil }^i{a}_{3,i}{\kappa }_{i,j,\nu }{t}^{j-\nu }=-\left(\sum _{i=0}^N(\sigma b_1{a}_{3,i}-\sigma b_2{a}_{8,i}){\mathrm{Ch}}_i^{*}\left(t\right)\right)+\left(\sum _{i=0}^N{a}_{2,i}{\mathrm{Ch}}_i^{*}\left(t\right)\right).\hfill \\ & \left(\sum _{i=0}^N{a}_{4,i}{\mathrm{Ch}}_i^{*}\left(t\right)\right)-b_4{\left(\sum _{i=0}^N{a}_{2,i}{\mathrm{Ch}}_i^{*}\left(t\right)\right)}^2-b_3\left(\sum _{i=0}^N{a}_{1,i}{\mathrm{Ch}}_i^{*}\left(t\right)\right)\left(\sum _{i=0}^N{a}_{5,i}{\mathrm{Ch}}_i^{*}\left(t\right)\right),\hfill \end{array}$$

(24)

$$\begin{array}{cc}\hfill \sum _{i=\lceil \nu \rceil }^m& \sum _{j=\lceil \nu \rceil }^i{a}_{4,i}{\kappa }_{i,j,\nu }{t}^{j-\nu }=-\left(\sum _{i=0}^N(\sigma a_{4,i}-0.5\sigma a_{9,i}){\mathrm{Ch}}_i^{*}\left(t\right)\right)-\left(\sum _{i=0}^N(a_{3,i}-a_{5,i}){\mathrm{Ch}}_i^{*}\left(t\right)\right).\hfill \\ & \left(\sum _{i=0}^N{a}_{2,i}{\mathrm{Ch}}_i^{*}\left(t\right)\right)+\left(\sum _{i=0}^N{a}_{4,i}{\mathrm{Ch}}_i^{*}\left(t\right)\right)\left(\sum _{i=0}^N{a}_{5,i}{\mathrm{Ch}}_i^{*}\left(t\right)\right),\hfill \end{array}$$

(25)

$$\begin{array}{cc}\hfill \sum _{i=\lceil \nu \rceil }^m& \sum _{j=\lceil \nu \rceil }^i{a}_{5,i}{\kappa }_{i,j,\nu }{t}^{j-\nu }=-\sigma b_5\left(\sum _{i=0}^N{a}_{5,i}{\mathrm{Ch}}_i^{*}\left(t\right)\right)+0.5{\left(\sum _{i=0}^N{a}_{2,i}{\mathrm{Ch}}_i^{*}\left(t\right)\right)}^2-0.5{\left(\sum _{i=0}^N{a}_{4,i}{\mathrm{Ch}}_i^{*}\left(t\right)\right)}^2,\hfill \end{array}$$

(26)

$$\begin{array}{cc}\hfill \sum _{i=\lceil \nu \rceil }^m& \sum _{j=\lceil \nu \rceil }^i{a}_{6,i}{\kappa }_{i,j,\nu }{t}^{j-\nu }=-b_6\left(\sum _{i=0}^N{a}_{6,i}{\mathrm{Ch}}_i^{*}\left(t\right)\right)+\left(\sum _{i=0}^N(a_{2,i}-a_{4,i}){\mathrm{Ch}}_i^{*}\left(t\right)\right)\left(\sum _{i=0}^N{a}_{9,i}{\mathrm{Ch}}_i^{*}\left(t\right)\right),\hfill \end{array}$$

(27)

$$\begin{array}{cc}\hfill \sum _{i=\lceil \nu \rceil }^m& \sum _{j=\lceil \nu \rceil }^i{a}_{7,i}{\kappa }_{i,j,\nu }{t}^{j-\nu }=-r\left(\sum _{i=0}^N{a}_{1,i}{\mathrm{Ch}}_i^{*}\left(t\right)\right)-b_1\left(\sum _{i=0}^N{a}_{7,i}{\mathrm{Ch}}_i^{*}\left(t\right)\right)\hfill \\ & +2\left(\sum _{i=0}^N{a}_{5,i}{\mathrm{Ch}}_i^{*}\left(t\right)\right)\left(\sum _{i=0}^N{a}_{8,i}{\mathrm{Ch}}_i^{*}\left(t\right)\right)-\left(\sum _{i=0}^N{a}_{4,i}{\mathrm{Ch}}_i^{*}\left(t\right)\right)\left(\sum _{i=0}^N{a}_{9,i}{\mathrm{Ch}}_i^{*}\left(t\right)\right),\hfill \end{array}$$

(28)

$$\begin{array}{cc}\hfill \sum _{i=\lceil \nu \rceil }^m& \sum _{j=\lceil \nu \rceil }^i{a}_{8,i}{\kappa }_{i,j,\nu }{t}^{j-\nu }=r\left(\sum _{i=0}^N{a}_{3,i}{\mathrm{Ch}}_i^{*}\left(t\right)\right)-b_1\left(\sum _{i=0}^N{a}_{8,i}{\mathrm{Ch}}_i^{*}\left(t\right)\right)\hfill \\ & -2\left(\sum _{i=0}^N{a}_{5,i}{\mathrm{Ch}}_i^{*}\left(t\right)\right)\left(\sum _{i=0}^N{a}_{7,i}{\mathrm{Ch}}_i^{*}\left(t\right)\right)+\left(\sum _{i=0}^N{a}_{2,i}{\mathrm{Ch}}_i^{*}\left(t\right)\right)\left(\sum _{i=0}^N{a}_{9,i}{\mathrm{Ch}}_i^{*}\left(t\right)\right),\hfill \end{array}$$

(29)

$$\begin{array}{cc}\hfill \sum _{i=\lceil \nu \rceil }^m& \sum _{j=\lceil \nu \rceil }^i{a}_{9,i}{\kappa }_{i,j,\nu }{t}^{j-\nu }=\left(\sum _{i=0}^N(-r{a}_{1,i}-b_1{a}_{7,i}){\mathrm{Ch}}_i^{*}\left(t\right)\right)-2\left(\sum _{i=0}^N{a}_{2,i}{\mathrm{Ch}}_i^{*}\left(t\right)\right).\hfill \\ & \left(\sum _{i=0}^N{a}_{6,i}{\mathrm{Ch}}_i^{*}\left(t\right)\right)+\left(\sum _{i=0}^N(2a_{6,i}+a_{7,i}){\mathrm{Ch}}_i^{*}\left(t\right)\right)\left(\sum _{i=0}^N{a}_{4,i}{\mathrm{Ch}}_i^{*}\left(t\right)\right)-\left(\sum _{i=0}^N{a}_{2,i}{\mathrm{Ch}}_i^{*}\left(t\right)\right).\hfill \end{array}$$

(30)

By collocating (22)–(30) at m points $t_s=\frac{s}{m-1}+1, s=1,2,\cdots ,m$, it will be reduced to an algebraic system with the coefficients $a_{k,i},k=1,2,\cdots ,9,i=0,2,\cdots ,m.$

Substituting Equation (21) into (20) and using the fact that ${\mathrm{Ch}}_i^{*}\left(0\right)={\mathrm{Ch}}_i^{*}$, the initial conditions (20) will be converted to

$$\sum _{i=0}^m{a}_{k,i}{\mathrm{Ch}}_i^{*}={\varphi }_k^0,k=1,2,\cdots ,9.$$

(31)

We solve the nonlinear system of Equations (22)–(31) for the unknowns $a_{k,i},$$k=1,2,\dots ,9,i=0,1,\cdots ,m$ by using the Newton iteration approach. This prompts us to develop the approximate solution via substitution in the form of (21).

5. Numerical Simulation and Discussion

5.1. Graphical and Tabular Findings

We will verify the accuracy and quality of the obtained scheme. We address (10)–(20) with various values of $\alpha$, m, and the parameter value r. In all computations, we take $\sigma =0.25$ and $a=1/2$ for computing the quantities $b_i,i=1,2,\cdots ,6,$ from Equation (19), and:

$${\varphi }_{1,0}={\varphi }_{3,0}={\varphi }_{9,0}=0.01,{\varphi }_{2,0}={\varphi }_{4,0}={\varphi }_{5,0}={\varphi }_{6,0}={\varphi }_{7,0}={\varphi }_{8,0}=0.0.$$

Because Reiterer et al. [21] observed that when $r>43.3$, the system exhibits hyperchaotic behavior and remains chaotic otherwise, in our work, we present both cases by solving the model under study (10)–(20) for $r\in [14.1,15.1]$ and $r=55$ in Figure 1, Figure 2, Figure 3, Figure 4, Figure 5 and Figure 6. Numerical simulation work was carried out with the help of the Mathematica software (MATHEMATICS 11) package.

The approximate solution with $\alpha =1.0,0.9,0.8,0.7,$$m=6,$ and $r=14.1$ is shown in Figure 1a–i. Figure 2a–i presents the approximate solution with $r=14.1,15.0,55.0,$$\alpha =0.95,$ and $m=7$. In Figure 3a–i, we give a comparison between our numerical results and those results obtained by the RK4 method at $\alpha =1$ with $m=6$ and $r=14.1$. The REF [22] of the solution is shown in Figure 4a–i at $\alpha =0.96$ and $r=14.1$ with $m=6,10$ in $[0,2]$. Figure 5 and Figure 6 present the phase projections on ${\varphi }_6-{\varphi }_7$ and ${\varphi }_6-{\varphi }_9$ with distinct values of r, with $m=8$ and $\alpha =0.99$.

In addition, to strongly prove and confirm the effectiveness of the given method, we present a comparison with a previously published work solving the same model by using the SCM with the help of other orthogonal polynomials, named Gegenbauer wavelet polynomials [23]. This comparison is presented in

Table 1. Values of the REF for the present method (first row) and the method in [23] (second row) with $r=14.1$.

t	${\mathit{\varphi }}_1$	${\mathit{\varphi }}_2$	${\mathit{\varphi }}_3$	${\mathit{\varphi }}_4$	${\mathit{\varphi }}_5$	${\mathit{\varphi }}_6$	${\mathit{\varphi }}_7$	${\mathit{\varphi }}_8$	${\mathit{\varphi }}_9$
0.0	2.52$\times {10}^{-7}$	4.02$\times {10}^{-6}$	6.85$\times {10}^{-7}$	2.29$\times {10}^{-8}$	3.90$\times {10}^{-7}$	0.82$\times {10}^{-7}$	0.73$\times {10}^{-6}$	5.68$\times {10}^{-7}$	4.85$\times {10}^{-7}$
0.0	3.65$\times {10}^{-6}$	5.96$\times {10}^{-5}$	4.65$\times {10}^{-6}$	0.85$\times {10}^{-7}$	5.97$\times {10}^{-5}$	9.65$\times {10}^{-6}$	4.05$\times {10}^{-5}$	6.90$\times {10}^{-6}$	5.02$\times {10}^{-6}$
0.2	4.55$\times {10}^{-8}$	3.62$\times {10}^{-7}$	3.95$\times {10}^{-6}$	3.35$\times {10}^{-6}$	1.25$\times {10}^{-7}$	7.50$\times {10}^{-8}$	1.90$\times {10}^{-8}$	7.30$\times {10}^{-7}$	0.60$\times {10}^{-7}$
0.2	5.01$\times {10}^{-7}$	1.50$\times {10}^{-6}$	2.38$\times {10}^{-5}$	1.95$\times {10}^{-5}$	7.65$\times {10}^{-6}$	6.05$\times {10}^{-6}$	5.62$\times {10}^{-6}$	6.95$\times {10}^{-6}$	6.02$\times {10}^{-6}$
0.4	6.26$\times {10}^{-6}$	0.05$\times {10}^{-6}$	1.05$\times {10}^{-8}$	4.50$\times {10}^{-7}$	8.98$\times {10}^{-6}$	5.66$\times {10}^{-7}$	7.58$\times {10}^{-8}$	4.28$\times {10}^{-7}$	3.95$\times {10}^{-7}$
0.4	7.95$\times {10}^{-5}$	7.99$\times {10}^{-5}$	3.95$\times {10}^{-7}$	6.65$\times {10}^{-5}$	0.05$\times {10}^{-5}$	3.80$\times {10}^{-6}$	6.95$\times {10}^{-6}$	3.98$\times {10}^{-6}$	7.65$\times {10}^{-6}$
0.6	7.28$\times {10}^{-7}$	7.98$\times {10}^{-7}$	5.26$\times {10}^{-7}$	7.28$\times {10}^{-7}$	8.96$\times {10}^{-8}$	2.58$\times {10}^{-6}$	5.29$\times {10}^{-7}$	2.24$\times {10}^{-8}$	4.96$\times {10}^{-7}$
0.6	6.65$\times {10}^{-6}$	6.29$\times {10}^{-6}$	7.95$\times {10}^{-6}$	8.32$\times {10}^{-7}$	8.85$\times {10}^{-6}$	0.29$\times {10}^{-5}$	8.95$\times {10}^{-5}$	2.24$\times {10}^{-7}$	8.53$\times {10}^{-5}$
0.8	4.95$\times {10}^{-7}$	6.68$\times {10}^{-7}$	5.02$\times {10}^{-5}$	7.95$\times {10}^{-7}$	9.95$\times {10}^{-6}$	8.02$\times {10}^{-8}$	4.05$\times {10}^{-7}$	3.03$\times {10}^{-8}$	1.95$\times {10}^{-7}$
0.8	3.75$\times {10}^{-6}$	5.50$\times {10}^{-6}$	4.32$\times {10}^{-4}$	5.26$\times {10}^{-6}$	0.25$\times {10}^{-6}$	7.95$\times {10}^{-7}$	6.50$\times {10}^{-6}$	4.74$\times {10}^{-7}$	3.62$\times {10}^{-6}$
1.0	1.65$\times {10}^{-7}$	2.32$\times {10}^{-8}$	0.29$\times {10}^{-7}$	3.26$\times {10}^{-8}$	1.65$\times {10}^{-7}$	4.00$\times {10}^{-7}$	4.00$\times {10}^{-6}$	5.65$\times {10}^{-7}$	6.85$\times {10}^{-8}$
1.0	2.01$\times {10}^{-6}$	3.85$\times {10}^{-7}$	3.96$\times {10}^{-6}$	2.95$\times {10}^{-7}$	2.75$\times {10}^{-5}$	3.90$\times {10}^{-6}$	2.80$\times {10}^{-5}$	4.85$\times {10}^{-5}$	3.96$\times {10}^{-7}$

and

Table 2. Values of the REF for the present method (first row) and the method in [23] (second row) with $r=55$.

t	${\mathit{\varphi }}_1$	${\mathit{\varphi }}_2$	${\mathit{\varphi }}_3$	${\mathit{\varphi }}_4$	${\mathit{\varphi }}_5$	${\mathit{\varphi }}_6$	${\mathit{\varphi }}_7$	${\mathit{\varphi }}_8$	${\mathit{\varphi }}_9$
0.0	3.95$\times {10}^{-6}$	4.02$\times {10}^{-5}$	0.85$\times {10}^{-7}$	1.78$\times {10}^{-7}$	5.42$\times {10}^{-6}$	9.01$\times {10}^{-7}$	4.32$\times {10}^{-5}$	2.85$\times {10}^{-6}$	3.85$\times {10}^{-7}$
0.0	6.65$\times {10}^{-5}$	6.82$\times {10}^{-4}$	7.33$\times {10}^{-5}$	0.75$\times {10}^{-6}$	2.65$\times {10}^{-4}$	4.50$\times {10}^{-6}$	3.65$\times {10}^{-5}$	1.85$\times {10}^{-5}$	2.99$\times {10}^{-6}$
0.2	5.85$\times {10}^{-7}$	5.85$\times {10}^{-6}$	1.65$\times {10}^{-5}$	0.26$\times {10}^{-5}$	8.96$\times {10}^{-6}$	4.85$\times {10}^{-7}$	3.65$\times {10}^{-7}$	6.95$\times {10}^{-6}$	1.68$\times {10}^{-6}$
0.2	4.65$\times {10}^{-6}$	1.95$\times {10}^{-6}$	2.95$\times {10}^{-4}$	5.65$\times {10}^{-5}$	9.00$\times {10}^{-5}$	7.65$\times {10}^{-5}$	0.75$\times {10}^{-6}$	4.32$\times {10}^{-6}$	2.09$\times {10}^{-5}$
0.4	5.11$\times {10}^{-6}$	3.95$\times {10}^{-5}$	3.65$\times {10}^{-7}$	9.95$\times {10}^{-6}$	0.85$\times {10}^{-5}$	1.01$\times {10}^{-7}$	4.32$\times {10}^{-7}$	0.75$\times {10}^{-5}$	1.00$\times {10}^{-5}$
0.4	0.02$\times {10}^{-5}$	9.26$\times {10}^{-4}$	6.96$\times {10}^{-6}$	2.85$\times {10}^{-4}$	1.95$\times {10}^{-5}$	0.85$\times {10}^{-6}$	7.98$\times {10}^{-5}$	6.98$\times {10}^{-5}$	4.35$\times {10}^{-5}$
0.6	2.02$\times {10}^{-7}$	3.98$\times {10}^{-6}$	4.25$\times {10}^{-6}$	5.85$\times {10}^{-7}$	6.65$\times {10}^{-5}$	1.50$\times {10}^{-6}$	0.95$\times {10}^{-6}$	2.28$\times {10}^{-7}$	3.66$\times {10}^{-7}$
0.6	1.12$\times {10}^{-6}$	2.54$\times {10}^{-5}$	0.29$\times {10}^{-5}$	5.58$\times {10}^{-6}$	1.74$\times {10}^{-4}$	7.00$\times {10}^{-5}$	6.62$\times {10}^{-4}$	5.85$\times {10}^{-6}$	4.97$\times {10}^{-5}$
0.8	3.32$\times {10}^{-7}$	3.75$\times {10}^{-6}$	4.85$\times {10}^{-5}$	5.78$\times {10}^{-5}$	0.12$\times {10}^{-6}$	1.85$\times {10}^{-7}$	6.50$\times {10}^{-7}$	6.52$\times {10}^{-6}$	4.29$\times {10}^{-6}$
0.8	4.85$\times {10}^{-6}$	4.95$\times {10}^{-6}$	3.95$\times {10}^{-4}$	7.96$\times {10}^{-5}$	3.96$\times {10}^{-5}$	5.22$\times {10}^{-5}$	1.52$\times {10}^{-6}$	4.85$\times {10}^{-6}$	0.85$\times {10}^{-5}$
1.0	3.85$\times {10}^{-6}$	4.68$\times {10}^{-5}$	1.24$\times {10}^{-7}$	8.98$\times {10}^{-7}$	5.85$\times {10}^{-6}$	2.05$\times {10}^{-7}$	3.65$\times {10}^{-5}$	8.98$\times {10}^{-6}$	6.65$\times {10}^{-7}$
1.0	5.01$\times {10}^{-5}$	2.96$\times {10}^{-4}$	3.95$\times {10}^{-5}$	4.36$\times {10}^{-6}$	3.29$\times {10}^{-4}$	0.55$\times {10}^{-6}$	0.01$\times {10}^{-5}$	8.85$\times {10}^{-5}$	2.01$\times {10}^{-6}$

with different values of the parameter $r=14.1$ (for the chaotic case) and $r=55$ (for the hyperchaotic case), with $\alpha =0.95$ and $m=10$ in each method, but with the same initial conditions and the same parameters as in Figure 2. We computed the residual error function for each method to show that our presented method with the ACPs is more accurate and computationally effective in solving the given system.

5.2. Discussion and Recommendations

In this subsection, we present a detailed explanation of the scientific and applied significance of the numerical results that we obtained in the previous subsection, in addition to clarifying the extent to which these results are affected and their dependence on some parameters affecting them, such as $\alpha$, m, and r, along with listing some recommendations related to the proposed method or the model itself. To show the importance of this scientific paper and what it presents as a new work, we present these discussions through the following points and items:

• The proposed method is demonstrated in Figure 1 to be an effective way to solve the proposed model in its fractional version in the Caputo sense.
• The results in Figure 2 show that the proposed technique is a reliable method to simulate both chaotic and hyperchaotic behaviors and is consistent with the previous study by Kouagou et al. [24].
• In Figure 3, we can note that there is excellent agreement between our technique and the RK4 method in this special case with the integer derivative $\alpha =1$, and this indicates that the proposed method is suitable.
• Through the results in Figure 4, we can increase the speed of the numerical computation by controlling the order of approximation m to improve the results as well as the efficiency of the given scheme.
• The obtained phase portraits in Figure 5 and Figure 6 are consistent with those of Kouagou et al. [24]. This demonstrates that this technique can manage high-dimensional chaotic systems.

Finally, from the graphical and tabular findings, we can see that the method used to solve the model under study in its fractional version using the Caputo operator is computationally accurate and suitable for solving complex dynamical systems with chaotic or hyperchaotic behavior. Additionally, we confirm that the numerical solution obtained from the given method depends on the values of $\alpha ,r,$ and m, as expected. Additionally, the suggested strategy significantly boosts the method’s effectiveness and results.

6. Conclusions

As we know very well, most of the numerical techniques now in use converge slowly for systems characterized by chaotic solutions, and this leads to inaccurate approximations. Though many of those methods can be highly accurate, they are less accurate for these kinds of models. As a result, we recommend the method described in this article. The results and the method’s efficiency are both markedly enhanced by this process.

In this study, a numerical simulation was performed by finding the numerical solutions for the nine-dimensional fractional Lorenz model for different values of $\alpha$, m, and the parameter r. To achieve this aim, we applied the SCM with ACPs. We first derive an approximate formula for the fractional derivative and use it to obtain the numerical scheme for solving the proposed model. This procedure creates a system of algebraic equations, and this system is solved by the Newton iteration method. To evaluate the method and measure its efficiency and accuracy, the residual error function was calculated, and a comparison was made with the RK4 method. From this comparison, we can confirm that the presented scheme is very suitable for the effective numerical study of this model. We can also control and reduce the error by increasing m, which increases the number of terms in the series. Additionally, we conclude that the Caputo operator used here is more suitable for describing the proposed model through the presented numerical simulations. Finally, the graphical and tabular findings show that the proposed technique is computationally accurate. We intend to deal with this model in the future, but on a larger scale, by generalizing this research to include additional types of polynomials or fractional derivatives.