Thème 4: résolution numérique des EDP linéaires

Problèmes statiques et dynamiques

Université Pierre et Marie CURIE

MNCS
Méthodes numériques et calcul scientifique

2 et 3 mars 2017

(UPMC) Thème 4 : Résolution numérique des EDP 2 et 3 mars 2017 1 / 57

Sommaire
1 Introduction
Généralités
Conditions aux limites
Problèmes statiques et dynamiques

(UPMC) Thème 4 : Résolution numérique des EDP 2 et 3 mars 2017 2 / 57

Sommaire
1 Introduction
Généralités
Conditions aux limites
Problèmes statiques et dynamiques
2 Exemple de problème statique linéaire : l’équation de Laplace
Discrétisation du problème
Reformulation vectorielle de la solution
Solution du problème

(UPMC) Thème 4 : Résolution numérique des EDP 2 et 3 mars 2017 2 / 57

Sommaire
1 Introduction
Généralités
Conditions aux limites
Problèmes statiques et dynamiques
2 Exemple de problème statique linéaire : l’équation de Laplace
Discrétisation du problème
Reformulation vectorielle de la solution
Solution du problème
3 Problèmes dynamiques
Discrétisation de l’espace et du temps
Algorithme explicite du premier ordre en temps
Algorithme implicite du premier ordre en temps
Algorithme de Crank–Nicolson

(UPMC) Thème 4 : Résolution numérique des EDP 2 et 3 mars 2017 2 / 57

Sommaire
1 Introduction
Généralités
Conditions aux limites
Problèmes statiques et dynamiques
2 Exemple de problème statique linéaire : l’équation de Laplace
Discrétisation du problème
Reformulation vectorielle de la solution
Solution du problème
3 Problèmes dynamiques
Discrétisation de l’espace et du temps
Algorithme explicite du premier ordre en temps
Algorithme implicite du premier ordre en temps
Algorithme de Crank–Nicolson
4 Équation des ondes

(UPMC) Thème 4 : Résolution numérique des EDP 2 et 3 mars 2017 2 / 57

 Introduction

Sommaire

1 Introduction
Généralités
Conditions aux limites
Problèmes statiques et dynamiques
2 Exemple de problème statique linéaire : l’équation de Laplace
Discrétisation du problème
Reformulation vectorielle de la solution
Solution du problème
3 Problèmes dynamiques
Discrétisation de l’espace et du temps
Algorithme explicite du premier ordre en temps
Algorithme implicite du premier ordre en temps
Algorithme de Crank–Nicolson
4 Équation des ondes

(UPMC) Thème 4 : Résolution numérique des EDP 2 et 3 mars 2017 3 / 57

 Introduction

Contexte

Équations différentielles ordinaires (EDO) : une variable indépendante

Équations aux dérivées partielles (EDP) : plusieurs variables
indépendantes
EDP omniprésentes en physique, chimie, sciences de la Terre, biologie :
mécanique des fluides, propagation des ondes, électromagnétisme,
phénomènes de diffusion...
Se limiter ici à des EDP linéaires régissant un champ scalaire u( #–
r , t).
Généralisations vectorielles possibles, et techniquement analogues.
Résolution numérique par la méthode des différences finies avec un
maillage régulier.
Mais il existe d’autres méthodes, éléments finis par exemple, mieux
adaptées dans le cas de domaines à géométrie complexe.

(UPMC) Thème 4 : Résolution numérique des EDP 2 et 3 mars 2017 4 / 57

 Introduction

Exemples d’EDP scalaires les plus connues

1 L’équation des ondes ou « de d’Alembert »
∂2u
− c2 ∆ u = 0 (1)
∂t2
2 L’équation de Schrödinger
∂u 1
i = − ∆u + V u (2)
∂t 2
3 L’équation de diffusion
∂u # –
= div[D grad u] = D ∆ u si D est uniforme (3)
∂t
4 Les équations de Poisson et de Laplace
∆u = ρ ∆u = 0 (4)
B Les équations de Navier-Stokes en mécanique des fluides sont en
général non-linéaires à cause du terme d’advection.
(UPMC) Thème 4 : Résolution numérique des EDP 2 et 3 mars 2017 5 / 57
 Introduction Généralités

Équations elliptiques/paraboliques/hyperboliques

Étude théorique des EDP : vaste domaine en mathématiques.

EDP linéaires (du second ordre) classées selon la forme des coefficients
devant les dérivées partielles :

∂2u ∂2u ∂2u

α + β + γ = ...
∂x2 ∂x∂y ∂y 2

β 2 − 4αγ > 0 hyperboliques : équation des ondes

β 2 − 4αγ = 0 paraboliques : diffusion ou Schrödinger
β 2 − 4αγ < 0 elliptiques : équations de Laplace, Poisson, Helmholtz

Propriétés mathématiques différentes

−→ méthodes numériques spécifiques.

(UPMC) Thème 4 : Résolution numérique des EDP 2 et 3 mars 2017 6 / 57

 Introduction Conditions aux limites

Conditions aux limites

L’équation est vérifiée dans un domaine D de l’espace.

Problème « bien posé » −→ conditions aux limites (CL).
Plusieurs types de conditions :
Dirichlet : u imposé sur le bord ∂D.
# –
Neumann : grad u. #–
n (flux) imposé sur ∂D.
Robin : plus général, combinaison linéaire des deux.
Cas général : conditions mixtes, i.e. différentes sur différentes parties
de ∂D.
Exemple de l’équation de la chaleur.
Dirichlet ⇔ parois isothermes (imposent la température)
Neumann ⇐ parois adiabatiques (transfert de chaleur nul)

(UPMC) Thème 4 : Résolution numérique des EDP 2 et 3 mars 2017 7 / 57

 Introduction Problèmes statiques et dynamiques

Problèmes statiques et dynamiques

Distinguer deux classes de problèmes selon la présence de la variable

temps :

Problèmes dynamiques : déterminer u( #–

r , t) dans D sur un intervalle
de temps + CL spatio-temporelles :
CI (conditions initiales) : u( #–
r , t = 0)
CL (conditions aux limites) : u( #– r = ·, t)

Problèmes statiques : déterminer u( #–

r ) en fonction des CL
(pas de variable temps explicite).
La solution statique peut correspondre à un régime stationnaire
atteint après un temps d’évolution assez long.

(UPMC) Thème 4 : Résolution numérique des EDP 2 et 3 mars 2017 8 / 57

 Introduction Problèmes statiques et dynamiques

Problèmes statiques

1 - Problèmes statiques : solution d’une EDP où le temps n’intervient pas.

Exemple : solution d’équilibre de l’équation de diffusion lorsque t → ∞.

Équation de Laplace !
∂2 ∂2
∆T = + T =0
∂x2 ∂y 2

(UPMC) Thème 4 : Résolution numérique des EDP 2 et 3 mars 2017 9 / 57

 Introduction Problèmes statiques et dynamiques

Problèmes dynamiques

2 - Problèmes dynamiques :

Deux méthodes simples seront exposées et analysées sur l’exemple de

l’équation de diffusion à 1-D :

∂u ∂2u
=D (5)
∂t ∂x2
avec diverses conditions initiales.

(UPMC) Thème 4 : Résolution numérique des EDP 2 et 3 mars 2017 10 / 57

 Un problème statique : l’équation de Laplace

Sommaire

1 Introduction
Généralités
Conditions aux limites
Problèmes statiques et dynamiques
2 Exemple de problème statique linéaire : l’équation de Laplace
Discrétisation du problème
Reformulation vectorielle de la solution
Solution du problème
3 Problèmes dynamiques
Discrétisation de l’espace et du temps
Algorithme explicite du premier ordre en temps
Algorithme implicite du premier ordre en temps
Algorithme de Crank–Nicolson
4 Équation des ondes

(UPMC) Thème 4 : Résolution numérique des EDP 2 et 3 mars 2017 11 / 57

 Un problème statique : l’équation de Laplace

Résolution de l’équation de Laplace dans un rectangle

Sur€le€bord€:€T(x,y)€connue y

Objectif : solution A l'intérieur T(x,y) inconnue

numérique approchée de vérifiant ∆T=0

∆T = 0 ∆x
à l’intérieur de D.

Simplification :
∆y
D rectangulaire,
de côtés ∆x et ∆y. x

Figure 1 – Géométrie du problème de

l’équation de Laplace : noter le
choix inhabituel de l’orientation des axes x et y.

(UPMC) Thème 4 : Résolution numérique des EDP 2 et 3 mars 2017 12 / 57

 Un problème statique : l’équation de Laplace

Exemple de carte d’isothermes

CL Dirichlet : températures imposées aux bords

˚
10C y
0 5 10 15 20 25
0
Figure 2 – Exemple de
solution de l’équation
14 13
12 11
10 ∆ T = 0, dans le cas où une
5 16
18
9
8
paroi est maintenue à 20◦ C,
20C
˚

7 une autre adjacente est à

10
6
5
10◦ C et les deux autres
4 parois sont à T = 0◦ C.
x 3
2
1
15 Les lignes indexées de 1 à 18
représentent les isothermes.
˚
0C
20

(UPMC) Thème 4 : Résolution numérique des EDP 2 et 3 mars 2017 13 / 57

 Un problème statique : l’équation de Laplace Discrétisation du problème

Maillage du domaine D
Discrétisation du domaine D : température évaluée en un nombre fini de
points. Le plus simple : grille de pas dx = dy = p
.'.'. (0,j' ).'.'.'.'.'.
1 .'.'.' j .'.'.'.'.'. nc y 0 j nc+1
0
1 y

.'.'.
∆x'=( nl+1)' p
.' .

p (i' ,j' )
i
i (i,j' )
(i,nc+1)
.'.'.'.'.'.

.'.'.'.'.'.
nl
p nl' +1
x x
∆y'=(nc+1)' p

Figure 3 – Point de grille ↔ (i, j). Figure 4 – Les bords du domaine

Intérieur du domaine : i = 1 → nl, (C.L.) sont numérotés : i = 0 pour le
j = 1 → nc. haut, i = nl + 1 pour le bas, j = 0 pour
la gauche, j = nc + 1 pour la droite.
(UPMC) Thème 4 : Résolution numérique des EDP 2 et 3 mars 2017 14 / 57
 Un problème statique : l’équation de Laplace Discrétisation du problème

Indexation 2D du domaine

À l’intérieur de la grille, on a :
nl = ∆x/p − 1 lignes horizontales.
nc = ∆y/p − 1 colonnes verticales
Chaque point est repéré par un couple d’indices (i, j) où le premier indice
note le numéro de la ligne et le second celui de la colonne.

On note Ti,j la température en ce point, de coordonnées

xi = i × p, yj = j × p.

Les inconnues du problème sont les n = nl × nc valeurs de Ti,j

où 1 6 i 6 nl, 1 6 j 6 nc.
Elles seront déduites des 2(nl + nc) valeurs de Ti,j imposées aux bords
i = 0 (nord) ou i = nl + 1 (sud)
j = 0 (ouest) ou j = nc + 1 (est).

(UPMC) Thème 4 : Résolution numérique des EDP 2 et 3 mars 2017 15 / 57

 Un problème statique : l’équation de Laplace Discrétisation du problème

Différences finies
Les dérivées partielles sont approchées par des différences finies, basées
sur des développements de Taylor. Le développement de Taylor :
∂T p2 ∂ 2 T p3 ∂ 3 T
T (x ± p) = T (x) ± p + +
 
± O p4 (6)
∂x 2 ∂x2 3! ∂x3
montre que :
∂2T
T (x + p) + T (x − p) = 2T (x) + p2 (x) +
 
O p4 (7)
∂x2
d’où :
∂2T T (x + p) + T (x − p) − 2 T (x)
(x) = + O p2
 
2 2
(8)
∂x p
ce qui conduit à :
∂2T Ti+1,j + Ti−1,j − 2 Ti,j
(x ) = +
 
i , yj O p2 (9)
∂x2 p2
dont la précision est du deuxième ordre en p.
(UPMC) Thème 4 : Résolution numérique des EDP 2 et 3 mars 2017 16 / 57
 Un problème statique : l’équation de Laplace Discrétisation du problème

Différences finies (suite)

Obtention alternative :

∂T p T (x ± p) − T (x)
 
x± ≈± (10)
∂x 2 p
qui conduit à :

∂2T 1 ∂T p ∂T p
    
(xi , yj ) ≈ xi + , yj − xi − , yj (11)
∂x2 p ∂x 2 ∂x 2
soit encore :
∂2T Ti+1,j + Ti−1,j − 2 Ti,j
(x ) = +
 
i , yj O p2 (12)
∂x2 p2

(UPMC) Thème 4 : Résolution numérique des EDP 2 et 3 mars 2017 17 / 57

 Un problème statique : l’équation de Laplace Discrétisation du problème

Différences finies : Laplacien 2D

On obtient ainsi le Laplacien 2D discrétisé à l’ordre 2 :

Ti+1,j + Ti−1,j + Ti,j+1 + Ti,j−1 − 4 Ti,j
∆ Ti,j = +
 
O p2 (13)
p2

Équation de Laplace ∆ T = 0 =⇒ système linéaire de N équations

reliant Ti,j aux 4 points voisins :

Ti+1,j + Ti−1,j + Ti,j+1 + Ti,j−1 − 4 Ti,j = 0 (14)

N = nc × nl valeurs à déterminer
à partir de 2(nc + nl) températures connues.
Certains points intérieurs ont des voisins sur les parois avec des
températures imposées qui vont constituer le second membre.
Noter que le pas p n’intervient plus.

(UPMC) Thème 4 : Résolution numérique des EDP 2 et 3 mars 2017 18 / 57

 Un problème statique : l’équation de Laplace Reformulation vectorielle de la solution

Réindexation 1-D du domaine (en RowMajor)

Ré-indexation 1D des points intérieurs : k = (i − 1)nc + j ⇒ vecteur Z.

Ti,j ⇐⇒ Zk 1 6 k 6 nl × nc

1N Nord j
y Figure 5 – Renumérotation des
1 2 nc points de grille pour représenter la
1W
solution sous la forme d’un vecteur :
2nc
Ouest

nc+1 2E le domaine est ré-indexé ligne par

ligne de la gauche vers la droite, puis
Est

i de haut en bas (comme une image

nl׀nc numérique) donc en RowMajor.
Les voisins au Nord et au Sud sont
k=(i€–1)nc+j
alors espacés de nc dans Z.
x 1S
Sud

(UPMC) Thème 4 : Résolution numérique des EDP 2 et 3 mars 2017 19 / 57

 Un problème statique : l’équation de Laplace Reformulation vectorielle de la solution

Formulation vectorielle

Pour un point sans contact avec les bords du domaine, l’équation (14)
s’écrit :
Zk−nc + Zk−1 − 4Zk + Zk+1 + Zk+nc = 0 (15)
Ne pas croire que le second membre du système linéaire est nul !
Pour un point en contact avec un bord, cette équation fait intervenir des
températures imposées : par exemple celles du bord Nord si i = 1.

T2N + Z1 − 4Z2 + Z3 + Z2+nc = 0 pour j = 2

=⇒ passer ces termes dans le second membre du système

Les CL de Dirichlet des bords créent le second membre
Chaque changement de membre d’une température imposée annule
un des coefficients des Zk0 .

(UPMC) Thème 4 : Résolution numérique des EDP 2 et 3 mars 2017 20 / 57

 Un problème statique : l’équation de Laplace Reformulation vectorielle de la solution

Matrice L2D

Réécrire le système d’équations sous la forme :

AZ=B (16)
où A = L2D est une matrice carrée N × N comportant
nl × nl blocs carrés de dimension nc × nc.
L2D est symétrique tridiagonale par blocs :
1 blocs de type L1D avec −4 au lieu de −2 sur la diagonale ;
2 blocs de type identité sur la sous-diagonale et la sur-diagonale.

Chaque ligne et chaque colonne contient 5 termes non-nuls,

sauf celles de numéro i × nc et i × nc + 1, où apparaissent les 0 des
contacts avec les bords latéraux,
et celles associées au premier et dernier bloc (bords haut et bas).

(UPMC) Thème 4 : Résolution numérique des EDP 2 et 3 mars 2017 21 / 57

 Un problème statique : l’équation de Laplace Reformulation vectorielle de la solution

Matrice L2D avec nl = 3 et nc = 4

 x
−4 1 0 0 1 0 0 0 0 0 0 0 

 1 −4 1 0 0 1 0 0 0 0 0 0 
 
 0 1 −4 1 0 0 1 0 0 0 0 0 
 
 0 0 1 −4 0 0 0 1 0 0 0 0 
 
 1 0 0 0 −4 1 0 0 1 0 0 0 
 
 0 1 0 0 1 −4 1 0 0 1 0 0  nl blocs
 
L2D =  0 0 1 0 0 1 −4 1 0 0 1 0  de taille
 
 0 0 0 1 0 0 1 −4 0 0 0 1  nc × nc
 
 
 0 0 0 0 1 0 0 0 −4 1 0 0 
 
 0 0 0 0 0 1 0 0 1 −4 1 0 
 
 0 0 0 0 0 0 1 0 0 1 −4 1 
 
 0 0 0 0 0 0 0 1 0 0 1 −4 

←−−nc−−→ ←−−nc−−→ ←−−nc−−→ y

(UPMC) Thème 4 : Résolution numérique des EDP 2 et 3 mars 2017 22 / 57

 Un problème statique : l’équation de Laplace Reformulation vectorielle de la solution

L2D tridiagonale par blocs : exemple nl = 5 et nc = 4

 −4 1 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
x

1 −4 1 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
 0 1 −4 1 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 
 0 0 1 −4 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0

 
 1 0 0 0 −4 1 0 0 1 0 0 0 0 0 0 0 0 0 0 0 
 0 1 0 0 1 −4 1 0 0 1 0 0 0 0 0 0 0 0 0 0 
 0 0 1 0 0 1 −4 1 0 0 1 0 0 0 0 0 0 0 0 0

 
 0 0 0 1 0 0 1 −4 0 0 0 1 0 0 0 0 0 0 0 0 
0 0 0 0 1 0 0 0 −4 1 0 0 1 0 0 0 0 0 0 0
 
 
 0 0 0 0 0 1 0 0 1 −4 1 0 0 1 0 0 0 0 0 0  nl blocs
L2D =  0 0 0 0 0 0 1 0 0 1 −4 1 0 0 1 0 0 0 0 0  de taille
0 0 0 0 0 0 0 1 0 0 1 −4 0 0 0 1 0 0 0 0
  nc × nc
 
 0 0 0 0 0 0 0 0 1 0 0 0 −4 1 0 0 1 0 0 0 
 0 0 0 0 0 0 0 0 0 1 0 0 1 −4 1 0 0 1 0 0 
0 0 0 0 0 0 0 0 0 0 1 0 0 1 −4 1 0 0 1 0
 
 
 0 0 0 0 0 0 0 0 0 0 0 1 0 0 1 −4 0 0 0 1 
 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 −4 1 0 0 
 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 1 −4 1 0

 
 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 1 −4 1 
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 1 −4 
y
←−−nc−−→ ←−−nc−−→ ←−−nc−−→ ←−−nc−−→ ←−−nc−−→

nl = 5 et nc = 4

(UPMC) Thème 4 : Résolution numérique des EDP 2 et 3 mars 2017 23 / 57

Sens de parcours et influence des bords (nl = 5 nc = 4)

 -4 1 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 x


 1 -4 1 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
 
 0 1 -4 1 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 
 
k =1  0 0 1 -4 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 
(1,1)  
 1 0 0 0 -4 1 0 0 1 0 0 0 0 0 0 0 0 0 0 0 
 
 0 1 0 0 1 -4 1 0 0 1 0 0 0 0 0 0 0 0 0 0 
 

nl blocs de taille nc × nc
 0 0 1 0 0 1 -4 1 0 0 1 0 0 0 0 0 0 0 0 0 
 0 0 0 1 0 0 1 -4 0 0 0 1 0 0 0 0 0 0 0 0 
 
 
 0 0 0 0 1 0 0 0 -4 1 0 0 1 0 0 0 0 0 0 0 
 
 0 0 0 0 0 1 0 0 1 -4 1 0 0 1 0 0 0 0 0 0 
 
 0 0 0 0 0 0 1 0 0 1 -4 1 0 0 1 0 0 0 0 0 
 0 0 0 0 0 0 0 1 0 0 1 -4 0 0 0 1 0 0 0 0 
 
 
 0 0 0 0 0 0 0 0 1 0 0 0 -4 1 0 0 1 0 0 0 
 
 0 0 0 0 0 0 0 0 0 1 0 0 1 -4 1 0 0 1 0 0 
 
 0 0 0 0 0 0 0 0 0 0 1 0 0 1 -4 1 0 0 1 0 
 0 0 0 0 0 0 0 0 0 0 0 1 0 0 1 -4 0 0 0 1 
 
 
 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 -4 1 0 0 
 
 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 1 -4 1 0 
 
 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 1 -4 1 
 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 1 -4 

y
←−nc−→ ←−nc−→ ←−nc−→ ←−nc−→ ← −nc−→
 Un problème statique : l’équation de Laplace Reformulation vectorielle de la solution

Le vecteur B résultant des CL

B est un vecteur creux à N = nl × nc composantes formé à partir des
2(nl + nc) températures aux bords. On peut le décomposer sous la forme :

 Nord   Sud   Ouest   Est 

T1N 0 T1W 0
   ..     .. 
 T2N   .   0   . 
 .   ..   .   
 .     .   
 .   .   .   0 
 Tnc N   0   0   T1E 
       
 0   0   T2W   0 
 .   ..
    . 
 .  . 
B=
    
 . .  0  . 
−  −  −   − 
   
 .. ..  ..
  
     
 .   .   .   0 
 0   0   0   T2E 
       
 0   T1S   Tnl W   0 
 .       . 
 .       . 
 .   T2S   0   . 
 .   .   .   
 .   .   .   
. . . 0
0 Tnc S 0 Tnl E

(UPMC) Thème 4 : Résolution numérique des EDP 2 et 3 mars 2017 25 / 57

 Un problème statique : l’équation de Laplace Solution du problème

Solution du système linéaire

Plusieurs méthodes permettent de résoudre le système linéaire (16),

L2D Z = B, par exemple :
1 Les méthodes les plus générales de résolution numérique des systèmes
linéaires, telles que la méthode d’élimination de Gauss-Jordan, la
décomposition « LU » 1 ;
2 Des méthodes adaptées au cas où A est une matrice « creuse », bande,
symétrique...
3 Des méthodes de relaxation qui s’apparentent, dans le cas d’une
équation telle que l’équation de Laplace, à la résolution temporelle de
l’équation de la chaleur.
4 L’utilisation des transformées de Fourier.

1. La factorisation « LU » exprime la matrice comme produit d’une matrice

triangulaire inférieure (Lower) par une matrice triangulaire supérieure (Upper).
(UPMC) Thème 4 : Résolution numérique des EDP 2 et 3 mars 2017 26 / 57
 Un problème statique : l’équation de Laplace Solution du problème

Résolution du système par méthode de relaxation

Méthode de relaxation : méthode itérative, basée sur la

décomposition en parties diagonale/non-diagonale : L2D = D + N
L2D Z = B équivaut à DZ = −(NZ − B).
Solution = point fixe de l’application Z 7→ −D−1 (NZ − B)
Calcul par récurrence en partant d’une valeur quelconque Z0 :

Zn+1 = −D−1 N Zn − B = N Zn − B /4

Convergence : Zn+1 − Z = N4 (Zn − Z) avec les valeurs absolues des valeurs

propres de N toutes inférieures à 4 (rayon spectral de N < 4).
Méthode aisée à coder et bien adaptée aux géométries complexes.

(UPMC) Thème 4 : Résolution numérique des EDP 2 et 3 mars 2017 27 / 57

 Un problème statique : l’équation de Laplace Solution du problème

Résolution du système linéaire à l’aide de Lapack

Méthodes générales (Gauss-Jordan, LU) ∝ N 3 opérations.

Ici : au plus 5 coefficients non nuls par ligne
=⇒ exploitation du caractère «creux» de la matrice : N opérations.
utilisation de Blas/Lapack :
Matrices générales : xGESV
−L2D définie positive : xPOSV
Symétrique : xSYSV
Matrice bande : xSBSV (amélioration significative)
Toutes ces propriétés : xPBSV

(UPMC) Thème 4 : Résolution numérique des EDP 2 et 3 mars 2017 28 / 57

 Un problème statique : l’équation de Laplace Solution du problème

Stockage bande
Fonctions pour matrices bandes =⇒ stockage spécifique de la matrice :
«stockage bande».
Ne fournir que les sur- et sous-diagonales significatives à la procédure.
Deux possibilités :
Écriture directe de la matrice contenant les diagonales
Écriture de la matrice complète puis « conversion ».
(en C : procédure matfull2band fournie)

ku=1 Matrice bande n=6 Stockage bande (kl+ku+1)×n

♣

kl=2 ♣

♦ ♥

♦ ♠

(UPMC) Thème 4 : Résolution numérique des EDP 2 et 3 mars 2017 29 / 57

 Problèmes dynamiques

Sommaire

1 Introduction
Généralités
Conditions aux limites
Problèmes statiques et dynamiques
2 Exemple de problème statique linéaire : l’équation de Laplace
Discrétisation du problème
Reformulation vectorielle de la solution
Solution du problème
3 Problèmes dynamiques
Discrétisation de l’espace et du temps
Algorithme explicite du premier ordre en temps
Algorithme implicite du premier ordre en temps
Algorithme de Crank–Nicolson
4 Équation des ondes

(UPMC) Thème 4 : Résolution numérique des EDP 2 et 3 mars 2017 30 / 57

 Problèmes dynamiques

Introduction aux problèmes dynamiques

Problèmes dynamiques :
Évolution temporelle à partir d’un état d’équilibre
Phénomènes de propagation

Caractéristiques :
Conditions initiales
Conditions aux limites (peuvent dépendre du temps)
=⇒ système forcé : transitoire puis régime permanent ν ≈ νforcé

Exemple : équation de la chaleur à une dimension spatiale pour simplifier.

(UPMC) Thème 4 : Résolution numérique des EDP 2 et 3 mars 2017 31 / 57

 Problèmes dynamiques

Équation de diffusion 1-D

Équation de diffusion à 1-D avec D = cste > 0 :

∂u ∂2u
=D 2 (17)
∂t ∂x

domaines spatial x ∈ [0, L], et temporel t ∈ [0, T ]

Conditions Initiales (CI) : u(x, 0) = u0 (x) ∀x

Conditions aux Limites (CL) : u(0, t) = ug (t) et u(L, t) = ud (t) ∀t

N.B. : espace et temps ont un rôle dissymétrique :

CL/CI
Degré différent des dérivées partielles

(UPMC) Thème 4 : Résolution numérique des EDP 2 et 3 mars 2017 32 / 57

 Problèmes dynamiques Discrétisation de l’espace et du temps

Discrétisation du domaine [O, L] × [0, T ]

t
L = (nx + 1)δx
n = nt + 1

T = (nt + 1)δt
Un
point intérieur
δt nx points intérieurs
unj inconnu
Cond. aux Limites
en x = 0 et x = L n=1
Cond. Initiale
n=0 x
en t = 0 (n = 0)
j=0 j=1 δx j = nx + 1

Figure 6 – Discrétisation du domaine [0, L] × [0, T ]

NB : visualisation ultérieure avec les 2 bords
⇒ à chaque temps, écrire Un (nx points)
plus les CL (2 points) : nx + 2 points

(UPMC) Thème 4 : Résolution numérique des EDP 2 et 3 mars 2017 33 / 57

 Problèmes dynamiques Discrétisation de l’espace et du temps

Discrétisation du domaine

Discrétiser les deux dimensions :

Espace : [0, L] est divisé en nx + 1 pas de largeur δx = nxL+1 .
nx points intérieurs + 2 bords
Temps : [0, T ] est subdivisé en nt + 1 pas de largeur δt = ntT+1 .

unj = u(x, t) pour x = j δx et t = n δt

CL =⇒ un0 ≡ ung et unnx +1 ≡ und pour 1 6 n 6 nt + 1

CI =⇒ u0j pour 1 6 j 6 nx .

Méthode de résolution fondamentalement itérative sur le temps :

unj + CL =⇒ un+1
j à l’intérieur du domaine.

(UPMC) Thème 4 : Résolution numérique des EDP 2 et 3 mars 2017 34 / 57

 Problèmes dynamiques Discrétisation de l’espace et du temps

Discrétisation des opérateurs : différences finies

Laplacien à l’ordre 2 (cf régimes stationnaires) :
unj+1 − 2unj + unj−1
∆u ≈
(δx)2
Dérivée ordre 1 par rapport au temps — au choix comme pour les EDO :
un+1 − unj
n
∂u

=
j
vers l’avant : (cf Euler forward)
∂t j δt
unj − un−1
n
∂u

=
j
vers l’arrière : (cf Euler backward)
∂t j δt
un+1 − un−1
n
∂u

=
j j
ou centrée :
∂t j 2δt

Ce choix détermine le caractère explicite ou implicite de la méthode.

(UPMC) Thème 4 : Résolution numérique des EDP 2 et 3 mars 2017 35 / 57

 Problèmes dynamiques Algorithme explicite du premier ordre en temps

Algorithme explicite du premier ordre en temps

Approche simpliste : dérivée vers l’avant

un+1 − unj unj+1 − 2unj + unj−1

 
=D pour 1 6 j 6 nx
j
(18)
δt (δx)2

Semblable à méthode d’Euler explicite (progressive).

Paramètre caractéristique sans dimension :

D δt
α= (19)
δx2

(UPMC) Thème 4 : Résolution numérique des EDP 2 et 3 mars 2017 36 / 57

 Problèmes dynamiques Algorithme explicite du premier ordre en temps

Système d’équations linéaires (algorithme explicite)

On peut alors écrire :

un+1
j = α unj+1 + (1 − 2α) unj + α unj−1 pour 1 6 j 6 nx (20)

Méthode explicite : un+1

j obtenus directement en fonctions des unk .

Pour j = 1 ou nx : unj−1 , unj+1 sur les bords

=⇒ fixés par conditions aux limites

Note : équivalent à méthode de relaxation si α = 1/2.

(UPMC) Thème 4 : Résolution numérique des EDP 2 et 3 mars 2017 37 / 57

 Problèmes dynamiques Algorithme explicite du premier ordre en temps

Schéma de progression explicite

t
Un+1
(n + 1)δt
Figure 7 – Schéma
illustrant la progression nδt
d’un pas de temps avec Un
l’algorithme explicite du
premier ordre en δt.

0 x
0 j−1 j j+1 nx + 1

un+1
j = α unj+1 + (1 − 2α) unj + α unj−1

(UPMC) Thème 4 : Résolution numérique des EDP 2 et 3 mars 2017 38 / 57

 Problèmes dynamiques Algorithme explicite du premier ordre en temps

Explicite : Formulation matricielle

Vecteur Un = (un1 , . . . , unnx ) des solutions à t = nδt
pour les nx points intérieurs
Le système d’équations linéaires prend la forme :
Un+1 = M(−α) Un + αVn ,
Vn = (ung , 0, . . . , 0, und ) : vecteur des conditions aux limites,
M(−α) = 11 + α L1D
où L1D = matrice tridiagonale représentant la dérivée seconde 1-D
   
−2 1 0 ··· 0 1−2α α 0··· 0
. .. .. .. 
1 .. .
  
 1 −2 .  α 1−2α α . 
. .. .. .. .. ..
   
L1D =
 0 .. . .

0  −→ M(−α) = 
 0 . . . 0 
 (21)
 .. . . ..
  .. .. ..

. . −2 . . 1−2α α
   
. 1 .
0 ··· 0 1 −2 0 ··· 0 α 1−2α

(UPMC) Thème 4 : Résolution numérique des EDP 2 et 3 mars 2017 39 / 57

 Problèmes dynamiques Algorithme explicite du premier ordre en temps

Explicite : stabilité
Problème important : instabilité de la méthode explicite
pour un pas de temps « un peu trop grand » (α > 0.5) :
0.3
1pas
2pas
0.25 4pas
6pas
0.2 α=1.2 8pas
10pas
0.15 12pas
14pas
0.1

0.05

-0.05

-0.1

-0.15
0 20 40 60 80 100 120 140

Figure 8 – Résolution explicite avec α = 1.2 > 0.5.

(UPMC) Thème 4 : Résolution numérique des EDP 2 et 3 mars 2017 40 / 57

 Problèmes dynamiques Algorithme explicite du premier ordre en temps

Explicite : analyse de stabilité de Fourier–von Neumann

Stabilité des modes propres :
u(x, t) = A(t) exp(ikx) ⇒ unj = A(n) eik j δx , (22)
dA/ dt = −Dk 2 A(t) A(t) = A(0) exp −Dk 2 t


solution si soit
L’algorithme explicite donne la version discrète de A(t) :

= 1 + α eik δx − 2 + e−ik δx
 
A(n+1)
A(n)
2
= 1 + α eik δx/2 − e−ik δx/2

(23)

= 1 − 4α sin2
 
k δx
2 .

→ suite géométrique A(n) = A(0) θn , de raison réelle θ = 1 − 4α sin2

 
k δx
2
L’algorithme explicite est stable si |θ| 6 1 pour tous les modes, soit

0 6 α 6 1/2 c-à-d δt 6 δx2 /2D

(UPMC) Thème 4 : Résolution numérique des EDP 2 et 3 mars 2017 41 / 57

 Problèmes dynamiques Algorithme explicite du premier ordre en temps

Stabilité de l’explicite : lien avec le spectre de L1D

La p-ième valeur propre de la matrice L1D , en dimension n, s’écrit en
fonction de β = π/2(n + 1) = πδx/2L :

− 4 sin2 pβ ; m = sin m p β
vecteur propre : u(p)



valeur propre :

On en déduit celles de la matrice M(−α) :

M(−α) = 11 + α L1D −→ 1 − 4α sin2 p πδx

2L .

Stabilité ⇐⇒ valeurs propres de M de module inférieur à 1 : ρ(M) < 1.

La relation de dispersion, iω = Dk 2 déduite de l’équation de la diffusion

est remplacée par :

2i sin ωδt
exp(i ωδt = 4α sin2 kδx




2 2 2 (24)

qui n’a de solution ω à partie imaginaire négative (stable) que si 4α 6 2.

Avant d’aller plus loin dans l’analyse numérique, voici quelques exemples.
(UPMC) Thème 4 : Résolution numérique des EDP 2 et 3 mars 2017 42 / 57
 Problèmes dynamiques Algorithme explicite du premier ordre en temps

Exemple : Diffusion d’une distribution en « marche »

1.0 Diffusion d'une marche

t =0

0.5 t→∞
u(x,t)

0.0

0.5

1.0
0 20 40 60 80 100 120
position : x

Figure 9 – Diffusion d’une marche avec des CL fixes +1 à gauche et −1 à

droite : u(x, t) tend vers une droite si t → ∞
(UPMC) Thème 4 : Résolution numérique des EDP 2 et 3 mars 2017 43 / 57
 Problèmes dynamiques Algorithme explicite du premier ordre en temps

Exemple : Diffusion d’une superposition de 2 modes

2.4 Rôle de la longueur d'onde sur la diffusion

2.2

2.0

1.8
u(x,t)

1.6

1.4

1.2

1.0

0.8
0 20 40 60 80 100 120
position : x

Figure 10 – Amortissement du mode 5 plus rapide que celui du mode 1

mais pour t → ∞, on obtiendrait une droite.
(UPMC) Thème 4 : Résolution numérique des EDP 2 et 3 mars 2017 44 / 57
 Problèmes dynamiques Algorithme implicite du premier ordre en temps

Solution au problème de stabilité : méthode rétrograde

Algorithme explicite (progressif) simple, mais instable sauf très petits pas
Stabiliser (comme Euler en EDO) en passant en implicite (rétrograde).
Algorithme implicite : dérivée temporelle « vers l’arrière ».

un+1 j+1 − 2uj

un+1 + uj−1
n+1 n+1
− unj
" #
=D
j
(25)
δt (δx)2

=⇒ système de nx équations linéaires à inverser :

j−1 + (1 + 2α) uj
−α un+1 j+1 = uj pour j ∈ [1, nx ]
n+1
− α un+1 n
(26)

sauf pour j = 1 et j = nx où interviennent les CL.

(UPMC) Thème 4 : Résolution numérique des EDP 2 et 3 mars 2017 45 / 57

 Problèmes dynamiques Algorithme implicite du premier ordre en temps

Schéma de progression implicite

ujn+1
Figure 11 – Schéma (n+1) δt
illustrant la
progression d’un pas n δt ujn
de temps avec
l’algorithme implicite
du premier ordre en δt.
x
j δx (nx+1) δx
0
(j −1) δx (j +1) δx

j−1 + (1 + 2α) uj
−α un+1 j+1 = uj
n+1
− α un+1 n

(UPMC) Thème 4 : Résolution numérique des EDP 2 et 3 mars 2017 46 / 57

 Problèmes dynamiques Algorithme implicite du premier ordre en temps

Formulation matricielle

De façon matricielle : M(α) Un+1 = Un + α Vn+1

 
1+2α -α 0 · · · 0
. . . . 
 -α . . . . . . .. 

. . .
 
M(α) = 11 − α L1D  0 .. .. .. 0 
=  (27)
 . . . . 
 . .. .. .. 
. -α
0 · · · 0 -α 1+2α

M tridiagonale mais M−1 matrice dense =⇒ un+1 j dépend de tous les

autres éléments, mais d’autant moins lorsqu’ils sont éloignés.
Lorsque α → 0, M(α)−1 → M(−α) : on retrouve l’algorithme explicite.

(UPMC) Thème 4 : Résolution numérique des EDP 2 et 3 mars 2017 47 / 57

 Problèmes dynamiques Algorithme implicite du premier ordre en temps

Stabilité de l’algorithme implicite

Analyse Fourier — Von Neumann : modes propres ∝ A(t) sin(kx).

Suite géométrique de raison :
1
θ=  61
1 + 4 α sin2

k δx
2

On obtient plus simplement ce résultat à partir des valeurs propres de

M(α)−1 = (11 − αL1D )−1 , inverses de 1 + 4α sin2 p 2L
πδx

.
Algorithme implicite inconditionnellement stable quelque soit δx, δt.
Toutefois, approximation du premier ordre en δt : description précise de la
solution transitoire nécessite un pas de temps petit.

(UPMC) Thème 4 : Résolution numérique des EDP 2 et 3 mars 2017 48 / 57

 Problèmes dynamiques Algorithme de Crank–Nicolson

Principe de l’algorithme de Crank–Nicolson

Objectif : accroître la précision =⇒ améliorer l’approximation en

∂
différences finies de ∂t .

Même stratégie que « Euler modifié » (cf. EDO) : moyenne des deux
algorithmes.

un+1 D (uj+1 − 2uj + uj−1 ) + (unj+1 − 2unj + unj−1 )

n+1 n+1 n+1
− unj
" #
=
j
δt 2 (δx)2

NB : différence finie temporelle inchangée, mais vue comme schéma centré

à l’instant milieu (n + 12 )δt, grâce à l’évaluation de la moyenne des
laplaciens discrets aux deux instants.

Algorithme de Crank–Nicolson du deuxième ordre en δt et δx.

(UPMC) Thème 4 : Résolution numérique des EDP 2 et 3 mars 2017 49 / 57

 Problèmes dynamiques Algorithme de Crank–Nicolson

Schéma de progression de Crank–Nicolson

ujn+1
Figure 12 – Schéma (n+1) δt
illustrant la progression n δt ujn
d’un pas de temps avec
l’algorithme de
Crank–Nicolson. x
j δx (nx+1) δx
0
(j −1) δx (j +1) δx

α α α α
− un+1 + (1 + α)un+1 − un+1 = unj−1 + (1 − α)unj + unj+1
2 j−1 j
2 j+1 2 2

(UPMC) Thème 4 : Résolution numérique des EDP 2 et 3 mars 2017 50 / 57

 Problèmes dynamiques Algorithme de Crank–Nicolson

Stabilité de l’algorithme de Crank–Nicolson

Etude de stabilité Von Neumann : décroissance géométrique des modes :

1 − 2 α sin2
 
k δx
θ=
2
(28)
1 + 2 α sin2
 
k δx
2

ce qui permet de montrer que cet algorithme est lui aussi

inconditionnellement stable (i.e. quel que soit le choix de δt).

θCN (α) = θexplicite ( α2 ) × θimplicite ( α2 )

(UPMC) Thème 4 : Résolution numérique des EDP 2 et 3 mars 2017 51 / 57

 Problèmes dynamiques Algorithme de Crank–Nicolson

Formulation matricielle

Système d’équations linéaires à résoudre (j ∈ [1, nx ]) :

α α α n α
j−1 + (1 + α)uj
− un+1 j+1 = uj−1 + (1 − α)unj + unj+1
n+1
− un+1
2 2 2 2
(29)
sauf pour j = 1 et j = nx où interviennent les CL.
De façon matricielle :

α −α αh n
   
=M Un + V + Vn+1
i
n+1
M U (30)
2 2 2

où M(α) = 11 − α L1D

(UPMC) Thème 4 : Résolution numérique des EDP 2 et 3 mars 2017 52 / 57

 Équation des ondes

Sommaire

1 Introduction
Généralités
Conditions aux limites
Problèmes statiques et dynamiques
2 Exemple de problème statique linéaire : l’équation de Laplace
Discrétisation du problème
Reformulation vectorielle de la solution
Solution du problème
3 Problèmes dynamiques
Discrétisation de l’espace et du temps
Algorithme explicite du premier ordre en temps
Algorithme implicite du premier ordre en temps
Algorithme de Crank–Nicolson
4 Équation des ondes

(UPMC) Thème 4 : Résolution numérique des EDP 2 et 3 mars 2017 53 / 57

 Équation des ondes

Équation des ondes : différences finies

Résolution numérique de l’équation de d’Alembert :

Discrétisation du domaine comme pour l’équation de la chaleur ;
CI arbitraires u(x, 0) ;
CL Dirichlet ou Neumann.
schéma de différences finies «naturel» :

∂2u 2
2∂ u
un+1 + un−1 − 2unj un + unj−1 − 2unj
= =
j j 2 j+1
c ⇒ c
∂t2 ∂x2 (δt)2 (δx)2
(31)
(deuxième ordre en temps et espace).
Note : symétrie entre temps et espace.

(UPMC) Thème 4 : Résolution numérique des EDP 2 et 3 mars 2017 54 / 57

 Équation des ondes

Algorithme explicite
δt
Poser α = c ⇒ équations explicites à trois niveaux de temps :
δx
un+1
j = α2 (unj+1 + unj−1 ) + 2(1 − α2 )unj − un−1
j . (32)

soit, sous forme matricielle :

Un+1 = 211 + α2 L1D Un − Un−1 + α2 L1D Vn

où V représente à nouveau les conditions aux limites sur les bords.

Ce schéma est appelé “leap-frog” en anglais, singulièrement traduit en

« saute-mouton ».

Ordre deux en temps ⇒ 2 CI : u(x, 0) et ∂u/∂t(x, 0),

ou 2 temps : u(x, 0) et u(x, δt).

(UPMC) Thème 4 : Résolution numérique des EDP 2 et 3 mars 2017 55 / 57

 Équation des ondes

Critère CFL

Relation de dispersion associée au schéma : sin2 ω δt

= α2 sin2 k δx



2 2
α 6 1 : solutions réelles (ie non divergentes).
α = 1 : relation de dispersion exacte (et unj disparaît).

⇒ Critère de Courant–Friedrichs–Lewy (CFL) :

α 6 1 ⇐⇒ méthode stable.
Valeur précise dépend de la méthode de discrétisation, dimension de
l’espace, éventuellement direction de propagation.
3 formulations équivalentes :
pas de temps δt < δx/c temps de propagation de l’onde ;
pas d’espace δx > c δt distance parcourue par l’onde ;
vitesse «diffusion information» δx/δt > c célérité de l’onde.

(UPMC) Thème 4 : Résolution numérique des EDP 2 et 3 mars 2017 56 / 57

 Équation des ondes

Schéma de progression et critère CFL

Le schéma explicite « saute-mouton » est désormais le suivant :

t

Figure 13 – Progression
d’un pas de temps dans un+1
j
l’algorithme explicite du (n + 1)δt
second ordre en temps et cδt
espace pour l’équation nδt
des ondes de d’Alembert
(n − 1)δt
un−1
Critère CFL : cδt < δx 1/c j δt

0 δx L x
Le critère CFL se traduit (dans ce cas simple 1D) par le fait que des points
intervenant dans la progression doivent être hors du cône de lumière, pour
laisser à l’onde « le temps de se propager ».

(UPMC) Thème 4 : Résolution numérique des EDP 2 et 3 mars 2017 57 / 57