Trois Méthodes, Trois Régions
Pour l'équation test dy/dx = λy, trois méthodes EDO explicites ont les régions de stabilité suivantes dans le plan complexe hλ :
Méthode d'Euler (premier ordre) : la région de stabilité est le disque |1 + hλ| ≤ 1, un cercle de rayon 1 centré à (-1, 0). Les valeurs réelles négatives de hλ doivent se trouver dans [-2, 0].
Runge-Kutta 2 (méthode du point milieu) (second ordre) : la région de stabilité est |1 + hλ + (hλ)²/2| ≤ 1. Plus grande que le disque d'Euler, mais toujours bornée.
Runge-Kutta 4 (quatrième ordre) : la région de stabilité satisfait |1 + hλ + (hλ)²/2 + (hλ)³/6 + (hλ)⁴/24| ≤ 1. Les valeurs réelles négatives de hλ s'étendent à environ -2,785. La région est sensiblement plus grande que celle d'Euler.
Euler Rétrograde (implicite) : la région de stabilité est l'ensemble du plan complexe sauf le disque |1 - hλ|⁻¹ > 1, équivalemment |1/(1-hλ)| ≤ 1. Pour λ dans le demi-plan gauche (Re(λ) < 0), c'est inconditionnellement stable — aucune contrainte sur h en matière de stabilité.
La Fonction d'Amplification
Pour toute méthode de Runge-Kutta, le facteur d'amplification par étape R(hλ) est une approximation polynomiale de e^(hλ) :
- Euler : R(z) = 1 + z (tronqué au degré 1)
- RK2 : R(z) = 1 + z + z²/2 (tronqué au degré 2)
- RK4 : R(z) = 1 + z + z²/2 + z³/6 + z⁴/24 (tronqué au degré 4)
La région de stabilité est {z ∈ ℂ : |R(z)| ≤ 1}. L'amplification de la solution vraie : |e^z| = e^(Re(z)). Pour Re(z) < 0 (EDO stable), la solution vraie décroît. La méthode numérique est stable si |R(z)| ≤ 1 — correspondant au comportement décroissant.
Valeurs Propres Purement Imaginaires : Systèmes Oscillatoires
De nombreux systèmes physiques ont des valeurs propres purement imaginaires : λ = iω (oscillations sans amortissement). Le système ressort-masse, la mécanique orbitale, la dynamique du pendule.
Pour λ = iω : hλ = ihω se trouve sur l'axe imaginaire.
Stabilité d'Euler sur l'axe imaginaire : |1 + ihω|² = 1 + (hω)² > 1 pour tout h > 0. Euler est instable pour toute taille de pas sur des valeurs propres purement imaginaires. L'« oscillation » calculée croît sans limites.
Stabilité de RK4 sur l'axe imaginaire : la région de stabilité s'étend à environ |hω| ≤ 2,83 sur l'axe imaginaire. Pour h suffisamment petit, RK4 gère les oscillations non amorties. Euler ne peut pas.
Cette géométrie explique pourquoi Euler échoue sur les systèmes conservatifs (ressort-masse, orbites, équations des ondes) même avec petit h, tandis que RK4 les gère bien.
La Géométrie des Problèmes Raides
Un système d'EDO raide a des valeurs propres de magnitudes très différentes. Le ratio de raideur : κ = max|Re(λᵢ)| / min|Re(λᵢ)| >> 1.
Pourquoi la raideur est coûteuse pour les solveurs explicites :
La stabilité exige h·max|λᵢ| ≤ C (où C dépend de la méthode). La valeur propre la plus négative fixe la limite.
La précision pour la dynamique lente exige h·min|λᵢ| ≥ ε (résoudre le mode le plus lent adéquatement).
Si κ est grand, ces deux exigences forcent un très petit h : assez petit pour la stabilité du mode rapide, assez grand pour échantillonner le mode lent. Le nombre d'étapes s'adapte à κ.
Image géométrique dans le spectre des valeurs propres : les valeurs propres du jacobien ∂f/∂y forment un ensemble de points dans le plan complexe. La région de stabilité d'un solveur explicite doit contenir tous les points h·λᵢ. Si les valeurs propres s'étendent de -1 à -1000, la région de stabilité doit couvrir une plage de 1000 le long de l'axe réel — exigeant un très petit h.
Solveurs implicites : la région de stabilité d'Euler rétrograde couvre l'ensemble du demi-plan gauche. Toutes les valeurs propres avec Re(λ) < 0 sont automatiquement à l'intérieur de la région de stabilité, quel que soit h. La contrainte sur h ne vient que de la précision, pas de la stabilité.
Ratio de Raideur et Coût
Considérez un réseau de réactions chimiques avec des réactions rapides (échelle de temps 10⁻⁶ s) et des réactions lentes (échelle de temps 1 s).
Ratio de raideur : κ = 10⁶ / 1 = 10⁶.
Avec RK4 (limite de stabilité h·|λ_max| ≤ 2,785) : h_max = 2,785 / 10⁶ ≈ 2,8 × 10⁻⁶ s.
Pour intégrer sur 10 s de temps de réaction : étapes = 10 / (2,8 × 10⁻⁶) ≈ 3,6 × 10⁶.
Avec Euler rétrograde (inconditionnellement stable) : h peut être choisi pour la précision des réactions lentes. h = 10⁻² s (100 échantillons sur l'échelle de 1 s). Étapes = 10 / 10⁻² = 1000.
Ratio de coût : 3,6 millions d'étapes explicites contre 1000 étapes implicites — un facteur de 3600. Chaque étape implicite nécessite de résoudre un système linéaire (le coût par étape est plus élevé), mais le coût total est beaucoup plus faible pour les problèmes très raides.
Pourquoi les Tubes n-Dimensionnels Ne Sont Pas Ce Que Vous Pensez
En 2D, un « tube » de rayon ε autour d'une courbe C est l'ensemble des points à distance ε de C. La section transversale est un cercle de rayon ε. Le volume du tube croît proportionnellement à sa longueur.
En n dimensions, la géométrie du tube change fondamentalement, en raison du phénomène du Chapitre 9 :
Le paradoxe des coins n-dimensionnels : dans l'espace n-dimensionnel, presque tout le volume d'un hypercube n-dimensionnel se trouve dans les coins — pas dans la région centrale. Quand n augmente, la fraction de volume à distance ε du centre tend vers zéro pour tout ε fixe.
Appliqué aux tubes de solution d'EDO :
En 2D : si la solution vraie passe par le centre d'un tube, la plupart des points proches se trouvent à proximité de la courbe. Les petites perturbations vous gardent près de la solution vraie.
En haute dimension : la plupart des points à l'intérieur de la boîte englobante du tube sont en fait loin de la courbe de solution vraie. Le « volume » du tube est dominé par les coins — des régions qui sont loin du centre dans plusieurs dimensions simultanément.
Conséquence pour la simulation : avec 28 EDO couplées (le problème d'interception de la Marine de Hamming), une perturbation de taille ε dans chaque dimension peut produire un déplacement total de ε√28 ≈ 5,3ε de la solution vraie. Le tube doit être compris en termes de norme L2 sur toutes les dimensions, pas seulement le déplacement maximal dans une dimension quelconque.
Stabilité en haute dimension : un système où chaque composante décroît indépendamment (chaque valeur propre a une partie réelle négative) peut toujours montrer des déplacements combinés importants car les erreurs des composantes s'ajoutent en norme L2. Le tube 28-dimensionnel n'est pas simplement 28 tubes 1-dimensionnels indépendants — la géométrie les couple.
De la Géométrie à la Conception
Les aperçus géométriques des Chapitres 18-20 s'assemblent comme un ensemble de principes de conception pour la simulation numérique :
Sélection de la taille du pas : h doit placer h·λ à l'intérieur de la région de stabilité pour chaque valeur propre. Pour les systèmes raides, les méthodes implicites suppriment la contrainte de stabilité, ne laissant que les exigences de précision.
Accumulation d'erreur en haute dimension : l'erreur globale est un vecteur dans l'espace n-dimensionnel. Sa norme croît comme √n fois l'erreur par composante. Les simulations haute-dimensionnelles nécessitent des exigences de précision par étape plus strictes.
Feedback comme stabilisateur : si la simulation incorpore un feedback (la sortie calculée influence les entrées ultérieures, comme dans un système de guidage), le feedback convergent amortit les erreurs. La simulation peut tolérer des entrées imprécises pour les quantités à l'intérieur de la boucle de feedback.
L'instabilité comme signal : pour les problèmes avec des champs de direction divergents, l'instabilité peut être exploitée : la direction de divergence porte l'information sur l'erreur de condition initiale, permettant un ajustement correctif.