USMA1Q — Méthodes Numériques¶
Séance 11 — Équations différentielles ordinaires II : systèmes et schémas mécaniques¶
Conservatoire National des Arts et Métiers
01 Avril 2026
nicholas-anton.collins-craft@enpc.fr
Un lien vers Google Colab qui peut faire tourner le notebook qui contient les exercices se trouve ici :
https://colab.research.google.com/github/nickcollins-craft/USMA1Q-Methodes-Numeriques/blob/main/
Plan du cours¶
| Jour | Sujet |
|---|---|
| 03 mars | Fondamentaux de la programmation en Python |
| 04 mars (matin) | Fonctions, listes approfondies, fichiers et graphiques |
| 04 mars (après-midi) | Calcul numérique avec NumPy |
| 11 mars (matin) | Tableaux en mémoire, matrices & précision numérique |
| 11 mars (après-midi) | Résolution des systèmes linéaires |
| 18 mars (matin) | Interpolation & recherche des racines |
| 18 mars (après-midi) | Systèmes d'équations non linéaires |
| 25 mars (matin) | Intégration numérique I |
| 25 mars (après-midi) | Intégration numérique II (Monte Carlo) |
| 01 avril (matin) | Équations différentielles ordinaires I |
| 01 avril (après-midi) | Équations différentielles ordinaires II |
| 08 avril | Examen |
Plan de la séance (indicatif)¶
| Horaire | Sujet |
|---|---|
| 0:00 – 0:55 | Systèmes d'EDO |
| 0:55 – 1:50 | Schémas d'intégration pour la mécanique |
| 1:50 – 2:00 | ☕ Pause |
| 2:00 – 3:00 | Examen blanc |
| 3:00 – 3:40 | Mini-projet encadré |
| 3:40 – 3:45 | Récapitulatif |
1 · Systèmes d'EDO¶
Rappel séance 10 : nous avons résolu des EDO scalaires $$\frac{\mathrm{d}y}{\mathrm{d}t} = f(t, y), \quad y(t_0) = y_0$$
Maintenant : que faire si l'inconnue est un vecteur $\mathbf{y}(t) \in \mathbb{R}^n$ ?
$$\frac{\mathrm{d}\mathbf{y}}{\mathrm{d}t} = \mathbf{f}(t, \mathbf{y}), \quad \mathbf{y}(t_0) = \mathbf{y}_0$$
Toutes les méthodes vues en séance 10 (Euler, RK4, solve_ivp) s'appliquent sans modification : il suffit que f renvoie un vecteur.
Exemples en science des matériaux :
- Cinétiques couplées (plusieurs phases en compétition)
- Mécanique des structures (déplacement + vitesse)
- Dynamique moléculaire simplifiée (positions + impulsions)
Réduction d'ordre : EDO d'ordre 2 → système du 1er ordre¶
Tout problème d'ordre supérieur se ramène à un système du premier ordre.
Exemple : l'oscillateur amorti $$m\ddot{x} + c\dot{x} + kx = 0, \quad x(0) = x_0,\; \dot{x}(0) = v_0$$
Changement de variables : $$y_1 = x \quad (\text{déplacement}), \qquad y_2 = \dot{x} \quad (\text{vitesse})$$
Système équivalent : $$\frac{\mathrm{d}}{\mathrm{d}t}\begin{pmatrix} y_1 \\ y_2 \end{pmatrix} = \begin{pmatrix} y_2 \\ \dfrac{-c\,y_2 - k\,y_1}{m} \end{pmatrix}$$
Le second membre $\mathbf{f}(t, \mathbf{y})$ est un vecteur à deux composantes.
Paramètres de l'oscillateur : taux d'amortissement $\zeta$¶
On réécrit l'équation sous forme adimensionnelle : $$\ddot{x} + 2\zeta\omega_0\dot{x} + \omega_0^2 x = 0$$
avec $\omega_0 = \sqrt{k/m}$ (pulsation propre) et $\zeta = c/(2\sqrt{mk})$ (taux d'amortissement).
| Régime | Valeur de $\zeta$ | Comportement |
|---|---|---|
| Sous-amorti | $\zeta < 1$ | Oscillations décroissantes |
| Critique | $\zeta = 1$ | Retour à zéro le plus rapide sans oscillation |
| Sur-amorti | $\zeta > 1$ | Retour exponentiel sans oscillation |
Implémentation avec solve_ivp¶
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
m, k = 1.0, 4.0 # masse [kg], raideur [N/m]
omega0 = np.sqrt(k / m) # pulsation propre [rad/s]
c = 0.5 # amortissement [N·s/m] → zeta = c/(2*sqrt(m*k))
def rhs_oscillateur(t, y):
# Second membre du système : y = [x, v].
x, v = y
dxdt = v
dvdt = (-c * v - k * x) / m
return [dxdt, dvdt]
# Conditions initiales
y0 = [1.5, 0.0] # x(0) = 1.5 m, v(0) = 0 m/s
sol = solve_ivp(rhs_oscillateur, [0, 20], y0, method='RK45', t_eval=np.linspace(0, 20, 1000), rtol=1e-9)
print("Composante y[0] = x(t), y[1] = v(t)")
print("Dernier point :", sol.t[-1], sol.y[:, -1])
Composante y[0] = x(t), y[1] = v(t) Dernier point : 20.0 [-0.00292262 -0.0186301 ]
Portraits de phase — les trois régimes¶
Le portrait de phase trace $v = \dot{x}$ en fonction de $x$ ; il révèle la structure géométrique de la solution sans regarder le temps.
- Sous-amorti ($\zeta < 1$) : spirale convergente vers l'origine
- Critique ($\zeta = 1$) : trajectoire directe, sans boucle
- Sur-amorti ($\zeta > 1$) : convergence plus lente, le long des directions propres
Oscillateur forcé — résonance¶
Ajoutons un forçage harmonique : $$m\ddot{x} + c\dot{x} + kx = F_0\sin(\omega t)$$
Le système augmenté reste de la forme $\mathbf{f}(t, \mathbf{y})$ :
$$\frac{\mathrm{d}}{\mathrm{d}t}\begin{pmatrix} y_1 \\ y_2 \end{pmatrix} = \begin{pmatrix} y_2 \\ \dfrac{F_0\sin(\omega t) - c\,y_2 - k\,y_1}{m} \end{pmatrix}$$
La résonance se produit quand $\omega \approx \omega_0$ : l'amplitude à l'état stationnaire est maximale.
F0, omega_force = 1.0, 2.0 # amplitude [N], pulsation de forçage [rad/s]
def rhs_force(t, y):
x, v = y
dvdt = (F0 * np.sin(omega_force * t) - c * v - k * x) / m
return [v, dvdt]
sol_force = solve_ivp(rhs_force, [0, 100], [0.0, 0.0], method='RK45', t_eval=np.linspace(0, 100, 5000))
À vous de faire les exercices 1.1 et 1.2 !¶
2 · Schémas d'intégration pour la mécanique¶
Contexte : en mécanique des structures, on discrétise spatialement (méthode des éléments finis) et on obtient un système semi-discret d'EDO du 2e ordre :
$$\mathbf{M}\ddot{\mathbf{u}} + \mathbf{C}\dot{\mathbf{u}} + \mathbf{K}\mathbf{u} = \mathbf{F}(t)$$
- $\mathbf{M}$ : matrice de masse
- $\mathbf{C}$ : matrice d'amortissement (souvent Rayleigh : $\mathbf{C} = \alpha\mathbf{M} + \beta\mathbf{K}$ pour $\alpha$ et $\beta$ constants scalaires).
- $\mathbf{K}$ : matrice de rigidité
- $\mathbf{F}(t)$ : chargement extérieur
Objectif : intégrer ce système directement (sans réduire à un système du 1er ordre) avec des schémas exploitant la structure physique du problème.
La méthode $\theta$¶
Forme unifiée qui englobe plusieurs schémas classiques.
$$\mathbf{u}_{n+1} = \mathbf{u}_n + h\,[(1-\theta)\mathbf{v}_n + \theta\,\mathbf{v}_{n+1}]$$ $$\mathbf{a}_{n+\theta}=\frac{\mathbf{v}_{n+1}-\mathbf{v}_{n}}{h}$$ $$\mathbf{M}\mathbf{a}_{n+\theta} + \mathbf{C}\mathbf{v}_{n+\theta} + \mathbf{K}\mathbf{u}_{n+\theta} = \mathbf{F}_{n+\theta}$$
| Valeur de $\theta$ | Schéma | Stabilité |
|---|---|---|
| $\theta = 0$ | Euler explicite | Conditionnellement stable |
| $\theta = 1/2$ | Point médian implicit | Inconditionnellement stable |
| $\theta = 1$ | Euler implicite | Inconditionnellement stable |
Pour $\theta \geq 1/2$ : inconditionnellement stable — le pas de temps n'est limité que par la précision souhaitée, pas par la stabilité. Si le système est linéaire, $\theta=1/2$ coincide avec la méthode de Crank–Nicolson.
Méthode de Crank–Nicolson¶
La dérivée est approximée par la moyenne des pentes au début et à la fin du pas.
$$\mathbf{u}_{n+1} = \mathbf{u}_n + \frac{h}{2}(\mathbf{v}_n + \mathbf{v}_{n+1})$$ $$\mathbf{M}\frac{\mathbf{v}_{n+1} - \mathbf{v}_n}{h} + \mathbf{C}\frac{\mathbf{v}_{n+1} + \mathbf{v}_n}{h} + \mathbf{K}\frac{\mathbf{u}_n + \mathbf{u}_{n+1}}{2}=\frac{\mathbf{F}_{n+1}+\mathbf{F}_{n}}{2}$$
Propriétés :
- Ordre 2 en temps (même ordre que RK2, mais structure adaptée au problème)
- Conserve exactement l'énergie pour un système non amorti ($\mathbf{C} = 0$) si le système est linéaire
- Inconditionnellement stable pour tout $h > 0$
N.B. la moyenne des valeurs au début et à la fin n'est pas nécessairement équivalent au valeur de point au milieu ! Crank–Nicolson et la méthode $\theta$ départ de deux principes différents, qui ne donne les mêmes résultats que dans le cas d'un système linéaire.
Méthode de Newmark-$\beta$¶
Standard industriel en dynamique des structures.
$$\mathbf{u}_{n+1} = \mathbf{u}_n + h\mathbf{v}_n + h^2\left[(\tfrac{1}{2}-\beta)\mathbf{a}_n + \beta\,\mathbf{a}_{n+1}\right]$$ $$\mathbf{v}_{n+1} = \mathbf{v}_n + h\left[(1-\gamma)\mathbf{a}_n + \gamma\,\mathbf{a}_{n+1}\right]$$ $$\mathbf{M}\mathbf{a}_{n+1} + \mathbf{C}\mathbf{v}_{n+1} + \mathbf{K}\mathbf{u}_{n+1} = \mathbf{F}_{n+1}$$
avec les paramètres $\beta \in [0, 1/2]$ et $\gamma \in [0, 1]$.
| $\beta$ | $\gamma$ | Schéma | Stabilité |
|---|---|---|---|
| $0$ | $1/2$ | Différences centrales (explicite) | Conditionnelle |
| $1/4$ | $1/2$ | Accélération moyenne (trapezoïdal) | Inconditionnelle |
| $1/6$ | $1/2$ | Accélération linéaire | Inconditionnelle |
Stabilité de Newmark — interprétation¶
Pour $\gamma = 1/2$ et $\beta \geq 1/4$ : inconditionnellement stable quelle que soit la fréquence du système.
Pour $\beta < 1/4$ (dont différences centrales $\beta = 0$) : conditionnellement stable, la condition est $h \leq 2/\omega_{\max}$ où $\omega_{\max}$ est la plus haute fréquence propre du système discrétisé.
En MEF : $\omega_{\max}$ est déterminé par le maillage — raffiner le maillage augmente $\omega_{\max}$ et force un plus petit pas de temps pour les schémas explicites. C'est pourquoi les codes de choc/impact utilisent souvent des schémas explicites mais avec un grand soin dans le choix du pas.
Implémentation de Newmark ($\beta$, $\gamma$)¶
def newmark_step(M, C, K, F_n, F_n1, u_n, v_n, a_n, h, beta=0.25, gamma=0.5):
# Un pas de Newmark-beta (système matriciel quelconque).
# Prédicteurs
u_pred = u_n + h * v_n + h**2 * (0.5 - beta) * a_n
v_pred = v_n + h * (1 - gamma) * a_n
# Résolution pour a_{n+1}: (M + gamma*h*C + beta*h^2*K) @ a_{n+1} = ...
K_eff = M + gamma * h * C + beta * h**2 * K
rhs = F_n1 - C @ v_pred - K @ u_pred
a_n1 = np.linalg.solve(K_eff, rhs)
# Correcteurs
v_n1 = v_pred + h * gamma * a_n1
u_n1 = u_pred + h**2 * beta * a_n1
return u_n1, v_n1, a_n1
# Exemple : oscillateur scalaire m=1, k=4, c=0 → M=[1], K=[4], C=[0]
import numpy as np
M, C, K = np.array([[1.0]]), np.array([[0.0]]), np.array([[4.0]])
h, T = 0.3, 20
N = int(T / h)
u, v = np.array([1.0]), np.array([0.0])
a = np.linalg.solve(M, -K @ u - C @ v)
traj = [u[0]]
for i in range(N):
u, v, a = newmark_step(M, C, K, np.array([0.]), np.array([0.]),
u, v, a, h)
traj.append(u[0])
print("Max |x|:", max(abs(x) for x in traj)) # doit rester proche de 1
Max |x|: 1.0
Méthode generalized-$\alpha$¶
Extension de Newmark permettant d'introduire une dissipation numérique contrôlée sur les hautes fréquences sans dégrader la précision sur les basses fréquences.
Paramétrisation par le rayon spectral à haute fréquence $\rho_\infty \in [0, 1]$ :
$$\alpha_m = \frac{2\rho_\infty - 1}{\rho_\infty + 1}, \quad \alpha_f = \frac{\rho_\infty}{\rho_\infty + 1}$$ $$\gamma = \tfrac{1}{2} - \alpha_m + \alpha_f, \quad \beta = \tfrac{1}{4}(1 - \alpha_m + \alpha_f)^2$$
| $\rho_\infty$ | Comportement |
|---|---|
| $1$ | Pas de dissipation (Newmark $\beta = 1/4$) |
| $0$ | Dissipation maximale (hautes fréquences annulées) |
Symplecticité — pourquoi ça compte¶
Un Hamiltonien décrit un système sans échange d'énergie avec l'extérieur : $$H(q, p) = \frac{p^2}{2m} + V(q) = \text{const le long des trajectoires}$$
Théorème de Liouville : le système hamiltonien préserve le volume dans l'espace des phases.
Un schéma numérique est symplectique si son application discrète préserve aussi (la version discrète de) cette structure.
Conséquence pratique :
- Les schémas non symplectiques (Euler explicite/implicite, RK4 ordinaire) font dériver l'énergie sur de longues intégrations.
- Les schémas symplectiques (Störmer–Verlet, Runge–Kutta partitionnés) conservent une énergie modifiée exactement → pas de dérive séculaire.
Schéma de Störmer–Verlet (leapfrog/saute-mouton)¶
Le plus simple des intégrateurs symplectiques du 2e ordre :
$$\mathbf{u}_{n+1} = \mathbf{u}_n + h\,\mathbf{v}_n + \frac{h^2}{2}\,\mathbf{a}_n$$ $$\mathbf{a}_{n+1} = \mathbf{M}^{-1}\left[\mathbf{F}(t_{n+1}) - \mathbf{K}\mathbf{u}_{n+1}\right]$$ $$\mathbf{v}_{n+1} = \mathbf{v}_n + \frac{h}{2}\,(\mathbf{a}_n + \mathbf{a}_{n+1})$$
Propriétés :
- Ordre 2 en temps
- Symplectique : préserve l'aire dans l'espace des phases
- Conditionnellement stable ($h < 2/\omega_0$), mais erreur d'énergie bornée (oscillations autour de la valeur exacte plutôt que dérive)
# Comparaison longue intégration : Euler explicite vs Störmer–Verlet
m, k = 1.0, 4.0
x0, v0 = 1.0, 0.0
E0 = 0.5 * k * x0**2 + 0.5 * m * v0**2
h, T = 0.3, 60
N = int(T / h)
t_arr = np.linspace(0, T, N + 1)
# Euler explicite
x_e, v_e = np.zeros(N + 1), np.zeros(N + 1)
x_e[0], v_e[0] = x0, v0
for n in range(N):
a = -k/m * x_e[n]
x_e[n + 1] = x_e[n] + h*v_e[n]
v_e[n + 1] = v_e[n] + h*a
# Störmer–Verlet
x_sv, v_sv = np.zeros(N + 1), np.zeros(N + 1)
x_sv[0], v_sv[0] = x0, v0
for n in range(N):
a_n = -k/m * x_sv[n]
x_sv[n + 1] = x_sv[n] + h*v_sv[n] + 0.5*h**2*a_n
a_n1 = -k/m * x_sv[n + 1]
v_sv[n + 1] = v_sv[n] + 0.5*h*(a_n + a_n1)
E_e = 0.5*k*x_e**2 + 0.5*m*v_e**2
E_sv = 0.5*k*x_sv**2 + 0.5*m*v_sv**2
print("Dérive énergie Euler :", (E_e[-1] - E0)/E0)
print("Dérive énergie Störmer :", (E_sv[-1] - E0)/E0)
Dérive énergie Euler : 5.102484264144263e+26 Dérive énergie Störmer : -0.03253404735298737
Non linéarité¶
Notons qu'on a fait des assomptions implicites que notre système $$\mathbf{M}\ddot{\mathbf{u}} + \mathbf{C}\dot{\mathbf{u}} + \mathbf{K}\mathbf{u} = \mathbf{F}(t)$$ est linéaire, et donc tout les systèmes peuvent être resolus en faisant un réarrangement des variables pour mettre la variable "de base" inconnue (e.g. $\mathbf{a}_{n+1}$, $\mathbf{v}_{n+1}$ etc., selon le schéma) dans le premier membre, et résoudre un système linéaire.
Tout les schémas marchent également quand on a des systèmes nonlinéaires e.g. $$\mathbf{M}(\mathbf{u})\ddot{\mathbf{u}} + \mathbf{C}(\mathbf{u},\dot{\mathbf{u}})\dot{\mathbf{u}} + \mathbf{K}(\mathbf{u},\dot{\mathbf{u}})\mathbf{u} = \mathbf{F}(t)$$ Même si c'est rare d'avoir un $\mathbf{M}$ non linéaire, $\mathbf{C}$ et surtout $\mathbf{K}$ le sont fréquemment (dans la plasticité pour exemple). Dans tel cas, le réarrangement des variables nous donne un système non linéaire de notre variable de base, et nous le résoudre avec des outils comme le Newton-Raphson.
Tableau récapitulatif des schémas¶
| Schéma | Ordre | Stabilité | Symplectique | Usage typique |
|---|---|---|---|---|
| Euler explicite | 1 | Cond. | Non | Pédagogie |
| Euler implicite | 1 | Incond. | Non | Problèmes raides |
| Point milieu implicit $\theta = 1/2$ | 2 | Incond. | Oui | Diffusion, chaleur, MEF |
| Newmark $\beta=1/4, \gamma=1/2$ | 2 | Incond. | Non | Dynamique des structures |
| Störmer–Verlet | 2 | Cond. | Oui | DM, hamiltoniens |
| Generalized-$\alpha$ | 2 | Incond. | Non | MEF avec dissipation |
À vous de faire les exercices 2.1 et 2.2 !¶
☕ Pause — 10 minutes¶
3 · Examen blanc¶
Il dure 1 heure — les notes (le cheat sheet) sont autorisées.
Objectifs :
- S'entraîner dans les conditions de l'examen final du 08 avril
- Identifier les points à retravailler
4 · Mini-projet — Analyse vibratoire d'une structure composite¶
Scénario : Le comportement dynamique d'un panneau de structure doit être analysé. On doit :
- Obtenir les velocités et les déplacements à partir des accelérations et les utiliser pour identifier la fréquence propre et le taux d'amortissement;
- Calculer la rigidé et l'utiliser pour modéliser un système avec une dégrée de liberté;
- Trouvez les fréquences propres d'un système avec deux dégrées de liberté et modéliser avec deux schémas d'intégration différentes.
Données fournies :
- Acceleration et temps de mésure.
Rendu :
Le soir (avant minuit) du 15 avril.
Points clés de la séance¶
| Concept | Message à retenir |
|---|---|
| Systèmes d'EDO | Réduire l'ordre → solve_ivp fonctionne sans modification |
| Énergie mécanique | Euler explicite gonfle, implicite dissipe, Störmer–Verlet et point au milieu implicite conservent |
| Newmark-$\beta$ | Standard MEF ; $\beta=1/4, \gamma=1/2$ inconditionnellement stable |
| Generalized-$\alpha$ | Dissipation contrôlée des hautes fréquences numériques |
| Symplecticité | Préservation de la structure hamiltonienne → pas de dérive séculaire |
Prochain rendez-vous : examen final le 08 avril (1h45) — toutes les séances 1–11 au programme.