Méthode des différences finies (encyclopédie mathématiques)

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.

En apparence, 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.

Toutefois, il convient de rester vigilant et critique sur les résultats obtenus tant que la seconde étape n’a pas été montrée en toute rigueur.

Sommaire

[modifier] Approximation des opérateurs par formules de Taylor

Article détaillé : différence finie.

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).

[modifier] Exemples d'approximations d'opérateurs

En un point $x$ et pour une valeur $h$ du pas de discrétisation tels que $u$ soit trois fois dérivable sur l’intervalle $[x - h, x + h]$ , la formule de Taylor-Young conduit aux deux relations :

$u( x + h ) = u( x ) + \sum_{ n = 1 }^3{ \frac{ h^n }{ n! } u^{(n)}( x ) } + h^3 \varepsilon_1( x, h )$

$u( x - h ) = u( x ) + \sum_{ n = 1 }^3{ \frac{ (-h)^n }{ n! } u^{(n)}( x ) } + h^3 \varepsilon_2( x, h )$

où les deux fonctions $ε i (x, h)$ convergent vers 0 avec $h$ . Par conséquent

$\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( x, h )$

$\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( x, h )$

correspondent à deux approximations de $u'(x)$ du 1^er ordre en $h$ .

En soustrayant les développements précédents, on obtient

$\frac{ u( x + h ) - u( x - h ) } { 2 h } = u'(x) + \frac{ h^2 }{ 6 } u^{(3)}( x ) + h^2 \varepsilon_3( x,h )$

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

Si $u$ est 4 fois dérivable et que les développements sont étendus à cet ordre, leur somme donne

$\frac{ u( x + h ) + u( x - h ) - 2 u( x ) } { h^2 } = u''(x) + \frac{ h^2 }{ 12 } u^{(4)}( x ) + h^2 \varepsilon_4( x,h )$

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

L’approximation précédente s’étend à des fonctions de plusieurs variables. En un point $(x, y)$ et pour une valeur $h$ du pas de discrétisation (le même dans les deux dimensions) tels que $u (x, y)$ soit 4 fois dérivable sur le rectangle $\scriptstyle [x-h, x+h]\times[y-h, y+h]$ , alors

$\frac{ u(x+h,y) + u(x-h,y) + u(x,y+h) + u(x,y-h) - 4 u(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_5(x,y,h )$

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

Remarque : 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.

[modifier] 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.

[modifier] Exemples de maillages

Pour une équation différentielle concernant une fonction d’une variable dont le domaine (dans $\R$ ) est l’intervalle $\displaystyle [0,1]$ , un maillage à pas constant est caractérisé par les $\displaystyle M+1$ nœuds $x_i = i h, 0 \leqslant i \leqslant M$ avec le pas $\displaystyle h = 1/M$ . Ce maillage comprend les deux points frontière $\displaystyle x_0$ et $\displaystyle 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 $\Omega \subset \R^2$ ) :

Si $\displaystyle \Omega$ est un rectangle $\displaystyle [0,1] \times [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 \leqslant i \leqslant M, 0 \leqslant j \leqslant N$ avec les pas $\displaystyle h = 1/M$ et $\displaystyle 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 où avec le pas . 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 $\partial \Omega$ le plus proche : dans ce cas, les irrégularités géographiques des nœuds frontières du maillage (observées lorsque $\displaystyle 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 $\partial \Omega$ . 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.

[modifier] Abaissement du degré de dérivation

Pour des raisons d'écriture algébrique et surtout pour l'étude a priori sur la convergence et la stabilité, il est parfois utile de reformuler le problème d’origine en un problème équivalent dont les ordres de dérivation sont inférieurs, ceci en introduisant des fonctions intermédiaires qui sont des dérivées ou dérivées partielles des fonctions du problème initial.

[modifier] Exemple

Le problème

$\displaystyle u''( x ) + \omega^2 u( x ) = 0 \,$ sur $\displaystyle [0,1]$ , avec $u( 0 ) = u_0, \, u'( 0 ) = u_1$

peut être reformulé de la manière suivante :

$U( x ) = \begin{pmatrix} \omega u( x ) \\ u'(x) \end{pmatrix}$

satisfaisant l'équation

$U'( x ) = A \, U(x), \, U_0 = \begin{pmatrix} \omega u_0 \\ u_1 \end{pmatrix}$

où

$A = \begin{pmatrix} 0 & \omega \\ - \omega & 0 \end{pmatrix}$

Le degré de dérivation est réduit au détriment du nombre de fonctions inconnues.

[modifier] 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’elle 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.

[modifier] Exemple de schéma numérique

Considérons le problème suivant :

$\displaystyle u'( x ) - \tau u( x ) = 0 \,$ sur $\displaystyle [0,1]$ avec $\displaystyle u( 0 ) = u_0.$

Ce problème reste académique dans la mesure où la solution exacte est connue, soit $\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 $\displaystyle h = 1/M$ , les inconnues $\displaystyle u_n$ reflétant $\displaystyle u( n h )$ sont liées par les relations

$\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

$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

$\frac{ u_1 - u_0 }{ h } - \tau u_0 = 0,$

$\frac{ u_{ n + 1 } - u_{ n - 1 } }{ 2 h } - \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^[3] :

Cette intuition se vérifie lorsque $\displaystyle \tau > 0$ .
Par contre, elle s’avère erronée lorsque $\displaystyle \tau$ est assez négatif (par exemple $\displaystyle \tau < -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 $\displaystyle u_1 - 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.

[modifier] 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 tendent 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 approximation 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.
Article détaillé : Consistance (mathématiques).
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.
Article détaillé : Stabilité d'un schéma numérique.

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^∞) : $\|u_h - u\|_\infty = h \max_i | u_i - u(ih) |.$
Convergence dans L² : $\|u_h - u\|_2 = h \sqrt{ \sum_i | u_i - u(ih) | ^2}.$
Convergence dans L¹ : $\|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.

[modifier] Exemple de convergence

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

$\ln (| \frac{u_n}{u( n h )}|) = n | \ln(1 + \tau h) - \tau h | = n O(h^2) = O(h)$

qui tend vers 0 lorsque $\displaystyle h$ tend vers 0, ceci uniformément en $n \leqslant 1/h.$

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

[modifier] Notes

↑ 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.
↑ 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.
↑ Il est possible de montrer que les deux schémas numériques sont uniformément convergents.

[modifier] Voir aussi

[modifier] Articles connexes

Différence finie
Consistance d'un schéma numérique
Stabilité d'un schéma numérique
Méthode des éléments finis
Méthode des volumes finis
Théorème de Lax
Théorème de Lax-Wendroff

Méthode des différences finies