Méthodes d'estimation et de vérification des erreurs : de l'erreur de troncature à la vérification du code

Reading Time: 9 minutes

Points à retenir clés

L’erreur de troncature est le fondement de toute vérification de code — Comprendre la dérivation de la série Taylor explique pourquoi le MMS fonctionne et pourquoi les taux de convergence correspondent aux prédictions théoriques.
La méthode des solutions manufacturées (MMS) est désormais automatisée – des frameworks comme Moose et FenicsX utilisent Sympy pour dériver symboliquement des termes sources, éliminant les erreurs de calcul manuel qui ont tourmenté les premiers efforts de vérification.
Vous n’avez besoin que de 3 à 4 raffinements de maillage – une fois que l’ordre de précision observé correspond à la prédiction théorique, le raffinement est un effort inutile. La vraie question est : quand la vérification est-elle « assez bonne » ?
Les termes d’erreur d’ordre de tête sont Diagnostic — la structure de l’extension de l’erreur de troncature vous indique exactement quelle approximation dérivée échoue, guidant votre sélection de schéma.

Ce que vous mesurez réellement

Avant de plonger dans les équations, clarifions ce que nous vérifions réellement. Dans la communauté de la simulation, Vérification et validation sont souvent confus, mais ils répondent à des questions fondamentalement différentes :

Vérification demande : « Est-ce que nous résolvons correctement les équations ? » C’est une question de cohérence interne. L’implémentation numérique résout-elle fidèlement le modèle mathématique, que ce modèle soit ou non décrivant la réalité ?
Validation demande : « Sommes-nous en train de résoudre les bonnes équations ? » Il compare les résultats de simulation à des expériences physiques ou à des données de référence.

Si la vue d’ensemble du cadre V&V du Post 182 vous a présenté une vue d’ensemble, cet article s’inscrit dans la machinerie mathématique qui rend la vérification rigoureuse. Considérez-le comme l’élément complémentaire de notre Guide des études de convergence, qui couvre le côté empirique. Ici, nous allons expliquer pourquoi les études de convergence fonctionnent et vous montrer comment les automatiser.

Framework mental : le compas d’erreur de troncature. Imaginez une erreur de troncature comme outil de navigation. L’extension de la série Taylor ne vous donne pas seulement une limite – elle révèle la direction dans laquelle votre schéma est biaisé. Une asymétrie de premier ordre dans les différences entre l’avenir ? C’est votre compas d’aiguille pointant vers la dérivée d’erreur dominante.

Erreur de troncature : fondement mathématique

La source la plus courante d’erreurs numériques dans les méthodes de différences finies, d’éléments finis et de volumes finis est Erreur de troncature – la différence entre l’opérateur différentiel exact et son approximation discrète.

Voici la chose que la plupart des manuels ne mettent pas suffisamment en avant : la dérivation des erreurs de troncature via la série Taylor n’est pas un exercice académique. C’est la base pratique pour chaque étude de raffinement de grille que vous aurez jamais menée. Voyons pourquoi.

Considérez le dérivé de temps d’Euler à l’envers :

$$frac{u(t) – u(t – delta t)}{delta t} = frac{du}{dt} + frac{1}{2} frac{d^2u}{dt2} delta t + o(delta t^2)$$

Le terme d’erreur principal est positif et proportionnel à Δt. Pour l’avant Euler, le signe se retourne. Pour une différence centrale dans l’espace :

$$frac{u(x + delta x) – u(x – delta x)}{2delta x} = frac{du}{dx} + o(delta delta x^2)$$

Remarquez la symétrie – les termes de l’ordre impair s’annulent, donnant une précision deuxième ordre. La symétrie paire/impaire de l’expansion de Taylor est la raison pour laquelle les différences centrales surpassent les schémas avant/arrière.

Une comparaison pratique

Schème	Formule	Erreur de troncature	Ordre
Différence de transfert	`(u(x+Δx) - u(x)) / Δx`	`−½ u'' Δx`	1er
Différence arriérée	`(u(x) - u(x-Δx)) / Δx`	`−½ u'' Δx`	1er
différence centrale	`(u(x+Δx) - u(x-Δx)) / 2Δx`	`+(1/6) u''' Δx²`	deuxième

Le tableau ci-dessus n’est pas seulement une comparaison, c’est un outil de diagnostic. Si votre étude de convergence montre le premier ordre lorsque le deuxième ordre était attendu, le terme d’erreur d’ordre principal vous indique exactement ce qui ne va pas.

L’expression de l’erreur de troncature

Pour une discrétisation générale, l’erreur de troncature prend la forme :

$$tau = C H^r$$

où h est le paramètre de discrétisation (taille de maillage ou pas de temps) et r est le taux de convergence. Cette expression est pas juste une limite asymptotique – c’est la quantité réelle que vous vérifiez de manière empirique. Chaque étude de raffinement de grille est fondamentalement une tentative de confirmer que l’expression d’erreur de troncature théorique correspond à l’erreur mesurée.

Ce lien entre le τ théorique et la mesure empirique des erreurs est ce qui rend la vérification du code rigoureuse plutôt qu’ondulée à la main.

De l’erreur de troncature à la convergence

Maintenant que nous avons l’expression de l’erreur de troncature, connectons-la au grand théorème de l’analyse numérique : le théorème d’équivalence de lax-richtmyer.

Le théorème stipule en anglais simple :

Pour les problèmes linéaires et bien posés, une discrétisation cohérente et stable est convergente.

Déballons chaque trimestre au niveau d’un praticien :

Consistent signifie que l’erreur de troncature passe à zéro sous la forme H → 0. Si votre schéma a τ = o(h²), il est cohérent.

stable signifie que les erreurs ne se développent pas de manière illimitée. La solution numérique reste limitée par rapport aux données initiales. Vous vérifiez la stabilité grâce à l’analyse de Von Neumann, aux méthodes énergétiques ou aux vérifications pratiques de l’état de la LCF.

Convergent signifie que la solution numérique s’approche de la solution exacte comme H → 0.

Le croquis de preuve du théorème (que vous devez connaître intuitivement, pas mémoriser) :

La cohérence garantit que l’opérateur discret se rapproche de l’opérateur continu.
La stabilité limite la propagation des erreurs à chaque pas de temps.
Ensemble, ils garantissent que l’erreur totale – erreur de troncature accumulée sur n = t/Δt étapes – reste bornée et converge vers zéro.

Le théorème de Lax-Richtmyer explique les pourquoi les exigences de stabilité existent. Vous ne pouvez pas simplement réduire H arbitrairement ; Vous devez également vous assurer que votre schéma est stable. Pour une intégration de temps explicite, c’est la condition CFL. Pour les schémas implicites, vous êtes généralement en sécurité, mais la tolérance à l’itération devient le nouveau problème de stabilité.

La méthode des solutions fabriquées

Maintenant, pour le gain pratique : comment vérifiez-vous que votre code résout correctement les équations ? L’approche la plus rigoureuse est la méthode des solutions manufacturées (MMS).

Le rapport Sandia 2000 de Salari et Park a établi le MMS comme norme de l’industrie, et il compte désormais plus de 555 citations. La raison ? Le MMS fonctionne pour les toutes PDE – linéaire, non linéaire, couplé, multiphysique – alors que les repères analytiques n’existent que pour les cas de test simples.

Le flux de travail MMS

La beauté du MMS est qu’il renverse le problème de vérification sur sa tête :

Choisir une solution manufacturée u_manufactured(x, t) — une fonction fluide arbitraire
Substituez à la PDE pour dériver le terme de forçage/source qui fait de u_manufacturé une solution exacte
Dériver les conditions initiales et limites de U_Fabriqué
Exécutez la simulation avec ces entrées modifiées
Comparer la solution calculée par rapport à U_Manufacture

Si la solution numérique correspond à la solution fabriquée dans les limites d’erreur attendues, votre code est vérifié.

Parcourons cela avec du code Python concret. Nous utiliserons deux des frameworks open source les plus largement adoptés : Moose et FenicsX.

Moose MMS : automatisation symbolique

Le module mms de Moose enveloppe SYMPY pour dériver automatiquement les fonctions de forçage. Voici comment vous configurez une étude de convergence spatiale pour une équation de diffusion :

import mms

# Define the PDE and manufactured solution
fs, ss = mms.evaluate("-div(grad(u))", "sin(2*pi*x)*sin(2*pi*y)")

# Print forcing function for MOOSE input file
mms.print_fparser(fs)

# Print exact solution and forcing function as MOOSE hit syntax
mms.print_hit(fs, "force")
mms.print_hit(ss, "exact")

La sortie vous indique exactement ce qu’il faut mettre dans votre fichier d’entrée .i :

8*pi^2*sin(2*x*pi)*sin(2*y*pi)
[force]
  type = ParsedFunction
  expression = '8*pi^2*sin(2*x*pi)*sin(2*y*pi)'
[]
[exact]
  type = ParsedFunction
  expression = 'sin(2*x*pi)*sin(2*pi*y)'
[]

Cette dérivation symbolique est cruciale. Pour une simple équation de diffusion 1D, vous pouvez dériver le terme source à la main. Pour Navier-Stokes ou l’élasticité avec des termes couplés, le calcul symbolique n’est pas facultatif – c’est le seul moyen d’éviter les erreurs.

Le fichier d’entrée Moose ressemble alors à ceci :

[Mesh]
  type = GeneratedMesh
  dim = 2
  nx = 8
  ny = 8
[]

[Kernels]
  [diff]
    type = ADDiffusion
    variable = u
  []
  [force]
    type = BodyForce
    variable = u
    function = force
  []
[]

[BCs]
  [all]
    type = FunctionDirichletBC
    variable = u
    function = exact
    boundary = 'left right top bottom'
[]
[]

[Postprocessors]
  [error]
    type = ElementL2Error
    function = exact
    variable = u
[]
[]

Ensuite, vous automatisez l’étude de convergence :

import mms

# Run 4 levels of refinement for both 1st and 2nd order elements
df1 = mms.run_spatial("diffusion_mms.i", 4, console=False)
df2 = mms.run_spatial("diffusion_mms.i", 4, "Variables/u/order=SECOND")

fig = mms.ConvergencePlot(xlabel="Element Size ($h$)", ylabel="$L_2$ Error")
fig.plot(df1, label="1st Order")
fig.plot(df2, label="2nd Order")
fig.save("convergence_plot.png")

Sur un tracé log-log, la pente de chaque ligne vous donne le taux de convergence observé. Pour les éléments du premier ordre, la pente doit s’approcher de 2. Pour les éléments du second ordre, elle doit s’approcher de 3.

FENICSX : études de convergence de Python Pure

FENICSX (le successeur de Fenics/Dolfin) fournit une interface Python tout aussi puissante. Voici comment calculer les normes d’erreur et les taux de convergence :

from dolfinx import default_scalar_type
from dolfinx.fem import (
    Expression, Function, functionspace,
    assemble_scalar, dirichletbc, form,
    locate_dofs_topological,
)
from dolfinx.fem.petsc import LinearProblem
from dolfinx.mesh import create_unit_square
from ufl import SpatialCoordinate, TestFunction, TrialFunction, div, dx, grad, inner
from mpi4py import MPI

import ufl
import numpy as np

def u_ex(mod):
    return lambda x: mod.cos(2 * mod.pi * x[0]) * mod.cos(2 * mod.pi * x[1])

u_numpy = u_ex(np)
u_ufl = u_ex(ufl)

def solve_poisson(N=10, degree=1):
    mesh = create_unit_square(MPI.COMM_WORLD, N, N)
    x = SpatialCoordinate(mesh)
    f = -div(grad(u_ufl(x)))
    V = functionspace(mesh, ("Lagrange", degree))
    u = TrialFunction(V)
    v = TestFunction(V)
    a = inner(grad(u), grad(v)) * dx
    L = f * v * dx
    u_bc = Function(V)
    u_bc.interpolate(u_numpy)
    facets = locate_entities_boundary(
        mesh, mesh.topology_dim - 1, lambda x: np.full(x.shape[1], True)
    )
    dofs = locate_dofs_topological(V, mesh.topology_dim - 1, facets)
    bcs = [dirichletbc(u_bc, dofs)]
    problem = LinearProblem(
        a, L, bcs=bcs,
        petsc_options={"ksp_type": "preonly", "pc_type": "lu"}
    )
    return problem.solve(), u_ufl(x)

L’idée clé ici est Calcul de norme d’erreur fiable. Lorsque l’erreur est faible, le calcul direct (u_ex - uh)^2 peut souffrir d’erreurs d’arrondi, car vous soustrayez deux nombres presque égaux. Le didacticiel FENICSX recommande d’interpoler d’abord les deux solutions dans un espace de fonction d’ordre supérieur :

def error_L2(uh, u_ex, degree_raise=3):
    degree = uh.function_space.ufl_element().degree
    family = uh.function_space.ufl_element().family_name
    mesh = uh.function_space.mesh
    
    # Create higher-order space for accurate subtraction
    W = functionspace(mesh, (family, degree + degree_raise))
    u_W = Function(W)
    u_W.interpolate(uh)
    
    u_ex_W = Function(W)
    u_ex_W.interpolate(u_ex)
    
    e_W = Function(W)
    e_W.x.array[:] = u_W.x.array - u_ex_W.x.array
    
    error = form(ufl.inner(e_W, e_W) * ufl.dx)
    error_global = mesh.comm.allreduce(assemble_scalar(error), op=MPI.SUM)
    return np.sqrt(error_global)

Ensuite, vous exécutez l’étude de convergence :

Ns = [4, 8, 16, 32, 64]
Es = np.zeros(len(Ns))
hs = np.zeros(len(Ns))

for i, N in enumerate(Ns):
    uh, u_ex = solve_poisson(N, degree=1)
    Es[i] = error_L2(uh, u_numpy)
    hs[i] = 1.0 / Ns[i]
    print(f"h: {hs[i]:.2e} Error: {Es[i]:.2e}")

# Compute observed convergence rates
rates = np.log(Es[1:] / Es[:-1]) / np.log(hs[1:] / hs[:-1])
print(f"Rates: {rates}")

La sortie des éléments du premier ordre indique les taux approchant de 2 :

Rates: [1.61 1.89 1.97 1.99]

Pour les éléments du second ordre, l’approche des taux 3. Il s’agit de la manifestation empirique de l’analyse théorique des erreurs de troncature dont nous avons discuté plus tôt.

Études pratiques de convergence : interpréter les résultats

Vous avez exécuté votre étude de raffinement de grille. Vous avez votre tracé de journal. Les pentes sont proches des valeurs théoriques. Mais qu’est-ce que cela vous dit réellement ?

La formule de commande observée

Pour les tailles de maillage h_i et h_{i-1} avec les erreurs correspondantes e_i et e_{i-1}, l’ordre observé est :

$$p approx frac{ln(e_{i-1} / e_i)}{ln(h_{i-1} / h_i)}$$

Si le rapport de raffinement de la grille est r ≈ 2 (commun dans les études), cela simplifie :

$$p approx log_2(e_{i-1} / e_i)$$

Lorsque votre ordre observé se rapproche de la prédiction théorique, vous avez confirmé que votre code fonctionne dans le régime asymptotique. C’est le jalon de la vérification.

Diplôme : une astuce pratique

Pour les normes d’erreur L2, le calcul direct peut masquer le vrai taux de convergence dû au tour d’horizon. La technique d’augmentation des degrés (interpolant dans un espace d’un degré plus avant la soustraction) est la solution de contournement standard dans FENICSX. Sans cela, vous pouvez signaler une convergence de second ordre lorsque le code converge réellement au troisième ordre.

Quand arrêter de vérifier

Voici la réalité pratique : la vérification coûte cher. Chaque raffinement de maillage double (ou quadruple pour la 2D) votre coût de calcul. Vous ne pouvez pas affiner à l’infini. Alors tu t’arrêtes quand ?

Le domaine a convergé vers trois heuristiques pratiques :

1. L’ordre observé correspond à la prédiction théorique. Lorsque vous avez plus de 3 améliorations de maillage et que l’ordre observé converge vers la valeur théorique dans la tolérance (généralement ±0,1), votre code est vérifié dans le régime asymptotique. Aucun raffinement supplémentaire n’est nécessaire à des fins de vérification.

2. L’indice de convergence de la grille (GCI) tombe en dessous de ~1 %. Le GCI quantifie la bande d’incertitude du raffinement du maillage. Lorsque GCI < 1 % de votre quantité d’intérêt, l’erreur numérique est négligeable par rapport à l’incertitude de modélisation.

3. L’erreur numérique tombe en dessous de l’incertitude physique. Si vos propriétés matérielles présentent une incertitude de 5 %, le raffinage jusqu’à ce que l’erreur numérique soit de 0,001 % est inutile. L’estimation d’erreur τ = C H^r vous indique où se trouve le plancher pratique.

Une liste de contrôle de décision pratique

[ ] Exécutez 3 à 4 raffinements de maillage avec un rapport de raffinement r ≈ 2
[ ] Calculer l’ordre observé P à partir de la pente du log-log
[ ] Vérifier que p correspond à la prédiction théorique dans ±0,1
[ ] Calculez GCI pour le maillage le plus fin
[ ] Comparez GCI à l’incertitude de modélisation (généralement 1 à 5 %)
[ ] Arrêtez quand : P est confirmé et GCI < ; 1 % du QI

Si vous utilisez MMS pour un code de production, cette liste de contrôle est suffisante. Pour les applications critiques pour la sécurité (comme la simulation de réacteurs nucléaires dans le manuel MPACT), des vérifications supplémentaires pour chaque terme dans l’équation régissant sont nécessaires.

Résumé

Le parcours de l’erreur de troncature à la vérification du code suit un chemin clair :

Série Taylor révèle l’expression d’erreur de troncature théorique τ = c h^r
Lax-Richtmyer Theorem garantit la convergence lorsque le schéma est cohérent et stable
MMS fournit le flux de travail pratique pour vérifier que votre code atteint le taux de convergence théorique
Études de raffinement de la grille confirment que les taux observés correspondent aux prédictions théoriques

L’idée clé qui lie cela : l’erreur de troncature n’est pas seulement un concept théorique, c’est la quantité mesurable que confirment vos études de raffinement de grille. Lorsque l’ordre observé fait correspondre la prédiction théorique à plus de 3 améliorations, votre code est vérifié.

Pour les prochaines étapes, consultez notre aperçu du cadre V&V pour le contexte généraliste et notre guide des études de convergence pour les recommandations pratiques de qualité de maillage qui complètent ces flux de travail de vérification.

Références

Salari, K. & Park, K. (2000). « Vérification de code par la méthode des solutions fabriquées. » Rapport Sandia National Laboratories SAND2000-0949. source
Oberkampf, W. & Roy, C. (2010). « Vérification et validation en informatique scientifique. » Cambridge University Press.
Langtangen, H.P. « Analyse des erreurs de troncature pour les méthodes de différences finies. » Méthodes numériques pour la documentation PDE.
Documentation MMS Moose. Framework Moose
Tutoriel de convergence Fenicsx. Dolfinx Dolfinx
Kindo, T. « Vérifier les simulations avec la méthode des solutions fabriquées. » Blog COMSOL. Source

Méthodes d’estimation et de vérification des erreurs : de l’erreur de troncature à la vérification du code