Problèmes sévères et solveurs rigides : lorsque les méthodes explicites échouent

Une simulation dépendante du temps peut sembler bien formulée et obliger toujours un solveur explicite à franchir des milliers de mesures extrêmement petites. Augmenter la commande d’Euler vers l’avant vers RK4 peut retarder le problème, mais cela ne le supprime pas toujours. Le système peut être rigide.

La rigidité apparaît lorsque la stabilité numérique nécessite un pas de temps beaucoup plus petit que l’étape nécessaire pour résoudre le comportement d’intérêt scientifique. Les composants rapides peuvent se décomposer presque immédiatement, tandis que des composants plus lents contrôlent la solution à long terme. Un intégrateur explicite doit rester stable pour les deux.

Les solveurs rigides implicites peuvent dépasser les modes de décroissance rapide car leurs régions de stabilité s’étendent beaucoup plus loin dans la moitié gauche du plan complexe. Le compromis est que chaque étape acceptée peut nécessiter une évaluation jacobienne, une factorisation matricielle, une algèbre linéaire itérative ou une solution non linéaire.

La bonne question n’est donc pas de savoir si les méthodes implicites sont toujours meilleures. Il s’agit de savoir si le coût de calcul d’une étape implicite est inférieur au coût des nombreuses petites étapes explicites nécessaires à la stabilité.

Le guide plus large Méthodes d’intégration temporelle pour les solveurs PDE compare les approches explicites, implicites et IMEX. Cet article se concentre spécifiquement sur la rigidité, les régions de stabilité, les jacobiens et les familles de solveurs rigides.

Ce que signifie la rigidité

Il n’y a pas de définition unique de la rigidité qui s’applique parfaitement à chaque système non linéaire, méthode numérique, tolérance et intervalle d’intégration.

Une description pratique utile est :

Un problème est rigoureux lorsqu’une méthode numérique doit prendre des mesures beaucoup plus petites que l’échelle requise pour la précision en raison des restrictions de stabilité.

Cette définition souligne que la rigidité dépend en partie du solveur. Un problème raide pour une méthode Runge-Kutta explicite peut être gérable pour une méthode implicite L-stable.

Cela dépend également de la sortie demandée. Un transitoire chimique rapide peut être scientifiquement important dans une étude et non pertinent dans une autre. S’il doit être résolu, une petite étape est une exigence de précision plutôt qu’une restriction inutile de rigidité.

Un système simple et rigide

Envisagez deux éléments en décomposition :

dy₁/dt = -y₁

dy₂/dt = -1000y₂

Les solutions exactes sont :

y₁(t) = y₁(0)e⁻ᵗ

y₂(t) = y₂(0)e⁻¹⁰⁰⁰ᵗ

Le deuxième composant disparaît beaucoup plus rapidement que le premier. Après son transitoire initial, la solution scientifiquement pertinente peut être contrôlée presque entièrement par y₁.

Le forward Euler appliqué à la deuxième équation n’est stable que lorsque :

|1 - 1000Δt| ≤ 1

Cela donne :

0 ≤ Δt ≤ 0.002

Le solveur doit continuer à utiliser une étape proche de cette échelle même après que y₂ est devenu négligeable. La composante lente pourrait être représentée avec précision avec des étapes beaucoup plus grandes, mais la stabilité explicite est contrôlée par le mode rapide.

Pourquoi le rapport de rigidité est incomplet

Pour un système linéaire :

dy/dt = Ay

Un rapport de rigidité couramment indiqué est :

κ =
max |Re(λᵢ)|
/
min |Re(λᵢ)|

Le rapport peut être informatif lorsque toutes les valeurs propres pertinentes ont des parties réelles négatives et que la matrice se comporte régulièrement. Une valeur importante suggère des taux de décroissance largement séparés.

Cependant, ce n’est pas un test universel. Il n’inclut pas l’intervalle d’intégration, la précision demandée, la condition initiale ou la direction dans laquelle la solution évolue. Il peut également devenir indéfini ou trompeur lorsque le Jacobien contient zéro valeurs propres, des modes de conservation ou un comportement fortement non normal.

Dans un système non linéaire, le Jacobien change le long de la trajectoire :

J(t) =
∂f
/
∂y
evaluated at y(t)

Un système peut donc se déplacer entre des régions rigides et non raides au cours d’une même simulation.

Indicateurs basés sur Jacobien

Les estimations des valeurs propres du Jacobien local sont utiles car les modes stables rapides sont associés à des valeurs propres qui ont de grandes parties réelles négatives.

Un indicateur moins cher utilise parfois la trace mise à l’échelle :

τ(J) =
trace(J) / d
=
div(f) / d

Sur un intervalle d’intégration avec une longueur caractéristique T, une valeur fortement négative de :

T trace(J) / d

Peut indiquer une contraction rapide de la phase-volume et des éventuelles restrictions de pas explicites.

Il s’agit d’un indicateur heuristique, pas une définition nécessaire et suffisante de la rigidité. Une trace fait la moyenne des parties réelles de la valeur propre et peut masquer une direction fortement raide lorsque d’autres directions la compensent.

Les systèmes non normaux peuvent également présenter une croissance transitoire qui n’est pas prédite de manière fiable par les valeurs propres ou les traces seules. Les normes logarithmiques, les diagnostics de solveur direct et les restrictions de taille de pas observées peuvent fournir des informations supplémentaires.

L’arrière-plan mathématique est résumé dans Aperçu des équations rigides. Une analyse plus approfondie de la raison pour laquelle la rigidité résiste à une simple définition basée sur des valeurs propres est fournie dans rigidité des ODE.

Signes pratiques de rigidité

Un modèle peut être raide lorsque :

• Une méthode explicite adaptative réduit à plusieurs reprises son pas malgré une solution fluide.
• Le solveur explicite prend de nombreuses étapes après qu’un transitoire rapide s’est déjà dégradé.
• Un petit changement de taille de pas fait passer le calcul de stable à explosif.
• Le Jacobien contient des modes fortement négatifs avec des échelles très différentes.
• Un solveur implicite effectue le même calcul avec beaucoup moins d’étapes acceptées.
• Les réactions chimiques, la diffusion, la relaxation ou les termes de pénalité introduisent des échelles de temps très courtes.

Ces symptômes ne doivent pas être utilisés pour rejeter d’éventuelles erreurs de codage ou de modélisation. Des signes, des unités, des conditions aux limites ou des Jacobiens incorrects peuvent créer un comportement qui ressemble à la rigidité.

L’équation du test de Dahlquist

La stabilité est généralement étudiée avec :

dy/dt = λy

La solution exacte après une étape est :

y(t + Δt) =
e^(λΔt)y(t)

Une méthode numérique produit :

yⁿ⁺¹ =
R(z)yⁿ

z = λΔt

La fonction R(z) est la fonction de stabilité. La méthode est absolument stable lorsque :

|R(z)| ≤ 1

L’ensemble de toutes ces valeurs forme la région de stabilité absolue.

Une introduction claire à l’équation de Dahlquist et à la stabilité A est disponible dans les Notes d’analyse numérique.

Pourquoi les méthodes explicites échouent

L’attaquant Euler a :

R(z) = 1 + z

Sa région de stabilité est un disque centré sur −1 avec un rayon. Le long de l’axe réel négatif :

−2 ≤ z ≤ 0

Le RK4 classique a une région plus grande, s’étendant à environ :

z ≈ −2.785

Il est toujours délimité. Aucune méthode Runge-Kutta explicite ne peut être stable.

Lorsqu’une valeur propre jacobienne est d’environ −1000, RK4 nécessite une étape inférieure à environ :

Δt ≲ 2.785 / 1000

Une méthode explicite d’ordre supérieur améliore la précision et peut agrandir l’intervalle stable, mais elle ne peut pas étendre la stabilité arbitrairement loin le long de l’axe réel négatif.

A-stabilité

Une méthode A-stable contient le demi-plan entier gauche dans sa région de stabilité :

Re(z) ≤ 0
⇒
|R(z)| ≤ 1

Euler à l’envers est A-stable parce que :

R(z) =
1 / (1 - z)

Pour les modes avec des parties réelles négatives, la méthode reste stable quelle que soit la taille du pas.

Cela ne signifie pas qu’un pas arbitrairement grand est exact. La stabilité A supprime une restriction de stabilité mais ne supprime pas l’erreur de troncature.

L-stabilité

Une méthode L-stable est A-stable et satisfait également :

R(z) → 0
as
z → -∞

Cette propriété amortit fortement les modes rapides non résolus. Euler à l’envers est L-stable.

Crank–Nicolson est A-stable mais pas stable. Sa fonction de stabilité est :

R(z) =
(1 + z/2)
/
(1 - z/2)

Comme z → -∞, le facteur approche −1. Un composant très rigide peut alterner signe au lieu de se décomposer rapidement. Cela peut créer des oscillations non physiques même si la méthode reste limitée.

La L-stabilité est particulièrement précieuse lorsque les transitoires rapides doivent disparaître plutôt que de rester visibles dans la solution calculée.

Le problème Prothero-Robinson

Le problème du test Prothero–Robinson a la forme générale :

dy/dt =
λ[y - g(t)]
+ g'(t)

Sa solution exacte est :

y(t) =
g(t)
+
[y(0) - g(0)]e^(λt)

Choisir :

g(t) = sin(t)

donne :

dy/dt =
λ[y - sin(t)]
+ cos(t)

et :

y(t) =
sin(t)
+
y(0)e^(λt)

Lorsque g(0)=0.

Pour un grand négatif λ, le transitoire exponentiel se désintègre rapidement, tandis que le forçage lisse sin(t) reste. Une méthode explicite peut encore être forcée de résoudre le mode de décroissance rapide pour la stabilité.

Le problème est également utilisé pour étudier la réduction de la commande rigide. Certaines méthodes ne parviennent pas à atteindre leur ordre de convergence nominale lorsque la rigidité, le forçage et les effets de limites ou de contraintes interagissent.

Comparaison Python avec Scipy

L’exemple suivant compare une méthode explicite avec deux solveurs rigides :

import numpy as np

from scipy.integrate import solve_ivp

stiffness = -1000.0
time_interval = (0.0, 10.0)
initial_state = [1.0]

def right_hand_side(time, state):
    value = state[0]

    return [
        stiffness
        * (value - np.sin(time))
        + np.cos(time)
    ]

def jacobian(time, state):
    return np.array([
        [stiffness]
    ])

def exact_solution(time):
    return (
        np.sin(time)
        + np.exp(stiffness * time)
    )

methods = [
    ("RK45", {}),
    ("Radau", {"jac": jacobian}),
    ("BDF", {"jac": jacobian}),
    ("LSODA", {})
]

for method, options in methods:
    solution = solve_ivp(
        fun=right_hand_side,
        t_span=time_interval,
        y0=initial_state,
        method=method,
        rtol=1e-7,
        atol=1e-9,
        **options
    )

    reference = exact_solution(
        solution.t
    )

    maximum_error = np.max(
        np.abs(
            solution.y[0]
            - reference
        )
    )

    print(method)
    print("  Success:", solution.success)
    print("  Accepted points:", solution.t.size)
    print("  Function evaluations:", solution.nfev)
    print("  Jacobian evaluations:", solution.njev)
    print("  LU factorizations:", solution.nlu)
    print("  Maximum error:", maximum_error)

Les décomptes exacts dépendent des tolérances, de la version du logiciel, de l’état initial, des informations de Jacobien et des paramètres du solveur. La comparaison importante est le coût total de calcul au même niveau d’erreur.

Familles de solveurs rigides

Aucun solveur unique n’est idéal pour chaque système rigide. Les considérations importantes comprennent :

• Dimension du système
• La parcimonie et la bande passante de Jacobien
• Coût d’évaluation du côté droit
• Disponibilité des dérivés analytiques
• tolérance demandée
• présence de contraintes ou d’une matrice de masse
• Si le problème change entre les phases rigides et non raides

Formules de différenciation en arrière

Les méthodes BDF sont des formules implicites en plusieurs étapes. Ils se rapprochent de la dérivée à l’aide de plusieurs valeurs de solution précédentes et évaluent l’équation gouvernante au nouveau moment.

Une formule BDF du premier ordre est Euler à l’envers :

(
    yⁿ⁺¹ - yⁿ
) / Δt
=
f(tⁿ⁺¹, yⁿ⁺¹)

BDF2 utilise deux valeurs précédentes :

(
    3yⁿ⁺¹
    - 4yⁿ
    + yⁿ⁻¹
) / (2Δt)
=
f(tⁿ⁺¹, yⁿ⁺¹)

Les implémentations d’ordre variable peuvent utiliser plusieurs commandes BDF et adapter à la fois la taille de la commande et la taille des pas. Les commandes un et deux sont stables. Les formules BDF d’ordre supérieur sacrifient la pleine stabilité A mais restent utiles pour de nombreux systèmes rigides.

Les méthodes BDF sont souvent efficaces pour les grands systèmes lorsque le Jacobien est clairsemé et qu’un solveur linéaire ou un préconditionneur parsemé degrés est disponible.

Les exemples incluent les BDF de Scipy, les méthodes ode15s de Matlab, les cadrans solaires CVODE BDF et les méthodes SCIML telles que QNDF et FBDF.

Méthodes de Rosenbrock

Les méthodes de Rosenbrock sont des méthodes Runge-Kutta linéairement implicites. Ils linéarisent les équations d’étape et résolvent les systèmes de la forme :

(
    I - γΔtJ
)kᵢ
=
known stage terms

Cela évite une solution non linéaire complète à chaque étape. La méthode nécessite toujours des informations jacobiennes ou une approximation utile.

Les méthodes de Rosenbrock sont souvent efficaces pour les petits et moyens raides, en particulier à des tolérances strictes. La même structure matricielle peut être réutilisée en plusieurs étapes.

Les exemples de l’écosystème du SCIML incluent :

• Rosenbrock23
• Rodas4
• Rodas5P

La Sciml Rosenbrock Documentation fournit des recommandations spécifiques à la méthode et des détails de mise en œuvre.

Méthodes Sdirk

Sdirk signifie Runge-Kutta implicite individuellement en diagonale. Chaque étape implicite utilise le même coefficient diagonal dans la matrice Runge-Kutta.

Cela produit des systèmes scéniques avec une structure répétée :

(
    I - γΔtJ
)kᵢ
=
stage residual

Une factorisation matricielle ou un préconditionneur peut donc être réutilisé à travers les étapes lorsque le Jacobien est maintenu fixe.

Les méthodes Sdirk sont des méthodes en une seule étape, ce qui les rend plus faciles à redémarrer après des événements ou des discontinuités que les formules BDF en plusieurs étapes.

Les exemples incluent les méthodes esdirk et kencarp. Certains fournissent des estimations de stabilité en L et d’erreurs intégrées.

Tr-BDF2

TR-BDF2 est une méthode composite en une seule étape. Il combine une étape de règle trapézoïdale avec une étape de différenciation vers l’arrière du second ordre.

Il ne doit pas être regroupé avec des solveurs BDF à plusieurs étapes ordinaires. Sa construction spéciale donne une précision de second ordre et de fortes propriétés de stabilité, notamment la stabilité en L pour le choix des paramètres standard.

TR-BDF2 est souvent une méthode raide à usage général pratique lorsqu’un amortissement modéré et robuste est plus important que l’ordre asymptotique très élevé.

Méthodes Runge-Kutta entièrement implicites

Firk signifie Runge–Kutta entièrement implicite. Les étapes sont couplées et nécessitent normalement la résolution d’un système non linéaire plus grand.

Les familles importantes comprennent :

• Méthodes Radau IIA
• Méthodes de Gauss-Legendre
• Méthodes lobatto

Les méthodes RADAU IIA peuvent fournir un ordre élevé, une stabilité A et une stabilité en L. Ils sont attractifs pour les problèmes de grande précision mais peuvent avoir un coût par étape plus élevé que les méthodes Rosenbrock, Sdirk ou BDF.

Le solveur Radau de Scipy est basé sur une méthode Radau IIa du cinquième ordre.

Méthodes explicites stabilisées

Les méthodes explicites stabilisées utilisent des polynômes de stabilité spécialement construits pour étendre l’intervalle stable le long de l’axe réel négatif.

Les exemples incluent :

• Méthodes RKC
• rock2
• rock4

Pour une méthode avec des étapes s, l’intervalle stable peut croître approximativement comme s². Cela rend l’approche utile pour les grands systèmes dominés par la diffusion où le spectre jacobien se situe principalement à proximité de l’axe réel négatif.

Ces méthodes ne sont pas stables et ne résolvent pas tous les types de rigidité. Ils peuvent ne pas convenir lorsque le système a des valeurs imaginaires substantielles, une rigidité de réaction hautement non linéaire, des contraintes algébriques différentielles ou un couplage difficile.

Le matériau EPFL sur les méthodes explicites stabilisées illustre comment les méthodes basées sur Chebyshev étendent le domaine de stabilité.

Méthodes IMEX

Lorsque la rigidité peut être séparée d’une dynamique non raide, une méthode implicite et explicite peut être plus efficace que de traiter implicitement le système complet :

dy/dt =
f_nonstiff(y)
+
f_stiff(y)

Le terme non raide est intégré explicitement, tandis que le terme rigide est intégré implicitement.

Les exemples incluent l’advection traitée explicitement et la diffusion traitée implicitement, ou le transport traité explicitement et les réactions chimiques traitées implicitement.

Les méthodes IMEX nécessitent une division d’opérateur physiquement et numériquement significative. Leur stabilité dépend de la paire complète explicite-implicite.

Commutation automatique de la rigidité

Certains algorithmes estiment la rigidité lors de l’intégration et passent entre des formules non raides et rigides.

LSODA bascule entre une méthode ADAMS et une méthode BDF. Il est disponible via Odepack, Scipy et plusieurs interfaces linguistiques.

Dans l’écosystème SCIML, une méthode automatique peut combiner un algorithme explicite avec un solveur rigide, par exemple :

AutoTsit5(
    Rosenbrock23()
)

La commutation automatique est utile lorsque la rigidité est inconnue ou n’apparaît que pendant une partie de la trajectoire. La détection n’est pas parfaite, de sorte que des problèmes difficiles peuvent toujours bénéficier de la sélection directe d’un solveur rigide.

Le Guide de solveur d’odes Sciml recommande les méthodes de commutation automatique lorsque la rigidité est incertaine.

Les questions jacobiennes

Pour une étape implicite, le résidu non linéaire peut être :

G(yⁿ⁺¹) =
yⁿ⁺¹
- yⁿ
- Δt f(tⁿ⁺¹, yⁿ⁺¹)

La méthode de Newton nécessite :

[
    I
    - ΔtJ
]δ
=
-G

où :

J =
∂f / ∂y

La solution linéaire peut dominer le temps d’exécution. Un solveur rigide fonctionnera mal s’il traite un grand Jacobien clairsemé comme une matrice dense.

Les performances peuvent s’améliorer grâce à :

• Un Jacobien analytique
• Différenciation automatique
• Un modèle de parcimonie déclaré
• Stockage de la matrice à bandes
• Factorisation directe clairsemée
• Méthodes Krylov
• Préconditionnement basé sur la physique
• Produits vectoriels Jacobian sans matrice

Le didacticiel SCIML sur résolution de grandes équations rigides montre comment la structure jacobienne et la sélection des solveurs linéaires affectent les performances.

Guide de sélection des solveurs

Pourquoi les seuils de dimension fixe sont trompeurs

Il est tentant de recommander les méthodes de Rosenbrock en dessous de 1 000 équations, les méthodes Sdirk jusqu’à 10 000 et les méthodes BDF au-dessus de cela. Ces seuils ne sont pas fiables.

Un problème dense de 500 équations peut être plus coûteux qu’un système clairsemé contenant des centaines de milliers d’équations. Les facteurs décisifs comprennent :

• Nombre d’entrées non nulles de Jacobien
• Remplissage lors de la factorisation
• Qualité de préconditionnement
• Coût de droite
• Tolérance requise
• Fréquence des mises à jour jacobiennes
• Disponibilité de produits sans matrice
• Évolutivité parallèle

Le choix du solveur doit être testé grâce à des comparaisons de précision du travail sur le problème réel.

L’exemple de propagation de la flamme

L’exemple de propagation de la flamme de Cleve Moler utilise l’équation scalaire non linéaire :

dy/dt =
y²(1 - y)

La valeur initiale est très faible. La solution reste proche de zéro pendant une longue période d’induction, passe rapidement, puis s’approche d’une seule.

Dans la comparaison MATLAB publiée, le solveur non rigide ode45 a suivi 3 040 étapes et 20 179 évaluations de fonctions. Le solveur rigide basé à Rosenbrock ode23s a suivi 99 étapes et 412 évaluations de fonctions.

Il s’agit d’environ trente fois moins d’étapes acceptées, et non d’un pour cent du nombre de pas d’origine.

L’exemple et sa représentation analytique Lambert w sont décrits dans équations différentielles raides par Cleve Moler.

Le résultat ne doit pas être interprété comme un avantage universel pour les méthodes implicites. Il montre à quel point l’efficacité du solveur peut changer lorsque la stabilité plutôt que la précision contrôle un intégrateur explicite.

Comparaison équitable des solveurs

Une référence juste doit comparer les méthodes avec une précision équivalente. Enregistrer :

• Erreur ou erreur globale dans une quantité d’intérêt
• Étapes acceptées et rejetées
• Évaluations de droite
• Évaluations jacobiennes
• Résolutions linéaires et factorisations
• Itérations non linéaires
• consommation de mémoire
• Heure totale de l’horloge murale

Un solveur prenant moins d’étapes peut encore être plus lent si chaque étape nécessite une factorisation coûteuse. Une méthode avec plus d’évaluations de fonctions peut gagner lorsque la fonction est bon marché et que l’algèbre linéaire concurrente est coûteuse.

erreurs courantes

• En supposant que chaque solution explicite divergente prouve la rigidité
• Traiter le taux de rigidité comme un diagnostic définitif
• Utiliser la trace jacobienne comme définition universelle de rigidité
• En supposant que chaque méthode implicite est A-stable
• En supposant que chaque méthode A-stable est L-stable
• Appeler TR-BDF2 une méthode BDF en plusieurs étapes ordinaire
• Développer Firk de manière incorrecte en tant que famille Runge-Kutta explicite
• Utilisation d’un Jacobien dense pour un grand système clairsemé
• Ne pas fournir un modèle de jacobien ou de parcimonie lorsque le solveur peut en utiliser un
• Sélection d’un solveur uniquement à partir du nombre d’équations
• Utilisation d’une méthode explicite stabilisée pour la rigidité en dehors de sa conception spectrale
• Comparaison des nombres d’étapes sans comparaison d’erreur et d’exécution
• Prendre des mesures implicites très importantes sans vérifier la précision temporelle

Un flux de travail pratique

1. Vérifiez les équations, les unités, les signes et les conditions initiales.
2. Exécutez une méthode explicite adaptative en tant que diagnostic.
3. Inspectez les étapes acceptées, les étapes rejetées et la taille minimale des étapes.
4. Estimez ou analysez le spectre jacobien lorsque cela est pratique.
5. Vérifiez si les composants rapides doivent être résolus physiquement.
6. Essayez un solveur rigide en L aux mêmes tolérances.
7. Fournir des informations jacobiennes et de la parcimonie.
8. Comparez le travail et la précision, non seulement les étapes acceptées.
9. Répétez avec des tolérances plus serrées.
10. Solveur de documents, tolérances, algèbre linéaire et critères de défaillance.

Guides connexes

• Méthodes d’intégration temporelle pour les solveurs PDE — Comparez les approches explicites, implicites et hybrides.
• Séparation de l’opérateur, fractionnement de Strang et schémas IMEX – Séparez les processus physiques rigides et non raides.
• études sur la qualité et la convergence du maillage – Déterminer si le raffinement spatial augmente la rigidité et modifie la solution.
• Vérification vs. Validation dans les simulations scientifiques — Convergence du solveur de places à l’intérieur Flux de travail de crédibilité.
• simulation & Projets de modélisation — Parcourez les ressources MatForge associées.

Lectures complémentaires

• Équation rigide : définitions, indicateurs et contexte historique
• Higham et Trefethen : rigidité des Odes
• L’IIT note une stabilité absolue et une stabilité A
• Matériau EPFL sur les méthodes explicites stabilisées
• Recommandations de solveur Sciml Ode
• MathWorks : équations différentielles rigides

Conclusion

La rigidité se produit lorsque la stabilité numérique force une méthode pour résoudre les composants rapides de plus près que la sortie souhaitée n’est requise. Les méthodes explicites de Runge-Kutta ont des régions de stabilité limitées, de sorte que les grands modes jacobiens négatifs peuvent imposer des étapes extrêmement petites.

Les méthodes implicites A-stables suppriment cette restriction de stabilité du demi-plan gauche. Les méthodes L-stables vont plus loin en amortissant fortement les modes rapides non résolus. Ces propriétés expliquent pourquoi les méthodes Backward Euler, BDF, Rosenbrock, Sdirk, TR-BDF2 et Radau peuvent surpasser les solveurs explicites sur les systèmes rigides.

Aucun diagnostic ou solveur ne fonctionne pour chaque problème. Les valeurs propres, les rapports de rigidité, la divergence à l’échelle et les restrictions d’étape observées fournissent toutes des informations partielles. La structure du jacobien et du solveur linéaire peut avoir plus d’importance que le nombre nominal d’équations.

Des méthodes explicites stabilisées sont utiles lorsque la rigidité est modérée et dominée par des valeurs propres d’axe négatif-réel. Les méthodes IMEX fonctionnent lorsque des opérateurs rigides et non rigides peuvent être séparés. Les méthodes de commutation automatique aident lorsque le régime change lors de l’intégration.

Le choix final doit être basé sur des tests de précision du travail. Le meilleur solveur rigide est la méthode qui atteint la précision scientifique requise avec le coût de calcul fiable le plus bas.