`Mathematical Problems in EngineeringVolume 2008 (2008), Article ID 698043, 10 pageshttp://dx.doi.org/10.1155/2008/698043`
Research Article

## Numerical Stability Test of Neutral Delay Differential Equations

Institute of Vibration Engineering Research, Nanjing University of Aeronautics and Astronautics, 210016 Nanjing, China

Received 24 October 2007; Revised 1 March 2008; Accepted 16 March 2008

Copyright © 2008 Z. H. Wang. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### 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.

#### 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 [19], which are different from retarded delay differential equations (RDDEs) that do not involve a time delay in the derivative of the highest order [1014]. 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., [1014]), 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 [1014]. 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 [19] 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 [1014], 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 [1719].

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.

#### References

1. G. Stépán and Z. Szabó, “Impact induced internal fatigue cracks,” in Proceedings of the ASME Design Engineering Technical Conferences (DETC '99), Las Vegas, Nev, USA, September 1999.
2. A. Bellen, N. Guglielmi, and A. E. Ruehli, “Methods for linear systems of circuit delay differential equations of neutral type,” IEEE Transactions on Circuits and Systems, vol. 46, no. 1, pp. 212–216, 1999.
3. A. G. Balanov, N. B. Janson, P. V. E. McClintock, R. W. Tucker, and C. H. T. Wang, “Bifurcation analysis of a neutral delay differential equation modelling the torsional motion of a driven drill-string,” Chaos, Solitons and Fractals, vol. 15, no. 2, pp. 381–394, 2003.
4. Z. N. Masoud, M. F. Daqaq, and N. A. Nayfeh, “Pendulation reduction on small ship-mounted telescopic cranes,” Journal of Vibration and Control, vol. 10, no. 8, pp. 1167–1179, 2004.
5. D. A. W. Barton, Dynamics and bifurcations of non-smooth delay equations, Ph.D. dissertation, University of Bristol, Bristol, UK, 2006.
6. Z. N. Masoud and A. H. Nayfeh, “Sway reduction on container cranes using delayed feedback controller,” Nonlinear Dynamics, vol. 34, no. 3-4, pp. 347–358, 2003.
7. Y. N. Kyrychko, K. B. Blyuss, A. Gonzalez-Buelga, S. J. Hogan, and D. J. Wagg, “Real-time dynamic substructuring in a coupled oscillator-pendulum system,” Proceedings of the Royal Society of London A, vol. 462, no. 2068, pp. 1271–1294, 2006.
8. Y. N. Kyrychko, S. J. Hogan, A. Gonzalez-Buelga, and D. J. Wagg, “Modelling real-time dynamic substructuring using partial delay differential equations,” Proceedings of the Royal Society of London A, vol. 463, no. 2082, pp. 1509–1523, 2007.
9. Z. N. Masoud, A. H. Nayfeh, and D. T. Mook, “Cargo pendulation reduction of ship-mounted cranes,” Nonlinear Dynamics, vol. 35, no. 3, pp. 299–311, 2004.
10. Y. Kuang, Delay Differential Equation with Applications in Population Dynamics, vol. 191 of Mathematics in Science and Engineering, Academic Press, Boston, Mass, USA, 1993.
11. Y. X. Qin, Y. Q. Liu, L. Wang, and Z. X. Zhen, Stability of Motion of Dynamical Systems with Time Lag, Science Press, Beijing, China, 2nd edition, 1989.
12. G. Stépán, Retarded Dynamical Systems: Stability and Characteristic Functions, vol. 210 of Pitman Research Notes in Mathematics Series, Longman Scientific & Technical, Harlow, UK, 1989.
13. S.-I. Niculescu, Delay Effects on Stability. A Robust Control Approach, vol. 269 of Lecture Notes in Control and Information Sciences, Springer, London, UK, 2001.
14. H. Y. Hu and Z. H. Wang, Dynamics of Controlled Mechanical Systems with Delayed Feedback, Springer, Berlin, Germany, 2002.
15. J. M. Krodkiewski and T. Jintanawan, “Stability improvement of periodic vibration of multi-degree-of-freedom systems by means of time-delay control,” in Proceedings of the International Conference on Vibration, Noise and Structural Dynamics, vol. 1, pp. 340–351, Venice, Italy, April, 1999.
16. Z. H. Wang and H. Y. Hu, “Calculation of the rightmost characteristic root of retarded time-delay systems via Lambert W function,” to appear in Journal of Sound and Vibration.
17. R. M. Corless, G. H. Gonnet, D. E. G. Hare, D. J. Jeffrey, and D. E. Knuth, “On the Lambert $W$ function,” Advances in Computational Mathematics, vol. 5, no. 4, pp. 329–359, 1996.
18. C. Hwang and Y.-C. Cheng, “A note on the use of the Lambert $W$ function in the stability analysis of time-delay systems,” Automatica, vol. 41, no. 11, pp. 1979–1985, 2005.
19. H. Shinozaki and T. Mori, “Robust stability analysis of linear time-delay systems by Lambert $W$ function: some extreme point results,” Automatica, vol. 42, no. 10, pp. 1791–1799, 2006.
20. M. Y. Fu, A. W. Olbrot, and M. P. Polis, “Robust stability for time-delay systems: the edge theorem and graphical tests,” IEEE Transactions on Automatic Control, vol. 34, no. 8, pp. 813–820, 1989.
21. M. Y. Fu, A. W. Olbrot, and M. P. Polis, “The edge theorem and graphical tests for robust stability of neutral time-delay systems,” Automatica, vol. 27, no. 4, pp. 739–741, 1991.
22. W. Michiels and T. Vyhlídal, “An eigenvalue based approach for the stabilization of linear time-delay systems of neutral type,” Automatica, vol. 41, no. 6, pp. 991–998, 2005.
23. K. Engelborghs, T. Luzyanina, and G. Samaey, “DDE-BIFTOOL v2.0: a Matlab package for the computation analysis of delay differential equations,” Tech. Rep. TW220, Department of Computer Science, Katholieke Universiteit Leuven, Leuven, Belgium, 2001.