Institute of Vibration Engineering Research, Nanjing University of Aeronautics and Astronautics, 210016 Nanjing, China
The stability of a delay differential equation can be investigated
on the basis of the root location of the characteristic function. Though a number of stability
criteria are available, they usually do not provide any information about the characteristic root
with maximal real part, which is useful in justifying the stability and in understanding the system
performances. Because the characteristic function is a transcendental function that has an infinite
number of roots with no closed form, the roots can be found out numerically only. While some
iterative methods work effectively in finding a root of a nonlinear equation for a properly chosen
initial guess, they do not work in finding the rightmost root directly from the characteristic function.
On the basis of Lambert W function, this paper presents an effective iterative algorithm for the calculation
of the rightmost roots of neutral delay differential equations so that the stability of the delay equations
can be determined directly, illustrated with
two examples.
1. Introduction
Many engineering systems can be
modeled as neutral delay differential equations (NDDEs) that
involve a time delay in the derivative of the highest order [1–9], which are different from retarded delay differential equations
(RDDEs) that do not involve a time delay in the derivative of the highest order
[10–14]. For example, a system,
which consists of a mass mounted on a linear spring to which a pendulum is
attached via a hinged massless rod, is used to predict the dynamic response of
structures to external forces using a set of actuators, and it is modeled as an
NDDE if the delay in actuators is taken into consideration [7]. While the RDDEs have been studied intensively
in the literature (see, e.g., [10–14]), the NDDEs have been investigated relatively few. Analysis shows that,
both RDDEs and NDDEs may exhibit very complicated nonlinear dynamics. For
example, a first-order autonomous DDE can exhibit chaotic motion, and a first-order autonomous
NDDE with a single delay
can even admit homoclinic snaking [5]. Complex behaviors of dynamical systems come out
after certain stationary solutions lose their stability, so stability analysis
plays a fundamental role in system dynamics. Usually, the stability analysis for
equilibriums of DDEs can be investigated on the basis of the method of
Lyapunov's function (al) including the LMI (linear matrix inequality) method,
or by means of the root location of the characteristic functions for
equilibriums [10–14]. In particular, the stability can
be studied on the basis of stability switches [7, 8, 10, 14], if the delay effect on the stability is addressed. In this
case, the delay interval is divided into a number of subintervals by the
critical values of delay for which the system changes its equilibriums from stable
status to unstable status, or from unstable status to stable status, as the
delay passes through the critical delays, and the system has the same stability
in each subinterval. If a DDE admits a number of stability switches, then the
system can be stabilized or destabilized by adjusting the delay value only. It
is worthy to note that, however, the stability may be very poor in some delay
intervals for which the system is asymptotically stable [14]. Thus, in practical applications, it
is required not only to know whether the time-delay system is asymptotically
stable, but also to know the stability margin. Therefore, a computational
algorithm for finding the
characteristic root with maximal real part (rightmost root for short) of an NDDE is preferable.
The
characteristic quasipolynomial
of a DDE has an
infinite number of roots that do not have closed form, and the roots can be
found numerically only. Though the famous Newton-Raphson method works effectively in finding a root of a nonlinear equation for a
properly chosen initial guess, it does not work in finding the rightmost root of an NDDE directly from the characteristic quasipolynomial. In [15], an iterative scheme was proposed for the
calculation of the
rightmost root of an RDDE, where the rightmost root was assigned to be the rightmost root of a
simplified polynomial resulted from the quasipolynomial in each step of the iteration. The problem is that the iterative sequence in [15] is frequently not convergent.
Recently, the author shows in [16] that the Newton-Raphson method or the Halley method
computes effectively the rightmost roots of RDDEs if Lambert function [17] is applied.
In
this paper, we are interested in the stability test of NDDEs,
whose
characteristic equations are assumed in the form
where is a constant, and the coefficients are polynomials in . The
systems discussed in [1–9] fall into the category of (1.1). If ,
then the trivial solution is asymptotically stable if and only if has roots with negative real parts only [10–14],
which is equivalent to ,
where
This may
not be the case of NDDEs. If ,
the condition “ has roots with negative real parts only” is not equivalent to , because the
infinite characteristic roots with negative real parts may accumulate on the imaginary axis as shown in
Section 2.1. If , has always roots with positive real part, so is unstable for all given .
Thus, only is assumed true in this paper, for which “ has roots with negative real parts only” is equivalent to .
The aim of this paper is to generalize the iterative method developed in [16] for
calculating the rightmost root of RDDEs to the stability test of NDDEs. The method will be briefly
introduced in Section 2, and two examples will be given in Section 3 to
demonstrate the efficiency of the proposed method. A few concluding remarks will be given in Section 4.
2. An Algorithm for Calculating the Rightmost Characteristic Root of Time-Delay Systems
A real corresponds to a real characteristic root or a
pair of complex conjugate characteristic roots. For simplicity, such
characteristic root(s) is called the rightmost root in this paper. In this section, the Newton-Raphson
method will be combined with a special function-Lambert function to find out
the rightmost root of .
2.1. An Explicit Stability Criterion for a First-Order RDDE
Let us consider a first-order
retarded delay differential equation described by
The characteristic equation
corresponding to the trivial solution is
On the basis of Lambert function,
the stability condition can be presented explicitly. In fact, Lambert
function is defined as the solution of a transcendental
equation
It has infinite branches, denoted by
respectively. is the unique branch that is analytic at the
origin and is called the principal branch. For more
details about Lambert function, it is referred to [17–19].
Now, if , then , and . Thus, the
characteristic roots can be expressed explicitly in terms of Lambert function
[17]
Moreover, it has
been proved in [19]: for
arbitrary ,
one has
where stands for the real part of . Thus, the rightmost root of (2.1) is ,
and the trivial solution, , of (2.1) is asymptotically stable if and
only if ,
namely,
Such a stability condition can be checked easily, because Maple, Matlab,
and Mathematica, the three popular mathematical softwares, provide a calculator
of Lambert function.
2.2. A Numerical Scheme for Neutral Delay Differential Equations
It is not
possible to gain an explicit form of the rightmost root, as done above, for
other delay differential equations. Thus, an iterative algorithm was proposed
in [16] for calculating the rightmost root of RDDEs whose characteristic
equation reads
where the
coefficients are polynomials. The main points of this
iterative method are summarized as follows.
2.2.1. Choice of the Initial Guess
A properly chosen initial value is
important in the applications of iterative methods. For our problem, one can
firstly chose freely a complex number and then refine it to be the rightmost root of
the following polynomial equation
where the
coefficients are assigned to fixed values. Note that the
notation here denotes the initial guess, rather than
the rightmost root given in Section 2.1.
2.2.2. Construction of the Algorithm
For certain fixed constants and , define
where is the principal branch of Lambert function,
and the constants are not
large to avoid numerical problems due to large factor. Then the Newton-Raphson
method is employed to find the rightmost root of ;
which has
quadratic convergence for unrepeated roots. Alternatively, Halley's method
can be used. This
algorithm has order 3 of convergence for unrepeated roots. The iteration is stopped
at step if
for a given tolerance .
2.2.3. Verification of the Computational Result
Due to (2.5), it is expected that resulted from (2.10) or (2.11) is the rightmost root of the
delay differential equation, namely, for any root of , one has
Equation (2.13) is guaranteed if the Nyquist
plot of passes through the origin of the complex plane
and the Nyquist plot of does not encounter the origin for very small . The method of Nyquist plot was originally proposed in [20] for RDDEs, and extended to
NDDEs in [21].
Such a scheme
works also for the quasipolynomial defined in (1.1) for NDDEs with .
2.3. Accumulation of the Characteristic Roots
The different
braches of Lambert function can be used to find different roots of . For
simplicity, let us calculate the roots of a first-order autonomous NDDE with
Let be the kth branch of Lambert function defined in Section 2.1, and
define
then all the
characteristic roots, computed by using the Newton-Raphson method or Halley's
method for , have
negative real parts, and they accumulate on the imaginary axis as shown in Table 1. As a result, the solution is not stable, though all the roots of have negative real parts.
Table 1: Numerical calculation of the roots of
in (
2.14)
with
.
The condition , required in
the proposed algorithm for finding the rightmost root, is not satisfied in this
example. The branch yields the leftmost root, rather than the
rightmost root.
3. Two Illustrative Examples
In this
section, the iterative method proposed in Section 2 will be applied to
calculate the rightmost roots of two NDDEs discussed in the literature.
3.1. A Neutral Delayed Oscillator
Let us
firstly consider the stability of a second-order NDDE [7] arising from
structure dynamics
where . The
characteristic equation is
Each root
must be a root of a certain branch of the following equation:
Then the
rightmost root can be found out from via the Newton-Raphson method or Halley's
method. For example, let , and
calculate the rightmost root for four special cases: .
To this end, one chose freely an initial guess,
say . Because the
simplified polynomial equation ,
corresponding to (2.8), has two complex roots
the initial
guess can be refined as for . The choice
of an initial guess with negative real part is also understandable from the
Nyquist plot. Due to Figure 1(a), where the Nyquist plot does not
encounter the origin of the complex plane, so the trivial solution is asymptotically stable, and consequently,
the rightmost root must have negative real part.
Figure 1: Graphical test for the
rightmost root of (
3.1): (a) the Nyquist plot of
; (b) the Nyquist plot of
.
With this , the 4th iteration
of the Newton-Raphson method gives . As shown in Figure 1(b), the Nyquist plot of passes through the origin, it follows that is the rightmost root for .
Similarly, starting from , the 4th
iteration, 6th iteration, and 4th iteration of the Newton-Raphson method give the
rightmost roots , , and for respectively. The numerical results are in
agreement with that obtained in [7].
Moreover, as shown in Figure 2, the curve of the
real parts of the rightmost roots with respect to the delay can be produced numerically
by means of the proposed algorithm, if the delay effect on the asymptotical stability
is addressed. From Figure 2, where the initial guess is taken as , we see that
the trivial solution of (3.1) exhibits two stability
switches in , occurs at and respectively, and the corresponding rightmost
roots can be found to be , respectively. This is the same result as that
obtained on the basis of stability switches [7]. In fact, gives . When , one has two
positive roots and , and the
corresponding minimal critical delays, determined from , are found
easily to be and respectively.
Figure 2: The solution
of (
3.1) is asymptotically stable if
, or it is
unstable if
.
From
Figure 2, we see also that at , the rightmost root has the smallest negative real part. This fact
indicates that a proper chosen delay value can improve the stability of an
NDDE.
3.2. An NDDE with Two Delays
Now, let us
consider an NDDE with two delays [22]
to show that
the method work also for NDDEs with multiple delays. The characteristic
equation is
As shown in Figure 3(a), the Nyquist plot of encounters the origin of the complex plane, so
the trivial solution of (3.5) is unstable, and the rightmost root
must have positive real part. To find out the rightmost root, let
where is the principal branch of Lambert function.
Then, choose freely an initial guess, say , and refine
it as the root of , namely,
replace it with . Then the third
iteration of the Newton-Raphson method gives Moreover, the Nyquist plot in Figure 3(b) shows that is the rightmost root, because the
Nyquist plot of does not encounter the origin. This result is
the same as the one obtained by using DDE-BIFTOOL [22, 23].
Figure 3: Graphical test for the rightmost root
of (
3.5):
(a) the Nyquist plot of
; (b) the
Nyquist plot of
.
Moreover, when a negative feedback
control is performed on to (3.5),
the unstable equilibrium is stabilized if one chooses, for example, ,
as shown in Figure 4.
Figure 4: The plot of the real part of the rightmost
root with respect to the feedback gain
for (
3.8).
4. Conclusions
In this
paper, the iterative method based on Lambert function for calculating the
rightmost roots of RDDEs is extended to the stability test of a kind of NDDEs
for which the asymptotical stability is guaranteed if all the characteristic
roots have negative real parts. Two illustrative examples show that the method works
effectively. The numerical scheme enables one not only to know whether the
time-delay system is asymptotically stable, but also to know the stability
margin. A rigorous mathematical treatment of the iterative method such as the
convergence of the iterative sequence, however, is not available in this paper
and is left for future consideration.
Though the investigation is made mainly for NDDEs
with fixed parameters, the proposed scheme does work for some NDDEs with a parameter
falling in a given interval. As shown in the first illustrative example from
structure dynamics, for example, the iterative method can produce a plot of the
real part of the rightmost root with respect to the delay, from which one can
easily determine for what value of delay the system is asymptotically stable,
and for what value of delay, the system is unstable. It reveals also that a
proper chosen delay value can improve the stability of an NDDE. In the second
illustrative example, an interval of the feedback gain is determined for
stabilizing the unstable equilibrium of the NDDE by using the iterative method.
Acknowledgment
This work
was supported by FANEDD of China under Grant no. 200430, and by NSF of China
under Grant no. 10532050.