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 :

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

Posons tout d’abord le vecteur X :

${\displaystyle 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 :

${\displaystyle 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 : ${\displaystyle X'=MX}$ Avec M une matrice n × n, dont les coefficients doivent être bien choisis. Tentons de les retrouver.

${\displaystyle {\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 :

${\displaystyle 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 :

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

${\displaystyle \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) :

${\displaystyle \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 :

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

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

Forme de la solution

Repartons de notre équation différentielle matricielle :

${\displaystyle X'=MX}$

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

${\displaystyle 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 ${\displaystyle P\in \mathrm {GL} _{n}\left(\mathbb {C} \right)}$ et ${\displaystyle A\in {\mathfrak {M}}_{n}\left(\mathbb {C} \right)}$. Alors

${\displaystyle \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 :

Théorème

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

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 :

${\displaystyle \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 :

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

É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 :

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

avec l’application matricielle :

${\displaystyle 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 :

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

avec ${\displaystyle \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 :

${\displaystyle \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 ${\displaystyle \left[\cdot ,\cdot \right]:A,B\mapsto AB-BA}$. Il s'agit d'une somme infinie.

Théorème de Magnus

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

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) :

Condition suffisante de convergence

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

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 :

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 ;

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 ${\displaystyle y'=ky}$, dont la solution analytique est ${\displaystyle y(t)=e^{kt}}$. La méthode d'Euler pour une telle équation donne :

${\displaystyle 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 à ${\displaystyle y_{n}=(1+hk)^{n}y_{0}}$, résultat partagé en fait par toutes les méthodes de Runge-Kutta. Posons la fonction ${\displaystyle \phi :z\mapsto 1+z}$, alors ${\displaystyle y_{n}=\phi (hk)^{n}y_{0}}$.

La fonction ${\displaystyle \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 :

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 ${\displaystyle \scriptstyle \{z\in \mathbb {C} |Re(z)<0\}}$ est compris dans sa région de stabilité.

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.