Solution operator approximations for characteristic roots of delay differential equations

Applied Numerical Mathematics 56 (2006) 305–317 www.elsevier.com/locate/apnum

Solution operator approximations for characteristic roots of delay differential equations Dimitri Breda Dipartimento di Matematica e Informatica, Università degli Studi di Udine, via delle Scienze 208, I-33100 Udine, Italy Available online 23 May 2005

Abstract In this paper a new method for the numerical computation of characteristic roots for linear autonomous systems of Delay Differential Equations (DDEs) is proposed. The new approach enlarges the class of methods recently developed (see [SIAM J. Numer. Anal. 40 (2002) 629; D. Breda, Methods for numerical computation of characteristic roots for delay differential equations: experimental comparison, in: BIOCOMP2002: Topics in Biomathematics and Related Computational Problems at the Beginning of the Third Millennium, Vietri, Italy, 2002, Sci. Math. Jpn. 58 (2) pp. 377–388; D. Breda, The infinitesimal generator approach for the computation of characteristic roots for delay differential equations using BDF methods, Research Report RR2/2002, Department of Mathematics and Computer Science, Università di Udine, Italy, 2002; IMA J. Numer. Anal. 24 (2004) 1; SIAM J. Sci. Comput. (2004), in press]) and in particular it is based on a Runge–Kutta (RK) time discretization of the solution operator associated with the system. Hence this paper revisits the Linear Multistep (LMS) approach presented in [SIAM J. Numer. Anal. 40 (2002) 629] for the multiple discrete delay case and moreover extends it to the distributed delay case. We prove that the method converges with the same order as the underlying RK scheme and illustrate this with some numerical tests that are also used to compare the method with other existing techniques. © 2005 IMACS. Published by Elsevier B.V. All rights reserved. Keywords: Solution operator semigroup; Eigenvalue problem; Runge–Kutta methods

E-mail address: [email protected] (D. Breda). 0168-9274/$30.00 © 2005 IMACS. Published by Elsevier B.V. All rights reserved. doi:10.1016/j.apnum.2005.04.010

306

D. Breda / Applied Numerical Mathematics 56 (2006) 305–317

1. Introduction The interest for Delay Differential Equations (DDEs) is still growing in many application fields, since the introduction of a delay into a model enriches its dynamics and allows a closer and more accurate description of the real-life phenomena. In this work we consider a system of m-dimensional linear (or linearized) DDEs with multiple discrete and distributed delays:   −τ  l−1 k   Ll y(t − τl ) + Ml (θ )y(t + θ ) dθ , t  0, (1) y (t) = L0 y(t) + l=1

−τl

where L0 , L1 , . . . , Lk ∈ Cm×m , 0 = τ0 < τ1 < · · · < τk = τ and Ml : [−τl , −τl−1 ] → Cm×m , l = 1, . . . , k, is sufficiently smooth. Delay systems such as (1) are particularly important in control theory, where the stability effects of delays is a crucial problem (see the recent survey [17]). In the exhaustive monograph [16] many graphical and analytical tests to this aim are reported, but they work only for particular classes of system (1). This motivates our interest in the numerical computation of the characteristic roots of (1), i.e. the (infinitely many) roots of   det (λ) = 0, (2) where (λ) = λI − L0 −

 k 

−τ  l−1

l=1

−τl

Ll e−λτl −



Ml (θ )eλθ dθ ,

λ ∈ C,

(3)

since it is well known that the zero solution of (1) is asymptotically stable if and only if these roots have strictly negative real part. Recently numerical approaches for characteristic roots computation have been proposed, which are based on the discretization of either the solution operator associated with (1) or the infinitesimal generator of the solution operator semigroup. We briefly recall (see [7]) that the solution operator T (t), t  0, associated with (1) is defined by T (t)ϕ = yt ,

ϕ ∈ X,

(4)

where X = C([−τ, 0], C ) endowed with the maximum norm, yt is the function m

yt (θ ) = y(t + θ ),

θ ∈ [−τ, 0],

and y is the solution of (1) with initial data ϕ ∈ X. The family {T (t)}t0 is a C0 -semigroup with infinitesimal generator A : D(A) ⊆ X → X given by Aϕ = ϕ  , with domain

ϕ ∈ D(A),

(5)

 



D(A) = ϕ ∈ X | ϕ ∈ X and ϕ (0) = L0 ϕ(0) + The two following important results [7,13]

 k 

−τ  l−1

l=1

−τl

Ll ϕ(−τl ) +

Ml (θ )ϕ(θ ) dθ

 .

D. Breda / Applied Numerical Mathematics 56 (2006) 305–317

  1 det (λ) = 0 ⇐⇒ λ = ln μ, t   det (λ) = 0 ⇐⇒ λ ∈ σ (A)

307

  μ ∈ σ T (t) \ {0}, (6)

where σ (·) denotes the spectrum, suggest the idea of turning the characteristic roots approximation problem into a corresponding eigenvalue problem for a suitable matrix discretization of either T (t) (i.e., solution operator approach) or A (i.e., infinitesimal generator approach). Engelborghs and Roose propose in [10,11] the solution operator approach via LMS time integration of system (1) without any distributed delay term. Their method computes approximations to the roots via a large, standard and sparse eigenvalue problem and it is implemented in the MATLAB package DDE-BIFTOOL for DDEs bifurcation analysis [8,9]. The distributed delay case is considered in [15]. The complete development of the infinitesimal generator approach first appears in [4,2] where a matrix approximation to A is obtained by discretizing the derivative in (5) using an RK or LMS method, respectively, ending up with a large and sparse eigenvalue problem as in [11]. Recently a similar approach has been developed via pseudospectral differencing methods in [5]. In this paper we revisit the solution operator approach applying particular RK methods for the time integration of (1) and extending to the general case the method proposed in [1] for the single delay case and prove convergence of the approximated roots to the exact ones. Indeed, in this work we apply Radau-IIA methods, but the approach can be developed also for general RK methods. The work is structured as follows. Section 2 describes the approach in the single delay case, extensions to the multiple and distributed delay case are presented in Section 3 and 4, respectively. Section 5 is about convergence analysis of the method. Numerical results illustrate the properties of the method and a comparison with other techniques in Section 6.

2. The single delay case Let us consider the system of DDEs (1) in the single delay case: y  (t) = L0 y(t) + L1 y(t − τ ),

t  0,

(7)

with initial data y(t) = ϕ(t), t ∈ [−τ, 0]. For fixed N , N positive integer, let us consider the constant stepsize mesh on the interval [−τ, 0] given by ⎧ Ω = {θj + ci h, j = 1, . . . , N, i = 1, . . . , s}, ⎪ ⎨ N θj = −j h, j = 1, . . . , N, τ ⎪ ⎩h = N, 0 < c1 < · · · < cs = 1, and replace the continuous space X by the discrete space XN = (Cm )ΩN ∼ = CmsN . Applying to (7) an s-stage RK method (A, b, c) of order p such that (1) 0 < c1 < · · · < cs = 1, (2) A is invertible, (3) b = (as1 , . . . , ass )T , where we take as points c1 , . . . , cs in ΩN the abscissae of the RK method, we can obtain an approximation x (1) ∈ XN to the solution y of (7) on the interval [h − τ, h]. Let us set x (δ) = ((x1(δ) )T , . . . , (xN(δ) )T )T ,

308

D. Breda / Applied Numerical Mathematics 56 (2006) 305–317

(δ) T (δ) T T δ = 0, 1, where xj(δ) = ((xj,1 ) , . . . , (xj,s ) ) is the ms-vector approximating the solution at the grid (δ) points δh + θj + ci h, j = 1, . . . , N , i = 1, . . . , s, of the discrete interval δh + ΩN , i.e. xj,i approximates (0) y(δh + θj + ci h) and in particular x corresponds to the discretization on ΩN of the initial data ϕ for (7). Within this setting, x (1) is given by  x1(1) = R(hL0 )(1s esT ⊗ Im )x1(0) + hR(hL0 )(A ⊗ L1 )xN(0) , j = 2, . . . , N, xj(1) = xj(0) −1 ,

where R(Z) = (Ims − A ⊗ Z)−1 ,

Z ∈ Cm×m ,

es = (0, . . . , 0, 1)T ∈ Rs and 1s = (1, . . . , 1)T ∈ Rs . We can thus write x (1) = SN x (0) ,

(8)

where SN is the msN × msN matrix ⎛ R(hL0 )(1s esT ⊗ Im ) 0 ⎜ Ims 0 ⎜ ⎜ 0 I ms SN = ⎜ .. .. ⎝ . . 0 0

... ... ... .. .

0 0 0 .. .

...

Ims

⎞ hR(hL0 )(A ⊗ L1 ) ⎟ 0 ⎟ ⎟. 0 ⎟ .. ⎠ .

(9)

0

It follows that (8) is the discrete counterpart of (4) at time t = h, i.e. SN is the discretization of T (h). Thus from (6) λ=

1 ln μ h

is an approximation to a characteristic root of (7) whenever μ = 0 is an eigenvalue of SN .

3. Extension to the multiple delay case As a generalization of (7) let us consider the system of DDEs: 

y (t) = L0 y(t) +

k 

Ll y(t − τl ),

t  0,

(10)

l=1

with initial data y(t) = ϕ(t), t ∈ [−τ, 0]. Application of the RK method (A, b, c) to (10), based on the same mesh ΩN , requires a suitable approximation zl,i to y(ci h − τl ), i = 1, . . . , s, l = 1, . . . , k − 1, since there is no reason why the points ci h − τl should belong to the discrete set ΩN or h + ΩN . Setting T T T , . . . , zl,s ) ∈ Cms , l = 1, . . . , k − 1, we get zl = (zl,1   k−1   (0)  T (1) (0) (11) (A ⊗ Ll )zl + h(A ⊗ Lk )xN . x1 = R(hL0 ) 1s es ⊗ Im x1 + h l=1

D. Breda / Applied Numerical Mathematics 56 (2006) 305–317

309

In order to preserve the convergence order p of the RK method, such an approximation zl,i can be constructed by centered Lagrange interpolation of p + 1 values at the nodes around the point ci h − τl . Therefore, for each of these points, we find the index αl,i := γ such that θγ  ci h − τl < θγ −1 and set p/2

zl,i =



lr (ci h − τl )xα(0) , l,i −r,s

(12)

r=1−p/2

where lr (·) is the rth Lagrange coefficient computed on the p + 1 grid points θαl,i −r , r = 1 − p/2 , . . . , p/2 , around ci h − τl (q denotes the smallest integer n such that n  q). The solution operator is now approximated by (9) with only the first block-row modified according to (12) and (11) and the characteristic roots approximations are recovered as in the single delay case. Remark 1. When the steplength h is too large or some discrete delays τl , l = 1, . . . , k − 1, are too close to the extremes 0 or τ it is not possible to use a centered interpolation. In this cases interpolation must be based on the p + 1 closest available node values. In any case, to preserve the order p of the RK method, one must always choose N  p + 1. 4. Extension to the distributed delay case As a further extension of (10) let us consider system (1) with initial data y(t) = ϕ(t), t ∈ [−τ, 0]. Before applying the RK method to (1), we substitute an appropriate quadrature rule for the integral term. Since we have to preserve the order p of the RK method and to maintain the computational effort as low as possible, the choice falls on a composite Gauss–Legendre quadrature formula based on p/2 nodes ξd ∈ [−1, 1] and weights ωd , d = 1, . . . , p/2 . This ensures an order of accuracy of at least p for the l−1 . quadrature on each delay interval [−τl , −τl−1 ], l = 1, . . . , k, with steplength Hl = τl −τ N Now, in a similar manner to the multiple delay case, the application of the RK method to (1), based on the same mesh ΩN , requires a suitable approximation wl,q,d,i of y(ci h + ψl,q,d ), i = 1, . . . , s, where ψl,q,d is the translation of ξd ∈ [−1, 1] into the interval [−τl + (q − 1)Hl , −τl + qHl ], l = 1, . . . , k, T T , . . . , wl,q,d,s )T ∈ Cms , we get q = 1, . . . , N . Setting wl,q,d = (wl,q,d,1  k−1   T  (0) (1) (A ⊗ Ll )zl + h(A ⊗ Lk )xN(0) x1 = R(hL0 ) 1s es ⊗ Im x1 + h l=1

+h

N k   l=1 q=1

 Hl   ωd A ⊗ Ml (ψl,q,d ) wl,q,d . 2 d=1 p/2

(13)

As in the previous section (Remark 1 included), for each point ci h + ψl,q,d we find the index βl,q,d,i := γ such that θγ  ci h + ψl,q,d < θγ −1 and set p/2

wl,q,d,i =



lr (ci h + ψl,q,d )xβ(0) . l,q,d,i −r,s

(14)

r=1−p/2

Once more the solution operator approximation (9) changes only in its first block-row according to (14) and (13) and the characteristic roots are again recovered by the logarithmic transformation of the eigenvalues of SN .

310

D. Breda / Applied Numerical Mathematics 56 (2006) 305–317

5. Convergence analysis We now develop an analysis of the convergence order of the computed roots to the exact ones. We first derive an expression for the discrete version of the characteristic equation (2), i.e. for the characteristic equation of the matrix SN . 5.1. The discrete characteristic equation We start from the eigenvalue problem for SN : μ is an eigenvalue of SN if and only if there exists x (0) ∈ XN \ {0} such that SN x (0) = μx (0) ,

(15)

which means by (9), (12)–(14)  s k−1   (0) (0) (0) bi L0 μx1,i + Ll μx1,s = x1,s + h

p/2

i=1

r=1−p/2

+

N k   l=1 q=1

l=1



(0) lr (ci h − τl )xα(0) + Lk xN,i l,i −r,s

p/2

p/2  Hl  ωd Ml (ψl,q,d ) lr (ci h + ψl,q,d )xβ(0) l,q,d,i −r,s 2 d=1 r=1−p/2

(16)

and (0) = xj(0) μxj,i −1,i ,

j = 2, . . . , N, i = 1, . . . , s.

(17)

(0) = μci +θj / h v, v ∈ Cm \ {0}, is a system of Now, observe that μ ∈ C is an eigenvalue of SN if and only if xj,i fundamental solutions of the system of difference equations (16) and (17). Since the characteristic roots (0) of (1) are approximated by λ = h1 ln μ, where μ is an eigenvalue of SN , substituting xj,i = eλ(ci h+θj ) v in (16) and (17) leads to  s  λh λ(cs h+θ1 ) λ(cs h+θ1 ) v=e v+h bi L0 eλh eλ(ci h+θ1 ) v e e i=1

+

k−1  l=1

+

p/2



Ll

lr (ci h − τl )eλ(cs h+θαl,i −r ) v + Lk eλ(ci h+θN ) v

r=1−p/2

p/2 N k   Hl  l=1 q=1

2



p/2

ωd Ml (ψl,q,d )

d=1



lr (ci h + ψl,q,d )e

λ(cs h+θβl,q,d,i −r )

r=1−p/2

and eλh eλ(ci h+θj ) v = eλ(ci h+θj −1 ) v,

j = 2, . . . , N, i = 1, . . . , s,

where the last is clearly true. Eq. (18) can be rewritten as  −τ  l−1 s k k     p+1  λh λci h −λτl λθ , e v=v+h bi e Ll e v+ Ml (θ )e v dθ + O h L0 v + i=1

l=1

l=1 −τ l

v

(18)

D. Breda / Applied Numerical Mathematics 56 (2006) 305–317

311

where the term O(hp+1 ) arises (in part) from the interpolation error in approximating the exponential function at the points ci h − τl , i = 1, . . . , s, l = 1, . . . , k − 1, and ci h + ψl,q,d , i = 1, . . . , s, l = 1, . . . , k, q = 1, . . . , N , d = 1, . . . , p/2 and (in part) from the error in the quadrature of the distributed term. Looking for nontrivial solutions, i.e. v ∈ Cm \ {0}, leads to the discrete characteristic equation   (19) det N (λ) = 0, where   −τ  l−1 k    eλh − 1 −λτl λθ Im − L0 − N (λ) = λ + Ml (θ )e dθ + O hp+1 Ll e ϕ(λ, h) l=1

(20)

−τl

with ϕ(λ, h) = λh

s 

bi eλci h ,

λ ∈ C.

i=1

5.2. Convergence order analysis Let us introduce for the exact and the approximated characteristic equations associated with (1) the functions d, dN : C → C given by   d(λ) = det (λ) and

  dN (λ) = det N (λ) ,

where (λ) and N (λ) are given by (3) and (20), respectively. Moreover let us denote by B(λ, ρ) the closed ball in C with center λ and radius ρ. Before giving the convergence result we first prove two useful lemmas concerning the functions d and dN above. Lemma 2. Let λ∗ ∈ C and ρ0 > 0. There exists h0 > 0 such that for h < h0 and for all λ ∈ B(λ∗ , ρ0 )   d(λ) − dN (λ)  C0 hp with C0 = C0 (λ∗ , ρ0 ) and p is the order of the RK method applied. Proof. Subtracting (3) from (20) we obtain  λh    e −1 − 1 Im + O hp+1 . N (λ) − (λ) = λ ϕ(λ, h) Now observe that

  ∞ s ∞      λh λk+1 hk+1 (k + 1) bi cik − 1 = dk+1 hk+1 ϕ(λ, h) − e − 1 = (k + 1)! k=p k=0 i=1

312

D. Breda / Applied Numerical Mathematics 56 (2006) 305–317

with dp+1 = 0 since the RK method has convergence order p and thus the following order conditions hold (see for instance [14]): s  1 , k = 0, . . . , p − 1. bi cik = k+1 i=1 Therefore for ρ0 > 0 there exists h0 > 0 such that for h < h0 and for all λ ∈ B(λ∗ , ρ0 )   N (λ) − (λ)  Khp where  ·  is a matrix norm on Cm×m and K is independent of h and λ. Since the derivative of det(X), X ∈ Cm×m , with respect to X is continuous the assertion follows. 2 Lemma 3. If the function d(λ) has a zero λ∗ ∈ C with multiplicity ν, then there exists ρ1 = ρ1 (λ∗ ) > 0 such that   d(λ) > C1 |λ − λ∗ |ν for λ ∈ B(λ∗ , ρ1 ) \ {λ∗ } where C1 = C1 (λ∗ ). Proof. Expanding d(λ) around λ∗ and taking into account the multiplicity ν we obtain   d (ν) (λ∗ ) (λ − λ∗ )ν + O |λ − λ∗ |ν+1 ν! with d (ν) (λ∗ ) = 0. The result follows. 2 d(λ) =

Now let us state the main result. Theorem 4. If (2) has a characteristic root λ∗ ∈ C with multiplicity ν, then there exists an h1 > 0 such that for h < h1 (19) has exactly ν roots λi , i = 1, . . . , ν (taking into account multiplicities) and   max |λ∗ − λi | = O hp/ν 1iν

where p is the order of the RK method applied. Proof. Let h1 > 0 such that h1  h0 and  1/ν C0 p/ν h1  min{ρ0 , ρ1 } C1 where h0 , C0 and ρ0 from Lemma 2 and C1 and ρ1 from Lemma 3. For h < h1 set  1/ν C0 hp/ν ρ= C1 Lemmas 2 and 3 yield     d(λ) > dN (λ) − d(λ) for |λ − λ∗ | = ρ. Hence, by Rouché’s Theorem (see for instance [6, Section 3, Theorem 3.8]), d(λ) and dN (λ) have exactly the same number of roots counted with their multiplicities in B(λ∗ , ρ). 2

D. Breda / Applied Numerical Mathematics 56 (2006) 305–317

313

In the next proposition we also show the nonexistence of “ghost roots”, i.e. if there exists a sequence of approximated characteristic roots that converges, then it converges to an exact characteristic root. Proposition 1. Let {ΩN (i) }i1 be a sequence of meshes on the interval [−τ, 0] such that N (i) → ∞ as i → ∞. If λ(i) → λ∗ as i → ∞, where λ(i) is a root of dN (i) (λ), then λ∗ is a root of d(λ). Proof. Since     ∗            d(λ ) = dN (i) λ(i) − d(λ∗ )  dN (i) λ(i) − d λ(i)  + d λ(i) − d(λ∗ ), the assertion follows by Lemma 2 and by continuity of d.

2

6. Numerical results The method proposed is implemented with s-stage Radau-IIA schemes of order p = 2s − 1 (s = 1, 2 and 3 are used), which satisfy the conditions given in Section 2. All the results concern the following systems and equations: Example 5.

  y  (t) = 2 − e−2 y(t) + y(t − 1).

(21)

Example 6 (from [10]).   y1 (t) = −0.5y1 (t) − tanh (y1 (t − 1.57)) + tanh (y2 (t − 0.2)), y2 (t) = −0.5y2 (t) + 2.34 tanh (y1 (t − 0.2)) − tanh (y2 (t − 1.57)).

(22)

Example 7 (from [12]). −0.1 −0.5   M1 y(t + θ ) dθ + M2 y(t + θ ) dθ, y (t) = L0 y(t) + L1 y(t − 1) + 

−0.3

(23)

−1

with coefficient matrices     −3 1 1 0 , L1 = , L0 = −24.646 −35.430 2.35553 2.00365     2 2.5 −1 0 , M2 = . M1 = 0 −0.5 0 −1 Eq. (21) has an exact rightmost root λ∗ = 2 while system (22), linearized around the steady state solution (y1∗ , y2∗ ) = (0, 0), has a rightmost root λ∗  0.3475. System (23) has a rightmost root λ∗  −1.2462. Fig. 1 refers to the single delay equation (21) and shows (left) the rightmost root convergence behavior of the proposed method (with p = 5), which agrees with the result stated in Theorem 4. Moreover the method is compared with other techniques, in particular the solution operator approach with LMS

314

D. Breda / Applied Numerical Mathematics 56 (2006) 305–317

Fig. 1. Rightmost root absolute error vs N (left) and computational time (seconds, right) vs required tolerance on the maximum absolute error TOL for system (21) using the solution operator approach based on RK (•) and BDF (∗) methods and the infinitesimal generator approach based on RK (◦) and pseudospectral differencing (+) methods.

Table 1 First ten rightmost roots value and modulus for system (22) Characteristic root

Modulus

0.34748172572643 −0.08116698520245 −0.43412304132027 ± i1.62752128215701 −0.82061576063772 ± i5.11180446782392 −1.20975894908637 ± i4.82695048752634 −1.24622687378869 ± i8.90021501181034

0.34748172572643 0.08116698520245 1.68442522507803 5.17725362947110 4.97624032015389 8.98704115253777

methods (5th-order BDF, see [1]) and the infinitesimal generator approach with RK (5th-order RadauIIA, see [4]) and pseudospectral differencing methods (based on Chebyshev points, see [5]). The linear (in a double-logarithmic scale context) convergence of the first three methods is clearly surpassed by the superlinear behavior of the pseudospectral technique for very low tolerances. The advantage is even more evident in terms of computational time (right). We then compare the solution operator approach based on RK methods and the infinitesimal generator approach based on pseudospectral differencing methods computing the first ten rightmost roots

D. Breda / Applied Numerical Mathematics 56 (2006) 305–317

315

Fig. 2. First ten rightmost roots of system (22).

Fig. 3. First ten rightmost roots absolute error vs N (left) and computational time (seconds, right) vs required tolerance on the maximum absolute error TOL for system (22) using the solution operator approach based on RK methods (1st row) and the infinitesimal generator approach based on psuedospectral differencing methods (2nd row).

316

D. Breda / Applied Numerical Mathematics 56 (2006) 305–317

Fig. 4. Rightmost root absolute error vs N (left) and computational time (seconds, right) vs required tolerance on the maximum absolute error TOL for system (23) using the solution operator approach based on RK methods (•: p = 5; ∗: p = 3; ◦: p = 1) and the infinitesimal generator approach based on psuedospectral differencing methods (+).

of system (22), listed in Table 1 and shown in Fig. 2. Fig. 3 (left) shows the error trend for each root (simple-real or complex-conjugate) in ascending order with respect to its modulus. This last plays a fundamental role in the sense that approximated roots of smaller modulus converge faster: for the pseudospectral approach this phenomenon clearly results from Theorem 4 in [5] and an analogous result can be deduced from the convergence analysis carried out in Section 5. Fig. 3 (right) demonstrates once more the superiority of the pseudospectral approach in terms of computational time with respect to the required tolerance for the maximum absolute error over all the computed roots. The above methods are also compared for computing the rightmost root of system (23). The convergence of order p = 2s − 1 is confirmed for all s-stage Radau-IIA methods even when distributed terms are concerned (see Fig. 4). The pseudospectral approach is still advantageous, especially in terms of computational time. Note that in the discrete delay case, if all the delays are commensurate, the computational time can be reduced by choosing a steplength that exactly divides the smallest delay interval. In fact this avoids the use of interpolation to recover the solution in the past since delayed grid points correspond to past grid points. Other results and comparisons can be found in [1,3] where a more detailed discussion on implementation issues is reported.

D. Breda / Applied Numerical Mathematics 56 (2006) 305–317

317

References [1] D. Breda, Methods for numerical computation of characteristic roots for delay differential equations: experimental comparison, in: BIOCOMP2002: Topics in Biomathematics and Related Computational Problems at the Beginning of the Third Millennium, Vietri, Italy, 2002, Sci. Math. Jpn., 58 (2), pp. 377–388. [2] D. Breda, The infinitesimal generator approach for the computation of characteristic roots for delay differential equations using BDF methods, Research Report RR2/2002, Department of Mathematics and Computer Science, Università di Udine, Italy, 2002. [3] D. Breda, Collection of numerical tests on the computation of characteristic roots for delay differential equations, Research Report RR2/2002, Department of Mathematics and Computer Science, Università di Udine, Italy, 2002. [4] D. Breda, S. Maset, R. Vermiglio, Computing the characteristic roots for delay differential equations, IMA J. Numer. Anal. 24 (1) (2004) 1–19. [5] D. Breda, S. Maset, R. Vermiglio, Pseudospectral differencing methods for characteristic roots of delay differential equations, SIAM J. Sci. Comput. (2004), in press. [6] J.B. Conway, Functions of One Complex Variable, second ed., Graduate Texts in Math., vol. 11, Springer, New York, 1978. [7] O. Diekmann, S.A. van Gils, S.M. Verduyn Lunel, H.O. Walther, Delay Equations—Functional, Complex and Nonlinear Analysis, Appl. Math. Sci., vol. 110, Springer, New York, 1995. [8] K. Engelborghs, T. Luzyanina, D. Roose, Numerical bifurcation analysis of delay differential equations using DDEBIFTOOL, ACM Trans. Math. Software 28 (1) (2000) 1–21. [9] K. Engelborghs, T. Luzyanina, G. Samaey, DDE-BIFTOOL v. 2.00: A Matlab package for bifurcation analysis of delay differential equations, Report TW330, Department of Computer Science, K. U. Leuven, Belgium, 2001. [10] K. Engelborghs, D. Roose, On stability of LMS methods and characteristic roots of delay differential equations, SIAM J. Numer. Anal. 40 (2) (2002) 629–650. [11] K. Engelborghs, D. Roose, Numerical computation of stability and detection of Hopf bifurcations of steady-state solutions of delay differential equations, Adv. Comput. Math. 10 (3–4) (1999) 271–289. [12] A. Fattouh, O. Sename, J.M. Dion, H∞ controller and observer design for linear systems with point and distributed delays, in: Proceedings of 2nd IFAC Workshop on Linear Time Delay Systems, Ancona, Italy, 2000. [13] J.K. Hale, S.M. Verduyn Lunel, Introduction to Functional Differential Equations, Appl. Math. Sci., vol. 99, Springer, New York, 1993. [14] J.D. Lambert, Numerical Methods for Ordinary Differential Systems, Wiley, Chichester, UK, 1997. [15] T. Luzyanina, K. Engelborghs, D. Roose, Computing stability of differential equations with bounded distributed delays, Numer. Algorithms 34 (1) (2003) 41–66. [16] S.I. Niculescu, Delay Effects on Stability: A Robust Control Approach, Lecture Notes in Control and Inform. Sci., vol. 269, Springer, New York, 2003. [17] J.P. Richard, Time-delay systems: an overview of some recent advances and open problems, Automatica 39 (2003) 1667– 1694.