This note follows on from the one devoted to first-order differential equations by considering several variables whose rates of change depend linearly on the levels of all the others. The new difficulty lies in this entanglement, and the strategy is to undo it by a change of basis: in a basis of eigenvectors of the transition matrix, the system reduces to scalar equations we know how to solve, at the cost of a detour through the complex numbers when the eigenvalues are not real. We then deal with non-diagonalisable matrices, with linear equations of order \(n\), which are systems in disguise, and we end with stability, with the theorem characterising the matrices whose solutions all converge to zero and the Routh–Hurwitz criterion. Numerical checks and figures are done in Python, and the solutions to the exercises are given in collapsible blocks.
First-order linear systems
Let \(a_{i,j}\) (for \(i=1,\dots,n\) and \(j=1,\dots,n\)) be real parameters and \(b_i(t)\) (for \(i=1,\dots,n\)) continuous functions from \(\mathbb R^+\) to \(\mathbb R\). The functions \(x_i(t)\), \(i=1,\dots,n\), defined on \(\mathbb R^+\) are solutions of the system of first-order linear differential equations:
\begin{equation*} \begin{cases} \dot x_1(t) = a_{1,1}\,x_1(t) + a_{1,2}\,x_2(t) + \dots + a_{1,n}\,x_n(t) + b_1(t)\\ \dot x_2(t) = a_{2,1}\,x_1(t) + a_{2,2}\,x_2(t) + \dots + a_{2,n}\,x_n(t) + b_2(t)\\ \quad\vdots\\ \dot x_n(t) = a_{n,1}\,x_1(t) + a_{n,2}\,x_2(t) + \dots + a_{n,n}\,x_n(t) + b_n(t) \end{cases} \end{equation*}if they are differentiable on \(\mathbb R^+\) and if, together with their derivatives \(\dot x_i(t)\), they satisfy these \(n\) equations for all \(t\) in \(\mathbb R^+\).
The system is linear because the levels \(x_j(t)\) only enter through linear combinations, and first-order because only first derivatives appear. The functions \(b_i(t)\) form the right-hand side; when they are all zero the system is said to be homogeneous.
Matrix notation
Let \(X(t)\equiv\bigl(x_1(t),\dots,x_n(t)\bigr)^\top\) denote the vector of variables, with values in \(\mathbb R^n\), \(\dot X(t)\) the vector of their derivatives and \(b(t)\equiv\bigl(b_1(t),\dots,b_n(t)\bigr)^\top\) the right-hand side vector. The system of the definition can then be written equivalently as:
\begin{equation*} \dot X(t) = A\,X(t) + b(t) \end{equation*}where the real square matrix \(A\), of dimension \(n\times n\), is defined by:
\begin{equation*} A = \begin{pmatrix} a_{1,1} & a_{1,2} & \dots & a_{1,n}\\ a_{2,1} & a_{2,2} & \dots & a_{2,n}\\ \vdots & \vdots & \ddots & \vdots\\ a_{n,1} & a_{n,2} & \dots & a_{n,n} \end{pmatrix} \end{equation*}We shall call \(A\) the transition matrix of the system. When \(n=1\) we recover the equation \(\dot x(t) = a\,x(t) + b(t)\) studied in the previous note, whose results we shall use freely.
We are interested in the real solutions of the system. Nevertheless, the search for these solutions will lead us to transform the system, and we shall sometimes have to compute complex solutions of the transformed system, that is, functions from \(\mathbb R^+\) to \(\mathbb C^n\). The following property shows that this detour raises no difficulty.
If \(A\) and \(b(t)\) are a real matrix and a real vector, then the complex function \(X(t) = X_1(t) + \mathrm i\,X_2(t)\), where \(X_1\) and \(X_2\) take values in \(\mathbb R^n\), is a solution of \(\dot X(t) = A\,X(t) + b(t)\) if and only if:
\begin{equation*} \begin{cases} \dot X_1(t) = A\,X_1(t) + b(t)\\ \dot X_2(t) = A\,X_2(t) \end{cases} \end{equation*}The real part of a complex solution is a real solution of the complete system, and its imaginary part a real solution of the homogeneous system.
Substituting \(X(t) = X_1(t) + \mathrm i\,X_2(t)\) into the system gives \(\dot X_1(t) + \mathrm i\,\dot X_2(t) = A\,X_1(t) + \mathrm i\,A\,X_2(t) + b(t)\). Since \(A\) and \(b(t)\) are real, \(A X_1(t) + b(t)\) and \(A X_2(t)\) are real, and two complex vectors are equal if and only if their real and imaginary parts are.
An equivalent transformation of the system
Solving the system directly is generally not possible, since the rate of change of each variable, \(\dot x_i(t)\), depends on the levels of all the variables. The idea is to consider new variables, independent linear combinations of the original ones, chosen so that the rate of change of each new variable depends only on its own level (and on the right-hand side). If we succeed, the system reduces to \(n\) scalar differential equations that we know how to solve.
More formally, let \(H\) be an invertible square matrix of dimension \(n\times n\), with entries in \(\mathbb C\), which we interpret as a change-of-basis matrix. Let:
\begin{equation*} Y(t) = H^{-1}X(t) \end{equation*}The vector \(Y(t)\) consists of \(n\) functions from \(\mathbb R^+\) to \(\mathbb C\), and by linearity of differentiation we have \(\dot Y(t) = H^{-1}\dot X(t)\). Substituting \(X(t) = H\,Y(t)\) into the system, then premultiplying by \(H^{-1}\), we get:
\begin{equation*} \begin{split} \dot X(t) = A\,X(t) + b(t) &\Leftrightarrow H\,\dot Y(t) = A\,H\,Y(t) + b(t)\\ &\Leftrightarrow \dot Y(t) = H^{-1}A\,H\,Y(t) + H^{-1}b(t) \end{split} \end{equation*}that is, writing \(\hat A \equiv H^{-1}AH\) and \(c(t)\equiv H^{-1}b(t)\):
\begin{equation*} \dot Y(t) = \hat A\,Y(t) + c(t) \end{equation*}This transformed system is equivalent to the original one: \(X(t)\) is a solution of the former if and only if \(Y(t) = H^{-1}X(t)\) is a solution of the latter, and the initial conditions correspond through \(Y(0) = H^{-1}X(0)\). The whole problem is to choose \(H\) so that \(\hat A\) is as simple as possible. The ideal case is when \(\hat A\) is diagonal; failing that, we shall see that one can always obtain a block triangular matrix, which is enough to solve the system step by step.
Reduction of the transition matrix
Diagonalisation
Let \(A\) be a real square matrix of dimension \(n\times n\) and \(I_n\) the identity matrix of the same dimension. The characteristic polynomial of \(A\) is the polynomial of degree \(n\) in \(\lambda\):
\begin{equation*} \chi_A(\lambda) = \left|A - \lambda I_n\right| \end{equation*}where \(|\cdot|\) denotes the determinant. The eigenvalues of \(A\) are the roots, real or complex, of \(\chi_A\), that is, the numbers \(\lambda\) such that \(|A-\lambda I_n| = 0\). The characteristic polynomial has at most \(n\) distinct roots. We denote by \(\lambda_k\), \(k=1,\dots,m\) (with \(m\leq n\)), these distinct roots and by \(r_k\) their multiplicities, with \(\sum_{k=1}^m r_k = n\). For each eigenvalue \(\lambda_k\), the matrix \(A-\lambda_kI_n\) is singular, so there exist non-zero vectors \(v\) such that:
\begin{equation*} (A-\lambda_kI_n)\,v = 0 \qquad\Leftrightarrow\qquad A\,v = \lambda_k v \end{equation*}Such a vector is an eigenvector associated with \(\lambda_k\). The set of solutions of this linear equation is a vector subspace of \(\mathbb C^n\), the eigenspace associated with \(\lambda_k\), whose dimension \(q_k\) satisfies \(1\leq q_k\leq r_k\).
The matrix \(A\) is diagonalisable if, for each eigenvalue \(\lambda_k\), one can find \(r_k\) linearly independent eigenvectors, that is, if \(q_k = r_k\) for all \(k=1,\dots,m\).
In this case, the \(n\) eigenvectors thus obtained are linearly independent and form the columns of an invertible change-of-basis matrix \(H\equiv(v_1,\dots,v_n)\). Since \(A\,v_i = \lambda_i v_i\) for each column, we have \(A\,H = H\,\hat A\) with:
\begin{equation*} \hat A = H^{-1}AH = \begin{pmatrix} \lambda_1 & 0 & \dots & 0\\ 0 & \lambda_2 & \ddots & \vdots\\ \vdots & \ddots & \ddots & 0\\ 0 & \dots & 0 & \lambda_n \end{pmatrix} \end{equation*}where each eigenvalue appears on the diagonal as many times as its multiplicity.
If the \(n\) eigenvalues of \(A\) are distinct, then \(A\) is diagonalisable.
Each eigenvalue has multiplicity \(r_k=1\) and has at least one eigenvector, so \(q_k = r_k\). It remains to check that eigenvectors associated with distinct eigenvalues are linearly independent. Suppose that \(v_1,\dots,v_p\), associated with distinct eigenvalues \(\lambda_1,\dots,\lambda_p\), satisfy \(\alpha_1v_1+\dots+\alpha_pv_p = 0\). Applying \(A-\lambda_pI_n\) to both sides, the last term vanishes and we get \(\alpha_1(\lambda_1-\lambda_p)v_1+\dots+\alpha_{p-1}(\lambda_{p-1}-\lambda_p)v_{p-1}=0\). Repeating the operation with \(\lambda_{p-1},\dots,\lambda_2\), only \(\alpha_1\prod_{j=2}^p(\lambda_1-\lambda_j)v_1 = 0\) remains, hence \(\alpha_1=0\) since the eigenvalues are distinct and \(v_1\neq0\). The same argument applied to each index shows that all the \(\alpha_i\) are zero.
If \(A\) is real, its characteristic polynomial has real coefficients and its complex roots come in conjugate pairs. If \(v\) is an eigenvector associated with \(\lambda\), then \(\bar v\) is an eigenvector associated with \(\bar\lambda\), as one sees by conjugating \(Av=\lambda v\). A real symmetric matrix has real eigenvalues and is always diagonalisable, with orthogonal eigenvectors.
Consider the real symmetric matrix:
\begin{equation*} A = \begin{pmatrix} 2 & 1\\ 1 & 2\end{pmatrix} \end{equation*}Its characteristic polynomial is:
\begin{equation*} \begin{split} \chi_A(\lambda) &= \begin{vmatrix} 2-\lambda & 1\\ 1 & 2-\lambda\end{vmatrix}\\ &= (2-\lambda)^2 - 1\\ &= (1-\lambda)(3-\lambda) \end{split} \end{equation*}The eigenvalues are therefore real and distinct: \(\lambda_1=1\) and \(\lambda_2=3\). The eigenvector \(v_1\) associated with \(\lambda_1\) must satisfy \((A-\lambda_1I_2)\,v_1 = 0\), that is:
\begin{equation*} \begin{pmatrix} 1 & 1\\ 1 & 1\end{pmatrix} \begin{pmatrix} v_{1,1}\\ v_{2,1}\end{pmatrix} = 0 \quad\Leftrightarrow\quad v_{2,1} = -v_{1,1} \end{equation*}An eigenvector is only defined up to a multiplicative factor; we set \(v_1 = (1,-1)^\top\). Likewise, the eigenvector \(v_2\) associated with \(\lambda_2\) satisfies \((A-\lambda_2I_2)\,v_2 = 0\):
\begin{equation*} \begin{pmatrix} -1 & 1\\ 1 & -1\end{pmatrix} \begin{pmatrix} v_{1,2}\\ v_{2,2}\end{pmatrix} = 0 \quad\Leftrightarrow\quad v_{2,2} = v_{1,2} \end{equation*}and we set \(v_2 = (1,1)^\top\). The change-of-basis matrix and its inverse are:
\begin{equation*} H = \begin{pmatrix} 1 & 1\\ -1 & 1\end{pmatrix} \qquad\text{and}\qquad H^{-1} = \frac{1}{2}\begin{pmatrix} 1 & -1\\ 1 & 1\end{pmatrix} \end{equation*}and one checks that:
\begin{equation*} H^{-1}AH = \begin{pmatrix} 1 & 0\\ 0 & 3\end{pmatrix} \end{equation*}
These computations can be checked numerically. The eig function of
numpy returns the eigenvalues and a matrix whose columns are normalised
eigenvectors (of unit norm), which changes nothing to the diagonalisation
since eigenvectors are only defined up to a factor.
A = np.array([[2, 1], [1, 2]])
lam, V = np.linalg.eig(A)
print(lam) # 3 and 1
print(V) # columns proportional to (1,1) and (1,-1)
H = np.array([[1, 1], [-1, 1]])
print(np.linalg.inv(H) @ A @ H) # diag(1, 3)
Reduced form of a non-diagonalisable matrix
What can be done when the transition matrix is not diagonalisable, that is, when for some eigenvalue \(\lambda_k\) of multiplicity \(r_k\) one only finds \(q_k < r_k\) linearly independent eigenvectors? We can no longer make \(\hat A\) diagonal, but we can make it block triangular, which will be enough to solve the system recursively.
Let \(A\) be a square matrix of dimension \(n\times n\), with real or complex entries, whose distinct eigenvalues are \(\lambda_1,\dots,\lambda_m\) with multiplicities \(r_1,\dots,r_m\). There exists an invertible change-of-basis matrix \(H\) such that \(\hat A = H^{-1}AH\) is block diagonal, with one block per distinct eigenvalue. The block associated with \(\lambda_k\) is an upper triangular matrix of dimension \(r_k\times r_k\) whose diagonal entries are all equal to \(\lambda_k\).
We take this result for granted1 and describe the construction of the columns of \(H\) associated with an eigenvalue \(\lambda_k\), writing \(N \equiv A-\lambda_kI_n\). We proceed in steps.
- Choose a basis \(v_1,\dots,v_{q_k}\) of the eigenspace, that is, \(q_k\) non-zero linearly independent vectors such that \(N v_i = 0\).
If \(q_k < r_k\), look for a vector \(w\), not belonging to the subspace spanned by the vectors already obtained, such that \(Nw\) belongs to this subspace:
\begin{equation*} (A-\lambda_kI_n)\,w = \alpha_1v_1+\dots+\alpha_{q_k}v_{q_k} \end{equation*}for coefficients \(\alpha_i\), not all zero, determined together with \(w\) (the linear system in \(w\) only has a solution for certain values of the \(\alpha_i\)). Such a \(w\) is called a generalised eigenvector.
- Add \(w\) to the list of vectors and, as long as \(r_k\) linearly independent vectors have not been obtained, repeat step 2, looking for a new vector whose image under \(N\) belongs to the subspace spanned by all the vectors already obtained.
One shows that this construction always yields \(r_k\) linearly independent vectors2. If these vectors are placed, in the order in which they were obtained, in the columns of \(H\), the corresponding block of \(H^{-1}AH\) is upper triangular: indeed \(A\,v_i = \lambda_kv_i\) for the eigenvectors, and for each generalised vector \(A\,w = \lambda_kw + \alpha_1v_1+\dots\) only involves \(w\) and vectors placed before it in \(H\).
Consider the real square matrix:
\begin{equation*} A = \begin{pmatrix} 3 & 3 & 1\\ -2 & -2 & -1\\ 2 & 3 & 2\end{pmatrix} \end{equation*}Its characteristic polynomial is \(\chi_A(\lambda) = (1-\lambda)^3\), which has a triple root \(\lambda=1\). The eigenvectors are the solutions of the system \((A-I_3)\,v = 0\):
\begin{equation*} \begin{cases} 2v_1 + 3v_2 + v_3 = 0\\ -2v_1 - 3v_2 - v_3 = 0\\ 2v_1 + 3v_2 + v_3 = 0 \end{cases} \quad\Leftrightarrow\quad 2v_1 + 3v_2 + v_3 = 0 \end{equation*}which has rank one. The eigenspace therefore has dimension two, there are only two linearly independent eigenvectors and the matrix \(A\) is not diagonalisable. We may choose:
\begin{equation*} v_1 = \begin{pmatrix} 0\\ 1\\ -3\end{pmatrix} \qquad\text{and}\qquad v_2 = \begin{pmatrix} 1\\ 0\\ -2\end{pmatrix} \end{equation*}The system \((A-I_3)\,w = \alpha_1v_1+\alpha_2v_2\) reads:
\begin{equation*} \begin{cases} 2w_1 + 3w_2 + w_3 = \alpha_2\\ -2w_1 - 3w_2 - w_3 = \alpha_1\\ 2w_1 + 3w_2 + w_3 = -3\alpha_1 - 2\alpha_2 \end{cases} \end{equation*}The first two equations are only compatible if \(\alpha_2 = -\alpha_1\), and the third is then identical to the first. Setting \(\alpha_1=1\) and \(\alpha_2 = -1\), we are left with \(2w_1+3w_2+w_3 = -1\) and may choose \(w = (0,0,-1)^\top\), which is not a linear combination of \(v_1\) and \(v_2\). The change-of-basis matrix and its inverse are:
\begin{equation*} H = \begin{pmatrix} 0 & 1 & 0\\ 1 & 0 & 0\\ -3 & -2 & -1\end{pmatrix} \qquad\text{and}\qquad H^{-1} = \begin{pmatrix} 0 & 1 & 0\\ 1 & 0 & 0\\ -2 & -3 & -1\end{pmatrix} \end{equation*}and we obtain:
\begin{equation*} \hat A = H^{-1}AH = \begin{pmatrix} 1 & 0 & 1\\ 0 & 1 & -1\\ 0 & 0 & 1\end{pmatrix} \end{equation*}The last column of \(\hat A\) is none other than the vector \((\alpha_1,\alpha_2,1)^\top\), in accordance with \(A\,w = w + v_1 - v_2\).
Numerically, the rank of \(A-I_3\) confirms that there are only two independent eigenvectors, and we recover the reduced form:
A = np.array([[3, 3, 1], [-2, -2, -1], [2, 3, 2]])
print(np.linalg.eigvals(A)) # 1, 1, 1 (up to rounding errors)
print(np.linalg.matrix_rank(A - np.eye(3))) # 1
H = np.array([[0, 1, 0], [1, 0, 0], [-3, -2, -1]])
print(np.linalg.inv(H) @ A @ H) # [[1, 0, 1], [0, 1, -1], [0, 0, 1]]
Solving the system
The transformed system
Let us return to the transformed system \(\dot Y(t) = \hat A\,Y(t) + c(t)\), where \(\hat A = H^{-1}AH\) is a diagonal or reduced form of \(A\). The entries of \(\hat A\) and the functions \(c_i(t)\) are complex in general. We shall need the solution of the scalar first-order equation with a complex coefficient.
Let \(\lambda\in\mathbb C\) and \(c(t)\) be a continuous function from \(\mathbb R^+\) to \(\mathbb C\). For any initial condition \(y(0) = y_0\in\mathbb C\), the differential equation \(\dot y(t) = \lambda\,y(t) + c(t)\) has a unique solution on \(\mathbb R^+\):
\begin{equation*} y(t) = e^{\lambda t}\left(y_0 + \int_0^t e^{-\lambda s}c(s)\,\mathrm ds\right) \end{equation*}This is the formula established for real coefficients in the previous note, and the proof is the same: the derivative of \(e^{-\lambda t}\) is \(-\lambda e^{-\lambda t}\) for any complex \(\lambda\), so that the equation is equivalent to \(\frac{\mathrm d}{\mathrm dt}\bigl(e^{-\lambda t}y(t)\bigr) = e^{-\lambda t}c(t)\), which we integrate between \(0\) and \(t\).
Two cases arise.
- If \(\hat A\) is diagonal, the transformed system reduces to \(n\) independent equations \(\dot y_i(t) = \lambda_iy_i(t) + c_i(t)\), which the previous property solves one at a time.
If \(\hat A\) is block diagonal, each block can be treated separately. Consider the block associated with the eigenvalue \(\lambda_k\) of multiplicity \(r_k\), and suppose it occupies rows \(i+1\) to \(i+r_k\) of \(\hat A\). The corresponding variables satisfy a subsystem of \(r_k\) equations, independent of the other \(n-r_k\) variables:
\begin{equation*} \begin{pmatrix} \dot y_{i+1}(t)\\ \dot y_{i+2}(t)\\ \vdots\\ \dot y_{i+r_k}(t)\end{pmatrix} = \begin{pmatrix} \lambda_k & \hat a_{k,1,2} & \dots & \hat a_{k,1,r_k}\\ 0 & \lambda_k & \dots & \hat a_{k,2,r_k}\\ \vdots & \ddots & \ddots & \vdots\\ 0 & \dots & 0 & \lambda_k \end{pmatrix} \begin{pmatrix} y_{i+1}(t)\\ y_{i+2}(t)\\ \vdots\\ y_{i+r_k}(t)\end{pmatrix} + \begin{pmatrix} c_{i+1}(t)\\ c_{i+2}(t)\\ \vdots\\ c_{i+r_k}(t)\end{pmatrix} \end{equation*}The block being upper triangular, it suffices to solve from the bottom up. The last equation, \(\dot y_{i+r_k}(t) = \lambda_ky_{i+r_k}(t) + c_{i+r_k}(t)\), is a scalar equation we know how to solve: given the initial condition \(y_{i+r_k}(0)\), it has a unique solution. Once this solution is obtained, we substitute it into the equation above, which in turn becomes a scalar first-order equation in \(y_{i+r_k-1}(t)\) with a known right-hand side. Climbing the triangle in this way, we obtain a unique solution for the \(r_k\) variables of the block, given their initial conditions.
In all cases, the transformed system has, for any initial vector \(Y_0\in\mathbb C^n\), one and only one solution \(Y(t)\) from \(\mathbb R^+\) to \(\mathbb C^n\) such that \(Y(0) = Y_0\).
Existence and uniqueness
If the functions \(b_i(t)\) of the right-hand side are continuous on \(\mathbb R^+\), then for any vector \(X_0\in\mathbb R^n\) the system \(\dot X(t) = A\,X(t) + b(t)\) has one and only one real solution on \(\mathbb R^+\) such that \(X(0) = X_0\). It is also the unique complex solution satisfying this initial condition, and it is equal to \(H\,Y(t)\), where \(Y(t)\) is the solution of the transformed system such that \(Y(0) = H^{-1}X_0\).
We have seen that the transformed system has one and only one complex solution \(Y(t)\) such that \(Y(0) = H^{-1}X_0\). Since the two systems are equivalent, \(X(t) = H\,Y(t)\) is the unique complex solution of the original system such that \(X(0) = X_0\). Write this solution as \(X(t) = X_1(t) + \mathrm i\,X_2(t)\). By the property established in the first section, the imaginary part satisfies \(\dot X_2(t) = A\,X_2(t)\) with, since \(X_0\) is real, \(X_2(0) = 0\). Now the zero function is a solution of this homogeneous system with this initial condition, and this solution is unique: \(X_2(t)\) is therefore zero at all times, and the unique complex solution is real.
An example with complex eigenvalues
We want to solve the system of differential equations:
\begin{equation*} \begin{cases} \dot x_1(t) = x_1(t) - x_2(t) + 2\\ \dot x_2(t) = x_1(t) + x_2(t) - 2t \end{cases} \end{equation*}The equivalent matrix form is \(\dot X(t) = A\,X(t) + b(t)\) with:
\begin{equation*} A = \begin{pmatrix} 1 & -1\\ 1 & 1\end{pmatrix} \qquad\text{and}\qquad b(t) = \begin{pmatrix} 2\\ -2t\end{pmatrix} \end{equation*}The characteristic polynomial is:
\begin{equation*} \chi_A(\lambda) = \begin{vmatrix} 1-\lambda & -1\\ 1 & 1-\lambda\end{vmatrix} = (1-\lambda)^2 + 1 \end{equation*}whose roots are \(\lambda_1 = 1-\mathrm i\) and \(\lambda_2 = 1+\mathrm i\). The eigenvalues are distinct, so the matrix is diagonalisable over \(\mathbb C\). The eigenvector \(v_1\) associated with \(\lambda_1\) satisfies:
\begin{equation*} \begin{pmatrix} \mathrm i & -1\\ 1 & \mathrm i\end{pmatrix} \begin{pmatrix} v_{1,1}\\ v_{2,1}\end{pmatrix} = 0 \quad\Leftrightarrow\quad \begin{cases} \mathrm i\,v_{1,1} - v_{2,1} = 0\\ v_{1,1} + \mathrm i\,v_{2,1} = 0 \end{cases} \end{equation*}The second equation is \(-\mathrm i\) times the first, so only one condition remains, \(v_{2,1} = \mathrm i\,v_{1,1}\), and we set \(v_1 = (-\mathrm i, 1)^\top\). The eigenvector associated with \(\lambda_2\) is the conjugate, \(v_2 = (\mathrm i, 1)^\top\), as one checks directly:
\begin{equation*} \begin{pmatrix} -\mathrm i & -1\\ 1 & -\mathrm i\end{pmatrix} \begin{pmatrix} \mathrm i\\ 1\end{pmatrix} = \begin{pmatrix} 1 - 1\\ \mathrm i-\mathrm i\end{pmatrix} = 0 \end{equation*}The change-of-basis matrix and its inverse are:
\begin{equation*} H = \begin{pmatrix} -\mathrm i & \mathrm i\\ 1 & 1\end{pmatrix} \qquad\text{and}\qquad H^{-1} = \frac{1}{2}\begin{pmatrix} \mathrm i & 1\\ -\mathrm i & 1\end{pmatrix} \end{equation*}Setting \(Y(t) = H^{-1}X(t)\) and \(c(t) = H^{-1}b(t)\) we get \(\dot Y(t) = \hat A\,Y(t) + c(t)\) with:
\begin{equation*} \hat A = \begin{pmatrix} 1-\mathrm i & 0\\ 0 & 1+\mathrm i\end{pmatrix} \qquad\text{and}\qquad c(t) = \begin{pmatrix} \mathrm i - t\\ -\mathrm i - t\end{pmatrix} \end{equation*}that is, equivalently:
\begin{equation*} \begin{cases} \dot y_1(t) = (1-\mathrm i)\,y_1(t) + \mathrm i - t\\ \dot y_2(t) = (1+\mathrm i)\,y_2(t) - \mathrm i - t \end{cases} \end{equation*}These two complex linear differential equations are independent. Rather than applying the integral formula, we proceed as in the real case: we look for a particular solution of each complete equation in the form of a polynomial \(\alpha + \beta t\), whose complex coefficients are identified by substitution into the equation, and then add the general solution of the homogeneous equation. For the first equation, \(\beta = (1-\mathrm i)(\alpha+\beta t) + \mathrm i - t\) for all \(t\) requires \((1-\mathrm i)\beta = 1\), that is \(\beta = \frac{1+\mathrm i}{2}\), and then \(\beta = (1-\mathrm i)\alpha + \mathrm i\), that is \(\alpha = \frac12\). The computation is identical for the second equation, and we obtain:
\begin{equation*} \begin{cases} y_1(t) = \gamma_1e^{(1-\mathrm i)t} + \frac{1+\mathrm i}{2}\,t + \frac12\\ y_2(t) = \gamma_2e^{(1+\mathrm i)t} + \frac{1-\mathrm i}{2}\,t + \frac12 \end{cases} \end{equation*}where the constants of integration are determined by \(Y(0) = H^{-1}X(0)\):
\begin{equation*} \begin{cases} \gamma_1 + \frac12 = y_1(0) = \frac{\mathrm i}{2}x_{0,1} + \frac12x_{0,2}\\ \gamma_2 + \frac12 = y_2(0) = -\frac{\mathrm i}{2}x_{0,1} + \frac12x_{0,2} \end{cases} \end{equation*}with \(X(0) = (x_{0,1},x_{0,2})^\top\). Note that \(y_2(t) = \overline{y_1(t)}\), which is no coincidence: the solution \(X(t)\) being real, the two coordinates of \(Y(t) = H^{-1}X(t)\) are conjugate since the rows of \(H^{-1}\) are. Returning to the original problem, the solution is \(X(t) = H\,Y(t)\), that is:
\begin{equation*} \begin{cases} x_1(t) = -\mathrm i\,y_1(t) + \mathrm i\,y_2(t)\\ x_2(t) = y_1(t) + y_2(t) \end{cases} \end{equation*}Substituting \(e^{(1\pm\mathrm i)t} = e^t(\cos t \pm \mathrm i\sin t)\) and \(\gamma_1 = \frac12\bigl(\mathrm i\,x_{0,1} + x_{0,2} - 1\bigr)\), \(\gamma_2 = \overline{\gamma_1}\), we obtain after simplification:
\begin{equation*} \begin{cases} x_1(t) = x_{0,1}\,e^t\cos t - (x_{0,2}-1)\,e^t\sin t + t\\ x_2(t) = x_{0,1}\,e^t\sin t + (x_{0,2}-1)\,e^t\cos t + t + 1 \end{cases} \end{equation*}One easily checks that these functions satisfy the system and the initial condition.
The particular solution \((t, t+1)^\top\) is polynomial, like the right-hand
side, and the solution of the homogeneous system is a combination of
\(e^t\cos t\) and \(e^t\sin t\): the real part of the eigenvalues, \(1\), governs
the amplitude and the imaginary part, \(\pm1\), the frequency of the
oscillations. The modulus of the homogeneous solution grows without bound,
the dynamics are unstable and the trajectories spiral away from the origin.
The figure below shows a few trajectories of the homogeneous
system in the \((x_1,x_2)\) plane, as well as the solution of the complete
system, on which a numerical solution obtained with solve_ivp has been
superimposed.
A = np.array([[1, -1], [1, 1]])
b = lambda t: np.array([2, -2*t])
x0 = np.array([1.0, 0.5])
t = np.linspace(0, 3, 301)
numerical = solve_ivp(lambda t, X: A @ X + b(t), (0, 3), x0, t_eval=t, rtol=1e-10, atol=1e-12)
x1 = x0[0]*np.exp(t)*np.cos(t) - (x0[1]-1)*np.exp(t)*np.sin(t) + t
x2 = x0[0]*np.exp(t)*np.sin(t) + (x0[1]-1)*np.exp(t)*np.cos(t) + t + 1
print(np.max(np.abs(numerical.y[0] - x1)), np.max(np.abs(numerical.y[1] - x2)))
Solve the system of differential equations:
\begin{equation*} \begin{cases} \dot x(t) = 3x(t) - 2y(t) + 1\\ \dot y(t) = x(t) + 1 \end{cases} \end{equation*}by writing it in matrix form, then transforming it in a basis of eigenvectors. Characterise the steady state and its stability.
Solution
Setting \(X(t) = (x(t),y(t))^\top\), the system reads \(\dot X(t) = A\,X(t) + b\) with:
\begin{equation*} A = \begin{pmatrix} 3 & -2\\ 1 & 0\end{pmatrix} \qquad\text{and}\qquad b = \begin{pmatrix} 1\\ 1\end{pmatrix} \end{equation*}The characteristic polynomial is \(\chi_A(\lambda) = (3-\lambda)(-\lambda) + 2 = \lambda^2 - 3\lambda + 2 = (\lambda-1)(\lambda-2)\), whose roots \(\lambda_1 = 1\) and \(\lambda_2 = 2\) are real and distinct. The eigenvectors satisfy \((A-I_2)\,v_1 = 0 \Leftrightarrow 2v_{1,1} - 2v_{2,1} = 0\) and \((A-2I_2)\,v_2 = 0 \Leftrightarrow v_{1,2} - 2v_{2,2} = 0\); we set \(v_1 = (1,1)^\top\) and \(v_2 = (2,1)^\top\), hence:
\begin{equation*} H = \begin{pmatrix} 1 & 2\\ 1 & 1\end{pmatrix} \qquad H^{-1} = \begin{pmatrix} -1 & 2\\ 1 & -1\end{pmatrix} \qquad \hat A = \begin{pmatrix} 1 & 0\\ 0 & 2\end{pmatrix} \qquad c = H^{-1}b = \begin{pmatrix} 1\\ 0\end{pmatrix} \end{equation*}The transformed system, \(\dot y_1(t) = y_1(t) + 1\) and \(\dot y_2(t) = 2y_2(t)\), has solution \(y_1(t) = \bigl(y_1(0)+1\bigr)e^t - 1\) and \(y_2(t) = y_2(0)\,e^{2t}\). Returning to the original variables, \(X(t) = H\,Y(t)\):
\begin{equation*} \begin{cases} x(t) = \bigl(y_1(0)+1\bigr)e^t + 2y_2(0)\,e^{2t} - 1\\ y(t) = \bigl(y_1(0)+1\bigr)e^t + y_2(0)\,e^{2t} - 1 \end{cases} \end{equation*}with \(y_1(0) = -x_0 + 2y_0\) and \(y_2(0) = x_0 - y_0\), where \((x_0,y_0)\) is the initial condition. The steady state is the constant solution \(X^{\star} = -A^{-1}b = (-1,-1)^\top\), and the solution can also be written:
\begin{equation*} X(t) = X^{\star} + (2y_0 - x_0 + 1)\,e^{t}\begin{pmatrix} 1\\ 1\end{pmatrix} + (x_0-y_0)\,e^{2t}\begin{pmatrix} 2\\ 1\end{pmatrix} \end{equation*}Both eigenvalues being strictly positive, the deviation from the steady state grows exponentially for any initial condition other than \(X^{\star}\): the steady state is an unstable node. In the phase plane, the trajectories leave it tangentially to the direction \(v_1\), associated with the smaller eigenvalue, before turning towards the direction \(v_2\), whose weight \(e^{2t}\) eventually dominates, as the figure below shows.
Characterisation of the solutions
The following two properties extend to the vector case the structure of the solutions of a linear equation: the general solution of the complete equation is the sum of a particular solution and of the general solution of the homogeneous equation, and the latter is a vector space.
If \(\Psi(t)\) is a particular solution of the system \(\dot X(t) = A\,X(t) + b(t)\), then \(X(t)\) is a solution of this system if and only if \(X(t) - \Psi(t)\) is a solution of the homogeneous system \(\dot Z(t) = A\,Z(t)\).
Since \(\Psi(t)\) is a solution of the complete equation, we have \(b(t) = \dot\Psi(t) - A\,\Psi(t)\). The function \(X(t)\) is a solution if and only if \(b(t) = \dot X(t) - A\,X(t)\), that is, if and only if \(\dot X(t) - \dot\Psi(t) = A\bigl(X(t) - \Psi(t)\bigr)\).
The set of real (respectively complex) solutions of the homogeneous linear system \(\dot X(t) = A\,X(t)\), \(t\in\mathbb R^+\), is a vector space over \(\mathbb R\) (respectively \(\mathbb C\)) of dimension \(n\).
By linearity, any linear combination of solutions of the homogeneous system is again a solution: the set of solutions is a vector space. For \(i=1,\dots,n\), let \(X_i(t)\) be the unique solution such that \(X_i(0) = U_i\), where \(U_i\) is the \(i\)-th vector of the canonical basis (whose \(j\)-th coordinate is \(1\) if \(j=i\) and \(0\) otherwise). The \(n\) functions \(X_i(t)\) are linearly independent, since their values at \(t=0\) are linearly independent vectors. Now let \(X_0 = (x_{1,0},\dots,x_{n,0})^\top\) be any initial condition. By linearity, \(\sum_{i=1}^nx_{i,0}X_i(t)\) is a solution of the homogeneous system, and it equals \(X_0\) at \(t=0\); by uniqueness, it is the unique solution such that \(X(0) = X_0\):
\begin{equation*} X(t) = \sum_{i=1}^n x_{i,0}\,X_i(t) \end{equation*}Every solution is therefore a linear combination of the \(X_i(t)\), which form a basis of the space of solutions3.
Let \(H = (v_1,\dots,v_n)\) be a change-of-basis matrix giving a reduced (or diagonal) form of the transition matrix \(A\). The general solution of the homogeneous system \(\dot X(t) = A\,X(t)\) is:
\begin{equation*} X(t) = \sum_{i=1}^n y_i(t)\,v_i \qquad t\in\mathbb R^+ \end{equation*}where \(y_1(t),\dots,y_n(t)\) are the coordinates of the solution \(Y(t)\) of the transformed system \(\dot Y(t) = \hat A\,Y(t)\) such that \(Y(0) = H^{-1}X(0)\). Each \(y_i(t)\) is of the form \(e^{\lambda_it}P_i(t)\), where \(\lambda_i\) is the eigenvalue associated with the column \(v_i\) and \(P_i\) a polynomial of degree strictly less than the multiplicity of \(\lambda_i\).
By definition, \(X(t) = H\,Y(t)\), which gives the first formula since the columns of \(H\) are the vectors \(v_i\). For the form of the \(y_i(t)\), consider the block of \(\hat A\) associated with the eigenvalue \(\lambda_k\) of multiplicity \(r_k\), and solve it from the bottom up as in the previous section, without right-hand side. The last equation gives \(y_{i+r_k}(t) = e^{\lambda_kt}P_{i+r_k}(t)\) with \(P_{i+r_k}\) constant, equal to \(y_{i+r_k}(0)\). Going up one row and substituting this solution, we obtain a scalar equation whose right-hand side is the product of \(e^{\lambda_kt}\) and a constant. By the integral formula, its solution is:
\begin{equation*} y_{i+r_k-1}(t) = e^{\lambda_kt}\left(y_{i+r_k-1}(0) + \hat a_{k,r_k-1,r_k}\int_0^t P_{i+r_k}(s)\,\mathrm ds\right) = e^{\lambda_kt}P_{i+r_k-1}(t) \end{equation*}where \(P_{i+r_k-1}\) is a polynomial of degree at most one, the integral of a polynomial being a polynomial of degree one unit higher. By induction going up the block, the right-hand side of row \(j\) is the product of \(e^{\lambda_kt}\) and a linear combination of the polynomials of the lower rows, and the solution is of the form \(e^{\lambda_kt}P_j(t)\) where the degree of \(P_j\) exceeds by at most one unit the largest degree of the lower rows. The first row of the block, reached after \(r_k-1\) steps, is therefore associated with a polynomial of degree at most \(r_k-1\). The argument applies in the same way to each eigenvalue.
When the transition matrix is diagonalisable, all the blocks have dimension one and the polynomials reduce to constants: the general solution is simply \(X(t) = \sum_{i=1}^n\gamma_ie^{\lambda_it}v_i\), where the constants \(\gamma_i\) are the coordinates of \(H^{-1}X(0)\). If \(A\) is real and \(\lambda = a + \mathrm ib\) is a complex eigenvalue with eigenvector \(v\), then \(\bar\lambda\) and \(\bar v\) also appear in the sum, with conjugate coefficients since the solution is real. The contribution of the pair is therefore \(2\,\mathrm{Re}\bigl(\gamma e^{\lambda t}v\bigr)\), a real combination of \(e^{at}\cos(bt)\) and \(e^{at}\sin(bt)\), as in the example of the previous section.
Consider the homogeneous system \(\dot X(t) = A\,X(t)\) with:
\begin{equation*} A = \begin{pmatrix} 1/2 & 1 & 1\\ 1 & 1/2 & 1\\ 1 & 1 & 1/2\end{pmatrix} \end{equation*}Show that the solution, given an initial condition \(X_0 = (x_{1,0},x_{2,0},x_{3,0})^\top\), is:
\begin{equation*} \begin{split} X(t) = {}& \frac{x_{1,0}+x_{2,0}-2x_{3,0}}{3}\,e^{-\frac t2}\begin{pmatrix} 1\\ 0\\ -1\end{pmatrix} + \frac{x_{1,0}+x_{3,0}-2x_{2,0}}{3}\,e^{-\frac t2}\begin{pmatrix} 1\\ -1\\ 0\end{pmatrix}\\ &+ \frac{x_{1,0}+x_{2,0}+x_{3,0}}{3}\,e^{\frac{5t}{2}}\begin{pmatrix} 1\\ 1\\ 1\end{pmatrix} \end{split} \end{equation*}for \(t\in\mathbb R^+\).
Solution
The matrix can be written \(A = J - \frac12I_3\), where \(J\) is the matrix all of whose entries are one. Since \(J\) has rank one, zero is a double eigenvalue of \(J\), and the trace of \(J\) being three, the third eigenvalue is \(3\). The eigenvalues of \(A\) are therefore \(-\frac12\) (double) and \(\frac52\). The matrix is real symmetric, hence diagonalisable: we must find two independent eigenvectors for \(-\frac12\), solutions of \(J\,v = 0\), that is \(v_1+v_2+v_3 = 0\), and one eigenvector for \(\frac52\), solution of \((J-3I_3)\,v = 0\):
\begin{equation*} v_1 = \begin{pmatrix} 1\\ 0\\ -1\end{pmatrix},\qquad v_2 = \begin{pmatrix} 1\\ -1\\ 0\end{pmatrix} \qquad\text{and}\qquad v_3 = \begin{pmatrix} 1\\ 1\\ 1\end{pmatrix} \end{equation*}By the previous remark, the general solution is:
\begin{equation*} X(t) = \alpha\,e^{-\frac t2}v_1 + \beta\,e^{-\frac t2}v_2 + \gamma\,e^{\frac{5t}{2}}v_3 \end{equation*}Evaluating this expression at \(t=0\) and identifying with \(X_0\), we get the system \(\alpha+\beta+\gamma = x_{1,0}\), \(-\beta+\gamma = x_{2,0}\), \(-\alpha+\gamma = x_{3,0}\). Summing the three equations gives \(3\gamma = x_{1,0}+x_{2,0}+x_{3,0}\), then \(\alpha = \gamma - x_{3,0}\) and \(\beta = \gamma - x_{2,0}\), which is the announced solution. The component along \(v_3\) dominates in the long run, unless \(x_{1,0}+x_{2,0}+x_{3,0} = 0\), in which case the solution converges to zero.
Linear equations of order \(n\)
A differential equation involving linearly the level of a function \(x\) and its first \(n\) derivatives is a linear differential equation of order \(n\). Normalising the coefficient of the \(n\)-th derivative, we write it in the form:
\begin{equation*} x^{(n)}(t) = a_1x^{(n-1)}(t) + a_2x^{(n-2)}(t) + \dots + a_nx(t),\qquad t\in\mathbb R^+ \end{equation*}where \(x^{(p)}\) denotes the \(p\)-th derivative of \(x\) and the \(a_i\) are real parameters. We assume \(a_n\neq0\), otherwise it is an equation of order \(n-1\) in \(\dot x\). We also assume, for the time being, that the equation has no right-hand side.
The equation of order \(n\) is equivalent to a system of \(n\) first-order equations. Setting \(x_1 = x^{(n-1)}\), \(x_2 = x^{(n-2)}\), …, \(x_n = x\), each variable is the derivative of the next one and the first satisfies the equation:
\begin{equation*} \begin{cases} \dot x_1(t) = a_1x_1(t) + a_2x_2(t) + \dots + a_nx_n(t)\\ \dot x_2(t) = x_1(t)\\ \quad\vdots\\ \dot x_n(t) = x_{n-1}(t) \end{cases} \end{equation*}that is, with \(X(t) = \bigl(x_1(t),\dots,x_n(t)\bigr)^\top\), the system \(\dot X(t) = A\,X(t)\) whose transition matrix is the companion matrix:
\begin{equation*} A = \begin{pmatrix} a_1 & a_2 & \dots & a_{n-1} & a_n\\ 1 & 0 & \dots & 0 & 0\\ 0 & 1 & \ddots & \vdots & \vdots\\ \vdots & \ddots & \ddots & 0 & 0\\ 0 & \dots & 0 & 1 & 0 \end{pmatrix} \end{equation*}The set of real (complex) solutions of the homogeneous linear differential equation of order \(n\) is a vector space over \(\mathbb R\) (\(\mathbb C\)) of dimension \(n\). The \(n\) initial conditions \(x(0)\), \(\dot x(0)\), …, \(x^{(n-1)}(0)\) determine a unique solution.
Immediate from the property on the dimension of the space of solutions of the equivalent homogeneous system, whose initial condition \(X(0)\) gathers precisely the \(n\) successive derivatives of \(x\) at zero.
The characteristic polynomial of the companion matrix is:
\begin{equation*} \chi_A(\lambda) = (-1)^n\bigl(\lambda^n - a_1\lambda^{n-1} - a_2\lambda^{n-2} - \dots - a_{n-1}\lambda - a_n\bigr) \end{equation*}Write \(D_n(\lambda) = |A-\lambda I_n|\) and argue by induction on \(n\). For \(n=1\), \(D_1(\lambda) = a_1 - \lambda\) and the formula holds. For \(n\geq2\), expand the determinant along the last column, whose only non-zero entries are \(a_n\) in the first row and \(-\lambda\) in the last row:
\begin{equation*} D_n(\lambda) = (-1)^{1+n}a_n\,M_{1,n} + (-1)^{2n}(-\lambda)\,M_{n,n} \end{equation*}where \(M_{i,j}\) is the minor obtained by deleting row \(i\) and column \(j\). The minor \(M_{n,n}\) is the determinant of the companion matrix of order \(n-1\) associated with \(a_1,\dots,a_{n-1}\), minus \(\lambda I_{n-1}\), so it is \(D_{n-1}(\lambda)\). The minor \(M_{1,n}\) is the determinant of the matrix formed by rows \(2\) to \(n\) and columns \(1\) to \(n-1\) of \(A-\lambda I_n\), which is upper triangular with \(1\)'s on the diagonal (and \(-\lambda\)'s just above): it equals \(1\). Thus \(D_n(\lambda) = (-1)^{n+1}a_n - \lambda D_{n-1}(\lambda)\), and substituting the induction hypothesis:
\begin{equation*} D_n(\lambda) = (-1)^{n+1}a_n + (-1)^n\bigl(\lambda^n - a_1\lambda^{n-1} - \dots - a_{n-1}\lambda\bigr) = (-1)^n\bigl(\lambda^n - a_1\lambda^{n-1} - \dots - a_{n-1}\lambda - a_n\bigr) \end{equation*}The eigenvalues of the companion matrix are therefore the roots of the polynomial \(\lambda^n - a_1\lambda^{n-1} - \dots - a_n\), which one obtains directly by looking for solutions of the form \(x(t) = e^{\lambda t}\) in the equation of order \(n\). We assume that this polynomial has \(m\leq n\) distinct roots \(\lambda_k\), with multiplicities \(r_k\) such that \(\sum_{k=1}^mr_k = n\).
The general solution of the homogeneous linear differential equation of order \(n\) is:
\begin{equation*} x(t) = \sum_{k=1}^m\sum_{p=0}^{r_k-1}\alpha_{k,p}\,t^p\,e^{\lambda_kt},\qquad t\in\mathbb R^+ \end{equation*}where the \(n\) constants \(\alpha_{k,p}\) are determined by the initial conditions.
Immediate by applying the theorem of the previous section to the equivalent system: \(x = x_n\) is a coordinate of \(X(t)\), hence a linear combination of the \(y_i(t) = e^{\lambda_it}P_i(t)\), and the degree of \(P_i\) is at most \(r_k-1\) when \(\lambda_i = \lambda_k\). The family of the \(n\) functions \(t^pe^{\lambda_kt}\) is linearly independent, and the space of solutions having dimension \(n\), it is a basis of it.
For \(n=2\), the equation \(\ddot x(t) = a_1\dot x(t) + a_2x(t)\) has characteristic polynomial \(\lambda^2 - a_1\lambda - a_2\). If the discriminant \(a_1^2+4a_2\) is strictly positive, the roots \(\lambda_1\) and \(\lambda_2\) are real and distinct and \(x(t) = \alpha_1e^{\lambda_1t} + \alpha_2e^{\lambda_2t}\). If it is zero, the double root \(\lambda = a_1/2\) gives \(x(t) = (\alpha_0 + \alpha_1t)\,e^{\lambda t}\). If it is strictly negative, the roots are conjugate, \(\lambda = a\pm\mathrm ib\) with \(a = a_1/2\) and \(b = \frac12\sqrt{-a_1^2-4a_2}\), and the real solution reads \(x(t) = e^{at}\bigl(\alpha_1\cos(bt) + \alpha_2\sin(bt)\bigr)\).
Let us now consider the equation with a right-hand side:
\begin{equation*} x^{(n)}(t) = a_1x^{(n-1)}(t) + \dots + a_nx(t) + b(t),\qquad t\in\mathbb R^+ \end{equation*}where \(b(t)\) is a continuous function from \(\mathbb R^+\) to \(\mathbb R\). Introducing a right-hand side raises no new problem: adding to a particular solution \(\psi(t)\) of the complete equation the general solution of the homogeneous equation, we obtain the general solution of the complete equation, since \(x(t)\) is a solution of the complete equation if and only if \(x(t)-\psi(t)\) is a solution of the homogeneous equation. It only remains to determine the \(n\) constants of integration from the initial conditions. For the usual right-hand sides, the form of the particular solution is known.
If the right-hand side is of the form \(b(t) = e^{\mu t}Q(t)\), where \(\mu\) is a real or complex number and \(Q\) a polynomial of degree \(q\), then the complete equation has a particular solution of the form \(e^{\mu t}R(t)\), with \(R\) a polynomial of degree \(q\), if \(\mu\) is not a root of the characteristic polynomial, and of the form \(t^re^{\mu t}R(t)\) if \(\mu\) is a root of order \(r\) of the characteristic polynomial.
We take this property for granted; the example of the previous section is an illustration with \(\mu = 0\): the right-hand side was polynomial of degree one, zero was not an eigenvalue, and the particular solution was polynomial of degree one. The coefficients of \(R\) are obtained by identification, substituting the postulated form into the equation.
Stability of a dynamic system
Definition and criterion
We consider the homogeneous system:
\begin{equation*} \dot X(t) = A\,X(t),\qquad t\in\mathbb R^+ \end{equation*}The real matrix \(A\) of dimension \(n\times n\) is said to be d-stable (differentially stable) if every solution \(X(t)\) of the homogeneous system tends to \(0\) as \(t\) tends to \(+\infty\). The zero steady state of the system is then said to be globally stable.
More generally, consider dynamics of the form \(\dot X(t) = A\,X(t) + b\), where the right-hand side \(b\) is a constant vector. A steady state is a constant solution \(X^{\star}\), that is, a point of \(\mathbb R^n\) satisfying \(A\,X^{\star} + b = 0\), i.e. \(X^{\star} = -A^{-1}b\) if \(A\) has full rank. The steady state is globally stable if, for any initial condition \(X_0\in\mathbb R^n\), the solution such that \(X(0) = X_0\) tends to \(X^{\star}\) as \(t\) tends to \(+\infty\). Since the deviation \(Z(t) = X(t) - X^{\star}\) satisfies \(\dot Z(t) = A\,Z(t)\), the steady state \(X^{\star}\) is globally stable if and only if \(A\) is d-stable.
The matrix \(A\) is d-stable if and only if all its eigenvalues have strictly negative real part.
We establish the necessary condition and the sufficient condition in turn.
Necessary condition. Suppose that \(A\) is d-stable. Let \(\lambda = a + \mathrm ib\) be an eigenvalue of \(A\) and \(v\in\mathbb C^n\) a non-zero associated eigenvector. Then \(Z(t) = e^{\lambda t}v\) is a complex solution of the system, since \(\dot Z(t) = \lambda e^{\lambda t}v = e^{\lambda t}A\,v = A\,Z(t)\). Its real and imaginary parts are real solutions, which tend to zero since \(A\) is d-stable, so \(Z(t)\) tends to zero. Let \(v_i\) be a non-zero coordinate of \(v\). The coordinate \(z_i(t) = v_ie^{\lambda t}\) tends to zero, which implies that \(|z_i(t)| = |v_i|\,e^{at}\) tends to zero as \(t\) tends to infinity, hence that \(a\) is strictly negative.
Sufficient condition. Suppose that all the eigenvalues \(\lambda_k = a_k + \mathrm ib_k\) of \(A\), with multiplicities \(r_k\), have strictly negative real parts \(a_k\) (\(k=1,\dots,m\)). Let \(X(t)\) be a solution of the system. By the theorem of the section on the characterisation of the solutions, each coordinate \(x_i(t)\) of \(X(t)\) is of the form \(\sum_{k=1}^me^{\lambda_kt}P_{i,k}(t)\), where the \(P_{i,k}\) are polynomials with complex coefficients, and:
\begin{equation*} |x_i(t)| = \left|\sum_{k=1}^me^{\lambda_kt}P_{i,k}(t)\right| \leq \sum_{k=1}^me^{a_kt}\,|P_{i,k}(t)| \end{equation*}Since \(a_k < 0\), each term \(e^{a_kt}|P_{i,k}(t)|\) tends to zero as \(t\) tends to infinity, the exponential dominating any polynomial. Thus \(X(t)\) tends to zero and the matrix \(A\) is d-stable.
Special cases
The previous criterion is easy to apply in low dimension.
- If \(n=1\), \(A = (a)\) is d-stable if and only if \(a < 0\).
- If \(n=2\), \(A = \begin{pmatrix} a & b\\ c & d\end{pmatrix}\) is d-stable if and only if \(\mathrm{tr}(A) = a + d < 0\) and \(|A| = ad - bc > 0\). Indeed, the trace is the sum of the eigenvalues and the determinant their product. If the eigenvalues are real, they are both strictly negative if and only if their sum is strictly negative and their product strictly positive. If they are complex conjugates, \(a\pm\mathrm ib\) with \(b\neq0\), their product \(a^2+b^2\) is always strictly positive and their sum \(2a\) is strictly negative if and only if their real part is. The figure below illustrates the three possible configurations when the determinant is non-zero.
If \(A\) is symmetric, its eigenvalues are real and it is d-stable if and only if it is negative definite, which by Sylvester's criterion reads:
\begin{equation*} (-1)^p\,|A_{1:p}| > 0\qquad\text{for all } 1\leq p\leq n \end{equation*}where \(A_{1:p}\) is the submatrix formed by the first \(p\) rows and columns of \(A\) (the leading principal minors alternate in sign, starting with \(a_{1,1} < 0\)).
Beyond dimension two, one can use the characteristic polynomial without computing its roots. The following proposition deals with the case \(n=3\).
For the polynomial \(\lambda^3 + a\lambda^2 + b\lambda + c\) to have all its roots with strictly negative real part, it is necessary and sufficient that its coefficients satisfy \(a > 0\), \(b > 0\), \(c > 0\) and \(ab > c\).
Let \(\lambda_1\), \(\lambda_2\) and \(\lambda_3\) be the three roots in \(\mathbb C\). Expanding \((\lambda-\lambda_1)(\lambda-\lambda_2)(\lambda-\lambda_3)\) and identifying, we have \(a = -(\lambda_1+\lambda_2+\lambda_3)\), \(b = \lambda_1\lambda_2+\lambda_1\lambda_3+\lambda_2\lambda_3\) and \(c = -\lambda_1\lambda_2\lambda_3\), whence the identity:
\begin{equation*} ab - c = -(\lambda_1+\lambda_2)(\lambda_2+\lambda_3)(\lambda_3+\lambda_1) \end{equation*}Necessary condition. If the three roots are real and strictly negative, we immediately have \(a > 0\), \(b > 0\), \(c > 0\), and each factor \(\lambda_i+\lambda_j\) being strictly negative, \(ab - c > 0\). If \(\lambda_1 = -r\) is real, with \(r > 0\), and \(\lambda_{2,3} = \alpha\pm\mathrm i\beta\) with \(\alpha < 0\), then \(a = r - 2\alpha > 0\), \(b = -2r\alpha + \alpha^2+\beta^2 > 0\), \(c = r(\alpha^2+\beta^2) > 0\) and, since \((\lambda_1+\lambda_2)(\lambda_1+\lambda_3) = |\lambda_1+\lambda_2|^2 > 0\) and \(\lambda_2+\lambda_3 = 2\alpha < 0\), again \(ab - c > 0\).
Sufficient condition. Suppose \(a > 0\), \(b > 0\), \(c > 0\) and \(ab > c\). A real polynomial of degree three has at least one real root \(\lambda_1\), and it is strictly negative since, the coefficients being all strictly positive, the polynomial is strictly positive for \(\lambda\geq0\). We factorise \(\lambda^3+a\lambda^2+b\lambda+c = (\lambda-\lambda_1)(\lambda^2+s\lambda+q)\), and by identification \(c = -\lambda_1q\), that is \(q = -c/\lambda_1 > 0\), then \(b = q - \lambda_1s\) and \(a = s - \lambda_1\). It follows that:
\begin{equation*} ab - c = (s-\lambda_1)(q-\lambda_1s) + \lambda_1q = s\bigl(q - \lambda_1s + \lambda_1^2\bigr) = s\bigl(b + \lambda_1^2\bigr) \end{equation*}hence \(s = (ab-c)/(b+\lambda_1^2) > 0\). The other two roots are those of the factor \(\lambda^2+s\lambda+q\), whose sum \(-s\) is strictly negative and whose product \(q\) is strictly positive: by the case \(n=2\), they have strictly negative real part.
The following two properties give necessary conditions for stability, easy to check, which allow one to rule out a matrix quickly.
If the real matrix \(A\) of dimension \(n\times n\) is d-stable, then its trace is strictly negative and its determinant has the sign of \((-1)^n\).
The trace is the sum of the eigenvalues and the determinant their product. The real eigenvalues are strictly negative, and the complex eigenvalues come in conjugate pairs \(a\pm\mathrm ib\) with \(a < 0\), whose sum \(2a\) is strictly negative and whose product \(a^2+b^2\) is strictly positive. The sum of all the eigenvalues is therefore strictly negative, and their product has the sign of \((-1)^{n_r}\), where \(n_r\) is the number of real eigenvalues, which has the parity of \(n\).
If the polynomial \(\lambda^n + a_1\lambda^{n-1} + \dots + a_{n-1}\lambda + a_n\) has all its roots with strictly negative real part, then all its coefficients are strictly positive.
The polynomial factorises as \((\lambda-\lambda_1)(\lambda-\lambda_2)\cdots(\lambda-\lambda_n)\). If \(\lambda_1 < 0\) is real, the coefficients of \(\lambda-\lambda_1\) are strictly positive. If \(\lambda_1 = a + \mathrm ib\) with \(a < 0\) and \(b\neq0\), we group the factor with that of the conjugate root: \((\lambda - a - \mathrm ib)(\lambda - a + \mathrm ib) = \lambda^2 - 2a\lambda + a^2+b^2\), whose coefficients are strictly positive. The polynomial is therefore the product of polynomials of degree one and two with strictly positive coefficients, and all its coefficients are strictly positive.
The characteristic polynomial of the companion matrix of the equation of order \(n\) is \((-1)^n(\lambda^n - a_1\lambda^{n-1} - \dots - a_n)\). The property applies to the monic polynomial \(\lambda^n - a_1\lambda^{n-1} - \dots - a_n\): a necessary condition for the stability of the equation of order \(n\) is that all the coefficients \(a_i\) be strictly negative. For \(n=2\), the condition \(a_1 < 0\) and \(a_2 < 0\) is also sufficient, since the trace of the companion matrix is \(a_1\) and its determinant \(-a_2\).
Diagonally dominant matrices
A square matrix \(A = (a_{i,j})\) of dimension \(n\times n\) with real entries is negative diagonally dominant if there exist \(n\) numbers \(d_i > 0\), \(i=1,\dots,n\), such that:
\begin{equation*} -d_i\,a_{i,i} > \sum_{j\neq i}d_j\,|a_{i,j}|\qquad\text{for all } 1\leq i\leq n \end{equation*}In particular, the diagonal entries of such a matrix are strictly negative, and each dominates, in absolute value, the other entries of its row weighted by the \(d_j\). This condition is easy to check, and it is enough to guarantee stability.
A negative diagonally dominant matrix is d-stable.
Let \(\lambda = \alpha + \mathrm i\beta\) be an eigenvalue of \(A\) and \(v\in\mathbb C^n\), non-zero, an associated eigenvector with coordinates \(v_j\). Choose an index \(k\) such that \(|v_k|/d_k = \max_j|v_j|/d_j\); we have \(v_k\neq0\) since \(v\) is non-zero. The \(k\)-th row of \(A\,v = \lambda v\) reads \((\lambda - a_{k,k})\,v_k = \sum_{j\neq k}a_{k,j}\,v_j\), whence, using \(|v_j|\leq d_j|v_k|/d_k\) for all \(j\):
\begin{equation*} \begin{split} |\lambda - a_{k,k}|\,|v_k| &\leq \sum_{j\neq k}|a_{k,j}|\,|v_j| \leq \frac{|v_k|}{d_k}\sum_{j\neq k}d_j\,|a_{k,j}|\\ &< \frac{|v_k|}{d_k}\bigl(-d_k\,a_{k,k}\bigr) = -a_{k,k}\,|v_k| \end{split} \end{equation*}where the last inequality is the definition of negative diagonal dominance. Dividing by \(|v_k| > 0\) gives \(|\lambda - a_{k,k}| < -a_{k,k}\), and since the real part of a complex number is bounded by its modulus, \(\alpha - a_{k,k} \leq |\lambda - a_{k,k}| < -a_{k,k}\), that is \(\alpha < 0\). All the eigenvalues have strictly negative real part and \(A\) is d-stable.
The Routh–Hurwitz criterion
The Routh–Hurwitz criterion generalises the proposition established for \(n=3\) by giving a necessary and sufficient condition on the coefficients of the characteristic polynomial, without computing the roots. Let \(p(\lambda) = \lambda^n + a_1\lambda^{n-1} + \dots + a_{n-1}\lambda + a_n\) be a monic polynomial with real coefficients4. Build the square matrix of dimension \(n\times n\):
\begin{equation*} R = \begin{pmatrix} a_1 & a_3 & a_5 & \dots & a_{2n-3} & a_{2n-1}\\ 1 & a_2 & a_4 & \dots & a_{2n-4} & a_{2n-2}\\ 0 & a_1 & a_3 & \dots & a_{2n-5} & a_{2n-3}\\ 0 & 1 & a_2 & \dots & a_{2n-6} & a_{2n-4}\\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots\\ 0 & 0 & 0 & \dots & a_{n-2} & a_n \end{pmatrix} \end{equation*}whose entry in row \(i\) and column \(j\) is \(a_{2i-j}\), with the convention \(a_0 = 1\) and \(a_k = 0\) if \(k < 0\) or \(k > n\). One can show that a necessary and sufficient condition for all the roots of \(p\) to have strictly negative real part is that the \(n\) leading principal minors of \(R\) (the determinants of the submatrices formed by the first \(p\) rows and columns, \(p=1,\dots,n\)) all be strictly positive. We take this result for granted, and content ourselves with checking that it recovers the conditions already established.
Write the Routh–Hurwitz matrix of the polynomial \(\lambda^2 + a_1\lambda + a_2\) and recover the stability condition of a two-dimensional system.
Solution
With \(n=2\), \(a_3 = 0\) and the matrix is:
\begin{equation*} R = \begin{pmatrix} a_1 & 0\\ 1 & a_2\end{pmatrix} \end{equation*}The leading principal minors are \(a_1\) and \(a_1a_2\). The roots have strictly negative real part if and only if:
\begin{equation*} \begin{cases} a_1 > 0\\ a_1a_2 > 0 \end{cases} \quad\Leftrightarrow\quad \begin{cases} a_1 > 0\\ a_2 > 0 \end{cases} \end{equation*}For a two-dimensional matrix \(A\), \(p(\lambda) = \lambda^2 - \mathrm{tr}(A)\,\lambda + |A|\), and we recover the condition \(\mathrm{tr}(A) < 0\) and \(|A| > 0\).
Write the Routh–Hurwitz matrix of the polynomial \(\lambda^3 + a_1\lambda^2 + a_2\lambda + a_3\) and recover the proposition established above for polynomials of degree three.
Solution
With \(n=3\), \(a_4 = a_5 = 0\) and the matrix is:
\begin{equation*} R = \begin{pmatrix} a_1 & a_3 & 0\\ 1 & a_2 & 0\\ 0 & a_1 & a_3\end{pmatrix} \end{equation*}The leading principal minors are \(a_1\), \(a_1a_2 - a_3\) and, expanding along the last column, \(a_3(a_1a_2 - a_3)\). The Routh–Hurwitz conditions therefore read \(a_1 > 0\), \(a_1a_2 > a_3\) and \(a_3(a_1a_2-a_3) > 0\), that is \(a_3 > 0\) given the second one. Finally \(a_1 > 0\), \(a_3 > 0\) and \(a_1a_2 > a_3\) imply \(a_2 > 0\). We recover exactly the conditions \(a_1 > 0\), \(a_2 > 0\), \(a_3 > 0\) and \(a_1a_2 > a_3\) of the proposition.
Footnotes:
It is a weak form of the Jordan reduction theorem, which further states that a suitable choice of the vectors leaves only \(0\) or \(1\) above the diagonal of each block. This refinement is not needed to solve the system.
The vectors obtained at step \(p\) lie in the kernel of \(N^p\), and the increasing sequence of kernels \(\ker N\subset\ker N^2\subset\dots\) reaches dimension \(r_k\).
See also chapter three of Philippe Michel, Cours de mathématiques pour économistes, Economica, 1989.
For a matrix \(A\), take \(p(\lambda) = (-1)^n\chi_A(\lambda)\), which has the same roots.