#+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-fr.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: fr
#+STARTUP: latexpreview
#+TITLE: Systèmes d'équations différentielles linéaires
#+DATE: Septembre 2026
#+AUTHOR: Stéphane Adjemian
#+EMAIL: stephane.adjemian@univ-lemans.fr
#+PROPERTY: header-args:python :python /tmp/blog-sedl/bin/python

#+BEGIN_QUOTE
Cette note prolonge celle consacrée aux
[[https://stephane-adjemian.fr/posts/equations-differentielles-du-premier-ordre/][équations différentielles du premier ordre]]
en considérant plusieurs variables dont les variations dépendent
linéairement des niveaux de toutes les autres. La difficulté nouvelle tient à
cet enchevêtrement, et la stratégie consiste à le défaire par un changement
de base\nbsp{}: dans une base de vecteurs propres de la matrice de transition,
le système se ramène à des équations scalaires que l'on sait résoudre, au
prix d'un détour par les nombres complexes lorsque les valeurs propres ne
sont pas réelles. On traite ensuite le cas des matrices non diagonalisables,
les équations linéaires d'ordre $n$, qui ne sont que des systèmes déguisés,
et l'on termine par la stabilité, avec le théorème caractérisant les
matrices dont toutes les solutions convergent vers zéro et le critère de
Routh–Hurwitz. Les vérifications numériques et les figures sont faites en
Python, et les corrigés des exercices sont donnés dans des blocs dépliables.
#+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 :exports none :results none
  import numpy as np
  import matplotlib.pyplot as plt
  from scipy.integrate import solve_ivp
#+end_src

* Systèmes linéaires du premier ordre
:PROPERTIES:
:CUSTOM_ID: systeme
:END:

#+BEGIN_definition
Soient $a_{i,j}$ (pour $i=1,\dots,n$ et $j=1,\dots,n$) des paramètres réels
et $b_i(t)$ (pour $i=1,\dots,n$) des fonctions continues de $\mathbb R^+$
dans $\mathbb R$. Les fonctions $x_i(t)$, $i=1,\dots,n$, définies sur
$\mathbb R^+$ sont solutions du /système d'équations différentielles
linéaires du premier ordre/\nbsp{}:
\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*}
si elles sont dérivables sur $\mathbb R^+$ et si elles vérifient, avec leurs
dérivées $\dot x_i(t)$, ces $n$ équations pour tout $t$ dans $\mathbb R^+$.
#+END_definition

Le système est linéaire parce que les niveaux $x_j(t)$ n'y interviennent que
par des combinaisons linéaires, et du premier ordre parce que seules les
dérivées premières apparaissent. Les fonctions $b_i(t)$ forment le /second
membre/\nbsp{}; lorsqu'elles sont toutes nulles on dit que le système est
/homogène/, ou sans second membre.

** Notation matricielle
:PROPERTIES:
:CUSTOM_ID: notation
:END:

On note $X(t)\equiv\bigl(x_1(t),\dots,x_n(t)\bigr)^\top$ le vecteur des
variables, à valeurs dans $\mathbb R^n$, $\dot X(t)$ le vecteur de leurs
dérivées et $b(t)\equiv\bigl(b_1(t),\dots,b_n(t)\bigr)^\top$ le vecteur du
second membre. Le système de la définition s'écrit alors de façon
équivalente\nbsp{}:

\begin{equation*}
\dot X(t) = A\,X(t) + b(t)
\end{equation*}
où la matrice réelle carrée $A$, de dimension $n\times n$, est définie
par\nbsp{}:

\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*}
On appellera $A$ la /matrice de transition/ du système. Lorsque $n=1$ on
retrouve l'équation $\dot x(t) = a\,x(t) + b(t)$ étudiée dans la
[[https://stephane-adjemian.fr/posts/equations-differentielles-du-premier-ordre/][note précédente]],
dont nous utiliserons librement les résultats.

#+BEGIN_remarque
On s'intéresse aux solutions réelles du système. Néanmoins, la recherche de
ces solutions nous conduira à transformer le système, et nous serons parfois
amenés à calculer les solutions complexes du système transformé, c'est-à-dire
des fonctions de $\mathbb R^+$ dans $\mathbb C^n$. La propriété suivante
montre que ce détour ne pose pas de difficulté.
#+END_remarque

#+BEGIN_property
Si $A$ et $b(t)$ sont une matrice et un vecteur réels, alors la fonction
complexe $X(t) = X_1(t) + \mathrm i\,X_2(t)$, où $X_1$ et $X_2$ sont à valeurs
dans $\mathbb R^n$, est solution de $\dot X(t) = A\,X(t) + b(t)$ si et
seulement si\nbsp{}:
\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*}
La partie réelle d'une solution complexe est une solution réelle du système
complet, et sa partie imaginaire une solution réelle du système homogène.
#+END_property

#+BEGIN_proof
En substituant $X(t) = X_1(t) + \mathrm i\,X_2(t)$ dans le système, il vient
\(\dot X_1(t) + \mathrm i\,\dot X_2(t) = A\,X_1(t) + \mathrm i\,A\,X_2(t) + b(t)\).
Puisque $A$ et $b(t)$ sont réels, $A X_1(t) + b(t)$ et $A X_2(t)$ sont réels,
et deux vecteurs complexes sont égaux si et seulement si leurs parties
réelles et imaginaires le sont.
#+END_proof

** Transformation équivalente du système
:PROPERTIES:
:CUSTOM_ID: transformation
:END:

Résoudre directement le système n'est généralement pas possible, car la
variation de chaque variable, $\dot x_i(t)$, dépend des niveaux de toutes les
variables. L'idée est de considérer de nouvelles variables, combinaisons
linéaires indépendantes des variables originales, choisies de sorte que la
variation de chaque nouvelle variable ne dépende que de son propre niveau (et
du second membre). Si l'on y parvient, le système se ramène à $n$ équations
différentielles scalaires que nous savons résoudre.

Plus formellement, soit $H$ une matrice carrée inversible de dimension
$n\times n$, à éléments dans $\mathbb C$, que nous interprétons comme une
matrice de changement de base. Posons\nbsp{}:

\begin{equation*}
Y(t) = H^{-1}X(t)
\end{equation*}
Le vecteur $Y(t)$ est composé de $n$ fonctions de $\mathbb R^+$ dans
$\mathbb C$, et par linéarité de la dérivation nous avons
$\dot Y(t) = H^{-1}\dot X(t)$. En substituant $X(t) = H\,Y(t)$ dans le
système, puis en prémultipliant par $H^{-1}$, il vient\nbsp{}:

\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*}

soit, en posant $\hat A \equiv H^{-1}AH$ et \(c(t)\equiv H^{-1}b(t)\)\nbsp{}:

\begin{equation*}
\dot Y(t) = \hat A\,Y(t) + c(t)
\end{equation*}
Ce /système transformé/ est équivalent au système initial\nbsp{}: $X(t)$ est
solution du premier si et seulement si $Y(t) = H^{-1}X(t)$ est solution du
second, et les conditions initiales se correspondent par $Y(0) = H^{-1}X(0)$.
Tout le problème est de choisir $H$ pour que $\hat A$ soit aussi simple que
possible. Le cas idéal est celui où $\hat A$ est diagonale\nbsp{}; à défaut, nous
verrons que l'on peut toujours obtenir une matrice triangulaire par blocs,
ce qui suffit à résoudre le système de proche en proche.

* Réduction de la matrice de transition
:PROPERTIES:
:CUSTOM_ID: reduction
:END:

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

Soit $A$ une matrice carrée réelle de dimension $n\times n$ et $I_n$ la
matrice identité de même dimension. Le /polynôme caractéristique/ de $A$ est
le polynôme de degré $n$ en \(\lambda\)\nbsp{}:

\begin{equation*}
\chi_A(\lambda) = \left|A - \lambda I_n\right|
\end{equation*}
où $|\cdot|$ désigne le déterminant. Les /valeurs propres/ de $A$ sont les
racines, réelles ou complexes, de $\chi_A$, c'est-à-dire les nombres
$\lambda$ tels que $|A-\lambda I_n| = 0$. Le polynôme caractéristique
possède au plus $n$ racines distinctes. On note $\lambda_k$, $k=1,\dots,m$
(avec $m\leq n$), ces racines distinctes et $r_k$ leurs ordres de
multiplicité, avec $\sum_{k=1}^m r_k = n$. Pour chaque valeur propre
$\lambda_k$, la matrice $A-\lambda_kI_n$ est singulière, et il existe donc
des vecteurs $v$ non nuls tels que\nbsp{}:

\begin{equation*}
(A-\lambda_kI_n)\,v = 0
\qquad\Leftrightarrow\qquad
A\,v = \lambda_k v
\end{equation*}
Un tel vecteur est un /vecteur propre/ associé à $\lambda_k$. L'ensemble des
solutions de cette équation linéaire est un sous-espace vectoriel de
$\mathbb C^n$, le /sous-espace propre/ associé à $\lambda_k$, dont la
dimension $q_k$ vérifie $1\leq q_k\leq r_k$.

#+BEGIN_definition
La matrice $A$ est /diagonalisable/ si, pour chaque valeur propre
$\lambda_k$, on peut trouver $r_k$ vecteurs propres linéairement
indépendants, c'est-à-dire si $q_k = r_k$ pour tout $k=1,\dots,m$.
#+END_definition

Dans ce cas, les $n$ vecteurs propres ainsi obtenus sont linéairement
indépendants et forment les colonnes d'une matrice de passage
$H\equiv(v_1,\dots,v_n)$ inversible. Puisque $A\,v_i = \lambda_i v_i$ pour
chaque colonne, on a $A\,H = H\,\hat A$ avec\nbsp{}:

\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*}
où chaque valeur propre apparaît sur la diagonale autant de fois que son
ordre de multiplicité.

#+BEGIN_property
Si les $n$ valeurs propres de $A$ sont distinctes, alors $A$ est
diagonalisable.
#+END_property

#+BEGIN_proof
Chaque valeur propre est de multiplicité $r_k=1$ et possède au moins un
vecteur propre, donc $q_k = r_k$. Il reste à vérifier que des vecteurs
propres associés à des valeurs propres distinctes sont linéairement
indépendants. Supposons que $v_1,\dots,v_p$, associés à des valeurs propres
distinctes $\lambda_1,\dots,\lambda_p$, vérifient
$\alpha_1v_1+\dots+\alpha_pv_p = 0$. En appliquant $A-\lambda_pI_n$ aux deux
membres, le dernier terme disparaît et l'on obtient
$\alpha_1(\lambda_1-\lambda_p)v_1+\dots+\alpha_{p-1}(\lambda_{p-1}-\lambda_p)v_{p-1}=0$.
En répétant l'opération avec $\lambda_{p-1},\dots,\lambda_2$ il ne reste que
$\alpha_1\prod_{j=2}^p(\lambda_1-\lambda_j)v_1 = 0$, d'où $\alpha_1=0$ puisque
les valeurs propres sont distinctes et $v_1\neq0$. Le même argument
appliqué à chaque indice montre que tous les $\alpha_i$ sont nuls.
#+END_proof

#+BEGIN_remarque
Si $A$ est réelle, son polynôme caractéristique est à coefficients réels et
ses racines complexes vont par paires conjuguées. Si $v$ est un vecteur
propre associé à $\lambda$, alors $\bar v$ est un vecteur propre associé à
$\bar\lambda$, comme on le voit en conjuguant $Av=\lambda v$. Une matrice
réelle symétrique a toutes ses valeurs propres réelles et est toujours
diagonalisable, avec des vecteurs propres orthogonaux.
#+END_remarque

#+BEGIN_exemple
Soit la matrice réelle symétrique\nbsp{}:
\begin{equation*}
A = \begin{pmatrix} 2 & 1\\ 1 & 2\end{pmatrix}
\end{equation*}
Le polynôme caractéristique associé est\nbsp{}:
\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*}
Les valeurs propres sont donc réelles et distinctes\nbsp{}: $\lambda_1=1$ et
$\lambda_2=3$. Le vecteur propre $v_1$ associé à $\lambda_1$ doit vérifier
$(A-\lambda_1I_2)\,v_1 = 0$, c'est-à-dire\nbsp{}:
\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*}
Un vecteur propre n'est défini qu'à un facteur multiplicatif près, on
posera $v_1 = (1,-1)^\top$. De même, le vecteur propre $v_2$ associé à
$\lambda_2$ vérifie \((A-\lambda_2I_2)\,v_2 = 0\)\nbsp{}:
\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*}
et l'on posera $v_2 = (1,1)^\top$. La matrice de passage et son inverse
sont\nbsp{}:
\begin{equation*}
H = \begin{pmatrix} 1 & 1\\ -1 & 1\end{pmatrix}
\qquad\text{et}\qquad
H^{-1} = \frac{1}{2}\begin{pmatrix} 1 & -1\\ 1 & 1\end{pmatrix}
\end{equation*}
et l'on vérifie que\nbsp{}:
\begin{equation*}
H^{-1}AH = \begin{pmatrix} 1 & 0\\ 0 & 3\end{pmatrix}
\end{equation*}
#+END_exemple

On peut vérifier ces calculs numériquement. La fonction ~eig~ de ~numpy~
renvoie les valeurs propres et une matrice dont les colonnes sont des
vecteurs propres normalisés (de norme un), ce qui ne change rien à la
diagonalisation puisque les vecteurs propres ne sont définis qu'à un facteur
près.

#+begin_src python :session sedl :exports code :results none
  A = np.array([[2, 1], [1, 2]])
  lam, V = np.linalg.eig(A)
  print(lam)                          # 3 et 1
  print(V)                            # colonnes proportionnelles à (1,1) et (1,-1)
  H = np.array([[1, 1], [-1, 1]])
  print(np.linalg.inv(H) @ A @ H)     # diag(1, 3)
#+end_src

** Forme réduite d'une matrice non diagonalisable
:PROPERTIES:
:CUSTOM_ID: forme-reduite
:END:

Que faire lorsque la matrice de transition n'est pas diagonalisable,
c'est-à-dire lorsque pour une valeur propre $\lambda_k$ de multiplicité
$r_k$ on ne trouve que $q_k < r_k$ vecteurs propres linéairement
indépendants\nbsp{}? On ne peut plus rendre $\hat A$ diagonale, mais on peut
la rendre triangulaire par blocs, ce qui suffira pour résoudre le système
de façon récursive.

#+BEGIN_property
Soit $A$ une matrice carrée de dimension $n\times n$, à éléments réels ou
complexes, dont les valeurs propres distinctes sont $\lambda_1,\dots,\lambda_m$
de multiplicités $r_1,\dots,r_m$. Il existe une matrice de passage $H$
inversible telle que $\hat A = H^{-1}AH$ soit diagonale par blocs, avec un
bloc par valeur propre distincte. Le bloc associé à $\lambda_k$ est une
matrice triangulaire supérieure de dimension $r_k\times r_k$ dont tous les
éléments diagonaux sont égaux à $\lambda_k$.
#+END_property

Nous admettons ce résultat[fn:: Il s'agit d'une forme faible du théorème de
réduction de Jordan, qui précise en outre qu'un choix convenable des
vecteurs permet de n'avoir que des $0$ ou des $1$ au-dessus de la diagonale
de chaque bloc. Cette précision n'est pas nécessaire pour résoudre le
système.] et décrivons la construction des colonnes de $H$ associées à une
valeur propre $\lambda_k$, en notant $N \equiv A-\lambda_kI_n$. On procède par
étapes.

1. On choisit une base $v_1,\dots,v_{q_k}$ du sous-espace propre, c'est-à-dire
   $q_k$ vecteurs non nuls linéairement indépendants tels que $N v_i = 0$.
2. Si $q_k < r_k$, on cherche un vecteur $w$, n'appartenant pas au sous-espace
   engendré par les vecteurs déjà obtenus, tel que $Nw$ appartienne à ce
   sous-espace\nbsp{}:
   \begin{equation*}
   (A-\lambda_kI_n)\,w = \alpha_1v_1+\dots+\alpha_{q_k}v_{q_k}
   \end{equation*}
   pour des coefficients $\alpha_i$ non tous nuls que l'on détermine en même
   temps que $w$ (le système linéaire en $w$ n'a de solution que pour
   certaines valeurs des $\alpha_i$). On dit que $w$ est un /vecteur propre
   généralisé/.
3. On ajoute $w$ à la liste des vecteurs et, tant que l'on n'a pas obtenu
   $r_k$ vecteurs linéairement indépendants, on recommence l'étape 2 en
   cherchant un nouveau vecteur dont l'image par $N$ appartient au
   sous-espace engendré par tous les vecteurs déjà obtenus.

On montre que cette construction aboutit toujours à $r_k$ vecteurs
linéairement indépendants[fn:: Les vecteurs obtenus à l'étape $p$ sont dans
le noyau de $N^p$, et la suite croissante des noyaux
$\ker N\subset\ker N^2\subset\dots$ atteint la dimension $r_k$.]. Si l'on
range ces vecteurs, dans l'ordre où ils ont été obtenus, dans les colonnes
de $H$, le bloc correspondant de $H^{-1}AH$ est triangulaire supérieur\nbsp{}:
en effet $A\,v_i = \lambda_kv_i$ pour les vecteurs propres, et pour chaque
vecteur généralisé $A\,w = \lambda_kw + \alpha_1v_1+\dots$ ne fait intervenir
que $w$ et des vecteurs placés avant lui dans $H$.

#+BEGIN_exemple
Soit la matrice carrée réelle\nbsp{}:
\begin{equation*}
A = \begin{pmatrix} 3 & 3 & 1\\ -2 & -2 & -1\\ 2 & 3 & 2\end{pmatrix}
\end{equation*}
Son polynôme caractéristique est $\chi_A(\lambda) = (1-\lambda)^3$, qui
admet une racine triple $\lambda=1$. Les vecteurs propres sont les solutions
du système \((A-I_3)\,v = 0\)\nbsp{}:
\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*}
qui est de rang un. Le sous-espace propre est donc de dimension deux, il
n'existe que deux vecteurs propres linéairement indépendants et la matrice
$A$ n'est pas diagonalisable. On peut choisir\nbsp{}:
\begin{equation*}
v_1 = \begin{pmatrix} 0\\ 1\\ -3\end{pmatrix}
\qquad\text{et}\qquad
v_2 = \begin{pmatrix} 1\\ 0\\ -2\end{pmatrix}
\end{equation*}
Le système $(A-I_3)\,w = \alpha_1v_1+\alpha_2v_2$ s'écrit\nbsp{}:
\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*}
Les deux premières équations ne sont compatibles que si $\alpha_2 = -\alpha_1$,
et la troisième est alors identique à la première. En posant $\alpha_1=1$ et
$\alpha_2 = -1$, il reste $2w_1+3w_2+w_3 = -1$ et l'on peut choisir
$w = (0,0,-1)^\top$, qui n'est pas combinaison linéaire de $v_1$ et $v_2$.
La matrice de passage et son inverse sont\nbsp{}:
\begin{equation*}
H = \begin{pmatrix} 0 & 1 & 0\\ 1 & 0 & 0\\ -3 & -2 & -1\end{pmatrix}
\qquad\text{et}\qquad
H^{-1} = \begin{pmatrix} 0 & 1 & 0\\ 1 & 0 & 0\\ -2 & -3 & -1\end{pmatrix}
\end{equation*}
et l'on obtient\nbsp{}:
\begin{equation*}
\hat A = H^{-1}AH = \begin{pmatrix} 1 & 0 & 1\\ 0 & 1 & -1\\ 0 & 0 & 1\end{pmatrix}
\end{equation*}
La dernière colonne de $\hat A$ n'est autre que le vecteur
$(\alpha_1,\alpha_2,1)^\top$, conformément à $A\,w = w + v_1 - v_2$.
#+END_exemple

Numériquement, le rang de $A-I_3$ confirme qu'il n'y a que deux vecteurs
propres indépendants, et l'on retrouve la forme réduite\nbsp{}:

#+begin_src python :session sedl :exports code :results none
  A = np.array([[3, 3, 1], [-2, -2, -1], [2, 3, 2]])
  print(np.linalg.eigvals(A))                    # 1, 1, 1 (aux erreurs d'arrondi près)
  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

* Résolution du système
:PROPERTIES:
:CUSTOM_ID: resolution
:END:

** Le système transformé
:PROPERTIES:
:CUSTOM_ID: systeme-transforme
:END:

Revenons au système transformé $\dot Y(t) = \hat A\,Y(t) + c(t)$, où
$\hat A = H^{-1}AH$ est une forme diagonale ou réduite de $A$. Les éléments
de $\hat A$ et les fonctions $c_i(t)$ sont a priori complexes. Nous aurons
besoin de la solution de l'équation scalaire du premier ordre avec un
coefficient complexe.

#+BEGIN_property
Soient $\lambda\in\mathbb C$ et $c(t)$ une fonction continue de $\mathbb R^+$
dans $\mathbb C$. Pour toute condition initiale $y(0) = y_0\in\mathbb C$,
l'équation différentielle $\dot y(t) = \lambda\,y(t) + c(t)$ admet une unique
solution sur \(\mathbb R^+\)\nbsp{}:
\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
C'est la formule établie pour les
[[https://stephane-adjemian.fr/posts/equations-differentielles-du-premier-ordre/#non-homogene-variable][coefficients réels]]
dans la note précédente, et la preuve est la même\nbsp{}: la dérivée de
$e^{-\lambda t}$ est $-\lambda e^{-\lambda t}$ pour tout $\lambda$ complexe, de sorte
que l'équation équivaut à $\frac{\mathrm d}{\mathrm dt}\bigl(e^{-\lambda t}y(t)\bigr) = e^{-\lambda t}c(t)$,
que l'on intègre entre $0$ et $t$.
#+END_proof

On distingue deux cas.

1. Si $\hat A$ est diagonale, le système transformé se réduit à $n$ équations
   indépendantes $\dot y_i(t) = \lambda_iy_i(t) + c_i(t)$, que la propriété
   précédente résout une à une.
2. Si $\hat A$ est diagonale par blocs, chaque bloc peut être traité
   séparément. Considérons le bloc associé à la valeur propre $\lambda_k$ de
   multiplicité $r_k$, et supposons qu'il occupe les lignes $i+1$ à $i+r_k$
   de $\hat A$. Les variables correspondantes vérifient un sous-système de
   $r_k$ équations, indépendant des $n-r_k$ autres variables\nbsp{}:
   \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*}
   Le bloc étant triangulaire supérieur, il suffit de résoudre par le bas.
   La dernière équation, $\dot y_{i+r_k}(t) = \lambda_ky_{i+r_k}(t) + c_{i+r_k}(t)$,
   est une équation scalaire que nous savons résoudre\nbsp{}: étant donnée la
   condition initiale $y_{i+r_k}(0)$, elle admet une unique solution. Une fois
   cette solution obtenue, on la substitue dans l'équation du dessus, qui
   devient à son tour une équation scalaire du premier ordre en
   $y_{i+r_k-1}(t)$ avec un second membre connu. En remontant ainsi le
   triangle, on obtient une solution unique pour les $r_k$ variables du bloc,
   étant données leurs conditions initiales.

Dans tous les cas, le système transformé admet, pour tout vecteur initial
$Y_0\in\mathbb C^n$, une solution $Y(t)$ et une seule de $\mathbb R^+$ dans
$\mathbb C^n$ telle que $Y(0) = Y_0$.

** Existence et unicité
:PROPERTIES:
:CUSTOM_ID: existence-unicite
:END:

#+BEGIN_theorem
Si les fonctions $b_i(t)$ du second membre sont continues sur $\mathbb R^+$,
alors pour tout vecteur $X_0\in\mathbb R^n$ le système
$\dot X(t) = A\,X(t) + b(t)$ admet une solution réelle et une seule sur
$\mathbb R^+$ telle que $X(0) = X_0$. C'est aussi l'unique solution complexe
vérifiant cette condition initiale, et elle est égale à $H\,Y(t)$, où $Y(t)$
est la solution du système transformé telle que $Y(0) = H^{-1}X_0$.
#+END_theorem

#+BEGIN_proof
On a vu que le système transformé admet une solution complexe et une seule
$Y(t)$ telle que $Y(0) = H^{-1}X_0$. Puisque les deux systèmes sont
équivalents, $X(t) = H\,Y(t)$ est l'unique solution complexe du système
initial telle que $X(0) = X_0$. Écrivons cette solution
$X(t) = X_1(t) + \mathrm i\,X_2(t)$. D'après la propriété établie dans la
[[#notation][première section]], la partie imaginaire vérifie
$\dot X_2(t) = A\,X_2(t)$ avec, puisque $X_0$ est réel, $X_2(0) = 0$. Or la
fonction nulle est solution de ce système homogène avec cette condition
initiale, et cette solution est unique\nbsp{}: $X_2(t)$ est donc nulle à tout
instant, et l'unique solution complexe est réelle.
#+END_proof

** Un exemple avec des valeurs propres complexes
:PROPERTIES:
:CUSTOM_ID: exemple-complexe
:END:

#+BEGIN_exemple
On veut résoudre le système d'équations différentielles\nbsp{}:
\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*}
L'écriture matricielle équivalente est $\dot X(t) = A\,X(t) + b(t)$ avec\nbsp{}:
\begin{equation*}
A = \begin{pmatrix} 1 & -1\\ 1 & 1\end{pmatrix}
\qquad\text{et}\qquad
b(t) = \begin{pmatrix} 2\\ -2t\end{pmatrix}
\end{equation*}
Le polynôme caractéristique est\nbsp{}:
\begin{equation*}
\chi_A(\lambda) = \begin{vmatrix} 1-\lambda & -1\\ 1 & 1-\lambda\end{vmatrix} = (1-\lambda)^2 + 1
\end{equation*}
dont les racines sont $\lambda_1 = 1-\mathrm i$ et $\lambda_2 = 1+\mathrm i$.
Les valeurs propres sont distinctes, la matrice est donc diagonalisable dans
$\mathbb C$. Le vecteur propre $v_1$ associé à $\lambda_1$ vérifie\nbsp{}:
\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*}
La seconde équation est égale à $-\mathrm i$ fois la première, il ne reste
donc qu'une condition, $v_{2,1} = \mathrm i\,v_{1,1}$, et l'on posera
$v_1 = (-\mathrm i, 1)^\top$. Le vecteur propre associé à $\lambda_2$ est le
conjugué, $v_2 = (\mathrm i, 1)^\top$, comme on le vérifie directement\nbsp{}:
\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*}
La matrice de changement de base et son inverse sont\nbsp{}:
\begin{equation*}
H = \begin{pmatrix} -\mathrm i & \mathrm i\\ 1 & 1\end{pmatrix}
\qquad\text{et}\qquad
H^{-1} = \frac{1}{2}\begin{pmatrix} \mathrm i & 1\\ -\mathrm i & 1\end{pmatrix}
\end{equation*}
En posant $Y(t) = H^{-1}X(t)$ et $c(t) = H^{-1}b(t)$ il vient
$\dot Y(t) = \hat A\,Y(t) + c(t)$ avec\nbsp{}:
\begin{equation*}
\hat A = \begin{pmatrix} 1-\mathrm i & 0\\ 0 & 1+\mathrm i\end{pmatrix}
\qquad\text{et}\qquad
c(t) = \begin{pmatrix} \mathrm i - t\\ -\mathrm i - t\end{pmatrix}
\end{equation*}
soit, de façon équivalente\nbsp{}:
\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*}
Ces deux équations différentielles linéaires complexes sont indépendantes.
Plutôt que d'appliquer la formule intégrale, on procède comme dans le cas
réel\nbsp{}: on cherche une solution particulière de chaque équation complète
sous la forme d'un polynôme $\alpha + \beta t$, dont les coefficients complexes
sont identifiés en substituant dans l'équation, puis on ajoute la solution
générale de l'équation sans second membre. Pour la première équation,
$\beta = (1-\mathrm i)(\alpha+\beta t) + \mathrm i - t$ pour tout $t$ impose
$(1-\mathrm i)\beta = 1$, soit $\beta = \frac{1+\mathrm i}{2}$, puis
$\beta = (1-\mathrm i)\alpha + \mathrm i$, soit $\alpha = \frac12$. Le calcul
est identique pour la seconde équation, et l'on obtient\nbsp{}:
\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*}
où les constantes d'intégration sont déterminées par
\(Y(0) = H^{-1}X(0)\)\nbsp{}:
\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*}
avec $X(0) = (x_{0,1},x_{0,2})^\top$. On remarque que $y_2(t) = \overline{y_1(t)}$,
ce qui n'est pas un hasard\nbsp{}: la solution $X(t)$ étant réelle, les deux
coordonnées de $Y(t) = H^{-1}X(t)$ sont conjuguées puisque les lignes de
$H^{-1}$ le sont. Revenons au problème initial, la solution est
$X(t) = H\,Y(t)$, soit\nbsp{}:
\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*}
En remplaçant $e^{(1\pm\mathrm i)t} = e^t(\cos t \pm \mathrm i\sin t)$ et
$\gamma_1 = \frac12\bigl(\mathrm i\,x_{0,1} + x_{0,2} - 1\bigr)$,
$\gamma_2 = \overline{\gamma_1}$, on obtient après simplification\nbsp{}:
\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*}
On vérifie sans peine que ces fonctions satisfont le système et la condition
initiale.
#+END_exemple

#+BEGIN_remarque
La solution particulière $(t, t+1)^\top$ est polynomiale, comme le second
membre, et la solution du système homogène est une combinaison de
$e^t\cos t$ et \(e^t\sin t\)\nbsp{}: la partie réelle des valeurs propres, $1$,
gouverne l'amplitude et la partie imaginaire, $\pm1$, la fréquence des
oscillations. Le module de la solution homogène croît sans limite, la
dynamique est instable et les trajectoires s'éloignent de l'origine en
spiralant. La figure [[fig:foyer][ci-dessous]] représente quelques
trajectoires du système homogène dans le plan $(x_1,x_2)$, ainsi que la
solution du système complet, où l'on a superposé une résolution numérique
obtenue avec ~solve_ivp~.
#+END_remarque

#+begin_src python :session sedl :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)
  numerique = 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(numerique.y[0] - x1)), np.max(np.abs(numerique.y[1] - x2)))
#+end_src

#+begin_src python :session sedl :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'Système homogène : foyer instable')
  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], numerique.y[0][::20], 'bo', markersize=3)
  ax2.plot(t[::20], numerique.y[1][::20], 'ro', markersize=3)
  ax2.set_xlabel(r'$t$')
  ax2.set_title(r'Système complet, $X(0)=(1, 0{,}5)$')
  ax2.legend()
  fig.tight_layout()
  fig.savefig("sedl-foyer.svg", transparent=True)
#+end_src

#+CAPTION: *À gauche, trajectoires du système homogène dans le plan de phase pour trois conditions initiales proches de l'origine. À droite, solution analytique du système complet (traits) et résolution numérique (points).*
#+LABEL: fig:foyer
[[file:sedl-foyer.svg]]

#+NAME: ex-noeud
#+BEGIN_exercice
Résolvez le système d'équations différentielles\nbsp{}:
\begin{equation*}
\begin{cases}
\dot x(t) = 3x(t) - 2y(t) + 1\\
\dot y(t) = x(t) + 1
\end{cases}
\end{equation*}
en le réécrivant sous forme matricielle, puis en le transformant dans une
base de vecteurs propres. Caractérisez l'état stationnaire et sa stabilité.
#+END_exercice

#+ATTR_HTML: :class corrige
#+BEGIN_details
#+BEGIN_summary
Corrigé
#+END_summary
En posant $X(t) = (x(t),y(t))^\top$, le système s'écrit
$\dot X(t) = A\,X(t) + b$ avec\nbsp{}:
\begin{equation*}
A = \begin{pmatrix} 3 & -2\\ 1 & 0\end{pmatrix}
\qquad\text{et}\qquad
b = \begin{pmatrix} 1\\ 1\end{pmatrix}
\end{equation*}
Le polynôme caractéristique est
$\chi_A(\lambda) = (3-\lambda)(-\lambda) + 2 = \lambda^2 - 3\lambda + 2 = (\lambda-1)(\lambda-2)$,
dont les racines $\lambda_1 = 1$ et $\lambda_2 = 2$ sont réelles et
distinctes. Les vecteurs propres vérifient
$(A-I_2)\,v_1 = 0 \Leftrightarrow 2v_{1,1} - 2v_{2,1} = 0$ et
\((A-2I_2)\,v_2 = 0 \Leftrightarrow v_{1,2} - 2v_{2,2} = 0\)\nbsp{}; on pose
$v_1 = (1,1)^\top$ et $v_2 = (2,1)^\top$, d'où\nbsp{}:
\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*}
Le système transformé, $\dot y_1(t) = y_1(t) + 1$ et $\dot y_2(t) = 2y_2(t)$,
a pour solution $y_1(t) = \bigl(y_1(0)+1\bigr)e^t - 1$ et
$y_2(t) = y_2(0)\,e^{2t}$. En revenant aux variables initiales,
\(X(t) = H\,Y(t)\)\nbsp{}:
\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*}
avec $y_1(0) = -x_0 + 2y_0$ et $y_2(0) = x_0 - y_0$, où $(x_0,y_0)$ est la
condition initiale. L'état stationnaire est la solution constante
$X^{\star} = -A^{-1}b = (-1,-1)^\top$, et la solution s'écrit encore\nbsp{}:
\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*}
Les deux valeurs propres étant strictement positives, l'écart à l'état
stationnaire croît exponentiellement pour toute condition initiale autre
que \(X^{\star}\)\nbsp{}: l'état stationnaire est un /nœud instable/. Dans le
plan de phase, les trajectoires s'en éloignent tangentiellement à la
direction $v_1$, associée à la plus petite valeur propre, avant de se
redresser vers la direction $v_2$, dont le poids $e^{2t}$ finit par dominer,
comme le montre la figure [[fig:noeud][ci-dessous]].
#+END_details

#+begin_src python :session sedl :exports none :results none
  fig, ax = plt.subplots(figsize=(6, 5))
  grille = np.linspace(-3.5, 1.5, 25)
  U, V = np.meshgrid(grille, grille)
  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'Nœud instable en $X^\star = (-1, -1)$')
  fig.tight_layout()
  fig.savefig("sedl-noeud.svg", transparent=True)
#+end_src

#+CAPTION: *Lignes de courant de l'exercice dans le plan de phase. Les droites en pointillés sont les directions propres : les trajectoires quittent l'état stationnaire tangentiellement à $v_1$ et se redressent vers $v_2$.*
#+LABEL: fig:noeud
[[file:sedl-noeud.svg]]

** Caractérisation des solutions
:PROPERTIES:
:CUSTOM_ID: caracterisation
:END:

Les deux propriétés suivantes généralisent au cas vectoriel la structure
des solutions d'une équation linéaire\nbsp{}: la solution générale de
l'équation complète est la somme d'une solution particulière et de la
solution générale de l'équation homogène, et cette dernière est un espace
vectoriel.

#+BEGIN_property
Si $\Psi(t)$ est une solution particulière du système
$\dot X(t) = A\,X(t) + b(t)$, alors $X(t)$ est solution de ce système si et
seulement si $X(t) - \Psi(t)$ est solution du système homogène
$\dot Z(t) = A\,Z(t)$.
#+END_property

#+BEGIN_proof
Puisque $\Psi(t)$ est une solution de l'équation complète, nous avons
$b(t) = \dot\Psi(t) - A\,\Psi(t)$. La fonction $X(t)$ est solution si et
seulement si $b(t) = \dot X(t) - A\,X(t)$, c'est-à-dire si et seulement si
$\dot X(t) - \dot\Psi(t) = A\bigl(X(t) - \Psi(t)\bigr)$.
#+END_proof

#+BEGIN_property
L'ensemble des solutions réelles (respectivement complexes) du système
linéaire homogène $\dot X(t) = A\,X(t)$, $t\in\mathbb R^+$, est un espace
vectoriel sur $\mathbb R$ (respectivement $\mathbb C$) de dimension $n$.
#+END_property

#+BEGIN_proof
Par linéarité, toute combinaison linéaire de solutions du système homogène
est encore une solution\nbsp{}: l'ensemble des solutions est un espace
vectoriel. Pour $i=1,\dots,n$, soit $X_i(t)$ l'unique solution telle que
$X_i(0) = U_i$, où $U_i$ est le \(i\)-ème vecteur de la base canonique (dont
la \(j\)-ème coordonnée vaut $1$ si $j=i$ et $0$ sinon). Les $n$ fonctions
$X_i(t)$ sont linéairement indépendantes, car leurs valeurs en $t=0$ sont des
vecteurs linéairement indépendants. Soit maintenant une condition initiale
$X_0 = (x_{1,0},\dots,x_{n,0})^\top$ quelconque. Par linéarité,
$\sum_{i=1}^nx_{i,0}X_i(t)$ est une solution du système homogène, et elle
vaut $X_0$ en \(t=0\)\nbsp{}; par unicité, c'est l'unique solution telle que
\(X(0) = X_0\)\nbsp{}:
\begin{equation*}
X(t) = \sum_{i=1}^n x_{i,0}\,X_i(t)
\end{equation*}
Toute solution est donc combinaison linéaire des $X_i(t)$, qui forment une
base de l'espace des solutions[fn:: Voir aussi le chapitre trois de Philippe
Michel, /Cours de mathématiques pour économistes/, Economica, 1989.].
#+END_proof

#+BEGIN_theorem
Soit $H = (v_1,\dots,v_n)$ une matrice de passage donnant une forme réduite
(ou diagonale) de la matrice de transition $A$. La solution générale du
système homogène $\dot X(t) = A\,X(t)$ est\nbsp{}:
\begin{equation*}
X(t) = \sum_{i=1}^n y_i(t)\,v_i \qquad t\in\mathbb R^+
\end{equation*}
où $y_1(t),\dots,y_n(t)$ sont les coordonnées de la solution $Y(t)$ du
système transformé $\dot Y(t) = \hat A\,Y(t)$ telle que $Y(0) = H^{-1}X(0)$.
Chaque $y_i(t)$ est de la forme $e^{\lambda_it}P_i(t)$, où $\lambda_i$ est la
valeur propre associée à la colonne $v_i$ et $P_i$ un polynôme de degré
strictement inférieur à l'ordre de multiplicité de $\lambda_i$.
#+END_theorem

#+BEGIN_proof
Par définition, $X(t) = H\,Y(t)$, ce qui donne la première formule puisque
les colonnes de $H$ sont les vecteurs $v_i$. Pour la forme des $y_i(t)$,
considérons le bloc de $\hat A$ associé à la valeur propre $\lambda_k$ de
multiplicité $r_k$, et résolvons-le par le bas comme dans la section
[[#systeme-transforme][précédente]], sans second membre. La dernière équation
donne $y_{i+r_k}(t) = e^{\lambda_kt}P_{i+r_k}(t)$ avec $P_{i+r_k}$ constant,
égal à $y_{i+r_k}(0)$. En remontant d'une ligne et en substituant cette
solution, on obtient une équation scalaire dont le second membre est le
produit de $e^{\lambda_kt}$ et d'une constante. D'après la formule intégrale,
sa solution est\nbsp{}:
\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*}
où $P_{i+r_k-1}$ est un polynôme de degré au plus un, l'intégrale d'un
polynôme étant un polynôme de degré supérieur d'une unité. Par récurrence
en remontant le bloc, le second membre de la ligne $j$ est le produit de
$e^{\lambda_kt}$ et d'une combinaison linéaire des polynômes des lignes
inférieures, et la solution est de la forme $e^{\lambda_kt}P_j(t)$ où le
degré de $P_j$ excède d'au plus une unité le plus grand degré des lignes
inférieures. La première ligne du bloc, atteinte après $r_k-1$ remontées,
est donc associée à un polynôme de degré au plus $r_k-1$. L'argument
s'applique de la même façon à chaque valeur propre.
#+END_proof

#+BEGIN_remarque
Dans le cas où la matrice de transition est diagonalisable, tous les blocs
sont de dimension un et les polynômes se réduisent à des constantes\nbsp{}: la
solution générale est simplement $X(t) = \sum_{i=1}^n\gamma_ie^{\lambda_it}v_i$,
où les constantes $\gamma_i$ sont les coordonnées de $H^{-1}X(0)$. Si $A$ est
réelle et si $\lambda = a + \mathrm ib$ est une valeur propre complexe de
vecteur propre $v$, alors $\bar\lambda$ et $\bar v$ apparaissent aussi dans
la somme, avec des coefficients conjugués puisque la solution est réelle.
La contribution de la paire est donc $2\,\mathrm{Re}\bigl(\gamma e^{\lambda t}v\bigr)$,
combinaison réelle de $e^{at}\cos(bt)$ et $e^{at}\sin(bt)$, comme dans
l'exemple de la [[#exemple-complexe][section précédente]].
#+END_remarque

#+NAME: ex-symetrique
#+BEGIN_exercice
Soit le système homogène $\dot X(t) = A\,X(t)$ avec\nbsp{}:
\begin{equation*}
A = \begin{pmatrix} 1/2 & 1 & 1\\ 1 & 1/2 & 1\\ 1 & 1 & 1/2\end{pmatrix}
\end{equation*}
Montrez que la solution, étant donnée une condition initiale
$X_0 = (x_{1,0},x_{2,0},x_{3,0})^\top$, est\nbsp{}:
\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*}
pour $t\in\mathbb R^+$.
#+END_exercice

#+ATTR_HTML: :class corrige
#+BEGIN_details
#+BEGIN_summary
Corrigé
#+END_summary
La matrice s'écrit $A = J - \frac12I_3$, où $J$ est la matrice dont tous les
éléments valent un. Puisque $J$ est de rang un, zéro est valeur propre double
de $J$, et la trace de $J$ étant égale à trois, la troisième valeur propre
est $3$. Les valeurs propres de $A$ sont donc $-\frac12$ (double) et
$\frac52$. La matrice est symétrique réelle, donc diagonalisable\nbsp{}: on
doit trouver deux vecteurs propres indépendants pour $-\frac12$, solutions
de $J\,v = 0$, c'est-à-dire $v_1+v_2+v_3 = 0$, et un vecteur propre pour
$\frac52$, solution de \((J-3I_3)\,v = 0\)\nbsp{}:
\begin{equation*}
v_1 = \begin{pmatrix} 1\\ 0\\ -1\end{pmatrix},\qquad
v_2 = \begin{pmatrix} 1\\ -1\\ 0\end{pmatrix}
\qquad\text{et}\qquad
v_3 = \begin{pmatrix} 1\\ 1\\ 1\end{pmatrix}
\end{equation*}
D'après la remarque précédente, la solution générale est\nbsp{}:
\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*}
En évaluant cette expression en $t=0$ et en identifiant avec $X_0$, on
obtient le système $\alpha+\beta+\gamma = x_{1,0}$, $-\beta+\gamma = x_{2,0}$,
$-\alpha+\gamma = x_{3,0}$. La somme des trois équations donne
$3\gamma = x_{1,0}+x_{2,0}+x_{3,0}$, puis $\alpha = \gamma - x_{3,0}$ et
$\beta = \gamma - x_{2,0}$, ce qui est la solution annoncée. La composante
selon $v_3$ domine à long terme, sauf si $x_{1,0}+x_{2,0}+x_{3,0} = 0$,
auquel cas la solution converge vers zéro.
#+END_details

* Équations linéaires d'ordre $n$
:PROPERTIES:
:CUSTOM_ID: ordre-n
:END:

Une équation différentielle faisant intervenir linéairement le niveau d'une
fonction $x$ et ses $n$ premières dérivées est une équation différentielle
linéaire d'ordre $n$. En normalisant le coefficient de la dérivée
\(n\)-ième, nous l'écrivons sous la forme\nbsp{}:

\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*}
où $x^{(p)}$ désigne la dérivée \(p\)-ième de $x$ et où les $a_i$ sont des
paramètres réels. On suppose $a_n\neq0$, sinon il s'agit d'une équation
d'ordre $n-1$ en $\dot x$. On suppose aussi, pour l'instant, que l'équation
n'a pas de second membre.

L'équation d'ordre $n$ est équivalente à un système de $n$ équations du
premier ordre. En posant $x_1 = x^{(n-1)}$, $x_2 = x^{(n-2)}$, \dots,
$x_n = x$, chaque variable est la dérivée de la suivante et la première
vérifie l'équation\nbsp{}:

\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*}
soit, avec $X(t) = \bigl(x_1(t),\dots,x_n(t)\bigr)^\top$, le système
$\dot X(t) = A\,X(t)$ dont la matrice de transition est la /matrice
compagnon/\nbsp{}:

\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
L'ensemble des solutions réelles (complexes) de l'équation différentielle
linéaire homogène d'ordre $n$ est un espace vectoriel sur $\mathbb R$
($\mathbb C$) de dimension $n$. Les $n$ conditions initiales $x(0)$,
$\dot x(0)$, \dots, $x^{(n-1)}(0)$ déterminent une unique solution.
#+END_property

#+BEGIN_proof
Direct par la propriété sur la dimension de l'espace des solutions du
système homogène équivalent, dont la condition initiale $X(0)$ regroupe
précisément les $n$ dérivées successives de $x$ en zéro.
#+END_proof

#+BEGIN_property
Le polynôme caractéristique de la matrice compagnon est\nbsp{}:
\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
Notons $D_n(\lambda) = |A-\lambda I_n|$ et raisonnons par récurrence sur $n$.
Pour $n=1$, $D_1(\lambda) = a_1 - \lambda$ et la formule est vérifiée.
Pour $n\geq2$, développons le déterminant selon la dernière colonne, dont
les seuls éléments non nuls sont $a_n$ en première ligne et $-\lambda$ en
dernière ligne\nbsp{}:
\begin{equation*}
D_n(\lambda) = (-1)^{1+n}a_n\,M_{1,n} + (-1)^{2n}(-\lambda)\,M_{n,n}
\end{equation*}
où $M_{i,j}$ est le mineur obtenu en supprimant la ligne $i$ et la colonne
$j$. Le mineur $M_{n,n}$ est le déterminant de la matrice compagnon d'ordre
$n-1$ associée à $a_1,\dots,a_{n-1}$, diminuée de $\lambda I_{n-1}$, c'est
donc $D_{n-1}(\lambda)$. Le mineur $M_{1,n}$ est le déterminant de la
matrice formée des lignes $2$ à $n$ et des colonnes $1$ à $n-1$ de
$A-\lambda I_n$, qui est triangulaire supérieure avec des $1$ sur la
diagonale (et des $-\lambda$ juste au-dessus)\nbsp{}: il vaut $1$. Ainsi
$D_n(\lambda) = (-1)^{n+1}a_n - \lambda D_{n-1}(\lambda)$, et en substituant
l'hypothèse de récurrence\nbsp{}:
\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

Les valeurs propres de la matrice compagnon sont donc les racines du
polynôme $\lambda^n - a_1\lambda^{n-1} - \dots - a_n$, que l'on obtient
directement en cherchant des solutions de la forme $x(t) = e^{\lambda t}$
dans l'équation d'ordre $n$. On suppose que ce polynôme possède $m\leq n$
racines distinctes $\lambda_k$, de multiplicités $r_k$ avec $\sum_{k=1}^mr_k = n$.

#+BEGIN_property
La solution générale de l'équation différentielle linéaire homogène d'ordre
$n$ est\nbsp{}:
\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*}
où les $n$ constantes $\alpha_{k,p}$ sont déterminées par les conditions
initiales.
#+END_property

#+BEGIN_proof
Direct en appliquant le théorème de la [[#caracterisation][section précédente]]
au système équivalent\nbsp{}: $x = x_n$ est une coordonnée de $X(t)$, donc une
combinaison linéaire des $y_i(t) = e^{\lambda_it}P_i(t)$, et le degré de
$P_i$ est au plus $r_k-1$ lorsque $\lambda_i = \lambda_k$. La famille des $n$
fonctions $t^pe^{\lambda_kt}$ est linéairement indépendante, et l'espace des
solutions étant de dimension $n$, elle en est une base.
#+END_proof

#+BEGIN_exemple
Pour $n=2$, l'équation $\ddot x(t) = a_1\dot x(t) + a_2x(t)$ a pour polynôme
caractéristique $\lambda^2 - a_1\lambda - a_2$. Si le discriminant
$a_1^2+4a_2$ est strictement positif, les racines $\lambda_1$ et $\lambda_2$
sont réelles distinctes et $x(t) = \alpha_1e^{\lambda_1t} + \alpha_2e^{\lambda_2t}$.
S'il est nul, la racine double $\lambda = a_1/2$ donne
$x(t) = (\alpha_0 + \alpha_1t)\,e^{\lambda t}$. S'il est strictement négatif,
les racines sont conjuguées, $\lambda = a\pm\mathrm ib$ avec $a = a_1/2$ et
$b = \frac12\sqrt{-a_1^2-4a_2}$, et la solution réelle s'écrit
$x(t) = e^{at}\bigl(\alpha_1\cos(bt) + \alpha_2\sin(bt)\bigr)$.
#+END_exemple

Considérons maintenant l'équation avec un second membre\nbsp{}:

\begin{equation*}
x^{(n)}(t) = a_1x^{(n-1)}(t) + \dots + a_nx(t) + b(t),\qquad t\in\mathbb R^+
\end{equation*}
où $b(t)$ est une fonction continue de $\mathbb R^+$ dans $\mathbb R$.
L'introduction d'un second membre ne pose pas de problème nouveau\nbsp{}: en
ajoutant à une solution particulière $\psi(t)$ de l'équation complète la
solution générale de l'équation homogène, on obtient la solution générale
de l'équation complète, puisque $x(t)$ est solution de l'équation complète
si et seulement si $x(t)-\psi(t)$ est solution de l'équation homogène. Il ne
reste plus qu'à déterminer les $n$ constantes d'intégration à l'aide des
conditions initiales. Pour les seconds membres usuels, on connaît la forme
de la solution particulière.

#+BEGIN_property
Si le second membre est de la forme $b(t) = e^{\mu t}Q(t)$, où $\mu$ est un
nombre réel ou complexe et $Q$ un polynôme de degré $q$, alors l'équation
complète admet une solution particulière de la forme $e^{\mu t}R(t)$, avec
$R$ un polynôme de degré $q$, si $\mu$ n'est pas racine du polynôme
caractéristique, et de la forme $t^re^{\mu t}R(t)$ si $\mu$ est racine
d'ordre $r$ du polynôme caractéristique.
#+END_property

Nous admettons cette propriété, dont l'exemple de la section
[[#exemple-complexe][précédente]] est une illustration avec \(\mu = 0\)\nbsp{}: le
second membre y était polynomial de degré un, zéro n'était pas valeur
propre, et la solution particulière était polynomiale de degré un. Les
coefficients de $R$ s'obtiennent par identification, en substituant la forme
postulée dans l'équation.

* Stabilité d'un système dynamique
:PROPERTIES:
:CUSTOM_ID: stabilite
:END:

** Définition et critère
:PROPERTIES:
:CUSTOM_ID: critere
:END:

On s'intéresse au système homogène\nbsp{}:

\begin{equation*}
\dot X(t) = A\,X(t),\qquad t\in\mathbb R^+
\end{equation*}
#+BEGIN_definition
On dit que la matrice réelle $A$ de dimension $n\times n$ est /d-stable/
(différentiellement stable) si toute solution $X(t)$ du système homogène
tend vers $0$ quand $t$ tend vers $+\infty$. On dit alors que l'état
stationnaire nul du système est /globalement stable/.
#+END_definition

Plus généralement, considérons une dynamique de la forme
$\dot X(t) = A\,X(t) + b$, où le second membre $b$ est un vecteur constant.
Un /état stationnaire/ est une solution constante $X^{\star}$, c'est-à-dire
un point de $\mathbb R^n$ vérifiant $A\,X^{\star} + b = 0$, soit
$X^{\star} = -A^{-1}b$ si $A$ est de plein rang. L'état stationnaire est
globalement stable si, pour toute condition initiale $X_0\in\mathbb R^n$, la
solution telle que $X(0) = X_0$ tend vers $X^{\star}$ quand $t$ tend vers
$+\infty$. Puisque l'écart $Z(t) = X(t) - X^{\star}$ vérifie
$\dot Z(t) = A\,Z(t)$, l'état stationnaire $X^{\star}$ est globalement stable
si et seulement si $A$ est d-stable.

#+BEGIN_theorem
Pour que la matrice $A$ soit d-stable, il faut et il suffit que toutes ses
valeurs propres soient de partie réelle strictement négative.
#+END_theorem

#+BEGIN_proof
On établit successivement la condition nécessaire et la condition
suffisante.

/Condition nécessaire./ On suppose que $A$ est d-stable. Soit
$\lambda = a + \mathrm ib$ une valeur propre de $A$ et $v\in\mathbb C^n$ un
vecteur propre associé, non nul. Alors $Z(t) = e^{\lambda t}v$ est une
solution complexe du système, car
$\dot Z(t) = \lambda e^{\lambda t}v = e^{\lambda t}A\,v = A\,Z(t)$. Ses parties
réelle et imaginaire sont des solutions réelles, qui tendent vers zéro
puisque $A$ est d-stable, donc $Z(t)$ tend vers zéro. Soit $v_i$ une
coordonnée non nulle de $v$. La coordonnée $z_i(t) = v_ie^{\lambda t}$ tend
vers zéro, ce qui implique que $|z_i(t)| = |v_i|\,e^{at}$ tend vers zéro
quand $t$ tend vers l'infini, donc que $a$ est strictement négatif.

/Condition suffisante./ On suppose que toutes les valeurs propres
$\lambda_k = a_k + \mathrm ib_k$ de $A$, de multiplicités $r_k$, sont de parties
réelles $a_k$ strictement négatives ($k=1,\dots,m$). Soit $X(t)$ une solution
du système. D'après le théorème de la section
[[#caracterisation][caractérisation des solutions]], chaque coordonnée
$x_i(t)$ de $X(t)$ est de la forme $\sum_{k=1}^me^{\lambda_kt}P_{i,k}(t)$, où
les $P_{i,k}$ sont des polynômes à coefficients complexes, et\nbsp{}:
\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*}
Puisque $a_k < 0$, chaque terme $e^{a_kt}|P_{i,k}(t)|$ tend vers zéro quand
$t$ tend vers l'infini, l'exponentielle l'emportant sur tout polynôme. Ainsi
$X(t)$ tend vers zéro et la matrice $A$ est d-stable.
#+END_proof

** Cas particuliers
:PROPERTIES:
:CUSTOM_ID: cas-particuliers
:END:

Le critère précédent est simple à mettre en œuvre en petite dimension.

- Si $n=1$, $A = (a)$ est d-stable si et seulement si $a < 0$.
- Si $n=2$, $A = \begin{pmatrix} a & b\\ c & d\end{pmatrix}$ est d-stable si
  et seulement si $\mathrm{tr}(A) = a + d < 0$ et $|A| = ad - bc > 0$. En
  effet, la trace est la somme des valeurs propres et le déterminant leur
  produit. Si les valeurs propres sont réelles, elles sont toutes deux
  strictement négatives si et seulement si leur somme est strictement
  négative et leur produit strictement positif. Si elles sont complexes
  conjuguées, $a\pm\mathrm ib$ avec $b\neq0$, leur produit $a^2+b^2$ est
  toujours strictement positif et leur somme $2a$ est strictement négative si
  et seulement si leur partie réelle l'est. La figure [[fig:stabilite][ci-dessous]]
  illustre les trois configurations possibles lorsque le déterminant est non
  nul.
- Si $A$ est symétrique, ses valeurs propres sont réelles et elle est
  d-stable si et seulement si elle est définie négative, ce qui d'après le
  critère de Sylvester s'écrit\nbsp{}:
  \begin{equation*}
  (-1)^p\,|A_{1:p}| > 0\qquad\text{pour tout } 1\leq p\leq n
  \end{equation*}
  où $A_{1:p}$ est la sous-matrice formée des $p$ premières lignes et
  colonnes de $A$ (les mineurs principaux dominants alternent en signe, en
  commençant par $a_{1,1} < 0$).

#+begin_src python :session sedl :exports none :results none
  fig, axes = plt.subplots(1, 3, figsize=(12, 4))
  cas = ((np.array([[-2, 1], [1, -2]]), r'Nœud stable : tr $=-4$, dét $=3$'),
         (np.array([[-0.5, -2], [2, -0.5]]), r'Foyer stable : tr $=-1$, dét $=4{,}25$'),
         (np.array([[1, 2], [2, 1]]), r'Col : tr $=2$, dét $=-3$'))
  grille = np.linspace(-2, 2, 20)
  U, V = np.meshgrid(grille, grille)
  for ax, (A, titre) in zip(axes, cas):
      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(titre)
  fig.tight_layout()
  fig.savefig("sedl-stabilite.svg", transparent=True)
#+end_src

#+CAPTION: *Portraits de phase de trois systèmes homogènes de dimension deux. À gauche et au centre la trace est négative et le déterminant positif, toutes les trajectoires convergent vers l'origine, sans oscillation (valeurs propres réelles) ou en spiralant (valeurs propres complexes). À droite le déterminant est négatif, les valeurs propres sont réelles et de signes opposés et l'origine est un col : seules les trajectoires issues de la direction propre stable convergent.*
#+LABEL: fig:stabilite
[[file:sedl-stabilite.svg]]

Au-delà de la dimension deux, on peut utiliser le polynôme caractéristique
sans calculer ses racines. La proposition suivante traite le cas $n=3$.

#+BEGIN_proposition
Pour que le polynôme $\lambda^3 + a\lambda^2 + b\lambda + c$ ait toutes ses
racines de partie réelle strictement négative, il faut et il suffit que ses
coefficients vérifient $a > 0$, $b > 0$, $c > 0$ et $ab > c$.
#+END_proposition

#+BEGIN_proof
Notons $\lambda_1$, $\lambda_2$ et $\lambda_3$ les trois racines dans
$\mathbb C$. En développant $(\lambda-\lambda_1)(\lambda-\lambda_2)(\lambda-\lambda_3)$
et en identifiant, on a $a = -(\lambda_1+\lambda_2+\lambda_3)$,
$b = \lambda_1\lambda_2+\lambda_1\lambda_3+\lambda_2\lambda_3$ et
$c = -\lambda_1\lambda_2\lambda_3$, d'où l'identité\nbsp{}:
\begin{equation*}
ab - c = -(\lambda_1+\lambda_2)(\lambda_2+\lambda_3)(\lambda_3+\lambda_1)
\end{equation*}
/Condition nécessaire./ Si les trois racines sont réelles et strictement
négatives, on a immédiatement $a > 0$, $b > 0$, $c > 0$, et chaque facteur
$\lambda_i+\lambda_j$ étant strictement négatif, $ab - c > 0$. Si
$\lambda_1 = -r$ est réelle, avec $r > 0$, et $\lambda_{2,3} = \alpha\pm\mathrm i\beta$
avec $\alpha < 0$, alors $a = r - 2\alpha > 0$,
$b = -2r\alpha + \alpha^2+\beta^2 > 0$, $c = r(\alpha^2+\beta^2) > 0$ et,
puisque $(\lambda_1+\lambda_2)(\lambda_1+\lambda_3) = |\lambda_1+\lambda_2|^2 > 0$
et $\lambda_2+\lambda_3 = 2\alpha < 0$, on a encore $ab - c > 0$.

/Condition suffisante./ Supposons $a > 0$, $b > 0$, $c > 0$ et $ab > c$. Un
polynôme réel de degré trois possède au moins une racine réelle
$\lambda_1$, et celle-ci est strictement négative puisque, les coefficients
étant tous strictement positifs, le polynôme est strictement positif pour
$\lambda\geq0$. On factorise
$\lambda^3+a\lambda^2+b\lambda+c = (\lambda-\lambda_1)(\lambda^2+s\lambda+q)$, et
par identification $c = -\lambda_1q$, soit $q = -c/\lambda_1 > 0$, puis
$b = q - \lambda_1s$ et $a = s - \lambda_1$. Il vient alors\nbsp{}:
\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*}
et donc $s = (ab-c)/(b+\lambda_1^2) > 0$. Les deux autres racines sont celles
du facteur $\lambda^2+s\lambda+q$, dont la somme $-s$ est strictement négative
et le produit $q$ strictement positif\nbsp{}: d'après le cas $n=2$, elles sont
de partie réelle strictement négative.
#+END_proof

Les deux propriétés suivantes donnent des conditions nécessaires de
stabilité, faciles à vérifier, qui permettent d'écarter rapidement une
matrice.

#+BEGIN_property
Si la matrice réelle $A$ de dimension $n\times n$ est d-stable, alors sa
trace est strictement négative et son déterminant est du signe de $(-1)^n$.
#+END_property

#+BEGIN_proof
La trace est la somme des valeurs propres et le déterminant leur produit.
Les valeurs propres réelles sont strictement négatives, et les valeurs
propres complexes vont par paires conjuguées $a\pm\mathrm ib$ avec $a < 0$,
dont la somme $2a$ est strictement négative et le produit $a^2+b^2$
strictement positif. La somme de toutes les valeurs propres est donc
strictement négative, et leur produit est du signe de $(-1)^{n_r}$, où $n_r$
est le nombre de valeurs propres réelles, qui a la parité de $n$.
#+END_proof

#+BEGIN_property
Si le polynôme $\lambda^n + a_1\lambda^{n-1} + \dots + a_{n-1}\lambda + a_n$ a
toutes ses racines de partie réelle strictement négative, alors tous ses
coefficients sont strictement positifs.
#+END_property

#+BEGIN_proof
Le polynôme se factorise en $(\lambda-\lambda_1)(\lambda-\lambda_2)\cdots(\lambda-\lambda_n)$.
Si $\lambda_1 < 0$ est réelle, les coefficients de $\lambda-\lambda_1$ sont
strictement positifs. Si $\lambda_1 = a + \mathrm ib$ avec $a < 0$ et
$b\neq0$, on regroupe le facteur avec celui de la racine conjuguée\nbsp{}:
$(\lambda - a - \mathrm ib)(\lambda - a + \mathrm ib) = \lambda^2 - 2a\lambda + a^2+b^2$,
dont les coefficients sont strictement positifs. Le polynôme est donc le
produit de polynômes du premier et du second degré à coefficients
strictement positifs, et tous ses coefficients sont strictement positifs.
#+END_proof

#+BEGIN_remarque
Le polynôme caractéristique de la matrice compagnon de l'équation d'ordre
$n$ est $(-1)^n(\lambda^n - a_1\lambda^{n-1} - \dots - a_n)$. La propriété
s'applique au polynôme unitaire \(\lambda^n - a_1\lambda^{n-1} - \dots - a_n\)\nbsp{}:
une condition nécessaire de stabilité de l'équation d'ordre $n$ est que tous
les coefficients $a_i$ soient strictement négatifs. Pour $n=2$, la condition
$a_1 < 0$ et $a_2 < 0$ est aussi suffisante, puisque la trace de la matrice
compagnon est $a_1$ et son déterminant $-a_2$.
#+END_remarque

** Matrices à diagonale dominante
:PROPERTIES:
:CUSTOM_ID: diagonale-dominante
:END:

#+BEGIN_definition
Une matrice carrée $A = (a_{i,j})$ de dimension $n\times n$ à éléments réels
est à /diagonale négative dominante/ s'il existe $n$ nombres $d_i > 0$,
$i=1,\dots,n$, tels que\nbsp{}:
\begin{equation*}
-d_i\,a_{i,i} > \sum_{j\neq i}d_j\,|a_{i,j}|\qquad\text{pour tout } 1\leq i\leq n
\end{equation*}
#+END_definition

En particulier, les éléments diagonaux d'une telle matrice sont strictement
négatifs, et chacun domine, en valeur absolue, les autres éléments de sa
ligne pondérés par les $d_j$. Cette condition est facile à vérifier, et elle
suffit à garantir la stabilité.

#+BEGIN_property
Une matrice à diagonale négative dominante est d-stable.
#+END_property

#+BEGIN_proof
Soit $\lambda = \alpha + \mathrm i\beta$ une valeur propre de $A$ et
$v\in\mathbb C^n$, non nul, un vecteur propre associé, de coordonnées $v_j$.
Choisissons un indice $k$ tel que \(|v_k|/d_k = \max_j|v_j|/d_j\)\nbsp{}; on a
$v_k\neq0$ puisque $v$ est non nul. La \(k\)-ième ligne de $A\,v = \lambda v$
s'écrit $(\lambda - a_{k,k})\,v_k = \sum_{j\neq k}a_{k,j}\,v_j$, d'où, en
utilisant $|v_j|\leq d_j|v_k|/d_k$ pour tout \(j\)\nbsp{}:
\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*}
où la dernière inégalité est la définition de la diagonale négative
dominante. En divisant par $|v_k| > 0$ il vient $|\lambda - a_{k,k}| < -a_{k,k}$,
et puisque la partie réelle d'un nombre complexe est majorée par son module,
$\alpha - a_{k,k} \leq |\lambda - a_{k,k}| < -a_{k,k}$, soit $\alpha < 0$. Toutes
les valeurs propres sont de partie réelle strictement négative et $A$ est
d-stable.
#+END_proof

** Le critère de Routh–Hurwitz
:PROPERTIES:
:CUSTOM_ID: routh-hurwitz
:END:

Le critère de Routh–Hurwitz généralise la proposition établie pour $n=3$ en
donnant une condition nécessaire et suffisante portant sur les coefficients
du polynôme caractéristique, sans calcul des racines. Soit
$p(\lambda) = \lambda^n + a_1\lambda^{n-1} + \dots + a_{n-1}\lambda + a_n$ un
polynôme unitaire à coefficients réels[fn:: Pour une matrice $A$, on prend
$p(\lambda) = (-1)^n\chi_A(\lambda)$, qui a les mêmes racines.]. On construit
la matrice carrée de dimension \(n\times n\)\nbsp{}:

\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*}
dont l'élément de la ligne $i$ et de la colonne $j$ est $a_{2i-j}$, avec la
convention $a_0 = 1$ et $a_k = 0$ si $k < 0$ ou $k > n$. On peut montrer
qu'une condition nécessaire et suffisante pour que toutes les racines de
$p$ soient de partie réelle strictement négative est que les $n$ mineurs
principaux dominants de $R$ (les déterminants des sous-matrices formées des
$p$ premières lignes et colonnes, $p=1,\dots,n$) soient tous strictement
positifs. Nous admettons ce résultat, et nous contentons de vérifier qu'il
redonne les conditions déjà établies.

#+NAME: ex-rh2
#+BEGIN_exercice
Écrivez la matrice de Routh–Hurwitz du polynôme $\lambda^2 + a_1\lambda + a_2$
et retrouvez la condition de stabilité d'un système de dimension deux.
#+END_exercice

#+ATTR_HTML: :class corrige
#+BEGIN_details
#+BEGIN_summary
Corrigé
#+END_summary
Avec $n=2$, $a_3 = 0$ et la matrice est\nbsp{}:
\begin{equation*}
R = \begin{pmatrix} a_1 & 0\\ 1 & a_2\end{pmatrix}
\end{equation*}
Les mineurs principaux dominants sont $a_1$ et $a_1a_2$. Les racines sont
de partie réelle strictement négative si et seulement si\nbsp{}:
\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*}
Pour une matrice $A$ de dimension deux, $p(\lambda) = \lambda^2 - \mathrm{tr}(A)\,\lambda + |A|$,
et l'on retrouve la condition $\mathrm{tr}(A) < 0$ et $|A| > 0$.
#+END_details

#+NAME: ex-rh3
#+BEGIN_exercice
Écrivez la matrice de Routh–Hurwitz du polynôme
$\lambda^3 + a_1\lambda^2 + a_2\lambda + a_3$ et retrouvez la proposition
établie plus haut pour les polynômes de degré trois.
#+END_exercice

#+ATTR_HTML: :class corrige
#+BEGIN_details
#+BEGIN_summary
Corrigé
#+END_summary
Avec $n=3$, $a_4 = a_5 = 0$ et la matrice est\nbsp{}:
\begin{equation*}
R = \begin{pmatrix} a_1 & a_3 & 0\\ 1 & a_2 & 0\\ 0 & a_1 & a_3\end{pmatrix}
\end{equation*}
Les mineurs principaux dominants sont $a_1$, $a_1a_2 - a_3$ et, en
développant selon la dernière colonne, $a_3(a_1a_2 - a_3)$. Les conditions
de Routh–Hurwitz s'écrivent donc $a_1 > 0$, $a_1a_2 > a_3$ et
$a_3(a_1a_2-a_3) > 0$, soit $a_3 > 0$ compte tenu de la deuxième. Enfin
$a_1 > 0$, $a_3 > 0$ et $a_1a_2 > a_3$ impliquent $a_2 > 0$. On retrouve
exactement les conditions $a_1 > 0$, $a_2 > 0$, $a_3 > 0$ et $a_1a_2 > a_3$ de
la proposition.
#+END_details