3 Linearisation and equilibria

3.2 Linear systems

Suppose \(x(t)\) satisfies the linear ODEs \(\dot{x}= Ax\), where \(A\) is a constant \(n \times n\) matrix and \(x \in \mathbb{R}^n\). If \(A\) has distinct eigenvalues \(\lambda_i\) with corresponding eigenvectors \( e_i\), then the general solution is a superposition of the eigenmodes. That is, \[ x(t)=\sum_{k=1}^n C_k\exp(\lambda_k t)e_k, \] where the \(C_k\) are constants determined from the initial condition \(x(0) = \sum C_k e_k\). This shows that eigenvalues and eigenvectors of \(A\) will be crucial to the understanding of the dynamics.

The eigen-pairs are also closely related to a particular coordinate transformation that simplifies the dynamics: \(x = [e_1\ e_2\ \cdots\ e_n]y\), or \(y=Ux\) with \(U=[e_1\ e_2\ \cdots\ e_n]^{-1}\), the inverse of the matrix formed by the eigenvectors. The the ODE for \(y\) becomes \(\dot{y} = U \dot{x} = UAx = UAU^{-1} y,\) i.e. \(y\) satisfies a linear ODE \(\dot{y}=UAU^{-1}y\).

The above choice of \(U=[e_1 \ e_2 \ \cdots \ e_n]^{-1}\) is particular in that \(UAU^{-1}\) is diagonal, and the system \(\dot{y} = UAU^{-1}y\) is essential \(n\) decoupled ODEs: \[ \dot{y}_1 = \lambda_1 y_1,\quad \dot{y}_2 = \lambda_2 y_2,\quad \cdots\quad \dot{y}_n = \lambda_n y_n. \] With this `natural' choice of transformation \(y=Ux\), the resulting system for \(y\) is called the normal form (depending only on the eigenvalues of \(A\)). We will work in the plane \(\mathbb{R}^2\) with real matrix \(A\), though extension to \(\mathbb{R}^n\) is not hard.

Now we discuss the behavior of the \(2\times 2\) system \(\dot{y}_1 = \lambda_1 y_1, \dot{y}_2 = \lambda_2y_2\) ( \(\lambda_1\) and \(lambda_2\) could be complex).

a) eigenvalues real and distinct: Suppose the eigenvalues \(\lambda_1\) and \(\lambda_2\) of \(A\) are real and distinct, then their corresponding eigenvectors \( e_1, e_2\) (assumed to be column vectors) are real and linearly independent. With the matrix \(U=[e_1, e_2]^{-1}\), we get \[ AU^{-1} = A[e_1\ e_2] = [Ae_1 \ Ae_2] = [\lambda_1 e_1 \ \lambda_2 e_2] =[e_1 \ e_2]\begin{bmatrix} \lambda_1 & \cr & \lambda_2\end{bmatrix} =U^{-1}\mbox{diag}(\lambda_1,\lambda_2). \] Left multiplying both sides with \(U\), we have \( UAU^{-1} = \mbox{diag}(\lambda_1,\lambda_2) \) as expected. As we shall see shortly, this transformation into normal form also makes it easier to understand the structure of the solutions, which depends on the signs of \(\lambda_1\) and \(\lambda_2\).

a i) \(\lambda_1< \lambda_2 < 0\): [stable node] In this case the linearisation in the normal form coordinates \(y^t=( u,v)\) is \begin{equation*} \dot{u}=\lambda_1 u, \quad \dot{v}=\lambda_2 v \end{equation*} with solutions \begin{equation*} u=u _0 e^{ \lambda_1 t}, \quad v=v _0 e^{ \lambda_2 t}. \end{equation*} Thus \((u,v) \to (0,0)\) as \(t \to \infty\) and both coordinate axis (\(u=0\) and \(v=0\)) are invariant. Moreover, if \(u_0, v_0\ne 0\) (i.e. off the coordinate axes) \( \frac{u}{u_0}=e ^{\lambda_1 t}=( e^{\lambda_2 t})^{\frac{\lambda_1}{\lambda_2}}; \qquad \frac{v}{v_0}=e ^{\lambda_2 t} \) and so \begin{equation}\label{parabola} \frac{u}{u_0}=\left( \frac{v}{v_0}\right)^{\frac{\lambda_1}{\lambda_2}}, \qquad \qquad \frac{\lambda_1}{\lambda_2}> 1, \end{equation} or equivalently \(uv^{-\lambda_1/\lambda_2}\) is a constant for points on the same trajectory. These are generalized parabolas, tangential to the \(v\)-axis at \((u,v)=(0,0).\)

Figure 3.2 Stable node in transformed \((u,v)\)-coordinates and in the original \((x_1,x_2)\)-coordinates.

This is called a stable node. In the original coordinates the \(u\)-axis corresponds to \(e_1\) and the \(v\)-axis to the \(e_2\) direction (you can see this from the transformation \(y=Ux\)), so the phase portrait is as shown in Figure 3.2. Thus in the original coordinates, lines corresponding to eigenvectors are invariant. Moreover almost all trajectories are tangential to \(e_2\) at \((0,0)\), i.e. tangential to eigenvector of eigenvalue with smallest modulus.

a ii) \(\lambda_1 > \lambda_2 > 0\): unstable node The phase portrait can be obtained using the same techniques as in the previous section. Indeed the manipulations are the same and the generalized parabola is also the same as changing the signs of both eigenvalues does not change the sign of their ratio. Another way of seeing the direct correspondence with the previous case is by reversing time. Set \(\tau = -t\) so \( \frac{\mathrm{d}}{\mathrm{d} t}=\frac{\mathrm{d}\tau}{\mathrm{d} t}\frac{\mathrm{d}}{\mathrm{d}\tau }=- \frac{\mathrm{d}}{\mathrm{d}\tau } \) and so if \begin{equation*} \frac{d}{dt} u = \lambda_1 u, \quad \frac{d}{dt}v = \lambda_2 v \end{equation*} then \begin{equation*} \frac{d}{d\tau} u = - \lambda_1 u, \quad \frac{d}{d\tau} v = - \lambda_2 v \end{equation*} which is the same as in case ai). Thus all we need to do is to change the direction of time, i.e. the arrows on the phase portraits to get the new phase portrait. This is called an unstable node, as shown in Figure 3.3.

Figure 3.3 Unstable node in transformed \((u,v)\)-coordinates and in the original \((x_1,x_2)\)-coordinates.

a iii) \(\lambda_1 < 0 < \lambda_2 \): saddle The analysis is as before but now \[\frac{u}{u_0}=(\frac{v}{v_0})^{\frac{\lambda_1}{\lambda_2}} \] is a generalized hyperbola as \(\lambda_1/\lambda_2<0\), as shown in Figure 3.4.

Figure 3.4 Stable node in transformed \((u,v)\)-coordinates and in the original \((x_1,x_2)\)-coordinates.

b) Complex conjugate eigenvalues \(\rho \pm i \omega, \omega \neq 0\) The eigenvectors \({z}_\pm\) are complex, and satisfy \begin{equation}\label{eq:cmplxeig} A{z}_\pm=(\rho\pm i\omega){z}_\pm. \end{equation} But we prefer to work with real quantities, and the first step is to take the real and imaginary parts of both sides of Equation \eqref{eq:cmplxeig} (only with \({z}_+\)), \[ A(\Re {z}_+ + i\Im {z}_+) =(\rho + i\omega) (\Re {z}_+ + i\Im {z}_+), \] or equivalently \[ A \Re {z}_+ = \rho \Re {z}_+ -\omega \Im {z}_+,\qquad A \Im {z}_+ = \rho \Im {z}_+ +\omega \Re {z}_+. \] To proceed, we take real and imaginary parts of the above eigenvector \({z}_\pm\) (remember that the real parts of \({z}_\pm\) are the same, and the imaginary parts only differ in their signs), forming the matrix \begin{equation*} U = \left[ \Re z \; , \; \Im z \right]^{-1} \qquad \mbox{or} \qquad U^{-1} = \left[ \Re z \; , \; \Im z \right]. \end{equation*} Then \begin{align*} A U^{-1} &= \left[ A \; \Re z , \; A \; \Im z \right] \\ &= \left[ \Re (A z )\; , \; \Im (A z ) \right], \qquad\qquad\qquad (\textrm{ as } A \text{ is real}) \\ &= \left[ \rho \Re ( z ) - \omega \Im ( z ) ,\ \rho \Im ( z ) + \omega \Im ( z ) \right] \cr &=[\Re {z}_+, \Im {z}_-] \begin{pmatrix} \rho & \omega \cr -\omega & \rho \end{pmatrix}. \end{align*} Thus we end up with \begin{equation*} \begin{pmatrix} \rho & \omega \\ - \omega & \rho \end{pmatrix} = UAU^{-1}, \end{equation*} where \(\begin{pmatrix} \rho & \omega \\ - \omega & \rho \end{pmatrix} \) is the complex normal form. In the new coordinates \(y=(u,v)^t=Ux\), the system becomes \(\dot{u} = \rho u + \omega v,\dot{v} = - \omega u + \rho v\). It is much easier to look at this system in the polar coordinates \(u = r \cos \theta, v=r \sin \theta\). Differentiating this new transform gives \begin{align*} \dot{u} & = \dot{r} \cos \theta - r \dot{\theta} \sin \theta = \rho r \cos \theta + \omega r \sin \theta \\ \dot{v} & = \dot{r} \sin \theta + r \dot{\theta} \cos \theta = - \omega r \cos \theta + \rho r \sin \theta . \end{align*} To eliminate \(\dot \theta\) to obtain an equation for \(\dot r\), multiply the first of these by \(\cos\theta\) and the second by \(\sin\theta\) and add to get \begin{equation}\label{requ} \dot{r} = \rho r, \quad \text{i.e.} \quad r = r_0 e^{\rho t}. \end{equation} Similarly to get the equation for \(\dot\theta\), multiply the first by \(\sin\theta\) and the second by \(\cos\theta\) and take the difference: \begin{equation*} \dot{\theta} = - \omega \quad\text{i.e.} \quad \theta =\theta_0 - \omega t, \end{equation*} which represents a clockwise rotation at constant angular velocity if \(\omega >0\). Using this to eliminate \(t\) from the equation for \(r\) shows that trajectories lie on spiral \(r=r_0 e^{\rho (\theta - \theta_0)/w}\).

b i) \(\rho>0\): unstable focus. In this case Equation (\ref{requ}) shows that solutions grow with time so trajectories spiral out of the origin. This is called a unstable focus clockwise if \(\omega >0\) counter-clockwise if \(\omega<0\). In the original coordinates, the phase portrait is a distorted spiral. To determine the direction of spiralling, we can consider the sign of \(\dot{x}_2\) on a horizontal line (where \(x_2=0\) through the stationary point or the sign of \(\dot{x}_1\) on a vertical line through the stationary point. If more detail is required the nullclines (see c(ii) below) indicate where solutions are flat or vertical as they move around the stationary point. This is called an unstable focus (see Figure 3.5).

Figure 3.5 Unstable focus in transformed \((u,v)\)-coordinates and in the original \((x_1,x_2)\)-coordinates.

b ii) \(\rho<0\) : stable focus In this case the \(\theta\( behaviour is the same as in the previous case but the radial velocity is towards zero. Solutions tend to the origin spiralling clockwise (if \(\omega>0\() as shown in Figure 3.6. In the original coordinates the solutions spiral inwards, with the direction given by consideration of the sign of \(\dot{x}_2\( (or \(\dot{x}_1\() on the horizontal line (resp. vertical line) through the stationary point. If more detail is required the nullclines (see c(ii) below) indicate where solutions are flat or vertical as they move around the stationary point. This is called a stable focus.

Figure 3.6 Stable focus in transformed \((u,v)\)-coordinates and in the original \((x_1,x_2)\)-coordinates.

b iii) \(\rho =0\): centre If \(\rho=0\) then \(\dot r=0\) and so \(r\) is constant -- solutions lie on circles in the transformed \((u,v)\)-coordinates with \(\theta\) changing at a constant rate, i.e. if \(r_0\ne 0\) then solutions are periodic with period \(\frac{2\pi}{|\omega |}\). This is called a centre, see Figure 3.7.

Figure 3.7 Centre in transformed \((u,v)\)-coordinates and in the original \((x_1,x_2)\)-coordinates.

Clearly if nonlinear terms are added then \(\dot r\) may no longer vanish, so we do not expect this type of behaviour to persist in typical nonlinear systems.

c) Repeated real roots \(\lambda\neq0.\) If the characteristic equation has two repeated roots \(\lambda=\lambda_1=\lambda_2\), then by Cayley-Hamiltonian Theorem (a matrix satisfies its own characteristic equation) \((A-\lambda I)^2=0\). Depending on the number of eigenvectors to the equation \((A-\lambda I){e}=0\), we have two cases (the equivalent two cases are \(A=\lambda I\) and \(A\neq \lambda I\).

c i) Repeated real roots \(\lambda\neq 0\): star Suppose there are two (linearly independent) eigenvectors \(e_1\) and \(e_2\) to the equation \((A-\lambda I){e} = {0}\), then \[ A[e_1, e_2] = [\lambda e_1, \lambda e_2] = \lambda [e_1, e_2]. \] Since \([e_1, e_2]\) is non-singular, we can right multiply \([e_1, e_2]^{-1}\) to the previous equation to get \begin{equation*} A = \lambda I = \begin{pmatrix} \lambda & 0 \\ 0 & \lambda \end{pmatrix}, \end{equation*} i.e. \(A\) is diagonal in any basis! Then using Equation (\ref{parabola}) with \(\lambda_1=\lambda_2\), \(\frac{x}{x_0} = \frac{y}{y_0}\) and so solutions lie on straight lines through the origin as shown in Figure 3.8. This is called a stable star if \(\lambda <0\) (so solutions tend to the origin) and an unstable star if \(\lambda >0\) (so solutions grow).

Figure 3.8 Stable stars and unstable stars

c ii): repeated roots \(\lambda\neq 0\): degenerate node. Suppose that there is only one eigenvalue \(e_1\) to the eigenvalue \(\lambda\) (although it is repeated), then we can find another vector \(e_2\) such that \( (A-\lambda I) e_2 = e_1\), or \(Ae_2 = \lambda e_2 + e_1\). Let \(U^{-1} = \left[ e_1; e_2 \right]\), the matrix with columns are the two vectors defined above. Then \begin{equation*} A U^{-1} = \left[ A e_1, A e_2 \right] = \left[ \lambda e_1, \lambda e_2 + e_1 \right] =[e_1,e_2]\begin{pmatrix} \lambda & 1\cr 0 & \lambda \end{pmatrix}. \end{equation*} Hence we get the \begin{equation*} \begin{pmatrix}\lambda & 1 \\ 0 & \lambda \end{pmatrix}= U A U^{-1}, \end{equation*} where the matrix \(\begin{pmatrix}\lambda & 1 \\ 0 & \lambda \end{pmatrix}\) is the normal form for this case of repeated roots. In the transformed coordinates \((u,v)\) defined as above, \begin{equation}\label{degnodeeqs} \dot{u} = \lambda u + v ,\quad \dot{v} = \lambda v. \end{equation} The second equation is easily solved to give \(v = v_0 e^{\lambda t}\). Substituting this into the first equation gives \(\dot u = \lambda u+v_0e^{\lambda t}\). Using the integrating factor \(e^{-\lambda t}\), we get \[ \frac{d}{dt} (u e^{- \lambda t}) = v_0. \] The integration of both sides lead to \(u e^{- \lambda t} - u_0 = v_0 t\) or \(u = u_0 e^{\lambda t} + v_0 t e^{ \lambda t}\). It is hard to analyse solutions directly from \eqref{degnodeeqs}. The phase portrait is given in Figure 3.9. They divide phase space into four regions according to the different combinations of the signs of \(\dot u\) and \(\dot v\). Bu considering the behaviour in each of these regions, we arrive at the phase portrait sketched. Of course, a rigorous justification takes more work, but this is enough to give a basic idea of the behaviour. The phase portrait of Figure 3.9 is an unstable degenerate node. For the case \(\lambda <0\) (a stable degenerate node), although the solution eventually converges to the origin, it may take a long excursion to finally move towards it.

Figure 3.9 Phase portrait for \(\dot{u}=\lambda u+v, \dot{v}=\lambda v\) with degenerate node (\(\lambda>0\)).

Direction of rotation for foci and centres. The sign of \(\omega\) in the normal form can always be chosen to be positive, but this might use a transformation that reverses the orientation of the plane, i.e. counter-clockwise rotation can be transformed into clockwise rotation. To determine the actual direction in an example either calculate the nullclines and the direction of the flow across the nullclines, or (and this is often easier) consider the direction of the flow on coordinate axes.

Example 3.3 Suppose that \( \dot x=-x-4y, \quad \dot y = x-y \) so the Jacobian matrix is \( \left(\begin{array}{cc} -1 & -4\\ 1& -1\end{array}\right) \) with characteristic equation \((s { + }1)^2+4=0\) so the eigenvalues are \(-1\pm i2\). The origin is therefore a stable node, but in which way does it rotate? Set \(x=0\) (the \(y\)-axis) and consider how solutions move across this line. On \(x=0\), \(\dot x=-4y\) and so if \(y>0\) then the motion is from right to left (as \(\dot x<0\) and if \(y<0\) (so \(\dot x>0\) the motion is from left to right. Thus the motion is counter-clockwise. It is often useful to indicate this on a diagram to make sure the figure is drawn appropriately.

Summary on the relation between the signs of eigenvalues and the behaviour of the underlying linear system. Given \(\dot{x} = Ax\), \(x \in \mathbb{R}^2\). Find the eigenvalues of \(A\) to characterise all the behaviours of the solution summarised as below.

  1. real, distinct:
  2. complex conjugate pair:
  3. other special cases:

For 2D systems, since the behaviours of the solutions, or equivalently of the roots depend only on the determinant and the trace of \(A\), they can be equally summarised as in Figure 3.10, where the parabola is the curve \(4\mbox{det} A = (\mbox{tr}A)^2\).

Figure 3.10 Different behaviours of the system \(\dot{x}=Ax\), depending on the determinant and the trace of \(A\).
The linearised systems for higher dimensions can be analysed in a similar way.