Équation différentielle linéaire/Généralisation, notation matricielle

Début de la boite de navigation du chapitre

fin de la boite de navigation du chapitre

Introduction

Dans les cas réels, et notamment en mécanique, il est courant d’avoir affaire à des équations différentielles d'ordre 2 ou plus. Le but de ce chapitre est de montrer comment, grâce à l'algèbre et à la notation matricielle, tout problème d'équation différentielle linéaire ordinaire à coefficients constants peut être ramené au cas « simple » d'une équation différentielle d'ordre un, traité dans le chapitre précédent.

En particulier, cela amène des résultats très généraux, comme l’existence de solutions, ou la dimension de l'espace des solutions, tout en proposant une méthode de résolution efficace.

Objet d'étude

Rappelons l’objet d'étude de ce chapitre :

Définition

Soient $a_{0},a_{1},a_{2},\ldots ,a_{n}$ et C des complexes, a_n étant non-nul. Soit f une fonction n fois dérivable sur $\scriptstyle \mathbb {R}$, à valeur dans $\scriptstyle \mathbb {R}$ ou $\scriptstyle \mathbb {C}$. On appelle équation différentielle linéaire ordinaire d'ordre n (à coefficients constants) la relation :

$\forall x\in \mathbb {R} ,\quad a_{n}f^{(n)}(x)+\ldots +a_{2}f^{(2)}(x)+a_{1}f^{(1)}(x)+a_{0}f(x)+C=0$

Puisque cela nous sera utile par la suite, introduisons dès maintenant quelques notations :

Notations

Les valeurs propres d'une matrice sont notées $\lambda _{i}$. Le vecteur propre assocé à $\lambda _{i}$ est noté ${\vec {\Lambda }}_{i}$. Le polynôme caractéristique d'une matrice M, défini comme le déterminant de $\left(M-X\cdot I_{n}\right)$, est noté $\chi _{M}$.

Mise en forme matricielle

Réduction de l’ordre

Grâce à la notation matricielle, nous allons réécrire cette relation sous la forme d'une simple équation du premier ordre. Pour cela, on va augmenter le nombre de dimensions. Nous ne travaillerons alors plus avec des scalaires, mais avec des vecteurs et des matrices.

Il est commode dans un premier temps de ne considérer que le cas homogène. La solution générale s'obtiendra aisément en ajoutant une solution particulière (constante) à la solution homogène. Pour simplifier encore davantage les calculs, divisons tous les coefficients par a_n (celui-ci étant non-nul) sans changer leur nom. On travaille donc sur l'équation différentielle :

$f^{(n)}+\ldots +a_{2}f^{(2)}+a_{1}f^{(1)}+a_{0}f=0$

Posons tout d’abord le vecteur X :

$X={\begin{pmatrix}f^{(n-1)}\\\vdots \\f^{(2)}\\f^{(1)}\\f\end{pmatrix}}$

Ce vecteur est dérivable, et on a bien sûr :

$X'={\begin{pmatrix}f^{(n)}\\\vdots \\f^{(3)}\\f^{(2)}\\f^{(1)}\end{pmatrix}}$

L'idée est alors de réécrire l'équation différentielle de départ, comme somme de X’ et de X : on aura ainsi une équation différentielle (matricielle) linéaire d'ordre un ! On cherche donc à écrire : $X'=MX$ Avec M une matrice n × n, dont les coefficients doivent être bien choisis. Tentons de les retrouver.

${\begin{pmatrix}f^{(n)}\\\vdots \\f^{(3)}\\f^{(2)}\\f^{(1)}\end{pmatrix}}={\begin{pmatrix}&&&\\&&?&&\\&&&\end{pmatrix}}{\begin{pmatrix}f^{(n-1)}\\\vdots \\f^{(2)}\\f^{(1)}\\f\end{pmatrix}}$

Travaillons ligne par ligne. La première ligne porte simplement les coefficients -a_n-1, ..., -a₀ dans cet ordre, de manière triviale. La seconde ligne peut être facilement résolue en mettant un coefficient 1 pour la dérivée (n+1)-ième, et 0 partout ailleurs. Il en est de même pour les autres. Finalement, on obtient :

$M={\begin{pmatrix}-a_{n-1}&\cdots &-a_{1}&-a_{0}\\1&\cdots &0&0\\\vdots &\ddots &\vdots &\vdots \\0&\cdots &1&0\end{pmatrix}}$

On peut remarquer le fait que M est une matrice compagnon.

Propriétés de la matrice M

Nous allons notamment nous intéresser aux valeurs et vecteurs propres de M. Le déterminant de cette matrice est le produit de ses valeurs propres, la trace est leur somme. Le déterminant et la trace de M sont respectivement :

$\det M=\prod _{i}\lambda _{i}=a_{0}\left(-1\right)^{n}$

$\mathrm {tr} \,M=\sum _{i}\lambda _{i}=-a_{n-1}$

Les valeurs propres sont les racines du polynôme caractéristique (au signe près) :

$\chi _{M}=(-1)^{n}[X^{n}+a_{n-1}X^{n-1}+\ldots +a_{2}X^{2}+a_{1}X+a_{0}]$

La résolution explicite, pour des équations d'ordre supérieur à 3, est difficile à la main, sauf cas particuliers. Néanmoins, l’utilisation de méthodes heuristiques (pour trouver une racine « triviale ») réduit le degré du problème — et il est toujours envisageable d’utiliser des logiciels de calcul formel (Mathematica, Matlab, Maple…), qui trouveront rapidement les racines.

En fait, M est toujours diagonalisable, que ses valeurs propres soient distinctes (auquel cas on peut diagonaliser M par une matrice de Vandermonde), ou pas. On a donc une matrice inversible P, telle que :

$P={\begin{pmatrix}{\vec {\Lambda }}_{n-1}&\cdots &{\vec {\Lambda }}_{1}&{\vec {\Lambda }}_{0}\end{pmatrix}}$

$PMP^{-1}=\Delta ={\begin{pmatrix}\lambda _{n-1}&&&\\&\ddots &&\\&&\lambda _{1}&\\&&&\lambda _{0}\end{pmatrix}}$

Résolution de l'équation différentielle

Forme de la solution

Repartons de notre équation différentielle matricielle :

$X'=MX$

Par analogie avec le cas simple décrit dans le chapitre précédent, une solution est de la forme :

$X=A\exp \left(Mx\right)$

La solution à l'équation de départ, non-homogène, s'obtient facilement en rajoutant une constante adaptée.

Mais qu'est-ce qu'une exponentielle de matrices ?

Exponentielle de matrices

Il existe deux manières d'approcher l'exponentielle de matrices, qui sont complémentaires : soit d'un point de vue géométrique, soit d'un point de vue analytique. Nous choisissons ici la première approche, probablement plus simple.

Rappelons tout d’abord une relation très importante :

Propriété

Soient $P\in \mathrm {GL} _{n}\left(\mathbb {C} \right)$ et $A\in {\mathfrak {M}}_{n}\left(\mathbb {C} \right)$. Alors

$\exp \left(P^{-1}AP\right)=P^{-1}\left(\exp A\right)P{\text{.}}$

Dans le cas d'une matrice diagonalisable, il est même possible d'expliciter complètement l'exponentielle de matrice, ce qui nous est fort utile ici. On a, en effet, le théorème suivant :

Début d’un théorème

Théorème

$\exp \left(Mx\right)=\sum _{i}e^{\lambda _{i}x}{\vec {\Lambda }}_{i}{\text{.}}$

Fin du théorème

La méthode alternative, mais qui ne permet pas toujours une résolution manuelle, est d’utiliser la série entière de l'exponentielle :

$\exp :X\mapsto e^{X}:=\sum _{k=0}^{\infty }{X^{k} \over k!}.$

Puisqu'on est capable de calculer toutes les puissances de M, l’utilisation de cette somme est envisageable. Elle est en réalité plus pratique dans les cas d'évaluation numérique : il n'est plus besoin de considérer les éventuelles valeurs propres ou vecteurs propres de la matrice.

Solution

On obtient, d'un seul coup, la solution pour la fonction f et pour toutes ses dérivées. Faisons rentrer le facteur A de la solution homogène dans l’expression développée de l'exponentielle, et prenons en compte la quantité constante non-homogène. La solution finale à l'équation différentielle ordinaire linéaire à coefficients constants d'ordre n de départ est alors :

$X\left(x\right)=\sum _{i}\left(A_{i}e^{\lambda _{i}x}+B_{i}\right){\vec {\Lambda }}_{i}$

Équations à coefficients non-constants

Échec des théorèmes d'ordre un

Le travail effectué jusqu'ici n'est valable que si les coefficients de l'équation différentielle sont des complexes. Peut-on généraliser ces résultats aux cas où il s'agit de fonctions de x ? La réponse, malheureusement, est non. Pire : il est impossible d'exprimer la solution, si solution il y a, par une formule finie (on parle de résolution en forme close). En fait, seules les équations d'ordre un (décrites dans le chapitre 1) et deux (décrites dans le chapitre 5) admettent des solutions analytiques.

En reprenant les notations du paragraphe précédent, l'équation différentielle homogène d'ordre n peut s'écrire sous la forme d'une équation matricielle d'ordre un :

$X'\left(x\right)+M\left(x\right)X\left(x\right)=0$

avec l’application matricielle :

$M:x\mapsto {\begin{pmatrix}a_{n-1}(x)&\cdots &a_{1}(x)&a_{0}(x)\\-1&\cdots &0&0\\\vdots &\ddots &\vdots &\vdots \\0&\cdots &-1&0\end{pmatrix}}$

Devrons-nous pour autant renoncer à leur résolution ?

Les séries de Magnus

Il existe un moyen de résoudre ces équations : les séries de Magnus. L’idée est que, par analogie avec le cas à coefficients constants, on va chercher X sous la forme d'une exponentielle de matrices :

$X(x)=e^{\Omega (x)}$

avec $\Omega$ une certaine application matricielle à déterminer.

Cette application est solution de quelques relations, qui permettent de déduire qu'elle peut être définie par une série, appelée série de Magnus :

$\Omega (t)=\int _{0}^{t}M(\tau )\,\mathrm {d} \tau -{\frac {1}{2}}\int _{0}^{t}\int _{0}^{\tau _{1}}[M(\tau _{2}),M(\tau _{1})]\,\mathrm {d} \tau _{2}\mathrm {d} \tau _{1}+{\frac {1}{4}}\int _{0}^{t}\int _{0}^{\tau _{1}}\int _{0}^{\tau _{2}}[[M(\tau _{3}),M(\tau _{2})],M(\tau _{1})]\,\mathrm {d} \tau _{3}\mathrm {d} \tau _{2}\mathrm {d} \tau _{1}+{\frac {1}{12}}\ldots$

avec $\left[\cdot ,\cdot \right]:A,B\mapsto AB-BA$. Il s'agit d'une somme infinie.

Début d’un théorème

Théorème de Magnus

Si $\Omega$ converge, alors $X(x)=e^{\Omega (x)}$ est solution de l'équation différentielle homogène à coefficients constants.

Fin du théorème

La question est, justement, de savoir quand cette série converge. Cette question n’est pas à l’heure actuelle (2007) complètement résolue. Toutefois, de nombreux résultats existent, parmi lesquels cette condition suffisante (mais pas nécessaire) :

Début d’un théorème

Condition suffisante de convergence

Si $\int _{0}^{x}\|M(t)\|_{2}\,\mathrm {d} t<\pi$, alors $\Omega (x)$ converge.

Fin du théorème

Propriétés

Stabilité analytique

La question de la stabilité, évoquée dans le chapitre précédent pour les équations du premier ordre, est en fait très générale.

Définition

Une équation différentielle est dite stable analytiquement lorsque toutes ses solutions éventuelles sont bornées et convergentes. Elle est dite instable dans le cas contraire.

La stabilité d'une équation différentielle linéaire ordinaire est liée aux valeurs propres de la matrice associée :

Début d’un théorème

Stabilité analytique d'une ODE linéaire

Une équation différentielle ordinaire linéaire est :

• analytiquement ω-stable si et seulement si, pour toute valeur propre λ, Re(λ) < 0 ;
• analytiquement instable dans le cas contraire ;
  • α-instable si et seulement si, il existe une valeur propre λ telle que Re(λ) > 0 ;

Fin du théorème

Remarques :

• dans le cas d'une équation ω-stable, le polynôme caractéristique du problème est un polynôme d'Hurwitz ;
• dans le cas d'une équation α-instable, il s'agit d'un polynôme de Schur. Les cas α-instables seraient ω-stables si on inversait le sens du temps ;
• ce théorème de stabilité se généralise, moyennant quelques conditions supplémentaires, aux cas d'équations différentielles non-linéaires ;
• le cas où toutes les valeurs propres sont de partie réelle nulle n’est pas analytiquement stable parce que les solutions obtenues oscillent et ne convergent donc pas — elles sont néanmoins bornées, et vérifient une autre forme de stabilité, moins restrictive, appelée « stabilité EBSB » (pour « entrée bornée, sortie bornée »).

Stabilité numérique

Une équation standard pour définir et évaluer la stabilité d'une méthode numérique est $y'=ky$, dont la solution analytique est $y(t)=e^{kt}$. La méthode d'Euler pour une telle équation donne :

$y_{n+1}=y_{n}+hf(t_{n},y_{n})=y_{n}+h(ky_{n})=y_{n}+hky_{n}=(1+hk)y_{n}.$

Introduisons une notion au travers de cette méthode, qui est en réalité plus générale, mais n'est ici qu’illustrative.

Définitions

Région de stabilité de la méthode d'Euler.

Ce qui se ramène à $y_{n}=(1+hk)^{n}y_{0}$, résultat partagé en fait par toutes les méthodes de Runge-Kutta. Posons la fonction $\phi :z\mapsto 1+z$, alors $y_{n}=\phi (hk)^{n}y_{0}$.

La fonction $\phi$ est appelée fonction de stabilité de la méthode d'Euler. Trivialement, la méthode donne un résultat borné si et seulement si cette fonction est de module inférieur ou égal à 1, et la limite est nulle si et seulement si ce module est inférieur strictement à 1. L'ensemble des complexes pour lesquels cette dernière condition est vérifiée est appelée région de stabilité.

On peut alors donner deux définitions de la stabilité numérique, par ordre croissant d'efficacité (et de rareté) :

• Une méthode d'intégration numérique qui, pour toute équation différentielle convergente, converge également, est dite A-stable.
• Une méthode d'intégration numérique dont la fonction de stabilité tend vers 0 lorsque son argument devient grand en module est dite L-stable.

Alors, on peut discuter de la stabilité d'une méthode, sachant le théorème suivant :

Début d’un théorème

Condition suffisante de A-stabilité

Une méthode d'intégration numérique est A-stable lorsque le demi-plan des complexes de partie négative $\scriptstyle \{z\in \mathbb {C} |Re(z)<0\}$ est compris dans sa région de stabilité.

Fin du théorème

Pour la méthode d'Euler, la région de stabilité se réduit à un cercle : il ne s'agit donc pas d'une méthode A-stable. Néanmoins, il est toujours possible (lorsque k < 0) de trouver un pas h tel que la fonction de stabilité soit de module inférieur à 1, et donc que l'intégration numérique converge.