Des méthodes numériques (Implicite, 2nd ordre EDO) - Partie 4

Nous apprenons comment et quand utiliser des méthodes numériques implicites pour résoudre les équations différentielles du second ordre.

Keywords: C, des méthodes numériques, maths, les simulations physiques

By Carmen Cincotti

April 11, 2022

Passons en revue : l’espace des phases
Équations différentielles ordinaires (EDOs) du second ordre
La suite
Des ressources (en français et anglais)

Je continue à bosser sur Pixar Physics Simulation Animation Course, en particulier l’article sur des méthodes numériques implicites. La semaine dernière, nous nous sommes introduits aux méthodes implicites du premier ordre. Cette semaine, nous travaillerons plus dur en étudiant des méthodes implicites du second ordre.

Cette leçon ouvrira la porte à la compréhension et à la création de simulations de tissu. Bien entendu, les simulations de tissu (tissés et tricotés, si vous voulez) consistent à nous servir des systèmes masse-ressorts. Contrairement à la simulation du système de particules que l’on a fait il y a quelques semaines, on utilisera un solveur d’Euler implicite du second ordre pour créer une telle simulation !

Passons en revue : l’espace des phases

Une simulation de tissu est juste comme une simulation d’un système de particules. Rappelons que des particules se représent en espace dans une configuration qui s’appelle espace des phases.

Chaque particule se représente dans l’espace par deux vecteurs de coordonnées généralisées (x, v) où x est un vecteur de position et v est un vecteur de vitesse. Le vecteur de l’espace des phases $p_n$ se construit des vecteurs x et y l’un après l’autre, ainsi il aura une taille finale de 2*n éléments (où n est la quantité de dimensions représentées dans notre système).

p_n = \begin{pmatrix} x_x \\ x_y \\ x_z \\ v_x \\ v_y \\ v_z \\ \end{pmatrix}

⚠️ Si vous êtes déjà perdu, il faut voir la discussion sur les particules où j’ai introduit ce sujet. Il faut que vous compreniez ce sujet avant de continuer !

Équations différentielles ordinaires (EDOs) du second ordre

La ressource de Pixar nous dit que :

La plupart des problèmes dynamiques sont écrits dans les termes d’une équation différentielle du second ordre :
$\bf \ddot x \textnormal{(t)} = f(x\textnormal{(t)}, \dot x\textnormal{(t)})$

où $\ddot x(t)$ est la dérivée de second ordre de la position x(t). Nous savons déjà que cette valeur est l’accélération. $\dot x(t)$ est la dérivée de la position… que l’on connait sous le nom de vitesse.

La résolution d’une équation différentielle du second ordre n’est pas idéale parce qu’elle est coûteuse et compliquée. Nous verrons que nous pouvons éviter ce calcul en approximant le terme du second ordre à l’aide de la série de Taylor dans un processus nommé linéarisation. En introduisant et en définissant la variable $v = \dot x$ , nous pouvons réécrire ces équations au premier ordre :

\frac{d}{dt}\begin{pmatrix} \bf x\textnormal{(t)} \\ \bf v\textnormal{(t)} \end{pmatrix} = \begin{pmatrix} \bf v\textnormal{(t)} \\ \bf f(x\textnormal{(t)}, \bf v\textnormal{(t)}) \end{pmatrix}

Simplifions ce système en rappelant que l’étape d’Euler implicite se définit comme suit :

x_{new} = x_0 + hf(x_{new}) \\ \because \Delta x = x_{new} - x_0 \\ \therefore \Delta x = h f(x_{new})

et que :

\Delta x = x(t_0 + h) - x(t_0) = \frac{d}{dt} x(t)

…on peut introduire des variables $x_0 = x(t_0)$ et $v_0 = v(t_0)$ pour enfin avoir une équation pour résoudre une étape d’Euler implicite :

\begin{pmatrix} \bf \Delta x \\ \bf \Delta v \end{pmatrix} = h \begin{pmatrix} \bf v_0 + \Delta v \\ \bf f(x_0 + \Delta x, v_0 + \Delta v) \end{pmatrix}

Linéarisation

Comme je l’ai dit, nous pouvons approximer $f(x_0 + \Delta x, v_0 + \Delta v)$ grâce aux séries de Taylor :

f(x_n + \Delta x, v_n + \Delta v) \approx f_n + \frac{\partial f}{\partial x} \Delta x + \frac{\partial f}{\partial v} \Delta v

💡 On peut faire une approximation d’une fonction réelle avec deux variables $f(x, y)$ en utilisant la série de Taylor comme suit :

f(x + \Delta x, y + \Delta y) \approx f(x,y) + \frac{\partial f(x,y)}{\partial x} \Delta x + \frac{\partial f(x,y)}{\partial y} \Delta y + ...

Vois cette ressource pour en savoir plus sur la série de Taylor et cette dérivation exacte.

Alors, l’équation devient :

\begin{pmatrix} \bf \Delta x \\ \bf \Delta v \end{pmatrix} = h \begin{pmatrix} \bf v_0 + \Delta v \\ \bf f_n + \frac{\partial f}{\partial x} \Delta x + \frac{\partial f}{\partial v} \Delta v \end{pmatrix}

Un aparté : dimensions

Il faut dire que les dimensions de $\Delta x$ et $\Delta v$ sont décidées par nous. Dans cet exemple, nous pouvons supposer que notre particule est définie par une position en 3D $\bf x \in \R^3$ et une vitesse $\bf v \in \R^3$ .

Cette distinction nous permet de mieux comprendre les dimensions de ce calcul. Supposons également dans ce cas que $f$ est en 3D. Une ressource que je vous recommande pour créer des simulations de tissus nous dit que $f$ est la force dans le système, mais je n’en suis pas si sûr à ce stade parce que je ne l’ai pas encore étudié. En tout cas, nous pouvons supposer que $\frac{\partial f}{\partial x} \in \R^{3x3}$

💡 Grâce au calcul vectoriel, on voit comment calculer la dérivation d’un vecteur en 3D :

\begin{bmatrix} f_x(x,y,z) \\ f_y(x,y,z) \\ f_z(x,y,z) \end{bmatrix} = \begin{bmatrix} \dfrac{\partial f_x}{\partial x} & \dfrac{\partial f_x}{\partial y} & \dfrac{\partial f_x}{\partial z} \\[1em] \dfrac{\partial f_y}{\partial x} & \dfrac{\partial f_y}{\partial y} & \dfrac{\partial f_y}{\partial z} \\[1em] \dfrac{\partial f_z}{\partial x} & \dfrac{\partial f_z}{\partial y} & \dfrac{\partial f_z}{\partial z} \\[1em] \end{bmatrix}

La dérivation est appelée la matrice jacobienne. On y reviendra avec plus de détails dans une partie ultérieure

En remplaçant $\Delta x$ dans l’équation de $\Delta v$ , nous arrivons au suivant :

\bf \Delta v = h (\bf f_0 + \frac{\partial f}{\partial x} \textnormal{h}(v_0 + \Delta v) + \frac{\partial f}{\partial v} \Delta v)

La forme linéaire

Nous sommes enfin arrivés au cœur d’une méthode implicite. On voit facilement que $\Delta v$ se met en chaque côté du signe égal. Il ne suffit pas de résoudre la valeur de la variable $\Delta v$ en la plaçant d’un seul côté. Il faut donc résoudre un système linéaire en utilisant la form linéaire $Ax = b$ .

💡 Un système linéaire peut être résolu en utilisant la forme :

Ax = b

où A est la matrice des coefficients, x est le vecteur solution et b est le vecteur constant. Prenons cet exemple avec un système linéaire comme suit :

3x + 2 = 2 \\ x - y = 4 \\ 5y + z = -1

En utilisant la forme $Ax = b$ , nous arrivons à :

\begin{bmatrix} 3 & 2 & 0 \\ 1 & -1 & 0 \\ 0 & 5 & 1 \\ \end{bmatrix} \begin{bmatrix} x \\ y \\ z \\ \end{bmatrix} = \begin{bmatrix} 2 \\ 4 \\ -1 \\ \end{bmatrix}

Après avoir arrangé nos termes dans la forme $Ax = b$ on arrive à :

\bf(I - \textnormal{h}\frac{\partial f}{\partial v} - \textnormal{h}^2\frac{\partial f}{\partial x})\Delta v = \textnormal{h}(f_0 + \textnormal{h}\frac{\partial f}{\partial x}v_0)

où I est la matrice identité ( $\bf I \in \R^{3x3}$ ).

L’étape suivante consiste bien sûr à résoudre $\Delta v$ . Je le ferai dans la prochaine série - une simulation de tissu !

La suite

Je continue à étudier la physique en créant une simulation de tissu !

Carmen's Graphics Blog