Les schémas aux différences finies

Comment résoudre une EDP

Nous avons abordé dans différentes pages de tangenteX plusieurs équations aux dérivées partielles (EDP) : l'équation de diffusion thermique, l'équation de Laplace, l'équation de Poisson. A chaque fois, nous avons cherché un moyen de résoudre numériquement ces différents types d'équations et à chaque fois, nous avons cherché le schéma numérique adéquate.

Il existe trois familles de schémas pour résoudre une EDP :

les schémas aux différences finies (les FDM) : ils conviennent pour les géométries simples. Ils sont adaptés aux problèmes de Dirichlet mais ils prennent difficilement en compte les conditions aux limites de type Von Neumann (conditions sur les dérivées). On les utilise beaucoup pour l'équation de la chaleur et toutes les EDP paraboliques.
les schémas aux éléments finis (les FEM) : ils conviennent particulièrement aux géométries complexes et sont très efficaces. Ils sont classiquement utilisés dans les logiciels de calcul multiphysiques professionnels. Pour nos besoins, ils sont complexes à mettre en oeuvre et sont couteux en temps de calcul et en espace mémoire. Python offre plusieurs librairies d'implémentation des schémas FEM : FiPy, SfePy et FEniCS. Elles sont d'un usage relativement complexe.
enfin, les schémas aux volumes finis. Ils sont basés sur les lois de conservation, ce qui explique qu'on les utilise beaucoup en mécanique des fluides. Ils permettent de traiter des géométries très complexes en 3D et prennent très bien en compte les conditions de Von Neumann.

Tous ces schémas sont basés sur le même principe : on discrétise la courbe, la surface ou le volume supposé continu sur lequel on veut intégrer l'EDP, c'est à dire qu'on le découpe en éléments, segments, surfaces ou volumes, petits mais pas infiniment petits. C'est l'étape de discrétisation.

Cela signifie qu'au lieu d'intégrer l'EDP comme une fonction continue on l'intègre comme une fonction discrète. Bien sur, pour ce faire, il va falloir trouver les méthodes qui nous permettent de mener une telle intégration. C'est l'objet des schémas aux différences, éléments ou volumes finis.

Dans cette page, nous nous limiterons aux différences finies, ce schéma étant le plus simple à comprendre et à mettre en oeuvre.

Quelques rappels

La dérivée

Il peut être utile de rappeler la définition de la dérivée première (d'ordre 1) d'une fonction d'une variable continue et dérivable sur un intervalle quelconque . C'est la limite, si elle existe de :
\( \lim_{h \rightarrow 0} \dfrac{f(x + h) - f(x)}{h} \)
On la note f'(x) ou plus généralement chez les physiciens \( \dfrac{df}{dx} \) et donc :
\( \dfrac{df}{dx} = \lim_{h \rightarrow 0} \dfrac{f(x + h) - f(x)}{h} \)

Vous notez qu'on établit donc une relation entre la dérivée f' d'une fonction et cette fonction.

Dans le cas d'une fonction g(x,y) à deux variables (par exemple), la définition de la dérivée par rapport à une des variables de cette fonction, la dérivée partielle, reste de la même forme. On a:
\( \dfrac{\partial g(x,y)}{\partial x} = \lim_{\Delta x \rightarrow 0} \dfrac{g(x + \Delta x,y ) - g(x,y)}{\Delta x} \)
Là aussi, on établit une relation entre la dérivée partielle de la fonction et la fonction elle-même.

Le théorème de Taylor-Young

Un théorème très utile de l'analyse va servir de base à la construction des méthodes des différences finies. Il s'agit du théorème de Taylor-Young, qui nous dit:
Soit I un intervalle ouvert non vide de \( \mathbb{R} \) et un point a quelconque dans cet intervalle. Soit une fonction f de I dans \( \mathbb{R} \) et un entier n positif ou nul. f est de classe \( C^n\) sur I. Alors il existe une fonction \( \epsilon (x) \) définie sur I, qui tend vers 0 quand x tend vers a, telle que, pour tout x de I :
\( f(x) = \sum_{k=0}^n \dfrac{f^{(k)}(a)}{k!} (x - a)^k + (x - a)^n \epsilon (x)\)
Cela veut dire qu'au voisinage d'un point quelconque a de l'intervalle, on peut exprimer une fonction f, qui répond aux conditions requises, en un polynome dont les facteurs sont les dérivées nième de la fonction en ce point a, et ceci en commettant une erreur qui est bornée par \( (x - a)^n \epsilon (x)\) quand x tend vers a. On démontre que si x tend vers a, alors \( \epsilon (x)\) tend vers 0. Dans les cours français, on note souvent ce reste \( O(x^n)\).

Prenons un exemple simple en posant f(x) = exp(x). Vous savez que la dérivée première de exp(x) est exp(x) et donc que la dérivée nième de exp(x) est exp(x). D'autre part, vous savez que exp(0) = 1. Si donc j'applique le théorème de Taylor-Young à exp(x) au voisinage de 0 (a = 0), j'obtiens:
\( e^x = \sum_{k=0}^n \dfrac{e^0}{k!} (x - 0)^k + (x - 0)^n \epsilon (x)\) soit \( e^x = \sum_{k=0}^n \dfrac{1}{k!} x^k + x^n \epsilon (x)\)
le dernier terme tendant vers 0.

Et bien sur, ce théorème est valable, sous les mêmes conditions, pour des fonctions à plusieurs variables !

Les développements limités

La notion de développement limité (DL) d'une fonction de classe \( C^n\) n'est pas abordée au lycée, et c'est pourtant un élément fondamental de l'analyse numérique et la base des schémas aux différences finies.

L'existence du développement limité (DL) d'une fonction f est un corollaire du théorème de Taylor-Young vu ci-dessus. On dit que f admet un DL d'ordre n au voisinage du point x0, qui s'exprime par :
\( f(x) = \sum_{k=0}^n \dfrac{f^{(k)}(x_0)}{k!} (x - x_0)^k + O_{x_0}((x - x_0)^n)\)
Par exemple, le DL d'ordre 2 de f(x) au voisinage de x0 sera :
\( f(x) = f(x_0) + f'(x_0)(x - x_0) + \dfrac{1}{2} f''(x_0) (x - x_0)^2 + O_{x_0}((x - x_0)^2)\)

Notions de stabilité en modélisation

Une des qualités essentielles d'un schéma numérique, c'est sa stabilité. Lorsqu'on emploie le terme "stabilité" en physique numérique, il est important de préciser de quoi on parle ! On peut rencontrer plusieurs types d'instabilité.

La première concerne la physique elle-même. Si le problème traité concerne un système chaotique, les solutions au problème seront instables et aucune opération sur le modèle mathématique et les méthodes numériques utilisées n'y changera rien! On parle alors de problème "chaotique", et il faudra se débrouiller pour tenir compte des résultats imprévisibles.

Un problème peut être "mal conditionné". C'est généralemnt du au choix d'un modèle mathématique non pertinent ou mal conçu. Le conditionnement d'un problème est une notion que l'on doit à Turing et qui traduit très schématiquement la petitesse de l'erreur commise lors des calculs numériques.

La méthode numérique choisie n'est pas adaptée. C'est le cas lorsque le schéma numérique choisi produit des erreurs numériques dues aux arrondis ou/et à la discrétisation. On en verra un exemple dans le cas du schéma d'Euler explicite.

Quoiqu'il en soit, si un problème est chaotique, on ne pourra produire de la stabilité ni avec le modèle mathématique, ni avec le schéma numérique. Si le problème est mal conditionné, aucun schéma numérique ne pourra y remédier. Et si le schéma numérique est instable, il faut en changer !

Les bases des méthodes aux différences finies

La discrétisation de l'espace et du temps

Les principes

En mathématique, et pas plus loin que ci-dessus dans l'énoncé du théorème de Taylor-Young, nous supposons habituellement les fonctions continues. Tous les outils de calcul différentiel "classique" sont basés sur ce critère. Notre problème est que la continuité n'est guère compatible avec le calcul numérique, et plus particulièrement avec le calcul numérique sur ordinateur.

Prenons le cas du calcul de la dérivée d'une fonction. Sa définition fait appel au concept de limite quand la valeur de h tend vers 0. Mathématiquement, on voit bien ce que cela veut dire: on trouvera toujours un réel aussi proche de zéro que l'on souhaite. Oui mais sur un ordinateur, les réels n'existent pas, ils sont "float". Et il existe un minorant et un majorant de l'ensemble des "float". Nous serons donc obligés de traduire "tend vers 0" par "plus grand ou égal au nombre le plus petit disponible sur mon ordinateur".

Il y a pire : sur un ordinateur, la continuité n'existe pas ! Appelons a le plus petit nombre "float" disponible sur mon ordinateur, -2,2*10^-308 sur mon PC un double précision en C. Construisons le nombre b = a + a. Si j'étais dans les réels, je pourrais trouver une infinité de nombres entre a et b. Sur mon ordinateur, il n'y a RIEN entre a et b.

Lorsque je manipule l'espace et le temps en physique numérique, ce ne sont jamais un espace et un temps continus. Ils sont découpés en segments plus ou moins minuscules dont seules les extrémités ont un sens. C'est l'opération de discrétisation de l'espace et du temps sur laquelle on va revenir.

Une remarque philosophique : il n'est pas impossible que la "réalité" de notre monde physique soit ainsi. L'espace est discret pour le physicien dans la mesure où nous ne disposons pas de physique applicable sous la dimension de Planck. Et rien ne prouve que le temps ne soit pas non plus "discret". Peut-être que finalement, les modèles discrets de la physique numérique ne sont pas aussi éloignés de la "réalité" que ça...

Comment procède-t-on ?

Imaginons que je veuille résoudre une EDP dépendante du temps, dans un espace de dimension 1, l'évolution de la température dans une barre fine par exemple. Mon espace de résolution de l'EDP sera donc de dimension 2 (en x et en t), un plan. Il me faut passer d'une surface continue, que je ne sais pas traiter avec un ordinateur à un dispostif discret. Pour réaliser ce dispositif, je vais mailler mon plan, c'est à dire tracer une grille sur ce plan, dont chaque maille aura pour dimension \( \Delta x\) et \( \Delta t\). Voilà comment cela se présente:

Sur cette grille, je repère chaque noeud par ses coordonnées d'espace (en abscisse) et de temps (en ordonnée), le couple \( (j \Delta x , n \Delta t) \), j et n variant entre 0 et la taille de la grille.

Le principe est de ne plus calculer la solution de l'EDP en tous points du domaine, ce que suppose la continuité, mais seulement à chaque noeud de la grille. Et c'est là qu'interviennent formule de Taylor et DL ! En chaque noeud, nous allons approximer la solution réelle de l'EDP en utilisant ces outils mathématiques.

Les principes des schémas de différences finies

Appelons u une fonction de classe \( C^4 \). Le théorème de Taylor nous permet d'écrire, avec h tendant vers 0 :
\( u(x + h) = u(x) + hu'(x) + \dfrac{h^2}{2}u''(x) + \dfrac{h^3}{4}u'''(x) + O(h4) \)
et son symétrique :
\( u(x - h) = u(x) - hu'(x) + \dfrac{h^2}{2}u''(x) - \dfrac{h^3}{4}u'''(x) + O(h4) \)

En combinant ces deux expressions, j'obtiens trois approximations de u'(x) qui me seront d'une grande utilité :
- l'approximation à droite (backward) : \( u'(x) = \dfrac{u(x+h) - u(x)}{h} + O(h) \),

- l'approximation à gauche (forward) : \( u'(x) = \dfrac{u(x) - u(x-h)}{h} + O(h) \),

- l'approximation centrée : \( u'(x) = \dfrac{u(x+h) - u(x-h)}{2h} + O(h^2) \), qui est plus précise car en \(O(h^2) \)

et l'approximation de u''(x) : \( u''(x) = \dfrac{u(x+h) - 2u(x) + u(x-h)}{h^2} + O(h^2) \).

Tous les schémas de DF sont basés sur ces quatre définitions, adaptées aux dérivées partielles. Pour une fonction u de deux variables x et y et classe \( C^4 \), on a pareillement:
- l'approximation à droite : \( \dfrac{\partial u}{\partial x} (x,y) = \dfrac{u(x + \Delta x,y) - u(x,y)}{ \Delta x} + O(\Delta x) \),

- l'approximation à gauche : \( \dfrac{\partial u}{\partial x} (x,y) = \dfrac{u(x,y) - u(x - \Delta x,y)}{ \Delta x} + O(\Delta x) \),

- l'approximation centrée : \( \dfrac{\partial u}{\partial x} (x,y) = \dfrac{u(x + \Delta x,y) - u(x - \Delta x,y)}{ 2 \Delta x} + O(\Delta x^2) \),

et l'approximation de la dérivée seconde : \( \dfrac{\partial^2 u}{\partial x^2} (x,y) = \dfrac{u(x + \Delta x,y) - 2 u(x,y) + u(x - \Delta x,y)}{\Delta x^2} + O(\Delta x^2) \).

qui sont bien sur transposables aux dérivées partielles par rapport à y.

Les différents schémas des DF

Posons le problème

Pour illustrer l'utilisation des schémas aux DF, je choisis une EDP que nous connaissons déjà : l'équation de diffusion thermique 1D, qui s'exprime par \( \dfrac{\partial T}{\partial t} = \dfrac{\partial^2 T}{\partial x^2} \) dans sa forme adimensionnée. Si vous voulez en savoir un peu plus sur cette équation, je vous invite à consulter la page qui lui est consacrée sur tangenteX.

Il s'agit de modéliser la diffusion thermique dans une barre fine et isolée, la température initiale de la barre est uniforme et fixée à Tinit et les extrémités étant maintenues à une température constante notée Text. Il s'agit d'un problème aux conditions aux limites de Dirichlet, pour une EDP parabolique :
\( \left\{ \begin{array}{lr} \dfrac{\partial T}{\partial t} = \dfrac{\partial^2 T}{\partial x^2} \; avec \; (x,t) \in ]0,1[ \times ]0,+\infty[ \\ T(x,0) = Tinit \\ T(0,t) = Text \\ T(1,t) = Text \end{array} \right. \)

Ce problème n'est pas chaotique et est bien conditionné. Reste donc la possibilité d'une instabilité liée au schéma numérique choisi... Je vous propose ci-dessous trois schémas numériques, sachant qu'il en existe bien d'autres (Lax, Lax-Wendroff, leap-frog, Dufort-Frankel, etc.), qui dérivent tous de ceux-ci.

Le schéma d'Euler explicite

Le schéma à construire est un schéma d'ordre 1 en temps et d'ordre 2 en espace, comme l'ordre des dérivées dans l'EDP. Je construis un maillage régulier, avec \(x_j = j \Delta x \) et \( t_n = n \Delta t \).

Pour discrétiser la dérivée partielle en temps du membre de gauche, j'utilise l'approximation à droite de cette dérivée, qui donne :
\( \dfrac{\partial T(x_j,t_n)}{\partial t} \approx \dfrac{T(x_j, t_n+\Delta t) - T(x_j,t_n)}{\Delta t} \)

Puis je discrètise la dérivée partielle du membre de droite d'ordre 2, par l'approximation centrée des dérivées secondes, ce qui me donne:
\( \dfrac{\partial^2 T(x_j,t_n)}{\partial x^2} \approx \dfrac{T(x_j + \Delta x,t_n) - 2T(x_j,t_n) + T(x_j - \Delta x,t_n)}{\Delta x^2} \)

D'où mon EDP discrétisée :
\( \dfrac{T(x_j, t_n+\Delta t) - T(x_j,t_n)}{\Delta t} = \dfrac{T(x_j + \Delta x,t_n) - 2T(x_j,t_n) + T(x_j - \Delta x,t_n)}{\Delta x^2} \)

Je vais simplifier l'écriture en posant \( T(x_j, t_n) = T_j^n \), ce qui me donne :
\( \dfrac{1}{\Delta t}(T_j^{n+1} - T_j^n) = \dfrac{1}{\Delta x^2}(T_{j+1}^n - 2 T_j^n + T_{j-1}^n) \)
ou encore en posant \( \lambda = \dfrac{\Delta t}{\Delta x^2}\) et regroupant les termes en indice n :

\( T_j^{n+1} = (1 - 2 \lambda)T_j^n + \lambda (T_{j+1}^n + T_{j-1}^n) \)

Schématiquement, on peut représenter ce schéma par :

Vous comprenez le nom du schéma : je calcule la température en un noeud et à un instant donné en utilisant la température à l'instant précédent pour les noeuds situés en dessous et de part et d'autre du noeud courant.

Pour la petite histoire, on nomme ce schéma le schéma FTCS (Forward Time Centred Space), car comme son nom l'indique il est approximé en avant pour le temps (approximation à droite) et centré pour l'espace.

Ce schéma est très simple, mais présente un inconvénient majeur: il est instable ! Il reste stable sous la condition que \( \Delta t \lt (\Delta x^2)/(2.0)\), condition s'appelle la condition de Courant-Friedrich-Levy, ou CFL. Autre défaut, il est relativement peu précis, car d'ordre \(O(\Delta t) \) sur le temps.

Ce schéma FTCS reste le plus simple à programmer mais son utilisation est délicate à cause de la condition CFL, dont il faut vérifier à chaque changement de pas sur x ou t qu'elle est bien vérifiée. Vous en trouverez une implémentation en Python dans la page sur l'équation de diffusion thermique.

Le schéma d'Euler implicite

Il est un peu plus compliqué que le schéma explicite mais il présente l'énorme avantage d'être intrinsèquement stable. Sa précision est la même que le schéma explicite. C'est une variante du schéma FTCS, qui consiste à calculer la dérivée temporelle à gauche (d'où le nom anglais BTCS) plutôt que de la calculer à droite comme dans le FTCS . On peut aussi calculer la dérivée spatiale au temps n+1 au lieu de la calculer au temps n.

Calculons la dérivée la dérivée temporelle à gauche, pour obtenir une relation implicite (d'où son nom français). J'obtiens l'équation discrète:

\( \dfrac{T_j^{n+1} - T_j^{n}}{\Delta t} = \dfrac{T_{j+1}^{n+1} - 2T_j^{n+1} + T_{j-1}^{n+1}}{\Delta x^2} \)
Le terme \( T_j^{n+1} \) ne dépend plus uniquement de \( T_j^n \) mais aussi de la valeur des noeuds voisins.

Afin de faciliter l'écriture du système linéaire, je vais réarranger mon équation comme ceci:

\( T_j^n = -\lambda T_{j-1}^{n+1} + (1 + 2 \lambda)T_j^{n+1} - \lambda T_{j+1}^{n+1} \)

Nous obtenons un système linéaire entre les \( T_j^{n+1} \) et \( T_j^n \). Il est de la forme \( T^n = A T^{n+1} \) et finalement \(T^{n+1} = A^{-1} T^n \) où la matrice A est la matrice tridiagonale des coefficients. Il se traite assez facilement avec les fonctions matricielles de Python.

Schématiquement, on peut représenter ce schéma par :

Je vous propose le script DFBCTS.py, qui utilise le schéma BCTS pour calculer une solution de l'équation de diffusion thermique, basé sur les scripts de calcul de la page sur l'équation de diffusion thermique.

En voici les instructions spécifiques, tout d'abord la construction de la matrice tridiagonale, qui utilise une fonction déjà vue :

A = MatDiag(Nx,-c,1 + 2*c)

puis l'inversion de cette matrice par la fonction inv() de scipy :

AI = inv(A)

Et enfin la résolution du système pour chaque pas de temps :

for t in range(0,Nt):

T[t+1,:] = dot(AI,T[t,:])

Nous obtenons une courbe identique à celle des autres méthodes, avec un avantage certain : pas besoin de réfléchir à la CFL lorsqu'on change de pas temporel ou spatial ! Attention toutefois à la consommation mémoire: 8 Go vivement conseillé !

Le schéma de Crank-Nicolson

Le schéma de Crank-Nicolson, qui date de 1947, mixte les deux approches explicites et implicites. La dérivée temporelle du membre de gauche est discrétisée par une approximation forward :
\( \dfrac{\partial T(x_j,t_n)}{\partial t} \approx \dfrac{T_j^{n+1} - T_j^n}{\Delta t} \)

La dérivée seconde spatiale est discrétisée par une valeur moyenne des approximations forward et backward :
\( \dfrac{\partial^2 T(x_j,t_n)}{\partial x^2} \approx \dfrac{1}{2 \Delta x^2}((T_{j+1}^{n+1} - 2T_j^{n+1} + T_{j-1}^{n+1}) + (T_{j+1}^n - 2T_j^n + T_{j-1}^n)) \)

Ce qui me donne l'équation discrètisée :
\( \dfrac{1}{\Delta t} (T_j^{n+1} - T_j^n) = \dfrac{1}{2 \Delta x^2}((T_{j+1}^{n+1} - 2T_j^{n+1} + T_{j-1}^{n+1}) + (T_{j+1}^n - 2T_j^n + T_{j-1}^n)) \)

En posant comme ci-dessus \( \lambda = \dfrac{\Delta t}{\Delta x^2}\), et en passant à droite de l'équation tous les termes connus d'indice n pour le temps, j'obtiens:

\( 2(\lambda + 1)T_j^{n+1} - T_{j+1}^{n+1} - T_{j-1}^{n+1} = 2(\lambda - 1)T_j^n + T_{j+1}^n + T_{j-1}^n \)

Cette équation discrète peut s'écrire sous la forme d'un système linéaire avec \( AT^{n+1} = BT^n \) où A et B sont les matrices diagonales des coefficients. La résolution de ce système linéaire se fait avec les fonctions de Python dot() et solve(). Vous en trouverez une implémentation en Python dans la page sur l'équation de diffusion thermique.

Schématiquement, on peut représenter ce schéma par :

Le schéma de Crank-Nicolson est plus précis que les schémas d'Euler, mais il reste en partie implicite et peut nécessiter un long temps de calcul. On peut le rendre entièrement explicite et on l'appelle dans ce cas le schéma de Heun. C'est un schéma d'ordre deux en temps, dont la condition de stabilité est que \( 0 \lt \Delta t \lt \dfrac{2}{\lambda}\). Je ne décrirai pas ce schéma ici.

Les scripts Python

Les scripts Python étudiés dans cette page sont disponibles dans le package DFPython.zip :

DFBTCS.py : résolution de l'équation de la chaleur par le schéma DF BTCS

Pour conclure

Pour la plupart de nos travaux, l'utilisation du schéma FTCS sera suffisante. Il n'est pas très précis, et instable. Mais à condition de respecter la CFL lorsqu'on définit les pas de discrétisation, il rend de grands services. Il est surtout facile à mettre en oeuvre. Il peut également servir à la résolution numérique des EDO. C'est d'ailleurs ce schéma qui est utilisé dans la méthode d'Euler présentée en TS pour la résolution numérique des EDO.

Une information pour les élèves de prépa : les méthodes aux différences finies, avec les schémas explicites et implicites d'Euler ont fait l'objet du problème du concours CCP option PC de 2015, avec une application de transfert thermique...