Résolution numérique d’équations différentielles

Équations différentielle scalaires d’ordre 1

On étudie ici des équations différentielles d’ordre 1 qui s’écrivent sous la forme :

$x'=f(x,t)$

où $f$ est une fonction de deux variables et où l’on cherche la fonction x de la variable $t$ définie et dérivable sur un intervalle $I$ de $\mathbb {R}$ vérifiant : $\left\lbrace {\begin{array}{l}\forall t\in I,(x(t),t)\in I\\\forall t\in I,x'(t)=f(x(t),t)\\\end{array}}\right.$

Afin d'obtenir une solution unique, nous allons adjoindre à cette équation différentielle une condition initiale sous la forme d’un couple $(x0,t0)\in I$ et chercher à résoudre le problème de Cauchy suivant : $\left\lbrace {\begin{array}{l}x'=f(x,t)\\x(t_{0})=x_{0}\\\end{array}}\right.$

Sous certaines conditions sur f que nous ne détaillerons pas, ce problème admet une unique solution, que nous allons chercher à déterminer numériquement.

Principe de la méthode d'Euler

La méthode d'Euler repose sur l'utilisation de l'approximation affine de la fonction. La dérivée $x'(t)={\frac {dx(t)}{dt}}$ est approximée par ${\frac {x_{n+1}-x_{n}}{t_{n+1}-t_{n}}}={\frac {x_{n+1}-x_{n}}{h}}$ où $h$ est le pas de calcul (suffisamment petit).

On peut donc écrire la relation de récurrence : $hx'(t_{n})=x_{n+1}-x_{n}\Rightarrow x_{n+1}=x_{n}+hx'(t_{n})\Rightarrow x_{n+1}=x_{n}+h.f(x_{n},t_{n})$

L’algorithme de résolution utilisera donc les deux relations de récurrence suivantes : $\left\lbrace {\begin{array}{l}t_{n+1}=t_{n}+h\\x_{n+1}=x_{n}+hf(x_{n},y_{n})\\\end{array}}\right.$

Exemple : chute libre

Considérons un objet de masse $M$ en chute libre verticale avec frottement fluide de coefficient $\lambda$ .

Sa trajectoire est repérée dans un référentiel galiléen avec l'axe ${\vec {z}}$ vers le bas.

Sa vitesse $v$ projetée sur l'axe ${\vec {z}}$ respecte l'équation de la dynamique suivante:

$M\mathrm {.} v{\text{'}}\left(t\right){=}M\mathrm {.} g-\mathrm {\lambda } v\left(t\right)$ (1)

Avec la condition initiale :

$v(0)=v_{0}$ (2)

L'équation (1) devient : $v{\text{'}}\left(t\right)=g-{\frac {\mathrm {\lambda } }{M}}v\left(t\right)$

Remarque : La solution exacte à ce problème est trivialement connue : $v\left(t\right)={\frac {m\mathrm {.} g}{\mathrm {\lambda } }}\left(1-{e}^{{\frac {-\mathrm {\lambda } }{m}}t}\right)$

Pour la solution approchée, on obtient l’algorithme suivant :

$\left\lbrace {\begin{array}{l}{t}_{n+1}={t}_{n}+h\\{v}_{n+1}={v}_{n}+h(g-{\frac {\lambda }{M}}{v}_{n})\\\end{array}}\right.$

Programmation python

On placera les valeurs ${t}_{n}$dans un tableau nommé t

De même on placera les valeurs ${v}_{n}$dans un tableau nommé v_euler

La boucle de calcul des approximations successives s'écrit :

1 for n in range (0, Nb_iterations):
2     t.append(t[n]+h)     # ajout d'une valeur dans le tableau recevant les valeurs de temps
3     v.append(v_euler[n]+h*(-g-v_euler[n]*lambda/m))    # ajout d'une valeur dans le tableau recevant les vitesses