Chapitre 3 : Méthode des différences finies (1D)

1. Problèmes stationnaires 1D
2. Problèmes instationnaires 1D
3. Consistance & stabilité
La méthode en 5 pas

PAS 1 : Définition du problème continu

PAS 2 : Introduction d’un maillage & valeurs nodales des champs

PAS 3 : Problème discrétisé, une équation algébrique par point

de maillage, approcher les dérivées avec des différences
finies de même ordre de précision partout.

PAS 4 : Ecriture matricielle du système linéaire

PAS 5 : Programmation du système linéaire & solution numérique

une fois (pb stationnaire), en boucle (pb instationnaire)

VALIDATION ?
1. Problèmes stationnaires 1D
1.1 Poisson 1D + CL Dirichlet

PAS 1 : Définition du problème continu

En premier lieu, on doit toujours proprement définir le problème continue, que l’on doit
résoudre. On commence avec le problème de Poisson en 1D, muni
de conditions limites de type Dirichlet (on impose la valeur de f sur les bords)

, sur les bords:

?
PAS 2 : Introduction d’un maillage & valeurs nodales des champs

Nous discrétisons l’espace en introduisant un maillage

maillage de M+1 pts

Les fonctions/champs sont représentées par des ensembles de valeurs nodales

discrétisées

M+1 valeurs nodales de g, connues

M+1 valeurs nodales de f, recherchées

Notation

Souvent, on travaille avec des maillages réguliers

pas d’espace
PAS 3 : Ecrire le problème discrétisé

La philosophie est simple : 1 équation algébrique par point de maillage

pour tout point de bord (point extérieur)

on écrit une représentation fidèle de la condition aux limites
= facile pour un condition de dirichlet

Ici donc 2 équations pour 2 points les bords

pour tout point intérieur Notation

on impose une version discrète de l’équation différentielle

On remplacera les dérivées par des différence finies. Pour un maillage générale ceci
donne (chapitre 1)

pour un maillage uniforme, on simplifie l’expression à

On obtiendra donc M-1 équations algébriques, pour M-1 points intérieurs

PAS 4 : Ecriture matricielle
On a obtenu donc M+1 eqs pour M+1 valeurs nodales. Ici, ces équations (maillage régulier)
seront

eq point 0

eq point 1

eq point 2

...

eq point i

...

eq point M-1

eq point M

que l’on peut assembler sous forme matricielle

 Notant ceci donne

eq point 0 1
eq point 1
eq point 2

... ... = ...

eq point M-1
eq point M 1

matrice taille (M+1) x (M+1) vecteur vecteur

= Un problème linéaire de la forme

PAS 5 : Programmation

Il suffit de définir la matrice A et le vecteur g, afin de pouvoir calculer la solution f à l’aide

de l’ordinateur. On utilisera f = A\g en Matlab et des bibliothèques d’algorithmes ailleurs.
VALIDATION ?

Si possible, on trouvera un test qui permet de retrouver des valeurs que l’on calcule
à la main. Ceci validera le programme, qui pourra ensuite être utilisé pour résoudre
d’autres problèmes
Le plus facile c’est de proposer une fonction f(x), par exemple

on calcule
alors

TEST : si on utilise cette fonction et ces constantes et pour

définir le vecteur colonne alors on devra trouver en le vecteur
colonne les valeurs nodales

Ce genre de test de validation permet d’évaluer comment l’erreur décroit avec

un nombre de points augmentant et de retrouver l’ordre du schéma numérique.
code poisson1d_dirichlet.m

+ validation
1.2 Poisson 1D + CL de Neuman

PAS 1 : Définition du problème continu

modif

CL:

PAS 2 : Introduction d’un maillage & valeurs nodales des champs

avec

PAS 3 : Problème discrétisé, une équation algébrique par point de maillage

modif
eq. point 0

précision
eq. points i d’ordre 2
partout

eq. point M
PAS 4 : Ecriture matricielle

eq point 0 modif
eq point 1
eq point 2

... ... ...

=

eq point M-1
eq point M

matrice taille (M+1) x (M+1) vecteur vecteur

PAS 5 : Programmation

Pareil
VALIDATION ?

Selon l’idée précédente on propose une fonction f(x), par exemple

on calcule

TEST : si on utilise cette fonction et ces constantes et pour

définir le vecteur colonne alors on devra trouver en le vecteur
colonne les valeurs nodales
code poisson1d_neuman.m

+ validation
1.3 Problème aux valeurs propres + CL Dirichlet

PAS 1 : Définition du problème continu

CL:

Pour quels existe-t-il des solutions f non-nulles

PAS 2 : Introduction d’un maillage & valeurs nodales des champs

avec

PAS 3 : Problème discrétisé, une équation algébrique par point de maillage

eq. point 0

eq. points i

eq. point M
PAS 4 : Ecriture matricielle

On cherche à écrire le problème discrétisé comme un problème aux

valeurs propres (généralisé) matriciel
Notons

1 0
1
... = ...
1
1 0

PAS 5 : Programmation

Il suffit de définir la matrice A et B, afin de pouvoir calculer les fonctions et les valeurs
propres à l’aide de l’ordinateur. On utilisera [F,LAMBDA]=eig(A,B) en Matlab
VALIDATION

Ici on connait la solution analytique

on renomme

Solution

Les conditions aux limites fixent

soit

Les fonctions propres sont donc

Est-ce qu’on les retrouve ?

code valeurspropres1d_general.m
+ validation
PAS 4’ : Ecriture matricielle, qui élimine les points de bord

Des conditions limites homogènes peuvent toujours être éliminées.

=
... ...

de la forme

PAS 5’ : Programmation

Il suffit de définir la matrice A et B, afin de pouvoir calculer les fonctions et les valeurs
propres à l’aide de l’ordinateur. On utilisera [F,LAMBDA]=eig(A) en Matlab
code valeurspropres1d_reduit.m
+ validation
1.4 Problème aux valeurs propres + CL différentes

PAS 1 : Définition du problème continu

CL:

Pour quels existe-t-il des solutions f non-nulles ?

PAS 2+3+4+5 ou 4’+5’ : A vous maintenant !

2. Problèmes instationnaires 1D
2.1 Advection 1D : schéma explicite UPWIND

PAS 1 : Définition du problème continu

1 CI : Solution ?
1 CL :

La fonction avance vers des valeurs de x croissantes, sans se déformer.

Peut-on retrouver ce comportement dans une solution numérique du
problème.
PAS 2 : Introduction d’un maillage, de temps discrets & des valeurs nodales des champs

maillage : M+1 points d’espace

avec uniforme

temps : N+1 temps

avec uniforme

PAS 3 : Ecrire le problème discrétisé

philosophie : 1 eq par point de maillage ET par temps discret

sur les points de bords, on impose les CL
sur les points intérieurs, on impose une traduction diff. finies de l’EDP
INITIALISATION : on commence par traduire la condition initiale qui fixe l’état à n=0

AVANCEMENT TEMPOREL : on cherche des équations algébriques pour tous les points,
qui permettent de trouver le champ à tous les temps suivants à partir du champ au(x) temps
précédent(s).
Sur le point de bord, à l’entrée du domaine (a>0), on dispose d’une condition
de Dirichlet qui fixe la valeur de la fonction à tout temps.

Sur les points intérieurs on impose une version discrétisée de l’EDP, ici

lié au choix du schéma

Dans le “schéma numérique” UPWIND, on discrétise la dérivée temporelle utilisant

une formule (F) et la dérivée spatiale utilisant une formule (B)

(F) (B)
Ceci donne

ou encore avec un nombre sans dimension

on comprend mieux
pourquoi.

En résumé nous avons donc pour tous les points et pour tout

Ceci est une formule récursive, qui permet de calculer les valeurs nodales du champ à
temps (n+1) à partir de celles du temps n précédent.

On appelle ce schéma explicite car la formule est explicite : on ne doit résoudre aucune
équation afin de trouver l’état suivant du champ.
PAS 4 : Ecriture matricielle

INITIALISATION

AVANCEMENT TEMPOREL Boucle sur

“Opérateur” ou matrice
d’évolution
PAS 5 : Programmation

Nous implémentons le précédent schéma dans le programme Matlab

testupwind.m.

et a est une vitesse donnée à l’entrée de la fonction.

code testupwind.m
+ validation

Le code retrouve bien la solution analytique si C<1

Mais instabilité numérique si C>1
A chaque pas de temps, on fait des erreurs d’approximation et ces erreurs peuvent
s’accumuler pour mener à un phénomène d’instabilité numérique.

Un code numérique doit être stable si on veut l’utiliser en pratique. Nous disposons
d’outils théoriques pour caractériser les domaines de stabilité d’un code = > Section 3
2.2 Diffusion 1D : schéma implicite CRANK-NICOLSON

PAS 1 : Définition du problème continu

PAS 2 : Maillage, temps discrets & discrétisation des champs

même qu’avant

PAS 3 : une équation par point de maillage et par temps

INITIALISATION : fixer l’état à n=0

AVANCEMENT TEMPOREL :
Sur les points de bords on écrit les conditions aux limites. Ici, pas de difficultés car on
a des CL de type Dirichlet. On les exprime à temps n+1

Pour tous les autres points intérieurs on essaie d’écrire une formule algébrique récursive,
qui est consistante avec l’EDP de diffusion. Ici le schéma de Crank-Nicolson dicte que

où est un nombre sans dimension important.

Pour comprendre ce schéma, on peut voir le processus de discrétisation en deux étapes

séparés (d’abord temps puis espace ou inversement) au tableau

Ce schéma est implicite, car on ne peut pas directement calculer le champ a temps n+1
à partir des valeurs nodales aux temps précédents.
 pour faire le lien avec la notation matricielle, il convient de réorganiser les termes ainsi

Tous les termes qui font appel aux valeurs nodales de f à temps n+1, sont mis à gauche
tous les autres à droite.

PAS 4 : écriture matricielle

INITIALISATION

AVANCEMENT TEMPOREL Boucle sur

(page suivante)
AVANCEMENT TEMPOREL Boucle sur

Si la matrice inverse est facile à calculer alors

opérateur d’évolution

Sinon, il faut résoudre un grand système linéaire à chaque pas de temps.

PAS 5 : Programmation

On programme un cas concret. Un matériau solide, initialement à température basse

est mise en contact avec une source chaude à temps t=0.

temps initial

temps ultérieurs
la chaleur diffuse
x
0 L
on souhaite simuler le transitoire. On choisit par exemple L=10 et

code testCNdiff.m

il n’y a pas de source de chaleur dans la barre.

3. Consistance & stabilité
3.1 Schémas numériques

Dans cette section on suppose qu’on résolve un problème instationnaire 1D tel que

(advection 1D)

(diffusion 1D)

(diff-adv 1D)

A l’aide de la méthode des différences finies. Comme d’habitude on note

pour le champ discret sur les points de maillage et les temps discrets. Dans toute cette
section, on travaille avec un pas d’espace et un pas de temps constant.
On peut imaginer une ‘infinité’ de schémas numériques qui permettent d’obtenir un
algorithme pour marcher en temps une solution. Quelques exemples

(FTCS adv)

(CN adv)

(I3N diff)

(autres exemples sur la feuille donnée et sur DOKEOS)

Chacun de ces schémas peut être plus au moins précis, plus au moins facile à utiliser
et plus ou moins stable en fonction des paramètres

La caractérisation de ces propriétés fait l’objet de cette section 3 du chapitre.

Les problèmes étudiés ici sont linéaires et 1D en espace. On peut alors toujours réécrire
les schémas sous une forme générale

2 niveaux
temporaux
(n+1 et n) (on ignore
les sources)
3 niveaux
temporaux
(n+1,n,n-1)

L’idée est ici de mettre à gauche tous les termes à temps n+1, à droite tous les termes
aux temps précédents (n,n-1,…). Chaque schéma se caractérise par un certain nombre
de niveaux temporaux et de points d’espace impliqués.

Un schéma est explicite si

ou

Dans ce cas, on obtient directement la solution à un temps n+1 (cf. section 2.1) à
partir de la solution aux temps précédents. Dans le cas opposé, le schéma est
implicite. La solution à temps n+1 se trouve comme la solution d’un système linéaire
(FTCS adv)

2 niveaux temporaux
3 points d’espace
explicite

(CN adv)

2 niveaux temporaux
3 points d’espace
implicite

(I3N diff)

3 niveaux temporaux
3 points d’espace
implicite
3.2 Consistance & ordre d’un schéma numérique

Il existe une multitude de schéma permettant de résoudre une même équations et

pour un schéma donné, il faut se convaincre que le schéma soit consistant, c.a.d.

Est-ce que le schéma représente bien l’équation exacte dans la limite ou

le pas d’espace et la pas de temps tendent vers zero, et partout
dans le domaine considéré ?

Pour tester la consistance, il suffit d’appliquer le schéma à une fonction continûment

différentiable . Ensuite on propose des développements limités

etc…

autour d’un seul endroit (ici ) et un seul temps (ici ) et on remplace ces expression
dans le schéma proposé

(comme en chapitre 1, mais avec des dérivées partielles en x et t)

Pour un schéma quelconque ceci résultera dans

( schéma ) = (equation exacte) + (résidu)

Un schéma ne sera que consistant si le résidu décroit à zéro pour et ce résidu

sera alors (le plus souvent) de la forme

(résidu) = = (terme) + (terme)

Ici p est un exposant qui indique l’ordre spatial du schéma, q indique l’ordre temporel.
On appelle min(p,q) l’ordre global du schéma.

(FTCS adv)

On remplace

….
 (equation exacte) (residu)

Le schéma est consistant, car on retrouve bien l’equation exacte dans la limite
. L’ordre temporel est 1, l’ordre spatial est 2, l’ordre global 1

D’autres exemples seront donnés en TD. Attention: pousser les DL suffisamment loin,
afin de pouvoir conclure sur l’ordre d’un schéma.

En théorie, la précision d’un schéma d’ordre supérieur est plus élevée, mais attention.
Ceci n’est que vrai dans la limite

exemple:

Parfois un schéma d’ordre 1, peut être plus précis qu’un schéma d’ordre sup à dx,dt fini
3.3 Convergence & stabilité d’un schéma numérique
3.3.a Définition : convergence d’un schéma

Dans la section 2 de ce chapitre, nous avons donné de nombreux exemples, dans

lesquels on rassemblait les équations des points discret sous forme matricielle. On
utilisera cette notation ici

solution numérique solution exacte (suffix ex)

Un schéma numérique sera dit convergent pour une norme donné (notée ), si
pour toute solution initiale et pour tout n
3.3.b. La consistance n’implique pas la convergence

Nous avons vu que le schéma donnant la boucle d’avancement temporel se laisse

écrire sous forme matricielle. Pour un schéma à deux niveaux temporaux, ceci donne

“Opérateur” ou matrice
d’évolution

soit

(1)
résidus
On peut écrire une équation pour soit-disant la solution exacte (suffix ex)

(2)

qui fait intervenir les résidus rencontrés dans l’étude de la consistance

On introduit l’erreur comme la différence entre la solution exacte et la solution numérique

En soustrayant les equations (2)-(1), on obtient une équation pour l’évolution cet l’erreur

Même quand l’erreur est zéro au départ, elle ne le restera pas forcément au cours des
itérations.

….
Si on introduit une norme ( notée ) on peut mesurer l’erreur de manière plus
globale et on peut aussi la majorer.

La consistance d’un schéma garantie que

mais ne garantie pas que

Cette dernière quantité peut même croître de manière incontrôlée lorsque

et en conséquence l’erreur ne décroitra pas forcément vers zéro:

Il faut donc imposer une condition supplémentaire pour assurer la convergence:

un schema numérique doit être stable.
3.3.c. Condition de stabilité

Un schéma numérique est stable pour la norme si pour tout vecteur , il existe une
constante indépendante de telle que

On voit ici que cette condition est associé à une durée maximale de la simulation T.
Avec cette relation, on simplifie (cf. slide précédent)

Utilisant maintenant la consistance du schéma

On trouve alors que le schéma converge

Théorème de Lax: Si un schéma consistant d’ordre (global) p est stable, l’erreur

tend vers 0 comme lorsque à fixe
3.4. Etude de stabilité par analyse de Fourier
3.4.a Introduction

L’étude de stabilité d’un schéma peut se faire de manière précise en étudiant le rayon
spectral de l’opérateur d’évolution. Ceci fait intervenir les notions de valeurs et
vecteur propres et sera un peu trop avancé pour ce cours (traité en M2).

Si on s’affranchit de l’effet des conditions aux limites, on peut utiliser l’analyse de

Fourier afin de trouver une condition de stabilité plus simple, qui s’exprime comme

On appelle le facteur d’amplification du schéma et k est ici un nombre

d’onde de l’espace de Fourier. Le calcul de ce facteur se fait directement à partir de
la définition du schéma.
3.4.b Quelques rappels sur l’analyse de Fourier

Soit une fonction réelle intégrable sur l’axe réel. On définit la transformée de
Fourier ici comme

La transformée de Fourier inverse est définie comme

Suite à l’utilisation de (au lieu de ) dans les exponentielles on a besoin d’un

facteur de normalisation (pour les anciens L3 d’Orsay)

Pour parler de stabilité, on aura besoin d’introduire une norme . On choisit ici la
norme que l’on définit comme

Une fonction qui à une norme < oo est carré-intégrable sur l’axe des réels:
On utilisera l’identité de Parséval

qui permet d’identifier la norme dans l’espace de Fourier à la norme dans l’espace
physique.
3.4.c. Facteur d’amplification d’un schéma à 2 niveaux
Soit un schéma numérique à 2 niveaux qui nous est donné, que l’on suppose de la
forme

Ici on ne spécifie pas comment varie m, mais vous pouvez bien imaginer que l’on
peut toujours écrire le schéma sous cette forme. On oublie l’existence de bords.
On construit maintenant des fonctions constantes par morceaux:

Pour cette fonction on peut, grâce au schéma numérique, écrire que

Toute dépendance de l’index j a disparue. On multiplie par et on intègre sur
l’axe réel afin de prendre la transformée de Fourier de l’équation. Par quelques manipula-
tiens basiques (propriété de translation)

<=>

Le facteur d’amplification est défini comme la quantité

Donnons un exemple, pour illustrer comment on procède en pratique

(FTCS adv)

On commence par réécrire le schéma sous forme générale

Ensuite on remplace

, ,

Ceci donne directement

Le facteur d’amplification devient

On constate
3.4.e. Stabilité

Nous avons déjà défini la norme d’une fonction. La norme d’un vecteur colonne
est la norme Euclidienne

C’est ce norme que l’on choisit afin de traduire la condition de stabilité (3.3.c)

On peut maintenant relier cette condition à une condition sur

Revenant sur la définition de ce facteur on avait

Avec la norme dans l’espace de Fourier, on peut immédiatement écrire

Application de l’identité de Parseval, permet de passer dans l’espace physique

(1)

et si on rappelle la définition des fonctions et on peut reconnaître

que

En conséquence, on peut simplifier l’équation (1) en

Si on boucle cette condition, on obtient donc

soit

Si on peut majorer

on aura la propriété de stabilité désirée

Stabilité & Théorème de von Neumann

Un schéma à 2 niveaux temporaux sera stable pour la norme s’il existe une
constante indépendante de telle que

Le théorème de von Neumann permet de relier cette stabilité à une condition sur le
facteur d’amplification du schéma

Un schéma à 2 niveaux est stable si et seulement si

avec le facteur d’amplification, une constante et à fixé

(Démo au tableau)
La notion de stabilité présentée ci-dessus, dépend du nombre de pas de temps
maximal: . Dans un cas concret, on aimerait avoir des garanties sur la
stabilité pour un nombre de pas de temps quelconque, d’où la notion de stabilité
asymptotique.

En TD, on ne s’intéressera qu’à la stabilité asymptotique.

Stabilité asymptotique & Critère de von Neumann

Un schéma à 2 niveau temporaux sera asymptotiquement stable pour la norme

s’il existe une constante indépendante de telle que

Le second théorème dit critère de von Neumann permet de relier cette stabilité asymp-
totique à une condition plus restrictive sur le facteur d’amplification du schéma

Un schéma à 2 niveaux est asymptotiquement stable si et seulement si

avec le facteur d’amplification et à fixé

(FTCS adv)

On commence par réécrire le schéma sous forme générale

Ensuite on remplace

, ,

Ceci donne directement

Le facteur d’amplification devient

On constate

Le schéma FTCS est inconditionnellement instable

3.4.f. Stabilité d’un schéma à 3 niveaux
Si nous avons un schéma à 3 niveaux temporaux, on procède de manière très similaire.
Soit un tel schéma

Suivant exactement la même procédure que précédemment, nous obtenons après la

transformée de Fourier

On divise cette équation par et on remarque que

Le facteur d’amplification est solution d’un polynôme d’ordre 2

Et il y a donc 2 solutions et
On peut montrer que l’extension suivante du critère de von Neumann s’applique

Un schéma à 3 niveaux est asymptotiquement stable si et seulement si

avec le facteur d’amplification et à fixé

Exemple en TD, consultez la feuille avec les propriétés de stabilité de quelques schémas