Cette semaine était chargée en maths ! J’ai continué à lire comme d’habitude Pixar Physics Simulation Animation Course, en particulier l’article sur des méthodes numériques implicites.

Il s’agit de résoudre des EDOs (équations différentielles ordinaires), considérées comme ‘raide’. On va voir pourquoi nous qualifions ces ODEs de ‘raide’ et comment nous devrions aborder la résolution de ces types d’ODE en utilisant des méthodes numériques implicites telles que la méthode Euler implicite.

⚠️ La ressource de Pixar nous conseille d’éviter une EDO raide en reformulant notre problème dynamique. Nous nous concentrerons sur les cas où cela n’est pas possible.

EDOs raides

Après avoir cherché une définition de ‘équation différentielle ordinaire raide’ sur Internet, j’ai fini par utiliser celle de Wikipédia qui, certes pour moi, n’est pas très suffisante, mais je vais quand même vous la donne :

Une équation différentielle raide est une équation différentielle dont la sensibilité aux paramètres va rendre difficile la résolution par des méthodes numériques explicites.

Bref : des EDOs raides sont difficiles à résoudre avec des méthodes numériques explicites (nous avons parlé des méthodes numériques explicites il y a quelques semaines).

Moi, je trouve que Pixar Physics Simulation Animation Course nous donne un exemple amusant sur lequel on va travailler pour mieux visualiser la définition Wikipédia. Imaginons que l’on veut simuler la dynamique d’une particule dont le changement de position se décrit par l’équation différentielle suivante :

où k est une grande constante positive. Si k est suffisamment grande, la particule ne bougera jamais trop loin de $y(t) = 0$ car $-ky(t)$ fera la revenir $y(t)$ à l’origine. On voit également que l’abscisse de la particule est plus libre de se déplacer, comme décrit par l’absence de cette grande constante positive k.

Dans cet exemple, cette EDO serait donc ‘raide’ si la variable k est si grande que la particule reste proche de la droite $y = 0$. Comme nous le verrons dans la partie suivante, les méthodes explicites auront du mal à simuler ce comportement attendu.

Le problème avec des méthodes explicites

Rappelons que l’équation de la méthode d’Euler explicite se décrit comme suit :

et dans cet exemple, $X_0 = \begin{pmatrix} x_0 \\ y_0 \end{pmatrix}$ et $\dot X(t) = \frac{d}{dt} \begin{pmatrix} x(t_0) \\ y(t_0) \end{pmatrix} = \begin{pmatrix} -x_0 \\ -ky_0 \end{pmatrix}$. Alors, après la substitution :

Ce qui produit la suivante :

Nous avons appris dans la discussion Partie 1 des méthodes numériques que pour trouver une solution, on doit prendre des petits pas dans le temps afin de s’approcher de la solution, car des pas trop grands pourront rendre le système instable. En réarrangeant les termes dans l’équation, on s’aperçoit rapidement que ces conditions doivent être remplies pour que la position de la particule décrite dans notre exemple se converge vers (0,0) :

🚨 Si nous prenons k égal à 1000 : le pas maximum permis est calculé comme suit : $1 - 1000h = -1 => h = 0,002$. Si on dépasse cette limite, l’ordonnée de la particule ne pourra jamais converger vers 0 !

Résolvez pour y en utilisant la méthode d’Euler explicite les paramètres h=0.0015 et k=1000.0 dans un solveur Euler explicite que j’ai écrit, nous pouvons visualiser l’ordonnée de la particule lorsqu’elle oscille et s’approche autour de zéro :

Voici comment l’ordonnée se converge, ou plutôt se diverge, qu’on utilise une valeur h > 0.002 (h = 0.003 et k = 1000.0):

OK, nous pourrions donc utiliser un très petit pas comme h = 0.0015 n’est-ce pas ? Peut-être… cependant, on se rendrait vite compte que la valeur d’abscisse bouge trop lentement, ce qui rendra notre simulation trop lente !

Pouvons-nous faire mieux ? Pourrions-nous utiliser une méthode plus stable pour se converger plus rapidement ? Oui - avec la méthode d’Euler implicite !

Euler implicite

La méthode d’Euler implicite est la plus simple des méthodes implicites à comprendre et c’est pourquoi je vais me concentrer dessus.

Voici l’équation d’Euler implicite :

Si vous aviez été quelqu’un d’observateur, vous auriez vu que la différence entre l’équation d’Euler implicite et celle d’Euler explicite est que le changement d’X se définit comme $h\dot X*{new}$ dans l’équation d’Euler implicite. Il faut donc résoudre ce changement en calculant $\dot X*{new}$… pour calculer enfin $X\_{new}$.

En fait, c’est comme marcher à reculons.

Autrement dit, c’est comme si on était à la position $X\_{new}$ et nous aimerions reculer pour arriver à la position $X_0$. Il faut dire qu’en anglais, Euler implicite s’appelle aussi ‘Backwards Euler’ pour cette raison.

C’est trop de résoudre $X*{new}$, d’autant plus que nous devons aussi résoudre $\dot X*{new}$… qui nous est aussi inconnue. Voyons s’il est possible de réécrire cela dans une équation résoluble.

On commence par fournir l’équation de la méthode d’Euler implicite :

où $\Delta X = X\_{new} - X_0$. Nous pouvons l’écrire comme suit :

En utilisant la dernière équation, on peut réécrire notre équation pour $\Delta X$ :

Nous pouvons tout multiplier par une matrice identité pour isoler notre variable $\Delta X$.

Après avoir isolé $\Delta X$, il est dorénavant possible de résoudre $X\_{new} = X_0 + \Delta X$. On va voir comment le faire avec du code en C dans la partie suivante.

Enfin, une solution et du code

Voici le code auquel je vais faire référence dans cette section.

Pour rappel, on veut faire converger cette EDO raide :

Et nous savons que $\frac{\partial}{\partial X}f(X(t))$ est la matrice jacobienne :

Revenons à l’équation $\Delta X = (\frac{1}{h}I - f'(X_0))^{-1} f(X_0)$. On calcule $\Delta X$ comme suit :

Ce qui est intéressant, c’est qu’on peut faire un pas avec n’importe quelle amplitude ! Si l’on calcule la limite de $\Delta X$ lorsque h s’approche d’infini, on verrait que :

Donc $X*{new} = X_0 + -(X_0) = \begin{pmatrix} 0 \\ 0 \end{pmatrix}$ lorsque _h* s’approche de l’infini !

Maintenant, nous pouvons utiliser le programme que j’ai écrit pour voir cet exemple. Après avoir navigué jusqu’au dossier où vous avez téléchargé le code, exécute la commande suivante :

où h=0.2, k=1000.0, x₀=100.0 et y₀=100.0. Vous devriez voir deux fichiers implicit.txt et explicit.txt qui viennent s’ajouter au dossier. Dans implicit.txt, on voit que y se converge vers 0 !

Mais dans explicit.txt, c’est une histoire complètement différente :

La suite

On verra comment résoudre des équations différentielles raides du second ordre ! À la prochaine !

Des ressources (en français et anglais)

• Wikipédia - Équation différentielle raide
• Large Steps in Cloth Simulation - CMU
• Implicit vs. Explicit Numerical Methods
• Developing a Matrix Library in C - UNC
• Approximating the Solution to the Cauchy Problem in C++
• Lecture 26 Implicit Euler’s Method - YouTube