Code
import numpy as np
import matplotlib.pyplot as plt1 Équations du mouvement, intégration numérique
Dans cette leçon, nous verrons comment résoudre numériquement les équations du mouvement issues du principe fondamental de la dynamique.
Note
Ce chapitre contient des exercices simples corrigés. Le code de correction peut être consulté en cliquant sur le texte indiquant “▸ Code”.
1.1 Seconde loi de Newton
On considère un corpuscule, ou une particule, de masse \(m\), position \(\vec{r}\) et vitesse \(\vec{v}\) définis au temps \(t\). Sous l’action d’une force \(\vec{f}\), la quantité de mouvement de la particule \(\vec{p} = m\vec{v}\) change selon la seconde loi de Newton : \[ \frac{\mathrm{d} \vec{p}}{\mathrm{d} t} = \vec{f}, \] ou, si l’on suppose la masse constante : \[ m\frac{\mathrm{d} \vec{v}}{\mathrm{d} t} = \vec{f}.\]
La seconde loi de Newton se généralise très bien à un système à \(N\) particules. On notera \(\{\vec{r}_i\}_{i=1}^N\) l’ensemble des positions, \(\{\vec{v}_i\}_{i=1}^N\) l’ensemble des vitesses, \(\{m_i\}_{i=1}^N\) l’ensemble des masses et \(\{\vec{f}_i\}_{i=1}^N\) l’ensemble des forces totales agissant sur chacune des particules. Nous verrons que le calcul de \(\vec{f}_i\) pour chaque particule est central dans les méthodes de simulations discrètes, puisque \(\vec{f}_i\) contrôle le mouvement de la particule \(i\). On peut alors écrire la seconde loi de Newton pour un système de \(N\) particules :
\[ m_i\frac{\mathrm{d}\vec{v}_i}{\mathrm{d}t} = \vec{f}_i. \tag{1.1}\]
1.2 Discrétisation
Le problème à trois corps est un problème de physique classique où trois particules interagissent par l’attraction gravitationnelle. Dans le cas général, ce problème n’est pas soluble analytiquement1, mais on peut le résoudre de manière approchée par des techniques numériques.
Les techniques de résolution numérique sont essentiellement basées sur un principe : celui de transformer une équation différentielle en un ensemble d’équations algébriques (qui ne font pas intervenir la dérivée). C’est le procédé de discrétisation. Cependant, ce procédé a un coût, celui de n’obtenir qu’une solution approchée du problème original. L’erreur commise, appelée erreur de discrétisation, est souvent inévitable, mais on peut la contrôler dans beaucoup de cas.
On va discuter dans cette leçon de différentes façons d’approximer la solution d’une équation différentielle, et d’une méthode acceptable pour la solution approchée de la seconde loi de Newton.
Partition du temps
La première étape du processus de discrétisation est de sélectionner des valeurs discrètes de la variable de notre équation différentielle, ici le temps. Dans le problème original \(t \in [0, +\infty)\). On défini donc un ensemble de \(M\) temps discrets : \[t_i = \Delta t\cdot i,\quad i \in [0, M-1],\] ainsi les temps \(t_i\) sont uniformément espacés de \(\Delta t\), que l’on appellera pas de temps. On verra que cette quantité contrôle bon nombre de propriétés numériques de notre système (les propriétés numériques sont uniquement dépendantes des choix de discrétisation, et sont à distinguer des propriétés physiques du système, qui dépendent des équations gouvernant le système).
1.3 Intégration en temps – premier ordre
Avant de discuter de l’équation de Newton, intéressons-nous à une équation différentielle plus simple :
\[ \dot{y} + y = 0,\quad y(0) = 1, \]
dont on connait la solution, \(y(t) = e^{-t}\). Posons l’équation pour les temps discrets \(t_i\). On introduit la notation \(y_i = y(t_i)\):
\[ \dot{y}_i + y_i = 0,\quad i\in[0, M-1]. \]
Schéma d’Euler explicite
Pour trouver les valeurs inconnues \(y_i\), il nous faut la dérivée de \(y\), que l’on ne connait pas. En revanche, on peut l’approximer par différence finie :
\[ \dot{y}_i \approx \frac{y(t_{i+1}) - y(t_i)}{t_{i+1} - t_i} = \frac{y_{i+1} - y_i}{\Delta t} \]
De toute évidence l’approximation devient exacte quand \(\Delta t \rightarrow 0\). On peut remplacer notre approximation dans l’équation différentielle :
\[ \frac{y_{i+1} - y_i}{\Delta t} + y_i = 0 \quad \Leftrightarrow \quad y_{i+1} = (1 - \Delta t)y_i,\quad i\in[0, M-1] \]
On a, par approximation de la dérivée en différences finies, transformé une équation différentielle en \(M\) équations algébriques, qui définissent en fait une relation de récurrence pour la suite \(\{y_i\}_{i=0}^{M-1}\). Puisque la relation de récurrence définit une suite géométrique, on peut écrire \(y_i\) directement :
\[ y_i = (1 - \Delta t)^i y_0,\] on remarque alors une propriété particulière : si \(\Delta t < 1\) la limite de la suite est nulle et la suite est monotone, comme la solution analytique, mais si \(2 > \Delta t > 1\) la suite oscille (si \(\Delta = 1\) la suite est nulle) et si \(\Delta t > 2\) la suite diverge. C’est l’exemple parfait d’une propriété numérique, ici que \(\Delta t = 1\) est un pas de temps critique qui conditionne la bonne marche de la résolution numérique. C’est une caractéristique du schéma d’Euler explicite qui est celui que l’on vient de mettre en place : il est conditionnellement stable.
Exercice
Passons à l’expérience numérique : on implémente le schéma d’Euler explicite et on montre l’effet du pas de temps. Pour cela, on reformule les équations algébriques ci-dessus sous forme d’algorithme :
- Résolution de l’équation différentielle discrétisée pour trouver la dérivée : \(\dot{y}_{i} \gets -y_i\)
- Mise-à-jour de la variable \(y_{i+1} \gets y_{i} + \Delta t \dot{y}_i\)
- \(i \gets i+1\)
Implémentez l’algorithme ci-dessus et essayez de reproduire la Figure 1.1, qui montre le résultat de calculs avec différents pas de temps \(\Delta t\). La solution de l’exercice est dans la cellule de code ci-dessous. Commentez les comportements montrées par la Figure 1.1.
Code
def euler_explicit(t_final, dt):
M = int(t_final / dt) + 1
y = np.zeros(M)
t = np.arange(M) * dt
y[0] = 1
for i in range(M-1):
y_dot = -y[i]
y[i+1] = y[i] + dt * y_dot
return t, y
t = np.linspace(0, 9, 200)
y_exact = np.exp(-t)
cmap = plt.get_cmap('copper')
fig, ax = plt.subplots()
ax.plot(t, y_exact, label="Solution", lw=4, color='k')
for dt in [0.5, 1, 1.5, 2.1]:
ax.plot(*euler_explicit(t.max(), dt),
label=f"$\\Delta t = {dt:.1f}$", color=cmap(dt / 2))
ax.legend()
ax.grid()
ax.set(xlabel="$t$",
ylabel="$y(t)$")
plt.show()
Schéma d’Euler implicite
Au lieu de poser l’équation au temps \(t_i\), pour lequel on connaît la valeur de \(y_i\) (puisqu’on la calcule à l’itération précédente), posons l’équation au temps \(t_{i+1}\) auquel on ne connaît pas encore \(y_{i+1}\) (c’est pour cela que le schéma est implicite) : \[ \dot{y}_{i+1} + y_{i+1} = 0,\] on doit à présent approximer \(\dot y_{i+1}\) par différences finies : \[ \dot{y}_{i+1} \approx \frac{y(t_{i+1}) - y(t_i)}{t_{i+1} - t_i} = \frac{y_{i+1} - y_i}{\Delta t}, \] en remplaçant dans l’équation différentielle : \[ \frac{y_{i+1} - y_i}{\Delta t} + y_{i+1} = 0 \quad \Leftrightarrow \quad y_{i+1} = \frac{1}{1 + \Delta t}y_i,\quad i\in[0, M-1] \] Au premier abord, on dirait une situation similaire au schéma précédent : on a une relation de récurrence qui définit une suite géométrique. En revanche, on voit que pour tout \(\Delta t > 0\), la raison de la suite \(1/(1 + \Delta t)\) est plus petite que 1. Le schéma que l’on vient de poser, appelé Euler implicite est donc inconditionnellement stable.
Exercice
Faisons la même expérience que précédemment, en formulant les équations sous forme d’algorithme. Comme on a une formulation implicite, l’écriture de l’algorithme est moins évidente. On doit partir de l’équation de différence finie pour \(\dot{y}_{i+1}\). En combinant avec l’équation différentielle, on peut écrire que :
\[\begin{align*} \dot{y}_{i+1} (= -y_{i+1}) & = \frac{y_{i+1} - y_i}{\Delta t}\\ \Leftrightarrow y_{i+1} & = y_i + \Delta t \dot{y}_{i+1}\\ \Leftrightarrow -\dot{y}_{i+1} &= y_i + \Delta t\dot{y}_{i+1}\\ \Leftrightarrow \dot{y}_{i+1} &= \frac{-1}{1 + \Delta t} y_i \end{align*}\]
On peut alors écrire l’agorithme suivant:
- \(\dot{y}_{i+1}\gets \frac{-1}{1 + \Delta t} y_i\)
- \(y_{i+1} \gets y_i + \Delta t \dot{y}_{i+1}\)
- \(i\gets i + 1\)
Implémentez l’algorithme ci-dessus et essayez de reproduire la Figure 1.2, qui montre des calculs par schéma d’Euler implicite pour différents pas de temps. Commentez la différence avec la Figure 1.1.
Code
def euler_implicit(t_final, dt):
M = int(t_final / dt) + 1
y = np.zeros(M)
t = np.arange(M) * dt
y[0] = 1
for i in range(M-1):
y_dot = -1 / (1 + dt) * y[i]
y[i+1] = y[i] + dt * y_dot
return t, y
t = np.linspace(0, 8, 200)
y_exact = np.exp(-t)
plt.plot(t, y_exact, label="Solution", lw=4, color='k')
for dt in [0.5, 1, 1.5, 2.1]:
plt.plot(*euler_implicit(t.max(), dt),
label=f"$\\Delta t = {dt:.1f}$", color=cmap(dt / 2))
plt.legend()
plt.xlabel("$t$")
plt.ylabel("$y(t)$")
plt.grid()
plt.show()
On observe, certes, que plus \(\Delta t\) est grand, moins la solution numérique est proche de la solution analytique, mais aucune solution numérique n’a de comportement oscillant comme précédemment. On voit également, contrairement à la figure d’Euler explicite, que la pente initiale est fausse pour les solutions approximées.
Améliorer la convergence
On a vu deux façons de mettre à jour une variable avec la méthode d’Euler :
- Euler explicite, qui utilise \(\dot y_{i}\) : \(y_{i+1} \gets y_i + \Delta t y_i\)
- Euler implicite, qui utilise \(\dot y_{i+1}\) : \(y_{i+1} \gets y_i + \Delta t y_{i+1}\)
Le schéma ci-dessus montre graphiquement ces deux façons de mettre à jour la variable :
Code
fig, ax = plt.subplots()
start = 1
end = 2
dt = 1
t = np.linspace(start - 0.1, end + 0.1, 200)
y = np.exp(-t)
annotate_kwargs = dict(xytext=(0, 1), textcoords='offset fontsize', fontsize=16)
ax.plot(t, y, color='k')
ax.plot([start, end], [np.exp(-start), np.exp(-end)], marker='o', color='r', ls='--')
ax.annotate('$y_i$', (start, np.exp(-start)), color='r', **annotate_kwargs)
ax.annotate('$y(t_{i+1})$', (end, np.exp(-end)), color='r', **annotate_kwargs)
# Forward euler update
fwd = (end, np.exp(-start) + dt * (-np.exp(-start)))
ax.axline((start, np.exp(-start)), slope=-np.exp(-start), color='blue', ls='--', lw=1.5)
ax.plot(fwd[0], fwd[1], marker='o', color='blue')
ax.annotate('$y_{i+1}^\\mathrm{explicite}$', fwd, color='b', **annotate_kwargs)
# Backward euler update
bwd = end, np.exp(-start) + dt * (-np.exp(-end))
ax.axline((end, np.exp(-end)), slope=-np.exp(-end), color='green', ls='--', lw=1.5)
ax.axline((start, np.exp(-start)), slope=-np.exp(-end), color='green', ls='--', lw=1.5)
ax.plot(bwd[0], bwd[1], marker='o', color='green')
ax.annotate('$y_{i+1}^\\mathrm{implicite}$', bwd, color='g', **annotate_kwargs)
ax.set(xlabel='$t$', ylabel="$y(t)$")
plt.show()
On remarque qu’aucune des pentes utilisées pour la mise-à-jour n’est satisfaisante, et que pour atteindre le point rouge sur la courbe noire, la solution analytique, il faut faire un pas avec une pente entre la pente initiale \(\dot y_{i}\) et la pente finale \(\dot y_{i+1}\). On propose donc le schéma de mise à jour suivant, avec une interpolation linéaire des deux pentes :
\[ y_{n+1} \gets y_i + \Delta t(\beta \dot y_i + (1 - \beta) \dot y_{i+1}),\quad \beta \in [0, 1]\]
Notons que pour \(\beta = 1\) on a le schéma explicite et \(\beta = 0\) le schéma implicite précédent, mais le schéma est implicite pour toute valeur de \(\beta\) non-nulle, car il faut calculer \(\dot y_{i+1}\). Résolvons pour l’inconnue \(\dot y_{i+1}\), on obtient :
\[ \dot y_{i+1} = \frac{-1}{1 + \Delta t(1 - \beta)}(y_i + \Delta t \beta \dot y_i) \]
On peut alors écrire un algorithme sous la forme d’un schéma de prédiction-correction :
- \(\dot y_i \gets -y_i\)
- \(y_\mathrm{pred} \gets y_i + \Delta t\beta \dot y_i\) (prédiction)
- \(\dot y_{i+1} \gets \frac{-1}{1 + \Delta t(1 - \beta)}y_\mathrm{pred}\)
- \(y_{i+1} \gets y_\mathrm{pred} + \Delta t(1-\beta)\dot y_{i+1}\) (correction)
- \(i \gets i+1\)
Important
La structure “prédiction \(\rightarrow\) mise-à-jour de la dérivée \(\rightarrow\) correction” est très utilisée par les schémas d’intégration en temps. C’est d’ailleurs cette structure que l’on utilisera pour résoudre l’Équation 1.1.
Exercice
Implémentez l’algorithme ci-dessus, et essayez de reproduire la Figure 1.4, qui compare les trois schémas d’intégrations en temps (\(\beta = 0\), \(\beta = 1\) et \(\beta = 0.5\)) pour deux valeurs de \(\Delta t\).
Code
def euler_partial(t_final, dt, beta):
M = int(t_final / dt) + 1
y = np.zeros(M)
t = np.arange(M) * dt
y[0] = 1
for i in range(M-1):
y_dot = -y[i]
y[i+1] = y[i] + dt * beta * y_dot
y_dot = -1 / (1 + dt * (1 - beta)) * y[i+1]
y[i+1] += dt * (1-beta) * y_dot
return t, y
fig, axs = plt.subplots(1, 2, figsize=(8, 4))
t = np.linspace(0, 8, 200)
y_exact = np.exp(-t)
for ax in axs:
ax.plot(t, y_exact, label="Solution", lw=4, color='k')
for ax, dt in zip(axs, [1.5, 0.5]):
ax.plot(*euler_implicit(t.max(), dt), label="$\\beta = 0$")
ax.plot(*euler_explicit(t.max(), dt), label="$\\beta = 1$")
ax.plot(*euler_partial(t.max(), dt, 0.5), label="$\\beta = 0.5$")
ax.legend()
ax.set(xlabel='$t$', ylabel='$y(t)$', title=f'$\\Delta t= {dt:.1f}$')
ax.grid()
fig.tight_layout()
plt.show()
On remarque que l’approximation avec \(\beta = 0.5\) est bien meilleures que les deux schémas précédents. En revanche, cela se fait au prix de deux désavantages :
- L’algorithme est implicite. Ce n’est pas très grave sur un système aussi simple, mais sur de grosses simulations le coût de calcul de \(\dot y_{i+1}\) peut être élevé.
- Malgré le fait que l’algorithme soit implicite, il est conditionnellement stable : selon la valeur de \(\beta\), \(\Delta t\) ne peut pas dépasser une certaine valeur
On a donc gagné en précision en ayant le même coût de calcul que le schéma d’Euler implicite mais sans son avantage principal d’être inconditionnellement stable. C’est à la personne faisant la simulation de choisir l’algorithme adapté !
1.4 Intégration en temps — second ordre
Les approximations par différences finies sont basées sur un développement limité de \(y(t + \Delta t) = y(t) + \dot{y}(t)\Delta t + O(\Delta t^2)\). Essayons d’appliquer ce concept à l’équation du second ordre suivantes :
\[ \ddot y + y = 0,\quad y(0) = 1,\ \dot y(0) = 0,\] qui a pour solution \(y(t) = \cos(t)\). Ici, l’état complet du système n’est plus uniquement déterminé par \(y(t)\), on a aussi besoin de connaître \(\dot y(t)\) pour avoir l’évolution du système.
Schéma d’Euler explicite
Faisons donc un développement limité au premier ordre de \(y\) et \(\dot y\):
\[ \begin{split} y(t + \Delta t) & = y(t) + \dot y(t)\Delta t + O(\Delta t^2),\\ \dot y(t + \Delta t) & = \dot y(t) + \ddot y(t)\Delta t + O(\Delta t^2). \end{split} \tag{1.2}\] Avec une notation indicielle \(y_{i+1}\) et \(\dot y_{i+1}\), on a le schéma d’intégration suivant :
\[\begin{align*} y_{i+1} &= y_i + \dot y_i \Delta t,\\ \dot y_{i+1} & = \dot y_i + \ddot y_i\Delta t = \dot y_i - y_i\Delta t, \end{align*}\] que l’on peut réécrire sous forme matricielle :
\[ \begin{pmatrix}y_{i+1} \\ \dot y_{i+1}\end{pmatrix} = \begin{pmatrix} 1 & \Delta t \\ -\Delta t & 1 \end{pmatrix}\cdot\begin{pmatrix} y_i \\ \dot y_i \end{pmatrix} \]
La stabilité du schéma est donc donnée par les propriétés spectrales de \(A(\Delta t) = \left(\begin{smallmatrix} 1 & \Delta t \\ -\Delta t & 1\end{smallmatrix}\right)\). Un calcul des valeurs propres donne :
\[\lambda_1 = 1 + \mathrm i\Delta t,\ \lambda_2 = 1 - \mathrm i\Delta t,\quad \mathrm{et}\quad |\lambda_1|^2 = |\lambda_2|^2 = 1 + \Delta t^2 > 1.\]
La matrice \(A(\Delta t)\) a donc pour effet de multiplier la norme du vecteur \((y_i, \dot y_i)\) par un facteur toujours plus grand que 1 : le schéma d’intégration est donc inconditionnellement instable !
Avertissement
La stabilité d’un schéma d’intégration est une propriété de la procédure de discrétisation et du système d’équation que l’on essaie de résoudre. Le schéma d’Euler explicite peut être stable pour un système et instable pour un autre !
Velocity-Verlet
Puisque le développement à l’ordre 1 n’est pas suffisant, essayons à un ordre supérieur :
\[ y(t + \Delta t) = y(t) + \dot y(t)\Delta t + \ddot y(t) \frac{\Delta t^2}{2} + \dddot y(t)\frac{\Delta t^3}{6} + O(\Delta t^4),\] on remarque que le signe de \(\Delta t\) n’intervient que dans les termes d’ordre impair : \[ y(t - \Delta t) = y(t) - \dot y(t)\Delta t + \ddot y(t) \frac{\Delta t^2}{2} - \dddot y(t)\frac{\Delta t^3}{6} + O(\Delta t^4),\] et en calculant la différence des deux équations, on a: \[ y(t + \Delta t) - y(t - \Delta t) = 2\dot y(t)\Delta t + O(\Delta t^3),\] on a donc gagné un ordre par rapport à l’expansion précédente (Équation 1.2) !
En revanche, le passage à la notation indicielle est moins aisé puisque l’on a trois temps \(t\), \(t + \Delta t\) et \(t - \Delta t\). Appelons les quantités définies à \(t + \Delta t\) “\(i + 1\)” et celles définies à \(t - \Delta t\) “\(i\)”. Le moment \(t\) se situe donc au milieu d’un intervalle de temps d’une longueur \(2\Delta t\). On notera alors les quantités définies à \(t\) “\(i + \frac{1}{2}\)”. On notera également \(\Delta T = 2\Delta t\) qui est la longueur de l’intervalle de temps entre \(i\) et \(i + 1\). On a alors :
\[\begin{split} \dot y_{i+\frac{1}{2}} & = \dot y_i + \ddot y_i \frac{\Delta T}{2},\\ y_{i+1} & = y_i + \dot y_{i + \frac{1}{2}}\Delta T,\\ \dot y_{i+1} & = \dot y_{i + \frac{1}{2}} + \ddot y_{i+1}\frac{\Delta T}{2}. \end{split} \tag{1.3}\]
La première équation est un développement limité explicite: \[\dot y(t + \frac{\Delta T}{2}) = \dot y(t) + \ddot y(t)\frac{\Delta T}{2} + O(\Delta T^2),\] la seconde d’un développement limité implicite : \[ \dot y(t + \frac{\Delta T}{2}) = \dot y((t + \Delta T) - \frac{\Delta T}{2}) = \dot y(t + \Delta T) - \ddot y(t + \Delta T)\frac{\Delta T}{2} + O(\Delta T^2).\]
On a obtenu l’algorithme Velocity-Verlet (1967) ! Les deux premières équations du système 1.3 sont l’étape de prédiction, et la troisième équation est l’étape de correction. Entre les deux, il faut calculer \(\ddot y_{i+1}\) qui ne dépend que de \(y_{i+1}\).
Exercice
Implémentez l’algorithme Velocity-Verlet pour l’équation linéaire du second ordre \(\ddot y + y = 0\) avec les conditions initiales \(y(0) = 1\), \(\dot y(0) = 0\) et comparez avec la solution analytique pour reproduire la Figure 1.5.
Code
def verlet_prediction(y, y_dot, y_dotdot, dt):
y_half = y_dot + y_dotdot * dt / 2
return y + y_half * dt, y_half
def verlet_correction(y_dot, y_dotdot, dt):
return y_dot + y_dotdot * dt / 2
def compute_ydotdot(y):
return -y
def verlet_step(y, y_dot, y_dotdot, dt):
y, y_dot = verlet_prediction(y, y_dot, y_dotdot, dt)
y_dotdot = compute_ydotdot(y)
y_dot = verlet_correction(y_dot, y_dotdot, dt)
return y, y_dot, y_dotdot
def verlet(t_final, dt):
M = int(t_final / dt) + 1
y = np.zeros(M)
y_dot = np.zeros_like(y)
t = np.arange(M) * dt
y[0] = 1
y_dot[0] = 0
y_dotdot = compute_ydotdot(y[0])
for i in range(M-1):
y[i+1], y_dot[i+1], y_dotdot = verlet_step(y[i], y_dot[i], y_dotdot, dt)
return t, y, y_dot
fig, ax = plt.subplots()
t = np.linspace(0, 4 * np.pi, 200)
y_analytical = np.cos(t)
ax.plot(t, y_analytical, 'k', lw=4, label='Solution')
for dt in [0.5, 1, 2.01]:
t_num, y, _ = verlet(t.max(), dt)
ax.plot(t_num, y, label=f'$\\Delta t = {dt:.2f}$', color=cmap(dt / 2))
ax.set(xlabel='$t$', ylabel='$y(t)$')
ax.legend(ncol=2)
plt.show()
On remarque ici que, comme pour le schéma explicite d’Euler, il y a un pas de temps critique (qui semble être proche de 2). On peut le calculer en faisant la même analyse que ci-dessus, à savoir exprimer le couple \((y_{i+1}, \dot y_{i+1})\) en fonction de \((y_i, \dot y_i)\) :
\[\begin{align*} y_{i+1} &= y_i + y_{i+1/2}\Delta T \\ & = y_i + (\dot y_i + \ddot y_i\frac{\Delta T}{2})\Delta T\\ & = y_i + (\dot y_i - y_i\frac{\Delta T}{2})\Delta T\\ & = (1-\frac{\Delta T^2}{2})y_i + \dot y_i\Delta T\\ \dot y_{i+1} & = \dot y_{i+1/2} + \ddot y_{i+1}\frac{\Delta T}{2}\\ & = \dot y_i + \ddot y_i\frac{\Delta T}{2} + \ddot y_{i+1}\frac{\Delta T}{2}\\ &= \dot y_i -y_i\frac{\Delta T}{2} - y_{i+1}\frac{\Delta T}{2}\\ \Leftrightarrow y_{i+1}\frac{\Delta T}{2} + \dot y_{i+1} & = -y_i \frac{\Delta T}{2} + \dot y_i \end{align*}\]
On peut donc écrire le système suivant :
\[ \underbrace{\begin{pmatrix} 1 & 0\\ \frac{\Delta T}{2} & 1\end{pmatrix}}_A\cdot \begin{pmatrix} y_{i+1}\\ \dot y_{i+1}\end{pmatrix} = \underbrace{\begin{pmatrix}1 - \frac{\Delta T^2}{2} & \Delta T \\ -\frac{\Delta T}{2} & 1\end{pmatrix}}_B \cdot \begin{pmatrix} y_{i}\\ \dot y_{i}\end{pmatrix}\]
Pour que le schéma soit stable, on souhaite que
\[ \begin{Vmatrix}y_{i+1} \\ \dot y_{i+1}\end{Vmatrix} \leq |R|\cdot \begin{Vmatrix} y_i \\ \dot y_i\end{Vmatrix},\] où \(R = A^{-1}B\). Il faut donc étudier les valeurs propres de la matrice \(A^{-1}B\) :
\[ R = A^{-1}B = \begin{pmatrix} 1 & 0\\ -\frac{\Delta T}{2} & 1\end{pmatrix}\cdot \begin{pmatrix}1 - \frac{\Delta T^2}{2} & \Delta T \\ -\frac{\Delta T}{2} & 1\end{pmatrix} = \begin{pmatrix}1 - \frac{\Delta T^2}{2} & \Delta T\\-\frac{\Delta T}{2}(2 - \frac{\Delta T^2}{2}) & 1 - \frac{\Delta T}{2}\end{pmatrix}.\]
En calculant le polynôme caractéristique, on obtient l’équation suivante pour les valeurs propres de \(R\): \[\left(1 - \frac{\Delta T^2}{2} - \lambda\right) + \frac{\Delta T^2}{2}\left(2- \frac{\Delta T^2}{2}\right) = 0\]
Ce qui donne les valeurs propres \(\lambda_1 = 1 - \frac{\Delta T^2}{2} + \mathrm i \frac{\Delta T^2}{2}(2 - \frac{\Delta T^2}{2})\) et \(\lambda_2 = 1 - \frac{\Delta T^2}{2} - \mathrm i \frac{\Delta T^2}{2}(2 - \frac{\Delta T^2}{2})\). Pour que les valeurs absolues des deux valeurs propres soient plus petites que 1, il faut que :
\[\begin{align*} |\lambda_{1,2}|^2 = 1 + \left(2 - \frac{\Delta T^2}{2}\right)\frac{\Delta T^2}{2}\left(\frac{\Delta T^2}{2}\left(2 - \frac{\Delta T^2}{2}\right) - 1\right) & \leq 1\\ \Rightarrow \left(2 - \frac{\Delta T^2}{2}\right)\left(\frac{\Delta T^2}{2}\left(2 - \frac{\Delta T^2}{2}\right) - 1\right) & \leq 0\\ \Rightarrow -\left(2 - \frac{\Delta T^2}{2}\right)\left(\frac{\Delta T^2}{2} - 1\right)^2 & \leq 0\\ \Rightarrow 2 - \frac{\Delta T^2}{2} & \geq 0\\ \Rightarrow |\Delta T| & \leq 2 \end{align*}\]
On a donc trouvé le pas de temps critique \(\Delta T_\mathrm{crit} = 2\).
Dimensionnalisation
L’équation que l’on a résolue est l’équation adimensionnalisée pour le système vibratoire masse-ressort sans amortissement. Avec une masse \(m\) et une raideur de ressort \(k\), l’équation s’écrit :
\[ m\ddot u(t) + ku(t) = 0 \quad \Leftrightarrow \quad \ddot u(t) + \omega^2 u(t) = 0,\text{ avec } \omega^2 = \frac{k}{m}.\]
Par changement de variable \(\hat t = \omega t\), on peut définir \(u(t) = y(\hat t)\) et \(\ddot u(t) = \omega^2 y(\hat t)\) et l’on obtient l’équation adimensionnalisée \(\ddot y + y = 0\). Le pas de temps critique \(\Delta T_\mathrm{crit}\) est donc lui aussi adimensionnalisé. Notons \(\Delta t_\mathrm{crit}\) le pas de temps critique dimensionnalisé (tel que \(\Delta T_\mathrm{crit} = \omega\Delta t_\mathrm{crit}\)), pour que Velocity-Verlet soit stable quand on résoud l’équation dimensionnalisée en \(u\), il faut que
\[ \Delta t \leq \Delta t_\mathrm{crit} = \frac{2}{\omega}\]
On voit donc que plus la fréquence propre \(\omega\) est grand, plus le pas de temps critique est faible.
1.5 Application à l’équation de Newton
Comme pour les équations scalaires vues ci-dessus, l’algorithme Velocity-Verlet consiste en la mise-à-jour des dérivées d’ordre 0 et 1 de la position, il se formule ainsi :
- \(\vec v_{i+\frac{1}{2}} \gets \vec v_i + \frac{\Delta t}{2}\left(\frac{\vec f_i}{m}\right)\)
- \(\vec r_{i+1} \gets \vec r_i + \Delta t \vec v_{i+\frac{1}{2}}\)
- Calcul de \(\vec f_{i+1}\)
- \(\vec v_{i+1} \gets \vec v_{i+\frac{1}{2}} + \frac{\Delta t}{2}\left(\frac{\vec f_{i+1}}{m}\right)\)
Ici, on a formulé l’algorithme pour une seule particule. La généralisation à \(N\) particules est immédiate : la trajectoire d’une particule ne dépend que de la force qui s’applique sur cette particule. On peut donc considérer que le mouvement des particules est complètement indépendant pour les étapes 1, 2 et 4, et \(\vec r\), \(\vec v\), \(\vec f\) sont la position, vitesse et force d’une particule parmi \(N\). En revanche, l’étape 3 nécessite de calculer la force s’appliquant sur chaque particule, qui peut donc dépendre de la position des autres particules. C’est ce que l’on détaillera dans le prochain chapitre.
Le TP 1 est une application de Velocity-Verlet à la balistique des particules uniformément accélérées.
1.6 Lectures conseillées
- Bonnet et Frangi (2006), chapitre 9
- Fiorelli Vilmart et Vilmart (2017)
- Molecular Dynamics, chapter 2
Bonnet, Marc, et Attilio Frangi. 2006. Analyse Des Solides Déformables Par La Méthode Des Éléments Finis. Mécanique. Palaiseau : [Paris]: Éditions de l’École polytechnique ; Diffusion, Ellipses.
Fiorelli Vilmart, Shaula, et Gilles Vilmart. 2017. « Computing the Long Term Evolution of the Solar System with Geometric Numerical Integrators ». Snapshots of Modern Mathematics from Oberwolfach;2017, 09. https://doi.org/10.14760/SNAP-2017-009-EN.
Verlet, Loup. 1967. « Computer "Experiments" on Classical Fluids. I. Thermodynamical Properties of Lennard-Jones Molecules ». Physical Review 159 (1): 98‑103. https://doi.org/10.1103/PhysRev.159.98.
Des solutions particulières existent, mais le cas général est un parfait exemple de système chaotique↩︎