APPLICATION OF LAMBERT W FUNCTION IN OSCILLATION THEORY

This paper is primarily concerned with analyzing some special properties of the characteristic equation of the linear delay differential equation of the first and odd order. The basic concepts and results keen to the oscillatory nature of solutions, as well as stability analysis are re-investigated using the properties of Lambert W function and compared with existing ones. A brief survey of the function origin and possible applications is presented together with a short behavior analysis.


INTRODUCTION
Best known today for proving the irrationality of π and introducing the concept of hyperbolic functions in trigonometry, the 18th-century Swiss mathematician, physicist, philosopher and astronomer Johann Heinrich Lambert considered in 1758 the series solution to the symmetric case of the trinomial equation x a −x b = (a−b)vx a+b .Six years later, the subsequent work of Euler on a more special case, which reduces to the form wa w = ax, resulted in some of the mathematical properties of this equation, where Lambert's initial contribution was explicitly emphasized.However, the importance and great applicability of this formulation was not appreciated and obtained no special notation until the last decade of the 20th century, when it was included in the computer algebra software Maple.
An influential work was done by Corless et al., as a result of a systematic literature search of the simplest nontrivial transcendental function in many isolated contributions collected in a compact overview [5], it was shown that the Lambert W function had been re-discovered at various moments in history.Moreover, a new discussion of complex branches, symbolic calculus and mainly the accurate implementation of the function has encouraged many researchers.Due to its simplicity, the wide range of applications associated with the Lambert W function arises in the modeling of diverse phenomena in various fields of science and engineering.In physics, it is utilized for the purpose of quantum theory, plasma physics and solar physics (see e.g.[8,17,18]).Applications in the field of theoretical computer science include the analysis of algorithms and counting search trees and graphs occuring in combinatorial applications as well as the iterated exponentiation [13].
From a strictly mathematical point of view, the most significant use of the Lambert W function is related to finding zeros of transcendental functions.Hence, without awareness of the Lambert W function, the location of zeros for such functions was determined either numerically or aproximately [12].When applying elementary analytical techniques, an explicit solution to various exponential or logarithmic polynomials can be formulated in terms of the Lambert W function, especially in the field of time-delay systems [1,14,[19][20][21][22].Moreover, an efficient implementation of this special function is available in various mathematicalengineering software tools.

LAMBERT W FUNCTION
We begin here by describing some basic notation and facts about the mentioned special function which we will use in later section.
The Lambert W function, henceforth also denoted by W (z), can be defined in analogy with the complex logarithm as the multivalued inverse of the function z = f (w) = we w , where w and z are assumed to be complex variables, satisfying the so-called defining equation z = W (z)e W (z) .Thus, the mentioned special function, also known as a Product Log function, which clearly indicates an extension of a natural logarithm, is an example of multivalued functions f (z), which are characterized, unlike single-valued-ones, by two or more distinct values for each value of z.The multivaluedness is represented by infinitely many branches indexed by k ∈ Z as W k 's, which restrict the original domain to a smaller subsets of the complex plane, on which our function is already single-valued and continuous.Using the Lagrange's Inverse Theorem, the principal branch (i.e., k = 0) of the Lambert W function can be represented by the following power series [6] with z ∈ C, while the other branches are defined in a series form given below: where ln k (z) represents the kth logarithm branch and the coefficients C lm can be expressed in terms of nonnegative Stirling numbers of first kind.Complete and detailed derivation is provided by Corless et al. in [5].
The branches of W k (z) are defined in detail in the complex plane (i.e., Re[W k (z)] → Im[W k (z)]), however, our interest is primarily focused on the real and imaginary parts of Recall the definition of the Lambert W function and denote their variables: z = f (w) = we w , w = f −1 (z) = W (z), where w = x + iy, z = u + iv.We have z = (x + iy)e x (cos y + i sin y).
(1) Separating real and imaginary parts of z, the following equations are obtained: For further needs, we restrict our interest to a case where v = 0 (it follows from the fact that the arguments of the Lambert W function used in the next section are always real), therefore ye x cos y + xe x sin y = 0, which implies x = −y cotgy or y = 0.Both cases will be discussed. - A real-valued function f (x) is strictly decreasing on (−∞, −1), strictly increasing on (−1, ∞) and has a global minimum point at x = −1 of z = −1/e.A second derivative indicates an inflection point at x = −2 of z = −2e −2 .The function f (x) is not injective, hence the inverse function cannot be constructed on the whole interval.For this reason, when restricted on each of monotone intervals, two real branches of the Lambert W function are defined as follows: for z ∈ [−1/e, 0): the principal branch satisfying W (z) ≥ −1 and denoted by W 0 (x) is real and increasing.Similarly, the branch satisfying W −1 < −1 and denoted by W −1 (z) is real and decreasing from −∞ to −1, where it meets with W 0 (z) at the so-called branch point.for z ∈ [0, ∞), the only part on the principal branch is real, which implies exactly one real root of the defining equation, which is positive except for W (0) = 0. with z < −1/e, every W k (z) is complex-valued.
For our future needs, only the principal branch will be assumed.
The properties are shown in Figure 1.Similarly, the −1th branch, which is real-valued only for z ∈ [−1/e, 0] is illustrated (see Figure 2).The important values are highlighted by vertical dashed lines passing throught the x-axis at −1/e and −π/2.

DELAY DIFFERENTIAL EQUATIONS
Since a result of any change of the system state is not instantaneous, a delay is proposed to play an essential role in modeling in order to represent basically a time taken to complete some hidden processes.The related theoretical framework of delay differential equations as a special subclass of functional differential equations has received a great deal of attention since the first half of the 1950s, and has found applications in various fields of science and engineering (see e.g.[3,7,11,15] and references therein).Such equations take into account the explicit dependence on the past in the prediction of the future.As a consequence, much more complicated dynamics (when compared with ordinary differential equations) needs to be treated, even in the case of a single constant delay.The initial value is represented by a so-called history function, therefore, a finite dimmensional space which was sufficient for expression of the general solution of an ODE must necessarily be changed into the infinite dimmensional-one.
As preliminaries, let us consider the initial value problem (IVP) for the general functional differential equation In the sequel work, we will focus our interest on the case of linear autonomous and homogeneous equation x (t) = L (x t ), with a continuous linear mapping where The first analytical method, which is often proposed because of its intuitiveness when dealing with discrete delays on a finite and relatively small interval, is called the Method of Steps.The core of the method, firstly proposed by Bellman in [3], is focused on the sequential forward extension of an initial solution in the direction of increasing time, which yields to an explicit solution of the delay differential equation represented by corresponding system of ODEs.In view of the future comparison, we assume there with no loss of generality an exponential initial function φ (t) = e −at .
A simple substitution y(t) = e at x(t) can be used to obtain a simplified pure-delay equation y (t) = −be aτ y(t − τ) with pertaining constant initial function φ (t) = 1, whose solution can be easily deduced using the Method of Steps approach as follows where [p] denotes the integer part of p and α = be aτ .The solution is obtained in terms of the original variable It should be noted that a similar explicit solution is not possible if the time-dependent initial function is taken into account and generated chain of ordinary differential equations becomes more and more complicated.Of course, difficult or lengthy integral calculations can be performed by appropriate tools of Computer Algebra Systems.Apparently, time needed for the computation is inversely proportional to the delay value.On the contrary, the power of the Lambert W function in solving such equations has been recently proven and extended as much as possible (see [1,23,24] to name a few).As will be seen, the complexity of the calculation is delay-independent which makes analysis much more general for any initial function φ (t).
A formal solution similarity with an ordinary differential equation may be reached using the Lambert W function.
Restricting our interest to the nontrivial solution, the characteristic equation is obtained by substituting the ansatz in the form x = Ce λt into the equation (3) as follows: λ + a + be −λ τ = 0 Multiplying both sides by e λt , we can rewrite the equation as Then, a multiplication by the factor τe aτ leads to τ(λ + a)e (λ +a)τ = −τbe aτ Solving the equation in λ , we get a solution to the transcendental characteristic equation in terms of Lambert W function: The existence of countably infinite branches, which results from the multivalued nature of the Lambert W function has already been stated.Therefore, using the principle of superposition, the solution of a linear homogeneous initial value problem can be expressed as a linear combination of an infinite number of particular branches of the Lambert W function as below for all t ≥ 0, where C k ∈ C represents particular coefficients, which can be efficiently determined by an approximation method proposed by Asl and Ulsoy (see [1]) using the preshape interval [−τ, 0], divided into 2N + 1 divisions depending on the number of branches considered in the real calculations, on which the history function is evaluated.An attractive and novel approach in determining coefficients performed by Sun Yi [24] focuses on the comparison with the solution in s-domain with the following result: where Φ(•) is the Laplace transform of φ (•).Thus, the complete solution of the homogeneous initial problem is given by: where λ k denotes kth characteristic root produced by W k .
Example 3.1.Let us consider the following IVP Figure 3 illustrates the distribution of the corresponding characteristic roots.The complete analytical solution is expressed as follows The solution is approximated by the sum of the first 5 branches, i.e. when N = 2.A comparison with exact solution is given in the Figure 2.

QUALITATIVE ANALYSIS
Moreover, in addition to closed-form solution, it would be of great interest from a practical standpoint to investigate various techniques in order to analyze the qualitative impact of the delays and to better understand the behavior of systems.A typical phenomenon of such equations is the possible occurrence of solutions exhibiting oscillatory behavior as same as time-delay induced instabilities.
In recent years, there has been much research activity devoted to the oscillatory and asymptotic behavior of solutions of delay differential equations.The reader is referred to [10,16], and the references therein.By convention, the solution is called oscillatory if it has arbitrarily large zeros and otherwise, it is called nonoscillatory.The equation itself is said to be oscillatory if all its solutions are oscillatory.Moreover, the equation is said to be asymptotically stable if all its solutions tend to zero as t → ∞.
Let us return back to the Lambert W function.Several facts result from a simple analysis performed in the previous section.First, non-oscillatory solutions can only be produced by some of two real branches W 0 and W −1 .Since a necessary and sufficient condition in order that the zero solution of be asymptotically stable in the sense of Lyapunov is that all characteristic roots have negative real parts, the following lemma will be very useful in further considerations.
To start with, assume first the simplest pure-delay differential equation in the form: with corresponding characteristic equation λ + be −λ τ = 0 We have already known that the transcendental nature of the characteristic equation (multiple zeros on the imaginary axis) makes the analysis more complicated.However, it has some interesting properties, which may be observed using the lemmas above: -There is only a finite number of zeros of the characterstic equation with positive real part and the rightmost is obtained using the principal branch only, where the asymptotic behavior of the solution can be sufficiently observed.
-There exists a real positive constant β such that all the zeros λ k 's of the characteristic equation satisfy Re(λ k ) ≤ β .
We were looking for the exponentials in the form x(t) = e λt , where λ is assumed to be complex.Therefore, the real and imaginary part of the rightmost root λ 0 are easily determined as Re and the asymptotic behavior is readily seen from Figures 1 and 2.
The following well-known theorem (see e.g.[16]) is hence easily derived using the Lambert W function.Moreover, if 1 e < bτ < π 2 holds, every solution tends to zero as t → ∞.
Proof.It is well known that every solution of (8) oscillates if and only if its characteristic equation has no real roots (see e.g.[10]).W (z) is complex-valued for z < −1/e, which directly implies bτ > 1/e.Furthermore, z > − π 2 , i.e., bτ < π 2 indicates a negative real part of the rightmost root of the characteristic equation, which is necessary and sufficient for an asymptotically stable solution.
Remark 4.1.Asymptotic behavior of the solution of (8) depending on the combination of parameters can be stated in detail as follows: unstable and the solution is non-oscillatory.
Apparently, if τ = 0, the equation is always stable and nonoscillatory.
Consider now the case of the equation ( 3).The change of variable which has been shown in the previous section allows us to use already obtained criteria.Logically, a condition for an oscillatory and asymptotically stable solution should be given by τbe aτ ∈ (1/e, π/2).
As will be shown, it is not entirely accurate.Using the previous analysis of the Lambert W function, we derive in detail a sufficient and necessary condition in respect to parameter values.First, we have just known that the rightmost root of the characteristic equation is expressed by: Therefore, the parametrized curve determining the second boundary of asymptotically stable solution in (a, b)-space is given as follows: where τ is fixed.By taking limits Remark 4.2.The asymptotics of (3) depend on delay below and above the critical value (see Figure 5).The following holds: -If a + b > 0 and bτe aτ < 1/e, a > −1/τ: Re[λ 0 ] < 0, Im[λ 0 ] = 0 and the zero point is asymptotically stable (if a < −1/τ, x = 0 is unstable.) -If a + b > 0, bτe aτ > 1/e and τ ∈ (0, τ): every solution oscillates around zero with decaying amplitude of oscillations.
- Finally, in order to study the asymptotic behavior of the oscillating periods and the amplitudes of the oscillatory solutions, it is useful to note an interesting property of the Lambert W function, whose imaginary parts change very slightly for negative arguments.Figure 4 illustrates the imaginary parts of several branches, which asymptotically approach to the corresponding multiple of π as z → −∞.The real part of the Lambert function is unbounded, although the same property can be observed.Due to this fact, the period of the oscillation much less depends on the value of a and b (occuring only in the argument of the Lambert W function) in comparison with the high dependence on the delay value, and it becomes larger and larger with the delay increase.However, it is easily seen from ( 5) that the amplitude of oscillations given by the real part of λ also depends on a because of the separated term in (5).

ODD-ORDER CASE
Consider now a little extension of the Lambert W function approach on a higher order linear delay differential equation, on which a similar procedure can be applicable.The oscillatory behavior of all solutions is established for the following odd-order delay differential equation It is well known from [10] that a case of even n directly implies an oscillatory behavior for b ∈ R + ; that is why we restrict our interest to odd-order equations.Characteristic equation of ( 12) is given as follows The rightmost characteristic root is easily determined based on result of Lemma 2.1 as It is clear that (12)  Hence, the necessary and sufficient for oscillatory and asymptotically stable solution of (13) can be explicitly given.Moreover, it is useful to realize that the oscillation period decreases with increasing n, while the amplitude of oscillations increases.
In the following, we compare the obtained results with an existing one.Let us consider more special case In [2], Džurina and Baculíková established new criteria for oscillation of the following third-order functional differential equation Their method is focused on the comparison principle: assuming the existence of the auxiliary function ξ (t) ∈ C 1 ([t 0 , ∞) such that ξ (t) ≥ 0, ξ (t) > t, and η(t) = τ(ξ (ξ (t))); then, if the following first order delay differential equations oscillate, the equation ( 13) is oscillatory (for a detailed proof of theorem see).
In our case, eq.(i) reduces to y (t) + a (t−τ) 2 2 y(t − τ) = 0 and evidently, lim holds.Let the auxiliary function be defined by ξ (t) = t + c, c > 0. Then ξ (ξ (t)) = t +2c and τ(ξ (ξ (t))) = t +2c−τ, where τ(t) < t, hence c ∈ (0, τ 2 ).The following equation is given Finding the maximum value in b = τ 3 , we obtain the following criterion for oscillatory solutions of (13) b > 27 eτ 3 , which is evidently much stronger than necessary.Of course, this is understandable because of the generality of the studied equation.In principle, the criterion depends on the choice of the auxiliary function ξ (t), which cannot be generalized.The method approach in [2], as well as in other publications dealing with higher-order delay differential equations (see e.g.[9,16]), does not provide any criteria of asymptotic stability of oscillatory solutions, which can be obtained using the Lambert W function.
Let us illustrate the above knowledge in the following example.
Example 5.1.Consinder the equation The first aim is to derive an analytical solution for the particular equation.Then, oscillatory and asymptotic properties are examinated.Following [24], we can easily determine coefficient's expression for higher-order differential equations as same as for first-order ones.
Taking the Laplace transform of (14) as .
Then, the inverse substitution of Q(λ k ) yields to the expression of C k as: .
Complete analytical solution of equation ( 14) is then expressed as follows Remark 5.1.Following [4], the most general possible case of the solvable characteristic equation in terms of the Lambert W function, that is when a characteristic polynomial has as many roots as its degree, is as follows: (λ + a) n = −be −τλ , from which a characteristic roots are determined as follows In this paper, we have studied the oscillatory and asymptotic properties of a simple linear delay differential equation using the Lambert W function approach, which was briefly introduced.The explicit solution obtained by the Method of Steps was compared to the approximative solution given by several branches of the special function.Certain noticeable properties of the characteristic equation observed via Lambert W function were stated and used in subsequent analysis.Next, we have derived a critical delay value, which causes periodic oscillations around zero.Then, we have provided an extension on the odd-order case, where the complete analytical solution was derived and illustrated on the numerical example.The criteria on the asymptotically stable and oscillatory solution were established and compared with existing ones.The dependence of the oscillation period on the parameter values was also analyzed.

Lemma 2 . 1 .
A solution of the following generalized odd-degree exponential polynomial k ax+b = (cx + d) n , k > 0, a, c = 0 can be expressed in terms of the Lambert function Taking the nth root, we find k ax+b n = cx + d, while separation of the unknown to the right side of the equation yields 1 c k b/n = (x + d/c)e − ax n ln k .Performing straightforward mathematical calculation, one can easily obtain − a ln k cn e bc−ad nc = − a ln k n (x + d/c)e − a ln k n (x+d/c) .Hence, directly from the definition of the Lambert W function, we have W (− a ln k n

Fig. 4
Fig.4 Solution of (6) using the exact relation (solid) and the sum of first three branches of W function (dashed)

√ b 2 − a 2 Theorem 4 . 2 .
it is easily seen that the intersection point of the first and second boundary, that is (a, b) = (−1/τ, 1/τ), lies on the line representing our trivial solution.As y approaches π, a, b goes to infinity .In accordance with our previous conclusions, the curve intersects the b-axis at (0, π/2τ) -which indicates the purely oscillatory solution of the delay differential equation with no ODE component treated above.Combining the equations (9) and (10) one gets an expression of y: y = arccos(−a/b) Substitution into (1) yields z = −aτe aτ a b − 1 b arccos − a b e aτ b 2 − a 2 Since our argument is z = −τbe aτ , the necessary condition for the existence of a pair of pure imaginary roots of the characteristic equation is obtained as below: bτ = 1 b a 2 τ + arccos − a b b 2 − a 2 from which a critical time lag can explicitly be given by τ = arccos − a b Assuming a + b > 0, b > a, the effect of such a critical delay (henceforth denoted as τ) is to cause free (constant amplitude) oscillations around the zero point.The following theorem summarizes the above knowledge.Assume that τ ∈ R + and a, b ∈ R. Then the following statements are equivalent.(a) Every solution of (3) is oscillatory.(b) bτe aτ > 1/e.Moreover, if τ < arccos(− a b ) √ b 2 −a 2 holds, every solution tends to zero as t → ∞.

−λ k t dt 2λ 2 k 3 √ 0 . 5 / 2 )
− e −2λ k × e λ k t , with λ k = 3 2 W k (−3Moreover, based on the obtained criteria, we can easily compute that the solution is oscillatory and tends to zero as t → ∞.
a, k = 0, ±1, ±2, . . ., ∞ The oscillations of solutions of such characteristic equation are then guaranteed by the condition e τa b > n

Table 1
The values of λ k 's and C k 's for k = 0, ±1, ±2