Synchronization of Fractional Partial Difference Equations via Linear Methods

Abstract

Discrete fractional models with reaction-diffusion have gained significance in the scientific field in recent years, not only due to the need for numerical simulation but also due to the stated biological processes. In this paper, we investigate the problem of synchronization-control in a fractional discrete nonlinear bacterial culture reaction-diffusion model using the Caputo h-difference operator and a second-order central difference scheme and an L1 finite difference scheme after deriving the discrete fractional version of the well-known Degn–Harrison system and Lengyel–Epstein system. Using appropriate techniques and the direct Lyapunov method, the conditions for full synchronization are determined.Furthermore, this research shows that the L1 finite difference scheme and the second-order central difference scheme may successfully retain the properties of the related continuous system. The conclusions are proven throughout the paper using two major biological models, and numerical simulations are carried out to demonstrate the practical use of the recommended technique.

1. Introduction

One of the most essential components of dynamic system analysis is the construction of adequate functions identified as controllers to ensure synchronization. To understand how these systems achieve their distinctive synchronization behavior, a great variety of mathematical models have been suggested and studied. For continuous-space systems, mathematical modeling of oscillating biological or chemical media, for example, takes the form of reaction-diffusion equations. This type of model shows intricate dynamical structures such as bifurcations, spatial patterns, and turning instability. It has been demonstrated that reaction-diffusion systems, such as low-dimensional oscillators, may exhibit synchronization. For example, Mesdoui et al. [1] examined the synchronization of the Degn–Harrison reaction-diffusion system. Ref. [2] was concerned with the synchronization control of the Lengyel–Epstein reaction-diffusion system. Furthermore, the synchronization of the FitzHugh–Nagumo reaction-diffusion model using a particular control rule was detailed in [3]. Other works regarding this subject may be found in [4,5].

In real-world applications, fractional-order nonlinear equations are frequently employed to describe a wide range of physical phenomena [6,7,8,9]. Scientists are still fascinated by fractional calculus because of its numerous applications in physics, chemistry, biology, electronics, electrical engineering, mechanics, signal processing, and control [10,11,12,13,14,15]. As a result, in recent years, scholars have grown particularly interested in it. However, over the last decade, there has been a spike in attention to fractional reaction-diffusion systems, particularly on the topic of synchronization. For instance, in [16], a hybrid technique for synchronizing between two integer and fractional-order reaction-diffusion systems is proposed, with applications to particular chemical models. Moreover, In [17], the dynamics of the activator–inhibitor system known as the Gierer–Meinhardt system, which is utilized to describe the interactions of chemical and biological phenomena, was investigated.

The discrete form of fractional calculus is a novel approach with enormous potential applications in a variety of scientific and industrial fields. The application has attracted tremendous attention in the past few years (see [18,19,20,21,22,23,24,25]). The purpose of the latest investigations in this area is fundamental. Fractional difference equations, on the one hand, enhance classical differential equations. In addition, they provide for a feasible comparison of the behaviors of fractional difference and fractional differential equations.

Many physical phenomena rely on spatially discrete systems, often known as discrete reaction-diffusion systems. In fact, the discreteness and structure of the underlying spatial domain influence dynamical behavior significantly. Active PIN-induced transport across cell membranes, for instance, is required for auxin spreading across plant leaves [26]. Peierls–Nabarro barriers are often used to prevent tiny faults from propagating across discrete media, although they can be avoided by carefully modifying system characteristics [27]. Discrete reaction-diffusion systems are more similar to biological systems than continuous ones, and certain investigations on the behavior of such systems are particularly fascinating (see [25,28,29,30]). Nevertheless, fractional discrete reaction-diffusion equations have not yet been extensively studied [31]. A fractional discrete diffusion equation was presented by [32]. In [33], the chaotic behavior of a variable fractional diffusion equation on discontinuous time scales is examined. Clearly, there is a gap in our comprehension of the dynamics of such systems.

To the best of the authors’ knowledge, this is the first work dealing with the synchronization and control of discrete fractional reaction-diffusion systems. This has prompted us to investigate the issue of complete synchronization in coupled discrete fractional reaction-diffusion systems. With the help of the fractional Lyapunov approach, linear control laws for the discrete fractional reaction-diffusion Degn–Harrison system and Lengyel–Epstein systems have been proposed after driving the discrete version of the considered systems using the L1 finite difference scheme and the second-order central difference scheme. The following is how the paper is managed: Section 2 introduces some essential concepts and lemmas for discrete fractional calculus. Section 3 describes the models investigated in this study, which are the fractional discrete reaction-diffusion Lengyel–Epstein and Degn–Harrison systems, and presents a unique discrete temporal fractional reaction-diffusion system. Section 4 contains the discrete fractional Lengyel–Epstein reaction-diffusion system’s master-slave formulation, along with unique control rules and demonstrations of convergence based on an appropriate Lyapunov functional. Section 5 employs the same approach to drive the master-slave discrete fractional Degn–Harrison reaction-diffusion system, as well as control laws and proofs of convergence. In Section 6, control laws are derived analytically and numerically in two dimensions to achieve synchronization between the master-slave systems of the investigated models.

2. Preliminaries

This part starts with an overview of some of the topic’s primary concepts.

Definition 1

([21]). Assuming $\mathrm{x}:\mathbb{N}\to \mathbb{R},$ the forward difference operator Delta is expressed by

$$\mathrm{\Delta }\mathrm{x}\left(\mathfrak{i}\right)=\mathrm{x}(\mathfrak{i}+1)-\mathrm{x}\left(\mathfrak{i}\right),\mathfrak{i}\in \mathbb{N}.$$

Additionally, the operators ${\mathrm{\Delta }}^n; n=1; 2; 3; \dots ,$ are recursively determined by

$${\mathrm{\Delta }}^n\mathrm{x}\left(\mathfrak{i}\right)=\mathrm{\Delta }\left({\mathrm{\Delta }}^{n-1}\mathrm{x}\left(\mathfrak{i}\right)\right),\mathfrak{i}\in \mathbb{N}.$$

More specifically, the second-order difference operator of the function $\mathrm{x}\left(\mathfrak{i}\right)$ is provided by

$${\mathrm{\Delta }}^2\mathrm{x}\left(\mathfrak{i}\right)=\mathrm{x}(\mathfrak{i}+2)-2\mathrm{x}(\mathfrak{i}+1)+\mathrm{x}\left(\mathfrak{i}\right).$$

(1)

Theorem 1

([21]). Given two functions $\mathrm{x};\mathrm{y}:\begin{array}{c}\hfill {\mathbb{N}}_{\mathrm{a}}\end{array}\to \mathbb{R}$ and $\mathrm{a};\mathrm{b}\in \mathbb{N};\mathrm{a}<\begin{array}{c}\hfill \mathrm{b}\end{array};$ we have the following formulae for summation by parts:

$$\begin{array}{cc}\hfill & \sum _{\mathfrak{i}=\mathrm{a}}^{\mathrm{b}-1}\mathrm{x}\left(\mathfrak{i}\right)\mathrm{\Delta }\mathrm{y}\left(\mathfrak{i}\right)=\mathrm{x}\left(\mathfrak{i}\right)\mathrm{y}\left(\mathfrak{i}\right){|}_{\mathrm{a}}^{\mathrm{b}}-\sum _{\mathfrak{i}=\mathrm{a}}^{\mathrm{b}-1}\mathrm{y}(\mathfrak{i}+1)\mathrm{\Delta }\mathrm{x}\left(\mathfrak{i}\right),\hfill \\ \hfill & \sum _{\mathfrak{i}=\mathrm{a}}^{\mathrm{b}-1}\mathrm{x}(\mathfrak{i}+1)\mathrm{\Delta }\mathrm{y}\left(\mathfrak{i}\right)=\mathrm{x}\left(\mathfrak{i}\right)\mathrm{y}\left(\mathfrak{i}\right){|}_{\mathrm{a}}^{\mathrm{b}}-\sum _{\mathfrak{i}=\mathrm{a}}^{\mathrm{b}-1}\mathrm{y}\left(\mathfrak{i}\right)\mathrm{\Delta }\mathrm{x}\left(\mathfrak{i}\right).\hfill \end{array}$$

Definition 2

([22]). Let $\begin{array}{c}\hfill \mathrm{x}\end{array}\in {\left(\mathfrak{h}\mathbb{N}\right)}_{\mathrm{a}}\to \mathbb{R}$. The $\mathfrak{h}$-sum of the $\begin{array}{c}\hfill \zeta -\end{array}$th order for each $\zeta >0$ has been provided by

$${}_{\mathfrak{h}}{\mathrm{\Delta }}_{\mathrm{a}}^{-\zeta }\mathrm{x}\left(\mathfrak{t}\right)={\displaystyle \frac{\mathfrak{h}}{\mathrm{\Gamma }\left(\zeta \right)}}\begin{array}{c}\hfill {\displaystyle \sum _{s=\frac{\mathrm{a}}{\mathfrak{h}}}^{\frac{\mathfrak{t}}{\mathfrak{h}}-\zeta }}{(\mathfrak{t}-\sigma \left(s\mathfrak{h}\right))}_{\mathfrak{h}}^{(\zeta -1)}\end{array}\mathrm{x}\left(s\mathfrak{h}\right),$$

$$\sigma \left(s\mathfrak{h}\right)=(s+1)\mathfrak{h},\mathfrak{t}\in {\left(\mathfrak{h}\mathbb{N}\right)}_{\begin{array}{c}\hfill \mathrm{a}\end{array}+\zeta \mathfrak{h}}.$$

where $\begin{array}{c}\hfill \mathrm{a}\end{array}\in \mathbb{R}$ is the initial value and the $\mathfrak{h}$-falling factorial function is stated as

$${\mathfrak{t}}_{\mathfrak{h}}^{\left(\zeta \right)}={\mathfrak{h}}^{\zeta }{\displaystyle \frac{\mathrm{\Gamma }(\frac{\mathfrak{t}}{\mathfrak{h}}+1)}{\mathrm{\Gamma }(\frac{\mathfrak{t}}{\mathfrak{h}}+1-\zeta )}}.$$

while

$${\left(\mathfrak{h}\mathbb{N}\right)}_{\begin{array}{c}\hfill \mathrm{a}\end{array}+\zeta \mathfrak{h}}=\{\mathrm{a}+(1-\zeta )\mathfrak{h},\mathrm{a}+(2-\zeta )\mathfrak{h},...\}.$$

Definition 3

([22]). For $\begin{array}{c}\hfill \mathrm{x}\end{array}\left(\mathfrak{t}\right)$ given on ${\left(\mathfrak{h}\mathbb{N}\right)}_{\begin{array}{c}\hfill \mathrm{a}\end{array}}$ and a stated $0<\zeta <1$, the Caputo $\mathfrak{h}$-difference operator is supplied:

$${}_{\mathfrak{h}}^C{\mathrm{\Delta }}_{\mathrm{a}}^{\zeta }\mathrm{x}\left(\mathfrak{t}\right){=}_{\mathfrak{h}}{\mathrm{\Delta }}_{\mathrm{a}}^{-(n-\zeta )}{\mathrm{\Delta }}_{\mathfrak{h}}^n\mathrm{x}\left(\mathfrak{t}\right).$$

where ${\mathrm{\Delta }}_{\mathfrak{h}}^n\mathrm{x}\left(\mathfrak{t}\right)={\displaystyle \frac{\mathrm{x}(\mathfrak{t}+\mathfrak{h})-\mathrm{x}\left(\mathfrak{t}\right)}{\mathfrak{h}}}.$

The following are a few essential properties that have been used in this work.

Lemma 1

([22]). For $t\in {\left(h\mathbb{N}\right)}_{\mathrm{a}+\zeta \mathfrak{h}}$ and $0<\zeta \le 1,$ the following proprieties hold:

•  

$${}_{\mathfrak{h}}{\mathrm{\Delta }}_{\mathrm{a}+(1-\zeta )\mathfrak{h}}^{-\zeta }{}_{\mathfrak{h}}^C{\mathrm{\Delta }}^{\zeta }\begin{array}{c}\hfill \mathrm{x}\end{array}\left(\mathfrak{t}\right)=\mathrm{x}\left(\mathfrak{t}\right)-\mathrm{x}\left(\mathrm{a}\right).$$

• For a constant$\begin{array}{c}\hfill \mathrm{x}\end{array}$

$${}_{\mathfrak{h}}^C{\mathrm{\Delta }}^{\zeta }\mathrm{x}=0.$$

Lemma 2

([22]). For $t\in {\left(h\mathbb{N}\right)}_{\mathrm{a}+\zeta \mathfrak{h}},$ inequality (2) holds true.

$${}_{\mathfrak{h}}^C{\mathrm{\Delta }}_{\mathrm{a}}^{\zeta }{\mathrm{x}}^2\left(\mathfrak{t}\right)\le 2\mathrm{x}{(\mathfrak{t}+\zeta \mathfrak{h})}_{\mathfrak{h}}^C{\mathrm{\Delta }}_{\mathrm{a}}^{\zeta }\mathrm{x}\left(\mathfrak{t}\right),$$

(2)

where $0<\zeta \le 1.$

Considering the fractional-order difference system:

$${}_{\mathfrak{h}}^C{\mathrm{\Delta }}_{\mathrm{a}}^{\zeta }\mathrm{x}\left(\mathfrak{t}\right)=\mathrm{\Phi }(\mathfrak{t}+\mathfrak{h}\zeta ,\mathrm{x}(\mathfrak{t}+\mathfrak{h}\zeta )),\mathfrak{t}\in {\left(\mathfrak{h}\mathbb{N}\right)}_{\mathrm{a}+\zeta \mathfrak{h}},$$

(3)

Theorem 2

([22]). Suppose $\mathrm{x}=0$ is the equilibrium point of system (3). If a positively definite and decreasing scalar function $V(\mathfrak{t},\mathrm{x}\left(\mathfrak{t}\right))$ exists so that ${}_{\mathfrak{h}}^C{\mathrm{\Delta }}_{\mathrm{a}}^{\zeta }V(\mathfrak{t},\mathrm{x}\left(\mathfrak{t}\right))\le 0$, the equilibrium point is asymptotically stable.

3. Model Description

The models in question are now approximated using two well-known approaches. These discrete models are, to our knowledge, the first in the literature. Wu et al. [32] proposed an interesting discretization of the fractional reaction equation shown below.

$$\left\{\begin{array}{c}{\displaystyle \frac{\partial \varkappa }{\partial \mathfrak{t}}}=K\mathrm{\Delta }\varkappa ,\mathrm{x}\in \mathrm{\Omega },\mathfrak{t}>0,\hfill \\ {\partial }_{\varkappa }=0,\mathrm{x}\in \partial \mathrm{\Omega },\mathfrak{t}>0,\hfill \\ \varkappa (\mathrm{x},0)=\varkappa \left(\mathrm{x}\right)>0,\mathrm{x}\in \mathrm{\Omega }.\hfill \end{array}\right.$$

(4)

This equation represents a classical diffusion equation with the initial boundary conditions, 0 is the initial point, and K is the diffusion coefficient.

According to the structure of the model (4) and the discretization employed by Wu et al. [32,33]. Considering $\mathrm{x}\in [0,L]$, we obtain ${\mathrm{x}}_{i+1}={\mathrm{x}}_i+{∆}_{\mathrm{x}},i=0,...,m,$ and by applying the central difference formula for $\mathrm{x}$, ${\displaystyle \frac{{\partial }^{2}u(\mathrm{x},\mathfrak{t})}{\partial {\mathrm{x}}^2}}$ as well as ${\displaystyle \frac{{\partial }^{2}w(\mathrm{x},\mathfrak{t})}{\partial {\mathrm{x}}^2}}$ may be approximated as

$$\left\{\begin{array}{c}{\displaystyle \frac{{\partial }^{2}u(\mathrm{x},\mathfrak{t})}{\partial {\mathrm{x}}^2}}\approx {\displaystyle \frac{u_{i+1}\left(\mathfrak{t}\right)-2u_i\left(\mathfrak{t}\right)+u_{i-1}\left(\mathfrak{t}\right)}{{\mathrm{\Delta }}_{\mathrm{x}}^2}},\hfill \\ {\displaystyle \frac{{\partial }^{2}w(\mathrm{x},\mathfrak{t})}{\partial {\mathrm{x}}^2}}\approx {\displaystyle \frac{w_{i+1}\left(\mathfrak{t}\right)-2w_i\left(\mathfrak{t}\right)+w_{i-1}\left(\mathfrak{t}\right)}{{\mathrm{\Delta }}_{\mathrm{x}}^2}}.\hfill \end{array}\right.$$

With the aid of the description of the second-order difference operator of ${\begin{array}{c}\hfill u\end{array}}_i$ and $w_i$, we obtain:

$$\left\{\begin{array}{c}{\displaystyle \frac{{\partial }^{2}u(\mathrm{x},t)}{\partial {\mathrm{x}}^2}}\approx {\displaystyle \frac{{\mathrm{\Delta }}^2{u}_{i-1}\left(\mathfrak{t}\right)}{{\mathrm{\Delta }}_{\mathrm{x}}^2}},\hfill \\ {\displaystyle \frac{{\partial }^{2}w(\mathrm{x},t)}{\partial {\mathrm{x}}^2}}\approx {\displaystyle \frac{{\mathrm{\Delta }}^2{w}_{i-1}\left(\mathfrak{t}\right)}{{\mathrm{\Delta }}_{\mathrm{x}}^2}}.\hfill \end{array}\right.$$

As a result, we may identify the previously mentioned model by Wu et al. [32,33].

$${}_{\mathfrak{h}}^C{\mathrm{\Delta }}_{{\mathfrak{t}}_0}^{\zeta }{\varkappa }_i\left(\mathfrak{t}\right)={\displaystyle \frac{1}{{\mathrm{\Delta }}_{\mathrm{x}}^2}}{\mathrm{\Delta }}^2{\varkappa }_{i-1}(\mathfrak{t}+\mathfrak{h}\zeta ).$$

(5)

With the periodic boundary conditions

$${\varkappa }_0\left(\mathfrak{t}\right)={\varkappa }_m\left(\mathfrak{t}\right),{\varkappa }_1\left(\mathfrak{t}\right)={\varkappa }_{m+1}\left(\mathfrak{t}\right).$$

(6)

Moving on to the models in question, the Degn–Harrison model and the Lengyel–Epstein model, we present the discrete fractional version of each.

Mesdoui et al. [1] designed the reaction-diffusion model commonly referred to as the Degn–Harrison reaction-diffusion model, which is represented as

$$\left\{\begin{array}{c}{\displaystyle \frac{\partial u}{\partial \mathfrak{t}}}=k_1\mathrm{\Delta }u+\mathrm{a}-u-{\displaystyle \frac{wu}{1+q{u}^2}},\mathrm{x}\in \mathrm{\Omega },\mathfrak{t}>0,\hfill \\ {\displaystyle \frac{\partial w}{\partial \mathfrak{t}}}=k_2\mathrm{\Delta }w+\mathrm{b}-{\displaystyle \frac{wu}{1+q{u}^2}},\mathrm{x}\in \mathrm{\Omega },\mathfrak{t}>0,\hfill \\ {\partial }_u={\partial }_w=0,x\in \partial \mathrm{\Omega },\mathfrak{t}>0\hfill \\ \begin{array}{c}\hfill u(\mathrm{x},0)=u_0\left(\mathrm{x}\right)>0,w(\mathrm{x},0)=w_0\left(\mathrm{x}\right)>0,\mathrm{x}\in \mathrm{\Omega }.\end{array}\hfill \end{array}\right.$$

(7)

where $\mathrm{\Omega }$ is a bounded domain in ${\mathbb{R}}^n$, and $\partial \mathrm{\Omega }$ is a suitably smooth border, while $k_1$ and $k_2$ are the respective diffusion coefficients of the reacting substances $u$ and w, which are supposed to be positive constants throughout the reaction phase.The Laplace operator is given by $\mathrm{\Delta }={\sum }_{i=1}^n{\displaystyle \frac{{\partial }^2}{\partial {\mathrm{x}}_i^2}}$.

Because time-fractional systems have been widely explored by scholars, Mesdoui et al. [1] presented the following fractional-time Degn–Harrison reaction-diffusion system.

$$\left\{\begin{array}{c}{}_0^C{D}_{\mathfrak{t}}^{\delta }u-k_1\mathrm{\Delta }u=\mathrm{a}-u-{\displaystyle \frac{uw}{1+q{u}^2}},\hfill \\ {}_0^C{D}_{\mathfrak{t}}^{\delta }w-k_2\mathrm{\Delta }w=\mathrm{b}-{\displaystyle \frac{uw}{1+q{u}^2}}.\hfill \end{array}\right.$$

(8)

where $0<\delta \le 1$ is the fractional order, ${}_0^C{D}_{\mathfrak{t}}^{\delta }$ denotes the Caputo fractional derivative, $k_1,k_2$ and $\sigma$ are strictly positive constants with the same initial conditions and Neumann boundary conditions considered by Mesdoui et al. [34].

Following the discretization defined previously, we may now provide the discrete fractional reaction-diffusion Degn–Harrison system.

$$\left\{\begin{array}{c}{}_{\mathfrak{h}}^C{\mathrm{\Delta }}_{{\mathfrak{t}}_0}^{\zeta }{u}_i\left(\mathfrak{t}\right)={\displaystyle \frac{k_1}{{\mathrm{\Delta }}_{\mathrm{x}}^2}}{\mathrm{\Delta }}^2{u}_{i-1}(\mathfrak{t}+\mathfrak{h}\zeta )+\mathrm{a}-u_i(\mathfrak{t}+\mathfrak{h}\zeta )-{\displaystyle \frac{u_i(\mathfrak{t}+\mathfrak{h}\zeta )w_i(\mathfrak{t}+\mathfrak{h}\zeta )}{1+q{\left(u_i(\mathfrak{t}+\mathfrak{h}\zeta )\right)}^2}},\hfill \\ {}_{\mathfrak{h}}^C{\mathrm{\Delta }}_{{\mathfrak{t}}_0}^{\zeta }{w}_i\left(\mathfrak{t}\right)={\displaystyle \frac{k_2}{{\mathrm{\Delta }}_{\mathrm{x}}^2}}{\mathrm{\Delta }}^2{w}_{i-1}(\mathfrak{t}+\mathfrak{h}\zeta )+\mathrm{b}-{\displaystyle \frac{u_i(\mathfrak{t}+\mathfrak{h}\zeta )w_i(\mathfrak{t}+\mathfrak{h}\zeta )}{1+q{\left(u_i(\mathfrak{t}+\mathfrak{h}\zeta )\right)}^2}}.\hfill \end{array}\right.$$

(9)

With the periodic boundary conditions

$$\left\{\begin{array}{c}{u}_0\left(\mathfrak{t}\right)=u_m\left(\mathfrak{t}\right),u_1\left(\mathfrak{t}\right)=u_{m+1}\left(\mathfrak{t}\right),\hfill \\ w_0\left(\mathfrak{t}\right)=w_m\left(\mathfrak{t}\right),w_1\left(\mathfrak{t}\right)=w_{m+1}\left(\mathfrak{t}\right),\hfill \end{array}\right.$$

(10)

and the initial condition

$$u_i\left({\mathfrak{t}}_0\right)={\psi }_1\left({\mathrm{x}}_i\right)\ge 0,w_i\left({\mathfrak{t}}_0\right)={\psi }_2\left({\mathrm{x}}_i\right)\ge 0.$$

Regarding the remaining model, the Lengyel–Epstein reaction-diffusion system was provided as a simulation of the chlorite-iodide-malonic-acid chemical reaction (CIMA). Yi et al. [35] investigated a specific model described by:

$$\left\{\begin{array}{c}{\displaystyle \frac{\partial u}{\partial \mathfrak{t}}}=\mathrm{\Delta }u+\mathrm{a}-u-{\displaystyle \frac{uw}{1+u^2}},\mathrm{x}\in \mathrm{\Omega },\mathfrak{t}>0,\hfill \\ {\displaystyle \frac{\partial w}{\partial \mathfrak{t}}}=\sigma \left(d\mathrm{\Delta }w+\mathrm{b}\left(u-{\displaystyle \frac{uw}{1+u^2}}\right)\right),\mathrm{x}\in \mathrm{\Omega },\mathfrak{t}>0,\hfill \\ {\partial }_u={\partial }_w=0,\mathrm{x}\in \partial \mathrm{\Omega },\mathfrak{t}>0,\hfill \\ u(\mathrm{x},0)=u_0\left(\mathrm{x}\right)>0,w(\mathrm{x},0)=w_0\left(\mathrm{x}\right)>0,\mathrm{x}\in \mathrm{\Omega }.\hfill \end{array}\right.$$

(11)

where $\mathrm{\Omega }$ is a bounded domain in ${\mathbb{R}}^n$, with a properly smooth boundary $\partial \mathrm{\Omega }$. $u$ reflects the chemical concentration of the activator iodide, whereas w indicates the inhibitor chlorite at a point $\mathrm{x}\in \mathrm{\Omega }$, $\mathrm{a}$ and $\begin{array}{c}\hfill \mathrm{b}\end{array}$ are related to the supply concentration, d is the value of the ratio of the coefficient of diffusion, and $\sigma >0$ is an adjusting parameter determined by the amount of starch concentration.

Given that fractional systems have been thoroughly studied over the years, the next fractional Lengyel–Epstein system was investigated:

$$\left\{\begin{array}{c}{}_0^C{D}_{\mathfrak{t}}^{\delta }u-d_1\mathrm{\Delta }u=\mathrm{a}-u-{\displaystyle \frac{4uw}{1+u^2}},\hfill \\ {}_0^C{D}_{\mathfrak{t}}^{\delta }w-d_2\mathrm{\Delta }w=\sigma \mathrm{b}\left(u-{\displaystyle \frac{uw}{1+u^2}}\right),\hfill \end{array}\right.$$

(12)

where $d_1,d_2$ and $\sigma$ are constants that are positive and have similar initial conditions and Neumann boundaries.

We analyze the discrete fractional reaction-diffusion Lengyel–Epstein system (15) via the model (12) and the discretization described above.

$$\left\{\begin{array}{c}{}_{\mathfrak{h}}^C{\mathrm{\Delta }}_{{\mathfrak{t}}_0}^{\zeta }{u}_i\left(\mathfrak{t}\right)={\displaystyle \frac{d_1}{{\mathrm{\Delta }}_{\mathrm{x}}^2}}{\mathrm{\Delta }}^2{u}_{i-1}(\mathfrak{t}+\mathfrak{h}\zeta )+\mathrm{a}-u_i(\mathfrak{t}+\mathfrak{h}\zeta )-{\displaystyle \frac{4u_i(\mathfrak{t}+\mathfrak{h}\zeta )w_i(\mathfrak{t}+\mathfrak{h}\zeta )}{1+{\left(u_i(\mathfrak{t}+\mathfrak{h}\zeta )\right)}^2}},\hfill \\ {}_{\mathfrak{h}}^C{\mathrm{\Delta }}_{{\mathfrak{t}}_0}^{\zeta }{w}_i\left(\mathfrak{t}\right)={\displaystyle \frac{d_2}{{\mathrm{\Delta }}_{\mathrm{x}}^2}}{\mathrm{\Delta }}^2{w}_{i-1}(\mathfrak{t}+\mathfrak{h}\zeta )+\sigma \mathrm{b}\left(u_i(\mathfrak{t}+\mathfrak{h}\zeta )-{\displaystyle \frac{u_i(\mathfrak{t}+\mathfrak{h}\zeta )w_i(\mathfrak{t}+\mathfrak{h}\zeta )}{1+{\left(u_i(\mathfrak{t}+\mathfrak{h}\zeta )\right)}^2}}\right).\hfill \end{array}\right.$$

(13)

Using periodic boundary conditions:

$$\left\{\begin{array}{c}{u}_0\left(\mathfrak{t}\right)=u_m\left(\mathfrak{t}\right),u_1\left(\mathfrak{t}\right)=u_{m+1}\left(\mathfrak{t}\right),\hfill \\ w_0\left(\mathfrak{t}\right)=w_m\left(\mathfrak{t}\right),w_1\left(\mathfrak{t}\right)=w_{m+1}\left(\mathfrak{t}\right),\hfill \end{array}\right.$$

(14)

as well as the initial condition

$$u_i\left({\mathfrak{t}}_0\right)={\psi }_1\left({\mathrm{x}}_i\right)\ge 0,w_i\left({\mathfrak{t}}_0\right)={\psi }_2\left({\mathrm{x}}_i\right)\ge 0.$$

4. Synchronization of Discrete-Time Fractional Reaction-Diffusion Lengyel–Epstein System

The most typical method for testing synchronization is to employ a controller to have the slave system output duplicate the master system output in some similar way. In this part, we create a controller that minimizes the state difference between synchronized systems to zero, which is known to be complete synchronization. Let the discrete reaction-diffusion master system (13) and the slave system be

$$\left\{\begin{array}{c}{}_{\mathfrak{h}}^C{\mathrm{\Delta }}_{{\mathfrak{t}}_0}^{\zeta }{U}_i\left(t\right)={\displaystyle \frac{d_1}{{\mathrm{\Delta }}_{\mathrm{x}}^2}}{\mathrm{\Delta }}^2{U}_{i-1}(\mathfrak{t}+\mathfrak{h}\zeta )+\mathrm{a}-U_i(\mathfrak{t}+\mathfrak{h}\zeta )-{\displaystyle \frac{4U_i(\mathfrak{t}+\mathfrak{h}\zeta )W_i(\mathfrak{t}+\mathfrak{h}\zeta )}{1+{\left(U_i(\mathfrak{t}+\mathfrak{h}\zeta )\right)}^2}}+V_1\left(\mathfrak{t}\right),\hfill \\ {}_{\mathfrak{h}}^C{\mathrm{\Delta }}_{{\mathfrak{t}}_0}^{\zeta }{W}_i\left(\mathfrak{t}\right)={\displaystyle \frac{d_2}{{\mathrm{\Delta }}_{\mathrm{x}}^2}}{\mathrm{\Delta }}^2{W}_{i-1}(\mathfrak{t}+\mathfrak{h}\zeta )+\sigma \mathrm{b}\left(U_i(\mathfrak{t}+\mathfrak{h}\zeta )-{\displaystyle \frac{U_i(\mathfrak{t}+\mathfrak{h}\zeta )W_i(\mathfrak{t}+\mathfrak{h}\zeta )}{1+{\left(U_i(\mathfrak{t}+\mathfrak{h}\zeta )\right)}^2}}\right)+V_2\left(\mathfrak{t}\right).\hfill \end{array}\right.$$

(15)

With the periodic boundary conditions

$$\left\{\begin{array}{c}{U}_0\left(\mathfrak{t}\right)=U_m\left(\mathfrak{t}\right),U_1\left(\mathfrak{t}\right)=U_{m+1}\left(\mathfrak{t}\right),\hfill \\ W_0\left(\mathfrak{t}\right)=W_m\left(\mathfrak{t}\right),W_1\left(\mathfrak{t}\right)=W_{m+1}\left(\mathfrak{t}\right),\hfill \end{array}\right.$$

(16)

and the initial condition

$$U_i\left({\mathfrak{t}}_0\right)={\mathrm{\Phi }}_1\left({\mathrm{x}}_i\right)\ge 0,W_i\left({\mathfrak{t}}_0\right)={\mathrm{\Phi }}_2\left({\mathrm{x}}_i\right)\ge 0.$$

The purpose of synchronization is to reduce the error regarding the master and slave systems to zero, which is described as

$$(e_{1i},e_{2i})=(U_i-u_i,W_i-w_i).$$

(17)

In what follows, we will identify the linear controllers $V_1$ and $V_2$ that cause the error system solution to be 0 as $\mathfrak{t}$ approaches $+\infty$. In other words, to be able to accomplish complete synchronization within the master-slave systems (13)–(15), we examine the asymptotical stability of the zero solution of the synchronization error system described in (17).

First, According to Ouannas et al. [2], it is easy to verify that

$$\begin{array}{cc}\hfill & |{\displaystyle \frac{4U_i(\mathfrak{t}+\mathfrak{h}\zeta )W_i(\mathfrak{t}+\mathfrak{h}\zeta )}{1+{\left(U_i(\mathfrak{t}+\mathfrak{h}\zeta )\right)}^2}}-{\displaystyle \frac{4u_i(\mathfrak{t}+\mathfrak{h}\zeta )w_i(\mathfrak{t}+\mathfrak{h}\zeta )}{1+{\left(u_i(\mathfrak{t}+\mathfrak{h}\zeta )\right)}^2}}|\le |U_i-u_i|+4|W_i-w_i|,\hfill \\ \hfill & |{\displaystyle \frac{\sigma b{U}_i(\mathfrak{t}+\mathfrak{h}\zeta )W_i(\mathfrak{t}+\mathfrak{h}\zeta )}{1+{\left(U_i(\mathfrak{t}+\mathfrak{h}\zeta )\right)}^2}}-{\displaystyle \frac{\sigma \mathrm{b}{u}_i(\mathfrak{t}+\mathfrak{h}\zeta )w_i(\mathfrak{t}+\mathfrak{h}\zeta )}{1+{\left(u_i(\mathfrak{t}+\mathfrak{h}\zeta )\right)}^2}}|\le |U_i-u_i|+\sigma b|W_i-w_i|.\hfill \end{array}$$

Theorem 3.

If there is a control matrix $M={\left(m_{ij}\right)}_{2\times 2}$ that satisfies $1-m_1>0$ and $m_2-\sigma b>0$ the master-slave reaction-diffusion system identified in (13)–(15) is synchronized applying the linear control rule indicated below.

$$\left\{\begin{array}{c}{V}_1\left(\mathfrak{t}\right)=-\left(m_1+{\displaystyle \frac{29}{4}}\right)e_{1i}\left(\mathfrak{t}\right),\hfill \\ V_2\left(\mathfrak{t}\right)=-m_2{e}_{2i}-(\sigma b+1)e_{2i}\left(\mathfrak{t}\right).\hfill \end{array}\right.$$

(18)

Proof. 

When (18) is substituted into the error system described in (17), the result is

$$\left\{\begin{array}{c}{}_{\mathfrak{h}}^C{\mathrm{\Delta }}_{{\mathfrak{t}}_0}^{\zeta }{e}_{1i}\left(\mathfrak{t}\right)={\displaystyle \frac{d_1}{{\mathrm{\Delta }}_{\mathrm{x}}^2}}{\mathrm{\Delta }}^2{e}_{1,i-1}(\mathfrak{t}+\mathfrak{h}\zeta )-e_{1i}(\mathfrak{t}+h\zeta )\hfill \\ -({\displaystyle \frac{4U_i(\mathfrak{t}+\mathfrak{h}\zeta )W_i(\mathfrak{t}+\mathfrak{h}\zeta )}{1+{\left(U_i(\mathfrak{t}+\mathfrak{h}\zeta )\right)}^2}}-{\displaystyle \frac{4u_i(\mathfrak{t}+\mathfrak{h}\zeta )w_i(\mathfrak{t}+\mathfrak{h}\zeta )}{1+{\left(u_i(\mathfrak{t}+\mathfrak{h}\zeta )\right)}^2}})-\left(m_1+{\displaystyle \frac{29}{4}}\right)e_{1i}\left(\mathfrak{t}\right),\hfill \\ {}_{\mathfrak{h}}^C{\mathrm{\Delta }}_{{\mathfrak{t}}_0}^{\zeta }{e}_{2i}\left(\mathfrak{t}\right)={\displaystyle \frac{d_2}{{\mathrm{\Delta }}_{\mathrm{x}}^2}}{\mathrm{\Delta }}^2{e}_{2,i-1}(\mathfrak{t}+\mathfrak{h}\zeta )+\sigma b(e_{2i}(\mathfrak{t}+\mathfrak{h}\zeta )\hfill \\ -({\displaystyle \frac{U_i(\mathfrak{t}+\mathfrak{h}\zeta )W_i(\mathfrak{t}+\mathfrak{h}\zeta )}{1+{\left(U_i(\mathfrak{t}+\mathfrak{h}\zeta )\right)}^2}}-{\displaystyle \frac{u_i(\mathfrak{t}+\mathfrak{h}\zeta )w_i(\mathfrak{t}+\mathfrak{h}\zeta )}{1+{\left(u_i(\mathfrak{t}+\mathfrak{h}\zeta )\right)}^2}}))-m_2{e}_{2i}-(\sigma \mathrm{b}+1)e_{2i}\left(\mathfrak{t}\right).\hfill \end{array}\right.$$

(19)

Developing a Lyapunov function of the type

$$L\left(\mathfrak{t}\right)={\displaystyle \frac{1}{2}}\sum _{i=1}^m\left(e_{1i}^2\left(\mathfrak{t}\right)+e_{2i}^2\left(\mathfrak{t}\right)\right),$$

(20)

gives

$$\begin{array}{cc}\hfill {}_{\mathfrak{h}}^C{\mathrm{\Delta }}_{{\mathfrak{t}}_0}^{\zeta }L\left(\mathfrak{t}\right)& \le \sum _{i=1}^m{e}_{1i}{(\mathfrak{t}+\mathfrak{h}\zeta )}_h^C{\mathrm{\Delta }}_{{\mathfrak{t}}_0}^{\zeta }{e}_{1i}\left(\mathfrak{t}\right)+e_{2i}{(\mathfrak{t}+\mathfrak{h}\zeta )}_h^C{\mathrm{\Delta }}_{{\mathfrak{t}}_0}^{\zeta }{e}_{2i}\left(\mathfrak{t}\right),\hfill \\ \hfill & =\sum _{i=1}^m{e}_{1i}(\mathfrak{t}+h\zeta )[{\displaystyle \frac{d_1}{{\mathrm{\Delta }}_{\mathrm{x}}^2}}{\mathrm{\Delta }}^2{e}_{1,i-1}(\mathfrak{t}+\mathfrak{h}\zeta )-e_{1i}(\mathfrak{t}+\mathfrak{h}\zeta )-({\displaystyle \frac{4U_i(\mathfrak{t}+\mathfrak{h}\zeta )W_i(\mathfrak{t}+\mathfrak{h}\zeta )}{1+{\left(U_i(\mathfrak{t}+\mathfrak{h}\zeta )\right)}^2}}\hfill \\ \hfill & -{\displaystyle \frac{4u_i(\mathfrak{t}+\mathfrak{h}\zeta )w_i(\mathfrak{t}+\mathfrak{h}\zeta )}{1+{\left(u_i(\mathfrak{t}+\mathfrak{h}\zeta )\right)}^2}}-\left(m_1+{\displaystyle \frac{29}{4}}\right)e_{1i}]+\sum _{i=1}^m{e}_{2i}(\mathfrak{t}+\mathfrak{h}\zeta )[{\displaystyle \frac{d_2}{{\mathrm{\Delta }}_{\mathrm{x}}^2}}{\mathrm{\Delta }}^2{e}_{2,i-1}(\mathfrak{t}+\mathfrak{h}\zeta )\hfill \\ \hfill & +\sigma \mathrm{b}{e}_{2i}(\mathfrak{t}+\mathfrak{h}\zeta )-\left({\displaystyle \frac{\sigma b{U}_i(\mathfrak{t}+\mathfrak{h}\zeta )W_i(\mathfrak{t}+\mathfrak{h}\zeta )}{1+{\left(U_i(\mathfrak{t}+\mathfrak{h}\zeta )\right)}^2}}-{\displaystyle \frac{\sigma \mathrm{b}{u}_i(\mathfrak{t}+\mathfrak{h}\zeta )w_i(\mathfrak{t}+\mathfrak{h}\zeta )}{1+{\left(u_i(\mathfrak{t}+\mathfrak{h}\zeta )\right)}^2}}\right)\hfill \\ \hfill & -(m_2{e}_{2i}+(\sigma \mathrm{b}+1)e_{2i})],\hfill \\ \hfill & \le \sum _{i=1}^m{\displaystyle \frac{d_1}{{\mathrm{\Delta }}_{\mathrm{x}}^2}}{e}_{1i}(\mathfrak{t}+\mathfrak{h}\zeta ){\mathrm{\Delta }}^2{e}_{1,i-1}(\mathfrak{t}+\mathfrak{h}\zeta )+{\displaystyle \frac{d_2}{{\mathrm{\Delta }}_{\mathrm{x}}^2}}{e}_{2i}(\mathfrak{t}+\mathfrak{h}\zeta ){\mathrm{\Delta }}^2{e}_{2,i-1}(\mathfrak{t}+\mathfrak{h}\zeta )\hfill \\ \hfill & -\sum _{i=1}^m(1-m_1)e_{1i}^2\left(\mathfrak{t}\right)+\left|e_{1i}\left(\mathfrak{t}\right)\right|\left(|U_i-u_i|+4|W_i-w_i|\right)\hfill \\ \hfill & -{\displaystyle \frac{29}{4}}{e}_{1i}^2(\mathfrak{t}+\mathfrak{h}\zeta )+\sum _{i=1}^m(\sigma \mathrm{b}-m_2)e_{2i}^2(\mathfrak{t}+\mathfrak{h}\zeta )+|e_{2i}(\mathfrak{t}+\mathfrak{h}\zeta )\left|\right(|U_i(\mathfrak{t}+\mathfrak{h}\zeta )-u_i(\mathfrak{t}+\mathfrak{h}\zeta )|\hfill \\ \hfill & +\sigma b|W_i(\mathfrak{t}+\mathfrak{h}\zeta )-w_i(\mathfrak{t}+\mathfrak{h}\zeta )\left|\right)-(\sigma \mathrm{b}+1)e_{2i}^2(\mathfrak{t}+\mathfrak{h}\zeta ),\hfill \\ \hfill & \le {\displaystyle \frac{d_1}{{\mathrm{\Delta }}_{\mathrm{x}}^2}}(\mathrm{\Delta }{e}_{1,i-1}\mathrm{\Delta }{e}_{1,i-1}(\mathfrak{t}+\mathfrak{h}\zeta ){|}_{m+1}^1-\sum _{i=1}^m{\left(\mathrm{\Delta }{e}_{1,i-1}(\mathfrak{t}+\mathfrak{h}\zeta )\right)}^2)+{\displaystyle \frac{d_2}{{\mathrm{\Delta }}_x^2}}{\left(\mathrm{\Delta }{e}_{2,i-1}\mathrm{\Delta }{e}_{2,i-1}(\mathfrak{t}+\mathfrak{h}\zeta )\right|}_{m+1}^1\hfill \\ \hfill & -\sum _{i=1}^m{\left(\mathrm{\Delta }{e}_{2,i-1}(\mathfrak{t}+\mathfrak{h}\zeta )\right)}^2)-\sum _{i=1}^m(1-m_1)e_{1i}^2\left(\mathfrak{t}\right)+(\sigma \mathrm{b}-m_2)e_{2i}^2(\mathfrak{t}+\mathfrak{h}\zeta )\hfill \\ \hfill & +e_{1i}^2\left(t\right)+4|e_{1i}\left(\mathfrak{t}\right)\left|\right|e_{2i}\left(\mathfrak{t}\right)|-{\displaystyle \frac{29}{4}}{e}_{1i}^2(\mathfrak{t}+\mathfrak{h}\zeta )+|e_{1i}(\mathfrak{t}+\mathfrak{h}\zeta )\left|\right|e_{2i}(\mathfrak{t}+\mathfrak{h}\zeta )|+\sigma \mathrm{b}{e}_{2i}^2(\mathfrak{t}+\mathfrak{h}\zeta )\hfill \\ \hfill & -(\sigma \mathrm{b}+1)e_{2i}^2(\mathfrak{t}+\mathfrak{h}\zeta ),\hfill \\ \hfill & \le -\sum _{i=1}^m\left({\displaystyle \frac{d_1}{{\mathrm{\Delta }}_{\mathrm{x}}^2}}{\left(\mathrm{\Delta }{e}_{1,i-1}(\mathfrak{t}+\mathfrak{h}\zeta )\right)}^2+{\displaystyle \frac{d_2}{{\mathrm{\Delta }}_{\mathrm{x}}^2}}{\left(\mathrm{\Delta }{e}_{2,i-1}(\mathfrak{t}+\mathfrak{h}\zeta )\right)}^2\right)-\sum _{i=1}^m(1-m_1)e_{1i}^2\left(\mathfrak{t}\right)\hfill \\ \hfill & -\sum _{i=1}^m(m_2-\sigma \mathrm{b})e_{2i}^2\left(t\right)-\sum _{i=1}^m{\left({\displaystyle \frac{5}{2}}|e_{1i}(\mathfrak{t}+\mathfrak{h}\zeta )|-|e_{2i}(\mathfrak{t}+\mathfrak{h}\zeta )|\right)}^2<0.\hfill \end{array}$$

This means the global asymptotic stability of the error system’s zero solution (19), based on the Lyapunov stability theory presented in Theorem 2. As a result, the master-slave systems (13)–(15) are completely synchronized. □

5. Synchronization of Discrete Fractional Degn–Harrison Reaction-Diffusion Systems

We investigate the synchronization of the fractional discrete Degn–Harrison models using the master-slave formalism, in which the discrete fractional Degn–Harrison reaction-diffusion systems are linked in such a way that the slave system asymptotically matches the master system. In this scenario, we create controllers that cause the difference between the states of synchronized systems to converge to zero, indicating that the systems are fully synchronized. The slave system that is linked to the master system (13) may be expressed as

$$\left\{\begin{array}{c}{}_{\mathfrak{h}}^C{\mathrm{\Delta }}_{{\mathfrak{t}}_0}^{\zeta }{U}_i\left(\mathfrak{t}\right)={\displaystyle \frac{k_1}{{\mathrm{\Delta }}_{\mathrm{x}}^2}}{\mathrm{\Delta }}^2{U}_{i-1}(\mathfrak{t}+\mathfrak{h}\zeta )+\mathrm{a}-U_i(\mathfrak{t}+\mathfrak{h}\zeta )-{\displaystyle \frac{U_i(\mathfrak{t}+\mathfrak{h}\zeta )W_i(\mathfrak{t}+\mathfrak{h}\zeta )}{1+q{\left(U_i(\mathfrak{t}+h\zeta )\right)}^2}}+S_1\left(\mathfrak{t}\right),\hfill \\ {}_{\mathfrak{h}}^C{\mathrm{\Delta }}_{{\mathfrak{t}}_0}^{\zeta }{W}_i\left(\mathfrak{t}\right)={\displaystyle \frac{k_2}{{\mathrm{\Delta }}_{\mathrm{x}}^2}}{\mathrm{\Delta }}^2{W}_{i-1}(\mathfrak{t}+\mathfrak{h}\zeta )+\mathrm{b}-{\displaystyle \frac{U_i(\mathfrak{t}+\mathfrak{h}\zeta )W_i(\mathfrak{t}+\mathfrak{h}\zeta )}{1+q{\left(U_i(\mathfrak{t}+\mathfrak{h}\zeta )\right)}^2}}+S_2\left(\mathfrak{t}\right).\hfill \end{array}\right.$$

(21)

With the periodic boundary conditions

$$\left\{\begin{array}{c}{U}_0\left(\mathfrak{t}\right)=U_m\left(\mathfrak{t}\right),U_1\left(\mathfrak{t}\right)=U_{m+1}\left(\mathfrak{t}\right),\hfill \\ W_0\left(\mathfrak{t}\right)=W_m\left(\mathfrak{t}\right),W_1\left(\mathfrak{t}\right)=W_{m+1}\left(\mathfrak{t}\right),\hfill \end{array}\right.$$

(22)

and the initial condition

$$U_i\left({\mathfrak{t}}_0\right)={\mathrm{\Phi }}_1\left({\mathrm{x}}_i\right)\ge 0,W_i\left({\mathfrak{t}}_0\right)={\mathrm{\Phi }}_2\left({\mathrm{x}}_i\right)\ge 0.$$

The goal of this part is to identify a control $S_i$ to induce the synchronization errors $e_i(\mathrm{x},\mathfrak{t})=(e_i^1(\mathrm{x},\mathfrak{t}),e_i^2(\mathrm{x},\mathfrak{t}))$ described by

$$e_{1i}(\mathrm{x},\mathfrak{t})=U_i(\mathfrak{t},\mathrm{x})-u_i(\mathfrak{t},\mathrm{x}),e_{2i}(\mathrm{x},\mathfrak{t})=W_i(\mathfrak{t},\mathrm{x})-w_i(\mathfrak{t},\mathrm{x}).$$

(23)

where $(u_i(\mathrm{x},\mathfrak{t}),w_i(\mathrm{x},\mathfrak{t}))$ and $(U(\mathrm{x},\mathfrak{t}),\begin{array}{c}\hfill W\end{array}(\mathrm{x},\mathfrak{t}))$ are the solutions of systems (13) and (21) that converge to zero as $\mathfrak{t}$ approaches infinity.

The error system is given by

$$\left\{\begin{array}{c}{}_{\mathfrak{h}}^C{\mathrm{\Delta }}_{{\mathfrak{t}}_0}^{\zeta }{e}_{1i}\left(\mathfrak{t}\right)={\displaystyle \frac{k_1}{{\mathrm{\Delta }}_{\mathrm{x}}^2}}{\mathrm{\Delta }}^2{e}_{1,i-1}(\mathfrak{t}+\mathfrak{h}\zeta )+\mathrm{a}-e_{1i}(\mathfrak{t}+\mathfrak{h}\zeta )-{\displaystyle \frac{U_i(\mathfrak{t}+\mathfrak{h}\zeta )W_i(\mathfrak{t}+\mathfrak{h}\zeta )}{1+q{\left(U_i(\mathfrak{t}+\mathfrak{h}\zeta )\right)}^2}}\hfill \\ +{\displaystyle \frac{u_i(\mathfrak{t}+\mathfrak{h}\zeta )w_i(\mathfrak{t}+\mathfrak{h}\zeta )}{1+q{\left(u_i(\mathfrak{t}+\mathfrak{h}\zeta )\right)}^2}}+S_1\left(\mathfrak{t}\right),\hfill \\ {}_{\mathfrak{h}}^C{\mathrm{\Delta }}_{{\mathfrak{t}}_0}^{\zeta }{e}_{2i}\left(t\right)={\displaystyle \frac{k_2}{{\mathrm{\Delta }}_{\mathrm{x}}^2}}{\mathrm{\Delta }}^2{e}_{2,i-1}(\mathfrak{t}+\mathfrak{h}\zeta )+\mathrm{b}-{\displaystyle \frac{U_i(\mathfrak{t}+\mathfrak{h}\zeta )W_i(\mathfrak{t}+\mathfrak{h}\zeta )}{1+q{\left(U_i(\mathfrak{t}+\mathfrak{h}\zeta )\right)}^2}}\hfill \\ +{\displaystyle \frac{u_i(\mathfrak{t}+\mathfrak{h}\zeta )w_i(\mathfrak{t}+\mathfrak{h}\zeta )}{1+q{\left(u_i(\mathfrak{t}+\mathfrak{h}\zeta )\right)}^2}}+S_2\left(\mathfrak{t}\right).\hfill \end{array}\right.$$

(24)

The error system (24) may be seen to satisfy the periodic boundary conditions.

$$\left\{\begin{array}{c}{e}_{1,0}\left(\mathfrak{t}\right)=U_0\left(\mathfrak{t}\right)-u_0\left(\mathfrak{t}\right)=U_m\left(\mathfrak{t}\right)-u_m\left(\mathfrak{t}\right)=e_{1m}\left(\mathfrak{t}\right),\hfill \\ e_{1,1}\left(\mathfrak{t}\right)=U_1\left(\mathfrak{t}\right)-u_1\left(\mathfrak{t}\right)=U_{m+1}\left(\mathfrak{t}\right)-u_{m+1}\left(\mathfrak{t}\right)=e_{1,m+1}\left(\mathfrak{t}\right),\hfill \\ e_{2,0}\left(\mathfrak{t}\right)=W_0\left(\mathfrak{t}\right)-w_0\left(\mathfrak{t}\right)=W_m\left(\mathfrak{t}\right)-w_m\left(\mathfrak{t}\right)=e_{2m}\left(\mathfrak{t}\right),\hfill \\ e_{2,1}\left(\mathfrak{t}\right)=W_1\left(\mathfrak{t}\right)-w_1\left(\mathfrak{t}\right)=W_{m+1}\left(\mathfrak{t}\right)-w_{m+1}\left(\mathfrak{t}\right)=e_{2,m+1}\left(\mathfrak{t}\right).\hfill \end{array}\right.$$

(25)

Before proceeding to the synchronization of master-slave systems, consider the following lemma.

Lemma 3

([34]). The following inequality holds

$$|{\displaystyle \frac{u_i(\mathfrak{t}+\mathfrak{h}\zeta )w_i(\mathfrak{t}+\mathfrak{h}\zeta )}{1+q{\left(u_i(\mathfrak{t}+\mathfrak{h}\zeta )\right)}^2}}-{\displaystyle \frac{U_i(\mathfrak{t}+\mathfrak{h}\zeta )W_i(\mathfrak{t}+\mathfrak{h}\zeta )}{1+q{\left(U_i(\mathfrak{t}+\mathfrak{h}\zeta )\right)}^2}}|\le Q\left(\right|U_i-u_i|+|W_i-w_i\left|\right),$$

(26)

where

$$Q\ge max\left\{{\displaystyle \frac{5}{4}}k,{\displaystyle \frac{1}{2\sqrt{q}}}\right\},\left|{\begin{array}{c}\hfill W\end{array}}_i\right|<k.$$

The controllers $S_1$ and $S_2$ are determined in the following Theorem to establish synchronization between the systems provided in (13) and (21).

Theorem 4.

Under the following control law, the master system (1) and slave system (2) are completely synchronized.

$$\left\{\begin{array}{c}{S}_1\left(\mathfrak{t}\right)=(1-2Q)e_{1i}\left(\mathfrak{t}\right),\hfill \\ S_2\left(\mathfrak{t}\right)=-2Q{e}_{2i}\left(\mathfrak{t}\right).\hfill \end{array}\right.$$

(27)

Proof. 

By substituting the control described in the Theorem in the error system, we obtain

$$\left\{\begin{array}{c}{}_{\mathfrak{h}}^C{\mathrm{\Delta }}_{{\mathfrak{t}}_0}^{\zeta }{e}_{1i}\left(\mathfrak{t}\right)={\displaystyle \frac{k_1}{{\mathrm{\Delta }}_{\mathrm{x}}^2}}{\mathrm{\Delta }}^2{e}_{1,i-1}(\mathfrak{t}+\mathfrak{h}\zeta )-{\displaystyle \frac{U_i(\mathfrak{t}+\mathfrak{h}\zeta )W_i(\mathfrak{t}+\mathfrak{h}\zeta )}{1+q{\left(U_i(\mathfrak{t}+\mathfrak{h}\zeta )\right)}^2}}+{\displaystyle \frac{u_i(\mathfrak{t}+\mathfrak{h}\zeta )w_i(\mathfrak{t}+\mathfrak{h}\zeta )}{1+q{\left(u_i(\mathfrak{t}+\mathfrak{h}\zeta )\right)}^2}}\hfill \\ -2Q{e}_{1i}(\mathfrak{t}+\mathfrak{h}\zeta ),\hfill \\ {}_{\mathfrak{h}}^C{\mathrm{\Delta }}_{{\mathfrak{t}}_0}^{\zeta }{e}_{2i}\left(\mathfrak{t}\right)={\displaystyle \frac{k_2}{{\mathrm{\Delta }}_{\mathrm{x}}^2}}{\mathrm{\Delta }}^2{e}_{2,i-1}(\mathfrak{t}+\mathfrak{h}\zeta )-{\displaystyle \frac{U_i(\mathfrak{t}+\mathfrak{h}\zeta )W_i(\mathfrak{t}+\mathfrak{h}\zeta )}{1+q{\left(U_i(\mathfrak{t}+\mathfrak{h}\zeta )\right)}^2}}+{\displaystyle \frac{u_i(\mathfrak{t}+\mathfrak{h}\zeta )w_i(\mathfrak{t}+\mathfrak{h}\zeta )}{1+q{\left(u_i(\mathfrak{t}+\mathfrak{h}\zeta )\right)}^2}}\hfill \\ -2Q{e}_{2i}(\mathfrak{t}+\mathfrak{h}\zeta ).\hfill \end{array}\right.$$

(28)

Next, we design a Lyapunov function as

$$L\left(\mathfrak{t}\right)={\displaystyle \frac{1}{2}}\sum _{i=1}^m\left({\left(e_{1i}\left(\mathfrak{t}\right)\right)}^2+{\left(e_{2i}\left(\mathfrak{t}\right)\right)}^2\right),$$

(29)

then, we have

$$\begin{array}{cc}\hfill {}_{\mathfrak{h}}^C{\mathrm{\Delta }}_{{\mathfrak{t}}_0}^{\zeta }L\left(\mathfrak{t}\right)& ={\displaystyle \frac{1}{2}}{}_{\mathfrak{h}}^C{\mathrm{\Delta }}_{{\mathfrak{t}}_0}^{\zeta }\sum _{i=1}^m{\left(e_{1i}\left(t\right)\right)}^2+{\left(e_{2i}\left(\mathfrak{t}\right)\right)}^2,\hfill \\ \hfill & \le \sum _{i=1}^m{e}_{1i}(\mathfrak{t}+\mathfrak{h}\zeta ){}_h^C{\mathrm{\Delta }}_{{\mathfrak{t}}_0}^{\zeta }{e}_{1i}\left(\mathfrak{t}\right)+e_{2i}(\mathfrak{t}+\mathfrak{h}\zeta ){}_{\mathfrak{h}}^C{\mathrm{\Delta }}_{{\mathfrak{t}}_0}^{\zeta }{e}_{2i}\left(\mathfrak{t}\right),\hfill \\ \hfill & =\sum _{i=1}^m{e}_{1i}(\mathfrak{t}+\mathfrak{h}\zeta )({\displaystyle \frac{k_1}{{\mathrm{\Delta }}_{\mathrm{x}}^2}}{\mathrm{\Delta }}^2{e}_{1,i-1}(\mathfrak{t}+\mathfrak{h}\zeta )-{\displaystyle \frac{U_i(\mathfrak{t}+\mathfrak{h}\zeta )W_i(\mathfrak{t}+\mathfrak{h}\zeta )}{1+q{\left(U_i(\mathfrak{t}+\mathfrak{h}\zeta )\right)}^2}}+{\displaystyle \frac{u_i(\mathfrak{t}+\mathfrak{h}\zeta )w_i(\mathfrak{t}+\mathfrak{h}\zeta )}{1+q{\left(u_i(\mathfrak{t}+\mathfrak{h}\zeta )\right)}^2}}\hfill \\ \hfill & -2Q{e}_i^1(\mathfrak{t}+\mathfrak{h}\zeta ))+e_{2i}(\mathfrak{t}+\mathfrak{h}\zeta )({\displaystyle \frac{k_2}{{\mathrm{\Delta }}_{\mathrm{x}}^2}}{\mathrm{\Delta }}^2{e}_{2,i-1}(\mathfrak{t}+\mathfrak{h}\zeta )\hfill \\ \hfill & -{\displaystyle \frac{U_i(\mathfrak{t}+\mathfrak{h}\zeta )W_i(\mathfrak{t}+\mathfrak{h}\zeta )}{1+q{\left(U_i(\mathfrak{t}+\mathfrak{h}\zeta )\right)}^2}}+{\displaystyle \frac{u_i(\mathfrak{t}+\mathfrak{h}\zeta )w_i(\mathfrak{t}+\mathfrak{h}\zeta )}{1+q{\left(u_i(\mathfrak{t}+\mathfrak{h}\zeta )\right)}^2}}2Q{e}_i^2(\mathfrak{t}+\mathfrak{h}\zeta )),\hfill \\ \hfill & =\sum _{i=1}^m{\displaystyle \frac{k_1}{{\mathrm{\Delta }}_{\mathrm{x}}^2}}{e}_{1i}(\mathfrak{t}+\mathfrak{h}\zeta ){\mathrm{\Delta }}^2{e}_{1,i-1}(\mathfrak{t}+\mathfrak{h}\zeta )+{\displaystyle \frac{k_2}{{\mathrm{\Delta }}_{\mathrm{x}}^2}}{e}_{2i}(\mathfrak{t}+\mathfrak{h}\zeta ){\mathrm{\Delta }}^2{e}_{2,i-1}(\mathfrak{t}+\mathfrak{h}\zeta )-2Q({\left(e_{1i}\right)}^2+{\left(e_{2i}\right)}^2)\hfill \\ \hfill & +\left({\displaystyle \frac{U_i(\mathfrak{t}+\mathfrak{h}\zeta )W_i(\mathfrak{t}+\mathfrak{h}\zeta )}{1+q{\left(U_i(\mathfrak{t}+\mathfrak{h}\zeta )\right)}^2}}+{\displaystyle \frac{u_i(\mathfrak{t}+\mathfrak{h}\zeta )w_i(\mathfrak{t}+\mathfrak{h}\zeta )}{1+q{\left(u_i(\mathfrak{t}+\mathfrak{h}\zeta )\right)}^2}}\right)(e_{1i}+e_{2i}),\hfill \\ \hfill & \le {\displaystyle \frac{k_1}{{\mathrm{\Delta }}_{\mathrm{x}}^2}}\sum _{i=1}^m{e}_{1i}(\mathfrak{t}+\mathfrak{h}\zeta )-\left(\mathrm{\Delta }\mathrm{\Delta }{e}_{1,i-1}(\mathfrak{t}+\mathfrak{h}\zeta )\right)+{\displaystyle \frac{k_2}{{\mathrm{\Delta }}_{\mathrm{x}}^2}}\sum _{i=1}^m{e}_{2i}(\mathfrak{t}+\mathfrak{h}\zeta )-\left(\mathrm{\Delta }\mathrm{\Delta }{e}_{2,i-1}(t+\mathfrak{h}\zeta )\right)\hfill \\ \hfill & -2Q\sum _{i=1}^m({\left(e_{1i}\right)}^2+{\left(e_{2i}\right)}^2)+\sum _{i=1}^m|{\displaystyle \frac{U_i(\mathfrak{t}+\mathfrak{h}\zeta )W_i(\mathfrak{t}+\mathfrak{h}\zeta )}{1+q{\left(U_i(\mathfrak{t}+\mathfrak{h}\zeta )\right)}^2}}+{\displaystyle \frac{u_i(\mathfrak{t}+\mathfrak{h}\zeta )w_i(\mathfrak{t}+\mathfrak{h}\zeta )}{1+q{\left(u_i(\mathfrak{t}+\mathfrak{h}\zeta )\right)}^2}}|\left(\right|e_{1i}|+|e_{2i}\left|\right),\hfill \\ \hfill & \le {\displaystyle \frac{k_1}{{\mathrm{\Delta }}_{\mathrm{x}}^2}}\left(\mathrm{\Delta }{e}_{1,i-1}\mathrm{\Delta }{e}_{1,i-1}{(\mathfrak{t}+\mathfrak{h}\zeta )|}_{m+1}^1-\sum _{i=1}^m{\left(\mathrm{\Delta }{e}_{1,i-1}(\mathfrak{t}+\mathfrak{h}\zeta )\right)}^2\right)-2Q\sum _{i=1}^m({\left(e_{1i}\right)}^2+{\left(e_{2i}\right)}^2)\hfill \\ \hfill & +{\displaystyle \frac{k_2}{{\mathrm{\Delta }}_{\mathrm{x}}^2}}\left(\mathrm{\Delta }{e}_{2,i-1}\mathrm{\Delta }{e}_{2,i-1}{(\mathfrak{t}+\mathfrak{h}\zeta )|}_{m+1}^1-\sum _{i=1}^m{\left(\mathrm{\Delta }{e}_{2,i-1}(\mathfrak{t}+\mathfrak{h}\zeta )\right)}^2\right)+Q\sum _{i=1}^m\left(\right|e_{1i}|+|e_{2i}{\left|\right)}^2,\hfill \\ \hfill & \le -{\displaystyle \frac{k_1}{{\mathrm{\Delta }}_{\mathrm{x}}^2}}\sum _{i=1}^m{\left(\mathrm{\Delta }{e}_{1,i-1}(\mathfrak{t}+\mathfrak{h}\zeta )\right)}^2-{\displaystyle \frac{k_2}{{\mathrm{\Delta }}_{\mathrm{x}}^2}}\sum _{i=1}^m{\left(\mathrm{\Delta }{e}_{2,i-1}(\mathfrak{t}+\mathfrak{h}\zeta )\right)}^2-Q\sum _{i=1}^m\left(\right|e_{1i}|-|e_{2i}{\left|\right)}^2\le 0.\hfill \end{array}$$

According to the Lyapunov stability theory stated in Theorem 2, this implies the global asymptotic stability of the zero solution of the error system (24). Consequently, the master-slave systems (13) and (21) are completely synchronized. □

6. Numerical Simulation

To demonstrate and confirm the synchronization techniques proposed in the preceding section. We provide the following examples with numerical simulations:

Example 1.

We consider the master-slave systems (13)–(15) with the following parameters: $(\mathrm{a},\mathrm{b},\sigma ,d_1,d_2)=(10,1,2,1,1.5),N=100,\mathfrak{h}=1.5,\mathfrak{t}\in [0,150],$$\mathrm{x}\in [0,20]$, the boundary conditions $(u_0\left(\mathfrak{t}\right),w_0\left(\mathfrak{t}\right))=(1,5),(u_1\left(\mathfrak{t}\right),w_1\left(\mathfrak{t}\right))=(1,5)$ and with the following initial conditions.

$$\left\{\begin{array}{c}({\psi }_1\left({\mathrm{x}}_i\right),{\psi }_2\left({\mathrm{x}}_i\right))=(7+0.3sin\left(5\pi {\mathrm{x}}_i\right),7+0.6sin\left(5\pi {\mathrm{x}}_i\right)),\hfill \\ ({\mathrm{\Phi }}_1\left({\mathrm{x}}_i\right),{\mathrm{\Phi }}_2\left({\mathrm{x}}_i\right))=(7+0.2cos\left(5\pi {\mathrm{x}}_i\right),7+0.2cos\left(5\pi {\mathrm{x}}_i\right)).\hfill \end{array}\right.$$

First, the assumptions given in Theorem 3 is satisfied for controlling the master-slave discrete fractional reaction-diffusion systems using the following linear controllers.

$$M=\left(\begin{array}{c}0.90\\ 03\end{array}\right).$$

As a result, systems (13) and (15) are completely synchronized. We provide the numerical solution of the system (13) in (30). Moreover, Figure 1 and Figure 2 show the solutions $u_i$, $w_i$, $U_i$ and $W_i$, also, Figure 3 illustrate the time development of error system states $e_{1i}$ and $e_{2i}$ in this case.

$$\left\{\begin{array}{c}{u}_i\left(\mathrm{n}\mathfrak{h}\right)={\psi }_1\left({\mathrm{x}}_i\right)+{\displaystyle \frac{{\mathfrak{h}}^{\zeta }}{\mathrm{\Gamma }\left(\zeta \right)}}\sum _{\mathfrak{p}=1}^{\mathrm{n}}{\displaystyle \frac{\mathrm{\Gamma }(\mathrm{n}-\mathfrak{p}+\zeta )}{\mathrm{\Gamma }(\mathrm{n}-\mathfrak{p}+1)}}\times [{\displaystyle \frac{u_{i+1}\left((\mathfrak{p}-1)\mathfrak{h}\right)-2u_i\left((\mathfrak{p}-1)\mathfrak{h}\right)+u_{i-1}\left((\mathfrak{p}-1)\mathfrak{h}\right)}{{\mathrm{\Delta }}_{\mathrm{x}}^2}}\hfill \\ +\begin{array}{c}\hfill \mathrm{a}\end{array}-u_i\left((\mathfrak{p}-1)\mathfrak{h}\right)-{\displaystyle \frac{4u_i\left((\mathfrak{p}-1)\mathfrak{h}\right)w_i\left((\mathfrak{p}-1)\mathfrak{h}\right)}{1+{\left(u_i\left((\mathfrak{p}-1)\mathfrak{h}\right)\right)}^2}}],\hfill \\ w_i\left(\mathrm{n}\mathfrak{h}\right)={\psi }_2\left({\mathrm{x}}_i\right)+{\displaystyle \frac{{\mathfrak{h}}^{\zeta }}{\mathrm{\Gamma }\left(\zeta \right)}}\sum _{\mathfrak{p}=1}^{\mathrm{n}}{\displaystyle \frac{\mathrm{\Gamma }(\mathrm{n}-\mathfrak{p}+\zeta )}{\mathrm{\Gamma }(\mathrm{n}-\mathfrak{p}+1)}}\times [{\displaystyle \frac{w_{i+1}\left((\mathfrak{p}-1)\mathfrak{h}\right)-2w_i\left((\mathfrak{p}-1)\mathfrak{h}\right)+w_{i-1}\left((\mathfrak{p}-1)\mathfrak{h}\right)}{{\mathrm{\Delta }}_{\mathrm{x}}^2}}\hfill \\ +\sigma \mathrm{b}\left(u_i\left((\mathfrak{p}-1)\mathfrak{h}\right)-{\displaystyle \frac{u_i\left((\mathfrak{p}-1)\mathfrak{h}\right)w_i\left((\mathfrak{p}-1)\mathfrak{h}\right)}{1+{\left(u_i\left((\mathfrak{p}-1)\mathfrak{h}\right)\right)}^2}}\right)],\hfill \\ 1\le i\le m,\hfill \\ \mathrm{n}>0.\hfill \end{array}\right.$$

(30)

Example 2.

To keep track of the system’s (13) performance, we alter the system’s parameters and order, taking into account parameter values: $(\mathrm{a},\mathrm{b},q,k_1,k_2,\zeta )=(1.2371,0.1,9,3,2,0.35)$ and

$$\left\{\begin{array}{c}{\psi }_1\left({\mathrm{x}}_i\right)=0.2(3+0.1cos\left(0.5{\mathrm{x}}_i\right)),\hfill \\ {\psi }_2\left({\mathrm{x}}_i\right)=0.2(4+0.3sin\left(0.2{\mathrm{x}}_i\right)).\hfill \end{array}\right.$$

(31)

Additionally, we set

$$\left\{\begin{array}{c}{\mathrm{\Phi }}_1\left({\mathrm{x}}_i\right)=0.7sin\left(0.3{\mathrm{x}}_i\right)),\hfill \\ {\mathrm{\Phi }}_2\left({\mathrm{x}}_i\right)=0.5cos\left(0.3{\mathrm{x}}_i\right)).\hfill \end{array}\right.$$

(32)

With the periodic conditions $(u_0\left(\mathfrak{t}\right),w_0\left(\mathfrak{t}\right))=(3,1),(u_1\left(\mathfrak{t}\right),w_1\left(\mathfrak{t}\right))=(3,1)$ and $(U_0\left(\mathfrak{t}\right),W_0\left(\mathfrak{t}\right))=(4,2),(U_1\left(\mathfrak{t}\right),W_1\left(\mathfrak{t}\right))=(4,2)$.

As a consequence of the numerical simulations, we can see that by adding appropriate controllers as shown in (27), the dynamics of (13) and (21) are synchronized, and the zero constant state of the synchronization error system expressed in (28) is asymptotically stable. Figure 4 and Figure 5 are numerical simulations of the master-slave systems under the considered parameters. Furthermore, Figure 6 indicates that the system’s zero steady-state is asymptotically stable, moreover, Figure 7 shows the same results in the $2D$ spatial domain.

7. Conclusions

In this study, we present a unique version of the Degn–Harrison reaction-diffusion systems and the Lengyel–Epstein reaction-diffusion systems that depend on the Caputo h-difference operator. We developed unique approaches for investigating synchronization in a spatiotemporal model of nonlinear bacterial colonies. First, for complete synchronization, suitable control schemes for synchronization are presented. The results of synchronization are based on Lyapunov theory and the master-slave formulation. To demonstrate the efficacy as well as the validity of the suggested synchronization schemas, numerical simulations of discrete time-fractional order Degn–Harrison systems and Lengyel–Epstein systems are provided. In the future, our plan is to further investigate bacterial colonies and the related reaction-diffusion synchronization phenomena, with the aim of developing sensor-based applications.