Pour décrire un processus d’évolution, ou le profil d’une ligne d’eau, par exemple, on est souvent amené à résoudre une équation différentielle ordinaire (EDO) du premier ordre. Cette équation écrit comment varie une fonction, en un point donné (un instant ou un point de l’espace), connaissant la valeur de cette fonction mathématique, le problème à résoudre s’écrit :
$$ \left\{ \begin{array}{rcr} \frac{dy}{dt} & = & f(y,t) \\ y(t=t_0) & = & y_0 \\ \end{array} \right. $$
où $\frac{dy}{dt}$ désigne la dérivée par rapport à t de la fonction $y$ (qui dépend de la variable $t$) ; la variable $y_0$ est appelée la condition à la limite ; elle conditionne la solution finale de l’équation.
Comme souvent on ne connait pas de solution analytique de ce problème, on va utiliser des méthodes approchées pour estimer la solution. On fait donc une discrétisation de la variable $t$. On note ainsi $\Delta t$ le pas de discrétisation, et on résout le problème aux points $t_0$, $t_1=t_0+\Delta t$, $t_2=t_0+2\Delta t$, ..., $t_n=t_0+n\Delta t$ où $n$ est un entier.
Méthode d’Euler explicite
Cette méthode est la plus intuitive ; elle consiste à considérer que, d’un point $t_i$ au point $t_{i+1}$, la fonction évolue linéairement, avec une trajectoire qui est celle qu’on peut calculer au point $t_i$.
Le problème se résout donc de la façon suivante :
on connait la fonction $f$, un point $t_i$ où on connait $y_i$
on peut donc calculer $y’_i=f(y,t)$
on estime alors la valeur de $y$ au point $t_{i+1}=t_i+\Delta t$ :
$$ y_{i+1}\simeq y_i + y’_i \Delta t $$
on peut alors itérer (résoudre pas à pas) pour passer au point suivant. Le problème est initialisé en partant de $t_0$ où on connait $y_0$ (condition à la limite).
On sent bien que ce schéma pourra donner de bons résultats uniquement si $\Delta t$ n’est pas trop grand. Des valeurs de $\Delta t$ trop grandes peuvent donner des résultats complètement faux, conduisant à des interprétations physiques erronées. Son intérêt est toutefois sa simplicité, et il s’implémente facilement sur un tableau.
Exemple d’application : processus exponentiel
Considérons le problème (simple) suivant :
$$ \left\{ \begin{array}{rcr}\label{eq:exp} \frac{dy}{dt} & = & -a y \\ y(t=t_0) & = & y_0 \\ \end{array} \right. $$
On a donc ici $f(y,t)=-ay$. La solution analytique se résout facilement, donnant $y(t)=y_0 \exp\left(-a(t-t_0)\right)$.
On peut résoudre le problème par la méthode d’Euler :
on choisit $\Delta t$ (par exemple, $\Delta t=1$)
calculer $ y_1=y_0 - a y_0 \Delta t$
calculer $ y_2=y_1 - a y_1 \Delta t$ etc.
On constate que la résolution n’est pas très précise ; ceci est lié au pas de calcul trop grand compte tenu de la méthode choisie et de l’équation à résoudre.
Schéma de Runge-Kutta d’ordre 4
Le schéma de Runge-Kutta d’ordre 4 est basée sur une approximation de la dérivée à un ordre supérieur (ordre 4). Le principe de la discrétisation est le même, mais on va faire quelques calculs supplémentaires pour approcher la dérivée :
on connait la fonction $f$, un point $t_i$ où on connait $y_i$
on peut donc calculer $f_1=y’_i=f(y_i,t_i)$ (cf. méthode d’Euler=pente au point ($t_i,y_i$))
on calcule $f_2=f(y_i+\frac12 \Delta t f_1,t_i+\frac12 \Delta t)$ (valeur estimée au milieu de l’intervalle, avec la pente prise en $t_i$)
on calcule $f_3=f(y_i+\frac12 \Delta t f_2,t_i+\frac12 \Delta t)$ (valeur estimée au milieu de l’intervalle $t_{i+1/2}$, avec la pente prise en $t_{i+1/2}$)
on calcule $f_4=f(y_i+ \Delta t f_3,t_i+ \Delta t)$ (valeur estimée en $t_{i+1}$, avec la pente prise en $t_{i+1/2}$)
on a alors
$$ y_{i+1}\simeq y_1 + \frac16 \Delta t(f_1+2f_2+2f_3+f_4) $$
on peut alors itérer (résoudre pas à pas) pour passer au point suivant. Le problème est initialisé en partant de $t_0$ où on connait $y_0$ (condition à la limite).
On voit clairement que la méthode est beaucoup plus précise. Même avec un pas de calcul beaucoup plus élevé, on approche correctement la solution, alors que la méthode d’Euler donne des résultats très éloignés de la solution exacte. On remarque aussi que, entre les points de discrétisation, on a approché la solution par des segments de droite (bien qu’on puisse avoir une interpolation plus évoluée si on a besoin de connaître des valeurs intermédiaires).
Le programme peut s’écrire de la façon suivante, avec le logiciel de calcul scientifique Scilab :