2.4 Linear oscillators
2.5.2 Stability analysis at equilibrium
y0(t) = Ay(t)
we try a solution of the form y(t) = eλtv, where v ∈ R2is a constant. Plugging into (2.25), we get λveλt = Aveλt.
We find that y(t) = eλtv is a solution of (2.25) if and only if
Av = λv. (2.27)
That is, λ is the eigenvalue and v is the corresponding eigenvector. The eigenvalue λ satisfies the following characteristic equation
det (λI − A) = 0.
In two dimension, this is
λ2− T λ + D = 0, where
T = a + d, D = ad − bc are the trace and determinant of A, respectively. The eigenvalues are
λ1= T +√
T2− 4D
2 , λ2= T −√
T2− 4D
2 .
There are three possibilities for the eigenvalues:
• Case 1: T2− 4D > 0. Then λ16= λ2and are real.
• Case 2: T2− 4D < 0. Then λ1, λ2and are complex conjugate.
• Case 3: T2− 4D = 0. Then λ1is a double root.
Case 1.
λi are real and there are two independent real eigenvectors v1 and v2. The corresponding two independent solutions are
y1 = eλ1tv1, y2 = eλ2tv2. A general solution has the form
y(t) = C1y1(t) + C2y2(t) If the initial data is y0, then
C1v1+ C2v2= y0. Let T be the matrix (v1, v2). Then
C1 C2
= T−1y0.
The 0 state is an equilibrium. Its behavior is determined by the sign of the eigenvalues:
• λ1, λ2 < 0: all solutions tend to 0 as t → ∞. We call 0 state a sink. It is a stable equilibrium.
• λ1, λ2 > 0: all solutions tend to infinity as t → ∞. In fact, all solutions tend to the 0 state as t → −∞. We call 0 state a source. It is an unstable equilibrium.
• λ1· λ2< 0. Let us take λ1< 0 and λ2 > 0 as an example for explanation. A general solution has the form
y(t) = C1eλ1tv1+ C2eλ2tv2
We have eλ1t→ 0 and eλ2t → ∞ as t → ∞. Hence if y(0) ∈ Ms := {γv1, γ ∈ R}, then the corresponding C2 = 0, and y(t) → 0 as t → ∞. We call the line Msa stable manifold.
2.5. 2 × 2 LINEAR SYSTEMS 55 On the other hand, if y(0) ∈ Mu := {γv2, γ ∈ R}, then the corresponding C1 = 0 and y(t) → 0 as t → −∞. We call the line Mu an unstable manifold. For any other y0, the corresponding y(t) has the following asymptotics:
y(t) → v1-axis, as t → −∞, y(t) → v2-axis, as t → +∞.
That is, all solutions approach the stable manifold as t → ∞ and the unstable manifold as t → −∞. The 0 state is the intersection of the stable and unstable manifolds. It is called a saddle point.
• λ1 = 0 and λ2 6= 0. In this case, a general solution has the form: y(t) = C1v1+ C2eλ2tv2. The equilibrium {¯y|A¯y = 0} is a line: {C1v1|C1 ∈ R}. When λ2 < 0, then all solutions approach C1v1. This means that the line C1v1is a stable line.
Examples
1. Consider the matrix
A =
1 1 4 1
. The corresponding characteristic equation is
det (λI − A) = (λ − 1)2− 4 = 0.
Hence, the two eigenvalues are
λ1= 3, λ2= −1.
The eigenvector v1corresponding to λ1 = 3 satisfies (A − λ1I)v1= 0.
This gives
v1 =
1 2
. Similarly, the eigenvector corresponding to λ2 = −1 is
v2 =
1
−2
.
2. Consider
y0 = Ay, A =
8 −11 6 −9
.
The eigenvalues of A are roots of the characteristic equation det (λI − A) = 0. This yields
The line in the direction of v1 is a stable manifold, whereas the line in v2 direction is a unstable manifold. The origin is a saddle point.
3. Consider
λiare complex conjugate.
λ1= α + iω, λ2= α − iω.
Since A is real-valued, the corresponding eigenvectors are also complex conjugate:
w1 = u + iv, w2= u − iv.
We have two independent complex-valued solutions: z1 = eλ1tw1and z2 = eλ2tw2.
Since our equation (2.25) has real coefficients, its real-valued solution can be obtained by taking the real part (or pure imaginary part ) of the complex solution. In fact, suppose z(t) = x(t) + iy(t) is a complex solution of the real-value ODE (2.25). Then
d
dt(x(t) + iy(t)) = A (x(t) + iy(t)) .
By taking the real part and the imaginary part, using the fact that A is real, we obtain dx
dt = Ax(t), dy
dt = Ay(t)
2.5. 2 × 2 LINEAR SYSTEMS 57 Hence, both the real part and the imaginary part of z(t) satisfy the equation.
Now, let us take the real part and the imaginary part of one of the above solution:
z1(t) = eαt(cos ωt + i sin ωt) (u + iv) Its real part and imaginary part are respectively
y1(t) = eαt(cos ωtu − sin ωtv) y2(t) = eαt(sin ωtu + cos ωtv)
The other solution z2is the complex conjugate of z1. We extract the same real solutions from z2. You may wonder now whether u and v are independent. Indeed, if v = cu for some c ∈ R, then
A(u + iv) = λ1(u + iv) gives
A(1 + ic)u = λ1(1 + ic)u Au = λ1u = (α + iω)u This yields
Au = αu, and ωu = 0,
because A is real. This implies ω = 0 if u 6= 0. This contradicts to that the eigenvalue λ1 has nontrivial imaginary part.
From the independence of u and v, we conclude that y1 and y2 are also independent, and constitute a real basis in the solution space S. A general solution is given by
y(t) = C1y1(t) + C2y2(t), where Ciare real and are determined by the initial data:
C1y1(0) + C2y2(0) = y0, C1u + C2v = y0. Let T be the matrix (u, v), then
C1 C2
= T−1y0. We may use another parameters to represent the solution.
y(t) = C1eαt(cos ωtu − sin ωtv) + C2eαt(sin ωtu + cos ωtv)
= eαt((C1cos ωt + C2sin ωt)u + (C2cos ωt − C1sin ωt)v)
= Aeαt(cos(ωt − ω0)u + sin(ωt − ω0)v) , where (C1, C2) = A(cos ω0, sin ω0).
There are three cases for the structure of the solutions.
• α = 0: The eigenvalues are pure imaginary. All solutions are ellipses.
• α < 0: The solution are spirals and tend to 0 as t → ∞. The 0 state is a spiral sink.
• α > 0: The solution are spirals and tend to 0 as t → −∞. The 0 state is a spiral source.
Example
7)/2. The corresponding eigenvectors are v1 = real parts and imaginary parts. They are
y1 = et/2(cos(ωt)u − sin(ωt)w) , y2 = et/2(sin(ωt)u + cos(ωt)w) where ω =√
7/2.
Case 3.
λ1 = λ2are real and only one eigenvector. Let us see examples first to get some intuition.
Examples eigenvector. The y2component satisfies single equation
y02 = ry2.
We obtain y2(t) = C2ertand By plugging this into the first equation y10 = ry1+ C2ert,
we find y1(t) = C2tertis a special solution. The general solution of y1is y1(t) = C2tert+ C1ert. We can express these general solutions in vector form:
y(t) = C1ert
2.5. 2 × 2 LINEAR SYSTEMS 59 has a double root λ = 2. The corresponding eigenvector satisfies
(A − 2I)v = 0 This yields a solution, called v1:
v1 =
1
−1
.
This is the only eigenvector. The solution e2tv1 is a solution of the ODE. To find the other independent solution, we expect that there is a resonant solution te2tin the direction of v1. Unfortunately, te2tv1is not a solution unless v1 = 0. Therefore, we try to another
y(t) = te2tv1+ eµtv2,
Now, we find two solutions
y1 = e2tv1
y2 = te2tv1+ e2tv2.
Remark. The double root case can be thought as a limiting case where the second eigenvalues λ2 → λ1and v2 → v1. In this case,
1 λ2− λ1
eλ2tv2− eλ1tv1
is also a solution. This is equivalent to differentiate eλtv
in λ at λ1, where v(λ) is the eigenvector corresponding to λ. This derivative is teλ1tv1+ eλ1t∂v
∂λ
The new vector∂v∂λ is denoted by v2. By plugging teλ1tv1+ eλ1tv2into the equation, we obtain v2
should satisfies
(A − λ1I)v2= v1.
Let us return to study general equation: y0 = Ay with eigenvalues of A being a double root (i.e. λ1 = λ2). Let us call the corresponding eigenvector v1. Our goal is to find another vector v2
such that
(A − λ1I)v2= v1. The characteristic equation is
p(λ) := det(λI − A) = (λ − λ1)2 = 0.
The Caley-Hamilton theorem states that A satisfies the matrix equation:
p(A) = 0.
This can be seen from the following argument. Let Q(λ) be the adjugate matrix of A − λI, i.e.
Q(λ) =
d − λ −b
−c a − λ
=
d −b
−c a
− λI.
This adjugate matrix commutes with A (check by yourself). Further, (A − λI)Q(λ) = Q(λ)(A − λI) = p(λ)I.
This is a polynomial in λ with matrix coefficients. The coefficients commute with A. When we plug λ = A, we immediately get p(A) = 0.
Now, suppose λ1is a double root of A. We can find an eigenvector v1such that Av1 = λ1v1.
2.5. 2 × 2 LINEAR SYSTEMS 61 We cannot find two independent eigenvectors corresponding to λ1unless A = λ1I (why?). Yet, we can find another vector v2such that
(A − λ1I)v2= v1.
The solvability of v2 comes from the follows. Let Nkbe the null space of (A − λ1I)k, k = 1, 2.
We have the following mapping
N2 A−λ−→ N1I 1 A−λ−→ {0}.1I
We have seen that the only eigenvector is v1. Thus, N1 =< v1 >, the span of v1. In the map A − λ1I : N2 → N1, the domain space is N2 = R2 from Caley-Hamilton theorem. We have seen that the kernel is N1, which has dimension 1. Therefore the range space has dimension 1. Here, we use a theorem of linear map: the sum of the dimensions of range and kernel spaces equals the dimension of the domain space. We conclude that the range (A − λ1I)N2has to be N1. Therefore, there exists a v2 ∈ N2such that
(A − λ1I)v2= v1. The matrix A, as represented in the basis v1 and v2, has the form
A[v1, v2] = [v1, v2]
λ1 1 0 λ1
This is called the Jordan canonical form of A. We can find two solutions from this form:
y1(t) = eλ1tv1,
y2(t) = teλ1tv1+ eλ1tv2
You can check the Wronskian W [y1, y2](t) 6= 0. Thus, y1 and y2 form a fundamental solution.
The general solution has the form
y(t) = C1y1(t) + C2y2(t).
The stability of the 0 state (called the critical state) relies on the sign of λ1. We have
• λ1 < 0: the 0 state is a stable equilibrium.
• λ1 > 0: the 0 state is an unstable equilibrium.
• λ1 = 0: the general solution reads
y(t) = C2tv2+ C1v1, which tends to ∞ as t → ∞. Therefore, the 0 state is “unstable.”
Summary of stability analysis We can plot a stability diagram on the plane of the two parameters T and D, the trace and the determinant of A. The eigenvalues of A are
λ1= T +√
T2− 4D
2 , λ2= T −√
T2− 4D
2 .
Let ∆ := T2− 4D. On the T -D plane, the parabola ∆ = 0, the line D = 0 and the line T = 0 partition the plane into the following regions. The status of the origin is as the follows.
• ∆ > 0, D < 0, the origin is a saddle point.
• ∆ > 0, D > 0, T > 0, the origin is an unstable node (source).
• ∆ > 0, D > 0, T < 0, the origin is an stable node (sink).
• ∆ < 0, T < 0, the origin is a stable spiral point.
• ∆ < 0, T > 0, the origin is an unstable spiral point.
• ∆ < 0, T = 0, the origin is an stable circular point.
• ∆ = 0, T < 0, the origin is a stable node.
• ∆ = 0, T > 0, the origin is an unstable node.
• ∆ = 0, T = 0, the origin is an unstable node.
Homeworks.
1. Consider A =
1 −2 3 −4
. Find the exact solution of y0 = Ay and analyze the stability of 0.
2. Consider A =
3 6
−1 −2
. Find the exact solution of y0 = Ay and analyze the stability of 0.
3. Solve the circuit system
I V
0
−RL1 −L1
1
C −CR1
2
I V
and analyze the stability of the 0 state.
Chapter 3
Nonlinear systems in two dimensions
3.1 Biological models
3.1.1 Lotka-Volterra system