ABSTRACT PERTURBATIONS OF STABILITY MATRICES WITH APPLICATIONS TO LOTKA-VOLTERRA MODELS OF ECOSYSTEMS BY Gary William Harrison A system of differential equations x’ = Ax + g(x), where g(x) = 0(HXH), has an asymptotically stable equilibrium point at the origin if A is a stability matrix, that is, if all the eigenvalues of A have negative real part. In theory this condition has been used to determine the stability of Lotka-Volterra models of ecosystems, but in practice it breaks down because A can only be estimated. Methods are develOped to establish that a matrix is stable when it is only known approximately. Specifically, bounds on e are found for the perturbed matrix A + 63 to be stable when A is stable. When B is a rank one matrix, necessary and sufficient bounds are derived from the Routh-Hurwitz stability criterion. For a general matrix B, sufficient bounds are derived T from the Lyapunov equation A S + SA = -Q. Perturbations eB satisfying these latter bounds are shown to form an Gary William Harrison open convex set, leading to a method to compute open convex neighborhoods of a stable matrix A which con- tain only stable matrices. Similar methods yield sufficient conditions for the bordered matrix [AT b] c d to be stable when A is stable. (1+1) (1+1)N = Iterative methods of the form NTS + S NTS(1) + S(l)N + wIATS(l)+S(l)A+Q) are investigated to solve the Lyapunov equation. They can converge only if (1+1) 3 Nxm + MAX“) the vector scheme Nx +q) converges, and no faster. Sufficient conditions are found for con- vergence of several methods corresponding to different choices of N, and examples are given showing that a block Seidel type method is usually most efficient and better than previously known methods of solution for large matrices. The implications for stability of Lotka- Volterra and other ecosystem models are discussed, and conditions are found for stability to be preserved when a new species is added. A domain of attraction for the equilibrium point, which in some cases is the entire positive quadrant, is found using Lyapunov functions. When coefficients in the model are time varying and the derivative of the variations is small enough, solutions are shown to follow close to a moving critical point. If in addition the coefficients are periodic, there is a periodic solution. PERTURBATIONS OF STABILITY MATRICES WITH APPLICATIONS TO LOTKA-VOLTERRA MODELS OF ECOSYSTEMS BY Gary William Harrison A THESIS Submitted to Michigan State University in partial fulfillment of the requirements for the degree of DOCTOR OF PHILOSOPHY Department of Mathematics 1975 To Andrea, Martha, and Bryce ii ACKNOWLEDGMENT S I wish to thank those who have contributed to this thesis. Dr. Robert Rosen introduced me to bio- mathematics and continued to give me helpful suggestions and references even after leaving Michigan State University. Dr. J. Sutherland Frame took over as my advisor when Dr. Rosen left and gave me invaluable guidance during the develOpment and writing of the thesis. Dr. P.K. WOng, Dr. Mary Winter, and Dr. M. Tomber completed my doctoral committee. Many others have provided re- ferences or contributed through the classes they taught. Less tangible but no less real has been the support of family and friends, especially my wife, Andrea and my children Martha and Bryce. I might have written a thesis under a different advisor, but without Andrea's cooperation any thesis would have been impossible. I also wish to thank my mother and father, my brother, Evan, and my office mate Mike McGrew. My research was supported financially by the National Science Foundation and by the Department of Mathematics at Michigan State University. iii TABLE OF CONTENTS Chapter Page I. INTRODUCTION AND BACKGROUND 1 1.1 Introduction . . . . . . . . . . . . . . l 1.2 Lotka-Volterra Model of an Ecosystem . 4 1.3 Description of Problem . . . . . . . . . 6 1.4 Notation . . . . . . . . . . . . . 9 1.5 Background . . . . . . . . . . . . . . . 11 II. APPLYING LYAPUNOV'S EQUATION 21 2.1 The Lyapunov Equation . . . . . . . . . 21 2.2 The Lyapunov Equation and the Symmetric Part of a Matrix . . . . . . . . . . 26 2.3 Stability Preserving Perturbations by Lyapunov's Equation . . . . . . . . . . 32 2.4 Convexity of S-Permissible Perturbations . . . . . . . . . . . . . 40 2.5 Bordering a Stability Matrix . . . . . . 46 III. ITERATIVE SOLUTIONS OF THE LYAPUNOV EQUATION 50 3.1 General Iterative Procedures . . . . . . 51 3.2 Rate of Convergence . . . . . . . . . . 53 3.3 Practical Iterative Procedures . . . . . 56 3.4 Simple Iteration . . . . . . . . . . . . 60 3.5 JacObi Iteration . . . . . . . . . . . . 61 iv Chapter Page 3.6 Seidel Iteration . . . . . . . . . . . . 65 3.7 Block Seidel Iteration . . . . . . . . . 66 3.8 Examples and Comparisons . . . . . . . . 70 IV. APPLYING THE ROUTH—HURWITZ STABILITY CRITERION 74 4.1 The Stability Criterion of Routh and Hurwitz . . . . . . . . . . . . . . 74 4.2 Applications to Perturbations . . . . . 77 4.3 Examples . . . . . . . . . . . . . . . . 81 V. APPLICATIONS TO LOTKA-VOLTERRA AND OTHER MODELS OF ECOSYSTEMS 87 5.1 The Community Matrix . . . . . . . . . . 87 5.2 D—Stable Volterra Matrices . . . . . . . 97 5.3 Addition of a New Species to an Ecosystem . . . . . . . . . . . . . . . 100 5.4 Domain of Attraction by Lyapunov Functions . . . . . . . . . . . . . . . 104 5.5 Time Varying Coefficients in the Lotka— Volterra Model . . . . . . . . . . . . . 109 CONCLUSIONS 120 BIBLIOGRAPHY 126 Figure 1. LIST OF FIGURES Graph of 11+wxi12 versus w, where In <...< 12 < xl<<() are real eigen- values of N -1 A. Graph of ll—wiiwynl2 versus w, where _-I_-iyn are the eigenvalues of D'YLR with 0 g y. Solution and critical values of equation (83) with w = 2V. Solution and critical values of equation (83) with w = w. Solution and critical values of equation (83) with w = 2. vi Page 55 64 115 116 117 CHAPTER I INTRODUCTION AND BACKGROUND 1.1. Introduction. When is an ecosystem stable? This is one of the most frequently asked questions in ecology today. In fact, the question of stability is one of the central questions in any mathematical analysis of a dynamical system. Let x(t) be the vector of "state variables", i.e. the variables which characterize the system, and assume (1) x’(t) = f(x,t) . Definition 1.1. If f is a function of x only, the system defined by (l) is said to be autonomous. Definition 1.2. A vector XE is called an equilibrium point or a critical point of (1) if f(xE,t) = O for all t. We note that an equilibrium point XE is a solution of (1). In this thesis we will use the following definition of stability introduced by Lyapunov [25]: Definition 1.3. An equilibrium point XE is stable if for every 6 > O and every t .2 0, there is O a 6 such that any solution of (1) with ”x(to) - XE” < 5 satisfies “x(t) - XE” < e for all t'Z to. Definition 1.4. An equilibrium point XE is asymptotically stable if it is stable and for every t0.2 0 there is a 6 such that ”x(to) - XE” < 5 implies x(t) converges to xE as t goes to infinity. The stability of a general solution x(t) of (1) is that of the equilibrium point 0 under the change of variables y(t) = x(t) - x(t). One of the standard techniques to prove stability is given by the following theorem, found in most advanced texts on differential equations. (See for example, Hahn [12], Hale [13], or Rosen [33].) Definition 1.5. A function g(x) is said to be 0(HxH), written g(x) = 0(HxH), if for every 6 > 0 there is a 5 such that “X“ < 6 implies “g(x)” < eHXH. Theorem 1.6 (Lyapunov, 1893). Let (2) x’=Ax+g(x,t), where A is a matrix, and g(x,t) = o(HxH). Then if the 0 vector is asymptotically stable for (3) x' = Ax the 0 vector is asymptotically stable for (2). Theorem 1.6 can be generalized to have nonconstant matrices A(t) in (2) and (3) by introducing the con- cept of uniform asymptotic stability, but we will not need this generality. Theorem 1.7. The system (3) is asymptotically stable if and only if all the eigenvalues of A have negative real part. Theorem 1.7 is our motivation to study matrices whose eigenvalues have negative real part. To avoid repetition, we make the following definition: Definition 1.8. A is a stability matrix (or a matrix A is said to be stable) if all the eigenvalues of A have negative real part. Lyapunov gave an alternative method to prove stability, known as his second or direct method. The function V in Theorem 1.9 is called a Lyapunov function. For a proof see Hahn [12] or Yoshizawa [44]. Theorem 1.9. Let 0 be an equilibrium point of (l) and assume there is a differentiable function V(x,t) defined on a neighborhood of O in the domain of f, satisfying (i) V(O,t) EEO: (ii) a(x) g V(x,t), for some continuous positive definite function a(x), (iii) v’(x(t) ,t) g 0. Then 0 is a stable equilibrium point of (1). If f(x,t) is bounded when x belongs to a compact set and (iii) is replaced by (iv) V'(x(t),t) g_—c(x), for some continuous positive definite function c(x), then 0 is asymptotically stable equilibrium point of (1). l;g, Lotka-Volterra Model of an Ecosystem. To study the stability of an ecosystem, we will use the Lotka-volterra model introduced for a two species system by Lotka in 1925 [24], and independently for the general n—species system by Volterra in 1926 [40]. (See also volterra [41].) Let pi(t) be the population (or, as is often preferred, the pOpulation density) of the ith species at time t in an ecosystem with n species. Then the growth rate p{(t) is given by n (4) p; = p.1(ri + .2) aijpj) j—l where ri is the intrinsic growth rate and dij is the interaction coefficient measuring the effect of species j on the growth of species i. If species j preys on is negative and a.. species i, 31 aij p031tive. If species 1 and j compete for the same resource dij and aji are both negative. The terms aii represent self-interaction or self-competition. They were omitted by both Lotka and Volterra in their original formulation, but there is good biological reason to include them and assume that they are less than or equal to zero. We will discuss this system and the significance of the coeffi- cients more fully in Chapter 5. We also refer the reader to Rosen [33], who discusses the relationship of (4) to the "law of mass action" of chemistry. We assume that the coefficients in (4) are such that there is an equilibrium point pE satisfying 0 i = l,2,...,n E (5) ri + ZEQijpj with p? > O for all j. Clearly, negative values for pi have no meaning. There are other equilibrium points where pi = O for some i and (5) is satisfied for the other values of i, but a zero value for pi corresponds to extinction for that species. To investigate the stability of pE we make the change of co-ordinates (6a) x = p - pE. The stability properties of pE are the same as those of the zero vector in n (6b) x; = Z) pgd..x. + _Z} x.d..x., i = l,2,...,n. The zero vector is asymptotically stable in (6) if the matrix _ _ E (7) A — [ai.]. aij — pioij is a stability matrix. Several authors (May [28], Rosen [33], and Strobeck [36]) have discussed the stability of pE by examining the eigenvalues of A. It would appear that for an actual ecosystem the mathematical analysis of stability would be rather straight- forward, since several numerical techniques for finding the eigenvalues of a matrix are well developed. But this presupposes an exact knowledge of the coefficients r. 1 and aij' whereas the biologist will know the signs of the aij but will only be able to give rough estimates for their magnitudes. This research started with the question: If the entries of A are known only approximately, can we still guarantee that 1k is a stability matrix? 1.3. Description of the PrOblem. To put the prdblem more precisely, if A is a given stability matrix and B is a given matrix, we want to know for what values of e, A + 63 is a stability matrix. More generally, let Eij denote the matrix with 1 in the (i,j) position and zeroes elsewhere. we would like to find an open region about the origin in nz-space, such that if the nz-vector (eij) is in the region, A + Z) 6 i,j would like to be able to perturb the entries of A in- ijEij is a stability matrix. That is, we dependently, not just in a prescribed ratio as is the case when the perturbation has the form CB. While most perturbation studies only deal with changes that are small compared with the original values, we will want to investigate perturbations that are of the same magnitude as the original matrix. In fact, we seek bounds on 6 that are as large as possible. In Chapter II, we find sufficient conditions for the perturbed matrices to be stable in both of the above cases, by using the criterion of Lyapunov for a stability matrix. Part of the significance of these results is that the bounds on e and on eij are not just theoretical, but methods are developed to compute them. We also find sufficient conditions for the bordered matrix [’AT b:] to c d be stable when A is stable. The methods of Chapter II depend on the solution S of the Lyapunov matrix equation (8) ATS + SA = -Q for a given symmetric, positive definite matrix Q. In Chapter III, we develOp iterative procedures to solve (8), which appear to be superior to existing procedures if the matrix A is large. In Chapter IV, we use the Routh-Hurwitz criterion for stability matrices to find computable, necessary and sufficient bounds on e for A + 63 to be a stability matrix, if B is a rank one matrix. In Chapter V, we return to an examination of the stability of the Lotka-Volterra system (4). Besides answering our original questions from the methods of Chapters II and IV, we use the criterion for a bordered matrix to be stable to find conditions for an ecosystem to remain stable when a new species is added, and use Lyapunov functions based on the solution S of equation (8) to establish a domain Of attraction for the equilibrium point pE, that is, a neighborhood of pE such that any solution starting in this neighborhood tends to pE. Finally, we consider not just constant, but also time varying perturbations B(t) and show that if the deriva- tive of B(t) is small enough, the solutions stay close to a certain time-varying critical vector c(t). If, in addition, B(t) is periodic of period T, there exists a periodic solution of period T. All the material in this thesis after section 2.1 is the original work of the author unless explicitly stated otherwise. l;4, Notation. Throughout this thesis we will use capital letters, such as A, to denote matrices and the corresponding subscripted lower case letters aij to denote the (i,j)-entry of A. Alternately A = [aij] will mean A is the matrix with aij as its (i,j) entry. All matrices are assumed to be real unless stated otherwise. We will use column vectors, which will be denoted by lower case letters, such as v, and the sub- scripted letter vi will denote the ith component of the vector v. we will also use the notation v = (vi) to define v as the column vector with ith component vi. Eigenvalues of matrices will generally be denoted by lower case Greek letters. Since real matrices can have complex eigenvalues and eigenvectors, we will not assume that any eigenvalue or eigenvector is real unless so stated. If unspecified, matrices will be assumed to be n x n and vectors to have dimension n. The transpose of a vector or a matrix will be denoted by VT or AT, and the conjugate transpose by v* (or A* should A be specified as a complex matrix). The symbol v’ will indicate the derivative of v 'with respect to time. D = Diag[d1...dn] will mean that D is a diagonal matrix With entries dii = di' dij = 0 1f 1 # 3. When _ 10 _ v is a vector, D(v) will denote the diagonal matrix Diag[v1...vn]. The identity matrix will be denoted by I. If unspecified it will be n x n. To specify a different dimension, say m x m, we will use a subscript Im. For two vectors u and v (or two matrices A and B) u << v (or A << B) will mean that each entry of v - u (of B - A) is positive. The spectral radius of a matrix A, the largest absolute value of an eigenvalue of A, will be denoted by 0(A). The symbol HvH will denote any vector norm. we will often use the vector norms defined by n (9.1) Hle - '2 )Vi) i=1 (9.2) llvllz = flv (9.3) [M], = 8111) \Vil Any matrix norm, HA”, satisfying (10) HAxll g llAllllxll. where “Ax“ and “x” are a vector norm, is said to be consistent with that vector norm. We will make use of the matrix norms ll n (11.1) HAHl = Hjfip=1 HAVHl = mgx igfi laijl' 1 _ (11.2) “All2 = sup llAvll2 = /o(ATA), (Ivnfl n (11.3) “A” = H fip ”Av“on = max Z) :aijl ° v co=1 i j=1 which are consistent with the corresponding vector norms. It is well known that for any consistent matrix norm, 0(A) g ”A . (See, for example, Isaacson and Keller [18] for a proof of this and equations (ll.1)—(11.3).) Whenever we use the symbol HA , it will stand for a consistent matrix norm. 1L5. Background. The question of the location of the eigenvalues of a matrix has often been examined, although the results are not usually stated in terms of perturbations of a stability matrix which preserve its stability. We now review some of these results and ShOW'hOW they are applicable to our problem. We assume throughout that the matrix A is a stability matrix. 1.5.1. Continuity of Eigenvalues. It is well known that the eigenvalues of a matrix are continuous functions of its entries. Thus we know that there is an open region in n2 space containing A which contains only stability matrices. 12 g:§;g, Sensitivity_Analysis. This approach to studying the effect of perturbations, based on the deri- vative of the eigenvalue, is often used in engineering applications. (See Porter and Crossley [31] or Tomovic [38]). Let M(e) be a matrix function of 6 having eigenvalue 1(6) and corresponding right and left eigen- vectors x(e) and yT(€). Differentiating (12) M(€)x(e) = x(€)x(€) with respect to e, multiplying on the left by yT(e), and solving for X'(€). yields . fiwx (13) x (e) = T Y (€)X(€) T In particular, if M = A + 8B and x = x(O), y = yT(0) are the right and left eigenvectors of A corresponding to x: x(0)l T (14) V(o) = 5%35. y X Setting B = Eij yields .X. (15) __a_).L__—_.X_1__l . 8a.. T 1] y X The eigenvalues are most sensitive to changes in the entries of A for which 'SEA— is largest. For small enough 6 iJ' the eigenvalues of A + SE can be approximated by 1(0) + ex’(0). But we do not really know x'(e) for 13 € # 0, since we do not know x(e) and y(e) for e # 0. Without this knowledge we do not have an error bound for the approximation to x(e), so that we cannot say for certain that 1(6) has not moved far enough away from 1(0) to cross the imaginary axis into the right half plane. 1.5.3. Characteristic EquatiOn. Since the eigen— values of A satisfy the characteristic equation _ n n—l _ (16) det[A—)I] — x + Cl). +...+ Cn-l" + cn — 0, it is necessary for all the coefficients ci to be posi- tive in order for all the eigenvalues to have negative real part. If all the eigenvalues of A are known to be real (for example if A is symmetric) this is both a necessary and sufficient condition for A to be a stability matrix. The Routh-Hurwitz criterion discussed in Chapter IV is a strengthening of this condition to make it necessary and sufficient for general matrices. In particular, we note that it is necessary that tr A = -c1 = sum of the eigenvalues be negative and that det A = (-l)ncn = product of the eigenvalues have the same sign as (-l)n for A to be stable. 1.5.4. Two By Two Matrices. If A is a 2 x 2 matrix it is necessary and sufficient for A to be stable that tr A be negative and det A be positive, which is useful for constructing counterexamples. Furthermore, the 14 eigenvalues of A are given by all+a 'a —a (17) 1 = ————-33-:n/4-ll§—33)2 + a 2 12a21 ° Perturbations of A which decrease the diagonal entries of A usually are stabilizing, but not always, since -3 1 changing ] is stable, but [:g i] is not. In this case from -2 to -4 decreased the average a11+322 . 2 ' . /a11‘a22 2 distance between thehn 2V' 2 ) + alza21 . of the off diagonal entries which increase alza21 (symmetric a11 value of the eigenvalues, but increased the Perturbations perturbations) move the eigenvalues further apart and hence are destabilizing. Perturbations which decrease alza21 (skew-symmetric perturbations) move the eigenvalues together. When becomes negative enough to make the quantity a12a21 under the radical zero, the eigenvalues meet and further reduction of alza21 yields complex conjugate eigenvalues. Thus skew symmetric perturbations of a 2 x 2 matrix with negative trace tend to be stabilizing or at least not destabilizing, though they may lead to imaginary components in the eigenvalues, which correspond to oscillations in the solutions of equation (3). We will see evidence that these Observations are also true for larger matrices. A precise mathematical formulation of this principle for general matrices would be valuable. 15 1.5.5. Gershgorin's Theorem and Generalizations. Gershgorin's theorem, that the eigenvalues x of a matrix R lie within the union of the circles defined by either (18) ix - riil S.§;irij! or (19) ix " rii! S.lerji| is well known. Thus if R is diagonally dominant and has negative diagonal entries, it is a stability matrix and any perturbation which preserves this characteristic pre- serves its stability. Generalizations of Gershgorin's theorem may be formulated (see Householder [17] or Wilkinson [43]) by writing a matrix R as E + C, where, for some non- singular matrix P, D = P-lEP is diagonal. (Here we drop our assumption that D and P be real.) Then D + P-ICP is similar to R, so that for any eigenvalue 1 of R, there is an eigenvector x such that (20) Dx + P-lCPX = 1x. It follows that (21) 1 g H (1)-11)"1 P'lcpll and (22) min (211 _ i) g HP’lcpll g lip-1)) up) (cu. If E is the diagonal part of R, and P = I, using the norms H H1 and H H, in (21) yields equations (18) and 16 (19) respectively. Various other choices of E and the norm used give various generalizations of Gershgorin's theorem. This can be applied to perturbations of a stability matrix A in two ways. First, if A is diagonalizable, letting E = A and C = 88 shows that A + 63 is stable for any 6 satisfying (23) eHP-lBPH g_min Reldiil. i Letting E = A and C = Z) eijEij shows that i,j A + Z: €..E.. is stable as long as i,j 1] 13 (24) H E: eijEin g ——:fL-——-min Reidii" i,j HP IIHPH i If A is normal then P is unitary and HPH2 = HP-l“2 = 1, but if A has nearly parallel eigenvectors, which generally happens when A has two almost equal eigenvalues, then the condition number HP-lHHPH may be very large. In this case,computation of D and P will be difficult and may be unstable. (Unstable is used here in the sense of numerical analysis, that small round-off errors will cause large errors in the final answer.) Second, suppose that A = E + C, that E is stable, and that (25) HP—lCPH g min Refdiif. l 17 Then A is stable and any perturbation of A which preserves (25), preserves its stability. It should be noted, however, that perturbations which change E, not only change the values of dii but also change the matrix P. 1.5.6. Symmetric and Skew Symmetric Parts. Any real matrix A can be written as (26) A=M+K where M = %(AT+A) is symmetric and K = g(A-AT) is skew symmetric. It is well known that the real parts of the eigenvalues of A lie between the smallest and largest eigenvalues of M, and that their imaginary parts lie be- tween the conjugate pair of pure imaginary eigenvalues of K of largest magnitude. Hence for M to be negative definite is a sufficient condition for A to be stable. If A is normal, the real parts of the eigenvalues of A are the eigenvalues of M and their imaginary parts are the eigenvalues of K, so that M negative definite is a necessary and sufficient condition for a normal matrix A to be stable. The stability criterion of Lyapunov discussed in Chapter II is a weakening of this condition to make it both sufficient and necessary for general matrices. The next lemma is a special case of a well known inequality for the eigenvalues of the sum of two symmetric matrices. (See Wilkinson [43].) 18 Lemma 1.10. Let M and V be symmetric matrices, u,v, and 1 be the largest eigenvalues of M,V, and M + V respectively. Then (27) ).g u + v. Proof: * * * (28) 1 = sup 2E--—(;‘I~{1-‘f\—7-—)—-‘)igsup X*MX + sup X*VX = u + v. x x x x x x x x x Corollary 1.11. Let A satisfy (26), B = V + L with V symmetric and L skew symmetric, and u and v be the largest eigenvalues of M and V, respectively. Assume u < 0. Then A + 6B is stable for all e such that u+ €v<0. Proof: The symmetric part of A + 68 is M + 8V, which by the Lemma has largest eigenvalue x(e) g u + ev. This simple Corollary is often a very practical way to establish that a perturbation preserves stability, and is the starting point for the methods of Chapter II. 1.5.7. Sign Stable Matrices. Quirk and Ruppert [32] (see also Maybee and Quirk [29p introduced the con- cept of sign stability in connection with economic problems and proved the following theorem. Definition. A matrix A is sign—stable if every matrix C such that cij is negative, zero, or positive whenever aij is negative, zero, or positive, respectively, is stable. 19 Theorem 1.12. An indecomposable matrix A is sign—stable if and only if A satisfies the conditions (i) aii g_o for every i, aii < 0 for some i. (ii) aijaji g_o for every i # j. (iii) For any sequence of 3 or more indices i,j,k,...,r (with i # j # k #...# r), the product a = 0. o la: 0.0a I 1] 3k r1 (iv) There exists a nonzero term in the expansion of det A. Of course a perturbation of a sign-stable matrix which preserves the sign pattern of the entries preserves the stability. Since the signs of the dij of equation (5) are easily determined but their magnitudes are not, we will make use of sign stability in our applications. 1.5.8. Norm of (A+1)(A—I)'1. Let C: (A+1)(A-1)‘1. Then each eigenvalue Y of C is given by _ 1.11 where 1 is an eigenvalue of A. The linear fractional transformation 1 4-%§%v maps the left half plane onto the interior of the unit circle. Therefore A is stable if and only if 0(0) < 1. Since 0(0) g_HcH for any consistent matrix norm, ”C“ < l is a sufficient condition for A to be stable. But it is difficult to apply this condition to the stability 20 of the perturbed matrix A + eB, because of the difficulty of expressing H(A+eB+I)(A+eB-I)_1H as a function of 6. 1.5.9. Canonical Forms. Stability of a matrix A can be determined by transforming it to one of the various canonical forms which display the eigenvalues. Or one can use the Schwarz canonical form [34] (see also Barnett and Storey [4]), which does not display the eigenvalues, but has all positive entries on the subdiagonal if and only if A is stable. Canonical forms are not well suited to perturbation problems, since the transformations in- volve long computational processes which obscure the re- lationship of the canonical form of A to that of A + 63. CHAPTER II APPLYING LYAPUNOV'S EQUATION gyl, The Lyapunov Equation. To study the preservation of the stability of a matrix A under perturbations we need a necessary and sufficient con- dition for A to be stable. The most useful of these is the following theorem due to Lyapunov: Theorem 2.1 (Lyapunov). If there exist symmetric positive definite matrices S and Q such that (1) ATS + SA = -Q then A is a stability matrix. Conversely if A is a stability matrix, then for any positive definite symmetric matrix 0. there is a unique positive definite symmetric matrix S which satisfies equation (1). Lyapunov's equation (1) has many uses and is discussed in many places in the literature. Bellman [5] and Barnett and Storey [4] both have detailed expositions. In this section we review the pertinent results from the literature. Section 2.2 contains original contributions to the theory of Lyapunov's equation, and sections 2.3, 21 22 2.4, and 2.5 are new applications of Lyapunov's equation, using it to study perturbations and bordering of stability matrices. Throughout this chapter, 8 and Q 'will always stand for symmetric positive definite matrices unless noted otherwise. That equation (1) is sufficient for A to be stable is closely related to the stability of the linear differential equation (2) x' = Ax. If we form the positive definite quadratic form xTSx, then (3) (xTSx)’ = xTATSx + xTASx = -xTQx, and xTSx is a positive definite Lyapunov function for (2) with a negative time derivative. By Theorem 1.9, this implies that the zero solution of (2) is asymptotically stable, which can happen only if all the eigenvalues of A have negative real part, i.e., only if A is a stability matrix. A proof that a positive definite solution S for equation (1) exists when A is a stability matrix can be found in Bellman [5], who gives the integral representation T (4) s = I" eA t Q eAt dt. 0 Olga Taussky [37] has given an algebraic proof. 23 To better understand equation (1), we examine the Operator A which maps the space of n X n matrices to the space of n x n matrices and is defined by (5) Ax = ATX + XA. It is easily seen that A is linear, and that it maps symmetric matrices to symmetric matriCes and skew— symmetric matrices to skew symmetric matrices. It is also apparent that equation (1) is linear in A, so that (6) GA + BB = dA + BB for any matrices A and B and scalars d and B. For an n X n matrix X, let XV denote the n2 x 1 column vector formed by taking X row by row. Then for any two matrices B and C, (BXC) = (BXCT) V XV' where X denotes the Kronecker product. It follows that - _ T _ T T (7) (AX)V - (A X+XA)V -— ((A x1) + (IxA ))XV. If 1 [ all a12 ... a1n 1 a21 a22 ... a2n (8) A = . ! ‘ an1 an2 °° annJ 24 then , r-AT+a I a I 11 21 anll T T 3121 A +a22I . anZI (9) (A XI) + (IXA ) = a I a I AT+a II i. 1n 2n "' nn J It is readily seen that if A is diagonal, triangular, symmetric, or skew symmetric, then (ATXI) + (IXAT) has the same property. Theorem 2.2. If a matrix A has eigenvalues 11...xn, then the eigenvalues of the linear Operator A are precisely the sums of all pairs of eigenvalues xi + xj. Proof: Let RAR_1 be upper triangular, so that the eigenvalues of A are the diagonal elements of RAR-l. Then T-1 T_1 T T T T (10) (R XR )[(A x1) + (IxA )](R XR) -1 -l = RT ATRT X I + I X RT ATRT which is lower triangular. The eigenvalues of A are -1 -l the diagonal entries of RT ATRT X I + I X RT ATRT which are the sums of pairs of diagonal entries of RAR-l. In fact if x and y are eigenvectors of AT, T so that ATx = 11x and yTA = y 12, then xyT is the 25 "eigenvector" of A corresponding to 11 + 12, since (11) A(xyT) AT(xyT) + (xyT)A = xlxyT + xyTx2 T (11+12)Xy - Corollary 2.3. If no two eigenvalues 11,12 (not necessarily distinct) of A satisfy 11 + 12 = 0, then the matrix equation ATX + XA = Y has a unique solution X for every matrix Y. Proof: The hypothesis guarantees that the linear operator A has no zero eigenvalues, hence it is non— singular. Corollary 2.4. If no two eigenvalues 11,12 (not necessarily distinct) of A satisfy 11 + 12 = 0, and Q is symmetric, then the solution S to ATS + SA = -Q is symmetric. Proof: (12) ATST + STA = -QT = -Q. By the uniqueness of solutions, ST = S. It is clear that any stability matrix satisfies the hypothesis of the preceding two corollaries, but any singular matrix or matrix with a conjugate pair of pure imaginary eigenvalues fails to do so. Li 26 243, The Lyapgnov Equations and the Symmetric Part of a Matrix. It does not seem to be widely recognized that there is a close connection between the Lyapunov equation (1) and negative definiteness of the symmetric part of a matrix. Equation (1) shows that a matrix A is stable if and only if for some positive definite symmetric matrix S, the matrix SA has negative definite symmetric part. The following corollary to Lyapunov's theorem expresses the relation between stability matrices and nega- tive definite matrices in another way. Corollary 2.5. A is a stability matrix if and only if A is similar to a matrix with negative definite symmetric part. Proof: If A is similar to a matrix with negative definite symmetric part it is clear that A is a stability matrix. If A is a stability matrix, equation (1) holds. Let S = RTR, R nonsingular. Multiplying (l) on the left -1 by RT and on the right by R-1 we have -1 -l _ -l (13) RT ATRT + RAR l = -RT QR - T"1 -1 Thus —R QR , which is negative definite, is twice the symmetric part of RAR-l. 27 If S is a symmetric positive definite matrix, we can form an inner product (x,y)S = x*Sy. We call a matrix B S—symmetric if (x,By')S = (Bx,y)s and S— skew-symmetric if (x,By)S = -(Bx,y)S. It is readily verified that V is S-symmetric if and only if SV is symmetric, and that an S-symmetric matrix has all real eigenvalues {vj} with (x,Vx) * (14) min V- = min W-S— = min w j J x ' S x x Sx (x,Vx) * (15) max v- = max —T;-§7§-= max §;§!§ j 3 x ' S x x Sx Likewise W is S—skew-symmetric if and only if SW is skew symmetric and an S-skew-symmetric matrix has pure imaginary eigenvalues {ting} with (x,Wx) * . . . SWX (16) min w. = min ————-T—-= min X* j x (X’X S x x Sx (x,Wx) * x SWX (17) max w. = max —————7—-= max * j x (X'X S x x Sx Lemma 2.6. Let S be a symmetric positive definite matrix. Any matrix A can be written as a sum of an S- symmetric matrix V and an S—skew-symmetric matrix W. If v1 and Vn are the smallest and largest eigenvalues of V, respectively, and iwl and imn, wi'g 0, are the smallest and largest (in absolute value) eigenvalues of W, respectively, then any eigenvalue 1 of A satisfies 28 (18) v1 g_Re(1) g_vn (l9) uh g]Im(xfl 3 Mn. Proof: Let V = %(S—1ATS+A) and W = %(-S—1ATS+A). Furthermore, for some x (20) 1x = Ax = Vx + wx from which it follows that * * XSVX xswx (21) X="'*——+T—- x Sx x Sx * * From the symmetry, x Sx and x SVx are real, and SW * skew-symmetric implies that x SWX is pure imaginary. Equations (18) and (19) now follow from (l4)-(l7). This leads us naturally to the following theorem, first proved in separate parts by VOgt [39], and Barnett and Storey [4] using different methods. Theorem 2.7. Let A be a matrix with eigenvalues [1k] and S a positive,definite symmetric matrix. Let ATS + SA = M and -ATS + SA = K. Then V = %S-1M ‘has real eigenvalues {v1 g_v2 g,..g_vn} and 1 . l (22) '5 min * V1 _<_ Re(1k) S Vn -- '2- max * for every k. 29 S—lK has pure imaginary eigenvalues NIH Likewise W = {3:1ij 0 g (01 _<_ (1)2 _<_...g “In and * 'k l . x Kx l x RX (23) - min —;——- = u) S_Im(x ) Sim = — max * 2 x x Sx l k n 2 x x Sx for every k. 2529;; V is the S—symmetric part of A and W is the S-skew symmetric part of A. The inequalities in (22) and (23) now follow from (18) and (19). The equalities in (22) and (23) follow from (14)—(17) since x*svx = l l * * * 5 x Mx and x SWX = 5 x Kx. When M = -Q with Q positive definite, Theorem 2.7 gives us another proof of the sufficiency of the Lyapunov condition for A to be stable. Corollary 2.8. Let ATS + SA = -I, A have eigenvalues with real parts o1 g,..g_an < 0, and Cl and on be the smallest and largest eigenvalues of S. respectively. Then 1 1 (24) OIS-ES-ESOnSHSH- Egggf; The smallest and largest eigenvalues of - % S"1 are - §%- and - 3%év respectively. By 1 n Theorem 2.7 with M = —I l l (25) —-2-5-l-_<_al_gang-20n. Equation (24) follows immediately, on g_HsH being true for any consistent norm. 30 The following stability concepts have been introduced in economic theory: Definition 2.9. A matrix A is D—stable if for any matrix D = diag[dl,d2,...,dn], DA is stable if and only if di > 0 for every i. Definition. A matrix A is S-stable if for any symmetric matrix S, SA is stable if and only if S is positive definite. Clearly S-stable a D—stable = stable. Arrow and McManus [2] have shown the following: Theorem 2.10. If AT + A is negative definite then A is S-stable. We now show that Theorem 2.10 is equivalent to Lyapunov's theorem. Suppose Lyapunov's theorem is true and AT + A is negative definite. If SA is stable, (26) (ATs)s'1 + s‘1(SA) = AT + A, and by Lyapunov's theorem S"1 is positive definite, which implies S is positive definite. Conversely if S is positive definite, so is 5‘1 and (26) is Lyapunov's equation for SA with -Q = AT + A. 31 On the other hand suppose the theorem of Arrow and McManus is true. If equation (1) holds with S and Q positive definite, then VSA is stable for any posi— tive definite V. Taking V = 8—1 shows that A is stable. Conversely, if A is stable we know that for any Q a solution S to equation (1) exists. We need only show that S is positive definite. But SA satisfies the hypothesis of Arrow and McManus. S-1(SA) = A is stable implies 8-1, and hence S, is positive definite. Corollary 2.11. Assume S and Q are symmetric, positive definite, ATS + SA = -Q, and V is symmetric and commutes with S. Then VA and AV are stable if and only if V is positive definite. Proof: Since V and S commute, VS-l is symmetric. VA = VS-lSA, which by Arrow and McManus is stable if and only if VS"l is positive definite, which happens if and only if V is positive definite. Corollary 2.12. If there is a positive definite diagonal matrix N and positive definite Q such that ATN + NA = -Q, then A is D-stable. Proof: Since any diagonal matrix D commutes with N, by the preceding corollary, DA is stable if and only if D is positive definite. 32 Definition 2.13. A matrix A is D-negative- definite if there exists a diagonal matrix D such that DA has negative definite symmetric part. We will have applications of D-stability and D- negative definiteness in Chapter V. 2;§, Stability Preserving Perturbationsjby Lyapunov's Equation. We now come to our main purpose for studying the Lyapunov equation, namely finding sufficient conditions for a perturbed matrix A + 68 to be stable when A is known to be stable. The basic idea is given by the following lemma, which we will call the Fundamental Lemma. Lemma 2.14. Let ATS + SA = -Q, with S and Q symmetric positive definite matrices. If -Q + €(BTS+SB) is negative definite, then A + EB is stable. Proof: (27) (A+€B)TS + S(A+€B) = -Q + c(BTS+SB). The right hand side of (27) is negative-definite. Thus by Lyapunov's theorem. A.+-eB is stable. The next three theorems illustrate how the fundamental lemma can be applied. 33 Theorem 2.15. Let ATS + SA = -Q, S and Q symmetric positive-definite. Let v be the largest eigenvalue of BTS + SB and q be the smallest eigen- value of Q. Then A + 63 is stable for all e > 0 such that (28) 6v - q < 0. Proof: By Lemma 1.10 the largest eigenvalue of -Q + €(BTS+SB) is less than 6V - q < 0. Thus -Q + c(BTS+SB) is negative-definite. The theorem now follows from Lemma 2.14. Theorem 2.16. Let ATS + SA = -Q, S and Q symmetric negative-definite. If 1 (29) c(IlBTH + IIBH) _<.-—-—_—-' HSIIHQ 1“ then A + 68 is stable. Proof: Let q be the smallest eigenvalue of Q and v the largest eigenvalue of BTS + SB. Then 1 _ sq- HQ 1ll Thus 6V - q is negative and A + 6B is stable by Theorem (30) ev g ellBTs + SBH g a(llBTll + NB“) IISH g 2.15. We note that if we use the 2-norm, then HBTH = “B“, and the condition (29) for A + SE to be stable becomes 34 1 (31) SHBHZ S leSHZHQ‘lllz Using the l-norm or the m—norm, HBTHl = HBHO, and condition (29) becomes 1 1 (32) ”Bill + ”EH 2 -1 = —1 °' IISII,IIQ ll, IISIIIIIQ II1 Theorem 2.17. Let S and Q be symmetric positive- definite matrices, ATS + SA = -Q, and nl.S.n2 S:~°S.nn be the eigenvalues of Q-1(BTS+SB). Then A + SE is stable: 1 l . (33a) for all --< 6 < ——, If n < 0 < nn ; n1 nn 1 1 . (33b) for all '3; < 6 < a, If fll < Uh S_O ; l . (33c) for all -a<< e <-fi;, 1f o.g n1 < nn. T T Proof: (A+€B) S + S(A+eB) = -Q + €(B S+SB), and A + 8B is stable for all e such that V(e) = —Q + €(BTS+SB) is negative definite. Since V(O) is negative definite, and the eigenvalues of V(e) are all real continuous functions of e. we may increase 6 and still have V(€) negative definite until det V(e) = 0. But (34) o = det V(e) = det(-Q+e(BTS+SB)) is equivalent to (35) o = det(% I-Q-1(BTS+SB)), 35 which implies % is an eigenvalue of Q—1(BTS+SB). The smallest positive 6 for which this occurs is % = fin unless nn'g 0, in which case V(e) is negative definite for all positive 6. Likewise V(e) remains negative definite as we decrease e from zero until % = hi, or if “1.2 0, V(e) is negative definite for all negative 6. In the preceding three theorems, it is often computationally convenient to let Q be the identity matrix. In this case the conditions of Theorems 2.15 and 2.17 for A + 68 to be stable both reduce to l l (36) -—'< e < -- Vl Vn where v1 is the smallest and Vn the largest eigenvalue of BTS + SB. The condition of Theorem 2.16 reduces to (37) c(llBTll+llBll) g 17%” . A stronger result, however, will often be obtained by using the criterion implied by the first inequality in (30): (38) eHBTS + SE” g_1. For some types of perturbation matrices B, the eigenvalues of BTS + SB are particularly easy to find. For example if B has rank one, BTS and SB have rank one and their sum has rank at most 2. Thus there are at most 2 nonzero eigenvalues of BTS + SB; they are easily found, since their sum is the trace of BTS + SB and their 36 product is the sum of all the 2 X 2 principal minors of BTS + SB. Assume that only one column of A is perturbed. which without loss of generality we may assume to be the last one, making B a rank one matrix of the form (39) B = [0.1)] where 0 denotes an n X (n—l) zero matrix and b denotes an n-dimensional column vector. We partition S into (40) s = [ MT] S ‘ . . T . where M is an n - l by n matrix and s 18 an n- dimensional row vector. Then 0 Mb (41) B S + SB = T wb M 2sTb The nonzero eigenvalues satisfy (42) v2 — 2sTbv - bTMMb = 0 or (43) v = sTb iflsTbV + bTMMb Noting that the quantity under the radical in (43) is just bTSZb, and that a similar argument can be carried out for a perturbation of other columns, we have the following theorem by applying equation (36) with 6 = l: 37 Let ATS + SA = let the jth Theorem 2.18. -1, column of A be perturbed by adding a column vector b, T and let sj be the jth row of S. The perturbed matrix is stable if (44) sgb +./bTSZb < 1. . . T T 2 Equation (43) shows that the quantity sjb + b S b will always be positive. Now let us assume that only the last row of A is perturbed so that B is the rank one matrix. 0 T (45) B = b = (bl’b , .,b ) T 2 n B ; and partition S into S = [M s] where s is the column vector (81,52....,sn)T. Now BTS + SB = bsT + st ] ZSlbl slb2+52bl slbn+snb1 _ . 52bl+slb2 2s2b2 ... $2bn+snb2 (46) —. snbl+slbn 25nbn J which is again rank two since st and bsT are both rank one. The two nonzero eigenvalues satisfy (47) v2-(ZZsb)v-Z(sb—sb)2=0 k k k k z z k k<£ 38 (48) _.\)=Z‘sb i/(Zsb)2+ 2(sb-s )2 kkk kkk k l___,_.._,____ . l] 13 2 Sji “/1? Ski Alternately, we could return to the left inequality in equation (36). Equation (50) follows from (52) and (54) since 3 . = s. for ever k. k1 1k y If ATS + SA = -Q, but Q is not the identity, Theorem 2.15 shows that we may replace the l in equations (44), (49), and (50) by q, the smallest eigenvalue of Q. 40 2.4. Convexity of S—Permissible Perturbations. Definition 2.21. Let ATS + SA be negative- definite. A matrix B is an S—permissible perturbation of A if (A+B)TS + S(A+B) is negative definite. In other words, a matrix B is an S—permissible perturbation of A if A + B is stable by the fundamental lemma, 2.14. Lemma 2.22. For a given symmetric, positive- definite matrix S, the set of matrices A such that ATS + SA is negative definite is a convex cone. Prggf: .Assume ATS + SA and BTS + SB are negative definite. Then for any a > 0, (dAT)S + S(dA) =o(ATS + SA) is negative definite. If 6,6 > o, o + B = 1, then (oA+53)Ts + S(dA+BB) = a(ATS+SA) + B(BTS+SB) is negative definite. Theorem 2.23. Let ATS + SA be negative definite. The set of S-permissible perturbations of A is an open convex set. Proof: Let B and C be S-permissible pertur- bations of A. If B,Y > O and B + Y = 1, then (55) A + BB + NC = S(A+B) + y(A+C). Convexity now follows from Lemma 2.22. 41 Since (A+B)TS + S(A+B) is negative definite, and eigenvalues are continuous functions of the entries of the matrix, there is an e such that if M is symmetric and HMH < e, (A+B)Ts + S(A+B) + M is negative definite. Let 6 = EHéHD' For any matrix R with HRH < 5 and ”RT” < 5, HRTS + SRH < 8, so that (A+B+R)TS + S(A+B+R) is negative definite. Thus the set of S— permissible perturbations is Open. This theorem is important because it allows us to answer our second question and compute open (and convex) neighborhoods of the origin which contain only stability preserving perturbations. Thus we can show that stability is preserved when the entries of A are perturbed in- dependently, and not just in a fixed ratio. Equivalently. a matrix can be shown to be stable when its entries are only known to lie in certain intervals: Corollary 2.24. Let ATS + SA = —I, and let 8.. > 0 for every i,j and Z} 8.. = 1. Then any 1] —- i,j 1] matrix R whose entries satisfy 6.. (56) a. < + $3_lw__ B . + lJ.WM . < r.. a.. 13 V/’ 2 13 13 2 sij — Ejsik s.j + Ezsik l for every i,j is stable. Proof: R can be eXpressed as 42 A+B, and from (50) B is a convex combination of S—permissible per- turbations. Example 2.25. sections 2.3 and 2.4 with the matrix r -10.0 0 O O O The solution to places is 1" .077 - -.004 .057 - -.015 -.010 - .003 .001 O -9.0 1.0 O O O T A S + .004 .088 .026 .064 .045 .031 .012 -7.0 —1.0 -2.0 O 5.0 O 0 SA = .057 .026 .265 .053 .080 .023 .005 -2.0 -5.0 O -2.0 2.0 -4.0 O A, given by O 0 0 0 -10.0 0 -8.0 -3.0 -l.0 0 0 -l.0 5.0 3.0 We illustrate the techniques of —4.0 -4.o -O.5.A rounded to 3 decimal .015 .064 .053 .232 .067 .098 .012 .010 .045 .080 .067 .474 .032 .048 .003 .031 .023 .098 .032 .295 .030 .001 .012 -.005 —.012 .048 .030 .377 43 It is not immediately obvious that S is positive definite or that A is stable, but if we take D = diag[0.2,0.2,0.2,0.2,0.5,0.2,0.4] then D + s is diagonally dominant and hence positive-definite and AT(D+S) + (D+S)A is (P -5 0 0 -0.4 -o 2 o o 0 0 —4.6 0 0 0 0 0 -O.4 O -1.8 O 0.5 O O —O.2 O O -1.8 -O.6 O 2 O 0 0 0.5 -0.6 -2.0 0 0 0 0 0 0.2 0 -1.4 0.4 L 0 O 0 0 O 0.4 —1.4_J which is diagonally dominant with negative diagonal entries and hence negative-definite. This implies that A is stable, which implies that S is positive definite. The upper and lower bounds for S-permissible perturbations bij computed from equation (50) are given by , r5.726 11.223 6.468 12.032 11.433 9.939 10.092 7.957 4.671 10.027 5.262 12.426 6.385 7.276 2.887 3.796 1.805 4.232 2.704 3.748 3.516 3.845 2.951 4.515 1.973 4.821 2.685 3.808 2.079 2.242 1.749 2.358 1.036 2.177 1.854 3.163 2.881 3.406 2.415 3.519 1.636 2.892 _2.608 2.539 2.652 2.538 2.324 2.429 1°317_, 44 (.48.64 9.812 24.623 8.920 9.281 10.571 10.404 7.697 26.853 6.596 16.313 5.852 10.556 8.780 4.304 3.172 40.321 2.921 4.787 3.206 3.398 L==- 3.460 4.758 3.052 23.672 2.927 5.662 3.490 1.994 1.864 2.434 1.791 58.621 1.911 2.257 ] 3.194 3.506 2.953 4.583 .2.873 47.147 3.490 7 2.628 2.701 2.584 2.702 2.986 2.837 216.374; Thus the (i,j) element of A may be changed from aij ij + bij for any zij < bij < uij while leav1ng the other entries fixed, and stability is preserved. For to a example, A will be stable if -216.87 g a77 g + .817 and the other entries are as given. Using the convexity of S-permissible perturbations, we know that A + B is stable for any B whose entries b.. satisfy 13 57 ..L.. b.. 3..u.., .. 0, .. = 1. ( ) 613 l] < 1] < #13 1) Bl] > 123 Bl] For example A + B will be stable for any B, l 1 . _ 2 _ l 49 L << B << 49 U. Or by taking 667 - §. B77 — 3, we see that we can add —2.327 < b67 < 1.928 to a67 and —72.124 < b77 < 0.439 to a while leaving the other 77 entries unchanged and preserve stability. 45 By comparison, HSHon = 0.757, so that by equation (37) A + B is stable if HBHl + H81!” 3 1.32. 7 = 2, b77 = 1, and all other entries of B are zero, A + 6B is stable for Equation (28) shows that if b6 -2.86 < e < .814. As another example of how the convexity of the S-permissible perturbation of a matrix can be applied, we prove the following theorem. Theorem 2.26. If M = [mij diagonal entries and the nonzero off diagonal entries ] has negative satisfy 1 ) ( (58) ‘mij' < zaij'mii for some dij with Z: dij = 1, then M is a stability i,j i741 matrix. Proof- Let A = diag[m ] and S = diag[ l ]. Since the m.. are all negative, S is positive definite 11 and ATS + SA==-I. Letting B = M - A, M ‘will be stable if B is an S—permissible perturbation. Let Eij denote the matrix with l in position m.. ~ (i,j) and zeros elsewhere. By equation (50), all Eij ij will be an S-permiSSible permutation if 1 = -2m.., 0 I 1 1 (59) 2mii="'§."."so S~s 11 ij ll 46 since all off diagonal entries of S are zero. Equation (59) is true for all nonzero off diagonal entries m. 13' m. C of M ‘by hypothesis. Thus B = Zldij(alJ-Eij) is a ij convex combination of S-permissible perturbations and hence is S-permissible itself. Example 2.27. The matrix F -1.0 0.2 O M = 0 -1.0 0.2 4.4 0 -3.0 )..- 4 . . l is seen to be stable by taking 012 — 023 — é and a _ 3 31-4 245, Bordering a Stability Matrix. Another important way to change a matrix is to border the matrix with a new row and column. Thus if we know that an n X n matrix A is stable, we seek conditions for the matrix cT d A b [. :] to be stable, where b and c are n X 1 column vectors and d is a scalar. We may treat this as a type of perturbation of A and study the stability of the new matrix by the methods of this chapter. Let ATS + SA = -Q, q be the smallest eigenvalue of Q, and r and k be arbitrary positive scalars. We will choose Optimal values for r and k later. 47 Then T A c r8 0 . r8 0 A b | + bT d c) k ; (3 k cT d (60) —rQ (er+kc) rQ 0 (er+kc)T 2kd c3 rq 0 er+kc + T (er+kc) 2kd+rq r0 (3 The smallest eigenvalue of is rq, so that C) rq the left hand side of (60) is negative definite if the 0 er+kc largest eigenvalue of is less (r8b+kc)T 2kd+rq than rq, that is if 2kd + rL+J(2kd+rq)2 + 4(rsb+ko)T(rsb+kc) (61) 2 < rq Isolating the radical, squaring both sides, and collecting terms we have (62) 2kdrq + (er+kc)T(er+kc) < 0. Since k,r, and q are all positive, this is equivalent to (63) d < _.J; (r T 2 k T T 2q k b S b + 2c Sb + E c c). The quantity in parenthesis is positive, and the least restrictive condition on d will be found by choosing k and r to minimize this quantity. We note that it is only the ratio of E which matters. The 48 function f(x) = a; + b + cx is readily seen to have a minimum at x =./ g by taking first and second deriva- . . k bTSZb . tives. Putting E = ——Er—- in (63), we have the following c c theorem: Theorem 2.28. Let ATS + SA = —Q, S and Q symmetric and positive definite, with q the smallest eigenvalue of Q. Then [AT b is stable if c d (64) d < - {-11- ((,/"cchTszb + cTSb). We note that the right hand side of (64) is always less than or equal to zero, being zero when c = —dSb, for some scalar a. Equation (64) is clearly only a sufficient condition for [AT b to be stable. Examples Lc d ‘ which are stable but have d > 0 are easy to construct. It is, however, the best result that can be r8 0 ] obtained by using a matrix of the form [0 k in equation (60). A better result might be possible using rS x [ T :1, but this makes the right hand side of (60) x k so complicated that bounds on the eigenvalues are difficult to find. If we wish to find a better (possibly positive) upper bound for d that guarantees [AT b1] to be stable, c d_) 49 could choose a (10 satisfying (64), solve the apunov equation for A b A b {:Q 0] s + s = CT do 1 1 CT d0 0 -l ‘'s o asing 1 as an initial guess in the iterative ‘3 20 0 methods discussed in the next chapter, and compute how far d 0 may be perturbed while preserving stability by equation (50) . Example 2.29. Let A be the same as in example 2.25, b? = [o,0,0,0,-1,—3,0], S the solution to ATS + SA = —I, and CT = [0,0,0,0,2,4,0]. Equation (64) becomes (1 < -(\/:’Tc \ATSZb + cTSb) - — -(4.3OXO.96-3.31) = —O.99. If we let bT = [0,0,0,0,-3,-1,0] and cT = nicho,0,+4,+2,0]; then d < -(4.47X1.47-6.34) -O.24. is a sufficient condition for the bordered A b matrix [ J to be stable. .cT d In this latter case, taking d = -0.24, computing the solution to the Lyapunov equation for the bordered matrix, and applying equation (50) shows that the bordered matrix is stable for d < 0.119. ‘1 1‘1““ CHAPTER III IFTERATIVE SOLUTION OF THE LYAPUNOV EQUATION The methods of the previous chapter depend on solvdxx; the Lyapunov matrix equation (1) ATS + SA = —Q. .As noted in Section 2.1, the Operator of A on the space n X n matrices defined by (2) As = ATS + SA is a linear operator. The major difficulty in solving equation (1) is the size of the system. Since S, Q - . 2 and AS are n X n matrices, there are n equations in n2 unknowns, and even using the symmetry of Q and S there are % n(n+1) equations in % n(n+1) unknowns. To solve such a system by Gaussian elimination would 1 l 3 l 2 2 l _ 1. 6 take §(§ n(n+1)) + (2 n(n+1)) + 3(5 n(n+1)) — EZ-n l 5 3 4 13 3 4 2 n . . . +.§ n + § n + 24 n + 3 n + 3 multiplications and divisions. (See Westlake [42], p.100.) Several authors have given methods to solve epxnfion (1) based on the characteristic polynomial of A. (See Brickley and McNamee [6], Frame [8], Jameson [l9]J But the best methods to find characteristic 50 51 equations require 11 matrix multiplications, and then these methods require another n matrix multiplications to get S after finding the characteristic polynomial. for a total of Zn4 multiplications. There are other methods (see Smith [35], and Kalman and Bertram [20]) which require first transforming A to a special form, usually a difficult and potentially an unstable procedure. A survey of several of these methods is found in Barnett and Storey [4], Chapter VI. Iterative procedures are often used to solve large systems of linear equations. One interative method to solve (l) is given by Barnett and Storey [4], but they report that it converges too slowly to be useful in many cases. This chapter adapts some of the well known iterative procedures for solving ordinary systems of linear equations to the solution of the matrix equation (1) and evaluates their effectiveness for this equation. Because of the nature of the matrices in the Lotka-VOlterra system, we will be especially interested in the effective- ness of these procedures for matrices A with large skew symmetric components . 3.1. General Iterative Procedures. Equation (1) may be written in the Operator notation as (3) ' 52 Let N be an operator which is easy to invert, and w any real number. Then (3) is equivalent to l- 1- - (4) 6 NS — Q NS + AS + Q. (5) S=S+wN-1(AS+Q). --l- -—1 (6) S = (I 2+wN A)S +‘wN Q. n The choice of the number w, called the overrelaxation parameter, will be discussed later. An alternate formu— lation is to let A = P - N, so that 1— 1. - - (4a) 6 NS = (Q - 1)NS + PS + Q (5a) S = (l-w)S + wN”1 (PS+Q) (6a) S = ((l-w)I 2+wN-1P)S + wN-lQ n We will find both forms useful. The iterative scheme derived from (5), (5b) S(1+1) = 8(1) (1) + wN—l (AS +Q) converges to the solution of (3), for any initial choice 8(0), if and only if the eigenvalues [11) of (INN-1 A) satisfy 'u' < 1. But (6) [u] = [l + w1:1 is an eigenvalue of N—IA} so that (4b) converges if and only if the eigevalues 1 of N -l A satisfy (7) :1 + W1] < l. Lemma 3.1. Let N Then the eigenvalues of {I'll-i have negative real part if and only if there is a positive real number w such that equation (5b) converges to the Solution of (l) for any initial matrix 3“”. Proof: Let [)‘k = xk+iyk} be the sczat of eigen- values of 51-117., and choose 0 < w < min 7:33-2- . This k xk+yk is possible when all the x.k are negative. For every k, . 2 2 2 (8) {1 + ka52 = l + Zka + (xk+yk)w is a convex quadratic function of w, which has value 1 -2 for w = O and w = T7. Therefore [1 + ka! < l for Xk+yk -2 -2 (9) O < w < min ——:k— -—:k—- 2 2 2 2' k Xk+yk Xk+yk On the other hand, suppose for some xk__>_o. Then l+w)\k= %k = Xk + iYk' l + ka + iyk is an eigenvalue of I 2 + WET '15. and has absolute value _>_ 1 for any n choice of w > O. 3.2. Rate of Convergence. For the convergence scheme (Sb), p(w) = C(I 2+wfi-11-x) is approximately the n factor by which the norm of the error is reduced with each iteration. If M = I 2 + wfi -15., the rate of n convergence is usually defined as !log 0 (M) '1 = -log 0 (M). (See Isaacson and Keller [18].) and A be any linear operators. Ikflmh; 54 Let {x11} be the set of eigenvalues of the operator 171-117).. To get the best rate of convergence, we want to choose w to minimize 0(I+wi\-I-l'I-\), which is equivalent to choosing w to minimize (10) p2 (w) = maxll + ka'lz = max(l+2ka+(x.fi+y}2<)w2) . k k For a fixed k, H (11) min'l+wx'2=1_ =__l‘._._ k 2 2 2 2 w Xka Xk+yk . . . ‘Xk . . and this minimum occurs for w = 2 2. Since g}, Xk+yk min p2(w) > minEl + ka!2, convergence will be slow for w w even the best choice of w if there is an eigenvalue >‘k = X'k + iyk with yk much larger than xk. If all the eigenvalues of 13-113 kn <°°'< X2< XII-<0! are real with then Figure2 1 shows that the best choice for w is where ‘1 + wxn =ll + wxllz , that isfor w:— X1+Xn° (An algebraic proof is in Isaacson and Keller [18].) For this choice of w, we have l _ >\ ..er 4 2 l (12) P(W) = (1+2). ( )+X2 ( ‘) = '. I. l xlfl‘n 1 (lehTZ' MH‘n Convergence will be slow if ”‘1! is very small compared to Mn.” that is if the condition number {7‘9- 1 171-15., is large. of I N ’“ I A: I E. . I I I I { I / / I "' I I I I I 5-0 _.’.___:2 :2. :3: :2. >‘l+)h )‘n )2 Figure 1. Graph of :l + Wkilz versus w, where kn <...< x2 < x1 < O are real eigen- values of fi”¥l§. Convergence occurs for W’ which makes :1 + wxigz < l for all i, that is for O <'w < 51% - kn w Convergence is most rapid for that which minimizes maxll + winZ, that i -2 Xl-I- x2 is for w = 56 In practice the determination of the best choice for w is complicated not only by complex eigenvalues, but also by a lack of knowledge of the exact location of the eigenvalues. 3.3. Practical Iterative Procedures. In practice, equation (5b) is used in the form (13a) P(i) = Ts(i) + S(i)A + Q (13b) NTv(i) + V(i)N = P(i) (13c) S(i+l) = 8(1) + wv(i). where N is a matrix for which it is easy to solve equation (13b) for V(l). Since we know that K = AT x I + I x AT and fi = NT x I + I x NT have eigenvalues that are sums of pairs of eigenvalues of A and N respectively, we can find criteria for con- vergence in terms of the eigenvalues or other properties of the matrices A and N. The following theorem is the most general in this direction. Theorem 3.2. A necessary condition for the matrix iterative scheme, (14) NTS(1+1) + S(l+l)N = NTS(1) + S(l)N + w(ATS(l)+S(1)A+Q) to converge for all initial values S(O), is that the vector iterative scheme (15) Nx(i+1) = Nx(i) + wa(i) + wq 57 converges for all initial values x(O) and the rate of convergence of the matrix scheme can be no faster than that of the vector scheme. Proof: In operator form, (14) is S(1+l) = (I 2+wfi'4lA)S(1) + wfi'ALQ and (15) is equiva— n X(i+l) (i) lent to (I+wN-1A)x + wN—lq. Let u be the eigenvalue of (I+wN-1A) with largest absolute value. Then (15) converges if and only if (u) < l and the rate of convergence is —log!u}. Since u is an eigenvalue of I + wN-lA, det[uI - (I+wN-1A)] = O, which implies det[(u—l)N - wA]==O, so that O is an eigenvalue of (u—l)N - wA. The eigenvalues of the Operator (Urlyfif- WA = (u—l)fi — WA are the sums of all pairs of eigenvalues of (u~l)N — wA. In particular 0 = O + 0 will be an eigenvalue. This implies det[(url)fi - wA] = O, which implies det[uI 2 — (I 2+wfi-4LA)] = 0. Thus u is an n n eigenvalue of I 2 + wI-{T -15., and p = O[I 2 + wI-{I -15.] 2 n n Eu). But (14) converges only if 1.2 p.2 I“), and the rate of convergence is -log p g_-log}u:. This completes the proof. Convergence of (15), however, is not sufficient for the convergence of (14), as seen in the following example: 58 _ —1 4 _ 1 0 Example 3.3.~ Let A - [_1 3] and N - [O 4] -1 ’1 4 I . 1 then N A = ..l 2 , which has eigenvalues - §.i i./15/8. 4 4 Since they have negative real part, by Lemma 3.1 there is a positive real number w such that (15) converges. But I I ’ 3-2 F]. -1 O ‘ q E 4 2 o -1 -. A = AT x I + I x AT = I I 4 O 2 -l I O 4 4 6 J g g 0 5 O O and N = NT x I + I x NT = i - O O 5 O o o o 8 J T 1 1 . g-l -§ --2- O : I I I 4 2 1 f __1_ I E S 0'3 1 S 5 5 g . 1 1 3 ‘ I 0 2' i ZIJ which has eigenvalues - $.i i\/l§/8, as predicted by the proof of Theorem 3.2, and a double eigenvalue 2/5. Since there are eigenvalues with both positive and negative real part, again by Lemma 3.1 there is no positive real ‘w which will make (14) converge. Even if we restrict ourselves to Operating only on symmetric matrices S, (14) will not converge. If Ts s ‘ ___ _ s = 11 12 , (N 1A(S)V is given by ale S22 59 .. , ,_ .. ‘1 '% "21' 0 S11 I f5?1 g O “E 512 g 0 g ’31" S21 0 % é 43‘: .. I{"322 -,- Since S is symmetric, we can reduce the dimensions to a 3 x 3 system by identifying s12 and $21. This requires us to add columns 2 and 3 and delete row 2 or row 3 in the matrix N-JLA, which then collapses to {-1 —1 0“ if} 2-} 5 5 5 _o 1 3.4 . . . l . 2 This matrix also has eigenvalues — §.i 1‘/15/8 and §° In other words, there are eigenvalues of the Operator N-lA with both positive and negative real part that correspond to symmetric "eigenvectors". An important consideration for using any itera- tive procedure is a good initial guess for the unknown. (0) One possibility is to take 8 = N'élQ, another is to 0) let S( be the identity matrix. 60 A convenient choice for Q is the identity matrix. If A is normal, A = M + K where M - %(AT+A) and K = %(-AT+A) and M-1 and K commute. Then (16) ATP-g W1) + (é Elna. = gnarl - KM—l + m'lm + Mflx] = -I. Thus -(AT-I-A)-'l = - M-l NIH can be used as an initial guess for S when A is close to normal and Q = 1. Our experience Shows, however, that when A is close to the form -D + K with D diagonal and K taking S(O) = % D'-l skew symmetric, is easier and works as well. AS mentioned, to actually use the iterative scheme (5b) (or equivalently (13)) we must find matrices N for which N is easy to invert, that is for which it is easy to solve (13b) for V(l). Since N = NT x I + I x NT is diagonal or triangular when N is diagonal or triangu- lar, respectively, these are the most natural choices for N. 3.4. Simple Iteration. % I is the identity Operator on n x n matrices, since % IS + S(% I) = S. Choosing N = % I gives the simple iterative scheme (17) S(1+l) = 8(1) + w(ATS(1)+S(l)A+Q). Criteria for convergence is given by the following corollary to Lemma 3.1. Corollary 3.4. Let A be a stability matrix. Then fluue is a positive real number w such that the 61 iterative scheme (17) converges to the Solution of ATS + SA = —Q, for any initial guess S(O). Proof: N = % I, N is the identity Operator, N'- A is just A, and (5b) reduces to the iterative scheme (17). But the eigenvalues of A are sums of pairs of eigenvalues of A, and hence all have negative real part. The conclusion follows from Lemma 3.1. If the ratio of the largest real part to the smallest real part of the eigenvalues of A is large, or if A has an eigenvalue with large imaginary part, the same will be true of A. In Section 3.2 we saw that the rate of convergence of (17) for such a matrix will be slow. Since each iteration of (17) takes n3 multi- plications, if it takes over 2n iterations to achieve the desired accuracy, the characteristic polynomial methods will be faster. 3.5. Jacobi Iteration. If N = diag[d1,...,dn], then equation (13b) componentwise is 13 13 3 ij p. I V. . = —]:_1—— 1] di+dj When we use the negatives of the diagonal entries of A for N in equations (13), the result is JacObi iteration ‘with overrelaxation applied to the operator A. Component— wise this can be written as 62 (2+1) (1) n (1.) n (1) (l9) sij (l-w)sij - W(kEi akiskj + kEi sik akj k#i ‘k#j + qij)/(aii+ajj) where 35;) is the value of the (i,j)-entry of S after the 2th iteration. This can be done with n3 + n2 multiplications per iteration. A sufficient condition for the Jacobi method with w = l to converge for ordinary vector equations Ax = y is that A be diagonally dominant. This is also a sufficient condition for the convergence of (19), since AT x I + I x AT is diagonally dominant whenever A is. For matrices which are exactly or nearly of the form A = -D + K ‘with D diagonal and K skew symmetric, which we expect to encounter often in a Lotka—volterra system, the following theorem is useful: Theorem 3.5. Let A = —D + K, K skew symmetric, D = diag[d1,...,dn], with di > 0 Vi. Then there is a real w > 0, such that i) --l + wD AS(1) '--l S O which makes (20) convergent. The eigenvalues of the Operator I + wD.ALA have the form 1 - w :_iwyk. Since for any fixed w, :1 — wiiwy!2 = (1+y2)w2-2w + l is strictly increasing as y increases, mixil — w i_iwyk!2 = :1 - w i_iwyn]2, and the best rate of convergence is obtained for that w which minimizes (l+yi)w2-2w + 1, that is for w = 1'2 . (See Figure 2) For this w. 031 + Wfi-élil = 1+yn //:;T"' :1 _ w 1'. iwyn! = \/ n2 , which makes the rate of convergence 1+y ____ n /y2 y2 slog \/ 2 t — 5 log 2 1+yn l+yn 64 Figure 2. _i -___.._-_.., . . , , 1 2 2 2 2 2 2 2 1+yn 1+yn 1+y2 1+yl Graph of :1 — w i_iwyn!2 versus w, where iyr1 are the eigenvalues of D'élK with O g_y1 g y2 <'°2§ Yn' Convergence occurs if )1 - w i_iwyn)2 < l for all n, that is if w < Best rate of convergence 1+yn is for w which minimizes max11 — w i_iwy )2, that is for w =-——j§. n n 1+yn W 65 Although this theorem guarantees convergence, it also shows that if the eigenvalues Of D'élK are much larger than 1 in absolute values, convergence will be too slow to be practical. 3.6. Seidel Iteration. A third possibility is to take N to be the upper (or lower) triangular part of A. Equation (13b) is then equivalent to the system of equations 1 j 21 Z} .v . + v. . = .. i = 1...n, ' = i...n, ( ) k=l aki kj kg: lkakj pij 3 which can be readily solved by taking the equations in the order i = l, j = 1...n; i = 2, j = 2...n; up to i - n, j = n, and utilizing the fact that vkj = ij' Using this choice for N in equation (5) and (13) is equivalent to Seidel Iteration applied to the operator A. Componentwise this can be written analogously to equation (19): '-1 n (“1) _ (2) _ 1 (1+1) (1.) (22) sij (l w)sij w(k=1 akiskj +k E11 akiskj 5'1 (1+1) “ u.) + Z) s. + 23 s a k1 1k akj kj+11k k3 + qij)/(a11 ajj) i: 1...n, j=i...n. Again this can be done with n3 + n2 multiplications per iteration. It is usually more efficient than the JacObi 66 method and is actually easier to program, since the old entries are replaced by the new entries as soon as they are computed, instead of requiring separate storage. Like the Jacobi method, the Seidel method with w = 1 will converge if A is diagonally dominant. The vector Seidel method is also known to converge when A is symmetric positive (negative) definite (see Fox [7], p.193), and A symmetric positive (negative) definite implies A is symmetric positive (negative) definite, so that this result also carries over to the matrix iterative procedure of (22). If A = (L + D - LT) where L is strictly lower triangular, D is diagonal, and T L + D + L is positive (negative) definite, the Seidel method can be shown to converge by mimicking Fox's proof. 3.7. Block Seidel Method. Iterative schemes from equation (5b) work best when N-J' is approximately A'él, but this is not the case for the Jacobi or Seidel schemes when A has large off diagonal entries and especially when the off diagonal part is nearly skew symmetric. In the latter case these schemes try to approximate the inverse of a matrix having eigenvalues with large imaginary components by inverting a matrix with real eigenvalues. The block Seidel method described below allows us to pull some of these large off diagonal entries into the part which is inverted, giving much more rapid convergence. 67 Let the n x n matrices A, S, and Q be partitioned into m x m blocks, [ A11 A12 "' A1r 1 rs11 S12 °°° S1r A21 A22 ... A2r 821 822 ... S2r A = S = {LArl Ar2 Arr J Sr1 Sr2 Srr (23) f‘ .1' Q11 Q12 er5 Q21 Q22 er Q = er Qr2 er where r = 3' For simplicity we assume that m divides n. It should be noted that since S and Q are symme- tric, Sii' Qii are symmetric and Sji = Sfj, jS = ng. Equation (1) written out block by block becomes r‘ T r (24) (k3; SjkAki) + RE: SikAkj = “Qij or (25) AT.s.. + S..A.. = — E) A:.Sk. - E; s. . - Q... ii ij ij 3] k=1 1 j k=1 lkAkj ij k¢i k¢j If m is small, equation (25) can be solved for the m2 entries in block Sij directly, by simply solving the 2 . m vector equation 68 T T _ (26) (Aii x I + I x Ajj)(sij)v — YV where Y is the right hand side of (25). This allows us to do a block Seidel iteration with overrelaxation completely analogous to equation (22) with single entries replaced by matrix blocks: T T _ T (1+1) r T (L) (27a) Aiivij + VijAjj - :23: Akiskj + _Z; Aki j k—i+1 +1 r + :2: Si]: )Ak.+ Sigfhxk. +Qi j +k= j+l j S(1+1) _ (z) _ (27b) Sij — (1 -w)Si ij wvij' The equations are solved in the order 1 = 1, j = l...r; i = 2, j = 2...r;...; i = r, j = r: and each time a block 8.. is found, 8.. is set equal to ST. if 13 31 13 To obtain the maximum benefit from this procedure it is best to first permute the rows and columns of ‘A to bring the largest possible Off diagonal elements into the diagonal blocks, so that instead of solving equation (1) directly one solves (28) PTATP(PTSP) + (PTSP)PTAP = —PTQP for PTSP, where P is a permutation matrix. On the computer this can be done by leaving A, S, and Q fixed but permuting the indices of the rows and columns. 69 This makes possible the Option of using a sequence of permutations to bring different off diagonal elements into the diagonal blocks on successive iterations. The number of multiplications for each iteration of the block Seidel method depends on the block size m. The right hand side of equation (27a) takes 2(3 - l) multiplications, and solving for the m2 unknowns in by Gaussian elimination takes % m6 + m4 + g m2 multiplications and divisions. Equation (27b) takes an Vij additional m2 multiplications. This must be done for (9 + 1) blocks, for a total of m S 2 l S l 3 2 n + 5 5 m + m + § + (6 m + 2 m — 3 plications. For m = 2, this is n + each of the % 3 l l 4 Malt! )n + g)n multi- 2 m 11 n + lg-n. T 2 This number could be reduced somewhat by taking advantage of the existing zeroes in the matrices Aii x I + I x A? 33" At the cost of storing an additional :§L(m2n2 + m3n) real numbers, the matrices Afi x I + I x Ajj could be stored in factored form after the first iteration, so that later iterations would require only m4 + m2 multi- plications and divisions to solve for Vij' In spite of the extra work involved, the examples Which follow Show that convergence is enough faster for the block Seidel method than for the Seidel or Jacobi method to make it the preferred method. As n gets 70 larger the number of additional operations becomes less significant, and for smaller matrices it would be better to solve them directly or use a characteristic polynomial method than an iterative method. 3.8. Examples and Comparisons. The methods described were tested on three different matrices and the results are given in the tables below. In each case the best overrelaxation parameter w was found somewhat experimentally and the result for the best w is given. It was observed that as w is increased the rate of convergence increases to a certain point, after which it rapidly decreases until the procedure becomes divergent. In the examples below D stands for the negative of the diagonal of A, w is the overrelaxation para- meter, and HR“ is the square root of the sum of squares of the entries of the residual ATS + SA + Q. In each case the procedure was run for 20 iterations unless the norm squared of the correction matrix N'41(ATS+SA+Q) -10 became less than 10 before that. Since a matrix A of Odd dimension cannot be divided evenly into 2 x 2 blocks, the block Seidel method was actually performed on the bordered matrix A O [O _1] ’ 71 Example 3.6. Simple iteration Jacobi Seidel Block * Seidel "-5 l 2 O l 1 -1 -4 7 —l O A = -—1— —2 -7 -6 0 0 10 0 1 0 -6 —4 L—l O O 4 -5 J rub.of (o) itera- Q S tions HRH T6diag[5,5,12,9.10] I .584 20 8.28xlO—2 I %’D-1 .333 20 1.12x10'2 I %D"1 .300 20 2.87x10"4 I %D-1 1.0 6 7.67x10‘8 *Two by two blocks were used. 72 Example 3.7. (0) no. of .4 Q S w iterations HR” Simple iterationI D i I 1 0.1 ” diverges I 4 I 5 * I Jacobi I I I i 0.125 20 1.32 I . I . _ D i I g 0.125 40 6.32XlO l 3 a . 1 -11 I -1 Seldel I 5D I 0.15 20 1.76XlO I , * I _ a Block Seidel I [ %I)]': 0.8 . lO 1.4lxlODS I * Two by two blocks were used and the rows and columns were first permuted to come in the order (1.5.3.4.2)- 73 Example 3.8. Simple iteration Jacobi Seidel Block * Seidel Block Seidel with sequence of 3 permuta- tions** -1O 0 —7 -2 O O O O -9 —1 -5 O O O 5 l -2 O -10 O O l 5 O -2 —8 —3 O 0 0 5 2 -1 0 —4 O O O 4 O -1 -4 0 0 0 0 5 3 -0.sj iterations -(AT+A)-l ] 0.02 60 11.9 ' T -1 ] diverges 3 -(A +A) ; 0.02 60 slowly . i 1 - s . - :5 D 1 ] diverges for all w > 1X10 4 I f; 13'1 0.4 20 6.16><10"1 I1 1 i :5 D" ] 0.3 20 2.53x10"1 7 I i q *TWO by two blocks were used and the rows and columns were first permuted to come in the order (1 3 2 4 5 7 6 8). 'k in the orders (1 3 2 4 5 7 6 8). (1 2 3 S 4 6 7 8). *The 3 permutations used put the rows and columns (1 8 2 3 4 5 6 7) and F... ‘1. . CHAPTER IV APPLYING THE ROUTH-HURWITZ STABILITY CRITERION In this chapter we use the criterion of Routh and Hurwitz to establish bounds on the size of 6 for which the perturbed matrix A + €B will be stable when A is stable. This method usually gives better bounds on the size of S than the methods of Chapter II, when only a single entry of A is perturbed, but is not applicable to all perturbations B. It also requires calculation of the characteristic polynomial with its associated difficulties, and does not yield any result like the convexity of the set of S- permissible perturbations that would allow us to com- bine the bounds for perturbations of individual entries to establish stability when more than one entry of A is perturbed. 4.1. The Stabilipy Criterion of Routh and Hurwitz. Let A be a real n by n matrix and let (1) p(X) = aoxn + boxn'l + alxn’z + b “'3 +... be the characteristic polynomial of A. The Hurwitz matrix H is defined to be the n x n square matrix 74 75 [b0 b1 b2 b3 I a0 a1 a2 a3 0 b0 b1 b2 . H = 0 a0 a1 a2 0 O b0 b1 0 0 a0 a1 . ”L ] ’15 where each row is terminated with as many zeroes as necessary to make H n x n. ii The following theorem, first proved by Hurwitz in 1895, can be found in Gantmacher [9]. Theorem 4.1 (Hurwitz, 1895). Assume aO > 0. All the roots of the polynomial p(X) (or equivalently all the eigenvalues of the matrix A) have negative real part if and only if all the leading principal minors of the matrix H are positive. Let us reduce the matrix H to upper triangular form by elementary row operations as follows: Multiply a rows 1.3.5,... by - ~59- and add them to rows 2,4,6,-°°. 0 Call the entry in position (2,2) c0, multiply rows b 2,4,6,... by — 39- and add them to rows 3,5,7,---. 0 Call the entry in position (3,3) d multiply rows 0! 76 c 3,5,7 by - 59- and add them to rows 4,6,8,°°°. Con- O tinue in this manner until the matrix is upper triangular. We note that none of these Operations change the determinant of H or any of its leading principal minors. The first, second, third,..., up to the nth principal minor of H are called the Hurwitz numbers D1,D2,D3 ... Dn' The numbers bO'CO'dO on the diagonal of the triangularized matrix are called the Routh numbers. Clearly D1 = b0, D2 = bOcO' D3 = bOcOdO, etc. This shows that the following criterion for stability, first given by Routh in 1877, is equivalent to Hurwitz' criterion. Theorem 4.2 (Routh, 1877). Assume aO > 0. All the roots of the polynomial p(x) (or equivalently all the eigenvalues of the matrix A) have negative real part if and only if the Routh numbers b d a0: O’CO' Op 0.. are all positive. ‘An independent proof of this theorem can also be found in Gantmacher, as well as a proof of the following theorem due to Orlando: Theorem 4.3 (Orlando, 1911). Let x1,x2,...,xn be the roots of p(l) . Then the (n—l)th Hurwitz number is given by (2) Dn-l = (-1) 2 a““1 n (xi+xk). 77 and the nth Hurwitz number Dn = det H is given by n(n+1) _ 2 n lgi O with det H(€O) = O, and A + SE is stable for all positive 6. The left hand inequality in (4) is established similarly by decreasing e from 0 until det H(€) is zero. This completes the proof. The characteristic polynomial of A + 6B, det(lI—(A+€B)), will in general be a polynomial in both A and 6. Unfortunately it will not be linear in 6, as the hypothesis Of Theorem 4.5 requires, except in special cases. J.S. Frame, in private communication, showed that if B has rank 1, the characteristic poly— nomial will be linear in a. The next lemma generalizes this result to perturbations B of arbitrary rank. Lemma 4.6. The degree in e of the character- istic polynomial of A + SE is less than or equal to the rank of B. Proof: Let m be the rank of B. Then there exist n x m matrices X and Y of rank m, such that B = XYT. Since _ _ _ F (7) XIn A X = XIn A €XY X In 0 SYT I _ 0 I €YT I 80 xIn—A x ' (8) det = det(lIn-A—€B). €YT I n1- Since 6 Ioccurs in only m rows of the matrix on the left in (7) and (8), the highest power of e that can In occurIWhen the determinant is expanded is e . The degree in e of the characteristic poly- nomial may be strictly less than the rank of B, as can be seen by taking A upper triangular and B strictly upper triangular. Suppose that B has all zero entries for one row or column, which for ease of notation we take to be the last column. Then B has rank 1, +eb a11 a12 "° a1n 1n a21 a22 ... a2n+€b2n (9) A + 68 = . an1 an2 .. ann+€bnn J and 81 I )X-all ‘a12 "' ”ain a21 X’a22 "' ’azn det(XI-(A+€B)) = det .anl an2 x—ann_‘ a21 x'azz "' 'b2n + E det ' Lan1 an2 ... -bnn j Thus the characteristic polynomial can easily be written as p(l) + €q(l) and the Hurwitz matrix as Hp + €Hq. 4.3. Examples. Example 4.7. We illustrate this technique by computing bounds on 8 such that the matrix —5 —10 0 e] I8 -1 -4 oi (ll) A(€) = ' is stxflale. A = A(O) is stable, since it is sign stable. The characteristic polynomial of A is 4 3 (12) p(x) = I + 91 + 113;.2 + 3191 + 390 and tine characteristic polynomial of A(e) is 4 3 (13) p().) + €q()\) = I. + 9). + 113).2 + 319). + 390 + e(-16). 82 The Hurwitz matrix for A(e) is ?F9 319 I I 1 113 390 H + 8H = I P q 9 319 L 1 113 390.1 (14) , 0 0 0 ] + 6 O O O -16 [O O O O O We form leHq by elementary row operations on the augmented matrix. r'9 319 0 0 I 0 <3 '0 <3 I 1 113 390 0 ; 0 0 ~16 0 1 0 9 319 0 ; 0 0 0 0 I_0 1 113 390 ' 0 0 0 —16.J First we reduce the left half to upper triangular form with the Routh numbers on the diagonal: r‘9 319 0 0 I 0 c) 0 0 I 0 77.56 390 0 I 0 0 -16 0 0 0 273.74 0 f 0 0 +1.86 0 [0 0 0 390 I 0 0 -.527 —l6_] 83 From this we see that the (n-l)th Hurwitz number for A(e) is Dn-l= 9 X 77.56 x (273.74+€ 1.86). Orlando's theorem shows that A(e) has a pair of conjugate pure imaginary eigenvalues when 6 = 33% = -l47.l7. The determinant of A(e) is (390—168), showing that A(€) 390 has a zero eigenvalue when 6 = T6_ = 24.375. Completing the computation we find [0 0 8.54 0 O O —.241 O —1 (15) H H = p q 0 O 147A17 O 0 0 -.0013 -—4iL—— I.. 24.375 4 . . l l which has nonzero eigenvalues W and - 224—3—75- Thus A(€) is indeed stable for (16) -147.17 < e < 24.375. Example 4.8. It will not always be possible to compute Dn-l as a linear function of e as in the previous example. For example, if we perturb the (2,4) element of the same matrix A, r —5 -10 0 0 I (17) A(e) = (F 84 the characteristic polynomial is 4 3 1 + 91 + 11312 + 3191 + 390 (18) 9(1) + €Q(l) + €(-2 x-lO) . and (19) H + 6H = + e P q 0 9 319 I 0 0 -2 0 I. l 113 390‘ O 0 0 -IO‘] Dn-l is now a more complicated function of 6, but [.0 -4.92><10"l 4.57 0 I -1 O 7.61x10-3 -1.44><10—l O (20) H H = _5 _3 p q ' 0 9.43x10 3.07x10 0 [ 0 1.88x10‘5 1.18x10‘3 —1.994 which has eigenvalues -1.99,O,1.OleO—3, and 9.67XlO-3. A(e) is stable for _ -1 _:l____ (21) —103 _ 3 < e < -1.99 — 0.502. 9.67x10" When we can apply Theorem 4.5 we get the true bounds on e: for A + SE to be stable, since A + eB has zero or pure imaginary eigenvalues when 6: is equal to the upper or lower bound. Our method based on the Lyapunov equation gives us sufficient bounds on e for A + SE to be stable, but may be much more restrictive than necessary. This point is readily seen in the follow— ing simple example. 85 -l O] 10 —1 We know that Example 4.9. Let A = [ -l e . . . l [10 _1] Will be stable if and only if e < Tau The solution to ATS + SA = -2I, is 51 5 5 I—! Equation (SO) of Chapter II says that we have stability for (23) -.0432 = 2 < e < 2 = .0356. 5 -\/512+52 5 + 5124.52 much more restrictive than actually necessary. But 2 0 1 0 0 (24) H +€H =[' }+6[ P ‘1 ..1 1 ,- ..0 —103! —l O O . and H H = [ ] has eigenvalues O, -10, and p q 0 —lO equation (4) tells us that we have stability for (25) -co<€<—]:" lO ' which are the true bounds. The eigenvalues U1 and UZ of leHq in (4) are usually the easiest eigenvalues to compute numerically. If the characteristic polynomial is easy to calculate and the perturbation B has rank one (in particular if only one entry of the matrix is being perturbed), the method Of this chapter is preferable to that of Chapter II. But if we wish to perturb entries in several different 86 rows and columns, especially if we wish to perturb them independently, the Lyapunov method is more suitable. r..- _‘ EU” 3 CHAPTER V APPLICATIONS TO LOTKA-VOLTERRA AND OTHER MODELS OF ECOSYSTEMS 5.1. The Community Matrix. We now return to our original problem, to investigate the stability of an ecosystem whose population levels are at equili- brium. We use the Lotka-Volterra model introduced in Chapter I, and investigate the stability of an equili- brium point in the first quadrant for the system of differential equations (1) pi = pi(ri + n 3 Hid Q m H n H s 3 Changing to matrix and vector notation (1) becomes (2) p’=D(p)[r+Gp]: ]. where, as before, D(p) is the matrix diag[pl.p2.---.pn We will call the matrix a = [Qij] the Volterra Matrix. The process of linearizing (2) about an equilibrium point as done in the introduction, may be written in this notation as follows: Find pE such that (3) r + GpE = 0. We assume that pE >> O. Letting (4) x=p-pE. We obtain 87 88 (5) x’ = D(p)dx = D(pE)ax + D(x)ax. Levins [23] has appropriately named the matrix _ E _ E (6) A — D(p )6: - [pialj]' the community matrix. Mathematicians will recognize it as the JacObian of (1) evaluated at pE. Since Ax is the linear part of (5) and (7) D(x)6.'x = o(x), the zero solution of (5) and hence the equilibrium point pE of (2) is asymptotically stable by Theorem 1.6 if A is a stability matrix. In this chapter we investi- gate the implications of our studies about stability preserving perturbations of A (and the insights we have gained about the nature of stability matrices) for stability of ecosystems. First we note some of the strengths and weaknesses of our approach to the problem of stability. The major weakness is that we can only prove local stability, not global stability. That is, we know that if the pertur— bation x(O) of the initial populations away from the equilibrium point is small enough, the pOpulations will return to their equilibrium values, but we do not know that any initial pOpulation will move toward the equili- brium values. we will later find some sufficient conditions for x(O) to be small enough. On the other hand we must remember that local stability is a necessary condition for global stability. 89 The major strength of our work is its applicability not only to Lotka-Volterra models, but also to other models of ecosystems, and in fact to general autonomous dynamical systems. For if (8) p’ = HP) and (9) f(pE) = 0. then letting x = p - pE we have (10) X' = Ax + g(x) ‘where A is the Jacobian matrix given by (1].) a.. = and g(x) = O(x). May [28] has reviewed some of the forms commonly used for the function f. We could also study age dependent and sex dependent effects by treating age and sex classes as different species. Gilpin and Justice [11] have shown that an arbitrary ecosystem model which they write in the form (12) pi = piRi(p) i = 1...n, can be approximated near an equilibrium point pE such that (13) R(pE) = 0 by a Lotka-Volterra model, which they (and many other authors) write in the form introduced by Gause [lO], r. n 2 _ __Zl_._ _ (15) Vii = 1. Although they give a detailed geometrical description of their procedure, what they have really done is to approximate R(p) by a first order Taylor series about E p . The relationships between (1), (12), and (14) are given by a . (16) a.. = —El-(pE) = a .Y.. i l] apj il 13 n 6R n E (17) ri — - Z) 3——-(pE)pE = Z? 913p j=1 93' 3 3=1 3 r. l _— (18) k ‘ aii 1 But if one wishes to study the stability of the equilibrium point pE by approximating (12) by a Lotka-Volterra model (1) and then forming the community matrix (6), he might as well take the Jacobian (11) directly, since if (19) fl (p) = piRi (p)! using (16) and (13) gives 5f. 8R. 1 E _ E i E = E - - (20) 5Pj (p ) — pi SE;'(p ) piaij for 1 ¢ 3 . 91 M. BR. 1 E _ E E 1 E _ E 5P1 (p) — Ri(p ) + pi Tpi (p ) — piaii . (21) and the Jacobian matrix calculated from (11) is the same as the community matrix calculated from (6). The constants ri, ki’ and Qij in (1) and (14) are usually given the following biological interpretations: ri is the intrinsic rate of growth of species i, i.e. births minus deaths per individual per unit of time when the population is near zero; ki is the carrying-capacity of the environment for species 1, i.e. the number of individuals of species i which causes the rate of growth to decline to zero due to overcrowding even in the absence of other species; and dij is the change in the growth rate of species i per each individual of species j, so that jg: aijpj is the total change (either damping or acceleration) of the growth rate due to species interaction. It is often better to measure population levels in terms of units of biomass (that is, the total mass of all members of a species) instead of numbers Of individuals. When this is done, the word individual is replaced by the phrase unit of biomass throughout the above paragraph. In equation (12) the intrinsic growth rate for species i would be Ri(O), and the carrying capacity 92 would be the value of pi which makes Ri(O...O pi O...O)==O. As Gilpin and Justice point out, these are not necessarily the same values as ri and ki obtained when approxi- mating (12) by (14). This is because, although (22) Ri(p) m ri + §Tdijpj near the equilibrium point pE, it may not be a good approximation away from pE, so that it may not be true that (23) Ri(0) ksri (24) R. (O...O k. O...O) :8 r. + O...k. = O. i i l 11 i Gilpin and Justice therefore state that when a general model of the form (12) is approximated by a Lotka-Volterra model the coefficients no longer have any biological significance, but what we must do is interpret them in terms of what happens near the equilibrium point. Thus we could reinterpret ri as the intrinsic growth rate at equilibrium, which is just balanced by the dampingk n (or acceleration) Z) d..p. when p = pE, and __.= j=l 13 3 ri as the number Of additional individuals the ii environment would support if the levels of the other species are changed so as to decrease the damping by one unit. We return to the community matrix A given by (6) or (11) and investigate the biological significance 93 of the entries of A. Whether the matrix A comes from a Lotka-Volterra model or a more complicated function f, the signs of the entries aij and aji correspond to the types of species interaction: (—+) species j preys on species i, (or is a parasite of species i), (--) both species compete for the same resource, (++) symbiosis, (00) no interaction, (+O) commensalism, and (—O) amensalism. Because of inefficiencies in energy use and transfer, the gains to a predator in a predator prey interaction are never as great as the loss to the prey. Thus if j preys on i, we can assume |a..[ > Ia..'. ij ' 31‘ The diagonal entries aii indicate the amount of self interaction: aii negative means that growth of species i is self damped near the equilibruim by some type of inner competition for resources, aii zero means that growth of species 1 is only limited through interactions with other species. Positive aii in the Lotka-Volterra model is biologically unrealistic, since it would imply positive Qii which would mean that in the absence of other species, species i grows without bounds and even faster than exponentially. The trophic level of a species is the number of steps in the food chain between it and the abiotic nutrients. Thus green plants are in the first trophic level, herbivores in the second, carnivores which eat 94 herbivores in the third, etc. (Some species, however, may be in more than one trophic level.) It is Often assumed that at the first trophic level the species exhibit self damping since plant growth is limited by abiotic factors, but that there is no self damping at higher trophic levels, Since herbivores and carnivores are limited only by the existence of prey and by being preyed upon. Holling [16], however, has found that many predators do exhibit self damping be— havior if their numbers become too large. The existence of some self damping is crucial for stability, because A cannot have all negative eigenvalues unless it has a negative trace. Of course if the aii are sufficiently negative, i.e. if there is enough self damping, then A is almost certainly stable. We have also seen that large symmetric off diagonal entries (competition and symbiosis), tend to be destabi- lizing, and that large skew symmetric entries (predator- prey) tend to have a neutral effect or else to be stabilizing in the sense that they pull the real parts of the eigenvalues closer together. This leads us to the conclusion that even though the system may be un- stable at one trOphic level. coupling it through predator- prey interactions with a stable trophic level may bring 95 about stability. (May [28] has calculated the conditions for this to happen in a 3 species two trophic level system.) The matrix A in Example 2.25 (pp.42 and 73) example of a possible community matrix for an ecosystem with four trophic levels, 2 plant species, 2 herbivores, 2 carnivores, and one top carnivore. This example shows no direct competition within any trophic level. It would be reasonable to add competition terms between the plants by making both a12 and a21 negative. The matrix A(O) in Example 4.7 (p.81) is another example of a community matrix, representing a simple food chain with only one species at each trophic level. Alternatively, it might be considered as the result of lumping all the species together in each trOphic level. Our original goal was to take such a community matrix, with the entries known only within certain error bounds, and devise methods to decide whether any matrix within these error bounds gives a stable system. We can now do this using the methods of Chapters II and IV, if we are willing to do enough numerical computation. Unfortunately, the author found no readily available examples where the community matrix for a natural ecosystem has been estimated. Most of the empirical work in ecology has been at the micro-ecology level, studies of 96 interactions between a few competing species and a few laboratory studies of predator-prey relations. In spite of declarations by many of the leaders of the discipline about the importance of studying ecosystems as complete systems and the dangers of looking at only isolated parts (see for example Holling [15] and Hardin [14]), empirical work on a macroecology scale is still in the pioneering stage. We will let examples 2.25, 2.29, 4.7, and 4.8 serve as illustrations Of what the mathematics developed in this thesis makes possible once the coeffi- cients of a community matrix have been estimated. Some recent theoretical work has been done which attempts to compute the interaction coefficients in a strictly competitive community (no predator—prey interactions) from the a measurement of "niche overlap." (See McArthur [26] and [27].) More work is clearly needed on methods to measure the strength of predator-prey in— teractions and their effect on birth and death rates. It should be remembered that the community matrix A is obtained by multiplying each row of the original matrix a = [dij] by the corresponding equilibrium values p?, which in turn depend on the growth rates ri. We Observe in nature that due to energy losses as one moves up the food chain, the equilibrium values for lower trophic levels, measured in units of biomass, are 97 generally greater than those of higher trOphic levels. A careful study of conditions for a Lotka-Volterra system to have a set of positive equilibrium values [pg], and the effect of these values on the community matrix needs to be undertaken. The equilibrium value factors in the community matrix complicate the effect that perturbations of the coefficients in (l) have on the community matrix. If growth rate vector r is perturbed to r + Ar and the Volterra matrix a is perturbed to a + afi, then the new equilibrium values pE + ApE satisfy E E (25) r + Ar + (6+Afl)(p +Ap ) = O. and the new community matrix A + B satisfies E E ,, _ E E (26) A + B = D(p +I_\.p Ham...) — D(p Mama) + D(Ap )(6HM). Thus (27) APE = -(6+M)'I(Ar+MpE) and the perturbation to the community matrix is (28) B = D(pE) m + D(ApE) (mm. 5.2. D-Stable Volterra Matrices. If the Volterra matrix a is D-stable (Def. 2.9), then the community matrix A will be stable for any vector of growth rates r which produces positive equilibrium 98 ‘values. Of course negative equilibrium pOpulations are unrealistic. In some systems the Volterra matrix a is clearly D—stable. If we consider only a trophic level (or simple food chain) model, the matrix a will have the sign pattern 'which will be sign—stable and thus D-stable if any of the diagonal terms are strictly negative. Since the first level can be assumed to exhibit self damping, all < O, and the community matrix A is stable. Another special form of the matrix a, used by Kerner [21], leads to D-stability. The basic assumption is that when species j preys on species i, the gain to species j is Y.d.., where —d.. is the loss to J l] 1] species i, and that the proportionality constant Yj depends only on the efficiency of species j in utili— zing its food, and not on Species i. This leads to the conclusion that equation (1) can be written in the form ’ (29) P1 " pi(ri+)‘i 2311345” (100 = -a0 I ' 31 13 99 for which the community matrix A has (i,j) entry _ E (30) aij — pilidij. Under the assumption that all the diagonal terms Qii are zero, Kerner shows that n E 1 (31) 0 = _ZI (pi(t)-pi 1n pi(t)%;7 i=1 1 is a constant of motion for this system and uses it to prove that there are neutrally stable oscillations about the equilibrium point and to build a "statistical mechanics" for the system. Of course if the diagonal entries are all negative, (32) D = diag 5' pixi gives (33) ATD + DA = diag[dii], the community matrix is stable, and the system has oscillations spiraling in toward the equilibrium point. These results are often criticized as being "fragile" since they are based on the probably invalid assumption that a.. = -d... Neither Kerner's school nor his critics 31 13 have realized that they can get the same conclusions with the following weaker hypothesis. If a is D-negative-definite (Def. 2.13), which will be the case if a has the form (34) a = B(M+K) I" 'u_ .: I; 100 'with E positive diagonal, M symmetric, negative definite, and K skew symmetric, then a will be D- stable and any positive equilibrium point will be asymptotically stable. (Kerner's formulation (29) assumes that M is either 0 or a negative diagonal matrix.) Thus we can tolerate competition terms and departures from Skew-symmetric in the predator-prey terms which are not too large compared to the self damping terms, and still have asymptotic stability for any realistic set of growth rates. As the symmetric terms become larger, we may still have asymptotic stability for particular equilibrium values, which means for particular growth rates {r1}, but not for all of them. Finally as the competition terms become too large, the eigenvalues are pushed far enough apart that some of them become positive, the system becomes unstable and some species become extinct. This is expressed in ecology as the competitive exclusion principle: that two Species competing for exactly the same resources cannot co—exist indefinitely. §;33 Addition of a New Species to an Ecosystem. Suppose we have an ecosystem modeled by equation (1) with stable community matrix given by equation (6), and add another species. We seek conditions for the new en- larged system to be stable. The new community matrix lOl ‘fiill be the old community matrix bordered by a new row and column plus a perturbation of the old community matrix due to the change in the equilibrium values. The enlarged system can be written as (35) pf = pi(r' + i = l...n 1 1 ). jpj + bipn+l 1W 0'. . 1 p p n+1 = pn+l(rn+l + j ijj + dpn+l)° W Thus the Volterra matrix a is replaced by the bordered a b matrix [- T :1] where b is the vector (b _ ...bn)T »C d 1 and cT is the vector (cl...cn). Let pE + ApE be the vector of new equilibrium values for species l,...,n and pg+l the equilibrium value for the added species. and let r be the vector (r1....,rn)T. Then pE+ApE a b _l [r E = - CT d (r -pn+1 n+1 (36a) C _ _ _ _ 1 fl czl+lanJk]' -la]b rr e e l T -l l L. -E C a a rn+l where h (36b) e = d — cTa’lb. . E -1 . . Since p = -d r, (36) implies r +CTpE E _ l _ T -l _ _ n+1 (37a) pn+1 — e(rn C a r) _ +1 d-cTa-lb 102 _ l T -l — _ E _ (37b) Ap _ a(r +l-c a r)« 1b _ _pn+la lb. Our first requirement for a stable ecosystem is that the new equilibrium values be positive, i.e. T E r -c p E n+1 (38a) p = - _ > O “*1 d—cTa 1b E E_ E —]1 (38b) p + Ap — p — pn+ld >> O The new community matrix is now seen to be (39) D E T E \ - E E' 1 . E - E E — p +Ap )[4 bJJ= ,A-pn+lD(a lbm D(p -pn+la lb)b pn+1 E ch d d pn+1 L pn+1 If S is the solution to the Lyapunov equation T (40) A S + SA = -I, we see from equations (28) and (64) of Chapter II that sufficient conditions for stability are (41) Vp§+l — l < 0, v the largest eigenvalue of -[aTD(a"lb)s + SD (CI—me'] pE (42a) d < Mngl V‘/;T chES2 Eb + c TSEb), l_\)‘pn+1 where (42b) E = D(pE—pfi+la‘lb). 103 Equations (38b) and (41) lead to the surprising conclusion that we are most likely to have stability in the expanded system, if the vector b, which gives the effects of the added species on the old species, is directly proportional to the vector r, which gives the intrinsic growth rates of the old species. For if b = fir with B > 0, then E _ E E which shows that the inequality (38b) is satisfied, i and the perturbation to A is (44) -p§+lD(a—lb)a = Bpfi+lA: which shows that v = -B in (41) and (42). This may be hard to accomplish in nature. since the usual assumption that the plants have a positive growth rate would imply that the added species should be preyed upon by the plants! The difficulty of adding a new species to an ecosystem and still maintaining stability, supports the theme stressed by Robert May, that stability is not en- hanced by increased complexity. In fact, our results are another addition to the list of mathematical studies indicating that increasing the number of species makes it more difficult to have stability. (See May [28].) ‘ 104 Perhaps the reason that many field ecologists hold the Opposite View, that complexity increases stability, is that they consider stability to be the lack of large fluctuations or oscillations over a short time period, 'which is something very different from Lyapunov asymptotic stability, which allows large oscillations as long as they are damped over a long time period. 5.4. Domain of Attraction by Lyapunov Functions When pE is asymptotically stable, solutions which start close enough to pE approach pE in the limit. We now want to investigate how close is close enough. Definition 5.1. The domain of attraction of a point pE for the differential equation (8) is the set of all points pO such that the solution p(t) with initial value pO satisfies lim p(t) = pE. t4ao Expositions of the following theorem, due to LaSalle [22], can be found in Hahn [12], and Yoshizawa [44]. Theorem 5.2 (LaSalle, 1960). Given the differential equation (45) y’ = F(y). let V(y) be a function with first order partial derivatives, V(y) 2 O, and V’(y(t)) g_o in a region R = (y: V(y) < a}, and let M be the largest invariant 105 subset of the set {y:V'(y) = 0}. Then every solution of (45) which starts in R approaches M as t 4 a. Corollary 5.3. Assume that in equation (45), F(O) = 0. Assume that V(y) has first order derivatives, V(y) > O and V'(y) < O for y # O on the set R = {y:V(y) < a}, and V(O) = V'(O) = 0. Then 0 is an asymptotically stable equilibrium point for (45) and R is contained in the domain of attraction of 0. We apply this corollary to get an estimate of the domain of attraction for the Lotka-Volterra system. The matrix norm used in the following theorem can be n ”1. n HZ. or n n, Theorem 5.4. Let pE be the equilibrium point of equation (1) given by equation (3). Let A defined by equation (6) be a stability matrix, (46) B = GD(pE) = D'1(pE>AD(pE). Q a positive definite matrix with smallest eigenvalue q, and S the solution to (47) BTS + SB = -Q with smallest eigenvalue 0. Then pE is asymptotically stable and any pOpulation vector p such that the vector y defined by ____i._ (48) Yi ' E 106 satisfies (49) ’ T ‘1/5 q S < . W Y (HBTH+HBH)HsH is in the domain of attraction of pE. Proof: For any solution p(t) of equation (1), let x(t) be given by equation (4), and y(t) by (48). Then (50) y(t) = D’1(pE)x(t) and from (5) (51) y’(t) = D’1(pE>x’>any 2 a ‘ > max 2 +iB , and from Corollary - 1' — k “k k 1 1 2.8 [Is]! 2 - E'd" and o _ .. 33—. Thus n l 1 (llBTH+HBH) ”SH 2-“- S, - ‘Vé1a1)3 m:x!ak+inL/21all ! Hence if B is ill-conditioned (7521' large) or has ' 1 large imaginary parts to its eigenvalues the inequality in (59) will be very restrictive. If, however, a is D-negative-definite, the next theorem shows that the entire positive quadrant is in the domain of attraction of pE. The Lyapunov function used was first applied to Lotka-Volterra systems by Aiken and Lapidus [1], who were following the school of thought of Kerner and of Montroll and Goel [30]. They applied it, however, only when a was a diagonal plus skew symmetric matrix, failing to realize once again that it gives the same results with much weaker hypothesis. 108 Theorem 5.5. Let the matrix a = [Qij] in equation (1) and (2) be such that D4 has negative definite (negative semidefinite) symmetric part for some matrix D = diag[dl...dn], di > O for all i. Then for any growth rates r such that the equili- brium vector pE defined by equation (3) satisfies (56) pE >> 0, pE is asymptotically stable (stable), and its domain of attraction includes all pOpulations p with (57) p >> 0. Proof: Let p(t) be any solution of (l) with initial value pO satisfying (57). Define p. (58) V(p) = Eldi(pi-p?-p§ ln(;%»- 1 Since (59) z — l — ln(z) _>_ 0 with equality only if z 1, V(p) 2.0 for any p >> 0 with equality only if p = pE. Now , . E p1,”) (60) V (p(t)) = :Ldi(pi(t)-pi m) and p; = Xi. so that by (5)(or equation (6) of Chapter I) (61) V’(t) = Zldi((p§+xi) Zldi.x. — p? Zlaijxj) 1 J J J J . _ ~ . _ T (62) V (t) — §’§;Xidiaijxj — X Dflx. 109 Stability when DC has negative semidefinite symmetric part now follows from Lyapunov's Theorem 1.9. Asymptotic stability and the domain of attraction when D6 has negative definite symmetric part follow from Corollary 5.3 to LaSalle's theorem since any point P >> O is in the set [p:V(p) < a} for some a. We note that the paper by Aiken and Lapidus contains an error. They claim asymptotic stability when a is a diagonal plus a skew symmetric matrix with non-positive diagonal entries but only one nonzero diagonal entry. But that only makes the symmetric part of a negative-semidefinite, proving stability but not asymptotic stability. We also note that this generalizes a result of McArthur [27], who arrived at the same conclusion when Dd was symmetric and negative—definite for some posi- tive diagonal matrix D, by using -XTDax as a Lyapunov function. But this limits him to a restricted set of competition equations, since predator-prey or unbalanced competition interactions introduce a skew symmetric component to a. 5.5. Time Varying Coefficients in the Lotka- Volterra Model. Certainly it is more realistic to expect the growth rates ri and even the interaction coefficients dij to vary with time instead of being 110 constant. To investigate how this affects stability, we add time varying perturbations to the coefficients in the Lotka-Volterra model. Equation (1) becomes n (63) pi = pi(ri+pi(t) + iii (Gij+Bij(t))Pj)- Let p(t) = (pi(t)) and B(t) = [Bij(t)]. There is no longer a fixed equilibrium point as in equation (1), but we define a time varying critical vector c(t) by (64) (r+p (t)) + (d+z3(t))c(t) = 0. (We may think of c(t) as a moving equilibrium point.) Let (65) z(t) = p(t) — c(t). so that n (66) z{(t) = Ci(t)[iéi(aij+flij(t))zi(t)1 n + zi(t)[iEi(aij+Bij(t))zi(t)J + c{(t) or in matrix notation (67) z’(t) = D(c(t))[61+0(t)]z(t) + D(z(t))[a+/3(t)]z(t) + c’(t), (68) c’(t) = —[a+a(t)1‘1[p'(t)-5'(t)c(t)]. we will show that if c'(t) is small enough, p(t) will stay close to the critical vector c(t), or in other 111 words that solutions z(t) of (67) will be bounded, the size of the bound depending on the size of c’(t) and the eigenvalues of a + B(t). Equation (67) contains a linear term, a quadratic term, and a forcing term c’(t). We express the linear part as (69) z’(t) = [A+B(t)]z(t) where (70) A + B(t) = D(c(t) ) (a+a(t)). Equation (70) does not define A and B(t) uniquely. We could define them uniquely by defining A, say by equation (6), and letting (70) define B(t), but it is not clear what is the best choice for A. We require A to be a stability matrix, so that we can solve the Lyapunov equation (40) for a positive—definite symmetric matrix S. Then if B(t) satisfies the con- ditions of Theorems 2.15, 2.16, or 2.17 for all t, zTSz will be a Lyapunov function for (69), the linear part of (67) will be asymptotically stable by Theorem 1.9, and this choice for S in equation (71) below will make q positive, as will be required in Theorems 5.7 and 5.8. Ideally we would like to find the positive—definite symmetric matrix S which makes q in (71) as large as possible, but it is not clear how to do so. 112 Lemma 5.6. Let S be a positive definite symmetric matrix, (71) -q = sup[1argest eigenvalue of [ATS+SA+BT(t)S+SB(t)]], t (72) d = sup 2HsH2HG+B(t)H2. t (73) L = sup 2HSc’(t)H§, t and o the smallest eigenvalue of S. Then 3 1 T 2 T (74) (zT(t)SZ(t))’ g.— % zTSz + d(§7§§) + L(57§£) Proof: From equations 67 to 71, (75) (zT(t)SZ(t))’ = zT(ATS+SA+BT(t)S+SB(t))z + 22T(G+B(t))TD(z)Sz + ZZTSc'. (76) zT(ATS+SA+BT(t)S+SB(t))z g_—quz. By the Cauchy Schwarz inequality (77) zzT(a+B(t))TD(z)SZ g_2HzH2H(a+B(t))TD(z)sZH2 _<_ 211mm IIZIISIIZHMz) ((2)2): 1 _<_ 211mm) lizilsilz(sz>2 . since HD(z)H2 g HZHZ' Also 1 (78) ZZTSc' g 2HzH2HSc’ll2 S ZHSC'H2(sz)2. Since sz g_% zTSz, (74) follows. 113 Theorem 5.7. Let q,d,L, and o be as in the previous lemma. Assume that q is positive, 2 9_. (79) L < 4d and let (80) w = V/q2 — 4dL .. then for any solution z(t) of (67) with l 2 T /6 . . 2(0) 6 R = {z:(z 82) < ‘35-(q+w)}, z(t) remains in R for all t. and for any 6 > O the values of l (zT(t)Sz(t))2 are in the set [0,‘ég-(q-w) + e] for sufficiently large t. 1 Proof: Let W(t) = (zT(t)Sz(t))§. From (74) a 1 d 2 (81) W s 56-57; W O -9 L o W + 0172)' From the comparison principle (see Yoshizawa [44]). W(t) g_u(t) for all t, where u(t) is the solution to I _1 d 2 q L _ (82) u (t) — 5(3'37'2-11 - '6 u +3175). 11(0) — W(O)- I But from (79) and the quadratic formula, u is negative for all u between ‘ég (q :.w) and positive elsewhere. , Thus u(t) 4‘24;- (q-w) for all 0 _<_ u(O) g‘é-g- (q+W), which completes the proof. 114 Theorem 5.7 shows that if z(O) = p(O) - c(O) is in R, z(t) is bounded, and hence if c(t) is bounded, p(t) = c(t) + z(t) is bounded. It also shows that if c’(t) is small enough in comparison to q, the solution vector p(t) follows the critical vector c(t) quite closely, because EggE-(q-w) goes to zero as L = sip 2H8c’(t)H goes to zero. Thus as c’(t) decreases the difference between p(t) and c(t) also decreases. Figures 3, 4, and 5 illustrate these points. The two dimensional system (83) pi pl[lO + (-2.5+cos wt)p + (—5 — 2 sin um)p2] 1 I p2 = p2[—5 + 4pl + (—l+cos wt)p2] was solved numerically by fifth order Runge-Kutta and plotted along with the critical values cl(t) and c2(t) defined by (64), for three different values of w. As w is decreased q and d remain the same, but L decreases, resulting in smaller differences between pi(t) and Ci(t)' One can also observe the initial transient solution corresponding to the time period when p(t) - c(t), as measured by (zT(t)Sz(t))1/2, is de- creasing, followed by the steady state oscillations during the time period when zT(t)Sz(t) is in the set [0, q-w) + e]. 0 2d ( Population Levels 115 ' 0.0 018 1T6 , 7.4 31.2 410 Time Figure 3. Solution and critical values of equation (83) with w = 2N. Prey = p1(t) -"‘ Predator = p2(t)-———— c1(t)~-—- c2(t)— — 116 Population Levels Time Eigu£§_4. Solution and critical values of equation (83) with u)=7T. Prey = p1(t)o --- Predator = p2(t)-———— cl(t)‘—° C2(t)-—_ 117 0 I N50 o )' 0 3.2 I\ (I) ‘. H‘ E\ : .\ (E!) \D a “. \ 2 .4~ .3 - '.\ p .. C: y a. .3 '\ 4.) (U r-I :5 Q-( 0 04 Time Figure 5. Solution and critical values of equation (83) with u) = E' Prey = pl(t)o --' Predator = p2(t) ———— cl(t) '- ' c2(t) --— 118 It would certainly be reasonable to model many natural ecosystems with periodic coefficients in equation (63). If the fluctuations are small enough, we can prove the existence of a periodic solution: Theorem 5.8. Let q,d,L,o and w be as in Theorem 5.7. Assume that q is positive and that equation (79) holds. If pi(t) and Bij(t) in equation (63) are periodic functions of period T for all i and j, then there exists a periodic solution to (63) with p(O) — c(O) <2 O(q—w). Proof: Let F:Rn 4 Rn be defined by F(zo) = z(T), where z(t) is the solution to (67) with initial value z(O) = 20. From standard theorems of differential equations, F is well defined and continuous. Let H = [20 :.¢/zgSzO _IgC-(q-W)) H is closed and convex. Further F maps H 4 H, since 1 (zT(T)Sz(T))2 g_u(T), where u is the solution of /" (82) with u(O) = 20 3‘53- (q—w), which implies u(t) S91: (q—w) for all t. Schauder's fixed point theorem (see Hale [13]) implies that F has a fixed point, i.e. z(T) = 2(0) for some solution 2 of (67) starting in H. From (64), c(t) is periodic with period T and hence (67) is periodic with period T. From the uniqueness of 119 solutions with given initial values, it follows that E(t+T) = 2(t) for all t, and therefore 5(t) = 2(t) + c(t) is periodic with period T, which completes the proof. Putting Theorems 5.7 and 5.8 together, we see that there is not only a periodic solution p(t) to (63), but that we can expect it to be close to the critical vector c(t) if c’(t) is small. Finally, for any other solution p(t), let y(t) = p(t) - p(t), so that O pj‘:(t) T "' (84) Y1 ’ 51(t) Yi + 3 pi(t)(“ij+gij(t))yj 7‘ . .. .. + .j. y (a +£31J(t))yJ Let _ Ei’(t) (85) A(t) = D(p(t))[a + B(t)] + Diag _ . pi (t) I then if there exists a positive definite symmetric matrix S, such that AT(t)S + SA(t) is negative definite for all t, p(t) is asymptotically stable. CONCLUSIONS This thesis develops methods to show that a stability matrix remains stable under perturbations. The problem was motivated by the Lotka-VCIterra models of ecosystems, and the necessity for such a study to make these models useful in practice was pointed out, but it is clear that the topic is important for any large dynamical system whose parameters can only be estimated. This will usually be the case for models of biological systems. The results will also be im— portant for the field of structural stability, since most of this thesis could be described as a study of the structural stability of linear systems under linear perturbations. The summary of pertinent results from matrix theory in the intro— duction should provide starting points for further research on the subject. Sensitivity methods based on derivatives of eigenvalues (Section 1.5.2), generalizations of Gershgorin's theorem (Section 1.5.5), and the transformation A.--o(A+I)(A—I)—1 (Section 1.5.8), appear the most promising to yield more in— formation about the effect of perturbations on stability. 120 121 This thesis has exploited the two classical necessary and sufficient criteria for stability, those of Lyapunov and of Routh and Hurwitz. Of these, the Lyapunov criterion is the best suited to perturbation analysis. Theorems 2.15-2.19 and especially Corollary 2.20 give us sufficient bounds on the perturbations that are computationally feasible, as are the conditions for preserving stability when bordering a matrix in Theorem 2.28. The convexity of the S-permissible perturbations, Theorem 2.23, is especially useful, since it provides a method to establish Open regions in n2 space that contain only stability matrices. This is what we really need, not just stability under a single perturbation, if we are to establish the stability of a matrix whose entries are only known to lie within certain intervals. Open regions about a matrix A containing only stability matrices can also be established using Gershgorin type theorems (Section 1.5.5). In practice these two methods should complement each other, since the Gershgorin methods are computationally easiest and give the best results when the eigenvalues of A are spread apart, whereas the Lyapunov equation methods are computationally easiest and give the best results when the eigenvalues of A are close together. 122 Improvements in the Lyapunov equation methods might be obtained by investigating the best choice for Q in the equation (1) of Chapter II. Further investi- gation of criteria for BTS + SB — Q in equation (27) to be negative-definite would also be useful, especially when B has rank two or higher. . The Lyapunov equation has many other applications besides perturbation analysis. In fact, we have seen some examples in Chapter V. Other examples can be found in Barnett and Storey [3] and Chapter 5 of May [28]. Thus the improved understanding of the relationship of this equation to the symmetric part of A (Section 2.2), and the iterative procedures to solve it (Chapter III) are valuable in themselves. The block Seidel method appears to be a strong candidate for the best procedure to solve the Lyapunov equation for large matrices (say dimension 10 or higher). While block Seidel and other iterative procedures for linear equations are well known, to the author's knowledge they have not previously been adapted to the solution of the Lyapunov equation. Theorem 3.2, which tells us how well we can expect these pro- cedures to perform on the Lyapunov equation, is important not only for the iteration schemes studied in this thesis but also for any which might be studied in the future. 123 While the Routh-Hurwitz criterion is often used to prove stability of a particular matrix, Theorem 4.5. employing this criterion and Orlando's theorem to give necessary and sufficient upper and lower bounds for a rank one perturbation to preserve stability of a matrix, is new. Further research in this direction should in— clude extensions of this result to perturbation matrices of higher rank. Perturbations of at least rank 2 will be important for ecosystem analysis, since changes in the strength of an interaction between two species will generally change both aij and aji' After developing all this mathematical machinery, it is disappointing to not find any examples in the ecological literature to apply it to. Certainly an attempt to estimate the coefficients in the Lotka-Volterra model for a natural ecosystem would be valuable. Perhaps the results of this thesis, making it possible to evaluate the stability of the model with only estimates of the coefficients, will make such a study more attractive to the ecologists. This thesis has contributed to a general under- standing of the types of species interactions that contribute to stability. More specific contributions to mathematical ecology include the conditions for 124 maintaining the stability of the system when adding a new species, and estimating the domain of attraction for an asymptotically stable equilibrium point of a Lotka—Volterra model. We also weakened the hypothesis for the Lyapunov function of Aiken and Lapidus and thereby generalized results both of Aiken and Lapidus and of McArthur, showing that the domain of attraction is the entire first quadrant when the community matrix is D—negative definite. This leads us to hope that the general theorem on domains of attraction (Theorem 5.4) could be improved, and points out the need for further investigation of D—negative definiteness. We have also seen that D—stability is important for ecosystem models, and conjecture that the two may be equivalent. When are all the equilibrium values positive, and whether ecosystems are "trophic-level—stable", are other questions that merit additional research. The final section, where the coefficients of the general Lotka-Volterra model are made time varying, is the first investigation of this problem. We saw that the solutions follow a moving critical vector and the closeness of the solutions to the critical vector depends on the derivative of the critical vector, and that if the coefficients are periodic there is a periodic solution. 125 The existence of a periodic solution when the time variations of the coefficients are small enough follows from a general theorem of differential equations (see Hale [13]), but we have shown an explicit bound on the size of the variations that is small enough. Weakening the hypothesis of Theorems 5.7 and 5.8 and investigating the uniqueness and stability of the periodic solutions should be goals of additional research in this area. But this will probably require non—linear techniques, whereas the methods developed in this thesis are essentially based on linearization. There are still many unanswered questions about the differential equations prOposed fifty years ago by Lotka and Volterra as a model for ecosystems. There has recently been a great increase in interest and re- search about them, but most of it of a theoretical rather than an applied nature. This thesis has helped to answer some of the theoretical questions about stability, and added others to the list of unanswered ones. But its major contribution is to improve the usefulness of the Lotka-Volterra equations for analyzing the stability of actual ecosystems. BIBLIOGRAPHY 10. 11. 12. BIBLIOGRAPHY Aiken, R., and Lapidus, L. 1973, The stability of interacting populations, Int. J. Systems Sci. 4, 691—695. Arrow, K.J., and McManus, M. 1958, A note on dynamic stability, Econometrica 26, 448—454. Barnett, S., and Storey, C. 1968, Some applications of the Liapunov matrix equation, J. Inst. Math. App. 4, 33-42. Barnett, S., and Storey, C. 1970, Matrix methods in stability theory, New York: Barnes and Noble. Bellman, R. 1970, Introduction to Matrix Analysis, New York, McGraw—Hill. Brickley, W.G., and MacNamee, J. 1960, Matrix and other direct methods for the solution of systems of linear difference equations, Phil. Trans. Roy. Soc. London, Series A, 252, 69—131. Fox, L. 1964, An introduction to numerical linear algebra, Oxford: Clarendon Press. Frame, J.S. 1963, The rectangular matrix equation AY + B = YC, Notices Am. Math. Soc. 10, 566. Gantmacher, F.R. 1959, Applications of the theory of Matrices, New York: Interscience. Gause, G.F. 1934, The struggle for existence, Baltimore: Williams and Wilkins. Gilpin, M.E., and Justice, K. 1973, A note on non— linear competition models, Math. Biosciences, 17, 57-63. Hahn, W. 1967, Stability of motion, New York: Springer Verlag. 126 13. 14. 15. 16. 17. 18. 19. 20. 21. 22. 23. 24. 25. 127 Hale, J.K. 1969, Ordinary differential equations, New York: Interscience. Hardin, G. 1969, Not peace, but ecology. In Diversity and stability in ecological systems, Brookhaven Symposium in Biology, No. 22, Springfield VA: Nat. Bureau of Standard, U.S. Dept. of Commerce, p.151. Holling, C.S. 1969, Stability in ecological systems. In Diversity and stabilipy in ecological systems, 0p. cit., p.128. Holling, C.S. 1966, The functional response of invertebrate predators to prey density, Mem. Entomol. Soc. Canada, 48. Householder, A. 1953, Principles of numerical analysis, New York: McGraw—Hill. Isaacson, E., and Keller, H. 1966, Analysis of numerical methods, New York: John Wiley and Sons. Jameson, A. 1968, Solution of the equation AX + XB = C by inversion of an M x M or N x N matrix, SIAM J. Appl. Math. 16, 1020-1023. Kalman, R.E., and Bertram, J.E. 1960, Control system analysis and design via the "second method" of Liapunov, trans. A.S.M.E.J. Basic Engpg., 82D, 371—400. Kerner, E.H. 1957, A statistical mechanics of interacting biological species, Bull. Math. Biophys., 19, 121-146. LaSalle, J.P. 1960, The extent of asymptotic stability, Proc. Nat. Acad. Sci. U.S.A., 46, 363-365. Levins, R. 1968, Evolution in a changing environment, Princeton: Princeton University Press. Lotka, A. 1925, Elements ofyphysical biology, Baltimore: Williams and Wilkins. (Reissued as Elements of mathematical biology, Dover, 1956.) Lyapunov, A. 1907, Probleme general de la stabilité’ du mouvement, Annales de la Faculte’des science de Toulouse, second series, 9. (Reissued in Annals Math. Studies, 19, Princeton: Princeton University Press, 1947.) 26. 27. 28. 29. 30. 31. 32. 33. 34. 35. 36. 37. 38. 39. 128 MacArthur, R.H. 1969, Species packing, or what competition minimizes, Proc. Nat. Acad. Sci., 64. 1369-1375. MacArthur, R.H. 1970, Species packing and competitive equilibrium for many species, Theor. POp. Biol., 1, 1-11. May, R.M. 1973, Stability_and complexity in model ecosystems, Princeton: Princeton University Press. Maybee, J., and Quirk, J. 1969, Qualitative problems in matrix theory, SIAM Review, 11, 30-51. Montroll, E., and Goel, N. 1971, On the Volterra and other nonlinear models of interacting populations Rev. of Modern Phys., 43, No. 2. Porter, B., and Crossley, R. 1972, Model control. theory and applications, New York: Barnes and Noble. Quirk, J., and Ruppert, R. 1965, Qualitative economics and stability of equilibrium, Rev. Economic Studies, 32, 311-326. Rosen, R. 1970, Dynamicalygystem theory_in biology, New York: Interscience. Schwarz, H.R. 1956, Ein Verfahren zur Stabilitatsfrage bei Matrizen-Eigenverte-probleme, Z. Angew. Math. Phys., 7, 473-500. Smith, R.A. 1966, Matrix calculations for Liapunov quadratic forms, J. Diff. Eqps., 2, 208—217. Strobeck, C. 1973, N species competition, Ecology, 54, 650-654. Taussky, O. 1961, A remark on a theorem of Lyapunov, J. Math. Anal. Appl., 2, 105-107. Tomovic, R. 1972, General sensitivity theory, New York: American Elsevier. Vogt, W.G. 1965, Transient response from the Liapunov stability equation, Preprints Joint Automatic Control Conf., Paper V4, 23-30. 40. 41. 42. 43. 44. 129 Volterra, V. 1926, Variazioni e fluttuazioni del numero d'individua en specie animali conviventi. (Translation by M.E. Wells in R.N. Chapman's Animal Ecology, New York: McGraw—Hill, 1931, pp.409—448.) Volterra, V. 1931, Leqons sur la theorie mathematique de la lutte pour la vie, Paris: Gauthier-Villars. Westlake, J. 1968, A handbook of numerical matrix inversion and solution of linear eggations, New York: John Wiley and Sons. Wilkinson, J.H. 1965, The algebraic eigenvalue problem, Oxford: Clarendon Press. Yoshizawa, T. 1966, Stability theory by Liapunov's second method, Tokyo: The Mathematical Society of Japan.