Structured Pseudospectra for Polynomial Eigenvalue Problems, with Applications.
Pseudospectra associated with the standard and generalized eigenvalue problems have been widely investigated in recent years. We extend the usual definitions in two respects, by treating the polynomial eigenvalue problem and by allowing structured perturbations of a type arising in control theory. We explore connections between structured pseudospectra, structured backward errors, and structured stability radii. Two main approaches for computing pseudospectra are described. One is based on a transfer function and employs a generalized Schur decomposition of the companion form pencil. The other, specific to quadratic polynomials, finds a solvent of the associated quadratic matrix equation and thereby factorizes the quadratic $\lambda$-matrix. Possible approaches for large, sparse problems are also outlined. A collection of examples from vibrating systems, control theory, acoustics, and fluid mechanics is given to illustrate the techniques.
Introduction
. Pseudospectra are an established tool for gaining insight into
the sensitivity of the eigenvalues of a matrix to perturbations. Their use is widespread
with applications in areas such as fluid mechanics, Markov chains, and control theory.
Most of the existing work is for the standard eigenproblem, although attention has also
been given to matrix pencils [4], [23], [33], [40], [46]. The literature on pseudospectra
is large and growing. We refer to Trefethen [41], [42], [43] for thorough surveys of
pseudospectra and their computation for a single matrix; see also the Web site [3].
In this work we investigate pseudospectra for polynomial matrices (or #-matrices)
0: m. We first define the #-pseudospectrum and obtain a computationally
useful characterization. We examine the relation between the backward
error of an approximate eigenpair of the polynomial eigenvalue problem associated
with (1.1), the #-pseudospectrum, and the stability radius. We consider both unstructured
perturbations and structured perturbations of a type commonly used in
control theory.
Existing methods for the computation of pseudospectra in the case
standard and generalized eigenvalue problems) do not generalize straightforwardly to
matrix polynomials. We develop two techniques that allow e#cient computation for
m > 1. A transfer function approach employs the generalized Schur decomposition of
the mn - mn companion form pencil. For the quadratic case an alternative
# Received by the editors May 1, 2000; accepted for publication (in revised form) by M. Chu
February 7, 2001; published electronically June 8, 2001. This work was supported by Engineering
and Physical Sciences Research Council grant GR/L76532.
http://www.siam.org/journals/simax/23-1/37145.html
Department of Mathematics, University of Manchester, Manchester, M13 9PL, England
(ftisseur@ma.man.ac.uk, http://www.ma.man.ac.uk/-ftisseur/, higham@ma.man.ac.uk,
http://www.ma.man.ac.uk/-higham/). The work of the second author was supported by a Royal
Society Leverhulme Trust Senior Research Fellowship.
solvent approach computes a solvent of the associated quadratic matrix equation
thereby factorizes the quadratic #-matrix; it works all
the time with n - n matrices once the solvent has been obtained. We give a detailed
comparison of these approaches and also outline techniques that can be e#ciently
used when n is so large as to preclude factorizations.
In the last section, we illustrate our theory and techniques on applications from
vibrating systems, control theory, acoustics, and fluid mechanics.
2. Pseudospectra.
2.1. Definition. The polynomial eigenvalue problem is to find the solutions
(x, #) of
where P (#) is of the form (1.1). If x #= 0 then # is called an eigenvalue and x the
corresponding right eigenvector; y #= 0 is a left eigenvector if y # P
of eigenvalues of P is denoted by #(P ). When Am is nonsingular P has mn finite
eigenvalues, while if Am is singular P has infinite eigenvalues. Good references for
the theory of #-matrices are [8], [20], [21], [37].
Throughout this paper we assume that P has only finite eigenvalues (and pseu-
doeigenvalues); how to deal with infinite eigenvalues is described in [16].
For notational convenience, we introduce
We define the #-pseudospectrum of P by
with
Here the # k are nonnegative parameters that allow freedom in how perturbations are
measured-for example, in an absolute sense (# k # 1) or a relative sense
By setting # unperturbed. The norm,
here and throughout, is any subordinate matrix norm. Occasionally, we will specialize
to the norm # p subordinate to the H-older vector p-norm.
When reduces to the
standard definition of #-pseudospectrum of a single matrix:
with #A# .
It is well known [43] that (2.4) is equivalent to
In the following lemma, we provide a generalization of this equivalence for the #-
pseudospectrum of P .
Lemma 2.1.
STRUCTURED PSEUDOSPECTRA FOR POLYNOMIAL EIGENPROBLEMS 189
Proof. Let S denote the set on the right-hand side of the claimed equality. We
first show that # (P ) implies # S. If # is an eigenvalue of P this is immediate,
so we can assume that # is not an eigenvalue of P and hence that P (#) is nonsingular.
Since
is singular, we have
so that # S.
Now let # S. Again we can assume that
with #, so that #x# = 1.
Then there exists a matrix H with
Lem. 6.3]). Let
y
y
and
We now apportion E between the A k by defining
where for complex z we define
Then
and #A k # k #, 0: m. Hence # (P ).
The characterization of the #-pseudospectrum in Lemma 2.1 will be the basis of
our algorithms for computing pseudospectra.
We note that for is the root neighborhood of the polynomial P
introduced by Mosier [28], that is, the set of all polynomials obtained by elementwise
perturbations of P of size at most #. This set is also investigated by Toh and Trefethen
[38], who call it the #-pseudozero set.
2.2. Connection with backward error. A natural definition of the normwise
backward error of an approximate eigenpair (x, #) of (2.1) is
and the backward error for an approximate eigenvalue # is given by
#) := min
#(x, #).
190 FRANC-OISE TISSEUR AND NICHOLAS J. HIGHAM
By comparing the definitions (2.3) and (2.6) it is clear that the #-pseudospectrum can
be expressed in terms of the backward error of # as
The following lemma gives an explicit expression for #(x, #) and #). This lemma
generalizes results given in [36] for the 2-norm and earlier in [5], [10] for the generalized
eigenvalue problem.
Lemma 2.2. The normwise backward error #(x, #) is given for x #= 0 by
p(|#x#
If # is not an eigenvalue of P then
#) =p(|#P (#)
Proof. It is straightforward to show that the right-hand side of (2.8) is a lower
bound for #(x, #). That the lower bound is attained is proved using a construction
for #A k similar to that in the proof of Lemma 2.1. The expression (2.9) follows on
using the equality, for nonsingular C # C n-n , min x#=0 #Cx#x#C
We observe that the expressions (2.7) and (2.9) lead to another proof of Lemma 2.1.
2.3. Structured perturbations. We now suppose that P (#) is subject to
structured perturbations that can be expressed as
with . The matrices D
and E are fixed and assumed to be of full rank, and they define the structure of the
perturbations; # is an arbitrary matrix whose elements are the free parameters. Note
that #A 0 , . , #Am in (2.10) are linear functions of the parameters in #, but that
not all linear functions can be represented in this form. We choose this particular
structure for the perturbations because it is one commonly used in control theory
[17], [18], [30] and it leads to more tractable formulae than a fully general approach.
Note, for instance, that the system
(which leads to a polynomial eigenvalue problem with may be interpreted as
a closed loop system with unknown static linear output feedback #; see Figure 2.1.
Note that unstructured perturbations are represented by the special case of (2.10)
with
For notational convenience, we introduce
Corresponding to (2.10) we have the following definition of structured backward error
for an approximate eigenpair (x, #):
STRUCTURED PSEUDOSPECTRA FOR POLYNOMIAL EIGENPROBLEMS 191
Fig. 2.1. Closed loop system with unknown static linear output feedback #.
and the backward error for an approximate eigenvalue is
#(x, #; D,E).
In the next result we use a superscript "+" to denote the pseudo-inverse [9].
Lemma 2.3. The structured backward error #(x, #; D,E) in the Frobenius norm
is given by
#F
if the system
is consistent; otherwise # F (x, #; D,E) is infinite.
Proof. It is immediate that # F (x, #; D,E) is the Frobenius norm of the minimum
Frobenius norm solution to (2.13). The result follows from the fact that
is the solution of minimum Frobenius norm to the consistent system
sect. 3.4.8].
To gain some insight into the expression (2.12) we consider the case of unstructured
but weighted perturbations, as in (2.11) but with
The system (2.13) is now trivially consistent and (2.12) gives
using the fact that #ab . The expression (2.14) di#ers
from that for #(x, #) in (2.8) for the 2-norm only by having the 2-norm of the vector
rather than the 1-norm in the denominator.
Lemma 2.4. If # is not an eigenvalue of P (#) then the structured backward error
#; D,E) is given by
(2.
192 FRANC-OISE TISSEUR AND NICHOLAS J. HIGHAM
Proof. We have
min
The companion form of P (#P (#) is given by
where
. 0
I
-A
I
I
I
Am
and
#Am
# .
As
# .
Then, using the identity det(I both AB and
BA are defined [47, p. 54],
det(P (#P
using [45, Lem. 1],
we have
But it is easily verified that
D,
STRUCTURED PSEUDOSPECTRA FOR POLYNOMIAL EIGENPROBLEMS 193
so that
We define the structured #-pseudospectrum by
Analogously to the unstructured case, #
so from Lemma 2.4 we have
which is a generalization of a result of Hinrichsen and Kelb [17, Lem. 2.2] for the
#-pseudospectrum of a single matrix.
2.4. Connection between backward error and stability radius. In many
mathematical models (e.g., those of a dynamical system) it is required for stability
that a matrix has all its eigenvalues in a given open subset C g # of the complex
plane. Various stability radii have been defined that measure the ability of a matrix
to preserve its stability under perturbations.
We partition the complex plane C into two disjoint subsets C g and C b , with
# an open set.
Consider perturbations of the form in (2.10). Following Pappas and Hinrichsen [30]
and Genin and Van Dooren [7], we define the complex structured stability radius of
the #-matrix P with respect to the perturbation structure (D, E) and the partition
by
r
Let #C b be the boundary of C b . By continuity, we have
r
#C s-t { #(P (#) +D#E(#C b # }
#; D,E).
Thus we have expressed the stability radius as an infimum of the eigenvalue backward
error. Using Lemma 2.4 we obtain the following result.
Lemma 2.5. If # is not an eigenvalue of P then
r
and for unstructured perturbations and the p-norm we have
r
The result for the unstructured case in the second part of this lemma is also
obtained by Pappas and Hinrichsen [30, Cor. 2.4] and Genin and Van Dooren [7,
Thm. 2].
194 FRANC-OISE TISSEUR AND NICHOLAS J. HIGHAM
3. Computation of pseudospectra. In this section, we consider the computation
of # (P ), concentrating mainly on the 2-norm. We develop methods for unstructured
perturbations and show how they can be extended to structured perturbations
of the form in (2.10).
Lemma 2.1 shows that the boundary of # (P ) comprises points z for which the
scaled resolvent norm p(|z|)#P (z) Hence, as for pseudospectra of a
single matrix, we can obtain a graphical representation of the pseudospectra of a
polynomial eigenvalue problem by evaluating the scaled resolvent norm on a grid of
points z in the complex plane and sending the results to a contour plotter. We refer
to Trefethen [42] for a survey of the state of the art in computation of pseudospectra
of a single matrix.
The region of interest in the complex plane will usually be determined by the
underlying application or by prior knowledge of the spectrum of P . In the absence
of such information we can select a region guaranteed to enclose the spectrum. If
Am is nonsingular (so that all eigenvalues are finite) then by applying the result
(A)| #A#" to the companion form (2.16) we deduce that
for any p-norm. Alternatively, we could bound max j |# j (P )| by the largest absolute
value of a point in the numerical range of P [24], but computation of this number
is itself a nontrivial problem. For much more on bounding the eigenvalues of matrix
polynomials see [15].
For the 2-norm, #P (z) denotes the smallest
singular value. If the grid is # and # min is computed using the Golub-Reinsch
SVD algorithm then the whole computation requires roughly
which is prohibitively expensive for matrices of large dimension and a fine grid. Using
the fact that # min (P (z)) is the square root of # min (P (z) # P (z)), we can approximate
with the power iteration or Lanczos iteration applied to P (z)
In the case of a single matrix, Lui [25] introduced the idea of using the Schur form of
A in order to speed up the computation of # min Unfortunately,
for matrix polynomials of degree m # 2 no analogue of the Schur form exists (that
is, at most two general matrices can be simultaneously reduced to triangular form).
We therefore look for other ways to e#ciently evaluate or approximate #P (z) -1 # for
many di#erent z.
3.1. Transfer function approach. The idea of writing pseudospectra in terms
of transfer functions is not new. Simoncini and Gallopoulos [34] used a transfer function
framework to rewrite most of the techniques used to approximate #-pseudospectra
of large matrices, yielding interesting comparisons as well as better understanding of
the techniques. Hinrichsen and Kelb [17] investigated structured pseudospectra of a
single matrix with perturbations of the form in (2.10), and they expressed the structured
#-pseudospectrum in terms of a transfer function.
Consider the equation
STRUCTURED PSEUDOSPECTRA FOR POLYNOMIAL EIGENPROBLEMS 195
It can be rewritten as
wm
where F and G are defined in (2.16). Hence
-I
# u.
Since this equation holds for all u, it follows that
-I
# .
This equality can also be deduced from the theory of #-matrices [21, Thm. 14.2.1].
We have thus expressed the resolvent in terms of a transfer function.
In control theory, P (z) -1 corresponds to the transfer function of the linear time-invariant
multivariate system described by
I
Several algorithms have been proposed in the literature [22], [27] to compute transfer
functions at a large number of frequencies, most of them assuming that G = I. Our
objective is to e#ciently compute the norm of the transfer function, rather than to
compute the transfer function itself.
For structured perturbations we see from (2.17) that the transfer function P (z)
is replaced by
E(z)P (z)
-D
# .
All the methods described below for the dense case are directly applicable with obvious
changes.
We would like a factorization of F - zG that enables e#cient evaluation or application
of di#erent z. There are various possibilities, including,
when G is nonsingular,
196 FRANC-OISE TISSEUR AND NICHOLAS J. HIGHAM
is a Schur decomposition, with W unitary and T upper tri-
angular. However this approach is numerically unstable when G is ill conditioned. A
numerically stable reduction is obtained by computing the generalized Schur decom-
position
where W and Z are unitary and T and S are upper triangular. Then
-I
# .
Hence once the generalized Schur decomposition has been computed, we can compute
x at a cost of O((mn) 2 ) flops, since T - zS is triangular of
dimension mn. For the 2-norm we can therefore e#ciently approximate #P (z)
using inverse iteration or the inverse Lanczos iteration, that is, the power method or
the Lanczos method applied to P (z)
The cost of the computation breaks into two parts: the cost of the initial transformations
and the cost of the computations at each of the # 2 grid points. Assuming
that (3.3) is computed using the QZ algorithm [9, Sec. 7.7] and the average number
of power method or Lanczos iterations per grid point is k, the total cost is about
For the important special case eigenvalue problem), this cost is
Comparing with (3.1) we see that this method is a significant improvement over the
SVD-based approach for a su#ciently fine grid and a small degree m.
For the 2-norm note that, because of the two outer factors in (3.4), we cannot
discard the unitary matrices Z and W , unlike in the analogous expression for the
resolvent of a single matrix in the standard eigenproblem. For the 1- and #-norms
we can e#ciently estimate #P (z) using the algorithm of Higham and Tisseur [14],
which requires only the ability to multiply matrices by P (z) -1 and P (z) -# .
An alternative to the generalized Schur decomposition is the generalized
Hessenberg-triangular form, which di#ers from (3.3) in that one of T and S is upper
Hessenberg. The Hessenberg form is cheaper to compute but more expensive to work
with. It leads to a smaller overall flop count when k# 2 >
# 25mn.
3.2. Factorizing the quadratic polynomial. The transfer function-based method
of the previous section has the drawback that it factorizes matrices of dimension m
times those of the original polynomial matrix. We now describe another method,
particular to the quadratic case, that does not increase the size of the problem.
Suppose we can find a matrix S such that A 2 S 2 +A 1 S+A that is, a solvent
of the quadratic matrix equation A 2
If we compute the Schur decomposition
STRUCTURED PSEUDOSPECTRA FOR POLYNOMIAL EIGENPROBLEMS 197
and the generalized Schur decomposition
then
so a vector can be premultiplied by Q(z) -1 or its conjugate transpose in O(n 2 ) flops
for any z. Moreover, for the 2-norm we can drop the outer Q and W # factors in (3.7),
by unitary invariance, and hence we do not need to form W . For the 2-norm, the
total cost of this method is
where c S is the cost of computing a solvent and we have assumed that we precompute
Z. Comparing this flop count with (3.5) we see that the cost per grid point of the
solvent approach is much lower.
The success of this method depends on two things: the existence of solvents
and being able to compute one at a reasonable cost. Some su#cient conditions for
the existence of a solvent are summarized in [13]. In particular, for an overdamped
problem, one for which A 2 and A 1 are Hermitian positive definite, A 0 is Hermitian
positive semidefinite, and a solvent is
guaranteed to exist.
Various methods are available for computing solvents [12], [13]. One of the most
generally useful is Newton's method, optionally with exact line searches, which requires
a generalized Sylvester equation in n-n matrices to be solved on each iteration,
at a total cost of about 56n 3 flops per iteration. If Newton's method converges within
8 iterations or so, so that c S # 448n 3 flops, this approach is certainly competitive in
cost with the transfer function approach.
When there is a gap between the n largest and n smallest eigenvalues ordered by
modulus, as is the case for overdamped problems [20, Sec. 7.6], Bernoulli iteration is an
e#cient way of computing the dominant or minimal solvent S [13]. If t iterations are
needed for convergence to the dominant or minimal solvent then the cost of Bernoulli
iteration is about c flops. Bernoulli iteration converges only linearly, but
convergence is fast if the eigenvalue gap is large.
A third approach to computing a solvent is to use a Schur method from [13],
based on the following theorem. Let F and G be defined as in (2.16), so that
-A
# .
Theorem 3.1 (Higham and Kim [13]). All solvents of Q(X) are of the form
is a generalized Schur decomposition with Q and Z unitary and T and S upper tri-
angular, and where all matrices are partitioned as block 2 - 2 matrices with n - n
blocks.
The method consists of computing the generalized Schur decomposition (3.9) by
the QZ algorithm and then forming
11 . The generalized Schur decomposition
may need to be reordered in order to obtain a nonsingular Z 11 . Note that the
198 FRANC-OISE TISSEUR AND NICHOLAS J. HIGHAM
unitary factor Q does not need to be formed. For this method, c
where the constant r depends on the amount of reordering required. From (3.8), the
total cost is now
which is much more favorable than the cost (3.5) of the transfer function method.
For higher degree polynomials we can generalize this approach by attempting to
linear factors by recursively computing solvents. However, for degrees
greater than 2 classes of problem for which a factorization into linear factors exists
are less easily identified and the cost of Newton's method (for example) is much higher
than for
3.3. Large-scale computation. All the methods described above are intended
for small- to medium-scale problems for which Schur and other reductions are pos-
sible. For large, possibly sparse, problems, di#erent techniques are necessary. These
techniques can be classified into two categories: those that project to reduce the size
of the problem and then compute the pseudospectra of the reduced problem, and
those that approximate the norm of the resolvent directly.
3.3.1. Projection approach. For a single matrix, A, Toh and Trefethen [39]
and Wright and Trefethen [48] approximate the resolvent norm by the Arnoldi method;
that is, they approximate or by # min (
where Hm is the square Hessenberg matrix of dimension m # n obtained from the
Arnoldi process and
Hm is the matrix Hm augmented by an extra row. Simoncini
and Gallopoulos [34] show that a better but more costly approximation is obtained by
approximating #(A- zI)
is the orthonormal
basis generated during the Arnoldi process. These techniques are not applicable
to the polynomial eigenvalue problem of degree larger than one because of the lack of
a Schur form for the Arnoldi method to approximate.
A way of approximating #P (z) -1 # for all z is through a projection of P (z)
a lower dimensional subspace. Let V k be an n - k matrix with orthonormal columns.
We can apply one of the techniques described in the previous sections to compute
pseudospectra of the projected polynomial eigenvalue problem
A possible choice for V k is an orthonormal basis of k selected linearly independent
eigenvectors of P (#). In this case,
P (#) is the matrix representation of the projection
of P (#) onto the subspace spanned by the selected eigenvectors. The eigenvectors can
be chosen to correspond to parts of the spectrum of interest and can be computed
using the Arnoldi process on the companion form pencil (F, G) or directly on P (#)
with the Jacobi-Davidson method or its variants [26], [35]. In the latter case, the
during the Davidson process.
3.3.2. Direct approach. This approach consists of approximating #P (z)
at each grid point z. Techniques analogous to those used for single matrices can be
applied, such as the Lanczos method applied to P (z) # P (z) or its inverse. We refer
the reader to [42] for more details and further references.
4. Applications and numerical experiments. We give a selection of applications
of pseudospectra for polynomial eigenvalue problems, using them to illustrate
the performance of our methods for computing pseudospectra. All our examples are
for 2-norm pseudospectra.
STRUCTURED PSEUDOSPECTRA FOR POLYNOMIAL EIGENPROBLEMS 199
4.1. The wing problem. The first example is based on a quadratic polynomial
A 0 from [6, Sec. 10.11], with numerical values modified as in
[20, Sec. 5.3]. The eigenproblem for Q(#) arose from the analysis of the oscillations
of a wing in an airstream. The matrices are
17.6 1.28 2.89
1.28 0.824 0.413
7.66 2.45 2.1
# .
The left plot in Figure 4.1 shows the boundaries of #-pseudospectra with perturbations
measured in the absolute sense . The
eigenvalues are plotted as dots. Another way of approximating a pseudospectrum
is by random perturbations of the original matrices [41]. We generated 200 triples
of complex random normal perturbation matrices (#A 1 ,
1: 3. In the right plot of Figure 4.1 are superimposed as small dots the
eigenvalues of the perturbed polynomials # 2
The solid curve marks the boundary of the #-pseudospectrum for
pictures show that the pair of complex eigenvalues are much more
sensitive to perturbations than the other two complex pairs.
The eigenvalues of Q(#) are the same as those of the linearized problem A- #I,
where
-A
# .

Figure

4.2 shows boundaries of #-pseudospectra for this matrix, for the same # as in

Figure

4.1. Clearly, the #-pseudospectra of the linearized problem (4.1) do not give
useful information about the behavior of the eigensystem of Q(#) under perturbations.
This emphasizes the importance of defining and computing pseudospectra for the
quadratic eigenvalue problem in its original form.
4.2. Mass-spring system. We now consider the connected damped mass-spring
system illustrated in Figure 4.3. The ith mass of weight m i is connected to the (i+1)st
mass by a spring and a damper with constants k i and d i , respectively. The ith mass
is also connected to the ground by a spring and a damper with constants # i and # i ,
respectively. The vibration of this system is governed by a second-order di#erential
equation
d
dt
where the mass matrix diagonal, and the damping matrix
C and sti#ness matrix K are symmetric tridiagonal. The di#erential equation leads
to the quadratic eigenvalue problem
In our experiments, we took all the springs (respectively, dampers) to have the same
constant except the first and last, for which the constant
is 2# (respectively, 2# ), and we took m i # 1. Then
200 FRANC-OISE TISSEUR AND NICHOLAS J. HIGHAM
Fig. 4.1. Wing problem. Left: # (Q), for # [10 approximation to #-
pseudospectrum with
Fig. 4.2. Wing problem. # (A), for A in (4.1) with # [10
and the quadratic eigenvalue problem is overdamped. We take an
of freedom mass-spring system over a 100 - 100 grid. A plot of the pseudospectra is
given in Figure 4.4.
For this problem we compare all the methods described. In the solvent approach
exact line searches were used in Newton's method and no reordering was used in the
generalized Schur method. The solvents from the Bernoulli and Schur methods were
refined by one step of Newton's method. The Bernoulli iteration converged in 12
iterations while only 6 iterations were necessary for Newton's method. The Lanczos
inverse iteration converged after 3 iterations on average. In Table 4.1 we give the estimated
flop counts, using the formulae from section 3, together with execution times.
The computations were performed in MATLAB 6, which is an excellent environment
for investigating pseudospectra. While the precise times are not important, the con-
STRUCTURED PSEUDOSPECTRA FOR POLYNOMIAL EIGENPROBLEMS 201
d i-1
d n-1
Fig. 4.3. An n degree of freedom damped mass-spring system.
Fig. 4.4. Pseudospectra of a 250 degree of freedom damped mass-spring system on a 100 - 100
grid.
clusion is clear: in this example, the three solvent-based methods are much faster
than the SVD and transfer function methods. (The high speed of the SVD method
relative to its flop count is attributable to MATLAB's very e#cient svd function.)
4.3. Acoustic problem. Acoustic problems with damping can give rise to large
quadratic eigenvalue problems (4.2), where, again, M is the mass matrix, C is the
damping matrix, and K the sti#ness matrix. We give in Figure 4.5 the sparsity
pattern of the three matrices M , C, and K of order 107 arising from a model of a
speaker box [1]. These matrices are symmetric and the sparsity patterns of M and
K are identical. There is a large variation in the norms: #M#
We plot in Figure 4.6 pseudospectra with perturbations measured in both an
absolute sense a relative sense
together with pseudospectra of the corresponding standard eigenvalue
problem of the form (4.1). The eigenvalues are all pure imaginary and are marked by
dots on the plot. The two first plots are similar, both showing that the most sensitive

Table
Comparison in terms of flops and execution time of di#erent techniques.
Method Estimated cost in flops Execution time
Golub-Reinsch SVD 26747n 3 102 min
Transfer function 3408n 3 106 min
Solvent: Schur 1677n 3 37 min
Matrices M and K
Matrix C
Fig. 4.5. Sparsity patterns of the three 107 - 107 matrices M,C, and K of an acoustic problem.
eigenvalues are located at the extremities of the spectrum; the contour lines di#er
mainly around the zero eigenvalue. The last plot is very di#erent; clearly it is the
eigenvalues close to zero that are the most sensitive to perturbations of the standard
eigenproblem form.
We mention that for this problem we have been unable to compute a solvent.
4.4. Closed loop system. In multi-input and multioutput systems in control
theory the location of the eigenvalues of matrix polynomials determine the stability
of the system. Figure 4.7 shows a closed-loop system with feedback with gains 1 and
The associated matrix polynomial is given by
# .
We are interested in the values of # for which P (z) has all its eigenvalues inside the
unit circle. By direct calculation with det(P (z)), using the Routh array, for example,
it can be shown that P (z) has all its eigenvalues inside the unit circle if and only if
# < 0.875.
The matrix P (z) can be viewed as a perturbed matrix polynomial with structured
perturbations:
STRUCTURED PSEUDOSPECTRA FOR POLYNOMIAL EIGENPROBLEMS 203
-4000 -3000 -2000 -1000 0 1000 2000 3000 4000
-4000 -3000 -2000 -1000 0 1000 2000 3000 4000
-4000 -3000 -2000 -1000 0 1000 2000 3000 4000
Fig. 4.6. Acoustic problem, Perturbations measured in an absolute sense
(top left) and relative sense (top right). Pseudospectra of the equivalent standard eigenvalue problem
are shown at the bottom.
where
I
zI
z 2 I
# .
We show in Figure 4.8 the structured pseudospectra as defined by (2.17). The dashed
lines mark the unit circle. Since the outermost contour has value just
touches the unit circle, this picture confirms the value for the maximal # that we
obtained analytically.
4.5. The Orr-Sommerfeld equation. The Orr-Sommerfeld equation is a linearization
of the incompressible Navier-Stokes equations in which the perturbations
in velocity and pressure are assumed to take the form
y,
204 FRANC-OISE TISSEUR AND NICHOLAS J. HIGHAM
z
z
Fig. 4.7. Closed-loop system with feedback gains 1 and 1
10.60.20.20.60.80.36Fig. 4.8. Structured pseudospectra of a closed-loop system with one-parameter feedback.
where # is a wavenumber and # is a radian frequency. For a given Reynolds number
R, the Orr-Sommerfeld equation may be written
We consider plane Poiseuille flow between walls at
in the streamwise x direction, for which the boundary conditions are
For a given real value of R, the boundary conditions will be satisfied only for certain
combinations of values of # and #. Two cases are of interest.
Case 1. Temporal stability. If # is fixed and real, then (4.3) is linear in the
parameter # and corresponds to a generalized eigenvalue problem. The perturbations
STRUCTURED PSEUDOSPECTRA FOR POLYNOMIAL EIGENPROBLEMS 205
are periodic in x and grow or decay in time depending on the sign of the imaginary
part of #. This case has been studied with the help of pseudospectra by Reddy,
Schmid, and Henningson [32].
Case 2. Spatial stability. For most real flows, the perturbations are periodic in
time, which means that # is real. Then the sign of the imaginary part of # determines
whether the perturbations will grow or decay in space. In this case, the parameter
is #, which appears to the fourth power in (4.3), so we obtain a quartic polynomial
eigenvalue problem. Bridges and Morris [2] calculated the spectrum of (4.3) using a
finite Chebyshev series expansion of # combined with the Lanczos tau method and
they computed the spectrum of the quartic polynomial by two methods: the QR
algorithm applied to the corresponding standard eigenvalue problem in companion
form, and Bernoulli iteration applied to determine a minimal solvent and hence to
obtain the n eigenvalues of minimal modulus.
For our estimation of the pseudospectra of the Orr-Sommerfeld equation we use
a Chebyshev spectral discretization that combines an expansion in Chebyshev polynomials
and collocation at the Chebyshev points with explicit enforcement of the
boundary conditions. We are interested in the eigenvalues # that are the closest to
the real axis, and we need Im(#) > 0 for stability. The linear eigenvalue problem
(Case 1) has been solved by Orszag [29]. The critical neutral point corresponding
to # and # both real for minimum R was found at
the frequency our calculations we set R and # to these
values and we computed the modes #, taking which gives matrices of order
1. The first few modes are plotted in Figure 4.9. For the first mode we obtained
which compares favorably with the result of Orszag. Figure
4.10 shows the pseudospectra in a region around the first few modes on a 100-100
grid, with # since A 4 is the identity matrix and is not
subject to uncertainty. The plot shows that the first mode is very sensitive. Interest-
ingly, the second and subsequent modes are almost as sensitive, with perturbations of
in the matrix coe#cients being su#cient to move all these modes across
the real axis, making the flow unstable. The pseudospectra thus give a guide to the
accuracy with which computations must be carried out for the numerical approximations
to the modes to correctly determine the location of the modes. For more on the
interpretation of pseudospectra for this problem, see [32] and [44].
Again, for comparison we computed the pseudospectra of the corresponding standard
eigenvalue problem. The picture was qualitatively similar, but the contour levels
were several orders of magnitude smaller, thus not revealing the true sensitivity of the
problem.
206 FRANC-OISE TISSEUR AND NICHOLAS J. HIGHAM
Fig. 4.9. The first few modes of the spectrum of the Orr-Sommerfeld equation for
109.39.058.88.99.159.59.5Fig. 4.10. Pseudospectra of the Orr-Sommerfeld equation for
STRUCTURED PSEUDOSPECTRA FOR POLYNOMIAL EIGENPROBLEMS 207



--R



http://www.


Elementary Matrices and Some Applications to Dynamics and Di
Stability Radii of Polynomial Matrices
Matrix Polynomials
Matrix Computations
Structured backward error and condition of generalized eigenvalue problems
Accuracy and Stability of Numerical Algorithms
Solving a Quadratic Matrix Equation by Newton's Method with Exact Line Searches
Numerical analysis of a quadratic matrix equation
A block algorithm for matrix 1-norm estimation
Bounds for Eigenvalues of Matrix Polynomials
More on Pseudospectra for Polynomial Eigenvalue Problems and Applications in Control Theory
Spectral value sets: A graphical tool for robustness analysis
Real and complex stability radii: A survey
Numerical solution of matrix polynomial equations by Newton's method

The Theory of Matrices


Numerical range of matrix polynomials
Computation of pseudospectra by continuation
Locking and restarting quadratic eigenvalue solvers
A determinant identity and its application in evaluating frequency response matrices
Root neighborhoods of a polynomial
Accurate solution of the Orr-Sommerfeld stability equation
Robust stability of linear systems described by higher order
Generalized Inverse of Matrices and Its Applications
Pseudospectra of the Orr-Sommerfeld operator

Transfer functions and resolvent norm approximation of large matrices

Backward error and condition of polynomial eigenvalue problems
The quadratic eigenvalue problem
Pseudozeros of polynomials and pseudospectra of companion matrices
Calculation of pseudospectra by the Arnoldi iteration
Portraits Spectraux de Matrices: Un Outil d'Analyse de la Stabilit-e
Pseudospectra of matrices
Computation of pseudospectra
Spectra and pseudospectra
Hydrodynamic stability without eigenvalues
On stability radii of generalized eigenvalue problems
Pseudospectra for matrix pencils and stability of equilibria
The Algebraic Eigenvalue Problem

--TR

--CTR
Kui Du, Note on structured indefinite perturbations to Hermitian matrices, Journal of Computational and Applied Mathematics, v.202 n.2, p.258-265, May, 2007
Graillat, A note on structured pseudospectra, Journal of Computational and Applied Mathematics, v.191 n.1, p.68-76, 15 June 2006
Kui Du , Yimin Wei, Structured pseudospectra and structured sensitivity of eigenvalues, Journal of Computational and Applied Mathematics, v.197 n.2, p.502-519, 15 December 2006
Kirk Green , Thomas Wagenknecht, Pseudospectra and delay differential equations, Journal of Computational and Applied Mathematics, v.196 n.2, p.567-578, 15 November 2006
structured perturbations;matrix polynomial;solvent;orr-sommerfeld equation;stability radius;pseudospectrum;quadratic matrix equation;backward error;transfer function;lambda-matrix;polynomial eigenvalue problem