Méthode des différences finies

En analyse numérique, la méthode des différences finies est une technique courante de recherche de solutions approchées d'équations aux dérivées partielles qui consiste à résoudre un système de relations (schéma numérique) liant les valeurs des fonctions inconnues en certains points suffisamment proches les uns des autres.

Cette méthode apparaît comme étant la plus simple à mettre en œuvre car elle procède en deux étapes : d'une part la discrétisation par différences finies des opérateurs de dérivation/différentiation, d'autre part la convergence du schéma numérique ainsi obtenu lorsque la distance entre les points diminue.

Approximation des opérateurs par formules de Taylor

Une discrétisation des opérateurs différentiels (dérivées premières, secondes, etc., partielles ou non) peut être obtenue par les formules de Taylor.

La formulation de Taylor-Young est préférable dans son utilisation simple, la formulation de Taylor avec reste intégral de Laplace permet de mesurer les erreurs (cf. plus bas).

Exemples d'approximations d'opérateurs

Approximations centrées

Si u est définie sur un voisinage d'un point x et trois fois dérivable en x, la formule de Taylor-Young conduit aux deux relations :

${\displaystyle u(x+h)=u(x)+\sum _{n=1}^{3}{{\frac {h^{n}}{n!}}u^{(n)}(x)}+h^{3}\varepsilon _{1}(h)}$

${\displaystyle u(x-h)=u(x)+\sum _{n=1}^{3}{{\frac {(-h)^{n}}{n!}}u^{(n)}(x)}+h^{3}\varepsilon _{2}(h)}$

où les deux fonctions ε_i(h) convergent vers 0 avec h. Par conséquent

${\displaystyle {\frac {u(x+h)-u(x)}{h}}=u'(x)+{\frac {h}{2}}u''(x)+{\frac {h^{2}}{3!}}u^{(3)}(x)+h^{2}\varepsilon _{1}(h)}$

${\displaystyle {\frac {u(x-h)-u(x)}{h}}=-u'(x)+{\frac {h}{2}}u''(x)-{\frac {h^{2}}{3!}}u^{(3)}(x)+h^{2}\varepsilon _{2}(h)}$

correspondent à deux approximations de u'(x) du 1^er ordre en h. La variable scalaire h est appelée pas du schéma aux différences.

En soustrayant les développements précédents, ce qui revient à faire la moyenne des deux différences finies antérieure et postérieure à u(x), on obtient

${\displaystyle {\frac {u(x+h)-u(x-h)}{2h}}=u'(x)+{\frac {h^{2}}{6}}u^{(3)}(x)+h^{2}\varepsilon _{3}(h)}$

qui est une approximation de u'(x) du 2^e ordre en h.

Approximations décentrées

Décentrage aval

Si u est définie sur un voisinage à droite d'un point x et trois fois dérivable en x, la formule de Taylor-Young conduit à la relation :

${\displaystyle u(x+h)=u(x)+\sum _{n=1}^{3}{{\frac {h^{n}}{n!}}u^{(n)}(x)}+h^{3}\varepsilon _{1}(h)}$

où la fonction ε₁(h) converge vers 0 avec h. Par conséquent

${\displaystyle {\frac {u(x+h)-u(x)}{h}}=u'(x)+{\frac {h}{2}}u''(x)+{\frac {h^{2}}{3!}}u^{(3)}(x)+h^{2}\varepsilon _{1}(h)}$

correspond à une approximation de u'(x) du 1^er ordre en h.

En répétant l'opération pour un décentrage aval, en écrivant que :

${\displaystyle u''(x)={\frac {{\frac {u(x+2h)-u(x+h)}{h}}-{\frac {u(x+h)-u(x)}{h}}}{h}}+\varepsilon _{2}(h)}$

d'où

${\displaystyle {\frac {u(x+2h)-2u(x+h)+u(x)}{h^{2}}}=u''(x)+hu^{(3)}(x)+h\varepsilon _{3}(h)}$

qui est une approximation de u''(x) du 2^e ordre en h.

Formules étendues aux ordres successifs

En étendant la taille du stencil, il est possible de déterminer des différences finies d'ordres supérieurs par des méthodes similaires (augmentation de l’ordre dans la formule de Taylor et détermination d'une combinaison linéaire adaptée pour annuler les termes superflus).

Si u est définie sur un voisinage d'un point x et quatre fois dérivable en x, par extension de la formule de Taylor, on peut montrer que les schémas à cinq points

${\displaystyle u'(x)={\frac {u(x-2h)-8u(x-h)+8u(x+h)-u(x+2h)}{12h}}+\varepsilon _{4}(h)}$

${\displaystyle u''(x)={\frac {-u(x-2h)+16u(x-h)-30u(x)+16u(x+h)-u(x+2h)}{12h^{2}}}+\varepsilon _{5}(h)}$

sont des approximations des dérivées première et seconde d'ordre 4.

Extension aux fonctions à plusieurs variables

Décentrage amont

En un point (x, y) tel que u soit quatre fois dérivable sur un carré [0, x + a] × [0, y + a], on peut écrire

${\displaystyle {\frac {(u(x+2h,0)-2u(x+h,0)+u(x,0))+(u(0,y+2h)-2u(0,y+h)+u(0,y))}{h^{2}}}=\Delta u(x,y)+{\frac {h^{2}}{12}}(\partial _{x}^{4}u(x,y)+\partial _{y}^{4}u(x,y))+h^{2}\varepsilon _{6}(h)}$

qui est une approximation du Laplacien Δu(x , y) du 2^e ordre en h (cf. équation de Laplace et équation de Poisson).

Stencil centré

En un point (x, y) tel que u soit quatre fois dérivable sur un carré [x - a, x + a] × [y - a, y + a], on peut écrire

${\displaystyle {\frac {u(x+h,y)+u(x-h,y)+u(x,y+h)+u(x,y-h)-4u(x,y)}{h^{2}}}=\Delta u(x,y)+{\frac {h^{2}}{12}}(\partial _{x}^{4}u(x,y)+\partial _{y}^{4}u(x,y))+h^{2}\varepsilon _{7}(h)}$

qui est une approximation du Laplacien Δu(x , y) du 2^e ordre en h (cf. équation de Laplace et équation de Poisson).

La notion d'« ordre » mentionnée ci-dessus correspond à un concept de convergence locale de l’opérateur discrétisé. La convergence globale de la solution discrète est un concept très différent, même s’il existe un lien de parenté entre les deux.

Maillage

Pour la méthode des différences finies, un maillage est un ensemble de points isolés (appelés nœuds) situés dans le domaine de définition des fonctions assujetties aux équations aux dérivées partielles, une grille sur les seuls nœuds de laquelle sont définies les inconnues correspondant aux valeurs approximatives de ces fonctions^[1].

Le maillage comprend également des nœuds situés sur la frontière du domaine (ou au moins « proches » de cette frontière) afin de pouvoir imposer les conditions aux limites et/ou la condition initiale avec une précision suffisante.

A priori, la qualité première d’un maillage est de couvrir au mieux le domaine dans lequel il se développe, de limiter la distance entre chaque nœud et son plus proche voisin. Cependant, le maillage doit également permettre d’exprimer la formulation discrète des opérateurs de différentiation : pour cette raison, les nœuds du maillage sont le plus souvent situés sur une grille dont les directions principales sont les axes des variables.

On appelle pas du maillage la distance entre deux nœuds voisins situés sur une droite parallèle à l’un des axes. Dans ce sens, le pas est une notion à la fois locale et directionnelle. On parlera de pas global pour désigner le plus grand pas local, une notion qui reste directionnelle.

Bien qu’un pas constant soit le plus souvent retenu (sans qu’il pose de problème théorique pour la résolution), il est parfois judicieux d’introduire un pas variable qui sera choisi plus fin dans les zones où la solution exacte subit de plus fortes variations : cette astuce permet de réduire le nombre d’inconnues sans porter atteinte à la précision des résultats. Par contre, la formulation est un peu plus complexe car la discrétisation des opérateurs différentiels doit en tenir compte.

Exemples de maillages

Pour une équation différentielle concernant une fonction d’une variable dont le domaine (dans ${\displaystyle \mathbb {R} }$) est l’intervalle [0 ; 1], un maillage à pas constant est caractérisé par les M + 1 nœuds x_i = i h, 0 ≤ i ≤ M avec le pas h = 1/M. Ce maillage comprend les deux points frontière x₀ et x_M sur lesquels sont imposées d’éventuelles conditions aux limites.

Considérons une équation aux dérivées partielles concernant une fonction de deux variables (domaine ${\displaystyle \Omega \subset \mathbb {R} ^{2}}$) :

• Si Ω est un rectangle [0 ; 1]×[0 ; 1] (dont les côtés sont parallèles aux axes), un maillage issu d’une grille (x_i , y_j) = (i h, j k), 0 ≤ i ≤ M, 0 ≤ j ≤ N avec les pas h = 1/M et k = 1/N est une simple généralisation du cas précédent.
• Si Ω est un disque centré à l’origine et de rayon 1, on considère le maillage constitué par les nœuds d’une grille qui sont situés dans le disque, soit (x_i , y_j) ∈ Ω où (x_i , y_j) = (i h, j k) avec le pas h= 1/M. Pour imposer d’éventuelles conditions aux limites (par exemple celles de Dirichlet qui fixent la valeur de la fonction sur ∂Ω), les rares nœuds se situant exactement sur la frontière sont trop peu représentatifs. Il convient alors d’étendre la propriété « être sur la frontière » à d’autres nœuds qui en sont proches, en englobant par exemple tous les nœuds du maillage qui n’ont pas quatre voisins directs. Les valeurs aux limites à fixer en ces nouveaux nœuds frontières peuvent être définies de diverses manières :
  • En prenant la valeur du problème exact qui est imposée au point de ∂Ω le plus proche : dans ce cas, les irrégularités géographiques des nœuds frontières du maillage (observées lorsque h diminue) engendrent des perturbations de la solution discrète, au mieux des anomalies locales n’ayant aucun lien avec la solution exacte.
  • En considérant que les valeurs des nouveaux nœuds frontières sont des inconnues, mais en ajoutant des relations différentielles discrétisées supplémentaires reliant « naturellement » ces inconnues aux valeurs des nœuds voisins et à celles de certains points de ∂Ω. Si l’approche est un peu plus complexe dans sa mise en œuvre, elle réduit significativement le défaut de la précédente.

Schéma numérique

Un schéma numérique peut être défini comme la formulation algébrique d’un problème discret conçu à l’aide de la méthode des différences finies. La démarche comprend les étapes suivantes :

• Choisir les opérateurs discrets qui sont des approximations des opérateurs différentiels de la formulation exacte.
• Générer un maillage du domaine de définition en étant attentif aux nœuds frontières et à la manière de traduire les conditions aux limites.
• En se fondant sur les expressions issues des opérateurs discrets, établir les relations liant les valeurs des fonctions aux nœuds du maillage (les inconnues).
• S’assurer que l’ensemble des inconnues et des relations qui les relient constitue un problème numérique qui ne soit pas sur- ou sous-déterminé. Cette vérification est une condition minimale pour espérer trouver une solution, mais elle ne donne aucune garantie sur la convergence globale.

Une fois que le schéma numérique est établi et que le problème discret est formulé, il s’agit non seulement de le résoudre, mais encore de s’assurer que la solution discrète converge vers la solution exacte lorsque les pas du maillage tendent vers 0.

Pour certains schémas dits explicites, il est possible d’ordonner les inconnues de telle sorte que chacune d’elles puisse être déterminée récursivement à partir des précédentes qui sont supposées être déjà calculées (matrice triangulaire). Pour les schémas implicites, il est parfois possible d’éviter de résoudre l’ensemble du système de toutes les équations. C’est en particulier le cas pour un système évolutif dont l’état, caractérisé par des variables spatiales, est défini par des conditions initiales (t = 0), puis évolue progressivement au cours du temps : le schéma numérique reste explicite dans la variable temporelle et son caractère implicite ne concerne que les variables spatiales.

Dans tous les cas, chaque équation du schéma numérique ne concerne qu’un petit nombre d’inconnues^[2]. Dans un environnement linéaire, cette propriété conduit à formuler le problème discret à l’aide de matrices creuses et à en tirer profit pour le résoudre en utilisant des méthodes appropriées. Cet avantage est indéniable lorsque la taille du maillage dépasse le cadre d’une étude didactique.

La résolution des schémas numériques s’appuie en général sur des méthodes algébriques classiques. Cependant d’autres formulations équivalentes peuvent faire appel à des méthodes d'optimisation.

Exemple de schéma numérique

Considérons le problème suivant :

${\displaystyle u'(x)-\tau u(x)=0\,\forall x\in [0,1],\ u(0)=u_{0}.}$

Ce problème reste académique dans la mesure où la solution exacte est connue :

${\displaystyle u(x)=u_{0}\,e^{\tau x}.}$

Avec le schéma d'Euler explicite d'ordre 1 appliqué sur un maillage régulier de pas h = 1/M, les inconnues u_n reflétant u(n h) sont liées par les relations

${\displaystyle {\frac {u_{n+1}-u_{n}}{h}}-\tau u_{n}=0,\;\forall \,0\leqslant n<M.}$

Ce schéma conduit à la relation de récurrence

${\displaystyle u_{n+1}=(1+h\tau )u_{n}}$

dont la solution explicite est

${\displaystyle u_{n}=\left(1+{\frac {\tau }{M}}\right)^{n}u_{0}.}$

Une autre formulation obtenue à l’aide du schéma d’ordre 2 (sauf au nœud n = 1 pour lequel on conserve le schéma d'ordre 1) donne

${\displaystyle {\frac {u_{1}-u_{0}}{h}}-\tau u_{0}=0,}$

${\displaystyle {\frac {u_{n+1}-u_{n-1}}{2h}}-\tau u_{n}=0,\;\forall \,0<n<M.}$

Comme le premier, ce second schéma est explicite.

Il est très facile de déterminer numériquement les solutions de ces deux schémas pour les comparer à la solution exacte. Il semble légitime de s’attendre à de meilleurs résultats avec le second schéma puisque son ordre est supérieur à celui du premier (il est possible de montrer que les deux schémas numériques sont uniformément convergents) :

• cette intuition se vérifie lorsque τ > 0 ;
• par contre, elle s’avère erronée lorsque τ est assez négatif (par exemple τ < –3). En effet, le schéma d’ordre 1 qui est appliqué au premier nœud génère dans ce cas un écart initial u₁ - u(h) qui va engendrer des oscillations sur les écarts constatés aux nœuds suivants, oscillations qui s’amplifient jusqu’aux derniers nœuds : le schéma d'Euler engendre sans conteste une solution nettement meilleure.

Cette comparaison montre clairement qu’une bonne représentation des opérateurs différentiels n’est pas une condition suffisante pour obtenir un bon schéma numérique.

Convergence

La convergence d’un schéma numérique est une propriété théorique globale assurant que l’écart (au sens d’une norme) entre la solution approchée et la solution exacte tend vers 0 lorsque le pas de discrétisation tend vers 0 (ou lorsque chacun des pas globaux associés aux différentes directions tend vers 0).

La solution approchée d’un schéma numérique reste peu crédible tant que sa convergence n’a pas été montrée. Cette preuve est sans doute le point le plus délicat de la méthode des différences fines, en tout cas celui qui nécessite l’usage d’outils analytiques.

Il ne suffit pas de vérifier à l’aide d’exemples numériques concrets que le comportement de la solution discrète est conforme aux attentes pour s’assurer de la convergence. Par contre, de tels exemples peuvent aider à prouver le contraire.

Conceptuellement, les écarts entre solution approchée et solution exacte se manifestent par une combinaison de deux phénomènes :

• L’écart induit localement par les approximations inhérentes à l’opérateur discrétisé. C’est la notion de consistance ou de robustesse du schéma. Par exemple, un ordre d’approximation élevé (obtenu par le théorème de Taylor) n’est pas d’une grande valeur lorsque la solution exacte n’a pas la régularité requise.
• Dans le cas d’un schéma explicite, la propagation des écarts qui, au cours des étapes du calcul, se combinent aux écarts précédents et peuvent progressivement s’amplifier. C’est la notion de stabilité du schéma qui se conçoit également lorsqu’il est implicite.

Ces concepts ne considèrent pas les erreurs d’arrondi qui peuvent encore compliquer les choses, comme en témoigne la figure ci-dessous, obtenue avec un exemple concret :

Instabilité numérique d’un schéma appliqué dans un cas particulier où il est « théoriquement » convergent

La norme pour laquelle la convergence est étudiée doit rester indépendante des pas de discrétisation. Cependant, il est fréquent d’utiliser des normes apparentées à celles des espaces L^p. Pour une fonction d’une variable :

• convergence simple : convergence en tout point de l’espace de définition ;
• convergence uniforme (dans L^∞) : ${\displaystyle \|u_{h}-u\|_{\infty }=h\max _{i}|u_{i}-u(ih)|}$ ;
• convergence dans L² : ${\displaystyle \|u_{h}-u\|_{2}=h{\sqrt {\sum _{i}|u_{i}-u(ih)|^{2}}}}$ ;
• convergence dans L¹ : ${\displaystyle \|u_{h}-u\|_{1}=h\sum _{i}|u_{i}-u(ih)|}$.

Dans le cadre d’un problème évolutif avec condition initiale, le théorème de Lax précise rigoureusement les notions de consistance et de stabilité, la seconde étant une condition nécessaire et suffisante pour assurer la convergence.

Exemple de convergence

Dans le dernier exemple présenté ci-dessus pour lequel on connaît à la fois la solution exacte ${\displaystyle \displaystyle u(x)=u_{0}e^{\tau x}}$ et la solution approchée (schéma d'Euler) ${\displaystyle \displaystyle u_{n}=u_{0}(1+\tau h)^{n}}$, le rapport satisfait

${\displaystyle \ln \left(\left|{\frac {u_{n}}{u(nh)}}\right|\right)=n|\ln(1+\tau h)-\tau h|=n\,O(h^{2})=O(h)}$

qui tend vers 0 lorsque ${\displaystyle \displaystyle h{\ddot {a}}}$ tend vers 0, ceci uniformément pour ${\displaystyle n\leqslant 1/h.}$

Ainsi ${\displaystyle \displaystyle u(nh)-u_{n}}$ tend uniformément vers 0, ce qui prouve la convergence de ce schéma d'Euler dans la norme ${\displaystyle \|.\|_{\infty }.}$

Notes et références

1. On observe ici une différence avec la méthode des éléments finis où les fonctions inconnues sont définies en tout point des éléments du maillage (constitué de triangles, de rectangles, etc.). Cependant, puisqu’il est possible d’interpoler les résultats obtenus par différences finies, cette dissemblance n’est pas fondamentale.
2. Même lorsque la formulation du problème d’origine fait intervenir des intégrales, il est souvent avantageux d’introduire des fonctions supplémentaires pour les représenter.

Voir aussi

Articles connexes

• Méthode des volumes finis
• Théorème de Lax-Wendroff

Liens externes

• Liens vers des outils de discrétisation (téléchargeables)
• Formulations mathématiques et résolution numérique en mécanique
• Exemple d'utilisation dans la simulation 1D

• Portail de l'analyse
• Portail de la physique
• Portail du génie mécanique