#+OPTIONS:   H:3 num:nil toc:nil \n:nil @:t ::t |:t ^:nil -:t f:t *:t TeX:t LaTeX:t skip:t d:t tags:not-in-toc creator:t timestamp:nil author:nil title:nil html5-fancy:t
#+HTML_DOCTYPE: html5
#+HTML_HEAD: <link rel="stylesheet" type="text/css" href="../../../css/stylesheet-blog.css">
#+HTML_HEAD: <link rel="stylesheet" type="text/css" href="../../../css/theorems-en.css">
#+HTML_HEAD: <link rel="stylesheet" type="text/css" href="../../../fontawesome/css/all.css" />
#+HTML_HEAD: <link rel="stylesheet" href="../../../highlight/styles/base16/solarized-light.min.css">
#+HTML_HEAD: <script src="../../../highlight/highlight.min.js"></script>
#+HTML_HEAD: <script>hljs.highlightAll();</script>
#+LANGUAGE: en-GB
#+STARTUP: latexpreview
#+TITLE: Systems of linear differential equations
#+DATE: September 2026
#+AUTHOR: Stéphane Adjemian
#+EMAIL: stephane.adjemian@univ-lemans.fr
#+PROPERTY: header-args:python :python /tmp/blog-sedl/bin/python

#+BEGIN_QUOTE
This note follows on from the one devoted to
[[https://stephane-adjemian.fr/en/posts/first-order-differential-equations/][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.
#+END_QUOTE

\\
\\
\\

#+BEGIN_SRC bash :results silent :exports none :async t
  python3 -m venv /tmp/blog-sedl
  source /tmp/blog-sedl/bin/activate
  pip install numpy scipy matplotlib
#+END_SRC

#+begin_src python :session sedl-en :exports none :results none
  import numpy as np
  import matplotlib.pyplot as plt
  from scipy.integrate import solve_ivp
#+end_src

* First-order linear systems
:PROPERTIES:
:CUSTOM_ID: systeme
:END:

#+BEGIN_definition
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^+$.
#+END_definition

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
:PROPERTIES:
:CUSTOM_ID: notation
:END:

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
[[https://stephane-adjemian.fr/en/posts/first-order-differential-equations/][previous note]],
whose results we shall use freely.

#+BEGIN_remarque
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.
#+END_remarque

#+BEGIN_property
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.
#+END_property

#+BEGIN_proof
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.
#+END_proof

** An equivalent transformation of the system
:PROPERTIES:
:CUSTOM_ID: transformation
:END:

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
:PROPERTIES:
:CUSTOM_ID: reduction
:END:

** Diagonalisation
:PROPERTIES:
:CUSTOM_ID: diagonalisation
:END:

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$.

#+BEGIN_definition
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$.
#+END_definition

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.

#+BEGIN_property
If the $n$ eigenvalues of $A$ are distinct, then $A$ is diagonalisable.
#+END_property

#+BEGIN_proof
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.
#+END_proof

#+BEGIN_remarque
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.
#+END_remarque

#+BEGIN_exemple
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*}
#+END_exemple

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.

#+begin_src python :session sedl-en :exports code :results none
  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)
#+end_src

** Reduced form of a non-diagonalisable matrix
:PROPERTIES:
:CUSTOM_ID: forme-reduite
:END:

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.

#+BEGIN_property
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$.
#+END_property

We take this result for granted[fn:: 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.] 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.

1. 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$.
2. 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/.
3. 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
vectors[fn:: 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$.]. 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$.

#+BEGIN_exemple
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$.
#+END_exemple

Numerically, the rank of $A-I_3$ confirms that there are only two
independent eigenvectors, and we recover the reduced form:

#+begin_src python :session sedl-en :exports code :results none
  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]]
#+end_src

* Solving the system
:PROPERTIES:
:CUSTOM_ID: resolution
:END:

** The transformed system
:PROPERTIES:
:CUSTOM_ID: systeme-transforme
:END:

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.

#+BEGIN_property
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*}
#+END_property

#+BEGIN_proof
This is the formula established for
[[https://stephane-adjemian.fr/en/posts/first-order-differential-equations/#non-homogene-variable][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$.
#+END_proof

Two cases arise.

1. 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.
2. 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
:PROPERTIES:
:CUSTOM_ID: existence-unicite
:END:

#+BEGIN_theorem
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$.
#+END_theorem

#+BEGIN_proof
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
[[#notation][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.
#+END_proof

** An example with complex eigenvalues
:PROPERTIES:
:CUSTOM_ID: exemple-complexe
:END:

#+BEGIN_exemple
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.
#+END_exemple

#+BEGIN_remarque
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 [[fig:foyer][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.
#+END_remarque

#+begin_src python :session sedl-en :exports code :results none
  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)))
#+end_src

#+begin_src python :session sedl-en :exports none :results none
  fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(11, 4))
  s = np.linspace(0, 3.6, 400)
  for r, phi in ((0.05, 0), (0.05, 2*np.pi/3), (0.05, 4*np.pi/3)):
      z0 = np.array([r*np.cos(phi), r*np.sin(phi)])
      z1 = np.exp(s)*(z0[0]*np.cos(s) - z0[1]*np.sin(s))
      z2 = np.exp(s)*(z0[0]*np.sin(s) + z0[1]*np.cos(s))
      ax1.plot(z1, z2, 'b', linewidth=1)
      ax1.annotate('', xy=(z1[-1], z2[-1]), xytext=(z1[-2], z2[-2]),
                   arrowprops=dict(arrowstyle='->', color='b'))
  ax1.plot(0, 0, 'ro', markersize=4)
  ax1.set_xlim(-2, 2); ax1.set_ylim(-2, 2)
  ax1.set_xlabel(r'$x_1$'); ax1.set_ylabel(r'$x_2$')
  ax1.set_title(r'Homogeneous system: unstable focus')
  ax2.plot(t, x1, 'b', linewidth=1, label=r'$x_1(t)$')
  ax2.plot(t, x2, 'r', linewidth=1, label=r'$x_2(t)$')
  ax2.plot(t[::20], numerical.y[0][::20], 'bo', markersize=3)
  ax2.plot(t[::20], numerical.y[1][::20], 'ro', markersize=3)
  ax2.set_xlabel(r'$t$')
  ax2.set_title(r'Complete system, $X(0)=(1, 0.5)$')
  ax2.legend()
  fig.tight_layout()
  fig.savefig("sedl-foyer.svg", transparent=True)
#+end_src

#+CAPTION: *Left: trajectories of the homogeneous system in the phase plane for three initial conditions close to the origin. Right: analytical solution of the complete system (lines) and numerical solution (dots).*
#+LABEL: fig:foyer
[[file:sedl-foyer.svg]]

#+NAME: ex-noeud
#+BEGIN_exercice
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.
#+END_exercice

#+ATTR_HTML: :class corrige
#+BEGIN_details
#+BEGIN_summary
Solution
#+END_summary
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
[[fig:noeud][below]] shows.
#+END_details

#+begin_src python :session sedl-en :exports none :results none
  fig, ax = plt.subplots(figsize=(6, 5))
  grid = np.linspace(-3.5, 1.5, 25)
  U, V = np.meshgrid(grid, grid)
  ax.streamplot(U, V, 3*U - 2*V + 1, U + 1, color='b', linewidth=0.7, density=1.1, arrowsize=0.8)
  for v, style in (((1, 1), 'r--'), ((2, 1), 'g--')):
      ax.plot([-1 - 3*v[0], -1 + 3*v[0]], [-1 - 3*v[1], -1 + 3*v[1]], style, linewidth=1.2)
  ax.plot(-1, -1, 'ko', markersize=4)
  ax.text(0.35, 0.8, r'$v_1$', color='r', bbox=dict(facecolor='white', edgecolor='none', pad=1))
  ax.text(0.9, 0.1, r'$v_2$', color='g', bbox=dict(facecolor='white', edgecolor='none', pad=1))
  ax.set_xlim(-3.5, 1.5); ax.set_ylim(-3.5, 1.5)
  ax.set_xlabel(r'$x$'); ax.set_ylabel(r'$y$')
  ax.set_title(r'Unstable node at $X^\star = (-1, -1)$')
  fig.tight_layout()
  fig.savefig("sedl-noeud.svg", transparent=True)
#+end_src

#+CAPTION: *Streamlines of the exercise in the phase plane. The dashed lines are the eigendirections: trajectories leave the steady state tangentially to $v_1$ and turn towards $v_2$.*
#+LABEL: fig:noeud
[[file:sedl-noeud.svg]]

** Characterisation of the solutions
:PROPERTIES:
:CUSTOM_ID: caracterisation
:END:

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.

#+BEGIN_property
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)$.
#+END_property

#+BEGIN_proof
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)$.
#+END_proof

#+BEGIN_property
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$.
#+END_property

#+BEGIN_proof
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 solutions[fn:: See also chapter three of
Philippe Michel, /Cours de mathématiques pour économistes/, Economica,
1989.].
#+END_proof

#+BEGIN_theorem
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$.
#+END_theorem

#+BEGIN_proof
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
[[#systeme-transforme][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.
#+END_proof

#+BEGIN_remarque
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 [[#exemple-complexe][previous section]].
#+END_remarque

#+NAME: ex-symetrique
#+BEGIN_exercice
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^+$.
#+END_exercice

#+ATTR_HTML: :class corrige
#+BEGIN_details
#+BEGIN_summary
Solution
#+END_summary
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.
#+END_details

* Linear equations of order $n$
:PROPERTIES:
:CUSTOM_ID: ordre-n
:END:

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)}$, \dots, $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*}

#+BEGIN_property
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)$, \dots,
$x^{(n-1)}(0)$ determine a unique solution.
#+END_property

#+BEGIN_proof
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.
#+END_proof

#+BEGIN_property
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*}
#+END_property

#+BEGIN_proof
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*}
#+END_proof

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$.

#+BEGIN_property
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.
#+END_property

#+BEGIN_proof
Immediate by applying the theorem of the [[#caracterisation][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.
#+END_proof

#+BEGIN_exemple
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)$.
#+END_exemple

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.

#+BEGIN_property
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.
#+END_property

We take this property for granted; the example of the
[[#exemple-complexe][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
:PROPERTIES:
:CUSTOM_ID: stabilite
:END:

** Definition and criterion
:PROPERTIES:
:CUSTOM_ID: critere
:END:

We consider the homogeneous system:

\begin{equation*}
\dot X(t) = A\,X(t),\qquad t\in\mathbb R^+
\end{equation*}

#+BEGIN_definition
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/.
#+END_definition

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.

#+BEGIN_theorem
The matrix $A$ is d-stable if and only if all its eigenvalues have
strictly negative real part.
#+END_theorem

#+BEGIN_proof
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
[[#caracterisation][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.
#+END_proof

** Special cases
:PROPERTIES:
:CUSTOM_ID: cas-particuliers
:END:

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
  [[fig:stabilite][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$).

#+begin_src python :session sedl-en :exports none :results none
  fig, axes = plt.subplots(1, 3, figsize=(12, 4))
  cases = ((np.array([[-2, 1], [1, -2]]), r'Stable node: tr $=-4$, det $=3$'),
           (np.array([[-0.5, -2], [2, -0.5]]), r'Stable focus: tr $=-1$, det $=4.25$'),
           (np.array([[1, 2], [2, 1]]), r'Saddle: tr $=2$, det $=-3$'))
  grid = np.linspace(-2, 2, 20)
  U, V = np.meshgrid(grid, grid)
  for ax, (A, title) in zip(axes, cases):
      dU = A[0, 0]*U + A[0, 1]*V
      dV = A[1, 0]*U + A[1, 1]*V
      ax.streamplot(U, V, dU, dV, color='b', linewidth=0.7, density=0.9, arrowsize=0.8)
      ax.plot(0, 0, 'ro', markersize=4)
      ax.set_xlim(-2, 2); ax.set_ylim(-2, 2)
      ax.set_xlabel(r'$x_1$'); ax.set_ylabel(r'$x_2$')
      ax.set_title(title)
  fig.tight_layout()
  fig.savefig("sedl-stabilite.svg", transparent=True)
#+end_src

#+CAPTION: *Phase portraits of three two-dimensional homogeneous systems. On the left and in the centre the trace is negative and the determinant positive, all trajectories converge to the origin, without oscillation (real eigenvalues) or spiralling (complex eigenvalues). On the right the determinant is negative, the eigenvalues are real with opposite signs and the origin is a saddle: only the trajectories starting on the stable eigendirection converge.*
#+LABEL: fig:stabilite
[[file:sedl-stabilite.svg]]

Beyond dimension two, one can use the characteristic polynomial without
computing its roots. The following proposition deals with the case $n=3$.

#+BEGIN_proposition
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$.
#+END_proposition

#+BEGIN_proof
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.
#+END_proof

The following two properties give necessary conditions for stability,
easy to check, which allow one to rule out a matrix quickly.

#+BEGIN_property
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$.
#+END_property

#+BEGIN_proof
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$.
#+END_proof

#+BEGIN_property
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.
#+END_property

#+BEGIN_proof
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.
#+END_proof

#+BEGIN_remarque
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$.
#+END_remarque

** Diagonally dominant matrices
:PROPERTIES:
:CUSTOM_ID: diagonale-dominante
:END:

#+BEGIN_definition
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*}
#+END_definition

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.

#+BEGIN_property
A negative diagonally dominant matrix is d-stable.
#+END_property

#+BEGIN_proof
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.
#+END_proof

** The Routh–Hurwitz criterion
:PROPERTIES:
:CUSTOM_ID: routh-hurwitz
:END:

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 coefficients[fn:: For a matrix $A$, take
$p(\lambda) = (-1)^n\chi_A(\lambda)$, which has the same roots.]. 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.

#+NAME: ex-rh2
#+BEGIN_exercice
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.
#+END_exercice

#+ATTR_HTML: :class corrige
#+BEGIN_details
#+BEGIN_summary
Solution
#+END_summary
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$.
#+END_details

#+NAME: ex-rh3
#+BEGIN_exercice
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.
#+END_exercice

#+ATTR_HTML: :class corrige
#+BEGIN_details
#+BEGIN_summary
Solution
#+END_summary
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.
#+END_details