# Controlling Fractional Difference Equations Using Feedback.

Divya D. Joshi
Department of Physics  Rashtrasant Tukadoji Maharaj Nagpur University
Amravati Rd  Nagpur- 440033  Maharashtra  India  Sachin Bhalekar
School of Mathematics and Statics  University of Hyderabad  Lingampally CR Rao
Road  Hyderabad Central University Rd  Gachibowli  Hyderabad-500046  Telangana
India  Prashant M. Gade
Department of Physics  Rashtrasant Tukadoji Maharaj Nagpur University
Amravati Rd  Nagpur- 440033  Maharashtra  India

###### Abstract

One of the most popular methods of controlling dynamical systems is feedback.
It can be used without acquiring detailed knowledge of the underlying system.
In this work, we study the stability of fractional-order linear difference
equations under feedback. The stability results are derived for an arbitrary
feedback time $\tau$. We study the cases of $\tau=1$ and $\tau=2$ in further
detail. The extension to the stability of fixed points under feedback for
nonlinear fractional order difference equations with fixed points $x_{*}=0$ is
also carried out.

##  1 Introduction

Differential Equations have been used for modeling various phenomena in
natural sciences for a long time. This modeling helps us in understanding the
physical phenomena and controlling them if necessary. The differential
equations have been found useful in modeling plethora of systems. They range
from the spreading of diseases [42], emotional self-regulation in romantic
couples [44], economic development and growth [55] to tumor growth [18]. This
modeling is inadequate in certain systems and generalization is required.
Fractional differential equations are the generalization of differential
equations for systems with memory. The order of fractional differential
equations could be real or even complex. Even though fractional calculus has
been around for more than 300 years, fractional differential equations are
applied in real world situations only in the past few decades. We find the
first mention of fractional derivatives by Leibniz and L’Hospital as early as
1695 but the very first definition of fractional derivative was introduced by
Liouville and Riemann in the second half of the 19th century. Thereafter
several eminent mathematicians like Caputo, Hadamard, Grunwald, Letnikov,
Riesz, and others gave various definitions for fractional derivatives. This
development in the field of calculus opened the door for many other fields
where fractional differential equations are an essential part of mathematical
modeling.

The systems defined by fractional differential equations are said to be non-
local. The reason is that the future depends on its entire history. Thus,
systems governed by fractional differential equations have long-term memory.
Hence, fractional differential equations are ideal models for systems where
memory plays an important role. Such systems are found in diverse fields. In
the field of material science, fractional differential equations are used to
define viscoelastic materials and their order determines the amount of
viscosity and elasticity present in the materials [39]. Fractional order
epidemic models of various diseases like COVID-19 [56], Ebola [17], HIV [30],
influenza A (H1N1) [20] have shown promising results that helped to
understand, analyze and control the spread of these diseases. Seismological
studies have shown that fractional order intensity measures for probabilistic
seismic demand modeling applied to highway bridges have superior performance
compared to traditionally used intensity measures. This improved efficiency
and proficiency, at the same time, maintaining practicality and sufficiency
[40] for seismological applications. Several books have been devoted to
applications of fractional calculus in real-life applications [46] and recent
applications of fractional calculus in the field of science and engineering
can be found in [45].

Fractional order systems are different from integer order systems in several
respects. However, they show phenomena observed in integer order systems such
as chaotic or aperiodic behavior. In certain systems, chaos is not desired and
several control schemes are designed to control the chaos. Two of the most
popular control schemes are the Ott-Grebogi-Yorke method [34] and the Pyragas
method [37]. For the first method, we need to know the stable and unstable
manifold of the desired orbit. For fractional order systems, the presence of
such manifolds itself is a matter of debate [10]. The method that can be used
without any detailed knowledge of the systems is the feedback method suggested
by Pyragas. We study the possibility of controlling fixed points using
feedback and show that the method indeed works in the fractional case as well.

In the case of systems with delay $\tau$, the value of the current state
depends on its value $\tau$ steps back. Such systems have applications both in
modeling as well as control. Fractional differential equations with delay have
found applications in control theory. Control theory deals with the control of
dynamical systems, uses feedback in several cases and has wide applications.
Controllers may be designed using a feedback mechanism. Here, the output of
the system is fed back to the system through a controller to influence the
behavior of the system and give the desired output.

Stability analysis of delay differential equation with control is carried out
by Bhalekar [7]. Controllers have been designed for controlling the fractional
ordered systems with input delay [41, 53, 54, 31]. The stability of the
Cournot duopoly model with distributed time delays is discussed by Culda et
al. in [13]. The existence of chaotic behavior, stability, and its
synchronization of fractional ordered systems with delay have been studied for
several systems like the Ikeda system [25], logistic system [47] and Chen
system [14]. The ongoing scenario in the world demands newer and updated
knowledge of the spreading of diseases and their controllability. As discussed
earlier, fractional-order differential equations are used for modeling
epidemiological problems. Certain realistic situations demand the introduction
of a delay in resulting equations and could lead to better modeling and
prediction. A rich dynamical behavior is obtained for the infection model of
fractional order with delay [26].

The fractional differential equations with delay have been studied in various
contexts. Rihan et al. investigated the fractional order model of the
interactions of tumor cells and the immune systems with two different time
delays. The stability of the solutions was observed to have improved and the
model leads to various complex behaviors [38]. Alzahrani et al. studied the
adverse effect of untimely or delayed reporting of infectious diseases more
realistically by using fractional differential equations with delay [2].

Fractional difference equations are a relatively much less studied models.
Difference equations are the discretized versions of differential equations.
The study of fractional difference equations can be seen as an approach to
studying fractional differential equations by the finite difference method.
Furthermore, if an integer-order difference equation is generalized to a
fractional-order one, then it is able to model the memory properties in the
system. This is due to the non-local property of the fractional order
difference operator. In this generalization, all the values of the system from
an initial point are considered while evaluating the new value. These
equations demand fewer computational resources and are simpler to code. The
stability conditions for difference equations of fractional order and even
complex order have already been obtained by [43, 9, 24]. Stability results of
two-term fractional difference equations are proposed by Brandibur and Kaslik
[11]. Stability conditions are crucial for studying the control and
synchronization of the systems. For systems of complex order, we have
numerically investigated a fractional difference equation along with a delay
term which can be viewed as a controller [24]. We observed that the parameter
range over which chaos is obtained is reduced on the introduction of delay
term for the complex fractional Lozi map. This can be viewed as a control. One
of the simplest control is a case of a fixed point where the system gives
steady output. Though the above work indicates that feedback can be useful for
control in fractional order maps, we need rigorous analytic conditions for
practical applications.

The theory of fractional finite differences is initiated by Lubich [27] and
Miller and Ross [28]. The topic is further developed by Atici and coworkers
[4]. The stability analysis of these equations is presented in [12, 11]. The
chaotic systems of these kinds are employed in the image processing by
Abdeljawad et al. [1]. Atici and Sengul introduced the fractional calculus of
variations and derived the discrete Euler-Lagrange equation [5]. A tumor-
growth in cancer is modelled by using the nabla-fractional difference
equations in [3]. Ouannas, Batiha and Pham [35] proposed the theory and
applications of chaotic difference equations of fractional order. These
systems are recently used to model COVID-19 [16]. Fractional order Mandelbrot
set and Julia sets are studied in [15].

This paper gives the stability analysis for stable fixed points for systems
defined by fractional difference equations coupled with a delay term. The
exact analysis is carried out for the linear system and is extended to the
nonlinear maps. Stabilization or destabilization of fixed points of fractional
order maps with feedback is studied. Because this is a discrete system, the
control term can be added proportionally to the value at the previous time.
This is unlike the differential equations studied above where the control term
is within the integration and thus the control is also fractional order. This
is a simpler and more practical case where the control at a given time depends
on the value of variable $\tau$ steps back and not the entire history. We
obtain the stability bounds for this system. The analysis is essentially for
linear systems. However, we observe that the same analysis gives equally good
stability bounds for nonlinear maps with appropriate linearization near the
zero fixed point. In this work, we have considered the $h-$ difference
operator while defining the system. For simplicity, we have taken $h=1$.
Therefore, $\tau$ has integer values. Fractional values of $\tau$ can be
considered for $h\neq 1$.

The plan of the paper is as follows. We give essential definitions followed by
the model. We carry out stability analysis of the fractional order difference
equations with delay. We show that the stability conditions can be expressed
in an equivalent matrix form. We study the cases $\tau=1$ and $\tau=2$ in
greater detail. The case of large $\tau$ is also discussed briefly. Several
examples are given for linear and nonlinear systems to corroborate our results
with numerical evidence in various systems.

##  2 Preliminaries

In this section, we have given some basic definitions. Let $h>0$,
$a\in\mathbb{R}$, $(h\mathbb{N})_{a}=\\{a,a+h,a+2h,\ldots\\}$ and
$\mathbb{N_{\circ}}=\\{0,1,2,\ldots\\}$.

###### Definition 2.1

(see [29]) The Z-transform of a sequence $\\{y(n)\\}_{n=0}^{\infty}$ is a
complex function given by

$Y(z)=Z[y](z)=\sum_{k=0}^{\infty}y(k)z^{-k}$

where $z\in\mathbb{C}$ is a complex number for which the series converges
absolutely.

###### Definition 2.2

(see [19, 6]) Let $h>0,\;a\in\mathbb{R}$ and
$(h\mathbb{N})_{a}=\\{a,a+h,a+2h,\ldots\\}$. For a function
$x:(h\mathbb{N})_{a}\rightarrow\mathbb{C}$, the forward h-difference operator
is defined as

$(\Delta_{h}x)(t)=\frac{x(t+h)-x(t)}{h},$

where $t$ $\in(h\mathbb{N})_{a}$.

Throughout this article, we take $a=0$ and $h=1$. We write $\Delta$ for
$\Delta_{1}$. Now, we generalize the fractional order operators defined in
[[29, 19, 6]].

###### Definition 2.3

For a function $x:(h\mathbb{N})_{a}\rightarrow\mathbb{C}$ the fractional h-sum
of order $\alpha=u+\iota v\in\mathbb{C},u>0$ is given by

$(_{a}\Delta_{h}^{-\alpha}x)(t)=\frac{h^{\alpha}}{\Gamma(\alpha)}\sum_{s=0}^{n}\frac{\Gamma(\alpha+n-s)}{\Gamma(n-s+1)}x(a+sh),\\\
$

where, $t=a+(\alpha+n)h,\;n\in\mathbb{N_{\circ}}$.

For $h=1$ and $a=0$ , we have

$\displaystyle(\Delta^{-\alpha}x)(t)$ $\displaystyle=$
$\displaystyle\frac{1}{\Gamma(\alpha)}\sum_{s=0}^{n}\frac{\Gamma(\alpha+n-s)}{\Gamma(n-s+1)}x(s)$
$\displaystyle=$
$\displaystyle\sum_{s=0}^{n}\left(\begin{array}[]{c}n-s+\alpha-1\\\ n-s\\\
\end{array}\right)x(s).$

Here, we used the generalized binomial coefficient

$\left(\begin{array}[]{c}\mu\\\ \eta\\\
\end{array}\right)=\frac{\Gamma(\mu+1)}{\Gamma(\eta+1)\Gamma(\mu-\eta+1)},\\\
\;\mu,\eta\in\mathbb{C},\;\text{Re}(\mu)>0,\;\text{and Re}(\eta)>0.$

If $n$ $\in$ $\mathbb{N_{\circ}}$ then

$\displaystyle\left(\begin{array}[]{c}\mu\\\
n\end{array}\right)=\frac{(\mu+1)}{n!\Gamma(\mu-n+1)}=\frac{\mu(\mu-1)\ldots(\mu-n-1)}{n!}.$

###### Definition 2.4

For $n\in\mathbb{N_{\circ}}$ and $\alpha=u+\iota v\in\mathbb{C},u>0,$ we
define

$\displaystyle\tilde{\phi}_{\alpha}(n)=\left(\begin{array}[]{c}n+\alpha-1\\\
n\\\ \end{array}\right)=(-1)^{n}\left(\begin{array}[]{c}-\alpha\\\
n\end{array}\right).$

Note: The convolution $\tilde{\phi}_{\alpha}*x$ of the sequences
$\tilde{\phi}_{\alpha}$ and $x$ is defined as

$\left(\tilde{\phi}_{\alpha}*x\right)(n)=\sum_{s=0}^{n}\tilde{\phi}_{\alpha}(n-s)x(s)$
$\therefore(\Delta^{-\alpha}x)(n)=(\tilde{\phi}_{\alpha}*x)(n)\\\ .$
$\displaystyle\therefore
Z(\Delta^{-\alpha}x)(n)=Z\left(\tilde{\phi}(n)\right)Z(x(n))$
$\displaystyle=(1-z^{-1})^{-\alpha}X(z),$

where $X$ is $Z$ transform of $x$.

###### Property 2.1

(see [36]) The time-shifting property shows how a change in the discrete
function’s time domain alters the Z-domain.

$Z[x(k-n)]=z^{-n}X(z)$

Proof: From Definition 2.1 we have,

$X(z)=\sum_{k=0}^{\infty}x(k)z^{-k}.$

Consider $k-n=m$ i.e., $k=m+n$. Thus, we write the z-transform equation as

$Z[x(k-n)]=\sum_{k=0}^{\infty}x(k-n)z^{-k}\\\
=\sum_{m=0}^{\infty}x(m)z^{-(m+n)}\\\ =\sum_{m=0}^{\infty}x(m)z^{-m}z^{-n}\\\
=z^{-n}\sum_{m=0}^{\infty}x(m)z^{-m}\\\ =z^{-n}X(z)$

###### Lemma 2.1

For $\alpha\in\mathbb{C},\;\text{Re}(\alpha)>0$,

$Z(\tilde{\phi}_{\alpha}(t))=\frac{1}{(1-z^{-1})^{\alpha}}.$

Proof: We have,

$\displaystyle Z(\tilde{\phi}_{\alpha}(t))$ $\displaystyle=$
$\displaystyle\sum_{j=0}^{\infty}\tilde{\phi}_{\alpha}(j)z^{-j}$
$\displaystyle=$
$\displaystyle\sum_{j=0}^{\infty}\left(\begin{array}[]{c}j+\alpha-1\\\
j\end{array}\right)z^{-j}$ $\displaystyle=$
$\displaystyle\sum_{j=0}^{\infty}(-1)^{j}\left(\begin{array}[]{c}-\alpha\\\
j\end{array}\right)z^{-j}$ $\displaystyle=$
$\displaystyle(1-z^{-1})^{-\alpha}.$

by using Newton’s generalization of the Binomial Theorem. [32, 33].

##  3 Model

Consider the fractional order linear difference equation

$x(t)=x_{0}+\sum_{j=0}^{t-1}\frac{\Gamma(t-j+\alpha-1)}{\Gamma(\alpha)\Gamma(t-j)}(a-1)x(j).$
(6)

In this paper, we study the stability analysis of a fractional difference
equation (6) coupled with a delay term.

### Modeling Equation

Introducing the delay term in the equation as $bx(t-\tau)$, we get

$x(t)=bx(t-\tau)+x_{0}+\sum_{j=0}^{t-1}\frac{\Gamma(t-j+\alpha-1)}{\Gamma(\alpha)\Gamma(t-j)}\left((a-1)x(j)\right),$
(7)

where, $b\in\mathbb{R}$ and $a\in\mathbb{C}$.

$\therefore x(t+1)-bx(t+1-\tau)=x_{0}+(a-1)(\tilde{\phi}_{\alpha}*x)(t).$

Taking Z-transform on both sides, we get

$\displaystyle zX(z)-zx_{0}-bz^{(1-\tau)}X(z)$ $\displaystyle=$
$\displaystyle\frac{x_{0}}{1-z^{-1}}+\frac{(a-1)}{(1-z^{-1})^{\alpha}}X(z).$
$\displaystyle\therefore
X(z)\left(z-bz^{(1-\tau)}-\frac{(a-1)}{(1-z^{-1})^{\alpha}}\right)$
$\displaystyle=$ $\displaystyle zx_{0}+\frac{x_{0}}{1-z^{-1}},$ (8)

where $|z|<1$.

### Characteristic Equation

From (8), the characteristic equation of (7) is

$\left(z(1-z^{-1})^{\alpha}-b(1-z^{-1})^{\alpha}z^{(1-\tau)}-(a-1)\right)=0,$
(9)

where the condition $|z|<1$ should be satisfied. Putting $z=e^{(it)},$ we get,

$\displaystyle
e^{(it)}(1-e^{-it})^{\alpha}-b(1-e^{-it})^{\alpha}e^{(1-\tau)it}-(a-1)$
$\displaystyle=$ $\displaystyle 0,$ $\displaystyle\text{i.e.
}e^{(it)}(1-e^{-it})^{\alpha}-b(1-e^{-it})^{\alpha}e^{(1-\tau)it}+1$
$\displaystyle=$ $\displaystyle a.$ (10)

### Matrix Representation

Equation (7) can be represented equivalently as the following system

$\displaystyle x(t)$ $\displaystyle=$ $\displaystyle
x_{0}+\sum_{j=0}^{t-1}\frac{\Gamma(t-j+\alpha-1)}{\Gamma(\alpha)\Gamma(t-j)}\left((a-1)x(j)\right)+by(t),y(t)=x(t-1)\text{
for }\tau=1.$ (11) $\displaystyle\therefore x(t+1)$ $\displaystyle=$
$\displaystyle x_{0}+(a-1)(\tilde{\phi}_{\alpha}*x)(t)+by(t+1).$

Taking z-transform, we get

$zX(z)-zx_{0}=\frac{x_{0}}{1-z^{-1}}+\frac{(a-1)}{(1-z^{-1})^{\alpha}}X(z)+b(zY(z)-zy_{0}),Y(z)=\frac{X(z)}{z}+x(-1).$
$\therefore
z(1-z^{-1})^{\alpha}X(z)-(a-1)X(z)-bz(1-z^{-1})^{\alpha}Y(z)=x_{0}(1-z^{-1})^{\alpha-1}+z(1-z^{-1})^{\alpha}x_{0}-bzy_{0}(1-z^{-1})^{\alpha},$
$\displaystyle zY(z)-X(z)=zx(-1).$
$\displaystyle\therefore\begin{bmatrix}z(1-z^{-1})^{\alpha}-(a-1)&-bz(1-z^{-1})^{\alpha}\\\
-1&z\end{bmatrix}\begin{bmatrix}X(z)\\\ Y(z)\end{bmatrix}=0.$
$\displaystyle\therefore\begin{vmatrix}z(1-z^{-1})^{\alpha}-(a-1)&-bz(1-z^{-1})^{\alpha}\\\
-1&z\end{vmatrix}=0.$ $\displaystyle\text{i.e.
}z(1-z^{-1})^{\alpha}-b(1-z^{-1})^{\alpha}+1=a.$

This is the equation for the model with $\tau=1$. Similarly, for $\tau=2$ we
get,

$\displaystyle\begin{vmatrix}z(1-z^{-1})^{\alpha}-(a-1)&0&-bz(1-z^{-1})^{\alpha}\\\
-1&(z-1)+1&0\\\ 0&-1&(z-1)+1\end{vmatrix}=0.$

We can generalize and write $\tau+1$ dimensional determinant for delay $\tau$
as follows.

$\displaystyle\begin{vmatrix}z(1-z^{-1})^{\alpha}-(a-1)&0&0&$\ldots$&-bz(1-z^{-1})^{\alpha}\\\
-1&(z-1)+1&0&\ldots&0\\\ 0&-1&(z-1)+1&\ldots&0\\\
\vdots&\ddots&\ddots&\ldots&\vdots\\\ 0&\ldots&-1&(z-1)+1&0\\\
0&0&\ldots&-1&(z-1)+1\end{vmatrix}=0.$

### Boundary curve

Separating real and imaginary parts of equation (10), we get

$Re(a)=2^{\alpha}\left(\sin\left(\frac{t}{2}\right)\right)^{\alpha}\left(\cos\left(\frac{\alpha\pi}{2}+t\left(1-\frac{\alpha}{2}\right)\right)-b\cos\left(\frac{\alpha\pi}{2}+t\left(1-\tau-\frac{\alpha}{2}\right)\right)\right)+1,$
$Im(a)=2^{\alpha}\left(\sin\left(\frac{t}{2}\right)\right)^{\alpha}\left(\sin\left(\frac{\alpha\pi}{2}+t\left(1-\frac{\alpha}{2}\right)\right)-b\sin\left(\frac{\alpha\pi}{2}+t\left(1-\tau-\frac{\alpha}{2}\right)\right)\right).$

Therefore, the parametric representation of the boundary condition is

$\displaystyle\beta(t)=$
$\displaystyle\\{2^{\alpha}\left(\sin\left(\frac{t}{2}\right)\right)^{\alpha}\left(\cos\left(\frac{\alpha\pi}{2}+t\left(1-\frac{\alpha}{2}\right)\right)-b\cos\left(\frac{\alpha\pi}{2}+t\left(1-\tau-\frac{\alpha}{2}\right)\right)\right)+1,$
(12) $\displaystyle
2^{\alpha}\left(\sin\left(\frac{t}{2}\right)\right)^{\alpha}\left(\sin\left(\frac{\alpha\pi}{2}+t\left(1-\frac{\alpha}{2}\right)\right)-b\sin\left(\frac{\alpha\pi}{2}+t\left(1-\tau-\frac{\alpha}{2}\right)\right)\right)\\},$

for $t\in[0,2\pi]$ in the complex plane. If the complex number $a$ lies inside
this anticlockwise oriented simple closed curve $\beta(t)$ then the system (7)
will be asymptotically stable.

###  3.1 Stability Analysis with $\tau=1$

Figure 1: bifurcation regions in $b\alpha$-plane for $\tau=1$ with
representative stability diagrams. Here we have magnified the range on the
y-axis for representation purposes. (NOT TO SCALE)

###### Theorem 3.1

Consider the delayed fractional order equation (7) with $\tau=1$. Let
$b=g_{1}(\alpha),b=g_{2}(\alpha)$ and $b=g_{4}(\alpha)$ be the bifurcation
curves as shown in Figure 1 in $b\alpha$-plane which are the branches of the
implicit curve

$\displaystyle g(b,\alpha)$ $\displaystyle=$ $\displaystyle
b\alpha\cos\left(\frac{1}{2}(1+\alpha)\left(\pi-\arctan\left(\frac{R}{S}\right)\right)\right)$
(13)
$\displaystyle-(-1+\alpha)\cos\left(\frac{1}{2}\left(\pi(1+\alpha)-(-1+\alpha)\arctan\left(\frac{R}{S}\right)\right)\right)$
$\displaystyle+\sin\left(\frac{1}{2}\left(\pi\alpha-(-3+\alpha)\arctan\left(\frac{R}{S}\right)\right)\right),$

where,

$R=\frac{-1+\alpha+b\alpha-b\alpha^{2}-\sqrt{(1+b\alpha)^{2}(1-2\alpha+4b\alpha+\alpha^{2})}}{b\alpha},$

and

$\displaystyle S$ $\displaystyle=$
$\displaystyle-\frac{\sqrt{2}}{b\alpha}(-(1+2(-1+b)b^{2}\alpha^{3}+b^{2}\alpha^{4}+\alpha^{2}-3b^{2}\alpha^{2}-2\alpha+2b\alpha$
$\displaystyle+(1-\alpha-b\alpha+b\alpha^{2})(\sqrt{(1+b\alpha^{2})(1-2\alpha+4b\alpha+\alpha^{2})})))^{\frac{1}{2}}.$

Furthermore, define the curve $b=g_{3}(\alpha)$ as the straight line $b=-1$ in
the $b\alpha$-plane.

We have the following stability results:
1\. If $b\in(-\infty,g_{4}(\alpha))\cup(g_{1}(\alpha),\infty)$, then the
system (7) is unstable.
2\. If $g_{4}(\alpha)<b<g_{3}(\alpha)$, then the boundary curve (12) produces
three disjoint regions; two of which are stable and one is unstable.
3\. If $g_{3}(\alpha)<b<g_{2}(\alpha)$, then the boundary curve (12) generates
a bounded region which is a stable region for the system.
4\. If $g_{2}(\alpha)<b<g_{1}(\alpha)$, then the single stable region in case
3 gets divided into two regions; one is stable and the other is unstable.
5\. If $b=g_{2}(\alpha)$, then the curve (12) is a smooth curve and the inside
part is stable.
6\. If $b=g_{3}(\alpha)$, then $\beta(0)=\beta(\pi)$ and the stable region
gets divided into two parts.
7\. If $b=g_{1}(\alpha)$ or $b=g_{4}(\alpha)$, then the curve $\beta$ will
have cusps.
Proof: The boundary curve $\beta(t)$ defined by (12) is closed because the
initial point $\beta(0)$ has the same value as that of final point
$\beta(2\pi)$. If $\beta(t)$ is simple and its orientation is anticlockwise
(positive) then the bounded part inside this curve is the stable region for
the system (7). As the parameter values change, the simple closed curve
$\beta(t)$ transforms to the non-simple curve i.e. a curve containing multiple
points. Therefore, the region bounded by the branch of the non-simple curve
$\beta(t)$ with anticlockwise orientation is the required stable region for
the system (7).

It is observed that the formation of cusps in $\beta(t)$ is responsible for
the generation of multiple points (self-intersections). Thus, the parameter
values at which cusps are produced in $\beta(t)$ are the points of
bifurcation. At the cuspidal point $\beta(t_{0})$, the curve
$\beta(t)=(x(t),y(t))$ cannot have the derivative
$\beta^{\prime}(t)=\frac{y^{\prime}(t)}{x^{\prime}(t)}$. Here, $x(t)$ and
$y(t)$ are the coordinates of the curve $\beta(t)$. In this case,
$\lim_{t\rightarrow t_{0}}\beta^{\prime}(t)$ can take either $\frac{0}{0}$ or
$\frac{\infty}{\infty}$ form. It is observed that at initial point $t=0$, this
limit takes $\frac{\infty}{\infty}$ form. Therefore, the initial point
$\beta(0)$ is responsible for the bifurcations in some cases.

To check another possibility (i.e. $\frac{0}{0}$ form), we solve the equations
$x^{\prime}(t)=0$ and $y^{\prime}(t)=0$ simultaneously. We eliminate $t$
between these equations and get a relation in the parameters $b$ and $\alpha$.
This is carried out by squaring and adding the equations $x^{\prime}(t)=0$ and
$y^{\prime}(t)=0$. Now, we get the value of $t$ as functions of $b$ and
$\alpha$. Substituting this $t$ value in the equation $y^{\prime}(t)=0$, we
get the implicit expression (13).

This curve $g(b,\alpha)=0$ has 3 branches viz
$b=g_{1}(\alpha),b=g_{2}(\alpha)$ and $b=g_{4}(\alpha)$ as shown in Figure 1.
Note that, $b=g_{2}(\alpha)$ is the straight line $b=1$ in the
$b\alpha$-plane.

Furthermore, the curve $\beta(t)$ intersects itself at $\beta(0)$ and
$\beta(\pi)$ as the parameter $b$ passes through the value $-1$. At this
intersection point $\beta(0)=\beta(\pi)$, $\beta^{\prime}(t)$ does not exist
at $t=0,\pi$ and we have a bifurcation apart from cusps. Thus, $b=-1$ is also
a bifurcation curve in the $b-\alpha$ plane. We denote it by
$b=g_{3}(\alpha)$.

If the $b$ value is above the curve $b=g_{1}(\alpha)$, the orientation of the
boundary curve $\beta(t)$ becomes clockwise and hence the system becomes
unstable for the values of ’$a$’ in the bounded as well as unbounded regions
described by $\beta(t)$.

As $b$ passes through the curve $b=g_{1}(\alpha)$, the $\beta(t)$ becomes
multiple curves and generates stable and unstable bounded regions as shown in
Figure 1. If we keep on decreasing the value of $b$ the unstable part becomes
smaller and vanishes at $b=1$. Note that $\beta$ is smooth at $t=0$ for
$g_{1}(\alpha)<b\leq g_{2}(\alpha)$. If $g_{2}(\alpha)<b<g_{3}(\alpha)$ then
the $\beta(t)$ is a simple closed curve with an anticlockwise orientation.
Therefore, the region bounded by $\beta(t)$ is the stable region and the
unbounded part of the complex plane is the unstable region.

As we take the values of $b$ between the curves $b=g_{3}(\alpha)$ and
$b=g_{4}(\alpha)$, $\beta(t)$ produces two stable and one unstable bounded
regions as shown in the Figure 1. The size of the unstable region increases if
we decrease the value of $b$.

At the bifurcation curve $b=g_{4}(\alpha)$, the stable regions vanish
completely and the $\beta(t)$ becomes a simple closed curve with clockwise
orientation. Thus, the system is completely unstable for $b<g_{3}(\alpha)$.
This proves the result.

###  3.2 Illustrative Examples

In this section, we provide corroborating evidence for stability results
obtained in the above section for $\tau=1$ and various values of $b$. We study
various regions shown in Figure 1.

###### Example 3.1

Let us begin with the region that lies above $g_{1}(\alpha)$. We take $b=18.3$
and $\alpha=0.1$ (Figure 2(a)). Here, we observe a non-smooth curve $\beta(t)$
with cusps on both ends. The number $a=-2$ lies inside the unstable region
according to Figure 2(b) and $a=(0.4+0.4\iota)$ which lies outside the curve,
is unstable as well (Figure 2(c)). The system is unstable for any value of
parameter $a$.

(a) Stability curve with $\tau=1$, $\alpha=0.1$ and $b=18.3$

(b) Unstable solution for $a=-2$.

(c) Unstable solution for $a=(0.4+0.4\iota)$.

Figure 2: Example 1

###### Example 3.2

Now consider the region between $b=g_{1}(\alpha)$ and $b=g_{2}(\alpha)$.
Consider $b=2$ and $\alpha=0.5$. Here, the stability curve is divided into two
parts (see Figure 3(a)). Consider a system (7) with $a=-1.2$. It is inside the
stable region and we get a stable trajectory as seen in Figure 3(b) For
$a=0.2$, which lies in the unstable region we get an unstable solution (Figure
3(c)). The solutions for the values of $a$ lying outside the stability curve
were checked and found unstable.

(a) Stability curve with $\tau=1$, $\alpha=0.5$ and $b=2$

(b) Stable solution for $a=-1.2$.

(c) Unstable solution for $a=0.2$.

Figure 3: Example 2

###### Example 3.3

On the boundary $b=g_{2}(\alpha)$ which is described by $b=1$ for any
$\alpha$. Consider the stability curve for $\alpha=0.4$ and $b=1$(see Figure
4(a)). It generates a single region. The quantity $a=0.5$ inside the curve
leads to a stable solution. For $a=-0.7\iota$ outside the curve, an unstable
solution is observed.

(a) Stability curve with $\tau=1$, $\alpha=0.4$ and $b=1$

(b) Stable solution for $a=0.5$.

(c) Unstable solution for $a=-0.7\iota$.

Figure 4: Example 3

###### Example 3.4

For $b$ in the region between $b=g_{2}(\alpha)$ and $b=g_{3}(\alpha)$, a
single stable region (Figure 5(a)) is observed. For $\alpha=0.3$ and $b=0.1$,
the value $a=0.3$ is inside this region and gives a stable solution (Figure
5(b)). On the other hand, $a=-0.4$ lies outside the stable region and leads to
an unstable solution (Figure 5(c)).

(a) Stability curve with $\tau=1$, $\alpha=0.3$ and $b=0.1$

(b) Stable solution for $a=-0.3$.

(c) Unstable solution for $a=-0.4$.

Figure 5: Example 4

###### Example 3.5

Examples with values of $b$ on the boundary curves give limiting cases. For
example $b=g_{2}(\alpha)=1$ is the limiting value of $b$ at which a single
stable region is observed and for in the next region, we have one stable and
one unstable region inside the stability curve. ( See Example 3.3, 3.2).
Similarly for $b=g_{3}(\alpha)=-1$ the values for $t=0$ and $t=\pi$ touch each
on the real axis and we have two stable regions. For $b<g_{3}(\alpha)=-1$, we
have 3 regions one of which is unstable. Consider the stability curve for
$\alpha=0.7$ and $b=-1$ (see Figure 6(a)). Stable trajectories are observed
for $a=1.1+0.5\iota$ and $a=1.3-0.8\iota$ (Figure 6(b) and 6(c) respectively).
We checked the stability of the parameter values lying outside the curve and
verified that they are unstable.

(a) Stability curve with $\tau=1$, $\alpha=0.7$ and $b=-1$

(b) Stable solution for $a=1.1+0.5\iota$.

(c) Stable solution for $a=1.3-0.8\iota$.

Figure 6: Example 5

###### Example 3.6

Now we consider region lying between $b=g_{3}(\alpha)$ and $b=g_{4}(\alpha)$.
The stability curve for $\alpha=0.2$ and $b=-3$ is sketched as shown in Figure
7(a). The curve has three regions, two of which look identical due to symmetry
with respect to the real axis. The system with $a=(4+1.3\iota)$ has a stable
solution (see Figure 7(b)). (Similarly, $a=(4-1.3\iota)$ leads to a stable
trajectory.) Thus these two regions have stable solutions. On the other hand,
for $a$ in the third region (say, $a=2$), we observe unstable trajectories
(see Figure 7(c)). It was checked that the parameter values $a$ lying outside
the stability curve lead to unstable solutions.

(a) Stability curve with $\tau=1$, $\alpha=0.2$ and $b=-3$

(b) Stable solution for $a=4+1.3\iota$.

(c) Unstable solution for $a=2$.

Figure 7: Example 6

###### Example 3.7

The last region lies below the curve $b=g_{4}(\alpha)$. Consider $\alpha=0.5$
and $b=-2.2$, the curve $\beta$ (see Figure 8(a)) encloses a single region.
There are no stable values of $a$ inside or outside the region. Consider
$a=1.5$ and $a=2.8$. Both values lead to unstable trajectories (see Figure
8(b) and 8(c), respectively).

(a) Stability curve with $\tau=1$, $\alpha=0.5$ and $b=-2.2$

(b) Unstable solution for $a=1.5$.

(c) Unstable solution for $a=2.8$.

Figure 8: Example 7

##  4 Stability result for $\tau=2$

Figure 9: bifurcation regions in $b\alpha$-plane for $\tau=2$ with
representative stability diagrams

###### Theorem 4.1

Now consider the delayed fractional order equation (7) with $\tau=2$. The
dynamics is richer and leads to five boundary curves that demarcate
qualitatively different behaviors. Let $b=g_{j}(\alpha)$, $j=1,2,\cdots,5$ be
the bifurcation curves in the $b\alpha$-plane (cf. Figure 9) which are the
branches of the implicit curve.

$\displaystyle g(b,\alpha)=$
$\displaystyle\cos\left(\frac{1}{2}\left(\pi\alpha-(-3+\alpha)K\right)\right)-b(1+\alpha)\cos\left(\frac{1}{2}\left(\pi\alpha-(3+\alpha)K\right)\right)+b\sin\left(\frac{1}{2}(1+\alpha)\left(\pi-K\right)\right)$
(14)
$\displaystyle+(\alpha-1)\sin\left(\frac{1}{2}\left(\pi(1+\alpha)-(-1+\alpha)K\right)\right),$

where,
$K=\arccos\left(\frac{-2+2\alpha-\alpha^{2}+b^{2}(2+2\alpha+\alpha^{2})}{2(-1+\alpha+b^{2}(1+\alpha))}\right)$.
For any $\alpha\in(0,1)$,
$-1=g_{5}(\alpha)<g_{4}(\alpha)<g_{3}(\alpha)<g_{2}(\alpha)<g_{1}(\alpha)=1$.

The stability results for the system (7) with $\tau=2$ are as follows:
1\. If $b\in(-\infty,g_{5}(\alpha))\cup(g_{1}(\alpha),\infty)$, then the
system is unstable.
2\. If $g_{5}(\alpha)<b<g_{4}(\alpha)$, then the boundary curve $\beta(t)$
generates four regions; of which the outer three are unstable and the central
one is stable.
3\. If $g_{4}(\alpha)<b<g_{3}(\alpha)$, then the number of unstable regions
reduces to one while the number of the stable region remains the same, i.e.
one. The two unstable regions described in Case 2 merge with the stable region
and hence the size of the stable region, in this case, is larger than that of
in Case 2.
4\. If $g_{3}(\alpha)<b<g_{2}(\alpha)$, then the region produced by the
boundary curve is a single closed curve which is stable.
5\. If $g_{2}(\alpha)<b<g_{1}(\alpha)$, then the boundary curve produces three
regions, one stable and two unstable.
Now we consider the boundary cases:
1\. If $b=g_{5}(\alpha)=-1$, then the curve $\beta(t)$ has three regions and
all the regions are unstable.
2\. If $b=g_{4}(\alpha)$, then the curve $\beta(t)$ has two regions, a stable
region with two cusps on the right side and an unstable region.
3\. If $b=g_{3}(\alpha)$, then the curve $\beta(t)$ has a single stable region
with a cusp on the left side.
4\. If $b=g_{2}(\alpha)$, then the curve $\beta(t)$ has a single region that
is stable and has two cusps.
5\. If $b=g_{1}(\alpha)=1$, then the curve $\beta(t)$ has two regions and both
regions are unstable.

Proof:
As discussed in the proof of Theorem 3.1, as the parameter values change the
simple closed curve $\beta(t)=(x(t),y(t))$ possesses multiple points and
generates different bounded regions, in this case also. Among these, the
regions with the anticlockwise oriented boundary are the stable ones. The
multiple points are generated through the cusps.

To find the conditions for the $\frac{0}{0}$ form of $\lim_{t\rightarrow
t_{0}}\beta^{\prime}(t)$, we solve $x^{\prime}(t)=0$ and $y^{\prime}(t)=0$
simultaneously. Squaring and adding these equations, we get

$t=\arccos\left(\frac{2(\alpha-1)-\alpha^{2}+b^{2}(\alpha^{2}+2\alpha+2)}{2(b^{2}(1+\alpha)+\alpha-1)}\right)$

.

Substituting this value in the equation $x^{\prime}(t)=0$, we get the
expression for the existence of cusps in $\beta(t)$ as equation (14). This
implicit curve $g(b,\alpha)$ has 5 branches namely $b=g_{j}(\alpha)$,
$j=1,2,\cdots,5$ as shown in Figure 9. These bifurcation curves produce 6
different regions labeled as A, B, C, D, E, and F of $b-\alpha$ plane (see
Figure 9). We list our observations below:

If the point $(b,\alpha)$ belongs to the region A i.e., if $b>g_{1}(\alpha)$,
the boundary curve $\beta(t)$ is a simple closed curve with a clockwise
orientation. Thus, the regions inside as well as outside $\beta(t)$ are
unstable.

In the region B, we have $g_{2}(\alpha)<b<g_{1}(\alpha)$. Here, $\beta(t)$
becomes a multiple curve and has two unstable and one stable bounded region.
As $b$ decreases in B, the range of stable region goes on increasing.

If $g_{3}(\alpha)<b<g_{2}(\alpha)$ (i.e. region C), $\beta(t)$ becomes a
simple closed curve with anticlockwise orientation. Thus the inside part of
$\beta(t)$ is the stable region and the outside part is unstable. As the value
of $b$ goes on decreasing, the stable region gets stretched horizontally and a
cusp is formed for $b=g_{3}(\alpha)$.

As $b$ is further reduced i.e., when $g_{4}(\alpha)<b<g_{3}(\alpha)$ (i.e.
region D), the cusp on the boundary curve $\beta(t)$ evolves into a loop-
shaped structure of clockwise orientation making $\beta$ a multiple curve with
one stable and one unstable region. As the value of $b$ moves closer towards
$g_{4}(\alpha)$, the area of unstable regions increases. At the bifurcation
curve $b=g_{4}(\alpha)$, two cusps are observed on the boundary curve $\beta$
on the right side.

In region E ($g_{5}(\alpha)<b<g_{4}(\alpha)$), both of the cusps on the right
of $\beta$ evolve into two unstable regions along with the previous unstable
region on the left side. Thus, the boundary curve $\beta(t)$ has 3 unstable
regions and one stable region. Reducing the value of $b$ further results in
the reduction of the size of the stable region.

When $b$ reaches the bifurcation curve $b=g_{5}(\alpha)$, the stable region
disappears completely, leaving the boundary curve $\beta(t)$ with 3 unstable
regions. These regions are unstable as a result of their clockwise
orientation.

In region F, i.e. for $b<g_{5}(\alpha)$ the boundary curve has four bounded
unstable regions with a clockwise orientation. Thus, below the bifurcation
curve $g_{5}(\alpha)$, the system becomes completely unstable for all values
of $a$.

As in the case of $\tau=1$, we have checked all the above results numerically.

(a)

(b)

Figure 10: Boundary curves are shown for the system (12) with $\alpha=0.2$,
$\tau=20,21,40,41$, (a) $b=-0.5$ and (b) $b=0.5$.

##  5 Asymptotic limit of large delay

We observe that we can stabilize a larger range of $a$ values with $\tau=1$
compared to $\tau=2$. The reason is that there are more routes to instability
for larger $\tau$. We now consider the system (15) for higher $\tau$ values.
Let us take $\tau=20,21,40$ and $41$ for $\alpha=0.2$. The boundary curves
$\beta$ have been plotted for the values of $b=-0.5$ and $b=0.5$ in the
figures 10 (a) and (b) respectively. In both figures, the stable region is in
the center and is surrounded by the spring-like structure. This spring-like
structure has a clockwise orientation and thus it is unstable. On the real
axis, the upper bound is one. The lower bound is indicated intersection of the
stability curve with the real axis. The lower bound of real stable $a$ is at
$\beta(\pi)$ when $\tau$ is even and $b$ is positive or $\tau$ is odd and $b$
is negative. However, as we increase $\tau$, the lower bound approaches
$\beta(\pi)$. For a given $\tau$, the stable region depends slightly on
whether $\tau$ is even or odd and the sign of $b$. This dependence reduces as
we increase $\tau$. The upper bound is fixed at $t=0$ for the boundary curve
$\beta(t)$, the lower bound approaches $\beta(\pi)$. We find that for large
$\tau$, the stable range of $a$ does not depend on the sign of $b$ and is
given by $a=1+2^{\alpha}(|b|-1)$.

The $b-a$ curves plotted in the figures 15 and 16 can be used to determine the
values of feedback coefficient $b$ to get the maximum stable range. As $\tau$
increases the $b-a$ curve reaches a limiting region. As mentioned above, the
stable $b-a$ region is given by $a=1$ and $a=1+2^{\alpha}(|b|-1)$ as
$\tau\rightarrow\infty$.

In a control system, the feedback term can be selected according to the
requirements of the systems. We see that the maximum range of asymptotically
stable fixed points can be obtained in the case of $\tau=1$. So, one can use
$\tau=1$ to maximize the range of asymptotically stable fixed points of the
system. Large delay can be used if the system requires enhanced chaos.

##  6 Efficacy of control for large negative multipliers

The range of stable real value of $a$ ranges from $1-2^{\alpha}$ to 1 for
system (6). The real range for integer order difference equation is between -1
and 1. There is a reduction in the real range for fractional order maps. Now
$a=-7$ is outside stable range for any $0<\alpha\leq 1$. The stability of
system (7) depends on the feedback coefficient ‘$b$’ and the delay ‘$\tau$’.
For $b=0$ , we have system (6) without feedback. Figure 11 shows stability
curves for system (7) with $\alpha=0.25,\tau=1,b=0$ (in black) and $b=6$ (in
red). The point $a=-7$ lies outside the black stability curve. It is enclosed
within small stable region of the red curve. We analyze the stability of this
point in both cases. We iterate the systems for $t=500$ time-steps. For $b=0$
case, the trajectory diverges and the solution is unstable (see Figure 12(a)).
For $b=6$ case, the trajectory asymptotically converges to zero and the
solution is stable(see Figure 12(b)). Thus large negative multipliers can be
stabilized with feedback for $\tau=1$.

Figure 11: Stability regions for system (7) $\alpha=0.25$.

(a)

(b)

Figure 12: Stability analysis for system (7) with $\alpha=0.25$, $\tau=1$ and
$a=-7$.

For the system (7) with $\tau=2$, the effect is less dramatic. Figure 13 shows
the stability curves for systems (7) with $\alpha=0.5,\tau=2,b=0$ (in black)
and $b=-0.6$ (in blue). In the first case $b=0$, the point $a=-1.1$ lies
outside the black stability curve and has an unstable solution (see Figure
14(a)). In second case $b=-0.6$, the point $a=-1.1$ lies inside the blue
stability curve and has a stable solution as seen in Figure 14(b).

Figure 13: Stability regions for system (7) $\alpha=0.5$.

(a)

(b)

Figure 14: Stability analysis for system (7) with $\alpha=0.5$, $\tau=2$ and
$a=-1.1$.

Thus, it is possible to stabilize a larger range of parameter values for a
fractional order system by adding a feedback term with delay (system (7)). The
case of $\tau=1$ gives the maximum range of stability with appropriate $b$ and
is recommended for control.

##  7 Nonlinear Maps

Consider a fractional order map with a feedback control

$x(t)=x_{0}+bx(t-\tau)+\frac{1}{\Gamma(\alpha)}\sum_{j=1}^{t}\frac{\Gamma(t-j+\alpha)}{\Gamma(t-j+1)}[f(x(j-1))-x(j-1)].$
(15)

###### Definition 7.1

A steady-state solution of (15) is called an equilibrium point. Thus, if
$x_{*}$ is an equilibrium point of (15) then $x_{0}=x_{*}$ implies
$x(t)=x_{*}$ for all $t=1,2,\cdots$.

Note that, the only equilibrium point of (15) is $x_{*}=0$ when we have
$f(0)=0$. Therefore, we consider the systems (15) with equilibrium point
$x_{*}=0$.

In a neighborhood of the point $x_{*}=0$, we have $f(x)\approx
f(x_{*})+(x-x_{*})f^{\prime}(x_{*})=ax$, where $a=f^{\prime}(0)$. Thus, the
local stability properties of the nonlinear system (15) at $x_{*}=0$ are the
same as those of the linear system (7).

We consider only real values of $a$ and $b$, in this case. Therefore, the
stable region given by the curve $\beta(t)$ defined by (12) is bounded by the
curves $a=1$, $a=-2^{\alpha}\left(1-b(-1)^{\tau}\right)$ and the parametric
curve
$a(t)=2^{\alpha}\left(\sin\left(\frac{t}{2}\right)\right)^{\alpha}\left(\cos\left(\frac{\alpha\pi}{2}+t\left(1-\frac{\alpha}{2}\right)\right)-\sin\left(\frac{\alpha\pi}{2}+t\left(1-\frac{\alpha}{2}\right)\right)\cot\left(\frac{\alpha\pi}{2}+t\left(1-\tau-\frac{\alpha}{2}\right)\right)\right)+1,$
$b(t)=\frac{\sin\left(\frac{\alpha\pi}{2}+t\left(1-\frac{\alpha}{2}\right)\right)}{\sin\left(\frac{\alpha\pi}{2}+t\left(1-\tau-\frac{\alpha}{2}\right)\right)},$
$t\in[0,2\pi]$ in the $b-a$ plane. We denote these curves as $b-a$ curves.

Thus, the system (15) can be controlled if the point $(b,a)$ with
$a=f^{\prime}(0)$ lies inside the region bounded by these curves. We consider
the nonlinear systems (15) of fractional order $\alpha$ and delay $\tau$. The
fractional logistic map in this model is defined by the equation (15) with
$f(x)=\lambda x(1-x)$, where $\lambda$ is a parameter and
$a=\lambda(1-2x_{*})=\lambda$.
Similarly, the fractional cubic map is (15) with $f(x)=\beta
x^{3}+(1-\beta)x$, where $a=3\beta x_{*}^{2}+(1-\beta)=1-\beta$.

Figure 15 shows the stability regions of all these systems for the $\tau=1$
for various values of $\alpha$.

We iterate these maps for $T$ time-steps where $T$ is large and the
equilibrium or asymptotic equilibrium point is assumed to be stable if the
convergence is obtained within $\delta$. We know the boundaries of the stable
fixed point in the $b-a$-plane a for linear system. For $a=f^{\prime}(0)$, we
observe that the $(b,a)$ values for stable fixed point lie in the analytically
obtained bounds for both $\tau=1$ and $\tau=2$. For $\tau=1$, we have plotted
$(b,a)$ values in Figure 15 for $\alpha=0.25,0.5,0.75$. It is observed for all
the stable zero fixed points of the logistic, and cubic maps lie within the
stable region defined by the $b-a$ curves.

We find that this $b-a$ region indeed encloses the stability region and gives
precise bounds. Thus if the stable equilibrium point is known, the required
strength of feedback for which the equilibrium point is stabilized can be
found. The uncontrolled case corresponds to $b=0$ and if we want to stabilize
(or destabilize) the equilibrium point an appropriate value of $b$ can be used
and desired range of stability can be obtained. We have studied
$\alpha=0.25,0.5,0.75$ and $\tau=1$. For $\tau=1$, the stability region
decreases in size with an increasing $\alpha$.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 15: Stable region for fractional systems with $\tau=1$. The stable
fixed points of fractional systems lie within the region enclosed by the $b-a$
curve. We discard $T$ time-steps and the convergence is within $\delta$. For
$\alpha=0.25$, $T=6\times 10^{4}$ and $\delta=10^{-9}$ for logistic map and
$T=4\times 10^{4}$ and $\delta-5^{-8}$ for cubic map. In all other cases
$T=2\times 10^{4}$ and $\delta=10^{-5}$.

Figure 16 confirms that the stable fixed point of the fractional maps lies
within the curve enclosed by the $\alpha$ dependent $b-a$ curves. Similarly,
for $\tau=2$, stability regions can be obtained for various values of
$\alpha$. The stable regions for the zero fixed point for logistic, and cubic
map with fractional order $\alpha=0.5$ and $\tau=2$ lie within the stable
region derived for linear case (see Figure 16).

(a)

(b)

Figure 16: Stable region for fractional systems with $\alpha=0.5$, $\tau=2$.
The stable fixed points of the fractional systems lie within the region
enclosed by the $b-a$ curve. We discard $T$ time-steps and the convergence is
within $\delta$. For logistic map $T=10^{5}$ and $\delta=10^{-9}$ and for
cubic map $T=5\times 10^{4}$ and $\delta=5\times 10^{-7}$.

Now we consider the case of large $\tau$. We consider the fractional logistic
map defined by the equation (15) with $f(x)=\lambda x(1-x)$, where
$f^{\prime}(0)=\lambda=a$. We iterate this map for $T=6\times 10^{4}$ and the
equilibrium point is assumed to be asymptotically stable if convergence is
obtained within $\delta=10^{-7}$. Here, we consider $\alpha=0.75$ and
$\tau=100$. The triangular stability region encloses the asymptotically stable
fixed points for the fractional logistic map. This is the generalized stable
region for large $\tau$. From the Figure 17 it is clear that the range of
stable region does not increase for any nonzero value of $b$. It decreases
with larger $|b|$. Thus, the large delay is not useful in stabilizing the
fixed point for the fractional system (15) and the best results are obtained
with $\tau=1$. However, as mentioned above, larger $\tau$ can be used to
destabilize the stable system.

Figure 17: Stability region in $b-a$ plane for the system (15) with
$\alpha=0.75$ and large $\tau$. The asymptotic stable fixed points of the
fractional logistic map defined by (15) with $\alpha=0.75$ and $\tau=100$ lie
inside the generalized stability region.

The control term in this model is a delay term with coefficient ‘$b$’. In
figure 18 (a) and (c), an increment in the range of nonzero stable fixed
points can be seen. We observe that the proper selection of control parameter
$b$ results in an extension of the range of the stable fixed points of the
system. The maximum range of the nonzero stable fixed points can be found with
the help of the $b-a$ curve for that particular $\alpha$ value the of
fractional system. The range of the stability of nonzero fixed points is
extended in presence of feedback. Preliminary numerical investigations
indicate that the above conditions are necessary but not sufficient for the
stability of nonzero fixed points. Further analytic investigations are needed
to find the stability of nonzero fixed points. Figures 18 (b) and (d) show the
time series of the system.

(a)

(b)

(c)

(d)

Figure 18: Figure (a) shows the bifurcation diagram for a logistic map with
$\tau=1$ and $b=1.1,\alpha=0.5$ where only a stable fixed point is observed.
The bifurcation diagram without control is also shown for reference. Figure
(c) shows bifurcation diagram for $\tau=2$ and $b=-0.58,\alpha=0.5$. The
bifurcation diagram without control is also shown for reference. The range
over which a stable fixed point is observed in either case. Figure (b) shows
the stabilization of the unstable chaotic state for $\tau=1$, $\alpha=0.5$,
and $\lambda=3.3$ to the stabilized fixed point. Similarly, figure (d) shows
the stabilization of the unstable chaotic state to stabilized fixed point for
$\tau=2$, $\alpha=0.5$, and $\lambda=3.23$.

##  8 Discussion and Conclusion

Fractional order systems with delay have been investigated in many contexts.
In neural networks, numerous works have been done. In [48], stability
conditions, and the existence of Hopf bifurcation for fractional order BAM
neural network (FOBAMNN) with delay have been established. In [22], FOBAMNN
with four delays has been turned into two delays and a correlation between
stability and delay term has been studied. Investigation of fractional order
neural network with multiple leakage delay shows that both fractional order
and time delay are very important in controlling the transient behaviors of
the FONN devised in [21]. The stability of fractional order triangle multi-
delayed neural networks has been studied in [49]. A comparative study of
integer order and fractional order delayed BAM neural network shows increased
stability region [51]. In [50], sufficient conditions for different delays for
ensuring stability and generation of Hopf bifurcation have been demonstrated.
Stability and bifurcation of an isovalent version of a fractional-order stage-
structured predator–prey system have been investigated in [52]. Global
asymptotic stabilization of fractional-order memristor-based neural networks
(FMNNs) with time delay can be achieved by adjusting two groups of parameters
as illustrated in [23]. The analysis of fractional order Bloch equation with
delay shows behaviors ranging from damped oscillations to oscillations with
increasing amplitude for various values of delay [8]. In this work, we carry
out basic investigations in the context of fractional order diffference
equations.

Control of chaos in dynamical systems is an important aspect of studies in the
theory of dynamical systems from viewpoint of applications. In integer order
differential and difference equations, it is a well-studied problem both
experimentally and theoretically. One of the simplest control schemes is
feedback and the Pyragas method is one such method with feedback delay control
to stabilize the chaotic systems. In this work, we studied systems defined by
fractional difference equations coupled with a delay term. The delay term acts
as a control in this system. We give analytic conditions for the stability of
the fixed point of these systems for the arbitrary delay. More detailed
analysis is carried out for $\tau=1$ and $\tau=2$ followed by analysis for
asymptotic limit. A detailed analysis is carried out for $\tau=1$ and
$\tau=2$. We give all bifurcation curves. They are given by $g(b,\alpha)$. For
real maps, we give the generalized stable regions for $\tau=1$, $\tau=2$ and
in the asymptotic limit. Using $g(b,\alpha)$, we obtain the upper and lower
bound of $b$ for particular $\alpha$ for which the stable regions exist. The
stability region of a system depends on the fractional order $\alpha$, delay
$\tau$, and the control parameter $b$. For nonlinear map $f(x)$, the zero
fixed point with slope ’$a=f^{\prime}(0)$’ in the same $b-a$ range. The
boundary curve (12) may enclose stable and/or unstable region/s for multiple
combinations of $\alpha$, $\tau$ and $b$. The criterion for finding out the
stable and unstable region is also given using the orientation of the curve.
Finally, we have studied nonlinear systems with fixed point $x_{*}=0$. Thus
the results apply to a broad range of systems and this is a practical scheme
for control. Our analysis can be used to stabilize/destabilize the fractional
order difference system by introducing the appropriate delay term. The chaotic
systems can be stabilized by selecting the delay $\tau=1$ whereas the stable
system can be made unstable/chaotic by selecting a larger value for $\tau$.

##  9 Acknowledgement

PMG and DDJ thank DST-SERB for financial assistance (Ref. CRG/2020/003993).

## References

  * [1] Thabet Abdeljawad, Santo Banerjee, and Guo-Cheng Wu.  Discrete tempered fractional calculus for new chaotic systems with short memory and image encryption.  Optik, 218:163698, 2020.
  * [2] Faris Alzahrani, Oyoon Abdul Razzaq, Daniyal Ur Rehman, Najeeb Alam Khan, Ali Saleh Alshomrani, and Malik Zaka Ullah.  Repercussions of unreported populace on disease dynamics and its optimal control through system of fractional order delay differential equations.  Chaos, Solitons & Fractals, 158:111997, 2022.
  * [3] Ferhan M Atıcı, Mustafa Atıcı, Ngoc Nguyen, Tilekbek Zhoroev, and Gilbert Koch.  A study on discrete and discrete fractional pharmacokinetics-pharmacodynamics models for tumor growth and anti-cancer effects.  Computational and Mathematical Biophysics, 7(1):10–24, 2019.
  * [4] Ferhan M Atici and Paul W Eloe.  A transform method in discrete fractional calculus.  International Journal of Difference Equations, 2(2), 2007.
  * [5] Ferhan M Atıcı and Sevgi Şengül.  Modeling with fractional difference equations.  Journal of Mathematical Analysis and Applications, 369(1):1–9, 2010\.
  * [6] Nuno RO Bastos, Rui AC Ferreira, and Delfim FM Torres.  Discrete-time fractional variational problems.  Signal Processing, 91(3):513–524, 2011.
  * [7] Sachin Bhalekar.  Stability and bifurcation analysis of a generalized scalar delay differential equation.  Chaos: An Interdisciplinary Journal of Nonlinear Science, 26(8):084306, 2016.
  * [8] Sachin Bhalekar, Varsha Daftardar-Gejji, Dumitru Baleanu, and Richard Magin.  Fractional bloch equation with delay.  Computers & Mathematics with Applications, 61(5):1355–1365, 2011\.
  * [9] Sachin Bhalekar and Prashant M Gade.  Stability analysis of fixed point of fractional-order coupled map lattices.  Communications in Nonlinear Science and Numerical Simulation, page 106587, 2022.
  * [10] Sachin Bhalekar and Madhuri Patil.  Nonexistence of invariant manifolds in fractional-order dynamical systems.  Nonlinear Dynamics, 102(4):2417–2431, 2020.
  * [11] Oana Brandibur and Eva Kaslik.  Stability results for two-term fractional-order difference equations.  In Proceedings of the International Conference on Fractional Differentiation and its Applications (ICFDA’21), pages 225–230. Springer, 2022\.
  * [12] Jan Čermák, István Győri, and Ludĕk Nechvátal.  On explicit stability conditions for a linear fractional difference system.  Fractional Calculus and Applied Analysis, 18:651–672, 2015.
  * [13] Loredana Camelia Culda, Eva Kaslik, and Mihaela Neamţu.  Stability and bifurcations in a general cournot duopoly model with distributed time delays.  Chaos, Solitons & Fractals, 162:112424, 2022.
  * [14] Varsha Daftardar-Gejji, Sachin Bhalekar, and Prashant Gade.  Dynamics of fractional-ordered chen system with delay.  Pramana, 79(1):61–69, 2012.
  * [15] Marius-F Danca and Michal Fečkan.  Mandelbrot set and julia sets of fractional order.  Nonlinear Dynamics, pages 1–16, 2023.
  * [16] Noureddine Djenina, Adel Ouannas, Iqbal M Batiha, Giuseppe Grassi, Taki-Eddine Oussaeif, and Shaher Momani.  A novel fractional-order discrete sir model for predicting covid-19 behavior.  Mathematics, 10(13):2224, 2022.
  * [17] Mustafa Ali Dokuyucu and Hemen Dutta.  A fractional order model for ebola virus with the new caputo fractional derivative without singular kernel.  Chaos, Solitons & Fractals, 134:109717, 2020.
  * [18] Heiko Enderling and Mark AJ Chaplain.  Mathematical modeling of tumor growth and treatment.  Current pharmaceutical design, 20(30):4934–4940, 2014.
  * [19] Rui AC Ferreira and Delfim FM Torres.  Fractional h-difference equations arising from the calculus of variations.  Applicable Analysis and Discrete Mathematics, pages 110–121, 2011\.
  * [20] Gilberto González-Parra, Abraham J Arenas, and Benito M Chen-Charpentier.  A fractional order epidemic model for the simulation of outbreaks of influenza a (h1n1).  Mathematical methods in the Applied Sciences, 37(15):2218–2226, 2014.
  * [21] Chengdai Huang, Heng Liu, Xiangyun Shi, Xiaoping Chen, Min Xiao, Zhengxin Wang, and Jinde Cao.  Bifurcations in a fractional-order neural network with multiple leakage delays.  Neural Networks, 131:115–126, 2020.
  * [22] Chengdai Huang, Juan Wang, Xiaoping Chen, and Jinde Cao.  Bifurcations in a fractional-order bam neural network with four different delays.  Neural Networks, 141:344–354, 2021.
  * [23] Jia Jia, Xia Huang, Yuxia Li, Jinde Cao, and Ahmed Alsaedi.  Global stabilization of fractional-order memristor-based neural networks with time delay.  IEEE transactions on neural networks and learning systems, 31(3):997–1009, 2019.
  * [24] Divya D Joshi, Prashant M Gade, and Sachin Bhalekar.  Study of low-dimensional nonlinear fractional difference equations of complex order.  Chaos: An Interdisciplinary Journal of Nonlinear Science, 32(11):113101, 2022.
  * [25] Lu Jun-Guo.  Chaotic dynamics of the fractional-order ikeda delay system and its synchronization.  Chinese Physics, 15(2):301, 2006.
  * [26] V Preethi Latha, Fathalla A Rihan, Rajan Rakkiyappan, and G Velmurugan.  A fractional-order delay differential model for ebola infection and cd8+ t-cells response: stability analysis and hopf bifurcation.  International Journal of Biomathematics, 10(08):1750111, 2017.
  * [27] Ch Lubich.  Discretized fractional calculus.  SIAM Journal on Mathematical Analysis, 17(3):704–719, 1986.
  * [28] Kenneth S Miller and Bertram Ross.  Fractional difference calculus.  In Proceedings of the international symposium on univalent functions, fractional calculus and their applications, pages 139–152, 1988.
  * [29] Dorota Mozyrska and Małgorzata Wyrwas.  The z-transform method and delta type fractional difference operators.  Discrete Dynamics in Nature and Society, 2015, 2015.
  * [30] Parvaiz Ahmad Naik, Jian Zu, and Kolade M Owolabi.  Global dynamics of a fractional order model for the transmission of hiv epidemic with optimal control.  Chaos, Solitons & Fractals, 138:109826, 2020.
  * [31] R Joice Nirmala, K Balachandran, L Rodríguez-Germa, and JJ Trujillo.  Controllability of nonlinear fractional delay dynamical systems.  Reports on mathematical physics, 77(1):87–104, 2016.
  * [32] Ivan Niven.  Formal power series.  The American Mathematical Monthly, 76(8):871–889, 1969.
  * [33] AoPS Online.  Binomial theorem.  https://artofproblemsolving.com/wiki/index.php/Binomial_Theorem\\#Generalizations.
  * [34] Edward Ott, Celso Grebogi, and James A Yorke.  Controlling chaos.  Physical review letters, 64(11):1196, 1990.
  * [35] Adel Ouannas, Iqbal M Batiha, and Viet-Thanh Pham.  Fractional Discrete Chaos: Theories, Methods and Applications.  World Scientific, 2023.
  * [36] M.D. PETALE.  Z-Transform: (Theory & Solved Examples).  Engineering Mathematics. MANGESH DEVIDASRAO PETALE.
  * [37] Kestutis Pyragas.  Continuous control of chaos by self-controlling feedback.  Physics letters A, 170(6):421–428, 1992.
  * [38] FA Rihan and G Velmurugan.  Dynamics of fractional-order delay differential model for tumor-immune system.  Chaos, Solitons & Fractals, 132:109592, 2020.
  * [39] H Schiessel, R Metzler, A Blumen, and TF0921 Nonnenmacher.  Generalized viscoelastic models: their fractional equations with solutions.  Journal of physics A: Mathematical and General, 28(23):6567, 1995\.
  * [40] Abdollah Shafieezadeh, Karthik Ramanathan, Jamie E Padgett, and Reginald DesRoches.  Fractional order intensity measures for probabilistic seismic demand modeling applied to highway bridges.  Earthquake Engineering & Structural Dynamics, 41(3):391–409, 2012\.
  * [41] Amar Si-Ammour, Said Djennoune, and Maamar Bettayeb.  A sliding mode control for linear fractional systems with input and state delays.  Communications in Nonlinear Science and Numerical Simulation, 14(5):2310–2318, 2009.
  * [42] David Smith, Lang Moore, et al.  The sir model for spread of disease-the differential equation model.  Convergence, 2004.
  * [43] R Stanisławski and Krzysztof J Latawiec.  Stability analysis for discrete-time fractional-order lti state-space systems. part i: New necessary and sufficient conditions for the asymptotic stability.  Bulletin of the Polish Academy of Sciences: Technical Sciences, (2), 2013.
  * [44] Joel S Steele and Emilio Ferrer.  Latent differential equation modeling of self-regulatory and coregulatory affective processes.  Multivariate Behavioral Research, 46(6):956–984, 2011.
  * [45] HongGuang Sun, Yong Zhang, Dumitru Baleanu, Wen Chen, and YangQuan Chen.  A new collection of real world applications of fractional calculus in science and engineering.  Communications in Nonlinear Science and Numerical Simulation, 64:213–231, 2018.
  * [46] Vladimir V Uchaikin.  Fractional derivatives for physicists and engineers, volume 2.  Springer, 2013.
  * [47] Dongping Wang and Juebang Yu.  Chaos in the fractional order logistic delay system.  In 2008 International Conference on Communications, Circuits and Systems, pages 646–651. IEEE, 2008.
  * [48] Changjin Xu, Maoxin Liao, Peiluan Li, Ying Guo, and Zixin Liu.  Bifurcation properties for fractional order delayed bam neural networks.  Cognitive Computation, 13(2):322–356, 2021.
  * [49] Changjin Xu, Zixin Liu, Peiluan Li, Jinling Yan, and Lingyun Yao.  Bifurcation mechanism for fractional-order three-triangle multi-delayed neural networks.  Neural Processing Letters, pages 1–27, 2022.
  * [50] Changjin Xu, Dan Mu, Zixin Liu, Yicheng Pang, Maoxin Liao, and Chaouki Aouiti.  New insight into bifurcation of fractional-order 4d neural networks incorporating two different time delays.  Communications in Nonlinear Science and Numerical Simulation, 118:107043, 2023.
  * [51] Changjin Xu, Dan Mu, Zixin Liu, Yicheng Pang, Maoxin Liao, Peiluan Li, Lingyun Yao, and Qiwen Qin.  Comparative exploration on bifurcation behavior for integer-order and fractional-order delayed bam neural networks.  Nonlinear Analysis: Modelling and Control, 27:1–24, 2022.
  * [52] Changjin Xu, Wei Zhang, Chaouki Aouiti, Zixin Liu, and Lingyun Yao.  Bifurcation insight for a fractional-order stage-structured predator–prey system incorporating mixed time delays.  Mathematical Methods in the Applied Sciences.
  * [53] Chun Yin, Shou-ming Zhong, and Wu-fan Chen.  Design of sliding mode controller for a class of fractional-order chaotic systems.  Communications in Nonlinear Science and Numerical Simulation, 17(1):356–366, 2012.
  * [54] Zhongli You, Michal Fečkan, and JinRong Wang.  Relative controllability of fractional delay differential equations via delayed perturbation of mittag-leffler functions.  Journal of Computational and Applied Mathematics, 378:112939, 2020\.
  * [55] Wei-Bin Zhang.  Differential Equations, Bifurcations and Chaos in Economics, volume 68.  World Scientific Publishing Company, 2005.
  * [56] Zizhen Zhang, Anwar Zeb, Oluwaseun Francis Egbelowo, and Vedat Suat Erturk.  Dynamics of a fractional order mathematical model for covid-19 epidemic.  Advances in Difference Equations, 2020(1):1–16, 2020.