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

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.


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 B Luděk Nechvátal nechvatal@fme.vutbr.cz JanČermák cermak.j@fme.vutbr.cz Tomáš Kisela kisela@fme.vutbr.cz 1 Institute of Mathematics, Brno University of Technology, Technická 2896/2, 616 69 Brno, Czechia 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 where α, τ ∈ R + , λ ∈ R and 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 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 − τ ). ( 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][19][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 supporting software packages. In fact, this technique can result into actually effective stability and asymptotic criteria.
The paper is organized as follows. Section 2 recalls 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 4 presents 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. 5 conclude the paper.

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 x of (1) with the Caputo fractional derivative (and a generally complex λ) can be written using the Mittag-Leffler type function 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 is the solution of (1) satisfying 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 S the set of all roots of (2) having non-negative real parts (note that S must be a finite set), then, for a non-integer α, (5) can be rewritten as where c s are complex coefficients depending on α, τ , λ, φ, and j ∈ {−1, 0, . . . , α − 1} (the particular value of j depends on limit behavior of φ at t = 0). Notice that j − α < 0, i.e., the function t j−α 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] used the following chain of steps to express the roots of (2) via the Lambert function introduced as the solution of 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 W k , k ∈ Z. Neither of them can be expressed in terms of elementary functions. In particular, W 0 is called a principal branch. For any z ∈ C, (W 0 (z)) is between −π and π . The other branches are numbered so that (W k (z)) is between (2k − 2)π and (2k + 1)π while (W −k (z)) is between −(2k + 1)π and −(2k − 2)π for any z ∈ C and k = 1, 2, . . .. More precisely, the ranges of W ±k and W ±(k+1) , k = 0, 1, . . . , are separated by the curves and the ranges of W 1 and W −1 are separated by the half-line These separating curves correspond to the branch cuts in the z-plane defined as in the case of W 0 , and in the case of W k , k = 0. Conventionally, the branch cut (having the argument π in the z-plane) is mapped by W k on its upper boundary in the w-plane. Only the branches W 0 and W −1 take 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 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 s 0 .

Lemma 1
Let z ∈ C. Then W 0 (z) has the largest real part (W 0 (z)) among all the other real parts (W k (z)), k ∈ Z.
The original proof of Lemma 1 is 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.

Some advances on the Lambert W function
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 W 0 (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 ∈ R and z ∈ C, is it possible to characterize the properties (W 0 (z)) < p and (W 0 (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 ∈ R and z ∈ C, is it possible to similarly elaborate on the properties | (W 0 (z))| > q and | (W 0 (z))| = q? The following result yields an affirmative answer to these questions.
(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 1 provide 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, W 0 (|z|)] (a = −1 if |z| ≥ exp(−1) and a = W 0 (−|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 (W 0 (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 (W 0 (z)).
To illustrate this evaluation technique, we compute W 0 (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 p 0 = 0.5 and the stopping criterion taken as | p k+1 − p k | ≤ 10 −16 . The same method gives (z) = q * ≈ 0.3808 in 7 iterations with the initial value q 0 = 0.5 and the same precision as in the case of the real part.
In fact, the value p * + i q * 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 w exp(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 1 was 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 G z (x, y) = x sin(y) + y cos(y) − |z| sin(Arg(z)) exp(−x) and In view of (15) and (18), the couples (x k , y k ), where x k = (W k (z)), y k = (W k (z)), have to meet the relations G z (x, y) = 0 and y = ± f z (x), respectively. The following assertion specifies ordering of imaginary parts of the branches of the Lambert function. Proof For the sake of formal simplicity, we identify complex numbers w = x + i y with couples (x, y) ∈ R 2 . First, let z ∈ C \ {0} be such that 0 ≤ Arg(z) ≤ π and define sets S z j , j ∈ Z, as (note that the equation G z (x, y) = 0 has no solution for y = (2 j − 1)π , j = 1, 2, . . . , and for y = 2 jπ , j = −1, −2, . . . ). We wish to show that S z j is a part of the range of W k just when j = k.
Let k ≥ 1 be arbitrary. Then, by the definition of W k (see also Sect. 2), Let j be such that (x k , y k ) ∈ S z j . We distinguish the following cases with respect to j.
If j > k, then y k > (2k + 1)π which contradicts (23). Let j = k. The ranges of W k and W k+1 , k = 0, 1, . . . , are separated by the curve provided y = 2kπ . To estimate the y-coordinate of a point (x, y) ∈ S z k , we put u = x in γ k .
First, let 2kπ < y < (2k + 1)π . Then any point (x, y) of S z k together with the corresponding point (x, v) of γ k have to fulfill the formula due to x = −v cot(v) and (24). Since 0 ≤ Arg(z) ≤ π , the right-hand side of (25) is non-negative, hence, we have y cot(y) ≥ v cot(v) implying y ≤ v. In other words, any point (x, y) ∈ S z k is located below or on the curve γ k separating W k and W k+1 . Second, let (2k − 1)π < y ≤ 2kπ . Then the points (x, y) ∈ S z k lie below γ k trivially (note that the equation G z (x, y) = 0 has a unique solution x for y = 2kπ , 0 < Arg(z) < π, and has no solution for y = 2kπ , Arg(z) = 0 or Arg(z) = π ). On the other hand, any point of S z k lies above the upper bound (2k − 1)π of γ k−1 , hence, we have proven that S z k is contained in the range of W k . Finally, if j < k, then we can similarly verify that any point of S z j is already located below or on γ k−1 . Thus, to summarize the previous observations, S z j is contained in the range of W k (k = 1, 2, . . . ) just when j = k; otherwise, S z j and the range of W k are disjoint. Consequently, (2k − 1)π < y k < (2k + 1)π , hence, y k < y k+1 for all k = 1, 2, . . ..
The procedure can be analogously applied to the case −π < Arg(z) < 0 (definition of the sets S z j now differ by shifting the particular y-domains vertically down by π ).
Lemma 2 is useful also for ordering of the real parts of the branches of the Lambert function. In particular, it enables to prove the classical assertion of Lemma 1 in a more lucid way compared to the existing proof techniques.
Further, we have Thus, f z is decreasing on its domain up to "a small part" which is, however, lying within the range of W 0 (we again identify w = x + i y with (x, y) ∈ R 2 ). Indeed, if Similarly, we can show that, for exp(−1) ≥ |z| > 0, the graph of f z between the relevant points is contained in the range of W 0 as well.

Remark 3
We have actually proven monotony of the real parts of W k (z) with respect to k which is a slightly stronger result than that stated in Lemma 1. dependence on the parameter λ. We remind that here we restrict ourselves to the case α > 1 due to the reason mentioned in Remark 1.
We start with a basic stability criterion for (1). We introduce a (parametric) set α,τ 0 given as see Fig. 1. This set already appeared in [6,Thm. 2] as the asymptotic stability region for (1) with 0 < α < 1. As indicated in [7], the D-subdivision method (combined with some tools of fractional calculus) used in [6] is extendable also to the case α > 1. In the following assertion, we confirm validity of this stability result for α > 1. Contrary to the existing techniques, we are able to prove this result (as a consequence of Theorem 1(i)) in an almost elementary way.
Theorem 1 can be used in a more subtle way to bring a deeper insight into behavior of (1). More precisely, relations (10)-(13) enable to reveal a relationship between the position of the rightmost characteristic root s 0 and the value of λ. Then, by (6), such a relationship easily results into a precise asymptotic description of the solutions of (1) in the unstable case ( (s 0 ) > 0). Indeed, while in the asymptotically stable case ( (s 0 ) < 0) the decay rate of solutions is algebraic and independent of the particular value of s 0 , in the unstable case ( (s 0 ) > 0), the growth rate is exponential and governed by the real part of the rightmost characteristic root s 0 (it is well known that this is true also for integer values of α). In addition, the imaginary part of s 0 is related to the frequency characteristics describing an oscillatory behavior of (1).
Thus, for given u, v ≥ 0, we need to find a region of all complex λ such that the rightmost characteristic root s 0 is lying on or left to the line u + i ω, ω ∈ R, and on or above (below) the line ω + i v (the line ω − i v), ω ∈ R. On this account, similarly as in the proof of Theorem 2, we put z = τ λ 1/α /α, p = τ u/α, q = τ v/α and consider u ≥ 0, απ/τ > v > 0. Then, Theorem 1 immediately implies that or |λ| ≥ u α exp(τ u) and α arccos u exp(τ u/α) Thus, if we introduce the set α,τ u as a set of all λ ∈ C such that either (27) or (28) holds, and α,τ v as a set of all λ ∈ C such that (29) holds (we put α,τ v = ∅ for απ/τ ≤ v < π and α,τ v = C for v = 0), then we can rewrite our previous observations into the following assertion: Lemma 3 Let α > 1, τ > 0, u ≥ 0, π > v ≥ 0, λ ∈ C and let s 0 be given by (9). Then
(b) It is easy to check that α,τ 0 introduced in (26) coincides with int α,τ 0 . The part (i) of Lemma 3 immediately implies Corollary 1 Let u ≥ 0 be fixed. Then α,τ u is the set of all λ ∈ C such that x(t) = O(exp(ut)) as t → ∞ for any solution x of (1).
To obtain an actually effective (and non-improvable) asymptotic result for the solutions of (1), we have to look at the problem inversely. More precisely, for a given complex λ / ∈ α,τ 0 , we need to find (non-negative) real values u 0 , v 0 such that the rightmost root s 0 of (2) satisfies The way how to do it easily follows from Theorems 1(ii) and 1(iv), respectively, taking into account the above introduced substitutions z = τ λ 1/α /α, p = τ u/α, q = τ v/α. Then the corresponding relations (11) and (13) become |λ| ≥ u α exp(τ u) and α arccos u exp(τ u/α) and |Arg(λ)| > τv/α and v α respectively (see also (28) and (29)). Thus, the unique solution of (30) 2 defines the value u 0 , while the unique positive solution of (31) 2 defines the value v 0 (provided Arg(λ) > 0). If Arg(λ) = 0, then s 0 is real, and we set v 0 = 0. Notice that (30) 2 , defining the relation between the modulus and argument of λ that is explicit with respect to the argument, forms the boundary of α,τ u . Its lefthand side, considered as a function of |λ|, is continuous, increasing and unbounded on [u α exp(τ u), ∞), and its graph in the complex plane creates a Jordan curve symmetric with respect to the real axis, see Fig. 2.
Similarly, (31) 2 provides the relation between the modulus and argument of λ that is explicit with respect to the modulus. The equality (31) 2 is the boundary of α,τ v , and its left-hand side, as a function of |Arg(λ)|, is continuous on (τ v/α, π ] and unbounded in a right neighborhood of the point τ v/α (for |Arg(λ)| = π , it takes a value on the negative real axis). This implies that α,τ v is unbounded for any 0 < v < απ/τ and its boundary splits the complex plane into two parts, see Fig. 2. Now we are in a position to formulate a complete asymptotic description for the solutions of (1). and v, respectively (the scenario corresponds to α = 1.2 and τ = 1). The particular blue curves represent the set of all λ ∈ C such that the rightmost characteristic root s 0 of (2) satisfies (s 0 ) = u, and the particular orange curves represent the set of all λ ∈ C such that the rightmost characteristic root s 0 of (2) satisfies | (s 0 )| = v. As an example, the blueish curvilinear rectangles then represent the set of all λ ∈ C such that 0.5 < (s 0 ) < 0.75 and 1.25 < | (s 0 )| < 1.5 Theorem 3 Let α > 1, τ > 0 and λ ∈ C.
(i) If λ ∈ α,τ 0 , then, for any solution x of (1), Moreover the algebraic decay order 1 − α cannot be improved; (ii) If λ / ∈ α,τ 0 , then, for any solution x of (1), Proof (i) The property is a direct consequence of (6) as the set S is empty.
(ii) Let s 0 = u 0 + i v 0 be the rightmost characteristic root, and S 0 be the set of the remaining characteristic roots s with a non-negative real part (we remind that (s) < (s 0 ) for any s ∈ S 0 ). Then, if α is a non-integer, we can write (6) as where c = c s 0 is the complex constant from (6) corresponding to the rightmost characteristic root s 0 . If α is an integer, then the dominating role of s 0 in 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 for any solution x of (1), and the constant u 0 is non-improvable. Moreover, for large t, the roots of the real and imaginary parts of x tend to the roots of cos(v 0 t) and sin(v 0 t), respectively. In both the cases, the distance between the subsequent roots tends to π/v 0 . 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 (s 0 ) < 0. However, the situation changes in the limit case α = 1 when, in accordance with the first-order theory, the rightmost characteristic root s 0 determines 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.
The real parts of the solutions of (1) with two above specified sets of entries, along with the growth-rate functions exp(u 0 t), are depicted in Figs. 3 and 4. The graphs suggest that the modulus of constant c introduced in Theorem 3(ii) is less than one for λ = λ 1 , and greater than one for λ = λ 2 .
To illustrate behavior of the solutions x in better detail, Figs. 5 and 6 depict the ratio (x(t))/ exp(u 0 t) for λ 1 and λ 2 , respectively. The resulting functions are bounded, but do not tend to zero which is a consequence of non-improvability of the constant u 0 in (32).
As mentioned in Remark 6(a), the distance between the subsequent roots of (x(t)) tends to π/v 0 . Figs. 7 and 8 illustrate 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. Fig. 3 The real part of the solution x of (1) for α = 1.2, τ = 1 and λ 1 = −2 + i, along with the corresponding growth-rate functions ± exp(0.4721t) Fig. 4 The real part of the solution x of (1) for α = 1.2, τ = 1 and λ 2 = −3 + i 0.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 Fig. 5 The real part of the solution x of (1) for α = 1.2, τ = 1 and λ 1 = −2 + i, divided by its growth-rate function exp(0.4721t) Fig. 6 The real part of the solution x of (1) for α = 1.2, τ = 1 and λ 2 = −3 + i 0.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. Fig. 7 The distance between the subsequent roots of (x(t)) for α = 1.2, τ = 1 and λ 1 = −2 + i is tending to π/1.2321 Fig. 8 The distance between the subsequent roots of (x(t)) for α = 1.2, τ = 1 and λ 2 = −3 + i 0.1 is 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, 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.