EDO et problèmes raides

Partons d'un exemple simple...

L'équation différentielle

Commençons par examiner un problème très simple de premier ordre, que vous rencontrerez sans doute dans vos cours en prépa ou à la fac (il est super classique).

Considérons le problème de Cauchy suivant :

\( \left\{ \begin{array}{l} \dfrac{dy}{dt} = -150y + 30 \\ y(0) = 1/5  \end{array} \right. \)

pour t dans l'intervalle [0,1]
Ce problème admet une solution et une seule, qui est d'ailleurs très facile à calculer. Mais pour s'amuser, nous allons demander à Python, et plus particulièrement à son module sympy de faire le calcul à notre place.

Sa solution analytique

Posons d'abord un peu différemment notre équation différentielle pour respecter les contraintes de sympy :
\( \dfrac{dy}{dt} + 150y - 30 = 0 \)
toujours avec :
\( y(0) = 1/5 \).
Ouvrez maintenant une session Python dans votre IDE favori (pour moi Anaconda) et tapez les instructions suivantes dans la console Python :
from sympy import *
init_printing()
t = symbols("t")
y = Function("y")
edo = y(t).diff(t) + 150*y(t) - 30
dsolve(edo)

Si vous voulez un peu plus de détail sur sympy, allez sur cette page (pour les généralités sur sympy) ou sur celle-ci (pour la résolution des EDO avec sympy) de TangenteX.com.
Sympy nous indique la solution générale de notre EDO :

\( y(t) = \dfrac{C_1}{150} e^{-150t}+ \dfrac{1}{5} \)

avec C1 la constante d'intégration que nous allons déterminer avec la condition initiale.

La solution particulière pour y(0) = 1/5 s'obtient très facilement :

\( y(0) = \dfrac{C_1}{150} e^{-150*0}+ \dfrac{1}{5} \)  soit  \( \dfrac{1}{5} = \dfrac{C_1}{150}  + \dfrac{1}{5} \)

d'où C1 = 0, ce qui nous donne comme solution de notre EDO la fonction constante \( y(t) = \dfrac{1}{5}\).

Sa solution numérique par la méthode d'Euler

Rappel sur la méthode d'Euler

Nous avons déjà rencontré la méthode d'Euler sur TangenteX.com, par exemple sur cette page. Rappelons brièvement le principe de cette méthode.

Découpons le domaine d'intégration de l'équation différentielle, noté X par exemple, en plusieurs intervalles de largeur h. Le principe de la méthode d'Euler (ici la méthode explicite) est basé sur la définition approchée de la dérivée d'une fonction au point \( t_i \) de l'intervalle \( y'(t_i) = \dfrac{y(t_{i+1}) - y(t_i)}{h} \), relation qui nous permet de calculer explicitement par récurrence la valeur de \( y(t_{i+1}) \) en fonction de la valeur de \( y(t_i) \) et de h. Nous avons la relation de récurrence :

\( \left\{ \begin{array}{l} t_{i+1} = t_i + h \\ y(t_{i+1}) = y(t_i) + h*y'(t_i)  \end{array} \right. \)

Cette relation permet de déduire que les erreurs dues aux approximations du calcul de y à chaque intervalle se cumulent, ce qui rend cette méthode assez imprécise. Mais elle a le bon goût d'être très facile à coder...

Son implémentation en Python

Dans la page citée en référence, j'ai implémenté la méthode d'Euler en C. Voici ce que cela donne en Python :

def Euler(fonction, dom, y0, h):
   Y = [];y = y0   
   for t in dom:
        Y.append(y)
        y += h*fonction(y,t)
   return Y

A part les fioritures informatiques, une application pure et simple de la relation de récurrence vue plus haut.

L'influence de la valeur du pas sur la solution

Choisissons d'abord un pas h = 0.015. Cette valeur n'est pas prise au hasard, comme nous le verrons plus loin. Sur une même figure, nous tracerons la solution analytique (en noir) et la solution numérique, ici en vert la solution obtenue avec la méthode d'Euler.

Méthode d'Euler h = 0.015

Vous notez le comportement étrange de notre solution : elle diverge complétement ! Pourtant le problème est numériquement bien posé. En effet, si vous introduisez une petite perturbation dans la condition initiale y0, en modifiant la valeur de la variable epsilon, vous ne constaterez pas de grosse perturbation sur la solution, ce qui serait le signe d'un problème mal posé.

D'où vient donc le problème ? Utilisons la relation de récurrence de l'algorithme d'Euler pour exprimer la valeur de y à un point n quelconque de l'intervalle, en fonction de y0 et de n. On obtient :

\( y_n - \dfrac{1}{5} = (1 - 150h)^n (y_0 - \dfrac{1}{5}) \).

On constate immédiatement que la fonction diverge à l'infini si la valeur absolue de (1 - 150h) est supérieure à 1. Il faut donc qu'elle soit inférieure ou égale à 1, c'est à dire que \( h \le \dfrac{1}{75} \) pour assurer la convergence. Dans le calcul ci-dessus, le critère n'est pas rempli et notre solution diverge très rapidement.

Diminuons maintenant la valeur de h, en la fixant à 0.001 par exemple, pour respecter le critère sur h :

Méthode d'Euler h = 0.001

Notre solution se comporte beaucoup mieux et converge sans problème ! Avec la méthode d'Euler, pour parvenir à résoudre un problème raide, il faut fixer un pas h très petit, dans les limites fixées par la nature de l'équation. Cela peut augmenter considérablement le temps de calcul, ce qui n'est pas toujours acceptable, outre les autres défauts de la méthode d'Euler explicite.

C'est que l'on appelle un problème raide (ou "mal conditionné") : la convergence de la solution dépend de façon très importante de la valeur du pas de discrétisation de l'intervalle. Certaines méthodes numériques ne conviennent pas à la résolution d'EDO "raides", comme Euler explicite par exemple. Nous verrons que d'autres s'y prêtent mieux.

Sa solution par la méthode RK4

Rappel sur la méthode RK4

J'ai présenté la méthode de Runge-Kutta d'ordre 4 (dite RK4) et son implémentation en FORTRAN dans une autre page de TangenteX.com, je ne vais donc pas revenir dessus.

Ici, nous allons nous intéresser à l'influence de la méthode, qui prend en compte dans son algorithme les élements d'ordre plus élevé dans le développement limité de la dérivée, sur sa capacité à traiter des problèmes raides.

Son implémentation en Python

Voici une implémentation de la méthode RK4 en Python, que j'utiliserai dans le script ProblemeRaidePO.py :

def RK4(fonction, dom, y0, h):
    Y = []; y = y0    
    for t in dom:
        Y.append(y)   
        k1 = h*fonction(y,t)
        k2 = h*fonction(y + k1/2.0, t + h/2.0)
        k3 = h*fonction(y + k2/2.0, t + h/2.0)
        k4 = h*fonction(y + k3, t + h)
        y += (k1 + 2.0*k2 + 2.0*k3 + k4)/6.0
    return Y

Il s'agit ici aussi de l'implémentation immédiate de l'algoritme de la méthode RK4.

La solution

Utilisons maintenant la routine RK4 plutôt que la routine Euler pour résoudre notre équation différentielle, sur le même domaine et avec une valeur de h = 0.015, celle pour laquelle Euler donnait un résultat divergent. La courbe résultante est la suivante :

Méthode RK4 h = 0.015

Nous constatons que cette fois la solution numérique converge de manière très satisfaisante vers la solution analytique sur le domaine d'intégration. La méthode RK4 semble convenable pour résoudre des problèmes raides, au moins dans l'exemple ci-dessus.

En fait, elle peut être pris en défaut sur des domaines d'intégration plus larges, avec des équations un peu plus raides que notre exemple. La méthode est meilleure que celle d'Euler, mais ce n'est pas la panacée.

Comment traiter une EDO "raide"

Dans certains problèmes de modélisation, vous serez confrontés à une EDO "raide". Vous vous en apercevrez lorsque, utilisant la méthode d'Euler la plus simple, votre solution numérique divergera très rapidement, à moins d'utiliser des pas de calcul très petits et donc pénalisants. L'utilisation de la méthode RK4 permettra peut-être de traiter votre problème. Mais le plus raisonnable est d'utiliser la fonction odeint() de Python. Nous avons déjà utilisé cette fonction dans plusieurs pages.

Son appel est très simple. La solution est retournée dans la liste S :

S,infodict = odeint(EDO,y0,X, printmessg=True,full_output=True)

Attardons nous sur le second paramètre, infodict, retourné par la fonction parce que je lui ai demandé (le paramètre printmessg=True dans l'appel). Si j'imprime infodict, j'obtiens sur la console Python toute une série d'informations qui décrivent le déroulement de l'algorithme. Un élement particulier de cette liste m'intéresse, on le retrouve vers la fin :

'mused': array([1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
       2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
       2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2], dtype=int32)}

La documentation de odeint() m'indique que mused est un vecteur qui stocke pour chaque pas d'intégration la méthode d'intégration utilisée. 1 désigne la méthode d'Adams et 2 désigne la méthode BDF (pour Backward Differentiation Formulas ou méthode de différentiation rétrograde). Pour info, la méthode d'Adams est utilisée pour les problèmes non raides et la méthode BDF pour les problèmes raides.

Vous remarquez que odeint() débute l'intégration en utilisant la méthode d'Adams, mais dès qu'elle détecte la "raideur" du problème, elle passe à une autre méthode, la méthode BDF, adaptée aux problèmes raides. Imaginez que vous ayez à programmer un tel mécanisme ...

Sauf s'il s'agit d'un programme pédagogique, je vous conseille donc fortement d'utiliser odeint() ou même ode(), cette dernière vous permettant une grande souplesse dans l'usage de différentes méthodes d'intégration. Le conseil est aussi valable si vous programmez avec SciLab ou Maple, qui possèdent les mêmes fonctions d'intégration d'une EDO, d'ailleurs basées sur la même librairie FORTRAN ODEPACK qui existe depuis quelques décennies...

La courbe résultante est sans surprise, la solution est satisfaisante :

odeint() h = 0.015

Systèmes différentiels et problèmes raides

Les problèmes raides peuvent également concerner les systèmes d'équations différentielles. On en trouve souvent dans les modèles de cinétique chimique, de dynamique des populations ou de calcul de structures. Les problèmes raides apparaissent lorsque la dynamique du système n'est pas homogène, c'est à dire que l'on observe une évolution lente du système dans une partie du domaine d'intégration et une évolution très rapide dans une partie généralement petite du domaine. On perçoit intuitivement que le pas d'intégration ne peut pas être identique dans les deux cas : il peut être relativement grand dans la première partie (évolution lente) et doit être très petit dans la partie "rapide" pour garantir la même qualité d'intégration. Pensez à la courbe classique de titrage en chimie et à sa pente très forte autour du point d'équivalence.

Sans entrer dans le détail de la problématique, voyons ce que cela donne sur un exemple simple et hyperclassique.

Un exemple

Le système

Considérons le système différentiel linéaire à coefficients constants, tiré d'un problème de cinétique chimique, suivant :

\( \left\{ \begin{array}{l} \dfrac{dx}{dt} = -20x - 19y \\ \dfrac{dy}{dt} = -19x - 20y \\ x(0) = 2 \\ y(0) = 0  \end{array} \right. \)

que je vais intégrer sur le domaine [0,1].

Sa solution analytique
Jacobienne, valeurs propres et vecteurs propres

Le premier réflexe, dans le cas d'un système différentiel, c'est d'en vérifier la stabilité. La théorie n'est pas immédiate mais la technique est simple : il suffit de calculer la matrice de Jacobi (ou jacobienne) du système et de vérifier ses valeurs propres.

Par définition (voir par exemple ici), la matrice jacobienne du système est égale à :

\( J = \left( \begin{array}{rcl} \dfrac{-20x - 19y}{\partial x}&  \dfrac{-20x - 19y}{\partial y} \\ \dfrac{-19x - 20y}{\partial x}&  \dfrac{-19x - 20y}{\partial y} \end{array} \right) \) soit la matrice \( J = \left( \begin{array}{rcl} -20&  -19 \\ -19&  -20 \end{array} \right) \)

Je peux donc écrire notre système sous forme matricielle : \( \dfrac{dF}{dt} = JF \) où J est la matrice \( \left( \begin{array}{rcl} -20&  -19 \\ -19&  -20 \end{array} \right) \)

Python nous permet de calculer rapidement les valeurs propres et les vecteurs propres de cette jacobienne. Sur la console Python, tapez les instructions suivantes, d'abord la définition de la jacobienne :

J = [[-20,-19],[-19,-20]]

puis, après avoir importé le package nécessaire, le calcul des valeurs propres et vecteurs propres :

from scipy.linalg import eig

eig(J)

ce qui nous donne un premier tableau contenant les valeurs propres et un second contenant les vecteurs propres sous forme de matrice (les vecteurs propres sont les colonnes de la matrice) :

(array([ -1.+0.j, -39.+0.j]), array([[ 0.70710678,  0.70710678], [-0.70710678,  0.70710678]]))

Nous constatons que les valeurs propres sont réelles et négatives : \( \lambda_1 = -1 \) et \( \lambda_2 = -39 \) . Les vecteurs propres associés sont V1 = [ 0.70710678,  -0.70710678] et V2 = [0.70710678,  0.70710678] ou encore \( V1 = \left( \begin{array}{rcl} \dfrac{\sqrt 2}{2} \\ -\dfrac{\sqrt 2}{2} \end{array} \right) \) et \( V2 = \left( \begin{array}{rcl} \dfrac{\sqrt 2}{2} \\ \dfrac{\sqrt 2}{2} \end{array} \right) \).

D'après la théorie des systèmes dynamiques, notre système est stable car ses deux valeurs propres sont réelles et ont toutes les deux des parties réelles négatives. Oui mais ... On constate que ces deux valeurs sont significativement différentes. C'est l'indice d'un problème raide. Lorsque les valeurs propres sont proches en valeur absolue, pas de problème. Mais dès qu'elles s'éloignent, il faut suspecter un problème raide. Et plus elles s'éloignent et plus le problème est raide...

Donc, face à un système différentiel : déterminer la matrice jacobienne puis calculer ses valeurs propres. Les comparer. Si elles sont trop différentes, on a affaire à un problème raide.

La solution générale et particulière du système

La solution générale du système s'exprime par une combinaison linéaire d'exponentielles qui dépendent des valeurs propres et des vecteurs propres de la jacobienne du système. Dans notre cas, la solution générale s'écrit :

\( X\left[ \begin{array}{rcl} x(t) \\ y(t) \end{array} \right] = a e^{\lambda_1 t} V1 +  b e^{\lambda_2 t} V2 \)

soit encore :

\( X\left[ \begin{array}{rcl} x(t) \\ y(t) \end{array} \right] = a e^{\lambda_1 t} \left( \begin{array}{rcl} \dfrac{\sqrt 2}{2} \\ -\dfrac{\sqrt 2}{2} \end{array} \right) +  b e^{\lambda_2 t} \left( \begin{array}{rcl} \dfrac{\sqrt 2}{2} \\ \dfrac{\sqrt 2}{2} \end{array} \right) \)

Les paramètres a et b dépendent des conditions initiales, soit x(0) = 2 et y(0) = 0. Il suffit de faire la substitution dans le système et l'on obtient un système linéaire simple qui nous permet de calculer a et b :

\( \left\{ \begin{array}{l} x(0) = a \dfrac{\sqrt 2}{2} + b \dfrac{\sqrt 2}{2} = 2 \\ y(0) = -a \dfrac{\sqrt 2}{2} + b \dfrac{\sqrt 2}{2} = 0   \end{array} \right. \)

La solution particulière de notre système devient donc :

\( \left\{ \begin{array}{l} x(t) = e^{-t} + e^{-39t}\\ y(t) = -e^{-t} + e^{-39t}   \end{array} \right. \)

système que nous utiliserons pour calculer la solution analytique avec la fonction Python suivante :

def SolAna(dom,vp1,vp2):
    xi = [];yi = []
    for t in dom:
        xi.append(
exp(vp1*t)  + exp(vp2*t))
        yi.append(-exp(vp1*t) + exp(vp2*t))
    return (xi,yi)

Comparaison des différentes méthodes d'intégration

Intégration avec un pas h = 0.05

Lançons le script ProblemeRaideSys.py, en fixant la valeur du pas dx = 0.05. La figure ci-dessous illustre les solutions x(t) et y(t) pour différentes méthodes d'intégration : Euler, RK4 et odeint(), avec la solution analytique en noir :

On constate plusieurs choses :

Intégration avec un pas h = 0.01

Lançons le script ProblemeRaideSys.py, en fixant la valeur du pas dx = 0.01. Nous obtenons :

Nous constatons que :

Les scripts Python

Les scripts Python étudiés dans cette page sont disponibles dans le package ProblemeRaide.zip :

Contenu et design par Dominique Lefebvre - www.tangenteX.com septembre 2017   Licence Creative Commons   Contact :

Cette œuvre est mise à disposition selon les termes de la Licence Creative Commons Attribution - Pas d’Utilisation Commerciale - Pas de Modification 3.0 France.