The Lambert function method in qualitative analysis of fractional delay differential equations

Fractional Calculus and Applied Analysis (2023) 26:1545–1565 https://doi.org/10.1007/s13540-023-00176-x ORIGINAL PAPER The Lambert function method in qualitative analysis of fractional delay differential equations Jan ˇ Cermák1 ·Tomáš Kisela1 ·Ludˇek Nechvátal1 Received: 13 December 2022 / Revised: 19 May 2023 / Accepted: 25 May 2023 / Published online: 16 June 2023 © The Author(s) 2023 Abstract We discuss an analytical method for qualitative investigations of linear fractional delay differential equations. This method originates from the Lambert function technique that is traditionally used in stability analysis of ordinary delay differential equations. Contrary to the existing results based on such a technique, we show that the method can result into fully explicit stability criteria for a linear fractional delay differential equation, supported by a precise description of its asymptotics. As a by-product of our investigations, we also state alternate proofs of some classical assertions that are given in a more lucid form compared to the existing proofs. Keywords Fractional delay differential equation (primary) ·Lambert function · Stability ·Asymptotic behavior Mathematics Subject Classification (Primary) 34K37 ·33E30 ·33E12 ·34K20 · 34K25 1 Introduction The paper discusses an analytical method for qualitative investigations of fractional delay differential equations (FDDEs). These equations are currently very intensively studied due to their importance in various application areas, with a special emphasis to control theory. Indeed, presence of both the time lag as well as non-integer derivative BLudˇek Nechvátal nechv[email protected].cz Jan ˇ Cermák [email protected].cz Tomáš Kisela [email protected].cz 1Institute of Mathematics, Brno University of Technology, Technická 2896/2, 61669 Brno, Czechia 123 1546 J. ˇ Cermák et al. order as control or tunning parameters in studied models provides a very efficient tool for various control processes such as stabilization or destabilization of the particular solutions of these models (for a pioneering work in this direction we refer to [17]). Systematic investigations of FDDEs were initiated in the paper [9]. Here, stability properties of Dαx(t)=λx(t−τ), (1) where α, τ ∈R+,λ∈Rand Dαis a fractional differential operator, were analyzed using the fact that (1) is asymptotically stable (i.e., any its solution is eventually tending to zero) if and only if all the roots of the characteristic equation sα−λexp(−sτ) =0(2) have negative real parts. To explore such a location of characteristic roots with respect to the imaginary axis, the Lambert function technique was utilized. The essence of the method consists in a representation formula for characteristic roots in terms of appropriate branches of this multi-valued function (for some precisions concerning the correct use of the Lambert function technique in stability analysis of (1), we refer also to [12]). A certain general disadvantage of this approach consists in its (seeming) disability to provide stability criteria in an explicit form depending on entry parameters only (i.e., on α,λand τin the case of (1)). As the other papers on stability and asymptotic properties of (1) followed, the Lambert function method was replaced by some alternate classical tools of stability investigations (such as D–partition method or τ-decomposition method) modified to the fractional case. Using these approaches, effective and non-improvable stability conditions for (1), supported by some asymptotic bounds, were derived in [16](the case λ∈R,0<α<1), [6] (the case λ∈C,0<α<1), and partially also in [7] (the case λ∈C,α>0). Some of the mentioned results can be extended also to the case of a two–term FDDE Dαx(t)=μx(t)+λx(t−τ). (3) In this respect, we refer to [2,5,15] (the case μ, λ ∈R,0<α<1) and [8](the case μ, λ ∈R,1<α<2). Following the integer-order case (see, e.g., [1]), (1) and (3) may serve as test equations for numerical analysis of FDDEs. From this point of view, it is very important to describe their basic qualitative properties in the strongest possible form. Then, when analyzing appropriate numerical schemes applied to these test equations, the ability to keep the key qualitative properties of the underlying exact equations is of basic importance. For some other recent advances in qualitative theory of FDDEs, we refer, e.g., to [3,10,11,18–20,23]. Following the above outlines, the aim of this paper is twofold. First, we deepen the existing knowledge on some qualitative properties of (1) with the Caputo fractional derivative. Second, perhaps a more important aspect of the paper consists in the way how we aim to do it. We come back to the Lambert function method used in [9] and show that this approach can offer more than formulae depending on the use of 123 The Lambert function method in qualitative... 1547 supporting software packages. In fact, this technique can result into actually effective stability and asymptotic criteria. The paper is organized as follows. Section 2recalls some existing findings on (1) and essentials of the Lambert function theory. In Sect. 3, we explore the Lambert function method in details. In particular, we give an alternate proof of the classical assertion saying that the characteristic root generated by the principal branch of the Lambert function has the largest real part, and formulate a criterion that enables to localize values of the principal branch in the complex plane. Section 4presents applications of these results to (1). Here, we extend the existing stability criteria for (1) to arbitrary (positive) real values of α, and formulate sharp asymptotic estimates for the solutions of (1). Some final remarks in Sect. 5conclude the paper. 2 Basic mathematical background In this section, we summarize some known facts relevant to our next investigations. First, we recall a close relationship between stability and asymptotic properties of (1), and distribution of the characteristic roots of (2). Then, we recall some basics of the Lambert function and its use in stability analysis of FDDEs. It was shown in [7] that any solution xof (1) with the Caputo fractional derivative (and a generally complex λ) can be written using the Mittag-Leffler type function Gλ,τ α,β (t)=t/τ−1  j=0 λj(t−jτ)αj+β−1 (α j+β) ,α,β>0,(4) where · denotes the upper integer part. More precisely, if φis a continuous initial (complex-valued) function on [−τ,0],φ0=φ(0)and φj,j=1,...,α−1, are (complex) constants (considered when α>1), then x(t)=α−1  j=0 φjGλ,τ α, j+1(t)+λ0 −τ Gλ,τ α,α(t−τ−ξ)φ(ξ)dξ(5) is the solution of (1) satisfying x(t)=φ(t)for all t∈[−τ,0], and limt→0+x(j)(t)= φj,j=1,...,α−1. Based on some asymptotic results on (4), the solution (5) can be rewritten by the use of the characteristic roots having non-negative real parts. We recall that (2) admits countably many roots, and only a finite number of them is lying right to any line (s)=p,p∈R(throughout the paper, the symbol (z)and (z)stands for the real and imaginary part of z∈C, respectively). If we denote by Sthe set of all roots of (2) having non-negative real parts (note that Smust be a finite set), then, for a non-integer α,(5) can be rewritten as x(t)= s∈S csexp(st)+O(tj−α)as t→∞ (6) 123 1548 J. ˇ Cermák et al. where csare complex coefficients depending on α,τ,λ,φ, and j∈{−1,0,...,α− 1}(the particular value of jdepends on limit behavior of φat t=0). Notice that j−α<0, i.e., the function tj−αalways tends to zero. By (6), the roots of (2) play an essential role in qualitative behavior of the solutions of (1). Following the classical integer-order pattern, the authors in [9]usedthe following chain of steps sαexp(sτ) =λ→sexpτ αs=λ1 α→τ αsexpτ αs=τ αλ1 α(7) to express the roots of (2) via the Lambert function introduced as the solution of W(z)exp(W(z)) =z,z∈C.(8) Before we recall the root formula for (2) based on this special function, some of its basic properties might be collected. The Lambert function is a multi-valued function (except at z=0) with infinitely many (single-valued) branches Wk,k∈Z. Neither of them can be expressed in terms of elementary functions. In particular, W0is called a principal branch. For any z∈C,(W0(z)) is between −πand π. The other branches are numbered so that (Wk(z)) is between (2k−2)π and (2k+1)π while (W−k(z)) is between −(2k+1)π and −(2k−2)π for any z∈Cand k=1,2,....More precisely, the ranges of W±kand W±(k+1),k=0,1,..., are separated by the curves {w=x+iy∈C:x=−ycot(y), 2kπ<|y|<(2k+1)π} and the ranges of W1and W−1are separated by the half-line {w=x+iy∈C:−∞<x≤−1,y=0}. These separating curves correspond to the branch cuts in the z-plane defined as {z=ξ+iη∈C:−∞<ξ≤−exp(−1), η =0} in the case of W0, and {z=ξ+iη∈C:−∞<ξ≤0,η=0} in the case of Wk,k= 0. Conventionally, the branch cut (having the argument π in the z-plane) is mapped by Wkon its upper boundary in the w-plane. Only the branches W0and W−1take on real values for a real z∈[−exp(−1), ∞)and a real z∈[−exp(−1), 0), respectively. Further details on the Lambert function (including some historical remarks) can be found in [4], for other comments, see also [13] and [22]. Now, following (7), all the roots of (2) can be expressed in the form sk=α τWkτ αλ1 α,k∈Z.(9) 123 The Lambert function method in qualitative... 1549 By (6), a crucial role in analysis of (1) is played by the rightmost characteristic root (i.e., the root of (2) with the largest real part). The following classical assertion says that this root is just s0. Lemma 1 Let z ∈C. Then W0(z)has the largest real part (W0(z)) among all the other real parts (Wk(z)),k∈Z. The original proof of Lemma 1is pretty long (see [22]). As a by-product of our next procedures, we are going to present an alternate (and more simple) way how to prove this assertion. Remark 1 As pointed out in [12], the expression (9) is not quite correct for some complex values of λ. More precisely, (7) contains taking the 1/α-power which means that the roots given by (9) are identical to those of (2) only in the case |Arg(λ)|≤απ (we recall that −π<Arg(·)≤π). This inequality is satisfied trivially when α≥1 but makes a restriction when 0 <α<1. In other words, if |Arg(λ)|>απ, then the representation (9) can produce some superfluous roots that are actually not the true roots of (2). As an example, we can consider, e.g., the case λ=−1, α=1/2 and τ=1 when (2) has the rightmost root s0≈−0.4172 −i2.2651 (i.e., (1) is asymptotically stable) while (9) produces s0≈0.4263 >0. On this account, we discuss qualitative properties of (1)forα>1. Comments to the case 0 <α<1 are provided in the final section. 3SomeadvancesontheLambertWfunction This section contains several key results on the Lambert function which proved to be useful in qualitative investigations of (1). To obtain an actually effective and strong asymptotic description of the solutions of (1), we need to effectively localize the position of the rightmost characteristic root in the complex plane. More precisely, by (6) and (9), we need to derive effective expressions of the real and imaginary parts of W0(z)in terms of z. Thus, keeping in mind intended stability and asymptotic analysis of (1), we can pose the following problems: For given p ∈Rand z ∈C,is it possible to characterize the properties (W0(z)) < p and (W0(z)) =p directly in terms of z and p,i.e., without an evaluation of the principal branch of the Lambert function? Further, for given q ∈Rand z ∈C,is it possible to similarly elaborate on the properties |(W0(z))|>q and |(W0(z))|=q? The following result yields an affirmative answer to these questions. Theorem 1 Let p,q∈R,p>−1,0<q<π, and z ∈C,z= 0. Then (i) (W0(z)) < p if and only if either |z|<pexp(p)or |z|≥|p|exp(p)and arccos pexp(p) |z|+|z|2−p2exp(2p) exp(p)<|Arg(z)|; (10) 123 1550 J. ˇ Cermák et al. (ii) (W0(z)) =p if and only if |z|≥|p|exp(p)and arccos pexp(p) |z|+|z|2−p2exp(2p) exp(p)=|Arg(z)|; (11) (iii) |(W0(z))|>q if and only if |Arg(z)|>q and q sin(|Arg(z)|−q)expqcot(|Arg(z)|−q)<|z|; (12) (iv) |(W0(z))|=q if and only if |Arg(z)|>q and q sin(|Arg(z)|−q)expqcot(|Arg(z)|−q)=|z|.(13) Proof (i) We write z=|z|exp(iArg(z)) and put xk=(Wk(z)),yk=(Wk(z)) where Wk,k∈Zare particular branches of the Lambert function. Substitution into (8) yields exp(xk)(xkcos(yk)−yksin(yk)) =|z|cos(Arg(z)), (14) exp(xk)(xksin(yk)+ykcos(yk)) =|z|sin(Arg(z)). (15) If we solve (14)–(15) with respect to unknowns xkexp(xk)and ykexp(xk), then xkexp(xk)=|z|cos(Arg(z)−yk), (16) ykexp(xk)=|z|sin(Arg(z)−yk). (17) To show that x0=(W0)(z)) < pwhenever |z|<pexp(p), we consider (16) implying x0exp(x0)≤|z|<pexp(p). Then the monotony property of the function g(p)=pexp(p)on (−1,∞)actually implies x0<p. Now we assume that |z|≥|p|exp(p). Squaring and adding (16) and (17) we get |z|2=((xk)2+(yk)2)exp(2xk), i.e., |yk|=|z|2−(xk)2exp(2xk) exp(xk).(18) 123 The Lambert function method in qualitative... 1551 For the principal branch, it holds x0≥−y0cot(y0), |y0|<π, i.e., x0sin(y0)+ y0cos(y0)≥0 whenever y0≥0. Multiplying this by exp(x0)and using (15), one gets |z|sin(Arg(z)) =exp(x0)(x0sin(y0)+y0cos(y0)) ≥0 which implies Arg(z)≥0fory0≥0. If y0<0, the same argumentation leads to Arg(z)≤0, hence Arg(z)y0≥0, i.e., |Arg(z)−y0|≤π. Then (16) with k=0is equivalent to arccos(x0exp(x0)/|z|)=|Arg(z)−y0|.(19) Moreover, sign analysis of (17) with respect to Arg(z)y0≥0 yields |Arg(z)|≥|y0|, i.e., |Arg(z)−y0|=|Arg(z)|−|y0|.(20) Then, using (18), (19) and (20), we are able to set up an implicit dependence between x0=(W0(z)) and zin the form f(x0,z)=0 where fis defined via f(p,z)=arccos pexp(p) |z|−|Arg(z)|+|z|2−p2exp(2p) exp(p) for all p>−1 and z∈Csuch that |p|exp(p)≤|z|.Letzbe fixed. Then df dp(p,z)=−(2p+1)exp(3p)+|z|2exp(p) exp(2p)|z|2−p2exp(2p)≤− (p+1)2exp(p) |z|2−p2exp(2p)≤0, hence, fis decreasing in pif |p|exp(p)≤|z|. Therefore, f(p,z)< f(x0,z)=0 whenever p>x0=(W0(z)) and |p|exp(p)≤|z|. (ii) The property follows directly from the proof of (i) using the fact that f(p,z)=0 if and only if p=x0due to monotony of fwith respect to p. (iii) Since W0is symmetric in the sense W0(z)=W0(z)for all z∈Cexcept those lying on the branch cut along the negative real axis between −∞ and −exp(−1),it suffices to assume the case y0=(W0(z)) > q>0. We have already observed that Arg(z)≥y0. In addition, a stronger property holds, namely Arg(z)>y0. Indeed, possible equality Arg(z)=y0implies y0=0 (due to (17)) which contradicts the assumption y0>q>0. Hence, it must be 0 <Arg(z)−y0<πas well as 0<Arg(z)−q<π.Wedivide(16)by(17) and put k=0 to get x0=y0cot(Arg(z)−y0). (21) 123 1552 J. ˇ Cermák et al. Taking the logarithm of (17) with k=0, we also have x0=ln|z|sin(Arg(z)−y0)−ln(y0). (22) Combining (21) and (22), we arrive at ln|z|sin(Arg(z)−y0)−ln(y0)−y0cot(Arg(z)−y0)=0 representing again an implicit dependence, now between (W0(z)) and z. If we denote h(q,z)=ln|z|sin(Arg(z)−q)−ln(q)−qcot(Arg(z)−q), then we have dh dq(q,z)=−qsin(2(Arg(z)−q)) −sin2(Arg(z)−q)−q2 qsin2(Arg(z)−q). While the denominator is positive, the numerator N(q,z)=−qsin(2(Arg(z)−q)) −sin2(Arg(z)−q)−q2 is negative for each 0 ≤q≤Arg(z). Indeed, we have dN dq(q,z)=2q(cos2(Arg(z)−q)−1)≤0 which implies that N(·,z)is non-increasing and, together with N(0,z)= −sin2(Arg(z)) < 0, negative on [0,Arg(z)]. Consequently, h(·,z)is decreasing and therefore h(q,z)>h(y0,z)=0 whenever 0 <q<y0<Arg(z). Taking into account the above mentioned symmetry, we arrive (after some elementary algebra) at (12). (iv) The required property is again a consequence of monotony of the function h from the previous part.  Remark 2 (a) The properties (ii) and (iv) of Theorem 1provide a new tool for evaluations of the principal branch of the Lambert function. Let z= 0 be a fixed complex number. Then the left-hand side of (11) is decreasing for all p∈[a,W0(|z|)](a=−1 if |z|≥exp(−1)and a=W0(−|z|)if |z|<exp(−1)) from πto the zero value. Hence, (11) has a unique root p∗lying in this interval, and this root equals just (W0(z)). Similarly, the left-hand side of (13) is increasing for all q∈(0,Arg(z)) from the zero value to infinity, i.e., (13) admits a unique positive root q∗which is just (W0(z)). To illustrate this evaluation technique, we compute W0(z)for z=1 2+i√3 2. Then |z|=1, Arg(z)=π/3 and the standard Newton method returns (z)=p∗≈0.4843 in 5 iterations with the initial value p0=0.5 and the stopping criterion taken as |pk+1−pk|≤10−16. The same method gives (z)=q∗≈0.3808 in 7 iterations with the initial value q0=0.5 and the same precision as in the case of the real part. 123 The Lambert function method in qualitative... 1553 In fact, the value p∗+iq∗matches the value produced by the MATLAB command lambertw(1/2+sqrt(3)/2*1i) to all the 15 digits behind the decimal point. Standardly, the Newton or Halley method is applied directly to the equation wexp(w)− z=0 using the complex arithmetic. MATLAB employs the latter method with some advanced guess of the starting point. For computing the values of the Lambert function with arbitrary precision, we refer to the recent paper [14]. (b) Using a different approach, the property (i) of Theorem 1was also discussed in [21]. In the sequel, we clarify ordering of the real as well as imaginary parts of the particular branches of the Lambert function. This ordering may be useful in a deeper asymptotic analysis of (1), and, moreover, results into an alternate proof of Lemma 1. Following the proof of Theorem 1, we introduce the functions Gz(x,y)=xsin(y)+ycos(y)−|z|sin(Arg(z)) exp(−x)and fz(x)=|z|2exp(−2x)−x2. In view of (15) and (18), the couples (xk,yk), where xk=(Wk(z)),yk=(Wk(z)), have to meet the relations Gz(x,y)=0 and y=±fz(x), respectively. The following assertion specifies ordering of imaginary parts of the branches of the Lambert function. Lemma 2 Let z ∈C\{0}. Then (Wk(z)) ≤(Wk+1(z)) for all k ∈Z. In fact, all the inequalities are strict with the only exception: If z ∈[−exp(−1), 0), then we have (W−1(z)) =(W0(z)) =0. Proof For the sake of formal simplicity, we identify complex numbers w=x+iy with couples (x,y)∈R2.First,letz∈C\{0}be such that 0 ≤Arg(z)≤πand define sets Sz j,j∈Z,as Sz j={(x,y)∈R2:Gz(x,y)=0,(2j−1)π < y<(2j+1)π}for j=1,2,...; Sz j={(x,y)∈R2:Gz(x,y)=0,0≤y<π}for j=0; Sz j={(x,y)∈R2:Gz(x,y)=0,−2π<y≤0}for j=−1; Sz j={(x,y)∈R2:Gz(x,y)=0,2jπ<y<(2j+2)π}for j=−2,−3,... (note that the equation Gz(x,y)=0 has no solution for y=(2j−1)π,j=1,2,..., and for y=2jπ,j=−1,−2,...). We wish to show that Sz jis a part of the range of Wkjust when j=k. Let k≥1 be arbitrary. Then, by the definition of Wk(see also Sect. 2), (2k−2)π < yk<(2k+1)π. (23) Let jbe such that (xk,yk)∈Sz j. We distinguish the following cases with respect to j. 123 1560 J. ˇ Cermák et al. as t→∞, where c=cs0is the complex constant from (6) corresponding to the rightmost characteristic root s0.Ifαis an integer, then the dominating role of s0in asymptotic behavior of (1) is well known. In this case, the assertion of (ii) holds as well.  Remark 6 (a) The asymptotic formula from Theorem 3(ii) immediately implies x(t)=O(exp(u0t)) as t→∞ (32) for any solution xof (1), and the constant u0is non-improvable. Moreover, for large t, the roots of the real and imaginary parts of xtend to the roots of cos(v0t)and sin(v0t), respectively. In both the cases, the distance between the subsequent roots tends to π/v0. These properties are illustrated by Example 1. (b) The asymptotic behavior of (1) significantly depends on stability of (1). In particular, the exponential terms in (6) are vanishing in the asymptotically stable case (s0)<0. However, the situation changes in the limit case α=1 when, in accordance with the first-order theory, the rightmost characteristic root s0determines an exponential decay rate of the solutions also in the asymptotically stable case. Since the above argumentation can be extended to this problem as well, our results provide a contribution also to the corresponding classical first-order theory. Example 1 Let α=1.2, τ=1, and consider (1) along with the initial conditions φ(t)=1(−1≤t≤0), φ0=φ(0)=1, and φ1=limt→0+x(t)=0. We compare the corresponding (numerical) solutions of (1) for two distinct values of λ, namely λ1=−2+i and λ2=−3+i0.1. As indicated by Fig. 2, both the values λ1,λ2lie in the instability region. In particular, the real parts u0of the corresponding rightmost roots are approximately 0.4721 and 0.4917, and their imaginary parts v0are 1.2321 and 1.5844, respectively. The real parts of the solutions of (1) with two above specified sets of entries, along with the growth-rate functions exp(u0t), are depicted in Figs. 3and 4. The graphs suggest that the modulus of constant cintroduced in Theorem 3(ii) is less than one for λ=λ1, and greater than one for λ=λ2. To illustrate behavior of the solutions xin better detail, Figs. 5and 6depict the ratio (x(t))/ exp(u0t)for λ1and λ2, respectively. The resulting functions are bounded, but do not tend to zero which is a consequence of non-improvability of the constant u0in (32). As mentioned in Remark 6(a), the distance between the subsequent roots of (x(t)) tends to π/v0.Figs.7and 8illustrate this fact. We can see that while in the case of λ1 the convergence is rather fast and the distance seems to be somewhat stabilized around the seventh root, in the case of λ2, the stabilization occurs around the hundredth root. 5 Concluding remarks The aim of the paper was to develop the Lambert function theory, and then apply the obtained results in qualitative investigations of (1). Using this approach, we were able 123 The Lambert function method in qualitative... 1561 Fig. 3 The real part of the solution xof (1)forα=1.2, τ=1andλ1=−2+i, along with the corresponding growth-rate functions ±exp(0.4721t) Fig. 4 The real part of the solution xof (1)forα=1.2, τ=1andλ2=−3+i0.1, along with the corresponding growth-rate functions ±exp(0.4917t) to formulate a precise asymptotic description of the solutions of (1). Particularly, in addition to an algebraic decay rate of the solutions in the stable case (described in some earlier papers), we could observe an exponential growth of the solutions in the 123 1562 J. ˇ Cermák et al. Fig. 5 The real part of the solution xof (1)forα=1.2, τ=1andλ1=−2+i, divided by its growth-rate function exp(0.4721t) Fig. 6 The real part of the solution xof (1)forα=1.2, τ=1andλ2=−3+i0.1, divided by its growth-rate function exp(0.4917t) unstable case; the rate of this growth was determined as a (unique) real root of an auxiliary transcendental equation. 123 The Lambert function method in qualitative... 1563 Fig. 7 The distance between the subsequent roots of (x(t)) for α=1.2, τ=1andλ1=−2+i is tending to π/1.2321 Fig. 8 The distance between the subsequent roots of (x(t)) for α=1.2, τ=1andλ2=−3+i0.1is tending to π/1.5844 However, the impact of the presented results is not limited to the theory of FDDEs only. Our approach offers an alternate way how to prove (and also strengthen) some classical assertions of the Lambert function theory. Moreover, to the best of our knowledge, the derived asymptotic formulae are new also in the first-order case. Here, 123 1564 J. ˇ Cermák et al. contrary to the fractional case, our results can be applied also in the stable case where a (non-improvable) rate of exponential decay of the solutions can be determined. Since we have formulated our results for (1) with a complex coefficient λ, their extension to the vector case is nearly straightforward provided the eigenvalues of a (real) system matrix are simple. Regarding eigenvalues with higher multiplicities, some additional argumentation seems to be necessary. Based on related cases discussed in earlier papers, one can expect a slight modification of the solutions growth, but no impact on the asymptotic frequency. Our final remark concerns the case 0 <α<1 not involved among the assumptions of the assertions of Sect. 4. The procedure of computing the characteristic roots uses the law of exponents which is, in general, not valid for complex numbers. Thus, some superfluous roots of the characteristic equation may appear if 0 <α<1(as illustrated via a counterexample in Remark 1). In this case, our stability and asymptotic formulae remain basically true, but we cannot confirm their strictness. In particular, we cannot claim that the above described rate of exponential growth of solutions is non-improvable. Nevertheless, we conjecture that a more thorough analysis of the corresponding branches of a complex power can overcome this problem, and thus achieve the strict asymptotic results for all α>0. Such an analysis provides another possible topic for the next research. Acknowledgements The research has been supported by the Grant GA20-11846S of the Czech Science Foundation. Funding Open access publishing supported by the National Technical Library in Prague. Declarations Conflict of interest The authors declare that they have no conflict of interest. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. References 1. Bellen, A., Zennaro, M.: Numerical Methods For Delay Differential Equations. The Clarendon Press, Oxford University Press, New York, Numerical Mathematics and Scientific Computation (2003) 2. Bhalekar, S.: Stability analysis of a class of fractional delay differential equations. Pramana-J. Phys. 81(2), 215–224 (2013). https://doi.org/10.1007/s12043-013-0569-5 3. Bhalekar, S.: Stability and bifurcation analysis of a generalized scalar delay differential equation. Chaos 26, Article ID 084306, 7 pp. (2016). https://doi.org/10.1063/1.4958923 4. Corless, R.M., Gonnet, G.H., Hare, D.E.G., Jeffrey, D.J., Knuth, D.E.: On the Lambert Wfunction. Adv. Comput. Math. 5(4), 329–359 (1996). https://doi.org/10.1007/BF02124750 123 The Lambert function method in qualitative... 1565 5. ˇ Cermák, J., Došlá, Z., Kisela, T.: Fractional differential equations with a constant delay: Stability and asymptotics of solutions. Appl. Math. Comput. 298, 336–350 (2017). https://doi.org/10.1016/j.amc. 2016.11.016 6. ˇ Cermák, J., Horníˇcek, J., Kisela, T.: Stability regions for fractional differential systems with a time delay. Commun. Nonlinear Sci. Numer. Simul. 31(1–3), 108–123 (2016). https://doi.org/10.1016/j. cnsns.2015.07.008 7. ˇ Cermák, J., Kisela, T.: Oscillatory and asymptotic properties of fractional delay differential equations. Electron. J. Differ. Equ. 2019, Paper No. 33, 15 pp. (2019) 8. ˇ Cermák, J., Kisela, T.: Stabilization and destabilization of fractional oscillators via a delayed feedback control. Commun. Nonlinear Sci. Numer. Simul. 117, Article ID 106960, 16 pp. (2023). https://doi. org/10.1016/j.cnsns.2022.106960 9. Chen, Y., Moore, K.L.: Analytical stability bound for a class of delayed fractional-order dynamic systems. Nonlinear Dyn. 29, 191–200 (2002) 10. Daftardar-Gejji, V., Sukale, Y., Bhalekar, S.: Solving fractional delay differential equations: A new approach. Fract. Calc. Appl. Anal. 18, 400–418 (2015). https://doi.org/10.1515/fca-2015-0026 11. Garrappa, R., Kaslik, E.: On initial conditions for fractional delay differential equations. Commun. Nonlinear Sci. Numer. Simul. 90, Article ID 105359, 16 pp. (2020). https://doi.org/10.1016/j.cnsns. 2020.105359 12. Hwang, C., Cheng, Y.C.: A note on the use of the Lambert Wfunction in the stability analysis of timedelay systems. Automatica 41(11), 1979–1985 (2005). https://doi.org/10.1016/j.automatica.2005.05. 020 13. Jeffrey, D.J., Hare, D.E.G., Corless, R.M.: Unwinding the branches of the Lambert Wfunction. Math. Sci. 21(1), 1–7 (1996) 14. Johansson, F.: Computing the Lambert Wfunction in arbitrary-precision complex interval arithmetic. Numer. Algorithms 83, 221–242 (2020). https://doi.org/10.1007/s11075-019-00678-x 15. Kaslik, E., Sivasundaram, S.: Analytical and numerical methods for the stability analysis of linear fractional delay differential equations. J. Comput. Appl. Math. 236(16), 4027–4041 (2012). https:// doi.org/10.1016/j.cam.2012.03.010 16. Krol, K.: Asymptotic properties of fractional delay differential equations. Appl. Math. Comput. 218(5), 1515–1532 (2011). https://doi.org/10.1016/j.amc.2011.04.059 17. Lazarevi´c, M.P.: Finite time stability analysis of PDαfractional control of robotic time-delay systems. Mech. Res. Commun. 33, 269–279 (2006). https://doi.org/10.1016/j.mechrescom.2005.08.010 18. Liu, L., Dong, Q., Li, G.: Exact solutions of fractional oscillation systems with pure delay. Fract. Calc. Appl. Anal. 25, 1688–1712 (2022). https://doi.org/10.1007/s13540-022-00062-y 19. Li, M., Wang, J.R.: Finite time stability of fractional delay differential equations. Appl. Math. Lett. 64, 170–176 (2017). https://doi.org/10.1016/j.aml.2016.09.004 20. Medved’, M., Pospíšil, M.: On the existence and exponential stability for differential equations with multiple constant delays and nonlinearity depending on fractional substantial integrals. Electron. J. Qual. Theory Differ. Equ. 2019, Paper No. 43, 17 pp. (2019) 21. Nishiguchi, J.: On parameter dependence of exponential stability of equilibrium solutions in differential equations with a single constant delay. Discrete Contin. Dyn. Syst. 36(10), 5657–5679 (2016). https:// doi.org/10.3934/dcds.2016048 22. Shinozaki, H., Mori, T.: Robust stability analysis of linear time-delay systems by Lambert Wfunction: Some extreme point results. Automatica 42(10), 1791–1799 (2006). https://doi.org/10.1016/j. automatica.2006.05.008 23. Tuan, T.H., Trinh, H.: A linearized stability theorem for nonlinear delay fractional differential equations. IEEE Trans. Autom. Control 63(9), 3180–3186 (2018). https://doi.org/10.1109/TAC.2018.2791485 Publisher’s Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations. 123