Instability in linear cooperative systems of ordinary differential equations

It is well known that, contrary to the autonomous case, the stability/instability of solutions of nonautonomous linear ordinary differential equations $x' = A(t) x$ is in no relation to the sign of the real parts of the eigenvalues of $A(t)$. In particular, the real parts of all eigenvalues can be negative and bounded away from zero, nonetheless there is a solution of magnitude growing to infinity. In this paper we present a method of constructing examples of such systems when the matrices $A(t)$ have positive off-diagonal entries (strongly cooperative systems). We illustrate those examples both with interactive animations and analytically. The paper is written in such a way that it can be accessible to students with diverse mathematical backgrounds/skills.


Introduction
It is a well-known fact that for an autonomous system of linear ordinary differential equations (ODEs) where A is a constant n by n matrix with real entries, the zero solution is asymptotically stable if and only if the real parts of the eigenvalues of A are negative.
Unfortunately, for nonautonomous systems of linear ODEs there is no hope for a similar result. Indeed, one can find examples of systems (1.1) such that for all t all the eigenvalues of A(t) are negative but there is a solution of (1.1) whose norm tends to infinity as t → ∞. For a nice paper on that subject, see [2]. The authors of [2] start by taking a 2 by 2 matrix B having negative (real) eigenvalues with eigendirections close to each other. In such a case, although each solution of the autonomous system x ′ = Bx converges to zero exponentially, there are solutions whose magnitude increases during certain time intervals (indeed, it suffices to notice that there are x ∈ R 2 such that Bx · x > 0). Then they construct an autonomous system x ′ = A(t)x by rotating the system x ′ = Bx around the origin at such angular velocity that some solution is being kept, for sufficient amount of time, in the set where A(t)x · x > 0.
In biological applications there are additional conditions imposed on the structure of differential equations. For instance, in some bacterial populations there is switching between two states (dormant vs. active). It should be remarked here that, to be sure, "real life" models are usually nonlinear, but a linear system, x ′ = A(t)x, can serve as a first approximation. If we let x1 stand for the density of bacteria in the dormant stage and x2 for the density of bacteria in the active stage, a12(t) (resp. a21(t)) describes the transition rate from the active to the dormant state (resp. from the dormant into the active state) at time t. It is straightforward that a12(t) and a21(t) must be nonnegative for each t. See, for example, [4].
A celebrated theorem, usually attributed to Erich Kamke [3], but obtained earlier by Max Müller [5], when applied to linear systems having matrices with nonnegative off-diagonal elements states that if at an initial time all coordinates of the initial value are nonnegative then they remain so at later times. Another celebrated theorem, due to Oskar Perron and Georg Frobenius, establishes that at each t there is a real eigenvalue of A(t) corresponding to an eigenvector having all coordinates nonnegative. It is a direct consequence of those theorems that in dimension two the construction given in [2] cannot be extended to such systems: there, eigenvectors of the "frozen" matrices A(t) rotate at constant angular velocity around the origin, whereas here at each t one eigenvector of A(t) must lie in the first quadrant.
It appears that, perhaps besides one example given as an aside remark in 4.4 on p. 583 of [2], the problem of finding linear nonautonomous (in particular, periodic) ODE systems whose matrices have nonnegative off-diagonal entries with all eigenvalues negative and bounded away from zero but which allow solutions divergent to infinity has not been addressed.
In the present paper a method is given of constructing such twodimensional systems. The idea is the following: during the first half of period, the evolution of the system is governed by a (far from symmetric) constant matrix having its eigenvector with both coordinates positive (the principal eigenvector ) close to one coordinate axis, whereas during the second half of period, the evolution of the system is governed by another (far from symmetric) constant matrix having its principal eigenvector close to the other coordinate axis; around the half-period there occurs a very fast (in Section 3 instantaneous) change in matrices.
The paper is written in a reasonably self-contained way. It is assumed that the reader knows the standard facts from linear ordinary differential equations (transition matrices, etc.). However, the knowledge of timedependent (or even time-periodic) linear ordinary differential equations is not indispensable (except Section 4): As the systems considered are piecewise constant, knowing basic properties of the matrix exponentials should suffice.
As the material related to matrices with positive off-diagonal entries does not usually form a part of the curriculum, needed proofs are given by using only the knowledge in calculus and elementary algebra (Subsection 2.4).
In the "heuristics" parts (Subsections 2.3 and 2.6, and Subsubsection 3.2.1) much stress is laid on how the direction of the solution of an autonomous system of two linear differential equations changes in time. In particular, it is shown that its evolution is governed by some nonlinear ODE. All this serves as an illustration whether the model chosen by us should exhibit the required properties. Those parts are, however, independent of the main body of the paper.
After those geometric illustrations a "hard" analytical proof is given, requiring, however, only the knowledge of the form of the matrix exponential.
All sections are endowed with remarks putting the material in the perspective of more advanced studies. Section 4 requires more advanced knowledge (however, not reaching beyond the Gronwall inequality or matrix norms). It can serve as a basis for some undergraduate homework.

Systems of linear ODEs
Consider a system of two linear ODEs where we assume that A : J → R 2×2 is a continuous matrix function (J ⊂ R is an interval not reducing to a singleton, and R 2×2 denotes the set of real 2 by 2 matrices).
It is a standard result in the course in ODEs that for each t0 ∈ J and each x0 ∈ R 2 there exists a unique solution, x(·; t0, x0), of the initial value problem and that solution is defined on the whole of J.
Usually stress is laid on fundamental matrices: X(·) is a fundamental matrix solution of (2.1) if its columns form a basis of the vector space of solutions of (2.1). For our purposes, however, it is better to use the transition matrix , that is, a matrix function of two variables, X = X(t; s), s, t ∈ J, such that for any s ∈ J and any x0 ∈ R 2 , there holds If X(·) is a fundamental matrix solution, the transition matrix is given by the formula The transition matrix is unique. We mention here important properties of the transition matrix: 1. X(s; s) = I, for any s ∈ J, where I is the identity matrix; 2. X(u; s) = X(u; t)X(t; s), for any s, t, u ∈ J; 3. X −1 (t; s) = X(s; t), for any s, t ∈ J.

Systems of autonomous linear ODEs. The matrix e tA
In modern courses in ODEs, when considering systems of autonomous linear ordinary differential equations usually a matrix function t → e tA is introduced, where Occasionally, for typographical reasons we write exp(tA) instead of e tA . It is proved that the above series has convergence radius infinity, the function t → e tA is differentiable, and the relations • e 0·A = I, • e (s+t)A = e sA e tA , s, t ∈ R, hold. Consequently, the solution of the initial-value problem for a system of ordinary differential equations with time-independent matrix A, In other words, the transition matrix X(·; ·) for the system (2.2) is given by the formula X(t; s) = e (t−s)A , s, t ∈ R.
We will use in the sequel the following fact.
If AB = BA then e t(A+B) = e tA e tB = e tB e tA for all t ∈ R. (2.4) However, for general A, B the above equalities need not hold.

2.3
The action of e tA on the unit circle In this subsection we shall analyze how the radiuses and directions of solutions of the systemẋ = Ax change in time. In other words, we investigate the action of e tA on vectors in R 2 .
We start by introducing some notation.
Recall that we can represent x ∈ R 2 in polar coordinates, x = x [cos θ sin θ] ⊤ , where x = √ x · x and θ is the polar angle of x. However, as there is no easy formula for θ in terms of the Cartesian coordinates of x, we prefer to write any nonzero vector x ∈ R 2 in the form We denote by S the set of all vectors y ∈ R 2 with unit radius (norm). In other words, S is the unit circle.
Let x(t) be a nontrivial (that is, not equal constantly to zero) solution of x ′ = Ax. That is, x(t) = e tA x0 for some nonzero x0.

How does e tA act on the lengths of vectors?
As a warm-up we try to find an ordinary differential equation satisfied by x(t) . After some calculus we obtain (2.5)

How does e tA act on the directions of vectors?
The present subsubsection can be skipped, since it will be needed later only for heuristic considerations in Subsubsection 3.2.1. Let us find an ordinary differential equation that is satisfied by the direction of x(t). We differentiate d dt or, after putting y(t) : We can say that y(t) is a solution of a system of two (nonlinear) ordinary differential equations, written in the matrix form as Observe that for any y ∈ S the vector A − (Ay · y)I y is perpendicular to y. It follows that for a solution y(t) of (2.6) we have from which we can conclude that, if at an initial moment t0 the value y(t0) is equal to one then it is equal to one at any time. So, although system (2.6) is well defined for all y ∈ R 2 , we will consider it for y belonging to S only.
If y ∈ S is not an eigenvector of A then the nonzero vector A − (Ay · y)I y, perpendicular to y, points either clockwise or counterclockwise.
We will return later, in Subsection 2.6, to analyzing the action of e tA .

Matrices with positive off-diagonal entries -Their spectral properties
We write R 2 + for the set of all those x = [x1 x2] ⊤ such that x1 ≥ 0 and x2 ≥ 0, and R 2 ++ for the set of all those x = [x1 x2] ⊤ such that x1 > 0 and x2 > 0.
Let M stand for the family of real 2 × 2 matrices having their offdiagonal entries positive, and let P stand for the family of real 2 × 2 matrices having all entries positive.
Our first result is usually known as the Frobenius-Perron theorem. As we are in dimension two, we will give here an elementary proof of it.
Then the following holds: (i) A has two real eigenvalues (denoted λ2 < λ1).
(ii) An eigenvector u corresponding to λ1 can be taken to have its coordinates positive.
(iv) An eigenvector v corresponding to λ2 has its coordinates (nonzero and ) of opposite signs.
Proof. (i) The characteristic polynomial of A has the form Consequently A has two real eigenvalues, λ2 < λ1.
so, since u1 and u2 have the same sign, we must have λ1 > a11 and λ1 > a22.
The other inequality follows from the first one by the fact that λ1 + λ2 = a11 + a22.
The larger eigenvalue, λ1, of A ∈ M will be called the principal eigenvalue of A (sometimes the terms dominant, leading, or Perron eigenvalue are used). An eigenvector u of A pertaining to the principal eigenvalue will be called a principal eigenvector of A. When speaking of a principal eigenvector we always assume that both its coordinates are positive.
A principal eigenvector of length one is called normalized. A normalized principal eigenvector of a matrix in M is unique.

Strongly cooperative systems of linear ODEs
We say that a system of linear ODEs is strongly cooperative if for each t ∈ J the matrix A(t) belongs to M. In this subsection we assume that the matrix function A(·) is continuous. Now we will give the two-dimensional version of the Müller-Kamke theorem. It is formulated in the linear setting, however a closer inspection shows that its proof is rather nonlinear in the spirit.
Theorem 2.1. Assume that system (2.7) is strongly cooperative. Then X(t; s) ∈ P for any s < t, s, t ∈ I.
Proof. Fix an initial moment s ∈ J. We start by noting that the first column of the matrix X(t; s) is the value at time t of the solution [x1(t) x2(t)] ⊤ of system x ′ = A(t)x satisfying the initial condition x(s) = [1 0] ⊤ . It follows from the uniqueness of the initial value problem for linear systems of ordinary differential equations that for any t ∈ J both x1(t) and x2(t) cannot be simultaneously equal to zero.
We have x ′ 2 (s) = a21(s)x1(s)+a22(s)x2(s) > 0, consequently x2(t) > 0 for t sufficiently close to s, t > s, t ∈ J. By continuity, since x1(s) > 0, x1(t) > 0 for t sufficiently close to s, t ∈ J. At any rate, there exists τ > s such that x1(t) > 0 and x2(t) > 0 for all t ∈ (s, τ ). We claim that τ = sup J, that is, x1(t) > 0 and x2(t) > 0 for all t ∈ J, t > s. Indeed, suppose to the contrary that this is not so, that is, there exists θ > s such that x1(θ) ≤ 0 or x2(θ) ≤ 0. As the product of the functions x1(t) and x2(t) is continuous, it follows from the Intermediate Value Theorem that the set { t > s : x1(t)x2(t) = 0 } is nonempty. Specialize τ to be the greatest lower bound of this set, and assume, for definiteness, that x1(τ ) = 0 (therefore x2(τ ) > 0). τ cannot be equal to s, because we have already shown that x1(t) > 0 directly to the right of s. Consequently τ > s, so x1(t) > 0 for t < τ , t sufficiently close to τ , from which it follows that We have thus shown that the first column of the matrix X(t; s), has, for all t > s, positive entries. By applying a similar reasoning to the solution of system x ′ = A(t)x satisfying the initial condition x(s) = [0 1] ⊤ we show that the second column of the matrix X(t; s), has, for all t > s, positive entries, too.
The full strength of Theorem 2.1 will be needed only in Section 4. In the main part, Section 3, we have matrices independent of time. In such a case we can give an alternative proof, using the theory of matrix exponentials only.
Then e tA ∈ P for all t > 0.
Proof. Assume first that A ∈ P. Then t k A k ∈ P for all t > 0 and k ∈ N, consequently e tA ∈ P.
If A belongs only to M but not to P, we putÃ := aI + A, where a := 1 − min{a11, a22}. ThenÃ ∈ P and, by the previous paragraph, e tÃ ∈ P. As (aI)A = A(aI), there holds e tÃ = e at e tA (see Eq. (2.4)), from which it follows immediately that e tA ∈ P.
Remark. One could be tempted to use the approach applied in the proof of Theorem 2.2 in proving Theorem 2.1. But this is not so: the obstacle is that X(t; s) for system x ′ = A(t)x need not be equal to exp( t s A(τ ) dτ ).

2.6
The action of e tA on the unit circle, continued Again, the present subsection can be skipped, because it will be helpful only in heuristic considerations why we have chosen such an example.
We assume that A ∈ M. As e (t−s)A is the transition matrix of the system x ′ = Ax, Theorem 2.2 gives that e tA ∈ P for all t > 0. We introduce the following notation: S+ := S ∩ R 2 + , and S++ := S ∩ R 2 ++ . For y ∈ S+ we denote G(y) := A − (Ay · y)I y.
Recall that y ′ = G(y) is a system of ODEs satisfied by the directions of the solutions of x ′ = Ax (see Eq. (2.6)). We already know that G(y) is perpendicular to y, and that G(y) equals the zero vector if and only if y is the normalized principal eigenvector u of A. Otherwise, G(y) is a nonzero vector, pointing either clockwise or counterclockwise.
We check that We want to show that for any y ∈ S+ situated between [1 0] ⊤ and u, the vector G(y) points counterclockwise toward u, and for any y ∈ S+ situated between [0 1] ⊤ and u, the vector G(y) points clockwise toward u.
In order to prove that notice, first, that y ∈ S+ can be written in polar coordinates as [cos θ sin θ] ⊤ , where θ ∈ [0, π/2], and second, that
• For t > s the directions x(t)/ x(t) tend clockwise to u as t → ∞. This occurs when x(s)/ x(s) lies between [0 1] ⊤ and u.
• For t > s the directions x(t)/ x(t) tend counterclockwise to u as t → ∞. This occurs when x(s)/ x(s) lies between [1 0] ⊤ and u.
The bottom line is that at each time the directions of the solution tend toward the principal eigenvector.

Construction
In the present section we give a construction of a nonautonomous (piecewise constant) linear system x ′ = A(t)x of ODEs such that for each t ∈ R the larger eigenvalue of A(t) equals −1/2 but there is a solution not converging to zero as t → ∞.

Idea of the construction
We consider a system of linear ODEs with A(t) defined as where A (1) , A (2) are 2 by 2 matrices. Notice that A(t) has discontinuity points at integers. A solution of the system (3.1) is defined in the following way: It is a continuous function ξ : R → R 2 such that It is straightforward to see that for any s ∈ R and any x0 ∈ R 2 there exists a unique solution x(t; s, x0) of equation (3.1) satisfying the initial condition x(s) = x0. Further, we can write the transition matrix X(t; s) as X(t; s)x0 = x(t; s, x0).
The transition matrix has all the properties mentioned earlier, except that at t or s being integers its one-sided derivatives satisfy the suitable equalities.
We will not need, however, the above equality in its full generality.
From now on, we assume that A (1) and A (2) belong to M. Our program is to find two matrices, A (1) , A (2) ∈ M, such that their principal eigenvalues are negative, yet the set of those y ∈ S++ for which A (1) y · y > 0 and A (2) y · y > 0 is large. Indeed, then it is quite likely that for some solution x(t) the directions x(t)/ x(t) will be in that set for quite a large fraction of time (or, which would be the best, always), so the magnitude of that solution grows from time t = 0 to time t = 2 (and, by periodicity, it must grow to infinity as time goes to infinity).
Where to look for such matrices? Certainly not among symmetric (Hermitian) matrices, since for such matrices one can prove quite easily that, if the principal eigenvalue of A is negative then Ay · y < 0 for all nonzero y. So, a matrix should be far from symmetric.
If we find two matrices, A (1) , A (2) ∈ M, such that their principal eigenvalues are negative, nevertheless the principal eigenvalue of e A (2) e A (1) is larger than one, we take w(t) as a solution sought for: at t = 2n the magnitude of the solution goes exponentially, as n → ∞, to infinity.

Definition of A (1) and A (2)
We define parameterized by a parameter c > 0 (c will be taken to be large). Observe that c measures how far from symmetric the matrices A (1) and A (2) are. It is easy to see that the eigenvalues of the matrices A (1) and A (2) are −1/2 and −3/2.
Denote by u (1) the normalized principal eigenvector of A (1) , and u (2) the normalized principal eigenvector of A (2) ,

Why could the example be O.K.?
We want to show that, under the choice of the matrices A (1) and A (2) as in the previous subsection, it is very likely that there are plenty of solutions x(t) such that their length tends exponentially fast to infinity.
In order to do that, let us look at the set of those y ∈ S++ such that A (1) y · y > 0. As the matrix A (2) is the transpose of A (1) , that set will be equal to the set of those y ∈ S++ such that A (2) y · y > 0.
We have By writing y ∈ S++ in polar coordinates as [cos(θ) sin(θ)] ⊤ , θ ∈ (0, π/2), we obtain that After simple calculation we get A (1) y · y > 0 if and only if provided that c > 1 + √ 3 2 . Now let us apply the knowledge of how the directions of a solution change, as formulated in subsection 2.6. Assume that the initial value x(0) is situated somewhere between the principal eigenvectors for A (1) and A (2) . Recall that at each moment the direction tends toward the normalized principal eigenvector at that moment, from time t = 0 to time t = 1 the directions tend clockwise toward the normalized principal eigenvector u (1) of A (1) . They can leave the "red" set, but again from time t = 1 to time t = 2 they tend counterclockwise toward the normalized principal eigenvector u (2) of A (2) . By periodicity, the directions oscillate.

Analysis
The reasoning given in the previous subsubsection cannot be considered a formal proof. Now we give an analytical solution.
As for the instability it is enough to find one solution such that for some sequence of time moments its lengths tend to infinity.
How to look for such a solution? By Theorem 2.1, both matrices e A (2) and e A (1) belong to P, their product, that is, P , belongs to P, too.
By Proposition 2.1, there exists precisely one normalized principal eigenvector w of P pertaining to the principal eigenvalue, µ, of P .
So it is sufficient to check that the principal eigenvalue µ of P is larger than one.
Remark 3.2. The Floquet theory states that there is a decomposition where Q(t) is a time-periodic matrix function (with period 2) and R is a constant matrix. The eigenvalues of e tR are called characteristic multipliers of x ′ = A(t)x, and a ν ∈ C such that e 2ν is a characteristic multiplier is called a Floquet exponent of x ′ = A(t)x. Generally, Floquet exponents are not defined uniquely. But in our case µ is the (positive real) characteristic multiplier, larger than the other one, and its natural logarithm can be called the principal Floquet exponent of x ′ = A(t)x.
We proceed now to computing (or, rather, estimating from below) µ.
The exponential of A (1) is given by the formula .
One can find this formula by some Computer Algebra System. However, we prefer to give a more analytical explanation. We which is easily seen, by comparing the Maclaurin series expansions, to be equal to (cosh( t 2 ))I + (2 sinh( t 2 ))B. Similarly we have exp(tA (2) .

One has
.

Continuous matrix function
In contrast to the previous sections, the reader is now assumed to have experienced more exposure to "harder" mathematical thinking.
One could think that perhaps a phenomenon described above has something to do with the discontinuity at integer times. Results contained in the present section show that this is not so.
We start by recalling that, if for a matrix C = [cij ] 2 i,j=1 ∈ R 2×2 we denote its Euclidean norm as C =  Another fact is the Gronwall inequality: Then We proceed now to the construction. 2 )I. It is easily seen that the normalized principal eigenvector, us, ofB(s) is an eigenvector of B(s) pertaining to the eigenvaluẽ λ(s) − (λ(s) + 1 2 ) = − 1 2 . As B(s) belongs to M, − 1 2 must be therefore its principal eigenvalue.
Observe that where A (1) and A (2) are as in Section 3.
For ǫ ∈ (0, 1/4) we define a matrix function Aǫ : [0, 2] → M by the formula The function Aǫ is continuous, and the principal eigenvalue of A(t) is constantly equal to −1/2. We extend the matrix function Aǫ to the whole of R by periodicity (with period 2).
Denote by Xǫ(t; s) the transition matrix for the systemẋ = Aǫ(t)x: Xǫ(t; s)x0 is the solution of the initial value problem We want to show that the matrices Xǫ(2; 0) converge (entrywise), as ǫ → 0 + , to the matrix X(2; 0) as in Section 3. In fact, this is a special case of the continuous dependence of solutions of the initial value problem on parameters, as presented, for example, in [1, Chapter 17], but we prefer to give its (simple) proof here.
By simply integrating the relevant equations for the transition matrix we see that for any s, t ∈ R.
Consequently, with the help of standard estimates of integrals, together with (4.1), we obtain An application of the Gronwall inequality gives that Xǫ(t; 0) ≤ e M t and X(t; 0) ≤ e M t for t ∈ [0, 2]. Observe that it follows from the above inequality that, as ǫ → 0, all the entries of Xǫ(2; 0) converge to the corresponding entries of Xǫ(2; 0) = P . As, by Theorem 2.1, Xǫ(2; 0) belong to P, Proposition 2.1 implies that the principal eigenvalues, µǫ, of Xǫ(2; 0) converge to the principal eigenvalue, µ, of P , which is > 1. Consequently, for ǫ > 0 sufficiently close to zero the principal eigenvalue of Xǫ(2; 0) is > 1. It suffices now to take the solution ofẋ = Aǫ(t)x taking the value wǫ at t = 0, where wǫ is the normalized principal eigenvector of Xǫ(2; 0).

Smoother time dependence
Repeating an argument from Section 4 we can further approximate continuous matrix functions A(t) by matrix functions that are smooth, for example C 1 , or even C ∞ .

Extensions of results
Our construction, whether in Section 3 or in Section 4, apparently gives only one unstable solution.
In reality, however, one can prove, without much effort, more: • Not only w(k) = µ k for all k ∈ N, but also lim t→∞ ln w(t) t = µ.
• w(t) is by far not the only solution possessing the above property. Indeed, if x(t) denotes a nontrivial solution such that its initial value, x(0), is in R 2 + , then the directions x(t)/ x(t) converge, as t → ∞, to the directions w(t)/ w(t) , that is, (and the exponential rate of convergence is equal to half the natural logarithm of the second eigenvalue of the transition matrix X(2; 0)), from which it follows that lim t→∞ ln x(t) t = µ holds for such x(t), too.
One could mention here that the convergence of the directions is a phenomenon known in the literature on mathematical biology as asynchronous exponential growth. The above could be a material for undergraduate work.

Discrete counterparts
The example given in Section 3 can serve as an example in the discretetime case. Namely, we have a periodic sequence exp(A (1) ), exp(A (2) ), exp(A (1) ), exp(A (2) ), . . . of matrices in P, with spectral radii e −1/2 , but whose norms have exponential growth rate larger than 1. Such families model populations inhabiting various habitats, see, e.g., [6]. A biological interpretation of the results would be that, although the growth rate of the whole population is at each moment negative (namely, equal to −1/2), nevertheless the time-averaged growth rate, that is, ln µ, is positive.