45  Eigenvalues, eigenvectors

In the previous chapters, you’ve seen how linear dynamics, when unstable, lead trajectories off to infinity. Chapter 43 looked at some ways that nonlinearity can tame instability, as in the simple models of linear growth.

In this chapter, we return to linear dynamics to develop a quantitative theory of stability. Such theory is important in many engineering and design applications. For instance, a building exposed to earthquake risk can be economically designed to be strong specifically against the type of shaking produced by earthquakes. An electronic circuit can be designed to be sensitive to certain kinds of communication signals while still resisting noise or jamming.

Application area 45.1  

Application area 45.1 Cruise control

Cruise control is a feature of many modern automobiles. It enables a driver to set a desired speed for the car to follow automatically, thus reducing cognitive load on the driver. Early models of cruise control introduced in the 1970s (and still common in 2020) had an extremely simple logic. The speedometer constantly monitors the car’s speed $v$. The driver sets a desired speed, which we will call $v_0$. The cruise control actuator—the mechanism that works to maintain the speed at $v_0$ by, let’s say, adjusting the throttle to increase or reduce speed as appropriate—is designed to produce a throttle input proportional to the deviation in actual speed $v$ from the set speed $v_0$. $$\text{throttle-input}=-a(v-v_0).$$ You can see that when $v=v_0$, the throttle input from the system is zero. When $v_0<v$, as when starting to go downhill, the throttle input will be negative, reducing the speed.

A simple model is that the throttle input is proportional to ${\partial }_{t}v$: positive input increases speed, negative input decreases speed. The overall system is then described by a linear, first-order differential equation $${\partial }_{t}v=-a(v-v_0).$$ This equation has a fixed point at $v=v_0$ (as desired). To describe more simply the dynamics around the fixed point, we can define a new dynamical state $y=v-v_0$. Differentiating both sides with respect to time gives ${\partial }_{t}y={\partial }_{t}v$, so we can re-write the original differential equation as $${\partial }_{t}y=-ay$$ which has a solution $y(t)=A{e}^{-at}$ as we’ve seen before. For any positive coefficient $a$, the control system is stable: a deviation from $v_0$ will be followed by exponential decay back to $v_0$.

A modern cruise control system adds an important feature: the car slows down automatically if it gets too close to the car in front. Let’s call the distance to the leading car $x$ and define “too close” to mean closer than distance $x_0$. The car is traveling at velocity $v$ and the car in front at $v_f$. If $v_f<v$, the distance between the cars gets smaller: ${\partial }_{t}x=-(v-v_f)$. On the other hand, if $x_0<x$, the car is far enough away that the velocity can be increased: ${\partial }_{t}v=b(x-x_0)$. We can unite these two first-order differential equations into a single, second-order equation: $${\partial }_{t}x=-(v-v_f)⟹{\partial }_{tt}x=-{\partial }_{t}v=-b(x-x_0)$$ To solve this system, define $y\equiv x-x_0$. This implies ${\partial }_{tt}y={\partial }_{tt}x$, so the equation describing the control system becomes $${\partial }_{tt}y=-by.$$ The solution to this is $y(t)=A\sin \left(\sqrt{b}t\right)$ as you can confirm by substitution.

Result: The simple car-following control system oscillates. That is not a very good experience for the driver! By constructing a theory of stability, we may be able to figure out a re-design that will cause the oscillations to go away.

45.1 Vector solutions to linear differential equations

The form in which we have been writing linear differential equations in two state variables is $$\begin{array}{rc}{\partial }_{t}x& =ax+by\\ {\partial }_{t}y& =cx+dy.\end{array}$$

A key part of constructing a theory of stability is finding a set of mathematical ideas that enable us to view dynamics in a simpler way. The idea we will introduce here is thinking about the state and trajectory of a differential equation using vectors. We will work here with systems with a two-variable dynamical state, but the results apply just as well to higher dimensional states. That is important in applied work, where the systems being modeled are complicated with many state components.

We can re-write the linear differential equation using vector and matrix notation. Suppose that we collect the $x$ and $y$ components of the state into a vector, $$\stackrel{\rightharpoonup }{\mathbf{w}\mathbf{(}\mathbf{t}\mathbf{)}}=\left[\begin{array}{c}x(t)\\ y(t)\end{array}\right].$$ The differential equation, in terms of $\stackrel{\rightharpoonup }{\mathbf{w}\mathbf{(}\mathbf{t}\mathbf{)}}$ is $${\partial }_t\stackrel{\rightharpoonup }{\mathbf{w}\mathbf{(}\mathbf{t}\mathbf{)}}=\left[\begin{array}{cc}a& b\\ c& d\end{array}\right]\stackrel{\rightharpoonup }{\mathbf{w}\mathbf{(}\mathbf{t}\mathbf{)}}.$$ Now imagine that we pick two non-colinear vectors, $\stackrel{\rightharpoonup }{{\mathbf{u}}_{\mathbf{1}}}$ and $\stackrel{\rightharpoonup }{{\mathbf{u}}_{\mathbf{2}}}$ that span the state space. Since the vectors are assumed to span the state, any initial condition can be written as a linear combination of those two vectors: $$\stackrel{\rightharpoonup }{\mathbf{w}\mathbf{(}\mathbf{0}\mathbf{)}}=\left[\begin{array}{c}x(0)\\ y(0)\end{array}\right]=m_1\stackrel{\rightharpoonup }{{\mathbf{u}}_{\mathbf{1}}}+m_2\stackrel{\rightharpoonup }{{\mathbf{u}}_{\mathbf{2}}}.$$

For the moment, we won’t worry about how best to choose $\stackrel{\rightharpoonup }{{\mathbf{u}}_{\mathbf{1}}}$ and $\stackrel{\rightharpoonup }{{\mathbf{u}}_{\mathbf{2}}}$; any two vectors that are not colinear will do.

As a running example, we will work with the pair of first-order differential equations $$\begin{array}{rc}{\partial }_{t}x& =x+y\\ {\partial }_{t}y& =2x,\end{array}$$ which, in vector/matrix form are $${\partial }_t\stackrel{\rightharpoonup }{\mathbf{w}\mathbf{(}\mathbf{t}\mathbf{)}}=\left[\begin{array}{cc}1& 1\\ 2& 0\end{array}\right]\stackrel{\rightharpoonup }{\mathbf{w}\mathbf{(}\mathbf{t}\mathbf{)}}.$$ Imagine that we choose, arbitrarily, $$\stackrel{\rightharpoonup }{{\mathbf{u}}_{\mathbf{1}}}=\left[\begin{array}{r}1\\ -3\end{array}\right]\text{and}\stackrel{\rightharpoonup }{{\mathbf{u}}_{\mathbf{2}}}=\left[\begin{array}{r}1.0\\ 0\end{array}\right].$$

For the example, we will calculate a trajectory starting at the initial condition $\stackrel{\rightharpoonup }{\mathbf{w}\mathbf{(}\mathbf{0}\mathbf{)}}=\left[\begin{array}{r}-1.1\\ 2.1\end{array}\right]$:

traj <- integrateODE(dx ~ x + y, dy ~ 2*x, # dynamics
                     x=-1.1, y=2.1, # initial conditions
                     domain(t=0:2))
traj_plot(y(t) ~ x(t), traj)

## Warning: All aesthetics have length 1, but the data has 500 rows.
## ℹ Please consider using `annotate()` or provide this layer with data containing
##   a single row.
## All aesthetics have length 1, but the data has 500 rows.
## ℹ Please consider using `annotate()` or provide this layer with data containing
##   a single row.

Figure 45.1: The trajectory calculated starting at (-1.1, 2.1). The graph is annotated with the vectors $vec{u}_1$ and $vec{u}_2$ and the flow field.

The initial condition (marked “0” in Figure 45.1 is, like any other point in the state space, a linear combination of $\stackrel{\rightharpoonup }{{\mathbf{u}}_{\mathbf{1}}}$ and $\stackrel{\rightharpoonup }{{\mathbf{u}}_{\mathbf{2}}}$. We can find the scalar coefficients of the linear combination using any of the methods presented in Block 5, for instance the telescope method. We will illustrate with qr.solve():

1

2

3

M <- cbind(rbind(1,-3), rbind(1,0))

w0 <- rbind(-1.1, 2.1)

qr.solve(M, w0)

So, $\stackrel{\rightharpoonup }{\mathbf{w}\mathbf{(}\mathbf{0}\mathbf{)}}=-0.7\stackrel{\rightharpoonup }{{\mathbf{u}}_{\mathbf{1}}}-0.4\stackrel{\rightharpoonup }{{\mathbf{u}}_{\mathbf{2}}}$. Keep these scalar coefficients, $-0.7$ and $-0.4$ in mind for the next example.

We can use integrateODE() to find the solution starting at any initial condition. In particular, we can find the solution $\stackrel{\rightharpoonup }{{\mathbf{u}}_{\mathbf{1}}}$ as the initial condition and, similarly, using $\stackrel{\rightharpoonup }{{\mathbf{u}}_{\mathbf{2}}}$ as the initial condition.

traj_u1 <- integrateODE(
  dx ~ x + y, dy ~ 2*x, # dynamics
  x=1, y=-3,  #initial conditions
  domain(t = 0:2))
traj_u2 <- integrateODE(
  dx ~ x + y, dy ~ 2*x, #dynamics
  x=1, y= 0,  #initial conditions
  domain(t = 0:2))

Figure 45.2 shows these trajectories.

## Warning: All aesthetics have length 1, but the data has 500 rows.
## ℹ Please consider using `annotate()` or provide this layer with data containing
##   a single row.
## All aesthetics have length 1, but the data has 500 rows.
## ℹ Please consider using `annotate()` or provide this layer with data containing
##   a single row.

Figure 45.2: Trajectories from the initial conditions $\stackrel{\rightharpoonup }{{\mathbf{u}}_{\mathbf{1}}}$ and $\stackrel{\rightharpoonup }{{\mathbf{u}}_{\mathbf{2}}}$.

At first glance, the two trajectories $\stackrel{\rightharpoonup }{{\mathbf{u}}_{\mathbf{1}}\mathbf{(}\mathbf{t}\mathbf{)}}$ and $\stackrel{\rightharpoonup }{{\mathbf{u}}_{\mathbf{2}}\mathbf{(}\mathbf{t}\mathbf{)}}$ in Figure 45.2 that start from $\stackrel{\rightharpoonup }{{\mathbf{u}}_{\mathbf{1}}}$ and $\stackrel{\rightharpoonup }{{\mathbf{u}}_{\mathbf{2}}}$ might not look much like the trajectory in Figure 45.1 that starts from $\stackrel{\rightharpoonup }{\mathbf{w}\mathbf{(}\mathbf{0}\mathbf{)}}=-0.7\stackrel{\rightharpoonup }{{\mathbf{u}}_{\mathbf{1}}}-0.4\stackrel{\rightharpoonup }{{\mathbf{u}}_{\mathbf{2}}}$. But in fact there is a very simple relationship between the trajectories: $$\stackrel{\rightharpoonup }{\mathbf{w}\mathbf{(}\mathbf{t}\mathbf{)}}=-0.7\stackrel{\rightharpoonup }{{\mathbf{u}}_{\mathbf{1}}\mathbf{(}\mathbf{t}\mathbf{)}}-0.4\stackrel{\rightharpoonup }{{\mathbf{u}}_{\mathbf{2}}\mathbf{(}\mathbf{t}\mathbf{)}}.$$ To state the situation more generally, any solution to the differential equations can be written as a linear combination of the solutions starting at $\stackrel{\rightharpoonup }{{\mathbf{u}}_{\mathbf{1}}}$ and $\stackrel{\rightharpoonup }{{\mathbf{u}}_{\mathbf{2}}}$, regardless of how $\stackrel{\rightharpoonup }{{\mathbf{u}}_{\mathbf{1}}}$ and $\stackrel{\rightharpoonup }{{\mathbf{u}}_{\mathbf{2}}}$ were chosen.

We can see this algebraically. Since $\stackrel{\rightharpoonup }{{\mathbf{u}}_{\mathbf{1}}\mathbf{(}\mathbf{t}\mathbf{)}}$ and $\stackrel{\rightharpoonup }{{\mathbf{u}}_{\mathbf{2}}\mathbf{(}\mathbf{t}\mathbf{)}}$ are solutions to the linear differential equation, it must be that $${\partial }_t\stackrel{\rightharpoonup }{{\mathbf{u}}_{\mathbf{1}}\mathbf{(}\mathbf{t}\mathbf{)}}=\left[\begin{array}{cc}1& 1\\ 2& 0\end{array}\right]\stackrel{\rightharpoonup }{{\mathbf{u}}_{\mathbf{1}}\mathbf{(}\mathbf{t}\mathbf{)}}\text{and}{\partial }_t\stackrel{\rightharpoonup }{{\mathbf{u}}_{\mathbf{2}}\mathbf{(}\mathbf{t}\mathbf{)}}=\left[\begin{array}{cc}1& 1\\ 2& 0\end{array}\right]\stackrel{\rightharpoonup }{{\mathbf{u}}_{\mathbf{2}}\mathbf{(}\mathbf{t}\mathbf{)}}.$$ Taking a linear combination of these equations gives
$${\partial }_t[m_1\stackrel{\rightharpoonup }{{\mathbf{u}}_{\mathbf{1}}\mathbf{(}\mathbf{t}\mathbf{)}}+m_2\stackrel{\rightharpoonup }{{\mathbf{u}}_{\mathbf{2}}\mathbf{(}\mathbf{t}\mathbf{)}}]=\left[\begin{array}{cc}1& 1\\ 2& 0\end{array}\right][m_1\stackrel{\rightharpoonup }{{\mathbf{u}}_{\mathbf{1}}\mathbf{(}\mathbf{t}\mathbf{)}}+m_2\stackrel{\rightharpoonup }{{\mathbf{u}}_{\mathbf{2}}\mathbf{(}\mathbf{t}\mathbf{)}}]$$ The same will be true in general, that is, for the matrix $\left[\begin{array}{cc}a& b\\ c& d\end{array}\right]$.

45.2 Eigenvectors and eigenvalues

In the previous section, we saw that the solution to any linear differential equation starting at any initial condition can be written as a linear combination $m_1\stackrel{\rightharpoonup }{{\mathbf{u}}_{\mathbf{1}}\mathbf{(}\mathbf{t}\mathbf{)}}+m_2\stackrel{\rightharpoonup }{{\mathbf{u}}_{\mathbf{2}}\mathbf{(}\mathbf{t}\mathbf{)}}$, where $\stackrel{\rightharpoonup }{{\mathbf{u}}_{\mathbf{1}}\mathbf{(}\mathbf{t}\mathbf{)}}$ is the solution starting at an initial condition $\stackrel{\rightharpoonup }{{\mathbf{u}}_{\mathbf{1}}}$ and $\stackrel{\rightharpoonup }{{\mathbf{u}}_{\mathbf{2}}\mathbf{(}\mathbf{t}\mathbf{)}}$ is the solution starting at $\stackrel{\rightharpoonup }{{\mathbf{u}}_{\mathbf{2}}}$. It does not matter how $\stackrel{\rightharpoonup }{{\mathbf{u}}_{\mathbf{1}}}$ and $\stackrel{\rightharpoonup }{{\mathbf{u}}_{\mathbf{2}}}$ are chosen, so long as they are not colinear, that is, so long as they span the state space.

In this section, we demonstrate that there is a particular way of selecting $\stackrel{\rightharpoonup }{{\mathbf{u}}_{\mathbf{1}}}$ and $\stackrel{\rightharpoonup }{{\mathbf{u}}_{\mathbf{2}}}$ that makes the solutions $\stackrel{\rightharpoonup }{{\mathbf{u}}_{\mathbf{1}}\mathbf{(}\mathbf{t}\mathbf{)}}$ and $\stackrel{\rightharpoonup }{{\mathbf{u}}_{\mathbf{2}}\mathbf{(}\mathbf{t}\mathbf{)}}$ have a very simple, purely exponential format. The vectors to be chosen are the eigenvectors of the matrix $\left[\begin{array}{cc}a& b\\ c& d\end{array}\right]$. We will call these eigenvectors $\stackrel{\rightharpoonup }{{\mathbf{\Lambda }}_{\mathbf{1}}}$ and $\stackrel{\rightharpoonup }{{\mathbf{\Lambda }}_{\mathbf{2}}}$. (This use of the Greek letter $\mathrm{\Lambda }$ (capital “lambda”) and it is lower-case version $\lambda$, is conventional in mathematics, physics, and engineering. So it is worth learning to identify the letters.)

Our task in this section is to show how to compute the eigenvectors $\stackrel{\rightharpoonup }{{\mathbf{\Lambda }}_{\mathbf{1}}}$ and $\stackrel{\rightharpoonup }{{\mathbf{\Lambda }}_{\mathbf{2}}}$ and that the solutions $\stackrel{\rightharpoonup }{{\mathbf{\Lambda }}_{\mathbf{1}}\mathbf{(}\mathbf{t}\mathbf{)}}$ and $\stackrel{\rightharpoonup }{{\mathbf{\Lambda }}_{\mathbf{2}}\mathbf{(}\mathbf{t}\mathbf{)}}$ are in fact simple exponentials. Chapter 46 derives the formula for the eigenvectors. Here, we use the R function eigen() to do the calculations for us.

Try it! 45.1

Try it! 45.1 Calculating eigenvectors

Eigenvectors can be calculated using the R function eigen() applied to the abcd matrix that defines the linear differential equation.

For the system of first-order differential equations $$\begin{gathered} {\partial }_{t}x=x+y \\ {\partial }_{t}y=2x \end{gathered}$$ the matrix is, as we’ve seen, $$\left[\begin{array}{cc}1& 1\\ 2& 0\end{array}\right].$$

Carrying out the eigenvector calculation is straightforward:

1

2

M <- cbind(rbind(1,2), rbind(1,0))

eigen(M)

The eigenvectors are the two columns of the matrix labeled vectors returned by the calculation. Here, that is
$$\stackrel{\rightharpoonup }{{\mathbf{\Lambda }}_{\mathbf{1}}}=\left[\begin{array}{r}0.7071\\ 0.7071\end{array}\right]\text{and}\stackrel{\rightharpoonup }{{\mathbf{\Lambda }}_{\mathbf{2}}}=\left[\begin{array}{r}-0.4472\\ 0.8944\end{array}\right].$$

The calculation also produces eigenvalues. Here that is ${\lambda }_1=2$ and ${\lambda }_2=-1$.

We can see what’s special about $\stackrel{\rightharpoonup }{{\mathbf{\Lambda }}_{\mathbf{1}}}$ and $\stackrel{\rightharpoonup }{{\mathbf{\Lambda }}_{\mathbf{2}}}$ by plotting them along with the flow field, as in Figure 45.3.

Figure 45.3: The eigenvectors for the $\left[\begin{array}{cc}1& 1\\ 2& 0\end{array}\right]$ plotted along with the flow.

The eigenvectors mark the directions where the flow is either directly toward the fixed point or directly away from it. Here, the flow on the subspace of $\stackrel{\rightharpoonup }{{\mathbf{\Lambda }}_{\mathbf{1}}}$ is away from the fixed point, while the flow along the subspace of $\stackrel{\rightharpoonup }{{\mathbf{\Lambda }}_{\mathbf{2}}}$ is inward to the fixed point.

The consequence of this alignment of the flow with the eigenvectors is that the trajectory from any initial condition $m_1\stackrel{\rightharpoonup }{{\mathbf{\Lambda }}_{\mathbf{1}}}$ will have the form $m_1(t)\stackrel{\rightharpoonup }{{\mathbf{\Lambda }}_{\mathbf{1}}}$ and similarly for an initial condition $m_2(t)\stackrel{\rightharpoonup }{{\mathbf{\Lambda }}_{\mathbf{2}}}$.

As we did in the previous section, let’s calculate the trajectories $\stackrel{\rightharpoonup }{{\mathbf{\Lambda }}_{\mathbf{1}}\mathbf{(}\mathbf{t}\mathbf{)}}$ and $\stackrel{\rightharpoonup }{{\mathbf{\Lambda }}_{\mathbf{2}}\mathbf{(}\mathbf{t}\mathbf{)}}$ starting at the two eigenvectors and plot out the $y(t)$ component of the solution. Since we are anticipating an exponential form for the function, we use semi-log axes, where an exponential looks like a straight line.

traj_eigen1 <- integrateODE(
  dx ~ x + y, dy ~ 2*x,  # dynamics
  x = 0.7071, y = 0.7071, # initial conditions
  domain(t = 0:1) 
  )
## Solution containing functions x(t), y(t).
traj_eigen2 <- integrateODE(
  dx ~ x + y, dy ~ 2*x,  # dynamics
  x=-0.4472, y=0.8944,   # initial conditions
  domain(t = 0:1))
## Solution containing functions x(t), y(t).

traj_plot(y(t) ~ t, traj_eigen1, color="magenta") |>
 traj_plot(y(t) ~ t, traj_eigen2, color="blue") |> 
  gf_refine(scale_y_log10(
    breaks=c(0.3290, 0.7071, 0.8944, 5.2248)))

Figure 45.4: Two time series, one showing stability, the other instability.

We have marked the $y$ axis with the starting and ending values of each function, so that you can find the exponential parameter $k$ for each function.

$y_1(t)=0.7071e^{k_{1}t}$

$y_2(t)=0.8944e^{k_{2}t}$.

To find $k_1$ and $k_2$, plug in $t=1$ to the solution:

$y_1(1)=5.2248=0.7071e^{k_1}⟹k_1=2$

$y_2(1)=0.3290=0.8944e^{k_2}⟹k_2=-1$

Look back at the results from the eigen(M) calculation. These values for $k_1$ and $k_2$ are exactly the eigenvalues that were computed from the matrix M. In standard notation, rather than $k_1$ and $k_2$, the notation ${\lambda }_1=k_1$ and ${\lambda }_2=k_2$ is preferred. (Remember, $\lambda$ is the Greek letter “lambda” in it is lower-case form.) Every solution to the differential equation has the form $$m_1{e}^{{\lambda }_{1}t}\stackrel{\rightharpoonup }{{\mathbf{\Lambda }}_{\mathbf{1}}}+m_2{e}^{{\lambda }_{2}t}\stackrel{\rightharpoonup }{{\mathbf{\Lambda }}_{\mathbf{2}}}.$$ The scalar coefficients $m_1$ and $m_2$ can be found from the initial condition. The stability of the system depends only on ${\lambda }_1$ and ${\lambda }_2$. If either one of these is positive, then the system is unstable.