Identification of numerical solutions of a fractal-fractional divorce epidemic model of nonlinear systems via anti-divorce counseling

1.   Introduction

Chaos is very interesting nonlinear phenomenon and has applications in many areas such as biology, economics, signal generator design, secure communication, many other engineering systems and so on. Because a nonlinear system in the chaotic state is very sensitive to its initial condition and chaos causes often irregular behavior in practical systems, chaos is sometimes undesirable ^{[1,2,3,4,5,6,7]}. The 3-dimensional Lorenz model (3DLM) was the first chaotic system proposed in the literature ^[3]. Celikovsky and Chen ^[2,8,9] introduced generalized Lorenz canonical forms, covering a large class of 3-dimensional autonomous systems. Further, Zhang et al. ^[10] by combining the advantages of both integer-order and fractional-order complex chaotic systems, proposed a hybrid-order complex Lorenz system. Chaotic systems have been widely addressed both from their mathematical properties ^[11,12] and practical applications ^{[11,12,13,14,15]}.

The 3DLM reveals the dependence of the solutions on the initial conditions for chaotic situations. Higher order Lorenz models have been derived to further study the stability of solutions and paths to chaos. The 3DLM, was obtained from the Rayleigh-Benard convection equations, which combine dissipative, heating and nonlinear advection physical processes. An interesting topic is the investigation of modes in the 3DLM and its generalization to higher dimensional models by increasing the number of modes. Shen ^[16,17] generalized the 3DLM to the 5DLM by adding two additional Fourier modes. This led to a better understanding of the role of some variables in the stability of solutions. Therefore, the 5DLM allows the role of modes in the predictability of solutions to be investigated for further understanding of the variables that increase the stability of solutions, and to find analytical solutions for critical points. Although the role of the modes on the stability of solutions, as well as some dynamical properties of the 5DLM, have been studied, the calculation of bounds for the variables and the global dynamics of the 5DLM are still open and challenging problems. Shen also extended the 5DLM to a six-dimensional Lorenz model (6DLM) ^[18], seven-dimensional Lorenz model (7DLM) ^[19] and generalized Lorenz model (GLM) ^[20], in order to examine the impact of an additional mode and its accompanying heating term on solution stability.

Given a chaotic dynamical system, if its chaotic attractor is bounded in the phase space and the trajectories of the system remain in a bounded region of the phase space, then we say that "the chaotic dynamical system is bounded". Estimation of the ultimate bound set (UBS) and the positive invariant set (PIS) of a chaotic system plays a very important role in studying its dynamic behavior. Among the most important applications, one can mention their use in controlling and synchronizing chaotic systems ^{[21,22,23,24,25,26,27,28,29]}. In fact, the bounds are necessary for both theoretical studies of chaotic attractors and numerical search of attractors. If we can show that, under certain considerations, there exists a GEAS for a chaotic system, then we can conclude that the system cannot have periodic or quasi-periodic responses, equilibrium points, or hidden attractors, outside this set of attractors. This issue has a great application in controlling systems and preventing their possible problems. Leonov ^[21] derived the first results about global UBS for the Lorenz model. Subsequently, Swinnerton-Dyer ^[30] demonstrated that the bounds of the states of the Lorenz equations could be determined by using Lyapunov functions. Several researchers further developed the idea and computed the GEAS and PIS for different chaotic systems ^{[31,32,33,34]}.

To the best of our knowledge, the GEAS and UBS for the 5DLM have not been investigated yet. In the present work, by changing system parameters and conditions, we create different attractive sets. Also, we calculate a small attractive set only dependent on the system parameters. The results obtained from the UBS have been used in the synchronization and control of dynamical systems ^[35,36,37]. Due to the importance of minimizing the synchronization time, by applying a finite time control scheme, an efficient synchronization method is given based on the obtained ultimate bound ^[38,39]. By developing these techniques, we can also estimate the ultimate bound of fractional chaotic systems ^[40,41,42].

This article is organized into 6 sections. The dynamical behavior of the 5DLM, including phase portraits, bifurcation diagrams and Hamilton energy are given in Section 2. In Section 3, we introduce some preliminary definitions and GEAS of the system. In Section 4, we present a method to compute small bound for the 5DLM. Section 5 presents the finite time synchronization problem using the results obtained in Section 4. The main conclusions are given in Section 6.

2.   Five-dimensional Lorenz model

The Rayleigh-Benard model for 2-dimensional $(x, z),$ dissipative and forced convection is ^[3]:

According to the studies of Rayleigh ^[3] and Saltzman ^[4], the following equations were obtained

where, $x_1, x_2$ and $x_3$ are function of time alone. All the parameters and variables mentioned above are given in Table 1.

By making some changes and manipulations, the partial differential equations (2.1) and (2.2) are converted into ordinary differential equations. Thus, the 3DLM chaotic system is expressed as ^[3]:

where, $\sigma, r$ are the Prandtl number, normalized Rayleigh number or the heating parameter and $b = \frac{4}{1+a^2}.$ Shen et al. ^[16] extended the 3DLM to the five-dimensional LM (5DLM) by including two additional Fourier modes with two additional vertical wave numbers. They used the five Fourier modes and rewrote $\psi$ and $\theta$ as the following:

where,

An additional mode $M_4 = \sqrt{2}\sin(lx) \sin(3mz)$ is included to derive the 6DLM. Here, $l$ and $m$ are defined as $\frac{\pi a}{H}$ and $\frac{\pi}{H},$ representing the horizontal and vertical wave numbers, respectively, and $H$ is the domain height and $\frac{2H}{a}$ represents the domain width.

By coordinate transformation, the original equation can be reduced to the following five-dimensional nonlinear dynamics:

where $d = \frac{9+a^2}{1+a^2}.$

Numerical analysis shows that the dynamical behavior of (2.8) changes from steady-state to chaotic, with the increase of $r$.

Figure 1 depicts the bifurcation diagram when the parameters values $\sigma = 10, b = \frac{8}{3}$ and $d = \frac{19}{3}$ are fixed, and $r$ varies on the interval $[0,100].$ In fact, one can see that chaos occurs after $r > 42.5$. For the value of the parameters $\sigma = 10, b = \frac{8}{3}, d = \frac{19}{3}, r = 43$ and $r = 25$, Lyapunov exponents are shown in Figure 2. The values of Lyapunov exponents at $500$th second are $L_1 = 1.1281$, $L_2 = 0.0073$, $L_3 = -1.3786$, $L_4 = -1.3774$ and $L_5 = -1.6688$. It is easy to observe that if $r = 43,$ then system (2.8) has the positive largest Lyapunov exponent. Therefore, the system (2.8) can exhibit chaotic behaviors. When selecting parameters $\sigma = 10, b = \frac{8}{3}, r = 25$ and $d = \frac{19}{3}$, the values of Lyapunov exponents at $500$th second are $L_1 = -0.5246$, $L_2 = -0.5358$, $L_3 = -6.9952$, $L_4 = -2.0305$ and $L_5 = -4.1633$. All negative Lyapunov exponents indicate that the behavior of the system is non-chaotic.

The time responses of the system for $r = 25$ and $r = 43$ are depicted in Figure 3. We verify that the system with $r = 25$ produces a steady-state solution and when $r = 43$ the system is in a chaotic state. Figure 4 depicts the phase portraits of system (2.8) when $\sigma = 10, b = \frac{8}{3}, r = 25$ and $d = \frac{19}{3}$. Figure 5 shows its chaotic behavior with $\sigma = 10, b = \frac{8}{3}, r = 43$ and $d = \frac{19}{3}$.

In this section, the Hamilton energy for the 5-dimensional Lorenz model (5DLM) is investigated. The Hamilton energy plays a crucial role in the stability of dynamical systems ^[43]. By continuously pumping or releasing energy in the system, we are able to stabilize chaos. Furthermore, the relation between the Hamilton energy and different chaotic attractors of system (2.8) and the energy dependence on attractors are discussed. Calculation of Hamilton energy for high-order Lorenz systems including six-dimensional Lorenz model (6DLM) ^[18], seven-dimensional Lorenz model (7DLM) ^[19] and generalized Lorenz model (GLM) ^[20], can be an interesting topic due to their special physical nature.

Let us rewrite (2.8) in the form:

where $\nabla H$ stands for the gradient vector of a smooth energy function $H(X)$, $J(X)$ represents a skew-symmetric matrix, and $R(X)$ denotes a symmetric matrix.

The Hamilton energy function can be expressed as:

Using the Helmholtz's theorem ^[43], it can be represented by:

Thus, we have:

Furthermore, it rate of variation is:

Figure 6 illustrates that chaotic and steady-state require lower and higher Hamilton energy, respectively.

3.   Ultimate bound sets

In this section, we pursue the goal of proving the existence of GEAS for the chaotic system (2.8). We will first mention a few prerequisites and definitions.

Let $\mathfrak{X} = (x_1, x_2, x_3, x_4, x_5)^T$, and consider that $\mathfrak{X}(t, t_0, \mathfrak{X}_0)$ is the solution of system

that satisfies $\mathfrak{X}_0 = \mathfrak{X}(t_0, t_0, \mathfrak{X}_0)$, with $t_0\geq 0$ representing the initial time. Also, $g:\mathbb{R}^5\rightarrow \mathbb{R}^5$ and $\Psi\subset \mathbb{R}^5$ is a compact set. Let us define the distance between $\mathfrak{X}(t, t_0, \mathfrak{X}_0)$ and $\Psi$ as:

Denote $\Psi_{\gamma} = \{\mathfrak{X}\mid\rho(\mathfrak{X}, \Psi) < \gamma\}$. Thus, one gets $\Psi \subset \Psi_{\gamma}$.

Definition 3.1. ^[21]. Suppose that $\Psi \subset \mathbb{R}^5$ is a compact set. If for any $\mathfrak{X}_0 \in \mathbb{R}^5/ \Psi,$

then $\Psi$ is an UBS of (3.1). Moreover, if for any $\mathfrak{X}_0 \in \Psi$ and all $t\geq t_0$, $\mathfrak{X}(t, t_0, \mathfrak{X}_0) \in \Psi$, then $\Psi$ is the PIS for (3.1).

From the above, it is interpreted that having a GEAS for a system guarantees that the system is UBS.

Definition 3.2. ^[21]. Given a Lyapunov function $\mathscr{L}_{\gamma}(\mathfrak{X})$, if there exist constants $K_{\gamma} > 0$ and $s_{\gamma} > 0,$ such that

for $\mathscr{L}_{\gamma}(\mathfrak{X}_0) > K_{\gamma}$ and $\mathscr{L}_{\gamma}(\mathfrak{X}) > K_{\gamma}$, then $\Psi_{\gamma} = \{\mathfrak{X}|\mathscr{L}_{\gamma}(\mathfrak{X}(t)) \leq K_{\gamma}\}$ is a GEAS of (2.8). Moreover, if for any $\mathfrak{X}_0 \in \Psi_{\gamma}$ and all $t > t_0$, $\mathfrak{X}(t, t_0, \mathfrak{X}_0) \in \Psi_{\gamma}$, then $\Psi_{\gamma}$ is a PIS.

The next theorem introduces the GEAS for system (2.8).

Theorem 3.1. For any $\sigma > 0,$ $b > 0,$ $r > 0, d > 0,$ $\beta > 0$ and $\alpha > 0$, with

If $\mathscr{L}_{\alpha, \beta}(\mathfrak{X}(t))\geq K_{\alpha, \beta}, \; t\geq t_0,$ then

This indicates that

is a GEAS and PIS of the system.

Proof. Let us define

and consider the definite positive Lyapunov function

The derivative of $\mathscr{L}_{\alpha, \beta}$ is as follows

From (3.5) we have

therefore,

when $\mathscr{L}_{\alpha, \beta}(\mathfrak{X}(t))\geq K_{\alpha, \beta}$. Then, we have the equivalent

Thus, if $\mathscr{L}_{\alpha, \beta}(\mathfrak{X}(t))\geq K_{\alpha, \beta}, \; t\geq t_0$, the following inequality results

By calculating the limit, one has

Therefore, the ellipsoid

is the GEAS and PIS of system (2.8). This ends the proof. □

Different cases can be highlighted:

Case 1: For $\alpha = 1, \beta = 1,$ then

is the GEAS of (2.8). Fore $\sigma = 10, b = \frac{8}{3}, r = 25$ and $d = \frac{19}{3},$ it yields

Figure 7 illustrates the attractors of (2.8) in distinct spaces by $\Psi_{1, 1}.$ For $\sigma = 10, b = \frac{8}{3}, r = 43$ and $d = \frac{19}{3},$ results

Figure 8 shows the attractors in different spaces defined by $\Psi_{1, 1}.$

Case 2: Let us consider $\alpha = 1, \beta = 2.$ Thus, the set

is the GEAS of (2.8).

For $\sigma = 10, b = \frac{8}{3}, r = 25$ and $d = \frac{19}{3},$ we have

Case 3: Define $\alpha = 2, \beta = 1$. Then,

is the GEAS of (2.8).

For $\sigma = 10, b = \frac{8}{3}, r = 25$ and $d = \frac{19}{3},$ one has

4.   Accurate ultimate bound

In this section we derive a more accurate and smaller boundary set than that established by Theorem 3.1. We state the following theorem.

Theorem 4.1. If $\sigma > 0, b > 0, r > 0$ and $d > 0,$ then we have the following boundaries for system (2.8) variables:

where

Proof. According to the results obtained from Theorem 3.1, we have

It is clear from the above equations that the upper bound of $x_1, x_2, \dots, x_5,$ depends on the lower bound of $\frac{\eta}{\sqrt{\alpha}}$ and $\frac{\eta}{\sqrt{\beta}}$ :

Therefore, one can determine the bound for $x_1$ as shown in Eq (4.1). For variables $x_2, x_3, x_4$ and $x_5,$ there is the same term:

Let us take $\frac{\alpha}{\beta} = \frac{1}{N}$ with $N\in \mathbb{N}$. Therefore,

Since

we find a more limited boundary set for $x_2$ as shown in Eq (4.2), and by doing the same process one can obtain

To confirm the theoretical results, we fix $\sigma = 10, b = \frac{8}{3}, r = 25$ and $d = \frac{19}{3}.$ Figure 9, shows the estimated bounds for each state variable. Furthermore, Table 2 compares the bounds estimated by Theorems 3.1 and 4.1, showing the advantage of Theorem 4.1.

□

Remark 4.1. A noteworthy point in calculating the bound for system (2.8) based on Theorem 4.1 is that the estimated boundary, in addition to being smaller, is independent of parameters $\beta$ and $\alpha.$

5.   Finite-time synchronization of the 5DLM

This section addresses the finite time synchronization of the 5DLM. In order to achieve fast and reliable synchronization, we use the results obtained in the previous section for the ultimate bound of (2.8). Let us consider the 5DLM (2.8) as the master, and define the slave by:

where $y_1, y_2, y_3, y_4, y_5$ are state vectors. The control law $u_1, u_2, u_3, u_4, u_5$ is designed for the drive (2.8) and response (5.1) systems can reach synchronization in finite time.

Lemma 5.1. Suppose that $a_1, a_2, \cdots, a_n$ are real numbers and that $0 < s < 1$. Then, we have:

Lemma 5.2. Inequality $2ab \leq \epsilon a^2 +\frac{1}{\epsilon} b^2$ holds for all real numbers $\epsilon > 0, a > 0$ and $b > 0$.

Lemma 5.3. Let us assume that $\mathcal{L}(t)$ is a continuous and positive-definite function that satisfies:

Then, the system is exponentially finite-time stable, where $\lambda, \mu > 0,$ and $0 < \omega < 1$ are constants.

Lemma 5.4. Let us suppose that $\mathcal{L}(t)$ is a Lyapunov function satisfying Eq. (5.2). Then, it holds

and $\mathcal{L}(t) = 0, \; \; \forall t\geq T$. The finite time $T$ is

Definition 5.1. Let us consider $\delta(t) = [\delta_1(t)\; \delta_2(t)\; \delta_3(t)\; \delta_4(t)\; \delta_5(t)]^T$, and define $\delta_i(t) = y_i(t) -x_i(t)$ as the synchronization errors. If there exists a positive value $T^*$ such that

and $\parallel \delta(t) \parallel = 0,$ for $t\geq T^*,$ then the systems (2.8) and (5.1), achieve finite-time synchronization.

The next theorem states the exponential finite-time synchronization condition of systems (2.8) and (5.1).

Theorem 5.1. The systems (2.8) and (5.1) can achieve finite-time synchronization by the control law:

where $\lambda, \mu > 0,$ and $s\in(0, 1)$. Further, the systems are synchronized in the time

where $\omega = \frac{s+1}{2}$.

Proof. Theorem 4.1 provides an accurate estimate of the ultimate bound of variables of system (2.8). Let consider $\alpha = \beta = 1$, then $\eta = \sigma+r$. Therefore, we have

From (2.8) and (5.1), the following error dynamics results:

Let us consider the Lyapunov function:

The time domain derivative of $\mathcal{L}(t)$ along the trajectories of (5.8) and using control law (5.5) is given by

From Lemma 5.2, we have

Then using Lemma 5.3, we obtain

where $\lambda, \mu > 0, \omega = \frac{1+s}{2}$. This implies that

Therefore,

which leads to

and since $\parallel \delta(t) \parallel^2 = \mathcal{L}_1(t),$ we obtain

This, from Lemma 5.3 and Lemma 5.4, guarantees the exponential synchronization of systems (2.8) and (5.1) in finite time $T^*$.

To be more realistic, an approach is presented here that requires only one controller to implement synchronization. □

Corollary 5.1. When the control functions are chosen as

then, the drive system (2.8) is exponentially synchronized with the slave system (5.1).

Proof. The Lyapunov function is selected as

In view of (5.8) and (5.10), the derivative of $\mathcal{L}(t)$ is

According to Lemma 5.2, Lemma 5.3 and (5.9) we have

where $\lambda, \mu > 0, \omega = \frac{1+s}{2}$. Therefore, synchronization is achieved exponentially with the control law (5.10).

To show the effectiveness of the control proposed in theorem 5.1, we perform numerical simulations. Let us choose the initial conditions $x_1(0) = 1, \; x_2(0) = -1, \; x_3(0) = 2, \; x_4(0) = 1, \; x_5(0) = 1,$ $y_1(0) = 0.1, \; y_2(0) = 0.1, \; y_3(0) = 0.2, \; y_4(0) = -1, \; y_5(0) = 1,$ and other parameters as $\lambda = 1,$ $\mu = 1$ and $s = \frac{1}{3}$. Figure 10 depicts the time series of the drive system (2.8) and the response system (5.1) without input control, in which the goal of synchronization has not been achieved. Figure 11 shows the case of using the proposed control (5.5). We verify that system (5.1) exponentially synchronizes with the master system (2.8) within the guaranteed convergence time.

The synchronization error for different modes with different controllers are depicted in Figure 12. The noteworthy point in these figures is that the goal of synchronization is achieved faster by increasing the number of controllers.

□

6.   Conclusions

The global dynamics of the 5DLM, which was obtained by increasing two modes to the original Lorenz system, was analyzed. Phase portraits, bifurcation diagrams and GEAS were estimated. Due to the dependence of the GEAS on the free parameters, a new boundary for the variables was estimated, which is more accurate than the existing one. Also, by employing a finite time control scheme, a synchronization method was proposed based on the obtained ultimate bound sets. The corresponding boundedness was numerically verified to demonstrate the efficiency of the presented method.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Conflict of interest

The authors declare that there are no conflicts of interests regarding the publication of this article.