- Email: [email protected]
On exact and discretized stability of a linear fractional delay differential equation
Journal Pre-proof On exact and discretized stability of a linear fractional delay differential equation
ˇ Jan Cerm´ ak, Ludˇek Nechv´atal
PII: DOI: Reference:
S0893-9659(20)30089-6 https://doi.org/10.1016/j.aml.2020.106296 AML 106296
To appear in:
Applied Mathematics Letters
Received date : 29 November 2019 Revised date : 12 February 2020 Accepted date : 12 February 2020 ˇ Please cite this article as: J. Cerm´ ak and L. Nechv´atal, On exact and discretized stability of a linear fractional delay differential equation, Applied Mathematics Letters (2020), doi: https://doi.org/10.1016/j.aml.2020.106296. This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2020 Elsevier Ltd. All rights reserved.
Journal Pre-proof
On exact and discretized stability of a linear fractional delay differential equation ˇ Jan Cerm´ ak∗, Ludˇek Nechv´ atal
repro of
Institute of Mathematics, Brno University of Technology, Technick´ a 2, CZ-616 69 Brno, Czech Republic
Abstract
The paper discusses the problem of necessary and sufficient stability conditions for a test fractional delay differential equation and its discretization. First, we recall the existing condition for asymptotic stability of the exact equation and consider an appropriate fractional delay difference equation as its discrete counterpart. Then, using the Laplace transform method combined with the boundary locus technique, we derive asymptotic stability conditions in the discrete case as well. Since the studied fractional delay difference equation serves as a backward Euler discretization of the underlying differential equation, we discuss a related problem of numerical stability (with a negative conclusion). Also, as a by-product of our observations, a fractional analogue of the classical Levin–May stability condition is presented. Keywords: Fractional delay differential and difference equation, Asymptotic stability, Numerical stability 2010 MSC: 34K37, 34K20, 39A30, 65L20
1. Introduction and preliminaries
al P
We consider the fractional delay differential equation CD
α
y(t) = λy(t − τ ),
t ∈ (0, ∞),
(1)
where λ ∈ C is a non-zero scalar, τ ∈ R+ is a lag and the symbol C Dα stands for the Caputo differential operator of real order 0 < α < 2. If f : (0, ∞) → C is a function, then this operator is defined by the convolution product Z t α (m) (t) = pm−α−1 (t − ξ)f (m) (ξ) dξ, t > 0 C D f (t) = pm−α−1 ∗ f 0
Jou rn
where m = 1 if 0 < α < 1 and m = 2 if 1 < α < 2. Here, pµ (t) = tµ /Γ(µ + 1), µ ∈ R is the Taylor monomial. If α = 1, we set C D1 f = f 0 . Details on the Caputo approach and other related matters can be found, e.g., in [1] and [2]. In particular, we note that the main preference of the Caputo derivative with respect to other fractional derivatives (like the Riemann–Liouville one) consists in the fact that all initial conditions associated with the studied fractional equation involve classical integer-order derivatives only, having well-known physical meaning in practical applications. Nevertheless, based on existing results, type of fractional derivative operator seems to be insignificant with respect to stability analysis of basic types of fractional differential equations and their discretizations (we shall comment this matter also in the final section). Recently, an enormous development of the qualitative theory of fractional delay differential equations has been initiated. These equations comprise both time delays and fractional order derivative parameters, and thus present a strong tool for mathematical modeling in various areas, in particular, those related to control theory. As a pioneering paper on this topic we mention [3] where the advantage of fractional order derivative as an additional degree of freedom in time delay control of a robotic systems was studied just from the stability point of view (for other related problems see, e.g., [4] and [5]). The application potential has inspired research development especially in stability theory for fractional delay differential equations (see, e.g., [6], [7] and [8]) whose basic foundations are now, at least partially, established. On the contrary, the corresponding theory for their discrete counterparts is just at the edge of existing knowledge. In particular, stability criteria in this area are very rare (for introductory results see, e.g., [9] and [10]). Therefore, the aim of this paper is to answer one of the key stability problems related to fractional delay difference equations. Following these outlines, the paper is organized as follows. Section 2 recapitulates the basic stability criterion for a test one-term fractional delay differential equation with a complex coefficient. The core of the paper is ∗ Corresponding
author ˇ Email addresses: [email protected] (Jan Cerm´ ak), [email protected] (Ludˇ ek Nechv´ atal)
Preprint submitted to Applied Mathematics Letters
February 12, 2020
Journal Pre-proof
repro of
contained in Section 3 devoted to stability analysis of a fractional delay difference equation that forms a direct discrete counterpart to the studied exact equation. Here, we first introduce some notions and tools related to discrete fractional calculus, and then apply the Laplace transform method and boundary locus technique to describe the stability region for the investigated difference equation. Also, a fractional analogue of the famous Levin–May stability condition from the area of discrete population dynamics is presented here. Section 4 points out that, contrary to the integer-order case, the backward Euler discretization of the underlying fractional delay differential equations is not τ -stable. Finally, Section 5 contains some summarizing remarks. 2. Stability of the exact equation
It is well known that (1) is asymptotically stable (i.e., any solution of (1) is eventually tending to the zero solution) if and only if the associated characteristic equation sα − λ exp(−sτ ) = 0
(2)
has roots with negative real parts only. When 0 < α < 1, an explicit reformulation of this condition was performed (as a part of some more general observations) in [11] for λ real, and in [12] for λ complex. In the latter case, the application of the modified D-decomposition method resulted into the stability set α απ |Arg(z)| − απ/2 , |Arg(z)| > Σα = z ∈ C : |z| < τ τ 2 (where −π < Arg(·) ≤ π) describing the set of all complex values λ such that (1) is asymptotically stable. The extension of this result for 1 < α < 2 can be done by use of a similar technique (see also [13]). Thus, as a starting point for our next investigations, we have Theorem 1. Let 0 < α < 2 and τ > 0 be reals. Then (1) is asymptotically stable if and only if λ ∈ Σα τ. Remark 2. Letting α → 2− , the region Σα τ becomes empty for any positive real lag τ and the stability is lost.
al P
The main goal of this paper is to consider a discrete analogue of (1) and find an appropriate discretized stability set (more precisely, a discrete counterpart to Theorem 1). This will be done in the next section where some essentials on discrete fractional h-calculus are mentioned first. 3. Stability of the discretized equation
In what follows, we use the notation hZj = {jh, (j + 1)h, . . . }, where h ∈ R+ (discretization stepsize) and j ∈ Z. In addition, we denote hZ+ = hZ1 . If k ∈ Z+ and f : hZ1−k → C, then a delayed (shifted) function f ρk : hZ+ → C of f is defined as f ρk (hn) = f h(n − k) , n = 1, 2, . . . .
Jou rn
Any function f : hZj → C will be associated with a sequence (fn ) defined as fn = f (nh), n = j, j + 1, . . .. One of the possible discretizations of (1) consists of considering the equation on the uniform mesh hZ+ and replacing the fractional operator C Dα by the backward Caputo fractional difference operator C∇α h defined for 0 < α < 2 and a function f : hZ−1 → C as α C∇h fn
(m−α−1)
= (ph
∗ ∇hm f )n ,
n = 1, 2, . . .
with m = 1 if 0 < α ≤ 1 and m = 2 if 1 < α < 2. Here, ∇h1 fn = ∇h fn = h−1 (fn − fn−1 )
and
∇h2 fn = ∇h (∇h f )n = h−2 (fn − 2fn−1 + fn−2 ),
n = 1, 2, . . .
are the standard first and second backward differences, “∗” stands for a convolution product defined as (f ∗ g)n = h
n X
fj gn−j+1 ,
n = 1, 2, . . . ,
j=1
f, g ∈ hZ+ ,
(µ)
and the function ph : hZ+ → R represents a discrete analogue to function pµ . For any µ ∈ R, it is defined as −µ − 1 x x(x − 1)(x − 2) . . . (x − j + 1) (µ) (ph )n = hµ (−1)n−1 where = , x ∈ R, j ∈ Z+ n−1 j j! is the binomial coefficient. Note that when α = 1, then C∇α h f becomes just ∇h f . For other basic information on nabla fractional differences and related notions (see, e.g., [14]). 2
Journal Pre-proof
Assuming τ = kh for a suitable k ∈ Z+ , the discretization of (1) has the form α C∇h yn
= λynρk ,
n = 1, 2, . . . ,
(3)
repro of
λ ∈ C\{0}. Note that the solution y : hZ+ → C of (3) now depends on a k-tuple of initial values y1−k , y2−k , . . . , y0 (and on ∇h y0 if 1 < α < 2). It can be seen that we have a unique existence of the solution as (3) actually represents an explicit formula for the unknown yn . The question of asymptotic stability of (3) is answered by the following assertion. Theorem 3. Let 0 < α < 2, h > 0 be reals and let k > 0 be an integer. Then (3) is asymptotically stable if α λ ∈ Sh,k where α Sh,k =
z ∈ C : |z| <
απ 2α α kπ − |Arg(z)| cos , |Arg(z)| > . hα 2k − α 2
To prove Theorem 3, we employ the Laplace transform method combined with the boundary locus technique. The Laplace transform of a function f : hZ+ → C is defined by a power series centered at the point s0 = h−1 as Lh {f }(s) = h
∞ X j=1
fj (1 − hs)j−1
for those s ∈ C where the series is convergent.
The well-known fact states that if the power series converges at some s 6= h−1 , then there exists r > 0 such that it converges on disk D(h−1 , r) = {s ∈ C : |s − h−1 | < r}. Moreover, Lh {f } is determined uniquely, i.e., two different sequences cannot have the same Laplace image. The next proposition collects a few useful properties (for details, see [15] and the references therein). Proposition 4. Let f, g : hZ+ → C be functions such that Lh {f } and Lh {g} converge on D(h−1 , rf ) and D(h−1 , rg ), respectively. Then (i) Lh {f ∗ g}(s) = Lh {f }(s) · Lh {g}(s) on D(h−1 , r) where r = min{rf , rg }; (µ)
al P
(ii) Lh {ph }(s) = s−µ−1 on D(h−1 , h−1 ) for any µ ∈ R; (iii) Lh {∇h f }(s) = sLh {f }(s) − f0 provided f is defined on hZ0 . Using these three properties, we immediately have α α−1 Lh {C∇α f0 h f }(s) = s Lh {f }(s) − s
α α−1 or Lh {C∇α f0 − sα−2 ∇h f0 h f }(s) = s Lh {f }(s) − s
(4)
Jou rn
for 0 < α ≤ 1 or 1 < α < 2, respectively, and for any f : hZ−1 → C such that Lh {f } converges on some D(h−1 , r), r > 0 (note that the result formally agrees with the continuous case, i.e., when the standard Laplace transform is applied to the Caputo derivative of a function f : (0, ∞) → C). Finally, it is easy to check that Lh {f ρk }(s) = (1 − hs)k Lh {f }(s) + h
k X j=1
f1−j (1 − hs)k−j ,
(5)
hence, the formulae (4) and (5) yield the Laplace image of the solution y in the form k −1 X Lh {y}(s) = sα−1 y0 + sα−2 ∇h y0 + λh y1−j (1 − hs)k−j · sα − λ(1 − hs)k
(6)
j=1
where the term sα−2 ∇h y0 vanishes if 0 < α ≤ 1. The equation sα − λ(1 − hs)k = 0
(7)
is called a characteristic equation of (3) (and forms a discrete counterpart to (2)). As a simple consequence of the Cauchy–Hadamard theorem we can state: Proposition 5. Let f : hZ+ → C be a function such that Lh {f } converges on D(h−1 , r) with r > 0. (i) If r > h−1 , then f ∈ `1 ; (ii) If r < h−1 , then |f | is not bounded from above, i.e., lim supn→∞ |fn | = ∞.
The behavior in the remaining case when r = h−1 is covered by the following version of the Wiener theorem (see [16] and [15]): 3
Journal Pre-proof
Proposition 6. Let g ∈ `1 . Then there exists f ∈ `1 satisfying Lh {f }(s) · Lh {g}(s) = 1 if and only if inf |Lh {g}(s)| : s ∈ cl D(h−1 , h−1 ) > 0.
(8)
Lh {u}(s) = sα−1 y0 + sα−2 ∇h y0 + λh
k X j=1
repro of
Proof of Theorem 3. A classical result from complex analysis claims that the radius of convergence of a power series is given by the distance from the series’ center to the nearest singular point. Observe that Lh {y} cannot have the radius of convergence greater than h−1 because of the presence of power function sα in Lh {y} (the origin is a singular point). On the other hand, if (7) has a root in D(h−1 , h−1 ), then Lh {y} has the radius of convergence smaller than h−1 and y is not asymptotically stable according to Proposition 5. Hence, let us further investigate the case when (7) has no root in D(h−1 , h−1 ), i.e., r = h−1 . In view of (6), the solution y can be written as a convolution y = u ∗ v where y1−j (1 − hs)k−j
and
Lh {v}(s) = sα − λ(1 − hs)k
−1
.
Let us show that v ∈ `1 if (7) has no root on cl D(h−1 , h−1 ) . Defining an auxiliary function v˜ : hZ+ → C as (−α−1)
v˜ = ph
we have
− λδ
where
δn = h−1
for
n = k + 1 and
δn = 0 otherwise,
−1 = 1. Lh {˜ v }(s) · Lh {v}(s) = sα − λ(1 − hs)k · sα − λ(1 − hs)k (−α−1)
al P
Moreover, v˜ ∈ `1 because of the asymptotic equivalence (ph )n ∼ Cn−α−1 as n → ∞ for a suitable C ∈ R (in (−α−1) fact, the exact sum of ph can be calculated). Finally, Lh {˜ v } is continuous and never zero on cl D(h−1 , h−1 ) since (7) has no root here by the assumption. Hence, by the Weierstrass theorem, we have also inf |Lh {˜ y }(s)| = min |Lh {˜ y }(s)| > 0 on cl D(h−1 , h−1 ) and Proposition 6 gives v ∈ `1 . The function u decays more slowly than v because of the presence of the term sα−1 (if 0 < α < 1) or sα−2 (ν−1) (ν−1) (if 1 < α < 2). Indeed, s−ν (0 < ν < 1) represents the Laplace image of ph for which (ph )n ∼ Cnν−1 (ν−1) 1 as n → ∞ for a suitable C ∈ R, hence, ph is not in ` . However, it is in a larger space, namely in `p where (ν−1) (1 − ν)−1 < p < ∞. Thus, by the Young’s inequality for convolution, ph ∗ v is also in `p and, consequently, it converges to zero. The remaining terms in Lh {u} correspond to summable functions which means that their convolution with v is again summable, thus converging to zero. Overall, we have proved that |yn | → 0 as n → ∞ provided (7) has no root on cl D(h−1 , h−1 ) . Note that we do −1 −1 not know stability behavior when (7) has a root on ∂D(h , h ) (in this case, the property (8) of Proposition 6 is not satisfied). Now we turn our attention to discussing conditions guaranteeing that characteristic equation (7) has a root on ∂D(h−1 , h−1 ). Let s be such a root, i.e., s = 2h−1 cos(ψ) exp(iψ)
for a suitable − π/2 < ψ ≤ π/2.
(9)
Jou rn
This s can equivalently be expressed as
s = h−1 1 + exp(i2ψ) .
(10)
Since the term (1 − hs)k is never zero for s ∈ ∂D(h−1 , h−1 ), the characteristic equation (7) can be rewritten as λ = g(s) = sα (1 − hs)−k .
(11)
Substituting (9) to g, we arrive (taking into account (10)) at g ∗ (ψ) = g(s)|s=2h−1 cos(ψ) exp(iψ) = (−1)k 2α h−α cosα (ψ) exp i(2k − α)ψ ,
−π/2 < ψ ≤ π/2.
Function g ∗ takes the values on a pair of curves lying in the complex plane that are symmetric with respect to the real axis. Clearly, the modulus |g ∗ (ψ)| is proportional to the value cosα (ψ) and thus is increasing when |ψ| moves from π/2 towards zero (because of cos is decreasing on [0, π/2]). It means that each of these two curves do not have self-intersections. Moreover, it crosses the real axis if ψ = jπ/(2k − α), j = 0, 1, . . . , k − 1 and if ψ = π/2. The crossing, which is nearest to the origin, occurs for j = k − 1 and lies on the negative real half-axis (since (−1)k exp(i(k − 1)π) = (−1)k (−1)k exp(−iπ) = −1). Hence, the restriction (k − 1)π/(2k − α) ≤ |ψ| ≤ π/2 causes that g ∗ just represents a simple closed curve in the complex plane. Substitution ψ = (kπ − |ϕ|)/(2k − α) 4
Journal Pre-proof
transforms g ∗ to g ∗∗ (ϕ) = 2α h−α cosα (kπ − |ϕ|)/(2k − α) exp(iϕ),
απ/2 ≤ |ϕ| ≤ π.
repro of
∗∗ The set Γα h,k = {g (ϕ) ∈ C : απ/2 ≤ |ϕ| ≤ π} splits (as a simple closed curve) the complex plane to two disjoint parts – an interior and an exterior. It can be shown that the preimage of the interior under the function g appearing in (11) is located outside of cl D(h−1 , h−1 ) and, conversely, the preimage of the exterior under g is just D(h−1 , h−1 ) \ {h−1 }. In the latter case, the exterior of Γα h,k is actually the codomain of g over D(h−1 , h−1 ) \ {h−1 }. The proof of this step requires all the parts of the function g ∗ (not only that one due to the restriction (k − 1)π/(2k − α) ≤ |ψ| ≤ π/2) and utilizes the open mapping theorem together with some continuity arguments (some technical refinements are needed as g is not holomorphic at h−1 and at the non-positive real half-axis), details how to proceed can be found in [17]. The above facts immediately imply that any choice of −1 λ belonging to the interior of Γα , h−1 ) , hence, (3) is asymptotically h,k causes that there is no root on cl D(h stable according to the previous part. Contrary, taking λ lying in the exterior of Γα h,k , we have a root of (7) located in D(h−1 , h−1 ) and thus (3) is unstable.
Remark 7. The stability condition in Theorem 3 is formulated as a sufficient one. To prove its necessity, one α has to dispose with the situation on the stability boundary ∂Sh,k which seems to be a much more complicated matter compared to the continuous counterpart (we remind that outside the closure of this area we have the instability). Both the exact and discrete stability sets are depicted in Figure 1. =(z) −
2α α (k − 1)π cos hα 2k − α
απ 2
α
απ 2
al P
π − απ/2 τ
απ 2
Σα τ
<(z)
=(z)
2α (k − 1)π cosα hα 2k − α
α Sh,k
Σα τ α Sh,k
−
−
−
−
π − απ/2 τ
α
<(z)
−
απ 2
α Figure 1: Exact and discrete asymptotic stability sets Σα τ and Sh,k . The figure corresponds to h = 0.5, k = 2 (hence τ = 1), α = 0.8 (left) and α = 1.1 (right), the scales of the axes are equal.
Remark 8. If we consider λ real and put h = 1 (the standard normalized discrete case), then the stability α property λ ∈ S1,k is reduced to the condition
Jou rn
−2α cosα (k − 1)π/(2k − α) < λ < 0
that provides a fractional generalization of the classical Levin–May stability condition −2 cos (k − 1)π)/(2k − 1) < λ < 0
for the delay difference equation
∇1 yn = λyn−k ,
n = 1, 2, . . .
(λ being a real scalar and k a positive integer) appearing in some problems of discrete population biology (see [18]). 4. Numerical stability of the backward Euler discretization Equation (1) is the simplest test equation that enables a generalization of the concept of A-stability to fractional delay differential equations. In the first-order delayed case, the notion of τ -stability was introduced in the following sense (see, e.g., [19]): A Θ-method for delay differential equations is called τ -stable if its appropriate difference scheme is asymptotically stable for any constant stepsize h = τ /k (k being a positive integer) whenever a test exact equation is asymptotically stable. Considering (1) with arbitrary 0 < α < 2 as a test exact equation, one can extend the notion of τ -stability for fractional delay differential equations as well. 5
Journal Pre-proof
Since (3) can be viewed as the backward Euler method (when θ = 1) applied to (1), we use our results to α + discuss its τ -stability. In other words, we need to verify the property Σα τ ⊂ Sh,k for any h = τ /k, k ∈ Z . Using the above explicit descriptions of corresponding stability regions, this property is equivalent to the inequality d(ϕ) = cos (kπ − ϕ)/(2k − α) − (ϕ − απ/2)/(2k) ≥ 0, απ/2 ≤ ϕ ≤ π (12)
=(z)
0.5 S1,1
p − 3π/4
repro of
that should be satisfied for any positive integer k. Nevertheless, this is not true. Indeed, if, e.g., α = 0.5 and 0.5 τ = h = k = 1, then Figure 2 clearly demonstrates that Σ0.5 is not a part of S1,1 (analytically, (12) does not 1 hold because of d(απ/2) = 0, d0 (απ/2) > 0 and d(π) < 0, hence, d(ϕ) ≥ 0 in a right neighborhood of απ/2 and d(ϕ) < 0 in a left neighborhood of π).
Σ0.5 1
π 4
<(z)
√ − 2
−
π 4
0.5 , the former is not a subset of the latter. Figure 2: The stability sets Σ0.5 and S1,1 1
5. Concluding remarks
al P
Consequently, the backward Euler method is not τ -stable. Note that this conclusion contradicts the situation in the first-order case when the backward Euler method embodies very strong stability properties (in fact, it is τ -stable even when two-term delay differential equation with complex coefficients is taken as a test equation, see [20]).
Jou rn
We have analyzed exact and discretized stability of the fractional delay differential equation (1). The stability region for its discretization (3) has been described explicitly and compared with the exact stability region. This comparison resulted into the conclusion that the backward Euler method for fractional delay differential equations is not τ -stable. The mathematical essence of this paper consists in two aspects. First, it is a derivation of the characteristic equation (7). More precisely, we have shown that (3) is asymptotically stable if its complex parameter λ lies outside the disk D(h−1 , h−1 ) (and (3) is not stable provided λ lies inside this disk). Till now, this theoretical condition has been known in the differential case only (see, e.g., [21] and [22]). Secondly, Theorem 3 provides an α effective reformulation of this condition and describes the appropriate discrete stability region Sh,k . From this point of view, Theorem 3 plays the same role as the corresponding effective conditions known from the differential case (see [11], [12] and partially also [8]). Regarding very few criteria on stability of fractional delay difference equations, some particular results on finite-time stability and Mittag-Leffler stability presented in [9] and [10] do not answer this fundamental stability problem discussed in the previous sections. Keeping in mind the numerical aspect of the studied problem, we note that there are some other papers suggesting suitable numerical schemes for solving fractional delay differential equations and answering typical questions such as the error estimate or the convergence rate (see, e.g., [23]). However, due to advanced form of the proposed methods, the analytical description of stability regions (similar to that performed in Section 3) is missing. Since the backward Euler scheme is, at least in the case of purely delay fractional differential equations, (µ) easy to implement (we note that the binomial coefficient in the function ph can easily be calculated recursively), the full description of its stability region is important. In particular, the inequality (12) can be useful when determining the maximal stepsize h such that the appropriate numerical solution can retain the key qualitative property (convergence to zero) of the exact solution. Possible extensions of Theorem 3 can be directed in several ways. Of course, one can consider (3) with other types of fractional operator. Some existing results in this area (see, e.g., [22]) indicate that the stability region derived for the Caputo and Riemann–Liouville fractional operators remains the same – up to the situation on the stability boundary where dynamics of the same equation with different fractional operators exhibits different stability properties. Further, (1) extended by a non-delayed term, i.e., the equation CD
α
y(t) = λy(t − τ ) + νy(t), 6
t ∈ (0, ∞)
Journal Pre-proof
with complex λ and ν, serves as a delay feedback control model enabling, among others, stabilization and synchronization of significant solutions of various differential systems (see, e.g., [4]). Thus, considering the corresponding Euler discretization (i.e., (3) extended by a non-delayed term), one can study and answer many issues connected with fractional control of discrete dynamical systems. Finally, another possible future research connected with the conclusions of this paper is related to stability analysis of more advanced numerical schemes for fractional delay differential equations.
repro of
Acknowledgement. The research has been supported by the grant 17-03224S of the Czech Science Foundation. The authors are grateful to the anonymous referees for valuable comments and suggestions that helped to improve the content of the paper. References
Jou rn
al P
[1] K. Diethelm, The Analysis of Fractional Differential Equations: An Application-Oriented Exposition Using Differential Operators of Caputo Type, Springer, Berlin, 2010. [2] Y. Zhou, Basic Theory of Fractional Differential Equations, World Scientific, 2014. [3] M. Lazarevi´ c, Finite time stability analysis of P Dα fractional control of robot time-delay systems. Mech. Res. Commun. 33 (2006) 269–279. [4] M. Lazarevi´ c, Stability and stabilization of fractional order time delay systems. Scientific Technical Review 61 (2011) 31–44. [5] B. Tao, M. Xiao, Q. Sun, J. Cao, Hopf bifurcation analysis of a delayed fractional-order genetic regulatory network model, Neurocomputing 275 (2018) 677–686. [6] S. Bhalekar, Stability analysis of a class of fractional delay differential equations, Pramana-J. Phys. 81 (2013) 215–224. [7] E. Kaslik, S. Sivasundaram, Analytical and numerical methods for the stability analysis of linear fractional delay differential equations, J. Comput. Appl. Math. 236 (2012) 4027–4041. [8] X. Teng, Z. Wang, Stability switches of a class of fractional-delay systems with delay-dependent coefficients, J. Comput. Nonlinear Dynam. 13 (2018) 111005 (9 pages). [9] F. Du, B. Jia, Finite-time stability of a class of nonlinear fractional delay difference systems, Appl. Math. Lett. 98 (2019) 233–239. [10] G.-C. Wu, D. Baleanu, L.L. Huang, Novel Mittag-Leffler stability of linear fractional delay difference equations with impulse, Appl. Math. Lett. 82 (2018) 71–78. [11] K. Krol, Asymptotic properties of fractional delay differential equations, Appl. Math. Comput. 218 (2011) 1515–1532. ˇ [12] J. Cerm´ ak, J. Horn´ıˇ cek, T. Kisela, Stability regions for fractional differential systems with a time delay, Commun. Nonlinear Sci. Numer. Simul. 31 (2016) 108–123. ˇ [13] J. Cerm´ ak, T. Kisela, Oscillatory and asymptotic properties of fractional delay differential equations, Electron. J. Differ. Eq. 2019 (2019) (15 pages). [14] F.M. Atici, P.W. Eloe, Linear systems of fractional nabla difference equations, Rocky Mountain J. Math. 41 (2011) 353–370. ˇ [15] J. Cerm´ ak, T. Kisela, L. Nechv´ atal, Stability regions for linear fractional differential systems and their discretizations, Appl. Math. Comput. 219 (2013) 7012–7022. [16] C.A. Desoer, M. Vidyassager, Feedback Systems: Input–Output Properties, Academic Press, New York, 1975. ˇ [17] J. Cerm´ ak, I. Gy¨ ori, L. Nechv´ atal, On explicit stability conditions for a linear fractional difference system, Fract. Calc. Appl. Anal. 18 (2015) 651–672. [18] S.A. Levin, R. May, A note on difference delay equations, Theor. Popul. Biol. 9 (1976) 178–187. [19] N. Guglielmi, Delay dependent stability regions of Θ-methods for delay differential equations, IMA J. Numer. Anal. 18 (1998) 399–418. [20] S. Maset, Stability of Runge Kutta methods for linear delay differential equations, Numer. Math. 87 (2000) 355–371. [21] C.P. Li, F.R. Zhang, A survey on the stability of fractional differential equations, Eur. Phys. J. Special Topics 193 (2011), 27–47. [22] D. Qian, Ch. Li, R. P. Agarwal, P. J. Y. Wong, Stability analysis of fractional differential system with Riemann-Liouville derivative, Math. Comput. Model. 52 (2010), 862–874. [23] A. Jhinga, V. Daftardar-Gejji, A new numerical method for solving fractional delay differential equations, Comput. Appl. Math. 38 (2019) 166 (18 pages).
7
Journal Pre-proof
Jou rn
al P
repro of
Jan Čermák: Conceptualization, Investigation, Methodology, Writing – original draft, Writing – review & editing. Luděk Nechvátal: Conceptualization, Investigation, Methodology, Writing – original draft, Writing – review & editing.
ˇ Jan Cerm´ ak, Ludˇek Nechv´atal
PII: DOI: Reference:
S0893-9659(20)30089-6 https://doi.org/10.1016/j.aml.2020.106296 AML 106296
To appear in:
Applied Mathematics Letters
Received date : 29 November 2019 Revised date : 12 February 2020 Accepted date : 12 February 2020 ˇ Please cite this article as: J. Cerm´ ak and L. Nechv´atal, On exact and discretized stability of a linear fractional delay differential equation, Applied Mathematics Letters (2020), doi: https://doi.org/10.1016/j.aml.2020.106296. This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2020 Elsevier Ltd. All rights reserved.
Journal Pre-proof
On exact and discretized stability of a linear fractional delay differential equation ˇ Jan Cerm´ ak∗, Ludˇek Nechv´ atal
repro of
Institute of Mathematics, Brno University of Technology, Technick´ a 2, CZ-616 69 Brno, Czech Republic
Abstract
The paper discusses the problem of necessary and sufficient stability conditions for a test fractional delay differential equation and its discretization. First, we recall the existing condition for asymptotic stability of the exact equation and consider an appropriate fractional delay difference equation as its discrete counterpart. Then, using the Laplace transform method combined with the boundary locus technique, we derive asymptotic stability conditions in the discrete case as well. Since the studied fractional delay difference equation serves as a backward Euler discretization of the underlying differential equation, we discuss a related problem of numerical stability (with a negative conclusion). Also, as a by-product of our observations, a fractional analogue of the classical Levin–May stability condition is presented. Keywords: Fractional delay differential and difference equation, Asymptotic stability, Numerical stability 2010 MSC: 34K37, 34K20, 39A30, 65L20
1. Introduction and preliminaries
al P
We consider the fractional delay differential equation CD
α
y(t) = λy(t − τ ),
t ∈ (0, ∞),
(1)
where λ ∈ C is a non-zero scalar, τ ∈ R+ is a lag and the symbol C Dα stands for the Caputo differential operator of real order 0 < α < 2. If f : (0, ∞) → C is a function, then this operator is defined by the convolution product Z t α (m) (t) = pm−α−1 (t − ξ)f (m) (ξ) dξ, t > 0 C D f (t) = pm−α−1 ∗ f 0
Jou rn
where m = 1 if 0 < α < 1 and m = 2 if 1 < α < 2. Here, pµ (t) = tµ /Γ(µ + 1), µ ∈ R is the Taylor monomial. If α = 1, we set C D1 f = f 0 . Details on the Caputo approach and other related matters can be found, e.g., in [1] and [2]. In particular, we note that the main preference of the Caputo derivative with respect to other fractional derivatives (like the Riemann–Liouville one) consists in the fact that all initial conditions associated with the studied fractional equation involve classical integer-order derivatives only, having well-known physical meaning in practical applications. Nevertheless, based on existing results, type of fractional derivative operator seems to be insignificant with respect to stability analysis of basic types of fractional differential equations and their discretizations (we shall comment this matter also in the final section). Recently, an enormous development of the qualitative theory of fractional delay differential equations has been initiated. These equations comprise both time delays and fractional order derivative parameters, and thus present a strong tool for mathematical modeling in various areas, in particular, those related to control theory. As a pioneering paper on this topic we mention [3] where the advantage of fractional order derivative as an additional degree of freedom in time delay control of a robotic systems was studied just from the stability point of view (for other related problems see, e.g., [4] and [5]). The application potential has inspired research development especially in stability theory for fractional delay differential equations (see, e.g., [6], [7] and [8]) whose basic foundations are now, at least partially, established. On the contrary, the corresponding theory for their discrete counterparts is just at the edge of existing knowledge. In particular, stability criteria in this area are very rare (for introductory results see, e.g., [9] and [10]). Therefore, the aim of this paper is to answer one of the key stability problems related to fractional delay difference equations. Following these outlines, the paper is organized as follows. Section 2 recapitulates the basic stability criterion for a test one-term fractional delay differential equation with a complex coefficient. The core of the paper is ∗ Corresponding
author ˇ Email addresses: [email protected] (Jan Cerm´ ak), [email protected] (Ludˇ ek Nechv´ atal)
Preprint submitted to Applied Mathematics Letters
February 12, 2020
Journal Pre-proof
repro of
contained in Section 3 devoted to stability analysis of a fractional delay difference equation that forms a direct discrete counterpart to the studied exact equation. Here, we first introduce some notions and tools related to discrete fractional calculus, and then apply the Laplace transform method and boundary locus technique to describe the stability region for the investigated difference equation. Also, a fractional analogue of the famous Levin–May stability condition from the area of discrete population dynamics is presented here. Section 4 points out that, contrary to the integer-order case, the backward Euler discretization of the underlying fractional delay differential equations is not τ -stable. Finally, Section 5 contains some summarizing remarks. 2. Stability of the exact equation
It is well known that (1) is asymptotically stable (i.e., any solution of (1) is eventually tending to the zero solution) if and only if the associated characteristic equation sα − λ exp(−sτ ) = 0
(2)
has roots with negative real parts only. When 0 < α < 1, an explicit reformulation of this condition was performed (as a part of some more general observations) in [11] for λ real, and in [12] for λ complex. In the latter case, the application of the modified D-decomposition method resulted into the stability set α απ |Arg(z)| − απ/2 , |Arg(z)| > Σα = z ∈ C : |z| < τ τ 2 (where −π < Arg(·) ≤ π) describing the set of all complex values λ such that (1) is asymptotically stable. The extension of this result for 1 < α < 2 can be done by use of a similar technique (see also [13]). Thus, as a starting point for our next investigations, we have Theorem 1. Let 0 < α < 2 and τ > 0 be reals. Then (1) is asymptotically stable if and only if λ ∈ Σα τ. Remark 2. Letting α → 2− , the region Σα τ becomes empty for any positive real lag τ and the stability is lost.
al P
The main goal of this paper is to consider a discrete analogue of (1) and find an appropriate discretized stability set (more precisely, a discrete counterpart to Theorem 1). This will be done in the next section where some essentials on discrete fractional h-calculus are mentioned first. 3. Stability of the discretized equation
In what follows, we use the notation hZj = {jh, (j + 1)h, . . . }, where h ∈ R+ (discretization stepsize) and j ∈ Z. In addition, we denote hZ+ = hZ1 . If k ∈ Z+ and f : hZ1−k → C, then a delayed (shifted) function f ρk : hZ+ → C of f is defined as f ρk (hn) = f h(n − k) , n = 1, 2, . . . .
Jou rn
Any function f : hZj → C will be associated with a sequence (fn ) defined as fn = f (nh), n = j, j + 1, . . .. One of the possible discretizations of (1) consists of considering the equation on the uniform mesh hZ+ and replacing the fractional operator C Dα by the backward Caputo fractional difference operator C∇α h defined for 0 < α < 2 and a function f : hZ−1 → C as α C∇h fn
(m−α−1)
= (ph
∗ ∇hm f )n ,
n = 1, 2, . . .
with m = 1 if 0 < α ≤ 1 and m = 2 if 1 < α < 2. Here, ∇h1 fn = ∇h fn = h−1 (fn − fn−1 )
and
∇h2 fn = ∇h (∇h f )n = h−2 (fn − 2fn−1 + fn−2 ),
n = 1, 2, . . .
are the standard first and second backward differences, “∗” stands for a convolution product defined as (f ∗ g)n = h
n X
fj gn−j+1 ,
n = 1, 2, . . . ,
j=1
f, g ∈ hZ+ ,
(µ)
and the function ph : hZ+ → R represents a discrete analogue to function pµ . For any µ ∈ R, it is defined as −µ − 1 x x(x − 1)(x − 2) . . . (x − j + 1) (µ) (ph )n = hµ (−1)n−1 where = , x ∈ R, j ∈ Z+ n−1 j j! is the binomial coefficient. Note that when α = 1, then C∇α h f becomes just ∇h f . For other basic information on nabla fractional differences and related notions (see, e.g., [14]). 2
Journal Pre-proof
Assuming τ = kh for a suitable k ∈ Z+ , the discretization of (1) has the form α C∇h yn
= λynρk ,
n = 1, 2, . . . ,
(3)
repro of
λ ∈ C\{0}. Note that the solution y : hZ+ → C of (3) now depends on a k-tuple of initial values y1−k , y2−k , . . . , y0 (and on ∇h y0 if 1 < α < 2). It can be seen that we have a unique existence of the solution as (3) actually represents an explicit formula for the unknown yn . The question of asymptotic stability of (3) is answered by the following assertion. Theorem 3. Let 0 < α < 2, h > 0 be reals and let k > 0 be an integer. Then (3) is asymptotically stable if α λ ∈ Sh,k where α Sh,k =
z ∈ C : |z| <
απ 2α α kπ − |Arg(z)| cos , |Arg(z)| > . hα 2k − α 2
To prove Theorem 3, we employ the Laplace transform method combined with the boundary locus technique. The Laplace transform of a function f : hZ+ → C is defined by a power series centered at the point s0 = h−1 as Lh {f }(s) = h
∞ X j=1
fj (1 − hs)j−1
for those s ∈ C where the series is convergent.
The well-known fact states that if the power series converges at some s 6= h−1 , then there exists r > 0 such that it converges on disk D(h−1 , r) = {s ∈ C : |s − h−1 | < r}. Moreover, Lh {f } is determined uniquely, i.e., two different sequences cannot have the same Laplace image. The next proposition collects a few useful properties (for details, see [15] and the references therein). Proposition 4. Let f, g : hZ+ → C be functions such that Lh {f } and Lh {g} converge on D(h−1 , rf ) and D(h−1 , rg ), respectively. Then (i) Lh {f ∗ g}(s) = Lh {f }(s) · Lh {g}(s) on D(h−1 , r) where r = min{rf , rg }; (µ)
al P
(ii) Lh {ph }(s) = s−µ−1 on D(h−1 , h−1 ) for any µ ∈ R; (iii) Lh {∇h f }(s) = sLh {f }(s) − f0 provided f is defined on hZ0 . Using these three properties, we immediately have α α−1 Lh {C∇α f0 h f }(s) = s Lh {f }(s) − s
α α−1 or Lh {C∇α f0 − sα−2 ∇h f0 h f }(s) = s Lh {f }(s) − s
(4)
Jou rn
for 0 < α ≤ 1 or 1 < α < 2, respectively, and for any f : hZ−1 → C such that Lh {f } converges on some D(h−1 , r), r > 0 (note that the result formally agrees with the continuous case, i.e., when the standard Laplace transform is applied to the Caputo derivative of a function f : (0, ∞) → C). Finally, it is easy to check that Lh {f ρk }(s) = (1 − hs)k Lh {f }(s) + h
k X j=1
f1−j (1 − hs)k−j ,
(5)
hence, the formulae (4) and (5) yield the Laplace image of the solution y in the form k −1 X Lh {y}(s) = sα−1 y0 + sα−2 ∇h y0 + λh y1−j (1 − hs)k−j · sα − λ(1 − hs)k
(6)
j=1
where the term sα−2 ∇h y0 vanishes if 0 < α ≤ 1. The equation sα − λ(1 − hs)k = 0
(7)
is called a characteristic equation of (3) (and forms a discrete counterpart to (2)). As a simple consequence of the Cauchy–Hadamard theorem we can state: Proposition 5. Let f : hZ+ → C be a function such that Lh {f } converges on D(h−1 , r) with r > 0. (i) If r > h−1 , then f ∈ `1 ; (ii) If r < h−1 , then |f | is not bounded from above, i.e., lim supn→∞ |fn | = ∞.
The behavior in the remaining case when r = h−1 is covered by the following version of the Wiener theorem (see [16] and [15]): 3
Journal Pre-proof
Proposition 6. Let g ∈ `1 . Then there exists f ∈ `1 satisfying Lh {f }(s) · Lh {g}(s) = 1 if and only if inf |Lh {g}(s)| : s ∈ cl D(h−1 , h−1 ) > 0.
(8)
Lh {u}(s) = sα−1 y0 + sα−2 ∇h y0 + λh
k X j=1
repro of
Proof of Theorem 3. A classical result from complex analysis claims that the radius of convergence of a power series is given by the distance from the series’ center to the nearest singular point. Observe that Lh {y} cannot have the radius of convergence greater than h−1 because of the presence of power function sα in Lh {y} (the origin is a singular point). On the other hand, if (7) has a root in D(h−1 , h−1 ), then Lh {y} has the radius of convergence smaller than h−1 and y is not asymptotically stable according to Proposition 5. Hence, let us further investigate the case when (7) has no root in D(h−1 , h−1 ), i.e., r = h−1 . In view of (6), the solution y can be written as a convolution y = u ∗ v where y1−j (1 − hs)k−j
and
Lh {v}(s) = sα − λ(1 − hs)k
−1
.
Let us show that v ∈ `1 if (7) has no root on cl D(h−1 , h−1 ) . Defining an auxiliary function v˜ : hZ+ → C as (−α−1)
v˜ = ph
we have
− λδ
where
δn = h−1
for
n = k + 1 and
δn = 0 otherwise,
−1 = 1. Lh {˜ v }(s) · Lh {v}(s) = sα − λ(1 − hs)k · sα − λ(1 − hs)k (−α−1)
al P
Moreover, v˜ ∈ `1 because of the asymptotic equivalence (ph )n ∼ Cn−α−1 as n → ∞ for a suitable C ∈ R (in (−α−1) fact, the exact sum of ph can be calculated). Finally, Lh {˜ v } is continuous and never zero on cl D(h−1 , h−1 ) since (7) has no root here by the assumption. Hence, by the Weierstrass theorem, we have also inf |Lh {˜ y }(s)| = min |Lh {˜ y }(s)| > 0 on cl D(h−1 , h−1 ) and Proposition 6 gives v ∈ `1 . The function u decays more slowly than v because of the presence of the term sα−1 (if 0 < α < 1) or sα−2 (ν−1) (ν−1) (if 1 < α < 2). Indeed, s−ν (0 < ν < 1) represents the Laplace image of ph for which (ph )n ∼ Cnν−1 (ν−1) 1 as n → ∞ for a suitable C ∈ R, hence, ph is not in ` . However, it is in a larger space, namely in `p where (ν−1) (1 − ν)−1 < p < ∞. Thus, by the Young’s inequality for convolution, ph ∗ v is also in `p and, consequently, it converges to zero. The remaining terms in Lh {u} correspond to summable functions which means that their convolution with v is again summable, thus converging to zero. Overall, we have proved that |yn | → 0 as n → ∞ provided (7) has no root on cl D(h−1 , h−1 ) . Note that we do −1 −1 not know stability behavior when (7) has a root on ∂D(h , h ) (in this case, the property (8) of Proposition 6 is not satisfied). Now we turn our attention to discussing conditions guaranteeing that characteristic equation (7) has a root on ∂D(h−1 , h−1 ). Let s be such a root, i.e., s = 2h−1 cos(ψ) exp(iψ)
for a suitable − π/2 < ψ ≤ π/2.
(9)
Jou rn
This s can equivalently be expressed as
s = h−1 1 + exp(i2ψ) .
(10)
Since the term (1 − hs)k is never zero for s ∈ ∂D(h−1 , h−1 ), the characteristic equation (7) can be rewritten as λ = g(s) = sα (1 − hs)−k .
(11)
Substituting (9) to g, we arrive (taking into account (10)) at g ∗ (ψ) = g(s)|s=2h−1 cos(ψ) exp(iψ) = (−1)k 2α h−α cosα (ψ) exp i(2k − α)ψ ,
−π/2 < ψ ≤ π/2.
Function g ∗ takes the values on a pair of curves lying in the complex plane that are symmetric with respect to the real axis. Clearly, the modulus |g ∗ (ψ)| is proportional to the value cosα (ψ) and thus is increasing when |ψ| moves from π/2 towards zero (because of cos is decreasing on [0, π/2]). It means that each of these two curves do not have self-intersections. Moreover, it crosses the real axis if ψ = jπ/(2k − α), j = 0, 1, . . . , k − 1 and if ψ = π/2. The crossing, which is nearest to the origin, occurs for j = k − 1 and lies on the negative real half-axis (since (−1)k exp(i(k − 1)π) = (−1)k (−1)k exp(−iπ) = −1). Hence, the restriction (k − 1)π/(2k − α) ≤ |ψ| ≤ π/2 causes that g ∗ just represents a simple closed curve in the complex plane. Substitution ψ = (kπ − |ϕ|)/(2k − α) 4
Journal Pre-proof
transforms g ∗ to g ∗∗ (ϕ) = 2α h−α cosα (kπ − |ϕ|)/(2k − α) exp(iϕ),
απ/2 ≤ |ϕ| ≤ π.
repro of
∗∗ The set Γα h,k = {g (ϕ) ∈ C : απ/2 ≤ |ϕ| ≤ π} splits (as a simple closed curve) the complex plane to two disjoint parts – an interior and an exterior. It can be shown that the preimage of the interior under the function g appearing in (11) is located outside of cl D(h−1 , h−1 ) and, conversely, the preimage of the exterior under g is just D(h−1 , h−1 ) \ {h−1 }. In the latter case, the exterior of Γα h,k is actually the codomain of g over D(h−1 , h−1 ) \ {h−1 }. The proof of this step requires all the parts of the function g ∗ (not only that one due to the restriction (k − 1)π/(2k − α) ≤ |ψ| ≤ π/2) and utilizes the open mapping theorem together with some continuity arguments (some technical refinements are needed as g is not holomorphic at h−1 and at the non-positive real half-axis), details how to proceed can be found in [17]. The above facts immediately imply that any choice of −1 λ belonging to the interior of Γα , h−1 ) , hence, (3) is asymptotically h,k causes that there is no root on cl D(h stable according to the previous part. Contrary, taking λ lying in the exterior of Γα h,k , we have a root of (7) located in D(h−1 , h−1 ) and thus (3) is unstable.
Remark 7. The stability condition in Theorem 3 is formulated as a sufficient one. To prove its necessity, one α has to dispose with the situation on the stability boundary ∂Sh,k which seems to be a much more complicated matter compared to the continuous counterpart (we remind that outside the closure of this area we have the instability). Both the exact and discrete stability sets are depicted in Figure 1. =(z) −
2α α (k − 1)π cos hα 2k − α
απ 2
α
απ 2
al P
π − απ/2 τ
απ 2
Σα τ
<(z)
=(z)
2α (k − 1)π cosα hα 2k − α
α Sh,k
Σα τ α Sh,k
−
−
−
−
π − απ/2 τ
α
<(z)
−
απ 2
α Figure 1: Exact and discrete asymptotic stability sets Σα τ and Sh,k . The figure corresponds to h = 0.5, k = 2 (hence τ = 1), α = 0.8 (left) and α = 1.1 (right), the scales of the axes are equal.
Remark 8. If we consider λ real and put h = 1 (the standard normalized discrete case), then the stability α property λ ∈ S1,k is reduced to the condition
Jou rn
−2α cosα (k − 1)π/(2k − α) < λ < 0
that provides a fractional generalization of the classical Levin–May stability condition −2 cos (k − 1)π)/(2k − 1) < λ < 0
for the delay difference equation
∇1 yn = λyn−k ,
n = 1, 2, . . .
(λ being a real scalar and k a positive integer) appearing in some problems of discrete population biology (see [18]). 4. Numerical stability of the backward Euler discretization Equation (1) is the simplest test equation that enables a generalization of the concept of A-stability to fractional delay differential equations. In the first-order delayed case, the notion of τ -stability was introduced in the following sense (see, e.g., [19]): A Θ-method for delay differential equations is called τ -stable if its appropriate difference scheme is asymptotically stable for any constant stepsize h = τ /k (k being a positive integer) whenever a test exact equation is asymptotically stable. Considering (1) with arbitrary 0 < α < 2 as a test exact equation, one can extend the notion of τ -stability for fractional delay differential equations as well. 5
Journal Pre-proof
Since (3) can be viewed as the backward Euler method (when θ = 1) applied to (1), we use our results to α + discuss its τ -stability. In other words, we need to verify the property Σα τ ⊂ Sh,k for any h = τ /k, k ∈ Z . Using the above explicit descriptions of corresponding stability regions, this property is equivalent to the inequality d(ϕ) = cos (kπ − ϕ)/(2k − α) − (ϕ − απ/2)/(2k) ≥ 0, απ/2 ≤ ϕ ≤ π (12)
=(z)
0.5 S1,1
p − 3π/4
repro of
that should be satisfied for any positive integer k. Nevertheless, this is not true. Indeed, if, e.g., α = 0.5 and 0.5 τ = h = k = 1, then Figure 2 clearly demonstrates that Σ0.5 is not a part of S1,1 (analytically, (12) does not 1 hold because of d(απ/2) = 0, d0 (απ/2) > 0 and d(π) < 0, hence, d(ϕ) ≥ 0 in a right neighborhood of απ/2 and d(ϕ) < 0 in a left neighborhood of π).
Σ0.5 1
π 4
<(z)
√ − 2
−
π 4
0.5 , the former is not a subset of the latter. Figure 2: The stability sets Σ0.5 and S1,1 1
5. Concluding remarks
al P
Consequently, the backward Euler method is not τ -stable. Note that this conclusion contradicts the situation in the first-order case when the backward Euler method embodies very strong stability properties (in fact, it is τ -stable even when two-term delay differential equation with complex coefficients is taken as a test equation, see [20]).
Jou rn
We have analyzed exact and discretized stability of the fractional delay differential equation (1). The stability region for its discretization (3) has been described explicitly and compared with the exact stability region. This comparison resulted into the conclusion that the backward Euler method for fractional delay differential equations is not τ -stable. The mathematical essence of this paper consists in two aspects. First, it is a derivation of the characteristic equation (7). More precisely, we have shown that (3) is asymptotically stable if its complex parameter λ lies outside the disk D(h−1 , h−1 ) (and (3) is not stable provided λ lies inside this disk). Till now, this theoretical condition has been known in the differential case only (see, e.g., [21] and [22]). Secondly, Theorem 3 provides an α effective reformulation of this condition and describes the appropriate discrete stability region Sh,k . From this point of view, Theorem 3 plays the same role as the corresponding effective conditions known from the differential case (see [11], [12] and partially also [8]). Regarding very few criteria on stability of fractional delay difference equations, some particular results on finite-time stability and Mittag-Leffler stability presented in [9] and [10] do not answer this fundamental stability problem discussed in the previous sections. Keeping in mind the numerical aspect of the studied problem, we note that there are some other papers suggesting suitable numerical schemes for solving fractional delay differential equations and answering typical questions such as the error estimate or the convergence rate (see, e.g., [23]). However, due to advanced form of the proposed methods, the analytical description of stability regions (similar to that performed in Section 3) is missing. Since the backward Euler scheme is, at least in the case of purely delay fractional differential equations, (µ) easy to implement (we note that the binomial coefficient in the function ph can easily be calculated recursively), the full description of its stability region is important. In particular, the inequality (12) can be useful when determining the maximal stepsize h such that the appropriate numerical solution can retain the key qualitative property (convergence to zero) of the exact solution. Possible extensions of Theorem 3 can be directed in several ways. Of course, one can consider (3) with other types of fractional operator. Some existing results in this area (see, e.g., [22]) indicate that the stability region derived for the Caputo and Riemann–Liouville fractional operators remains the same – up to the situation on the stability boundary where dynamics of the same equation with different fractional operators exhibits different stability properties. Further, (1) extended by a non-delayed term, i.e., the equation CD
α
y(t) = λy(t − τ ) + νy(t), 6
t ∈ (0, ∞)
Journal Pre-proof
with complex λ and ν, serves as a delay feedback control model enabling, among others, stabilization and synchronization of significant solutions of various differential systems (see, e.g., [4]). Thus, considering the corresponding Euler discretization (i.e., (3) extended by a non-delayed term), one can study and answer many issues connected with fractional control of discrete dynamical systems. Finally, another possible future research connected with the conclusions of this paper is related to stability analysis of more advanced numerical schemes for fractional delay differential equations.
repro of
Acknowledgement. The research has been supported by the grant 17-03224S of the Czech Science Foundation. The authors are grateful to the anonymous referees for valuable comments and suggestions that helped to improve the content of the paper. References
Jou rn
al P
[1] K. Diethelm, The Analysis of Fractional Differential Equations: An Application-Oriented Exposition Using Differential Operators of Caputo Type, Springer, Berlin, 2010. [2] Y. Zhou, Basic Theory of Fractional Differential Equations, World Scientific, 2014. [3] M. Lazarevi´ c, Finite time stability analysis of P Dα fractional control of robot time-delay systems. Mech. Res. Commun. 33 (2006) 269–279. [4] M. Lazarevi´ c, Stability and stabilization of fractional order time delay systems. Scientific Technical Review 61 (2011) 31–44. [5] B. Tao, M. Xiao, Q. Sun, J. Cao, Hopf bifurcation analysis of a delayed fractional-order genetic regulatory network model, Neurocomputing 275 (2018) 677–686. [6] S. Bhalekar, Stability analysis of a class of fractional delay differential equations, Pramana-J. Phys. 81 (2013) 215–224. [7] E. Kaslik, S. Sivasundaram, Analytical and numerical methods for the stability analysis of linear fractional delay differential equations, J. Comput. Appl. Math. 236 (2012) 4027–4041. [8] X. Teng, Z. Wang, Stability switches of a class of fractional-delay systems with delay-dependent coefficients, J. Comput. Nonlinear Dynam. 13 (2018) 111005 (9 pages). [9] F. Du, B. Jia, Finite-time stability of a class of nonlinear fractional delay difference systems, Appl. Math. Lett. 98 (2019) 233–239. [10] G.-C. Wu, D. Baleanu, L.L. Huang, Novel Mittag-Leffler stability of linear fractional delay difference equations with impulse, Appl. Math. Lett. 82 (2018) 71–78. [11] K. Krol, Asymptotic properties of fractional delay differential equations, Appl. Math. Comput. 218 (2011) 1515–1532. ˇ [12] J. Cerm´ ak, J. Horn´ıˇ cek, T. Kisela, Stability regions for fractional differential systems with a time delay, Commun. Nonlinear Sci. Numer. Simul. 31 (2016) 108–123. ˇ [13] J. Cerm´ ak, T. Kisela, Oscillatory and asymptotic properties of fractional delay differential equations, Electron. J. Differ. Eq. 2019 (2019) (15 pages). [14] F.M. Atici, P.W. Eloe, Linear systems of fractional nabla difference equations, Rocky Mountain J. Math. 41 (2011) 353–370. ˇ [15] J. Cerm´ ak, T. Kisela, L. Nechv´ atal, Stability regions for linear fractional differential systems and their discretizations, Appl. Math. Comput. 219 (2013) 7012–7022. [16] C.A. Desoer, M. Vidyassager, Feedback Systems: Input–Output Properties, Academic Press, New York, 1975. ˇ [17] J. Cerm´ ak, I. Gy¨ ori, L. Nechv´ atal, On explicit stability conditions for a linear fractional difference system, Fract. Calc. Appl. Anal. 18 (2015) 651–672. [18] S.A. Levin, R. May, A note on difference delay equations, Theor. Popul. Biol. 9 (1976) 178–187. [19] N. Guglielmi, Delay dependent stability regions of Θ-methods for delay differential equations, IMA J. Numer. Anal. 18 (1998) 399–418. [20] S. Maset, Stability of Runge Kutta methods for linear delay differential equations, Numer. Math. 87 (2000) 355–371. [21] C.P. Li, F.R. Zhang, A survey on the stability of fractional differential equations, Eur. Phys. J. Special Topics 193 (2011), 27–47. [22] D. Qian, Ch. Li, R. P. Agarwal, P. J. Y. Wong, Stability analysis of fractional differential system with Riemann-Liouville derivative, Math. Comput. Model. 52 (2010), 862–874. [23] A. Jhinga, V. Daftardar-Gejji, A new numerical method for solving fractional delay differential equations, Comput. Appl. Math. 38 (2019) 166 (18 pages).
7
Journal Pre-proof
Jou rn
al P
repro of
Jan Čermák: Conceptualization, Investigation, Methodology, Writing – original draft, Writing – review & editing. Luděk Nechvátal: Conceptualization, Investigation, Methodology, Writing – original draft, Writing – review & editing.