# Numerical Stability Test of Neutral Delay Differential Equations

Authors: Z. H. Wang
Journal: Mathematical Problems in Engineering
 (2008)

---

## Abstract

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.


---

## Body

## 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 NDDEdirectly 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 W function [17] is applied.In
this paper, we are interested in the stability test of NDDEs,(1.1)x˙(t)+Cx˙(t−τ)=Ax(t)+Bx(t−τ),x∈ℝn,
whose
characteristic equations are assumed in the form
(1.2)Δ(λ):=(1−pe−λτ)λn+a1(e−λτ)λn−1+⋯+an(e−λτ)=0,
where p is a constant, and the coefficients aj(z) are polynomials in z. The
systems discussed in [1–9] fall into the category of (1.1). If p=0,
then the trivial solution x=0 is asymptotically stable if and only if Δ(λ) has roots with negative real parts only [10–14],
which is equivalent to α<0,
where
(1.3)α=max⁡{Reλ:Δ(λ)=0,λ∈ℂ}.
This may
not be the case of NDDEs. If |p|=1,
the condition “Δ(λ) has roots with negative real parts only” is not equivalent to α<0, because the
infinite characteristic roots with negative real parts may accumulate on the imaginary axis as shown in
Section 2.1. If |p|>1, Δ(λ) has always roots with positive real part, so x=0 is unstable for all given τ>0.
Thus, only |p|<1 is assumed true in this paper, for which “Δ(λ) has roots with negative real parts only” is equivalent to α<0.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 W 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(2.1)x˙(t)+ax(t)+bx(t−τ)=0.
The characteristic equation
corresponding to the trivial solution x=0 is
(2.2)Δ(λ):=λ+a+be−λτ=0.
On the basis of Lambert W function,
the stability condition can be presented explicitly. In fact, Lambert W
function w=W(z) is defined as the solution of a transcendental
equation
(2.3)wew=z,(z∈ℂ).
It has infinite branches, denoted by Wk(z),k=0,±1,±2,…,
respectively. W0(z) is the unique branch that is analytic at the
origin z=0, and is called the principal branch. For more
details about Lambert W function, it is referred to [17–19].Now, ifΔ(λ)=0, then (λ+a)eλτ=−b, and τ(λ+a)e(λ+a)τ=−bτeaτ. Thus, the
characteristic roots can be expressed explicitly in terms of Lambert W function
[17]
(2.4)λk=−a+Wk(−bτeaτ)τ,k=0,±1,±2,….
Moreover, it has
been proved in [19]: for
arbitrary z∈ℂ,
one has
(2.5)max⁡k=0,±1,±2,…ReWk(z)=ReW0(z),
where Rez stands for the real part of z∈ℂ. Thus, the rightmost root of (2.1) is λ0,
and the trivial solution, x=0, of (2.1) is asymptotically stable if and
only if Reλ0<0,
namely,
(2.6)ReW0(−bτeaτ)<τRea.
Such a stability condition can be checked easily, because Maple, Matlab,
and Mathematica, the three popular mathematical softwares, provide a calculator
of Lambert W 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
(2.7)Δ(λ):=λn+a1(e−λτ1,e−λτ2,…,e−λτm)λn−1+⋯+an(e−λτ1,e−λτ2,…,e−λτm)=0,
where the
coefficients aj(z1,z2,…,zm) 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λ0 and then refine it to be the rightmost root of
the following polynomial equation
(2.8)λn+a1(e−λ0τ1,e−λ0τ2,…,e−λ0τm)λn−1+⋯+an(e−λ0τ1,e−λ0τ2,…,e−λ0τm)=0,
where the
coefficients aj are assigned to fixed values. Note that the
notation λ0 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 constantsa>0 and b, define
(2.9)F(λ)=aλ+b−W0((aλ+b−Δ(λ))eaλ+b),
where W0(z) is the principal branch of Lambert W function,
and the constants a,b are not
large to avoid numerical problems due to large factor. Then the Newton-Raphson
method is employed to find the rightmost root of Δ(λ)=0;
(2.10)λi+1=λi−F(λi)F′(λi),(i=0,1,2,…);
which has
quadratic convergence for unrepeated roots. Alternatively, Halley's method
(2.11)λi+1=λi−(1−12F(λi)F"(λi)(F′(λi))2)−1F(λi)F′(λi),(i=0,1,2,…) can be used. This
algorithm has order 3 of convergence for unrepeated roots. The iteration is stopped
at step N if
(2.12)|λN−λN−1|<ε, for a given tolerance ε.


#### 2.2.3. Verification of the Computational Result

Due to (2.5), it is expected that λN resulted from (2.10) or (2.11) is the rightmost root of the
delay differential equation, namely, for any root λ of Δ(λ)=0, one has
(2.13)Reλ≤ReλN.
Equation (2.13) is guaranteed if the Nyquist
plot of Δ(iω+ReλN)/(1+iω)n passes through the origin of the complex plane
and the Nyquist plot of Δ(iω+ReλN+μ)/(1+iω)n does not encounter the origin for very small μ>0. 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 |p|<1.


### 2.3. Accumulation of the Characteristic Roots

The different
braches of LambertW function can be used to find different roots of Δ(λ)=0. For
simplicity, let us calculate the roots of a first-order autonomous NDDE x˙(t)=−x˙(t−τ)−ax(t) with
(2.14)Δ(λ)=λ+λe−λτ+a.
Let Wk(z) be the kth branch of Lambert W function defined in Section 2.1, and
define
(2.15)Fk(λ)=λ+a−Wk(−τλeaτ)/τ,(k=0,±1,±2,…),
then all the
characteristic roots, computed by using the Newton-Raphson method or Halley's
method for Fk(λ), have
negative real parts, and they accumulate on the imaginary axis as shown in Table 1. As a result, the solution x=0 is not stable, though all the roots of Δ(λ) have negative real parts.Table 1
Numerical calculation of the roots ofΔ(λ) in (2.14)
with τ=0.3,a=0.5.
BranchCharacteristic rootBranchCharacteristic root−200−0.2572×10−9−0.4178×104i1−0.3631×10−2+0.1063×102i−100−0.9644×10−9−0.2084×104i2−0.4200×10−3+0.3147×102i−20−0.2498×10−5−0.4084×103i3−0.1517×10−3+0.5239×102i−10−0.1052×10−4−0.1990×103i4−0.7747×10−4+0.7333×102i−5−0.4688×10−4−0.9427×102i5−0.4688×10−4+0.9427×102i−4−0.7747×10−4−0.7333×102i10−0.1052×10−4+0.1990×103i−3−0.1517×10−3−0.5239×102i20−0.2498×10−5+0.4084×103i−2−0.4200×10−3−0.3147×102i100−0.9644×10−9+0.2084×104i−1−0.3631×10−2−0.1063×102i200−0.2572×10−9+0.4178×104i0−0.2410500−0.1386×10−10+0.1046×105iThe condition|p|<1, required in
the proposed algorithm for finding the rightmost root, is not satisfied in this
example. The branch W0(z) yields the leftmost root, rather than the
rightmost root.


## 2.1. An Explicit Stability Criterion for a First-Order RDDE

Let us consider a first-order
retarded delay differential equation described by(2.1)x˙(t)+ax(t)+bx(t−τ)=0.
The characteristic equation
corresponding to the trivial solution x=0 is
(2.2)Δ(λ):=λ+a+be−λτ=0.
On the basis of Lambert W function,
the stability condition can be presented explicitly. In fact, Lambert W
function w=W(z) is defined as the solution of a transcendental
equation
(2.3)wew=z,(z∈ℂ).
It has infinite branches, denoted by Wk(z),k=0,±1,±2,…,
respectively. W0(z) is the unique branch that is analytic at the
origin z=0, and is called the principal branch. For more
details about Lambert W function, it is referred to [17–19].Now, ifΔ(λ)=0, then (λ+a)eλτ=−b, and τ(λ+a)e(λ+a)τ=−bτeaτ. Thus, the
characteristic roots can be expressed explicitly in terms of Lambert W function
[17]
(2.4)λk=−a+Wk(−bτeaτ)τ,k=0,±1,±2,….
Moreover, it has
been proved in [19]: for
arbitrary z∈ℂ,
one has
(2.5)max⁡k=0,±1,±2,…ReWk(z)=ReW0(z),
where Rez stands for the real part of z∈ℂ. Thus, the rightmost root of (2.1) is λ0,
and the trivial solution, x=0, of (2.1) is asymptotically stable if and
only if Reλ0<0,
namely,
(2.6)ReW0(−bτeaτ)<τRea.
Such a stability condition can be checked easily, because Maple, Matlab,
and Mathematica, the three popular mathematical softwares, provide a calculator
of Lambert W 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
(2.7)Δ(λ):=λn+a1(e−λτ1,e−λτ2,…,e−λτm)λn−1+⋯+an(e−λτ1,e−λτ2,…,e−λτm)=0,
where the
coefficients aj(z1,z2,…,zm) 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λ0 and then refine it to be the rightmost root of
the following polynomial equation
(2.8)λn+a1(e−λ0τ1,e−λ0τ2,…,e−λ0τm)λn−1+⋯+an(e−λ0τ1,e−λ0τ2,…,e−λ0τm)=0,
where the
coefficients aj are assigned to fixed values. Note that the
notation λ0 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 constantsa>0 and b, define
(2.9)F(λ)=aλ+b−W0((aλ+b−Δ(λ))eaλ+b),
where W0(z) is the principal branch of Lambert W function,
and the constants a,b are not
large to avoid numerical problems due to large factor. Then the Newton-Raphson
method is employed to find the rightmost root of Δ(λ)=0;
(2.10)λi+1=λi−F(λi)F′(λi),(i=0,1,2,…);
which has
quadratic convergence for unrepeated roots. Alternatively, Halley's method
(2.11)λi+1=λi−(1−12F(λi)F"(λi)(F′(λi))2)−1F(λi)F′(λi),(i=0,1,2,…) can be used. This
algorithm has order 3 of convergence for unrepeated roots. The iteration is stopped
at step N if
(2.12)|λN−λN−1|<ε, for a given tolerance ε.


### 2.2.3. Verification of the Computational Result

Due to (2.5), it is expected that λN resulted from (2.10) or (2.11) is the rightmost root of the
delay differential equation, namely, for any root λ of Δ(λ)=0, one has
(2.13)Reλ≤ReλN.
Equation (2.13) is guaranteed if the Nyquist
plot of Δ(iω+ReλN)/(1+iω)n passes through the origin of the complex plane
and the Nyquist plot of Δ(iω+ReλN+μ)/(1+iω)n does not encounter the origin for very small μ>0. 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 |p|<1.


## 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λ0 and then refine it to be the rightmost root of
the following polynomial equation
(2.8)λn+a1(e−λ0τ1,e−λ0τ2,…,e−λ0τm)λn−1+⋯+an(e−λ0τ1,e−λ0τ2,…,e−λ0τm)=0,
where the
coefficients aj are assigned to fixed values. Note that the
notation λ0 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 constantsa>0 and b, define
(2.9)F(λ)=aλ+b−W0((aλ+b−Δ(λ))eaλ+b),
where W0(z) is the principal branch of Lambert W function,
and the constants a,b are not
large to avoid numerical problems due to large factor. Then the Newton-Raphson
method is employed to find the rightmost root of Δ(λ)=0;
(2.10)λi+1=λi−F(λi)F′(λi),(i=0,1,2,…);
which has
quadratic convergence for unrepeated roots. Alternatively, Halley's method
(2.11)λi+1=λi−(1−12F(λi)F"(λi)(F′(λi))2)−1F(λi)F′(λi),(i=0,1,2,…) can be used. This
algorithm has order 3 of convergence for unrepeated roots. The iteration is stopped
at step N if
(2.12)|λN−λN−1|<ε, for a given tolerance ε.


## 2.2.3. Verification of the Computational Result

Due to (2.5), it is expected that λN resulted from (2.10) or (2.11) is the rightmost root of the
delay differential equation, namely, for any root λ of Δ(λ)=0, one has
(2.13)Reλ≤ReλN.
Equation (2.13) is guaranteed if the Nyquist
plot of Δ(iω+ReλN)/(1+iω)n passes through the origin of the complex plane
and the Nyquist plot of Δ(iω+ReλN+μ)/(1+iω)n does not encounter the origin for very small μ>0. 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 |p|<1.


## 2.3. Accumulation of the Characteristic Roots

The different
braches of LambertW function can be used to find different roots of Δ(λ)=0. For
simplicity, let us calculate the roots of a first-order autonomous NDDE x˙(t)=−x˙(t−τ)−ax(t) with
(2.14)Δ(λ)=λ+λe−λτ+a.
Let Wk(z) be the kth branch of Lambert W function defined in Section 2.1, and
define
(2.15)Fk(λ)=λ+a−Wk(−τλeaτ)/τ,(k=0,±1,±2,…),
then all the
characteristic roots, computed by using the Newton-Raphson method or Halley's
method for Fk(λ), have
negative real parts, and they accumulate on the imaginary axis as shown in Table 1. As a result, the solution x=0 is not stable, though all the roots of Δ(λ) have negative real parts.Table 1
Numerical calculation of the roots ofΔ(λ) in (2.14)
with τ=0.3,a=0.5.
BranchCharacteristic rootBranchCharacteristic root−200−0.2572×10−9−0.4178×104i1−0.3631×10−2+0.1063×102i−100−0.9644×10−9−0.2084×104i2−0.4200×10−3+0.3147×102i−20−0.2498×10−5−0.4084×103i3−0.1517×10−3+0.5239×102i−10−0.1052×10−4−0.1990×103i4−0.7747×10−4+0.7333×102i−5−0.4688×10−4−0.9427×102i5−0.4688×10−4+0.9427×102i−4−0.7747×10−4−0.7333×102i10−0.1052×10−4+0.1990×103i−3−0.1517×10−3−0.5239×102i20−0.2498×10−5+0.4084×103i−2−0.4200×10−3−0.3147×102i100−0.9644×10−9+0.2084×104i−1−0.3631×10−2−0.1063×102i200−0.2572×10−9+0.4178×104i0−0.2410500−0.1386×10−10+0.1046×105iThe condition|p|<1, required in
the proposed algorithm for finding the rightmost root, is not satisfied in this
example. The branch W0(z) yields the leftmost root, rather than the
rightmost root.


## 3. Two Illustrative Examples

In this
section, the iterative method proposed in Section2 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
(3.1)x¨(t)+2ξx˙(t)+x(t)+px¨(t−τ)=0,
where 0<p<1. The
characteristic equation is
(3.2)Δ(λ):=λ2+2ξλ+1+pλ2e−λτ=0.
Each root
must be a root of a certain branch of the following equation:
(3.3)Fk(λ):=2ξλ+1−Wk(−λ2(1+pe−λτ)e2ξλ+1)=0,(k=0,±1,±2,…).
Then the
rightmost root can be found out from F0(λ)=0 via the Newton-Raphson method or Halley's
method. For example, let p=0.2,ξ=0.05, and
calculate the rightmost root for four special cases: τ=1,5,7.5,11.To this end, one chose freely an initial guess,
sayλ0=1.0+3.0i. Because the
simplified polynomial equation λ2+2ξλ+1+pλ2e−λ0=0,
corresponding to (2.8), has two complex roots
(3.4)−0.5971×10−1+0.1036×10i,−0.4813×10−1+0.1038×10i,
the initial
guess can be refined as λ0=−0.4813×10−1+0.1038×10i for τ=1. 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 x=0 is asymptotically stable, and consequently,
the rightmost root must have negative real part.Graphical test for the
rightmost root of (3.1): (a) the Nyquist plot of Δ(iω)/(1+iω)2; (b) the Nyquist plot of Δ(iω−0.1155)/(1+iω)2.
(a)(b)With thisλ0, the 4th iteration
of the Newton-Raphson method gives λ4=−0.1155−0.9220i. As shown in Figure 1(b), the Nyquist plot of Δ(iω−0.1155)/(1+iω)2 passes through the origin, it follows that λ4=−0.1155−0.9220i is the rightmost root for τ=1.Similarly, starting fromλ0=−0.4813×10−1+0.1038×10i, the 4th
iteration, 6th iteration, and 4th iteration of the Newton-Raphson method give the
rightmost roots 0.3125×10−1+0.9810i, −0.6278×10−1−0.8692i, and 0.2546×10−1−0.9986i for τ=5,7.5,11, respectively. The numerical results are in
agreement with that obtained in [7].Moreover, as shown in Figure2, 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 λ0=−0.4813×10−1+0.1038×10i, we see that
the trivial solution x=0 of (3.1) exhibits two stability
switches in τ∈[0,7], occurs at τ10=3.2742, and τ20=6.1742, respectively, and the corresponding rightmost
roots can be found to be ±1.1032i, ±0.9252i, respectively. This is the same result as that
obtained on the basis of stability switches [7]. In fact, Δ(iω)=0 gives (1−p2)ω4+(4ξ2−2)ω2+1=0. When p=0.2,ξ=0.05, one has two
positive roots ω1=1.1032 and ω2=0.9252, and the
corresponding minimal critical delays, determined from Δ(iω)=0, are found
easily to be τ10=3.2742 and τ20=6.1742, respectively.Figure 2
The solutionx=0 of (3.1) is asymptotically stable if τ∈[0,3.2742)∪(6.1742,7], or it is
unstable if τ∈(3.2742,6.1742).From
Figure2, we see also that at τ=2.12,  the rightmost root −0.2201−0.1061×10i 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]
(3.5)ddt(x(t)−34x(t−1)+12x(t−2))=14x(t)+34x(t−1)
to show that
the method work also for NDDEs with multiple delays. The characteristic
equation is
(3.6)Δ(λ):=λ(1−34e−λ+12e−2λ)−14−34e−λ=0.
As shown in Figure 3(a), the Nyquist plot of Δ(iω)/(1+iω) encounters the origin of the complex plane, so
the trivial solution x=0 of (3.5) is unstable, and the rightmost root
must have positive real part. To find out the rightmost root, let
(3.7)F(λ)=λ−14−W0((λ−14−Δ(λ))eλ−1/4),
where W0(z) is the principal branch of Lambert W function.
Then, choose freely an initial guess, say λ0=1.0+2.0i, and refine
it as the root of λ(1−3e−λ0/4+e−2λ0/2)−1/4−3e−λ0/4=0, namely,
replace it with λ0=−0.9513×10−1−0.2924i. Then the third
iteration of the Newton-Raphson method gives λ3=0.7581. Moreover, the Nyquist plot in Figure 3(b) shows that λ3=0.7581 is the rightmost root, because the
Nyquist plot of Δ(iω+0.7581+0.001)/(1+iω) does not encounter the origin. This result is
the same as the one obtained by using DDE-BIFTOOL [22, 23].Graphical test for the rightmost root
of (3.5): (a) the Nyquist plot of Δ(iω)/(1+iω); (b) the
Nyquist plot of Δ(iω+0.7581+0.001)/(1+iω).
(a)(b)Moreover, when a negative feedback
control−ux(t) is performed on to (3.5),
(3.8)ddt(x(t)−34x(t−1)+12x(t−2))−14x(t)−34x(t−1)=−ux(t),
the unstable equilibrium is stabilized if one chooses, for example, 0.5<u<1.5,
as shown in Figure 4.Figure 4
The plot of the real part of the rightmost
root with respect to the feedback gainu for (3.8).


## 3.1. A Neutral Delayed Oscillator

Let us
firstly consider the stability of a second-order NDDE [7] arising from
structure dynamics
(3.1)x¨(t)+2ξx˙(t)+x(t)+px¨(t−τ)=0,
where 0<p<1. The
characteristic equation is
(3.2)Δ(λ):=λ2+2ξλ+1+pλ2e−λτ=0.
Each root
must be a root of a certain branch of the following equation:
(3.3)Fk(λ):=2ξλ+1−Wk(−λ2(1+pe−λτ)e2ξλ+1)=0,(k=0,±1,±2,…).
Then the
rightmost root can be found out from F0(λ)=0 via the Newton-Raphson method or Halley's
method. For example, let p=0.2,ξ=0.05, and
calculate the rightmost root for four special cases: τ=1,5,7.5,11.To this end, one chose freely an initial guess,
sayλ0=1.0+3.0i. Because the
simplified polynomial equation λ2+2ξλ+1+pλ2e−λ0=0,
corresponding to (2.8), has two complex roots
(3.4)−0.5971×10−1+0.1036×10i,−0.4813×10−1+0.1038×10i,
the initial
guess can be refined as λ0=−0.4813×10−1+0.1038×10i for τ=1. 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 x=0 is asymptotically stable, and consequently,
the rightmost root must have negative real part.Graphical test for the
rightmost root of (3.1): (a) the Nyquist plot of Δ(iω)/(1+iω)2; (b) the Nyquist plot of Δ(iω−0.1155)/(1+iω)2.
(a)(b)With thisλ0, the 4th iteration
of the Newton-Raphson method gives λ4=−0.1155−0.9220i. As shown in Figure 1(b), the Nyquist plot of Δ(iω−0.1155)/(1+iω)2 passes through the origin, it follows that λ4=−0.1155−0.9220i is the rightmost root for τ=1.Similarly, starting fromλ0=−0.4813×10−1+0.1038×10i, the 4th
iteration, 6th iteration, and 4th iteration of the Newton-Raphson method give the
rightmost roots 0.3125×10−1+0.9810i, −0.6278×10−1−0.8692i, and 0.2546×10−1−0.9986i for τ=5,7.5,11, respectively. The numerical results are in
agreement with that obtained in [7].Moreover, as shown in Figure2, 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 λ0=−0.4813×10−1+0.1038×10i, we see that
the trivial solution x=0 of (3.1) exhibits two stability
switches in τ∈[0,7], occurs at τ10=3.2742, and τ20=6.1742, respectively, and the corresponding rightmost
roots can be found to be ±1.1032i, ±0.9252i, respectively. This is the same result as that
obtained on the basis of stability switches [7]. In fact, Δ(iω)=0 gives (1−p2)ω4+(4ξ2−2)ω2+1=0. When p=0.2,ξ=0.05, one has two
positive roots ω1=1.1032 and ω2=0.9252, and the
corresponding minimal critical delays, determined from Δ(iω)=0, are found
easily to be τ10=3.2742 and τ20=6.1742, respectively.Figure 2
The solutionx=0 of (3.1) is asymptotically stable if τ∈[0,3.2742)∪(6.1742,7], or it is
unstable if τ∈(3.2742,6.1742).From
Figure2, we see also that at τ=2.12,  the rightmost root −0.2201−0.1061×10i 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]
(3.5)ddt(x(t)−34x(t−1)+12x(t−2))=14x(t)+34x(t−1)
to show that
the method work also for NDDEs with multiple delays. The characteristic
equation is
(3.6)Δ(λ):=λ(1−34e−λ+12e−2λ)−14−34e−λ=0.
As shown in Figure 3(a), the Nyquist plot of Δ(iω)/(1+iω) encounters the origin of the complex plane, so
the trivial solution x=0 of (3.5) is unstable, and the rightmost root
must have positive real part. To find out the rightmost root, let
(3.7)F(λ)=λ−14−W0((λ−14−Δ(λ))eλ−1/4),
where W0(z) is the principal branch of Lambert W function.
Then, choose freely an initial guess, say λ0=1.0+2.0i, and refine
it as the root of λ(1−3e−λ0/4+e−2λ0/2)−1/4−3e−λ0/4=0, namely,
replace it with λ0=−0.9513×10−1−0.2924i. Then the third
iteration of the Newton-Raphson method gives λ3=0.7581. Moreover, the Nyquist plot in Figure 3(b) shows that λ3=0.7581 is the rightmost root, because the
Nyquist plot of Δ(iω+0.7581+0.001)/(1+iω) does not encounter the origin. This result is
the same as the one obtained by using DDE-BIFTOOL [22, 23].Graphical test for the rightmost root
of (3.5): (a) the Nyquist plot of Δ(iω)/(1+iω); (b) the
Nyquist plot of Δ(iω+0.7581+0.001)/(1+iω).
(a)(b)Moreover, when a negative feedback
control−ux(t) is performed on to (3.5),
(3.8)ddt(x(t)−34x(t−1)+12x(t−2))−14x(t)−34x(t−1)=−ux(t),
the unstable equilibrium is stabilized if one chooses, for example, 0.5<u<1.5,
as shown in Figure 4.Figure 4
The plot of the real part of the rightmost
root with respect to the feedback gainu for (3.8).


## 4. Conclusions

In this
paper, the iterative method based on LambertW 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.



---

*Source: 698043-2008-04-14.xml*