PSEUDOSPECTRAL DIFFERENCING METHODS ... - Semantic Scholar
PSEUDOSPECTRAL DIFFERENCING METHODS FOR CHARACTERISTIC ROOTS OF DELAY DIFFERENTIAL EQUATIONS D. BREDA† §, S. MASET‡ AND R. VERMIGLIO† ¶ Abstract. In [BMV03]and [Bre02] the authors proposed to compute the characteristic roots of Delay Differential Equations (DDEs) with multiple discrete and distributed delays by approximating the derivative in the infinitesimal generator of the solution operator semigroup by Runge-Kutta (RK) and Linear Multistep (LMS) methods, respectively. In this work the same approach is proposed in a new version based on pseudospectral differencing techniques. We prove the “spectral accuracy” convergence behavior typical of pseudospectral schemes, as also illustrated by some numerical experiments. Key words. Pseudospectral methods, delay differential equations, characteristic roots. AMS subject classifications. 65J10, 65Q05, 34K06, 34K20.
1. Introduction. The interest for delay differential equations (DDEs) is still growing in many application fields, since the introduction of the delay in the model enriches its dynamics and allows a closer and more accurate description of the real-life phenomena. In this paper we consider a system of m-dimensional linear (or linearized) DDEs with multiple discrete and distributed delays:
(1.1)
y (t) = L0 y(t) +
k
0 Ll y(t − τl ) +
l=1
M (θ)y(t + θ)dθ,
t ≥ 0,
−τ
where L0 , L1 , . . . , Lk ∈ Cm×m , 0 = τ0 < τ1 < . . . < τk = τ and M : [−τ, 0] → Cm×m is a piecewise smooth function. Indeed in the implementation of the method we write the distributed term as 0 (1.2) −τ
−τ l−1 k M (θ)y(t + θ)dθ = Ml (θ)y(t + θ)dθ,
t ≥ 0,
l=1 −τ l
where every function Ml , l = 1, . . . , k, is smooth. Delay systems such as (1.1) are particularly important in control theory, where the stability effects of delays is a crucial problem (see the recent survey [Ric03]). In the exhaustive monograph [Nic00] many graphical and analytical tests to this aim are reported, but they work only for particular classes of system (1.1). This motivates our interest in the numerical computation of the characteristic roots of (1.1), i.e. the (infinitely many) roots of (1.3)
det(∆(λ)) = 0,
† Dipartimento di Matematica e Informatica, Universit` a degli Studi di Udine, Via delle Scienze 208, I-33100 Udine, Italy. ‡ Dipartimento di Matematica e Informatica, Universit` a degli Studi di Trieste, Via Valerio 12, I-34127 Trieste, Italy, ([email protected]). § ([email protected]). ¶ ([email protected]).
1
2
D. BREDA, S. MASET AND R. VERMIGLIO
where ∆(λ) = λI −
(1.4)
k
Ll e
−λτl
0 −
l=0
M (θ)eλθ dθ, λ ∈ C, −τ
since it is well known that the zero solution of (1.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 to (1.1) or the infinitesimal generator of the solution operator semigroup. We briefly recall that the solution operator T (t), t ≥ 0, associated to (1.1) is defined by T (t)ϕ = yt ,
ϕ ∈ X,
where X = C ([−τ, 0] , Cm ) endowed with the maximum norm, yt is the function yt (θ) = y(t + θ),
θ ∈ [−τ, 0],
and y is the solution of (1.1) with initial data ϕ ∈ X. The family {T (t)}t≥0 is a C0 -semigroup with infinitesimal generator A : D(A) ⊆ X → X given by Aϕ = ϕ , ϕ ∈ D (A) ,
(1.5) with domain (1.6) D (A) =
ϕ ∈ X | ϕ ∈ X and ϕ (0) =
k l=0
0 Ll ϕ(−τl ) +
M (θ)ϕ(θ)dθ −τ
.
So (1.1) can be restated as the following abstract Cauchy problem ([DGVLW95]): dyt t>0 dt = Ayt , . y0 = ϕ, The two following important results ([DGVLW95], [HVL93]) 1. det(∆(λ)) = 0 ⇔ λ = 1t ln µ, µ ∈ σ(T (t)) \ {0}; 2. det(∆(λ)) = 0 ⇔ λ ∈ σ(A); where σ(·) denotes the spectrum, suggest the idea to turn the characteristic roots approximation problem into a corresponding eigenvalue problem for suitable matrix discretization of either T (t) (i.e. solution operator approach) or A (i.e. infinitesimal generator approach). Engelborghs and Roose propose in [ER02] and [ER99] the solution operator approach via LMS time integration for system (1.1) without distributed delay term. Their method computes approximations to the roots from a large, standard and sparse eigenvalue problem and it is implemented in the MATLAB package DDE-BIFTOOL for DDEs bifurcation analysis ([ELR02] and [ELS01]). The distributed delay case is considered in [LER03] by using LMS methods and in [Bre03] by using RK methods. The complete development of the infinitesimal generator approach first appears in [BMV03] and [Bre02] where a matrix approximation to A is obtained discretizing the derivative in (1.5) by RK or LMS method, respectively, ending up with a large and sparse eigenvalue problem as in [ER99].
PSEUDOSPECTRAL METHODS FOR CHARACTERISTIC ROOTS
3
This work continues the investigation along this line and proposes the infinitesimal generator approach via pseudospectral differencing methods (see [Tre00] and the bibliography therein). This technique involves the exact differentiation of interpolants at selected sets of points. The resulting differentiation matrix is non-sparse, but we can take advantage of the well-known “spectral accuracy” to obtain very accurate approximation with small matrix dimension. This behavior represents in fact, for sufficiently small tolerance, the outstanding advantage of this method compared to the previously cited discretization schemes. The paper is structured as follows. Section 2 describes the pseudospectral approach. Section 3 is about convergence analysis of the method. In §4 we discuss about implementation issues. In §5 we perform numerical experiments on test delay systems of type (1.1) to illustrate the properties of the new method and to critically compare it with LMS and RK discretization techniques. 2. The pseudospectral approach. To simplify the notation, we rewrite the system of DDEs (1.1) as follows: y (t) = f (yt ) where (2.1)
f (ϕ) =
k l=0
0 Ll ϕ(−τl ) +
M (θ)ϕ(θ)dθ,
ϕ ∈ X.
−τ
Given a positive integer N , let us consider the mesh ΩN of N + 1 distinct points in [−τ, 0] ΩN = {θN,i , i = 0, 1, . . . , N } with 0 = θN,0 > θN,1 > . . . > θN,N ≥ −τ . We replace the continuous space X by the space XN of the discrete functions defined on the mesh ΩN , i.e. any ϕ ∈ X is discretized into the block-vector x ∈ XN of components xi = ϕ(θN,i ) ∈ Cm , i = 0, 1, . . . , N . Let LN x, x ∈ XN , be the unique Cm -valued interpolating polynomial of degree ≤ N with (LN x)(θN,i ) = xi , i = 0, 1, . . . , N . Thus we approximate the infinitesimal generator A by the matrix AN : XN → XN , called spectral differentiation matrix, defined as follows: (AN x)0 = fN (LN x) (2.2) (AN x)i = (LN x) (θN,i ), i = 1, 2, . . . , N where fN is an approximation of f in which the distributed delay integral term in (2.1) is substituted by a suitable interpolatory quadrature rule or fN is equal to f in the case of simple matrix function M for which the integral can be exactly computed. By using the Lagrange representation of LN x we obtain a0 a1 · · · aN d10 d11 · · · d1N m(N +1)×m(N +1) AN = . .. . . .. ∈ C .. . . . dN 0 dN 1 · · · dN N with (2.3)
aj = fN (lN,j (·)Im ) , j = 0, 1, . . . , N (θN,i )Im , i = 1, . . . , N, j = 0, 1, . . . , N dij = lN,j
4
D. BREDA, S. MASET AND R. VERMIGLIO
where lN,j , j = 0, 1, . . . , N , are the Lagrange coefficients relevant to the nodes in m ΩN and, by definition, fN (lN,j (·)Im ) = fN lN,j (·)e(i) i=1 with e(i) ’s the canonical vectors in Rm . Explicit expressions of the dij ’s for particular choices of ΩN (e.g. Chebyshev extremal points) can be found in [Tre00, §6]. So finally the original problem of the numerical computation of the characteristic roots of (1.1) can be turned into the eigenvalue problem for the matrix AN , i.e. eigenvalues of AN directly approximate the characteristic roots. How accurate are these approximations? In the next section we develop the convergence analysis. 3. Convergence analysis. The convergence analysis of the approximated roots to the exact ones is based on a comparison between the continuous (i.e. exact) characteristic equation (1.3) and its discrete (i.e. approximated) counterpart. Let us thus characterize the continuous and the discrete eigenvalue problems for A and AN , respectively. Let λ ∈ C and ϕ ∈ D (A) be such that Aϕ = λϕ, that means
ϕ (θ) = λϕ(θ), ϕ (0) = f (ϕ) .
θ ∈ [−τ, 0]
Hence λ is eigenvalue of A if and only if there exists u ∈ Cm , u = 0, such that λu = f eλ· u . ˆ By introducing the linear operator A(λ) : Cm → Cm given by ˆ A(λ)u = f eλ· u , u ∈ Cm (3.1) we get that λ is eigenvalue of A if and only if ˆ det λIm − A(λ) (3.2) =0 which is the continuous characteristic equation (1.3). In the same fashion, let x ∈ XN be such that AN x = λx, thus by (2.2) we get fN (LN x) = λx0 and (3.3)
(LN x) (θN,i ) = λxi ,
i = 1, 2, . . . , N.
By defining pN (·; λ, u), u ∈ Cm , as the collocation polynomial relevant to the nodes θN,1 , θN,2 , . . . , θN,N for the initial value problem y (θ) = λy(θ), θ ∈ [−τ, 0] (3.4) y(0) = u,
PSEUDOSPECTRAL METHODS FOR CHARACTERISTIC ROOTS
5
whose solution is the exponential y(·; λ, u) = eλ· u, one can easily derive that (3.3) is equivalent to (LN x)(·) = pN (·; λ, x0 ). Consequently λ is eigenvalue of AN if and only if there exists u ∈ Cm , u = 0, such that λu = fN (pN (·; λ, u)) . Finally, by introducing the linear operator AˆN (λ) : Cm → Cm given by (3.5)
AˆN (λ)u = fN (pN (·; λ, u)) ,
u ∈ Cm ,
we get that λ is eigenvalue of AN if and only if (3.6) det λIm − AˆN (λ) = 0 which is the discrete characteristic equation. Assumption 3.1. We assume to use the Chebyshev extremal nodes on [−τ, 0], i.e. π τ θN,i = (3.7) cos (i ) − 1 , i = 0, . . . , N. 2 N Let us denote by B(λ, ρ) the closed ball in C of center λ and radius ρ. In next lemma we give an upper bound of the error of the collocation polynomial for (3.4). Lemma 3.2. Let λ∗ ∈ C and ρ0 > 0. Under Assumption 3.1 there exists N0 ∈ N such that, for N ≥ N0 , for λ ∈ B(λ∗ , ρ0 ) and for u ∈ Cm , the collocation polynomial pN (·; λ, u) of (3.4) based on (3.7) exists and it is unique. Moreover it holds C0 max |pN (θ; λ, u) − y(θ; λ, u)| ≤ √ θ∈[−τ,0] N
C1 N
N |u|
where C0 = C0 (λ∗ , ρ0 ) and C1 = C1 (λ∗ , ρ0 ) = τ (|λ∗ | + ρ0 )e are constants independent of N . Proof. By introducing the (linear) integral Volterra operator Kλ : X → X as θ (Kλ ϕ)(θ) = λ
ϕ(s)ds,
ϕ ∈ X,
λ ∈ C,
θ ∈ [−τ, 0],
0
and the (linear and bounded) Lagrange interpolation operator LN −1 : X → X relevant to the nodes (3.7), we can rewrite (3.4) and its collocation equations, respectively, as the equations on the space X (3.8)
y = u + Kλ y
and (3.9)
pN = u + Kλ LN −1 pN
where y = y(·; λ, u), pN = pN (·; λ, u) and u denotes the function of X of constant value u. By subtracting (3.8) from (3.9) and setting eN = pN − y and rN = (LN −1 − I)y, we get the following equation for eN (3.10)
eN = Kλ LN −1 eN + Kλ rN .
6
D. BREDA, S. MASET AND R. VERMIGLIO
It can be directly shown that the solutions of (3.10) are the functions Kλ eˆN , where eˆN is a solution of eˆN = LN −1 Kλ eˆN + rN . Since I − Kλ is invertible and, by the known estimate for the Lagrange interpolation on Chebyshev nodes, LN −1 Kλ − Kλ → 0 as N → ∞, it follows that, for sufficiently large N (and independently of λ), I − LN −1 Kλ is invertible and (I − LN −1 Kλ )−1 ≤ 2 (I − Kλ )−1 . Hence eN = Kλ (I − LN −1 Kλ )−1 rN is the unique solution of (3.10) and as a consequence eN ≤ 2 Kλ (I − Kλ )−1 rN . Since rN is the error of the Lagrange interpolation of the exponential y(·; λ, u) = eλ· u relevant to the nodes (3.7), we conclude that rN ≤ max{1, e−(λ)τ }
(τ |λ|)N |u|. N!
The final assertion follows by application of the Stirling formula √ N ! ≥ 2πN (N/e)N . Now let us set the following notation relevant to the continuous and discrete characteristic equations (3.2) and (3.6) respectively: ˆ (3.11) d(λ) = det λIm − A(λ) dN (λ) = det λIm − AˆN (λ)
(3.12)
ˆ with A(λ) and AˆN (λ) defined in (3.1) and (3.5), respectively. The following assumption and lemmas are the next necessary steps to state the final theorem of convergence. Assumption 3.3. We assume that the function fN in (2.2) satisfies sup fN < +∞.
N ∈N
Lemma 3.4. Let λ∗ ∈ C and ρ0 > 0. Under Assumptions 3.1 and 3.3, there exists N0 ∈ N such that, for N ≥ N0 and λ ∈ B(λ∗ , ρ0 ), we have N C1 1 |d(λ) − dN (λ)| ≤ C2 εN + √ N N where εN =
sup λ∈B(λ∗ ,ρ0 )
|f (y(·; λ, u)) − fN (y(·; λ, u)) | , |u|
PSEUDOSPECTRAL METHODS FOR CHARACTERISTIC ROOTS
7
C1 = C1 (λ∗ , ρ0 ) as defined in Lemma 3.2 and C2 = C2 (λ∗ , ρ0 ) are suitable constants independent of N . Proof. By Lemma 3.2 : ˆ |A(λ)u − AˆN (λ)u| ≤ |f (y(·; λ, u)) − fN (y(·; λ, u)) | + + |fN (y(·; λ, u)) − fN (pN (·; λ, u)) | ≤ N C0 C1 |u| ≤ εN + sup fN √ N N N ∈N for N ≥ N0 and λ ∈ B(λ∗ , ρ0 ). Now, since the derivative of det(A) with respect to A ∈ Cm×m is continuous, the assertion easily follows. Lemma 3.5. If the function d(λ) has a zero λ∗ ∈ C with multiplicity ν, then there exists ρ1 = ρ1 (λ∗ ) > 0 such that |d(λ)| > C3 |λ − λ∗ |
ν
for λ ∈ B(λ∗ , ρ1 ) \ {λ∗ } where C3 = C3 (λ∗ ). Proof. Expanding d(λ) around λ∗ and taking into account the multiplicity ν we obtain d(λ) =
d(ν) (λ∗ ) ν+1 (λ − λ∗ )ν + O |λ − λ∗ | ν!
with d(ν) (λ∗ ) = 0. The result follows. We are now able to state the convergence result. Theorem 3.6 (convergence). Let λ∗ ∈ C be a zero of the function d(λ) with multiplicity ν, let ρ1 as in Lemma 3.5 and let (3.13)
ρN =
C2 C3
1/ν N 1/ν C1 1 εN + √ N N
with C1 , C2 as in Lemma 3.4 with ρ0 = ρ1 and C3 as in Lemma 3.5. Under Assumptions 3.1 and 3.3, if ρN ≤ ρ1 then the function dN (λ) has ν zeros λi , i = 1, . . . , ν, (taking into account multiplicities) such that (3.14)
max |λ∗ − λi | ≤ ρN .
1≤i≤ν
Proof. Since ρN ≤ ρ1 , Lemma 3.4 and Lemma 3.5 yield |d(λ)| > |dN (λ) − d(λ)| for λ − λ = ρN . Hence, by Rouch´e’s Theorem [Con78, §5, Theorem 3.8], d and dN have exactly the same number of zeros counted with their multiplicities in B(λ∗ , ρN ). It is important to point out that the spectral accuracy of the collocation polynomial stated in Lemma 3.2 is preserved in the behavior of the final error (3.14) in Theorem 3.6, if the distributed term in (2.2) can be exactly computed, i.e. εN = 0. On the other hand, we also maintain the spectral accuracy in the error (3.14) if we make the further assumptions that
8
D. BREDA, S. MASET AND R. VERMIGLIO
1. the function M in (1.1) is piecewise C ∞ and there exists C > 0 such that max M (q) (θ) ≤ C q , q = 0, 1, 2, . . . θ∈[−τ,0]
2. a piecewise interpolatory quadrature rule based on N + 1 distinct nodes is used for the distributed term. Finally, we observe that if the quadrature error is independent of N , i.e. εN ≤ TOL (e.g. machine precision), the error bound in theorem 3.6 is still valid with εN replaced by TOL and thus the convergence follows spectral accuracy down to the tolerance TOL. 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 3.7. Let {ΩN (i) }i≥1 be a sequence of meshes on the interval [−τ, 0] such that N (i) → ∞ and εN (i) → 0 for i → ∞. If λ(i) → λ∗ for 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 3.4 and by continuity of d. 4. Implementation issues. The application of the pseudospectral approach presented in this work allows two different alternatives. The first one uses a unique interpolating polynomial for the whole delay interval [−τ, 0]. The second one can be interpreted as a piecewise pseudospectral approach since it makes use of a different N -degree interpolating polynomial for each delay interval [−τl , −τl−1 ], l = 1, . . . , k. The discretized infinitesimal generator results in the m(1 + kN ) × m(1 + kN ) matrix (1) (1) (1) (k−1) (k) (k) a1 · · · aN a0 a1 · · · aN · · · aN (1) (1) (1) d10 d11 · · · d1N . .. . . .. . . . . . (1) (1) (1) d d · · · d NN N0 N1 . . AN = . .. . (k) (k) (k) d10 d11 · · · d1N .. .. . . . . .. . . (k) (k) (k) dN 0 dN 1 · · · dN N (l)
(l)
where, for l = 1, . . . , k, aj with j = 0, . . . , N , and dij with i = 1, . . . , N and j = 0, . . . , N , are the coefficients given by (2.3) and relevant to the interval [−τl , −τl−1 ]. (l) (l−1) Note that a0 = aN holds for l = 2, . . . , k. Despite the increased approximant matrix dimension, this second solution avoids the evaluation of Lagrange coefficients at the discrete delay points when computing the first row since these are always mesh points. Moreover, if we are in the case of applying a quadrature rule to the distributed delay terms (1.2), using the interpolation points as quadrature nodes, the piecewise alternative permits to avoid all Lagrange coefficients evaluation. This
PSEUDOSPECTRAL METHODS FOR CHARACTERISTIC ROOTS
9
advantage clearly increases as the number of delays increases and moreover sparsity can be exploited when k is large. The two alternatives in the pseudospectral approach have been implemented in two different MATLAB routines for the numerical computation of the characteristic roots of system of DDEs using Chebyshev extremal grids for both versions of the general case. The algorithms apply Clenshaw-Curtis quadrature (see [Tre00, §12]) for the distributed delay terms, which is based on Chebyshev extremal nodes. For stability analysis, one can be interested in computing either the first rrightmost roots or the roots λ in the right half-plane (λ) ≥ α for a given real α. In the first case, since the roots with smaller modulus are approximated faster and, in general, they are not the rightmost ones, it is advisable to try with N > r (so it is easier that the first r-rightmost roots are contained in the first r-roots ordered by modulus). In the second case, by using the following bound |λ| ≤
k
Ll e−ατl + τ B max{1, e−ατ } := K
l=0
for roots λ with (λ) ≥ α, where B is a uniform upper bound of M (θ) in [−τ, 0], by (3.13) and (3.14) with εN = 0 one obtains for simple roots N τ (K + ρ1 )e C2 1 ∗ √ |λ − λ| ≤ C3 N N 2 and by assuming C C3 1 and τ (K +ρ1 )e τ Ke one can roughly estimate the required N for a given tolerance. For a more detailed discussion on the above and other implementation issues we refer the interested reader to [Bre02] and [BMV03].
5. Numerical results. We present here a sequence of results on some systems of linear DDEs (1.1) in order to show the features of the pseudospectral approach and to compare it with RK and LMS methods proposed in [Bre02] and [BMV03], respectively. All the results concern the following systems and equations: Example 1 [ER02] y (t) = 2y(t) − ey(t − 1)
(5.1)
Example 2 [ER02] y1 (t) = −0.5y1 (t) − tanh (y1 (t − 1.57)) + tanh (y2 (t − 0.2)) (5.2) y2 (t) = −0.5y2 (t) + 2.34 tanh (y1 (t − 0.2)) − tanh (y2 (t − 1.57)) Example 3 [FSD00] (5.3)
y (t) = L0 y(t) + L1 y(t − 1) +
0 M (θ)y(t + θ)dθ −1
with coefficients matrices −3 1 L0 = , −24.646 −35.430
L1 =
1 0 2.35553 2.00365
,
10
D. BREDA, S. MASET AND R. VERMIGLIO
M (θ) =
1 −θ 1 −θ
,
θ ∈ [−τ, 0]
Example 4 (5.4)
−0.1
y (t) = L0 y(t) + L1 y(t − 1) +
−0.5
M1 y(t + θ)dθ + −0.3
with coefficients matrices −3 1 L0 = , −24.646 −35.430 M1 =
2 2.5 0 −0.5
−1
L1 =
1 0 2.35553 2.00365
,
M2 y(t + θ)dθ
M2 =
−1 0 0 −1
,
Example 5 [ER02] (5.5)
y (t) = L0 y(t) + L1 y(t − τ )
with tridiagonal coefficients m × m matrices −2 1 1 1−α L0 = 2 +α h β
and
1 .. .. . . + RIm .. .. . . 1 1 −2
−2 1 . . 1 1−α 1 .. .. L1 = 2 . . h β .. .. 1 1 −2
where α, β and R are real parameters and h = π/(m + 1) with m ∈ N. Equation (5.1) presents an exact rightmost double root λ = 1 while system (5.2), linearized around the steady state solution (y1∗ , y2∗ ) = (0, 0), exhibits a rightmost root λ 0.3475. Systems (5.3) and (5.4) present, respectively, a rightmost root λ −0.6873 and λ 1.2462. The latter is a modified version of the former: the distributed delay is split in two separate terms. Moreover the matrices M1 and M2 are independent of θ. This modification is introduced in order to exalt the advantages of the piecewise alternative approach with respect to the standard one. The last test equation (5.5) presents a complex-conjugate rightmost root λ −0.06362±1.007i and it arises from the discretization in space by a classical second order central difference scheme of a PDE with time delay in the diffusion term. Figure 5.1 shows the trend of the absolute error on the rightmost root with respect to the number of interpolation points N for the single delay equation (5.1): convergence is confirmed as predicted in Theorem 3.6. Observe that the method works also in the case of root with multiplicity ν > 1, but in this case the maximum tolerance
11
PSEUDOSPECTRAL METHODS FOR CHARACTERISTIC ROOTS 0
10
−2
10
−4
10
−6
10
−8
10
0
1
10
10
Fig. 5.1. Rightmost root absolute error vs number of grid points N for system (5.1).
0
10
0.25 −5
10
0.2 0.15
−10
10
0.1 0.05
−15
10
0
10
1
10
0 0 10
−5
10
−10
10
−15
10
0
10
0.25 −5
10
0.2 0.15
−10
10
0.1 0.05
−15
10
0
10
1
10
0 0 10
−5
10
−10
10
−15
10
Fig. 5.2. Rightmost root absolute error vs number of grid points N (1st column) and computational time (seconds, 2nd column) vs required tolerance T OL for system (5.4) with single (1st row) and piecewise (2nd row) pseudospectral approach.
reached is equal to the machine precision under νth root. The convergence behavior confirms the typical “spectral accuracy” of spectral methods. The different approaches, single and piecewise interpolation, are compared in Figure 5.2 for the computation of the rightmost root of system (5.4). Performances are slightly different in favor of the piecewise method: this is a consequence of avoiding the computation of Lagrange coefficients as explained in §4, which would be more evident in the case of systems with many delays. The approximant matrix dimension is four time larger for the piecewise approach: this disadvantage is recovered (in terms of computational time, see the 2nd column) by using sparse eigenvalue solvers (e.g. the function eigs of MATLAB ). The same analysis is carried out in Figure 5.3 for system (5.3) comparing the piecewise pseudospectral approach with LMS and RK methods of order 5 proposed in [Bre02] and [BMV03], respectively, both showing convergence of order 5 as expected. For what concerns the density of the delay interval discretization (i.e. the number N for the pseudospectral method and N = τ /h, h being the steplength, for the other methods), one can see in the 1st column that the RK approach is more accurate than the pseudospectral one up to a tolerance of about 10−9 , advantage that falls rapidly for smaller values due to the “spectral accuracy” of the pseudospectral
12
D. BREDA, S. MASET AND R. VERMIGLIO
0
10
0.4
−5
10
0.3 0.2
−10
10
0.1 −15
10
0
10
1
10
0 0 10
2
10
−5
10
−10
10
−15
10
0
10
0.4
−5
10
0.3 0.2
−10
10
0.1 −15
10
0
10
1
10
0 0 10
2
10
−5
10
−10
10
−15
10
0
10
0.4
−5
10
0.3 0.2
−10
10
0.1 −15
10
0
10
1
10
2
10
0 0 10
−5
10
−10
10
−15
10
Fig. 5.3. Rightmost root absolute error vs number of grid points N (1st column) and computational time (seconds, 2nd column) vs required tolerance T OL for system (5.3) with RK (1st row), LMS (2nd row) and pseudospectral (3rd row) approach.
technique. This behavior holds in general for “linearly” convergent approaches (LMS and RK) compared to “superlinearly” convergent ones (pseudospectral), the “crossing” tolerance value depending on the particular system of DDEs analyzed by the constants involved in the error (here, “linear” and “superlinear” convergence refer to a double-logarithmic plot context). Moreover we compare RK and pseudospectral approaches computing the first ten rightmost roots of system (5.2), listed in Table 5.1 and shown in Figure 5.4. The 1st column of Figure 5.5 shows the error trend for each root (simple-real or complexconjugate) in ascending order with respect to their modulus. This last plays a fundamental role in the sense that approximated roots of smaller modulus converge sooner: for the pseudospectral approach this phenomenon clearly results from Theorem 3.6 and an analogous result can be deducted from the convergence analysis carried out in [BMV03]. The 2nd column evidences 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. Table 5.1 First ten rightmost roots value and modulus for system (5.2). characteristic root 0.34748172572643 -0.08116698520245 -0.43412304132027±i1.62752128215701 -0.82061576063772±i5.11180446782392 -1.20975894908637±i4.82695048752634 -1.24622687378869±i8.90021501181034
modulus 0.34748172572643 0.08116698520245 1.68442522507803 5.17725362947110 4.97624032015389 8.98704115253777
13
PSEUDOSPECTRAL METHODS FOR CHARACTERISTIC ROOTS 10
8
6
4
ℑ(λ)
2
0
−2
−4
−6
−8
−10 −1.4
−1.2
−1
−0.8
−0.6
−0.4
−0.2
0
0.2
0.4
ℜ(λ)
Fig. 5.4. First ten rightmost roots of system (5.2).
0
10
1.2 1
−5
10
0.8 0.6 −10
10
0.4 0.2
−15
10
0
10
1
10
0 0 10
2
10
−5
10
−10
10
−15
10
0
10
1.2 1
−5
10
0.8 0.6 −10
10
0.4 0.2
−15
10
0
10
1
10
0 0 10
2
10
−5
10
−10
10
−15
10
Fig. 5.5. First ten rightmost roots absolute error vs number of grid points N (1st column) and computational time (seconds, 2nd column) vs required tolerance on the maximum absolute error T OL (seconds, 2nd column) for system (5.2) with RK (1st row) and pseudospectral (2nd row) approach.
Finally, to test the performance of our approach on time delay systems of higher order, we carry on the convergence analysis on the rightmost characteristic root of (5.5) with m = 20 and parameters α = 0.1, β = 2 and R = 0.08. Results are given in Figure 5.6. Note that in general space discretization of PDEs with delay in time leads to DDE systems which possess already a quite large dimension m. This means that the choice of N for the approximation of the characteristic roots turns out to be crucial in saving computational time. Let us underline that in the second case we adopt a sparse eigenvalue solver to exploit the fact that the entries of AN are (apart from the first block row) m × m diagonal blocks. 0
10
0.5
0.45
−2
10
0.4 −4
10
0.35 −6
10
0.3
−8
10
0.25
0.2
−10
10
0.15 −12
10
0.1 −14
10
0.05
−16
10
0
10
1
10
0 0 10
−5
10
−10
10
−15
10
Fig. 5.6. Rightmost root absolute error vs number of grid points N (1st column) and computational time (seconds, 2nd column) vs required tolerance T OL for system (5.5) with pseudospectral approach.
14
D. BREDA, S. MASET AND R. VERMIGLIO
To complete this test we use N = 100 to allow the computation of r = m(1+N ) = 2020 roots, which is the dimension of the discretized infinitesimal generator. The eigenvalues are shown in Figure 5.7 (left), while a zoom in Figure 5.7 (right) gives an idea of the particular conformation of the spectrum and of how clustered the roots are. 600
35
30
400
25 200
ℑ(λ)
ℑ(λ)
20 0
15 −200
10
−400
5
−600 −600
−500
−400
−300
−200
−100
0
0 −0.1
−0.095
ℜ(λ)
−0.09
−0.085
−0.08 ℜ(λ)
−0.075
−0.07
−0.065
−0.06
Fig. 5.7. Characteristic roots of system (5.5).
REFERENCES [Bre03] D. Breda, Solution operator approximation for delay differential equation characteristic roots computation via Runge-Kutta methods, special issue of Appl. Numer. Math. for the Third Conference on Volterra and Delay Equations dedicated to Professor Alan Feldstein’s 70th Birthday, to appear. [Bre02] D. Breda, The infinitesimal generator approach for the computation of characteristic roots for delay differential equations using BDF methods, Research Report UDMI RR17/2002 (2002), Dipartimento di Matematica e Informatica, Universit` a degli Studi di Udine, Italy. [BMV03] D. Breda, S. Maset and R. Vermiglio, Computing the characteristic roots for delay differential equations, IMA J. Numer. Anal., 24 (1) (2004), pp. 1-19. [Con78] J.B. Conway, Functions of One Complex Variable, Springer Verlag, GTM series n. 11, 2nd edition, New York, U.S.A., 1978. [DGVLW95] O. Diekmann, S.A. van Gils, S.M. Verduyn Lunel and H.O. Walther, Delay Equations - Functional, Complex and Nonlinear Analysis, Springer Verlag, AMS series n. 110, New York, U.S.A., 1995. [ELR02] K. Engelborghs, T. Luzyanina and D. Roose, Numerical bifurcation analysis of delay differential equations using DDE-BIFTOOL, ACM Trans. Math. Softw., 28 (1) (2002), pp. 1-21. [ELS01] K. Engelborghs, T. Luzyanina, G. Samaey, DDE-BIFTOOL v. 2.00: a Matlab package for bifurcation analysis of delay differential equations, Report TW330 (2001), Department of Computer Science, K. U. Leuven, Belgium. [ER02] K. Engelborghs and D. Roose, On Stability of LMS methods and characteristic roots of delay differential equations, SIAM J. Numer. Anal., 40 (2) (2002), pp. 629-650. [ER99] K. Engelborghs and 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), pp. 271-289. [FSD00] A. Fattouh, O. Sename and J.M. Dion, H∞ controller and observer design for linear systems with point and distributed delays, proceedings of 2nd IFAC workshop on Linear Time Delay Systems, 11-13 September 2000, Ancona, Italy. [HVL93] J.K. Hale and S.M. Verduyn Lunel, Introduction to Functional Differential Equations, Springer-Verlag, AMS series n. 99, New York, U.S.A., 1993. [LER03] T. Luzyanina, K. Engelborghs and D. Roose, Computing stability of differential equations with bounded distributed delays, Numer. Algorithms, 34 (1) (2003), pp. 41-66. [Nic00] S.I. Niculescu, Delay Effects on Stability: A Robust Control Approach, Monograph Springer, LNCIS n.269, New York, U.S.A., 2000. [Ric03] J.P. Richard, Time-delay systems: an overview of some recent advances and open problems, Automatica, 39 (2003), pp. 1667-1694.
PSEUDOSPECTRAL METHODS FOR CHARACTERISTIC ROOTS
15
[Tre00] L.N. Trefethen,Spectral methods in MATLAB, SIAM, Software - Environment - Tools series, Philadelphia, U.S.A., 2000.
1. Introduction. The interest for delay differential equations (DDEs) is still growing in many application fields, since the introduction of the delay in the model enriches its dynamics and allows a closer and more accurate description of the real-life phenomena. In this paper we consider a system of m-dimensional linear (or linearized) DDEs with multiple discrete and distributed delays:
(1.1)
y (t) = L0 y(t) +
k
0 Ll y(t − τl ) +
l=1
M (θ)y(t + θ)dθ,
t ≥ 0,
−τ
where L0 , L1 , . . . , Lk ∈ Cm×m , 0 = τ0 < τ1 < . . . < τk = τ and M : [−τ, 0] → Cm×m is a piecewise smooth function. Indeed in the implementation of the method we write the distributed term as 0 (1.2) −τ
−τ l−1 k M (θ)y(t + θ)dθ = Ml (θ)y(t + θ)dθ,
t ≥ 0,
l=1 −τ l
where every function Ml , l = 1, . . . , k, is smooth. Delay systems such as (1.1) are particularly important in control theory, where the stability effects of delays is a crucial problem (see the recent survey [Ric03]). In the exhaustive monograph [Nic00] many graphical and analytical tests to this aim are reported, but they work only for particular classes of system (1.1). This motivates our interest in the numerical computation of the characteristic roots of (1.1), i.e. the (infinitely many) roots of (1.3)
det(∆(λ)) = 0,
† Dipartimento di Matematica e Informatica, Universit` a degli Studi di Udine, Via delle Scienze 208, I-33100 Udine, Italy. ‡ Dipartimento di Matematica e Informatica, Universit` a degli Studi di Trieste, Via Valerio 12, I-34127 Trieste, Italy, ([email protected]). § ([email protected]). ¶ ([email protected]).
1
2
D. BREDA, S. MASET AND R. VERMIGLIO
where ∆(λ) = λI −
(1.4)
k
Ll e
−λτl
0 −
l=0
M (θ)eλθ dθ, λ ∈ C, −τ
since it is well known that the zero solution of (1.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 to (1.1) or the infinitesimal generator of the solution operator semigroup. We briefly recall that the solution operator T (t), t ≥ 0, associated to (1.1) is defined by T (t)ϕ = yt ,
ϕ ∈ X,
where X = C ([−τ, 0] , Cm ) endowed with the maximum norm, yt is the function yt (θ) = y(t + θ),
θ ∈ [−τ, 0],
and y is the solution of (1.1) with initial data ϕ ∈ X. The family {T (t)}t≥0 is a C0 -semigroup with infinitesimal generator A : D(A) ⊆ X → X given by Aϕ = ϕ , ϕ ∈ D (A) ,
(1.5) with domain (1.6) D (A) =
ϕ ∈ X | ϕ ∈ X and ϕ (0) =
k l=0
0 Ll ϕ(−τl ) +
M (θ)ϕ(θ)dθ −τ
.
So (1.1) can be restated as the following abstract Cauchy problem ([DGVLW95]): dyt t>0 dt = Ayt , . y0 = ϕ, The two following important results ([DGVLW95], [HVL93]) 1. det(∆(λ)) = 0 ⇔ λ = 1t ln µ, µ ∈ σ(T (t)) \ {0}; 2. det(∆(λ)) = 0 ⇔ λ ∈ σ(A); where σ(·) denotes the spectrum, suggest the idea to turn the characteristic roots approximation problem into a corresponding eigenvalue problem for suitable matrix discretization of either T (t) (i.e. solution operator approach) or A (i.e. infinitesimal generator approach). Engelborghs and Roose propose in [ER02] and [ER99] the solution operator approach via LMS time integration for system (1.1) without distributed delay term. Their method computes approximations to the roots from a large, standard and sparse eigenvalue problem and it is implemented in the MATLAB package DDE-BIFTOOL for DDEs bifurcation analysis ([ELR02] and [ELS01]). The distributed delay case is considered in [LER03] by using LMS methods and in [Bre03] by using RK methods. The complete development of the infinitesimal generator approach first appears in [BMV03] and [Bre02] where a matrix approximation to A is obtained discretizing the derivative in (1.5) by RK or LMS method, respectively, ending up with a large and sparse eigenvalue problem as in [ER99].
PSEUDOSPECTRAL METHODS FOR CHARACTERISTIC ROOTS
3
This work continues the investigation along this line and proposes the infinitesimal generator approach via pseudospectral differencing methods (see [Tre00] and the bibliography therein). This technique involves the exact differentiation of interpolants at selected sets of points. The resulting differentiation matrix is non-sparse, but we can take advantage of the well-known “spectral accuracy” to obtain very accurate approximation with small matrix dimension. This behavior represents in fact, for sufficiently small tolerance, the outstanding advantage of this method compared to the previously cited discretization schemes. The paper is structured as follows. Section 2 describes the pseudospectral approach. Section 3 is about convergence analysis of the method. In §4 we discuss about implementation issues. In §5 we perform numerical experiments on test delay systems of type (1.1) to illustrate the properties of the new method and to critically compare it with LMS and RK discretization techniques. 2. The pseudospectral approach. To simplify the notation, we rewrite the system of DDEs (1.1) as follows: y (t) = f (yt ) where (2.1)
f (ϕ) =
k l=0
0 Ll ϕ(−τl ) +
M (θ)ϕ(θ)dθ,
ϕ ∈ X.
−τ
Given a positive integer N , let us consider the mesh ΩN of N + 1 distinct points in [−τ, 0] ΩN = {θN,i , i = 0, 1, . . . , N } with 0 = θN,0 > θN,1 > . . . > θN,N ≥ −τ . We replace the continuous space X by the space XN of the discrete functions defined on the mesh ΩN , i.e. any ϕ ∈ X is discretized into the block-vector x ∈ XN of components xi = ϕ(θN,i ) ∈ Cm , i = 0, 1, . . . , N . Let LN x, x ∈ XN , be the unique Cm -valued interpolating polynomial of degree ≤ N with (LN x)(θN,i ) = xi , i = 0, 1, . . . , N . Thus we approximate the infinitesimal generator A by the matrix AN : XN → XN , called spectral differentiation matrix, defined as follows: (AN x)0 = fN (LN x) (2.2) (AN x)i = (LN x) (θN,i ), i = 1, 2, . . . , N where fN is an approximation of f in which the distributed delay integral term in (2.1) is substituted by a suitable interpolatory quadrature rule or fN is equal to f in the case of simple matrix function M for which the integral can be exactly computed. By using the Lagrange representation of LN x we obtain a0 a1 · · · aN d10 d11 · · · d1N m(N +1)×m(N +1) AN = . .. . . .. ∈ C .. . . . dN 0 dN 1 · · · dN N with (2.3)
aj = fN (lN,j (·)Im ) , j = 0, 1, . . . , N (θN,i )Im , i = 1, . . . , N, j = 0, 1, . . . , N dij = lN,j
4
D. BREDA, S. MASET AND R. VERMIGLIO
where lN,j , j = 0, 1, . . . , N , are the Lagrange coefficients relevant to the nodes in m ΩN and, by definition, fN (lN,j (·)Im ) = fN lN,j (·)e(i) i=1 with e(i) ’s the canonical vectors in Rm . Explicit expressions of the dij ’s for particular choices of ΩN (e.g. Chebyshev extremal points) can be found in [Tre00, §6]. So finally the original problem of the numerical computation of the characteristic roots of (1.1) can be turned into the eigenvalue problem for the matrix AN , i.e. eigenvalues of AN directly approximate the characteristic roots. How accurate are these approximations? In the next section we develop the convergence analysis. 3. Convergence analysis. The convergence analysis of the approximated roots to the exact ones is based on a comparison between the continuous (i.e. exact) characteristic equation (1.3) and its discrete (i.e. approximated) counterpart. Let us thus characterize the continuous and the discrete eigenvalue problems for A and AN , respectively. Let λ ∈ C and ϕ ∈ D (A) be such that Aϕ = λϕ, that means
ϕ (θ) = λϕ(θ), ϕ (0) = f (ϕ) .
θ ∈ [−τ, 0]
Hence λ is eigenvalue of A if and only if there exists u ∈ Cm , u = 0, such that λu = f eλ· u . ˆ By introducing the linear operator A(λ) : Cm → Cm given by ˆ A(λ)u = f eλ· u , u ∈ Cm (3.1) we get that λ is eigenvalue of A if and only if ˆ det λIm − A(λ) (3.2) =0 which is the continuous characteristic equation (1.3). In the same fashion, let x ∈ XN be such that AN x = λx, thus by (2.2) we get fN (LN x) = λx0 and (3.3)
(LN x) (θN,i ) = λxi ,
i = 1, 2, . . . , N.
By defining pN (·; λ, u), u ∈ Cm , as the collocation polynomial relevant to the nodes θN,1 , θN,2 , . . . , θN,N for the initial value problem y (θ) = λy(θ), θ ∈ [−τ, 0] (3.4) y(0) = u,
PSEUDOSPECTRAL METHODS FOR CHARACTERISTIC ROOTS
5
whose solution is the exponential y(·; λ, u) = eλ· u, one can easily derive that (3.3) is equivalent to (LN x)(·) = pN (·; λ, x0 ). Consequently λ is eigenvalue of AN if and only if there exists u ∈ Cm , u = 0, such that λu = fN (pN (·; λ, u)) . Finally, by introducing the linear operator AˆN (λ) : Cm → Cm given by (3.5)
AˆN (λ)u = fN (pN (·; λ, u)) ,
u ∈ Cm ,
we get that λ is eigenvalue of AN if and only if (3.6) det λIm − AˆN (λ) = 0 which is the discrete characteristic equation. Assumption 3.1. We assume to use the Chebyshev extremal nodes on [−τ, 0], i.e. π τ θN,i = (3.7) cos (i ) − 1 , i = 0, . . . , N. 2 N Let us denote by B(λ, ρ) the closed ball in C of center λ and radius ρ. In next lemma we give an upper bound of the error of the collocation polynomial for (3.4). Lemma 3.2. Let λ∗ ∈ C and ρ0 > 0. Under Assumption 3.1 there exists N0 ∈ N such that, for N ≥ N0 , for λ ∈ B(λ∗ , ρ0 ) and for u ∈ Cm , the collocation polynomial pN (·; λ, u) of (3.4) based on (3.7) exists and it is unique. Moreover it holds C0 max |pN (θ; λ, u) − y(θ; λ, u)| ≤ √ θ∈[−τ,0] N
C1 N
N |u|
where C0 = C0 (λ∗ , ρ0 ) and C1 = C1 (λ∗ , ρ0 ) = τ (|λ∗ | + ρ0 )e are constants independent of N . Proof. By introducing the (linear) integral Volterra operator Kλ : X → X as θ (Kλ ϕ)(θ) = λ
ϕ(s)ds,
ϕ ∈ X,
λ ∈ C,
θ ∈ [−τ, 0],
0
and the (linear and bounded) Lagrange interpolation operator LN −1 : X → X relevant to the nodes (3.7), we can rewrite (3.4) and its collocation equations, respectively, as the equations on the space X (3.8)
y = u + Kλ y
and (3.9)
pN = u + Kλ LN −1 pN
where y = y(·; λ, u), pN = pN (·; λ, u) and u denotes the function of X of constant value u. By subtracting (3.8) from (3.9) and setting eN = pN − y and rN = (LN −1 − I)y, we get the following equation for eN (3.10)
eN = Kλ LN −1 eN + Kλ rN .
6
D. BREDA, S. MASET AND R. VERMIGLIO
It can be directly shown that the solutions of (3.10) are the functions Kλ eˆN , where eˆN is a solution of eˆN = LN −1 Kλ eˆN + rN . Since I − Kλ is invertible and, by the known estimate for the Lagrange interpolation on Chebyshev nodes, LN −1 Kλ − Kλ → 0 as N → ∞, it follows that, for sufficiently large N (and independently of λ), I − LN −1 Kλ is invertible and (I − LN −1 Kλ )−1 ≤ 2 (I − Kλ )−1 . Hence eN = Kλ (I − LN −1 Kλ )−1 rN is the unique solution of (3.10) and as a consequence eN ≤ 2 Kλ (I − Kλ )−1 rN . Since rN is the error of the Lagrange interpolation of the exponential y(·; λ, u) = eλ· u relevant to the nodes (3.7), we conclude that rN ≤ max{1, e−(λ)τ }
(τ |λ|)N |u|. N!
The final assertion follows by application of the Stirling formula √ N ! ≥ 2πN (N/e)N . Now let us set the following notation relevant to the continuous and discrete characteristic equations (3.2) and (3.6) respectively: ˆ (3.11) d(λ) = det λIm − A(λ) dN (λ) = det λIm − AˆN (λ)
(3.12)
ˆ with A(λ) and AˆN (λ) defined in (3.1) and (3.5), respectively. The following assumption and lemmas are the next necessary steps to state the final theorem of convergence. Assumption 3.3. We assume that the function fN in (2.2) satisfies sup fN < +∞.
N ∈N
Lemma 3.4. Let λ∗ ∈ C and ρ0 > 0. Under Assumptions 3.1 and 3.3, there exists N0 ∈ N such that, for N ≥ N0 and λ ∈ B(λ∗ , ρ0 ), we have N C1 1 |d(λ) − dN (λ)| ≤ C2 εN + √ N N where εN =
sup λ∈B(λ∗ ,ρ0 )
|f (y(·; λ, u)) − fN (y(·; λ, u)) | , |u|
PSEUDOSPECTRAL METHODS FOR CHARACTERISTIC ROOTS
7
C1 = C1 (λ∗ , ρ0 ) as defined in Lemma 3.2 and C2 = C2 (λ∗ , ρ0 ) are suitable constants independent of N . Proof. By Lemma 3.2 : ˆ |A(λ)u − AˆN (λ)u| ≤ |f (y(·; λ, u)) − fN (y(·; λ, u)) | + + |fN (y(·; λ, u)) − fN (pN (·; λ, u)) | ≤ N C0 C1 |u| ≤ εN + sup fN √ N N N ∈N for N ≥ N0 and λ ∈ B(λ∗ , ρ0 ). Now, since the derivative of det(A) with respect to A ∈ Cm×m is continuous, the assertion easily follows. Lemma 3.5. If the function d(λ) has a zero λ∗ ∈ C with multiplicity ν, then there exists ρ1 = ρ1 (λ∗ ) > 0 such that |d(λ)| > C3 |λ − λ∗ |
ν
for λ ∈ B(λ∗ , ρ1 ) \ {λ∗ } where C3 = C3 (λ∗ ). Proof. Expanding d(λ) around λ∗ and taking into account the multiplicity ν we obtain d(λ) =
d(ν) (λ∗ ) ν+1 (λ − λ∗ )ν + O |λ − λ∗ | ν!
with d(ν) (λ∗ ) = 0. The result follows. We are now able to state the convergence result. Theorem 3.6 (convergence). Let λ∗ ∈ C be a zero of the function d(λ) with multiplicity ν, let ρ1 as in Lemma 3.5 and let (3.13)
ρN =
C2 C3
1/ν N 1/ν C1 1 εN + √ N N
with C1 , C2 as in Lemma 3.4 with ρ0 = ρ1 and C3 as in Lemma 3.5. Under Assumptions 3.1 and 3.3, if ρN ≤ ρ1 then the function dN (λ) has ν zeros λi , i = 1, . . . , ν, (taking into account multiplicities) such that (3.14)
max |λ∗ − λi | ≤ ρN .
1≤i≤ν
Proof. Since ρN ≤ ρ1 , Lemma 3.4 and Lemma 3.5 yield |d(λ)| > |dN (λ) − d(λ)| for λ − λ = ρN . Hence, by Rouch´e’s Theorem [Con78, §5, Theorem 3.8], d and dN have exactly the same number of zeros counted with their multiplicities in B(λ∗ , ρN ). It is important to point out that the spectral accuracy of the collocation polynomial stated in Lemma 3.2 is preserved in the behavior of the final error (3.14) in Theorem 3.6, if the distributed term in (2.2) can be exactly computed, i.e. εN = 0. On the other hand, we also maintain the spectral accuracy in the error (3.14) if we make the further assumptions that
8
D. BREDA, S. MASET AND R. VERMIGLIO
1. the function M in (1.1) is piecewise C ∞ and there exists C > 0 such that max M (q) (θ) ≤ C q , q = 0, 1, 2, . . . θ∈[−τ,0]
2. a piecewise interpolatory quadrature rule based on N + 1 distinct nodes is used for the distributed term. Finally, we observe that if the quadrature error is independent of N , i.e. εN ≤ TOL (e.g. machine precision), the error bound in theorem 3.6 is still valid with εN replaced by TOL and thus the convergence follows spectral accuracy down to the tolerance TOL. 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 3.7. Let {ΩN (i) }i≥1 be a sequence of meshes on the interval [−τ, 0] such that N (i) → ∞ and εN (i) → 0 for i → ∞. If λ(i) → λ∗ for 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 3.4 and by continuity of d. 4. Implementation issues. The application of the pseudospectral approach presented in this work allows two different alternatives. The first one uses a unique interpolating polynomial for the whole delay interval [−τ, 0]. The second one can be interpreted as a piecewise pseudospectral approach since it makes use of a different N -degree interpolating polynomial for each delay interval [−τl , −τl−1 ], l = 1, . . . , k. The discretized infinitesimal generator results in the m(1 + kN ) × m(1 + kN ) matrix (1) (1) (1) (k−1) (k) (k) a1 · · · aN a0 a1 · · · aN · · · aN (1) (1) (1) d10 d11 · · · d1N . .. . . .. . . . . . (1) (1) (1) d d · · · d NN N0 N1 . . AN = . .. . (k) (k) (k) d10 d11 · · · d1N .. .. . . . . .. . . (k) (k) (k) dN 0 dN 1 · · · dN N (l)
(l)
where, for l = 1, . . . , k, aj with j = 0, . . . , N , and dij with i = 1, . . . , N and j = 0, . . . , N , are the coefficients given by (2.3) and relevant to the interval [−τl , −τl−1 ]. (l) (l−1) Note that a0 = aN holds for l = 2, . . . , k. Despite the increased approximant matrix dimension, this second solution avoids the evaluation of Lagrange coefficients at the discrete delay points when computing the first row since these are always mesh points. Moreover, if we are in the case of applying a quadrature rule to the distributed delay terms (1.2), using the interpolation points as quadrature nodes, the piecewise alternative permits to avoid all Lagrange coefficients evaluation. This
PSEUDOSPECTRAL METHODS FOR CHARACTERISTIC ROOTS
9
advantage clearly increases as the number of delays increases and moreover sparsity can be exploited when k is large. The two alternatives in the pseudospectral approach have been implemented in two different MATLAB routines for the numerical computation of the characteristic roots of system of DDEs using Chebyshev extremal grids for both versions of the general case. The algorithms apply Clenshaw-Curtis quadrature (see [Tre00, §12]) for the distributed delay terms, which is based on Chebyshev extremal nodes. For stability analysis, one can be interested in computing either the first rrightmost roots or the roots λ in the right half-plane (λ) ≥ α for a given real α. In the first case, since the roots with smaller modulus are approximated faster and, in general, they are not the rightmost ones, it is advisable to try with N > r (so it is easier that the first r-rightmost roots are contained in the first r-roots ordered by modulus). In the second case, by using the following bound |λ| ≤
k
Ll e−ατl + τ B max{1, e−ατ } := K
l=0
for roots λ with (λ) ≥ α, where B is a uniform upper bound of M (θ) in [−τ, 0], by (3.13) and (3.14) with εN = 0 one obtains for simple roots N τ (K + ρ1 )e C2 1 ∗ √ |λ − λ| ≤ C3 N N 2 and by assuming C C3 1 and τ (K +ρ1 )e τ Ke one can roughly estimate the required N for a given tolerance. For a more detailed discussion on the above and other implementation issues we refer the interested reader to [Bre02] and [BMV03].
5. Numerical results. We present here a sequence of results on some systems of linear DDEs (1.1) in order to show the features of the pseudospectral approach and to compare it with RK and LMS methods proposed in [Bre02] and [BMV03], respectively. All the results concern the following systems and equations: Example 1 [ER02] y (t) = 2y(t) − ey(t − 1)
(5.1)
Example 2 [ER02] y1 (t) = −0.5y1 (t) − tanh (y1 (t − 1.57)) + tanh (y2 (t − 0.2)) (5.2) y2 (t) = −0.5y2 (t) + 2.34 tanh (y1 (t − 0.2)) − tanh (y2 (t − 1.57)) Example 3 [FSD00] (5.3)
y (t) = L0 y(t) + L1 y(t − 1) +
0 M (θ)y(t + θ)dθ −1
with coefficients matrices −3 1 L0 = , −24.646 −35.430
L1 =
1 0 2.35553 2.00365
,
10
D. BREDA, S. MASET AND R. VERMIGLIO
M (θ) =
1 −θ 1 −θ
,
θ ∈ [−τ, 0]
Example 4 (5.4)
−0.1
y (t) = L0 y(t) + L1 y(t − 1) +
−0.5
M1 y(t + θ)dθ + −0.3
with coefficients matrices −3 1 L0 = , −24.646 −35.430 M1 =
2 2.5 0 −0.5
−1
L1 =
1 0 2.35553 2.00365
,
M2 y(t + θ)dθ
M2 =
−1 0 0 −1
,
Example 5 [ER02] (5.5)
y (t) = L0 y(t) + L1 y(t − τ )
with tridiagonal coefficients m × m matrices −2 1 1 1−α L0 = 2 +α h β
and
1 .. .. . . + RIm .. .. . . 1 1 −2
−2 1 . . 1 1−α 1 .. .. L1 = 2 . . h β .. .. 1 1 −2
where α, β and R are real parameters and h = π/(m + 1) with m ∈ N. Equation (5.1) presents an exact rightmost double root λ = 1 while system (5.2), linearized around the steady state solution (y1∗ , y2∗ ) = (0, 0), exhibits a rightmost root λ 0.3475. Systems (5.3) and (5.4) present, respectively, a rightmost root λ −0.6873 and λ 1.2462. The latter is a modified version of the former: the distributed delay is split in two separate terms. Moreover the matrices M1 and M2 are independent of θ. This modification is introduced in order to exalt the advantages of the piecewise alternative approach with respect to the standard one. The last test equation (5.5) presents a complex-conjugate rightmost root λ −0.06362±1.007i and it arises from the discretization in space by a classical second order central difference scheme of a PDE with time delay in the diffusion term. Figure 5.1 shows the trend of the absolute error on the rightmost root with respect to the number of interpolation points N for the single delay equation (5.1): convergence is confirmed as predicted in Theorem 3.6. Observe that the method works also in the case of root with multiplicity ν > 1, but in this case the maximum tolerance
11
PSEUDOSPECTRAL METHODS FOR CHARACTERISTIC ROOTS 0
10
−2
10
−4
10
−6
10
−8
10
0
1
10
10
Fig. 5.1. Rightmost root absolute error vs number of grid points N for system (5.1).
0
10
0.25 −5
10
0.2 0.15
−10
10
0.1 0.05
−15
10
0
10
1
10
0 0 10
−5
10
−10
10
−15
10
0
10
0.25 −5
10
0.2 0.15
−10
10
0.1 0.05
−15
10
0
10
1
10
0 0 10
−5
10
−10
10
−15
10
Fig. 5.2. Rightmost root absolute error vs number of grid points N (1st column) and computational time (seconds, 2nd column) vs required tolerance T OL for system (5.4) with single (1st row) and piecewise (2nd row) pseudospectral approach.
reached is equal to the machine precision under νth root. The convergence behavior confirms the typical “spectral accuracy” of spectral methods. The different approaches, single and piecewise interpolation, are compared in Figure 5.2 for the computation of the rightmost root of system (5.4). Performances are slightly different in favor of the piecewise method: this is a consequence of avoiding the computation of Lagrange coefficients as explained in §4, which would be more evident in the case of systems with many delays. The approximant matrix dimension is four time larger for the piecewise approach: this disadvantage is recovered (in terms of computational time, see the 2nd column) by using sparse eigenvalue solvers (e.g. the function eigs of MATLAB ). The same analysis is carried out in Figure 5.3 for system (5.3) comparing the piecewise pseudospectral approach with LMS and RK methods of order 5 proposed in [Bre02] and [BMV03], respectively, both showing convergence of order 5 as expected. For what concerns the density of the delay interval discretization (i.e. the number N for the pseudospectral method and N = τ /h, h being the steplength, for the other methods), one can see in the 1st column that the RK approach is more accurate than the pseudospectral one up to a tolerance of about 10−9 , advantage that falls rapidly for smaller values due to the “spectral accuracy” of the pseudospectral
12
D. BREDA, S. MASET AND R. VERMIGLIO
0
10
0.4
−5
10
0.3 0.2
−10
10
0.1 −15
10
0
10
1
10
0 0 10
2
10
−5
10
−10
10
−15
10
0
10
0.4
−5
10
0.3 0.2
−10
10
0.1 −15
10
0
10
1
10
0 0 10
2
10
−5
10
−10
10
−15
10
0
10
0.4
−5
10
0.3 0.2
−10
10
0.1 −15
10
0
10
1
10
2
10
0 0 10
−5
10
−10
10
−15
10
Fig. 5.3. Rightmost root absolute error vs number of grid points N (1st column) and computational time (seconds, 2nd column) vs required tolerance T OL for system (5.3) with RK (1st row), LMS (2nd row) and pseudospectral (3rd row) approach.
technique. This behavior holds in general for “linearly” convergent approaches (LMS and RK) compared to “superlinearly” convergent ones (pseudospectral), the “crossing” tolerance value depending on the particular system of DDEs analyzed by the constants involved in the error (here, “linear” and “superlinear” convergence refer to a double-logarithmic plot context). Moreover we compare RK and pseudospectral approaches computing the first ten rightmost roots of system (5.2), listed in Table 5.1 and shown in Figure 5.4. The 1st column of Figure 5.5 shows the error trend for each root (simple-real or complexconjugate) in ascending order with respect to their modulus. This last plays a fundamental role in the sense that approximated roots of smaller modulus converge sooner: for the pseudospectral approach this phenomenon clearly results from Theorem 3.6 and an analogous result can be deducted from the convergence analysis carried out in [BMV03]. The 2nd column evidences 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. Table 5.1 First ten rightmost roots value and modulus for system (5.2). characteristic root 0.34748172572643 -0.08116698520245 -0.43412304132027±i1.62752128215701 -0.82061576063772±i5.11180446782392 -1.20975894908637±i4.82695048752634 -1.24622687378869±i8.90021501181034
modulus 0.34748172572643 0.08116698520245 1.68442522507803 5.17725362947110 4.97624032015389 8.98704115253777
13
PSEUDOSPECTRAL METHODS FOR CHARACTERISTIC ROOTS 10
8
6
4
ℑ(λ)
2
0
−2
−4
−6
−8
−10 −1.4
−1.2
−1
−0.8
−0.6
−0.4
−0.2
0
0.2
0.4
ℜ(λ)
Fig. 5.4. First ten rightmost roots of system (5.2).
0
10
1.2 1
−5
10
0.8 0.6 −10
10
0.4 0.2
−15
10
0
10
1
10
0 0 10
2
10
−5
10
−10
10
−15
10
0
10
1.2 1
−5
10
0.8 0.6 −10
10
0.4 0.2
−15
10
0
10
1
10
0 0 10
2
10
−5
10
−10
10
−15
10
Fig. 5.5. First ten rightmost roots absolute error vs number of grid points N (1st column) and computational time (seconds, 2nd column) vs required tolerance on the maximum absolute error T OL (seconds, 2nd column) for system (5.2) with RK (1st row) and pseudospectral (2nd row) approach.
Finally, to test the performance of our approach on time delay systems of higher order, we carry on the convergence analysis on the rightmost characteristic root of (5.5) with m = 20 and parameters α = 0.1, β = 2 and R = 0.08. Results are given in Figure 5.6. Note that in general space discretization of PDEs with delay in time leads to DDE systems which possess already a quite large dimension m. This means that the choice of N for the approximation of the characteristic roots turns out to be crucial in saving computational time. Let us underline that in the second case we adopt a sparse eigenvalue solver to exploit the fact that the entries of AN are (apart from the first block row) m × m diagonal blocks. 0
10
0.5
0.45
−2
10
0.4 −4
10
0.35 −6
10
0.3
−8
10
0.25
0.2
−10
10
0.15 −12
10
0.1 −14
10
0.05
−16
10
0
10
1
10
0 0 10
−5
10
−10
10
−15
10
Fig. 5.6. Rightmost root absolute error vs number of grid points N (1st column) and computational time (seconds, 2nd column) vs required tolerance T OL for system (5.5) with pseudospectral approach.
14
D. BREDA, S. MASET AND R. VERMIGLIO
To complete this test we use N = 100 to allow the computation of r = m(1+N ) = 2020 roots, which is the dimension of the discretized infinitesimal generator. The eigenvalues are shown in Figure 5.7 (left), while a zoom in Figure 5.7 (right) gives an idea of the particular conformation of the spectrum and of how clustered the roots are. 600
35
30
400
25 200
ℑ(λ)
ℑ(λ)
20 0
15 −200
10
−400
5
−600 −600
−500
−400
−300
−200
−100
0
0 −0.1
−0.095
ℜ(λ)
−0.09
−0.085
−0.08 ℜ(λ)
−0.075
−0.07
−0.065
−0.06
Fig. 5.7. Characteristic roots of system (5.5).
REFERENCES [Bre03] D. Breda, Solution operator approximation for delay differential equation characteristic roots computation via Runge-Kutta methods, special issue of Appl. Numer. Math. for the Third Conference on Volterra and Delay Equations dedicated to Professor Alan Feldstein’s 70th Birthday, to appear. [Bre02] D. Breda, The infinitesimal generator approach for the computation of characteristic roots for delay differential equations using BDF methods, Research Report UDMI RR17/2002 (2002), Dipartimento di Matematica e Informatica, Universit` a degli Studi di Udine, Italy. [BMV03] D. Breda, S. Maset and R. Vermiglio, Computing the characteristic roots for delay differential equations, IMA J. Numer. Anal., 24 (1) (2004), pp. 1-19. [Con78] J.B. Conway, Functions of One Complex Variable, Springer Verlag, GTM series n. 11, 2nd edition, New York, U.S.A., 1978. [DGVLW95] O. Diekmann, S.A. van Gils, S.M. Verduyn Lunel and H.O. Walther, Delay Equations - Functional, Complex and Nonlinear Analysis, Springer Verlag, AMS series n. 110, New York, U.S.A., 1995. [ELR02] K. Engelborghs, T. Luzyanina and D. Roose, Numerical bifurcation analysis of delay differential equations using DDE-BIFTOOL, ACM Trans. Math. Softw., 28 (1) (2002), pp. 1-21. [ELS01] K. Engelborghs, T. Luzyanina, G. Samaey, DDE-BIFTOOL v. 2.00: a Matlab package for bifurcation analysis of delay differential equations, Report TW330 (2001), Department of Computer Science, K. U. Leuven, Belgium. [ER02] K. Engelborghs and D. Roose, On Stability of LMS methods and characteristic roots of delay differential equations, SIAM J. Numer. Anal., 40 (2) (2002), pp. 629-650. [ER99] K. Engelborghs and 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), pp. 271-289. [FSD00] A. Fattouh, O. Sename and J.M. Dion, H∞ controller and observer design for linear systems with point and distributed delays, proceedings of 2nd IFAC workshop on Linear Time Delay Systems, 11-13 September 2000, Ancona, Italy. [HVL93] J.K. Hale and S.M. Verduyn Lunel, Introduction to Functional Differential Equations, Springer-Verlag, AMS series n. 99, New York, U.S.A., 1993. [LER03] T. Luzyanina, K. Engelborghs and D. Roose, Computing stability of differential equations with bounded distributed delays, Numer. Algorithms, 34 (1) (2003), pp. 41-66. [Nic00] S.I. Niculescu, Delay Effects on Stability: A Robust Control Approach, Monograph Springer, LNCIS n.269, New York, U.S.A., 2000. [Ric03] J.P. Richard, Time-delay systems: an overview of some recent advances and open problems, Automatica, 39 (2003), pp. 1667-1694.
PSEUDOSPECTRAL METHODS FOR CHARACTERISTIC ROOTS
15
[Tre00] L.N. Trefethen,Spectral methods in MATLAB, SIAM, Software - Environment - Tools series, Philadelphia, U.S.A., 2000.
