Sommaire

1 Introduction
Généralités
Thème 4: résolution numérique des EDP linéaires Conditions aux limites
Problèmes statiques et dynamiques 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
Université Pierre et Marie CURIE Solution du problème
MNCS
Méthodes numériques et calcul scientifique 3 Problèmes dynamiques
Discrétisation de l’espace et du temps
2 et 3 mars 2017 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 1 / 60 (UPMC) Thème 4 : Résolution numérique des EDP 2 et 3 mars 2017 2 / 60

Introduction Introduction

Contexte Exemples d’EDP scalaires les plus connues

Équations différentielles ordinaires (EDO) : une variable indépendante 1 L’équation des ondes ou « de d’Alembert »
Équations aux dérivées partielles (EDP) : plusieurs variables ∂2u
− c2 ∆ u = 0 (1)
indépendantes ∂t2
Les équations aux dérivées partielles (EDP) interviennent dans la 2 L’équation de Schrödinger
description de très nombreux problèmes de physique, chimie, sciences de la ∂u 1
i = − ∆u + V u (2)
Terre, biologie : mécanique des fluides, propagation des ondes, ∂t 2
électromagnétisme, phénomènes de diffusion... 3 L’équation de diffusion
Nous nous limitons ici à des EDP linéaires relatives à un champ scalaire ∂u # –
u( #–
r , t). Les généralisations vectorielles sont bien sûr possibles, et = div[D grad u] = D ∆ u si D est uniforme (3)
∂t
techniquement analogues. 4 Les équations de Poisson et de Laplace
Résolution numérique par la méthode des différences finies avec un
∆u = ρ ∆u = 0 (4)
maillage régulier.
Mais il existe d’autres méthodes, éléments finis par exemple, mieux B Les équations de Navier-Stokes en mécanique des fluides sont en
adaptées dans le cas de domaines à géométrie complexe. 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 3 / 60 (UPMC) Thème 4 : Résolution numérique des EDP 2 et 3 mars 2017 4 / 60
 Introduction Introduction

Équation des ondes Équation de Schrödinger

L’équation des ondes ou « de d’Alembert » :

L’équation de Schrödinger :
∂2u
− c2 ∆ u = 0 (5) ∂u 1
∂t2 i = − ∆u + V u (7)
∂t 2
où u = u( #–
r , t) est un champ scalaire fonction de la position #–
r et du
temps t, comme la pression d’un gaz dans le cas du son. qui décrit en mécanique quantique la fonction d’onde d’une particule
La constante c est la vitesse de propagation des ondes et ∆ le Laplacien massive dans un potentiel V .
scalaire : En remplaçant t par z, elle représente aussi la propagation selon z d’une
∂2 ∂2 onde harmonique à la limite paraxiale, dans un milieu de permittivité
∆ = + dans le plan en coordonnées cartésiennes (6)  = n2 = V + 1.
∂x2 ∂y 2

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

Introduction Introduction

Équation de la diffusion Équations de Poisson et de Laplace

L’équation de diffusion
Des cas particuliers sont les équations de Poisson et de Laplace, qui
∂u #– #– # – peuvent être vues comme des cas limites lorsque la dépendance en temps
+ div J = 0 avec J = −D grad u
∂t est supprimée, ou remplacée par une dépendance harmonique.
#– L’équation de Poisson
Ici J est le flux (de u) par unité de surface et D le coefficient de diffusion ∆u = ρ
(m2 s−1 ), ce qui donne, dans le cas où D est homogène : est en électrostatique celle qui régit le potentiel Φ( #–
r ) créé par la densité
#–
de charges ρ( r ).
∂u
= D ∆u (8) L’équation de Laplace est un cas particulier où la charge est nulle :
∂t
Cette équation décrit : ∆u = 0 (9)
les transferts de chaleur par conduction dans un milieu continu ; elle
On va notamment la retrouver comme solution stationnaire de l’équation
est alors appelée équation de la chaleur et u ≡ T et D ≡ λ/Cp .
de la chaleur.
la diffusion des particules d’un fluide dans un autre (équation de Fick).

(UPMC) Thème 4 : Résolution numérique des EDP 2 et 3 mars 2017 7 / 60 (UPMC) Thème 4 : Résolution numérique des EDP 2 et 3 mars 2017 8 / 60
 Introduction Généralités Introduction Conditions aux limites

Équations elliptiques/paraboliques/hyperboliques Conditions aux limites

L’étude théorique des EDP est un vaste domaine en mathématiques qui On doit se donner des « conditions au limites » (CL). L’équation étant
est hors du cadre de ce cours. vérifiée dans un domaine D de l’espace (ou espace-temps), on distingue
Les EDP linéaires du second ordre évoquées plus haut sont classées selon des conditions de 2 types :
la forme des coefficients figurant devant les dérivées partielles : Conditions de Dirichlet : on impose la valeur de u sur la bordure de
D. Dans le cas où on étudie un problème dépendant du temps, cela
∂2u ∂2u ∂2u
α 2
+β + γ 2 = ... inclut des conditions « initiales »
∂x ∂x∂y ∂y Conditions de Neumann : c’est la valeur de la dérivée normale
# –
∂u/∂n = grad u · #–n que l’on impose.
β 2 − 4αγ > 0 hyperboliques : équation des ondes
Robin : plus général, combinaison linéaire des deux.
β 2 − 4αγ = 0 paraboliques : diffusion ou Schrödinger
Dans le cas général, on a une combinaison des deux selon différentes
β 2 − 4αγ < 0 elliptiques : équations de Laplace, Poisson, Helmholtz
parties de la bordure.
Elles ont des propriétés mathématiques différentes, et leur résolution par À titre d’exemple pour l’équation de la chaleur :
des techniques numériques requiert parfois des méthodes distinctes. Des CL de Dirichlet correspondent à des parois isothermes, qui
Dans la pratique, nous nous concentrerons sur une méthode de résolution imposent leur température.
# –
approchée par des méthodes numériques, dans le cas principalement de Des CL de Neumann avec grad u · #– n = 0 annulent le transfert de
l’équation de la chaleur. chaleur : parois adiabatiques.
(UPMC) Thème 4 : Résolution numérique des EDP 2 et 3 mars 2017 9 / 60 (UPMC) Thème 4 : Résolution numérique des EDP 2 et 3 mars 2017 10 / 60

Introduction Problèmes statiques et dynamiques Introduction Problèmes statiques et dynamiques

Problèmes statiques et dynamiques Problèmes statiques

Selon que la variable temps est ou non présente, deux classes de problèmes Nous nous intéresserons en premier lieu aux problèmes statiques.
doivent être distinguées : Il s’agit de problèmes où l’on cherche la solution d’une EDP où le temps
Problèmes dynamiques : il s’agit de déterminer u( #–
r , t) dans un n’intervient pas. C’est par exemple la solution d’équilibre de l’équation de
domaine D, sur un certain intervalle de temps, en fonction des CL diffusion qui sera asymptotiquement atteinte quand t → ∞, et qui doit
spatio-temporelles : donc vérifier l’équation de Laplace.
CI (conditions initiales) : u( #–
r , t = 0) = · · · Dans ce cas u représente un champ de température T (x, y) (avec une
CL (conditions aux limites) : u( #– r = ·, t) = · · · origine arbitraire) et on a dans le plan en coordonnées cartésiennes :
Problèmes statiques : le temps n’intervient pas, il s’agit de déterminer
∂2T ∂2T
u( #–
r ) en fonction des seules CL ; cela correspond en général à un ∆T = + =0. (10)
∂x2 ∂y 2
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 11 / 60 (UPMC) Thème 4 : Résolution numérique des EDP 2 et 3 mars 2017 12 / 60
 Introduction Problèmes statiques et dynamiques Un problème statique : l’équation de Laplace

Problèmes dynamiques Résolution de l’équation de Laplace dans un rectangle

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

Puis, dans une deuxième partie, nous considérerons les problèmes

dynamiques. Nous cherchons une solu- A l'intérieur T(x,y) inconnue
vérifiant ∆T=0
Deux méthodes simples seront exposées et analysées sur l’exemple de tion numérique approchée de
∆x
l’équation de diffusion à 1-D : ∆ T = 0 à l’intérieur d’un
domaine D que pour simpli-
∂u ∂2u
=D (11) fier nous supposerons rectan-
∂t ∂x2 gulaire, de côtés ∆x et ∆y ∆y
avec diverses conditions initiales. (voir figure 1). 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 13 / 60 (UPMC) Thème 4 : Résolution numérique des EDP 2 et 3 mars 2017 14 / 60

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

Exemple de carte d’isothermes Discrétisation du domaine

Nous supposons en outre que les parois ont des températures imposées de
l’extérieur (CL de Dirichlet).
˚
10C y
0 5 10 15 20 25
0
Figure 2 – Exemple de Nous procédons à une discrétisation du domaine D : la température sera
solution de l’équation
calculée en un nombre fini de points. Pour simplifier, le maillage choisi est
14 13
12 11 ∆ T = 0, dans le cas où une
une grille de pas dx = dy = p comme représenté sur les figures 3 et 4.
10
5 16
18
9
8
paroi est maintenue à 20◦ C,
20C
˚

7 une autre adjacente est à On a nl = ∆x/p − 1 lignes et nc = ∆y/p − 1 colonnes.

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 15 / 60 (UPMC) Thème 4 : Résolution numérique des EDP 2 et 3 mars 2017 16 / 60
 Un problème statique : l’équation de Laplace Discrétisation du problème Un problème statique : l’équation de Laplace Discrétisation du problème

Maillage du domaine D Indexation 2D du domaine

1 .'.'.' j .'.'.'.'.'. nc y
1 (0,j' ).'.'.'.'.'. À l’intérieur de la grille, on a :
0 .'.'. j nc+1
0 nl = ∆x/p − 1 lignes horizontales.

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

p y
nc = ∆y/p − 1 colonnes verticales
i (i,j' )

.'.'.
i (i' ,j' )
.'.'.'.'.'.

Chaque point est repéré par un couple d’indices (i, j) où le premier indice
(i,nc+1)
note le numéro de la ligne et le second celui de la colonne.

.'.'.'.'.'.
nl
x p On note Ti,j la température en ce point, de coordonnées
∆y'=(nc+1)' p
nl' +1 xi = i × p, yj = j × p.
x
Figure 3 –
Les inconnues du problème sont les n = nl × nc valeurs de Ti,j
Chaque point de la grille est repéré
par un couple (i, j) où i repère la
Figure 4 – Les bords du domaine où 1 6 i 6 nl, 1 6 j 6 nc.
(C.L.) sont numérotés : i = 0 pour le
Elles seront déduites des 2(nl + nc) valeurs de Ti,j imposées aux bords
ligne et j la colonne. L’intérieur du haut, i = nl + 1 pour le bas, j = 0 pour
domaine est numéroté de i = 1 à nl la gauche, j = nc + 1 pour la droite. i = 0 (nord) ou i = nl + 1 (sud)
et de j = 1 à nc. j = 0 (ouest) ou j = nc + 1 (est).

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

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

Différences finies Différences finies (suite)

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 : On note que cette approximation peut s’obtenir de façon alternative à
partir de :
∂T p2 ∂ 2 T p3 ∂ 3 T
T (x ± p) = T (x) ± p + +
 
± O p4 (12)
∂x 2 ∂x2 3! ∂x3 ∂T p T (x ± p) − T (x)
 
x± ≈± (16)
montre que : ∂x 2 p
∂2T qui conduit à :
T (x + p) + T (x − p) = 2T (x) + p2 (x) +
 
O p4 (13)
∂x2
∂2T 1 ∂T p ∂T p
    
d’où : (xi , yj ) ≈ xi + , yj − xi − , yj (17)
∂2T T (x + p) + T (x − p) − 2 T (x) ∂x2 p ∂x 2 ∂x 2
(x) = +
 
O p2 (14)
∂x2 p2 soit encore :
ce qui conduit à :
∂2T Ti+1,j + Ti−1,j − 2 Ti,j
(xi , yj ) = + O p2
 
∂2T Ti+1,j + Ti−1,j − 2 Ti,j (18)
(x ) = +
 
∂x2 p 2
i , yj O p2 (15)
∂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 19 / 60 (UPMC) Thème 4 : Résolution numérique des EDP 2 et 3 mars 2017 20 / 60
 Un problème statique : l’équation de Laplace Discrétisation du problème Un problème statique : l’équation de Laplace Reformulation vectorielle de la solution

Différences finies : Laplacien 2D Réindexation 1-D du domaine (en RowMajor)

On ré-indexe les points à l’intérieur de la grille par un seul indice k :
On obtient ainsi le Laplacien 2D discrétisé à l’ordre 2 :
k = (i − 1)nc + j (21)
Ti+1,j + Ti−1,j + Ti,j+1 + Ti,j−1 − 4 Ti,j
∆ Ti,j = +
 
O p2 (19) qui fait correspondre (cf fig. 5) au tableau 2-D des Ti,j à calculer, le
p2
vecteur à N = nl × nc composantes Zk .
L’équation de Laplace se traduit alors par un système linéaire de N 1N Nord j
équations reliant les valeurs de Ti,j aux 4 points voisins : y Figure 5 – Renumérotation des
1 2 nc points de grille pour représenter la
1W
Ti+1,j + Ti−1,j + Ti,j+1 + Ti,j−1 − 4 Ti,j = 0 (20) solution sous la forme d’un vecteur :
2nc

Ouest
nc+1 2E le domaine est ré-indexé ligne par
où 2 × (nl + nc) températures sont connues (celles sur les points aux bords ligne de la gauche vers la droite, puis

Est
indiqués par des croix sur la figure 4) et N = nl × nc sont à déterminer. i de haut en bas (comme une image
nl׀nc numérique) donc en RowMajor.
Système non homogène car certains points intérieurs ont des voisins sur les
parois avec des températures imposées qui vont constituer le 2d membre. Les voisins au Nord et au Sud sont
alors espacés de nc dans Z.
Noter que le pas p n’intervient plus. x 1S k=(i€–1)nc+j
Sud
(UPMC) Thème 4 : Résolution numérique des EDP 2 et 3 mars 2017 21 / 60 (UPMC) Thème 4 : Résolution numérique des EDP 2 et 3 mars 2017 22 / 60

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

Formulation vectorielle Matrice L2D

En plaçant à droite du signe égal les températures imposées sur les bords,
Pour un point sans contact avec les bords du domaine, l’équation (20)
on constitue, à partir des CL, le second membre du système qui prend
s’écrit :
alors la forme :
Zk−nc + Zk−1 − 4Zk + Zk+1 + Zk+nc = 0 (22)
Ne pas croire que le second membre du système linéaire est nul ! AZ=B (23)
Pour un point en contact avec un bord, cette équation fait intervenir des où A = L2D est une matrice carrée N × N comportant
températures imposées : par exemple celles du bord Nord si i = 1. nl × nl blocs carrés de dimension nc × nc.
L2D est symétrique tridiagonale par blocs :
T2N + Z1 − 4Z2 + Z3 + Z2+nc = 0 pour j = 2
1 blocs de type L1D avec −4 au lieu de −2 sur la diagonale ;
=⇒ passer ces termes dans le second membre du système 2 blocs de type identité sur la sous-diagonale et la sur-diagonale.
Les CL de Dirichlet des bords créent le second membre Chaque ligne et chaque colonne contient 5 termes non-nuls comme (22) le
Chaque changement de membre d’une température imposée annule prévoit, sauf celles de numéro i × nc et i × nc + 1, où apparaissent les 0,
un des coefficients des Zk0 . (bords latéraux), ainsi que 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 23 / 60 (UPMC) Thème 4 : Résolution numérique des EDP 2 et 3 mars 2017 24 / 60
 Un problème statique : l’équation de Laplace Reformulation vectorielle de la solution Un problème statique : l’équation de Laplace Reformulation vectorielle de la solution

Matrice L2D avec nl = 3 et nc = 4 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

 x 1 −4 1 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0
−4 1 0 0 1 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



1 −4 1 0 0 1 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 −4 1 0 0 1 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 −4 0 0 0 1
 
 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
   

 1 0 0 0 −4 1 0 0 1 0 0 0 

  0 0 0 0 1 0 0 0 −4 1 0 0 1 0 0 0 0 0 0 0

 
 0 1 0 0 1 −4 1 0 0 1 0 0  nl blocs
  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
L2D =  0 0 1 0 0 1 −4 1 0 0 1 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 1 0 0 1 −4 0 0 0 1  nc × nc
   
  0 0 0 0 0 0 0 0 0 1 0 0 1 −4 1 0 0 1 0 0 
 
 0 0 0 0 1 0 0 0 −4 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 1 0 0 1 −4 1 0 
   
   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 1 0 0 1 −4 1 
   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 1 −4 
 
 0 0 0 0 0 0 0 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−−→ 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 25 / 60 (UPMC) Thème 4 : Résolution numérique des EDP 2 et 3 mars 2017 26 / 60

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

Sens de parcours et influence des bords (nl = 5 nc = 4) Le vecteur B résultant des CL


-4 1 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
x B est un vecteur creux à N = nl × nc composantes formé à partir des
1 -4 1 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0  2(nl + nc) températures aux bords. On peut le décomposer sous la forme :


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 
 
 Nord   Sud   Ouest   Est 

1 0 0 0 -4 1 0 0 1 0 0 0 0 0 0 0 0 0 0 0 
 
 T1N 0 T1W 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 
  T2N   .   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 0 1 0 0 1 -4 1 0 0 1 0 0 0 0 0 0 
     
  Tnc N  T1E 

 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 0 0 1 0 0 1 -4 0 0 0 1 0 0 0 0 
  0   0   T2W   0 

0 0 0 0 0 0 0 0 1 0 0 0 -4 1 0 0 1 0 0 0 
  .   ..
    . 
 .
 −  .. 

B=
      

 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 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 0 1 0 0 0 -4 1 0 0 
   0  T2E 
   0   0 
0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 1 -4 1 0 
         
  0   T1S   Tnl W   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 
  .       . 
←−nc−→ ←−nc−→ ←−nc−→ ←−nc−→ ←−nc−→
y  .   T2S   0   . 
 .   .   .   
 .   .   .   
. . . 0
i = 2, j = 1 =⇒ k = 5 alias bk=(i−1)nc+j
0 Tnc S 0 Tnl E
= −δi,1 Ti−1,j − δi,nl Ti+1,j − δj,1 Ti,j−1 − δj,nc Ti,j+1 (24)

(UPMC) Thème 4 : Résolution numérique des EDP 2 et 3 mars 2017 27 / 60 (UPMC) Thème 4 : Résolution numérique des EDP 2 et 3 mars 2017 28 / 60
 Un problème statique : l’équation de Laplace Solution du problème Un problème statique : l’équation de Laplace Solution du problème

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

La méthode de relaxation est une méthode itérative qui utilise
Plusieurs méthodes permettent de résoudre le système linéaire (23), exclusivement l’équation (20) satisfaite par les températures aux points de
L2D Z = B, par exemple : la grille pour calculer Ti j , en fonctions de ses 4 voisins.
1 Les méthodes les plus générales de résolution numérique des systèmes Cela revient à séparer L2D en sa partie diagonale D = −411 et sa partie non
linéaires, telles que la méthode d’élimination de Gauss-Jordan, la diagonale N. On recherche la solution de l’équation NZ − B = −DZ = 4Z
décomposition « LU » 1 ; en partant d’une valeur quelconque Z0 , et calculant par récurrence :
Des méthodes adaptées au cas où A est une matrice « creuse », bande, Zn+1 = N Zn − B /4 dont le point fixe est la solution Z recherchée


2
symétrique...
On se convainc aisément que cette méthode est convergente en écrivant :
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 N n
Zn+1 − Z = (Z − Z)
l’équation de la chaleur. 4
4 L’utilisation des transformées de Fourier. et en montrant que la valeur absolue des valeurs propres de N est toujours
inférieure à 4 (rayon spectral inférieur à 4).
1. La factorisation « LU » exprime la matrice comme produit d’une matrice
Cette méthode (et ses variantes plus rapidement convergentes) est aisée à
triangulaire inférieure (Lower) par une matrice triangulaire supérieure (Upper). coder et s’adapte aisément à des géométries plus complexes.
(UPMC) Thème 4 : Résolution numérique des EDP 2 et 3 mars 2017 29 / 60 (UPMC) Thème 4 : Résolution numérique des EDP 2 et 3 mars 2017 30 / 60

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

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

Les méthodes générales, comme Gauss-Jordan ou LU demandent un On peut améliorer l’efficacité en tirant partie des propriétés de la matrice à
nombre d’opérations en N 3 pour résoudre un système de N équations 2 . inverser
Dans le cas présent, avec au plus 5 coefficients non nuls par ligne, on Par exemple on peut observer que −L2D est définie positive, et
devrait exploiter le caractère « creux » de la matrice, et descendre à un recourir à la procédure xPOSV.
nombre d’opérations proportionnel à N . On peut aussi s’appuyer sur fait qu’elle est symétrique, et utiliser la
En TE, nous utiliserons la bibliothèque Blas/Lapack optimisée, et on procédure correspondante xSYSV.
pourra donc explorer des régimes de tailles relativement intéressants, de Plus utilement, on remarquera que −L2D est une « matrice bande »,
l’ordre de 1000 à 10000 points. et une amélioration significative sera obtenue avec la fonction xSBSV
adaptée aux matrices symétriques-bandes.
On commencera avec la fonction Lapack xGESV pour les matrices
générales. Si enfin on tient compte de toutes ces propriétés, on utilisera avec
profit la procédure xPBSV.

2. Pour une matrice symétrique positive la factorisation de Choleski permet une

résolution en ∼ N 2
(UPMC) Thème 4 : Résolution numérique des EDP 2 et 3 mars 2017 31 / 60 (UPMC) Thème 4 : Résolution numérique des EDP 2 et 3 mars 2017 32 / 60
 Un problème statique : l’équation de Laplace Solution du problème Problèmes dynamiques

Stockage bande Introduction aux problèmes dynamiques

Les deux dernières fonctions (xSBSV et xPBSV) requièrent un stockage
spécifique de la matrice, appelé « stockage bande », où on ne fournit que Comme indiqué dans l’introduction générale, les problèmes dynamiques
les sur/sous-diagonales significatives (cf le « guide Lapack »). sont ceux où l’on suit l’évolution temporelle à partir d’un état « hors
Cette matrice réduite peut être construite directement en écrivant d’équilibre », ou plus généralement où l’on s’intéresse à des phénomènes de
successivement les (sous)-diagonales pertinentes. propagation.
En C, pour simplifier son obtention à partir de la matrice complète aplatie,
on pourra recourir à la procédure matfull2band fournie. Ils seront donc caractérisés par des « conditions initiales », et des
« conditions aux limites » susceptibles de dépendre du temps. Dans ce
ku=1 Matrice bande n=6 Stockage bande (kl+ku+1)×n dernier cas, on pourra étudier un système « forcé », qui, comme pour les
♣

systèmes linéaires décrits par une EDO, donne lieu à un transitoire suivi
♣
kl=2
d’un régime permanent forcé à la fréquence d’excitation.
♠

♦ ♥
Comme précédemment, nous allons nous concentrer sur le problème de
l’équation de la chaleur, et pour simplifier, réduire la dimensionnalité de
♦ ♠ l’espace à une dimension.
♥

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

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

Équation de diffusion 1-D Discrétisation du domaine [O, L] × [0, T ]

On cherche la solution numérique u(x, t) de l’équation de diffusion à 1-D, t
où le coefficient de diffusion D, nécessairement positif, est supposé L = (nx + 1)δx
uniforme et constant : n = nt + 1
∂u ∂2u

T = (nt + 1)δt
=D 2 (25) Un
∂t ∂x
point intérieur
dans le domaine spatial x ∈ [0, L], et temporel t ∈ [0, T ], avec des unj inconnu
δt nx points intérieurs

conditions initiales (CI) : Cond. aux Limites

en x = 0 et x = L n=1
u(x, 0) = u0 (x) ∀x (26)
Cond. Initiale
n=0 x
en t = 0 (n = 0)
et des conditions aux limites (CL) : j=0 j=1 δx j = nx + 1

u(0, t) = ug (t) et u(L, t) = ud (t) ∀t (27) Figure 6 – Discrétisation du domaine [0, L] × [0, T ]
NB : visualisation ultérieure avec les 2 bords
Il est à noter que l’espace et le temps jouent ici des rôles dissymétriques, ⇒ à chaque temps, écrire Un (nx points)
pas seulement à cause des CL/CI mais surtout à cause du degré différent plus les CL (2 points) : nx + 2 points
des dérivées partielles.
(UPMC) Thème 4 : Résolution numérique des EDP 2 et 3 mars 2017 35 / 60 (UPMC) Thème 4 : Résolution numérique des EDP 2 et 3 mars 2017 36 / 60
 Problèmes dynamiques Discrétisation de l’espace et du temps Problèmes dynamiques Discrétisation de l’espace et du temps

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

Il faut discrétiser les deux dimensions : Pour le Laplacien, nous utilisons la même expression que pour le problème
Espace : [0, L] est divisé en nx + 1 pas de largeur δx = L statique, qui est précise au deuxième ordre en δx :
nx +1 .
nx points intérieurs + 2 bords unj+1 − 2unj + unj−1
Temps : [0, T ] est subdivisé en nt + 1 pas de largeur δt = T
nt +1 .
∆u ≈
(δx)2
On note unj , la valeur numérique de u(x, t) pour x = j δx et t = n δt.
Les CL fixent donc les un0 ≡ ung et les unnx +1 ≡ und pour Pour la dérivée du premier ordre par rapport au temps, nous devons choisir
1 6 n 6 nt + 1. une expression, qui peut être :
Les CI fixent quant à elles les u0j pour 1 6 j 6 nx . vers l’avant, impliquant (un+1
j − unj )/δt (cf Euler forward)
La méthode de résolution sera fondamentalement itérative, chaque pas de vers l’arrière, avec (unj − un−1
j )/δt (cf Euler backward)
temps nécessitant de calculer les nx valeurs de un+1
j aux points intérieurs ou centrée avec (un+1
j − un−1
j )/2δt
à l’aide de uj et des conditions aux limites ug et un+1
n n+1
d . Le choix que nous allons faire détermine l’ordre de l’évaluation en δt, mais
aussi et surtout 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 37 / 60 (UPMC) Thème 4 : Résolution numérique des EDP 2 et 3 mars 2017 38 / 60

Problèmes dynamiques Algorithme explicite du premier ordre en temps Problèmes dynamiques Algorithme explicite du premier ordre en temps

Algorithme explicite du premier ordre en temps Système d’équations linéaires (algorithme explicite)
L’algorithme le plus simple et le plus naturel consiste à utiliser la dérivée
vers l’avant et donc à remplacer l’équation différentielle (25) par l’équation On peut alors écrire :
aux différences finies :
un+1
j = α unj+1 + (1 − 2α) unj + α unj−1 pour 1 6 j 6 nx (30)
un+1 unj j+1 − 2uj +
 n n
− u unj−1

=D pour 1 6 j 6 nx
j
(28)
δt (δx)2 qui met clairement en évidence le caractère explicite puisque les un+1
j
s’obtiennent directement en fonction des valeurs des unk à l’instant initial
Du point de vue temporel, cette approche peut être vue comme une
du pas de temps considéré.
extension de la méthode d’Euler explicite (progressive) étudiée sur les
Lorsque l’un des unj−1 ou α unj+1 n’est pas intérieur, c’est à dire que j = 1
EDO.
ou j = nx , on retrouve une condition de bord, et un0 ≡ ung ou unnx +1 ≡ und
Nous faisons ainsi apparaître le paramètre caractéristique sans dimension :
apparaîtront comme un terme additionnel, avec un poids α.
D δt
α= (29) Notons que la méthode de relaxation évoquée dans la première partie, si on
δx2
l’applique à une dimension, est exactement équivalente au cas où α = 21 .
qui joue un rôle essentiel dans toute la discussion qui va suivre.
(UPMC) Thème 4 : Résolution numérique des EDP 2 et 3 mars 2017 39 / 60 (UPMC) Thème 4 : Résolution numérique des EDP 2 et 3 mars 2017 40 / 60
 Problèmes dynamiques Algorithme explicite du premier ordre en temps Problèmes dynamiques Algorithme explicite du premier ordre en temps

Schéma de progression explicite Explicite : Formulation matricielle

En considérant le vecteur Un = (un1 , . . . , unnx ) des solutions à t = nδt pour
Il est utile de représenter cette situation par un schéma de progression :
t les nx points intérieurs, le système d’équations linéaires prend la forme :
Un+1 Un+1 = M(−α) Un + αVn ,
(n + 1)δt
où Vn = (ung , 0, . . . , 0, und ) est le vecteur des conditions aux limites,
Figure 7 – Schéma
nδt et la matrice est
illustrant la progression
d’un pas de temps avec Un M(−α) = 11 + α L1D
l’algorithme explicite du où L1D est la matrice tridiagonale représentant le laplacien à une
premier ordre en δt. dimension :

0
   
−2 1 0 ··· 0 1−2α α 0··· 0
x
.. .. 
0 j−1 j+1 nx + 1
. ..
j 1 .. .
  
 1 −2 .  α 1−2α α . 
. .. .. .. .. ..
   
L1D = 0 .. . .

0  −→ M(−α) =  0 . . . 0 
 (31)
un+1 = α unj+1 + (1 − 2α) unj + α unj−1
 
j
 .. . . ..
  .. .. ..

. . −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 41 / 60 (UPMC) Thème 4 : Résolution numérique des EDP 2 et 3 mars 2017 42 / 60

Problèmes dynamiques Algorithme explicite du premier ordre en temps Problèmes dynamiques Algorithme explicite du premier ordre en temps

Explicite : stabilité Explicite : analyse de stabilité de Fourier–von Neumann

Un problème important auquel on est confronté dans ce type de résolution Étudions la stabilité d’une onde stationnaire :
est l’instabilité de la méthode explicite
u(x, t) = A(t) sin(kx) ⇒ unj = A(n) sin (kjδx) , (32)
pour un pas de temps « un peu trop grand » (α > 0.5) :
0.3
1pas
c’est une solution exacte pour les CL ug = ud = 0, si k = pπ/L (p entier)
et si A(t) vérifie dA/ dt = −Dk 2 A(t) donc A(t) = A(0) exp −Dk 2 t


2pas
0.25 4pas
6pas
0.2 α=1.2 8pas Or l’algorithme explicite (28) donne, après division par unj :
10pas
0.15 12pas

sin(k(j + 1)δx) sin(k(j − 1)δx)

14pas  
0.1 A(n+1) −A(n)
= A(n+1)
−1 = α −2+
0.05
A(n) A(n) sin(kjδx) sin(kjδx)
= α 2 cos(kδx)−2 = −4α sin2
 
kδx
 
0 en développant les numérateurs 2 .
-0.05
soit une suite géométrique A(n) = A0 θn de raison réelle θ = 1 − 4α sin2
 
kδx
-0.1 2
-0.15
L’algorithme explicite est stable si |θ| 6 1 quel que soit le mode p, soit
0 20 40 60 80 100 120 140

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

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 43 / 60 (UPMC) Thème 4 : Résolution numérique des EDP 2 et 3 mars 2017 44 / 60
 Problèmes dynamiques Algorithme explicite du premier ordre en temps Problèmes dynamiques Algorithme explicite du premier ordre en temps

Stabilité de l’explicite : lien avec le spectre de L1D Exemple : Diffusion d’une distribution en « marche »
La p-ième valeur propre de la matrice L1D , en dimension n, s’écrit en Diffusion d'une marche
1.0
fonction de β = π/2(n + 1) = πδx/2L :
t =0
− 4 sin 2
pβ ; u(p) = sin m p β



valeur propre : vecteur propre : m
0.5 t→∞
On en déduit celles de la matrice M(−α) :
M(−α) = 11 + α L1D −→ 1 − 4α sin2 p πδx

2L .

u(x,t)
0.0
La condition de stabilité traduit donc simplement le fait que ces valeurs
propres sont de module inférieur à 1, c’est-à-dire que son rayon spectral
est inférieur à 1. 0.5

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

est remplacée par : 1.0
0 20 40 60 80 100 120
position : x
2i sin ωδt
exp(i ωδt
= 4α sin2 kδx




2 2 2 (33)
qui n’a de solution ω à partie imaginaire négative (stable) que si 4α 6 2. Figure 9 – Diffusion d’une marche avec des CL fixes +1 à gauche et −1 à
droite : u(x, t) tend vers une droite si t → ∞
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 45 / 60 (UPMC) Thème 4 : Résolution numérique des EDP 2 et 3 mars 2017 46 / 60

Problèmes dynamiques Algorithme explicite du premier ordre en temps Problèmes dynamiques Algorithme implicite du premier ordre en temps

Exemple : Diffusion d’une superposition de 2 modes Solution au problème de stabilité : méthode rétrograde
Rôle de la longueur d'onde sur la diffusion L’algorithme explicite est simple mais la condition de stabilité oblige à faire
2.4
de très petits pas...
2.2 On a vu, avec les EDO, que ces problèmes peuvent être résolus en utilisant
2.0 un algorithme implicite. C’est ce que nous faisons ici en choisissant la
dérivée temporelle « vers l’arrière », introduisant les nouvelles équations :
1.8

un+1 j+1 − 2uj

un+1 + un+1
n+1
− unj
" #
u(x,t)

1.6
=D
j j−1
(34)
1.4 δt (δx)2
1.2
Les un+1
j (j ∈ [1, nx ]) s’obtiennent alors en fonction des unj
1.0 (j ∈ [0, nx + 1]) en résolvant le système de nx équations linéaires :
0.8
0 20 40 60 80 100 120
j−1 + (1 + 2α) uj
−α un+1 j+1 = uj pour j ∈ [1, nx ]
n+1
position : x − α un+1 n
(35)

Figure 10 – Amortissement du mode 5 plus rapide que celui du mode 1 sauf pour j = 1 et j = nx où interviennent les CL.
mais pour t → ∞, on obtiendrait une droite. Soulignons le parallèle avec la méthode d’Euler implicite (rétrograde).
(UPMC) Thème 4 : Résolution numérique des EDP 2 et 3 mars 2017 47 / 60 (UPMC) Thème 4 : Résolution numérique des EDP 2 et 3 mars 2017 48 / 60
 Problèmes dynamiques Algorithme implicite du premier ordre en temps Problèmes dynamiques Algorithme implicite du premier ordre en temps

Schéma de progression implicite Formulation matricielle

Cette méthode met en œuvre le schéma de progression suivant : Ce système peut être écrit de façon matricielle sous la forme :
t M(α) Un+1 = Un + α Vn+1

ujn+1 avec les mêmes notations que précédemment, et notamment :

Figure 11 – Schéma (n+1) δt  
illustrant la 1+2α -α 0 · · · 0
progression d’un pas n δt ujn . . . . 
 -α . . . . . . .. 

de temps avec  . . .

M(α) = 11 − α L1D  0 .. .. .. 0 
=  (36)
l’algorithme implicite
du premier ordre en δt.
x  . . . .
 . .. .. ..


. -α
j δx (nx+1) δx 0 · · · 0 -α 1+2α
0
(j −1) δx (j +1) δx Notons que si M est tridiagonale, son inverse est au contraire une matrice
dense. Ainsi un+1
j dépend de tous les autres éléments, mais d’autant moins
qu’ils sont plus éloignés. Bien sûr, lorsque α → 0, M(α)−1 → M(−α) ce
j−1 + (1 + 2α) uj
−α un+1 j+1 = uj
n+1
− α un+1 n
qui signifie que cette dépendance tend à disparaître et on retrouve alors
l’algorithme explicite.
(UPMC) Thème 4 : Résolution numérique des EDP 2 et 3 mars 2017 49 / 60 (UPMC) Thème 4 : Résolution numérique des EDP 2 et 3 mars 2017 50 / 60

Problèmes dynamiques Algorithme implicite du premier ordre en temps Problèmes dynamiques Algorithme de Crank–Nicolson

Stabilité de l’algorithme implicite Principe de l’algorithme de Crank–Nicolson

On cherche à nouveau les modes propres définis par (32).

Si on reporte leur expression dans la nouvelle équation aux différences Pour accroître la précision, il faut donc améliorer l’approximation en
finies (34), on obtient la raison de la suite géométrique : différences finies du terme de dérivée partielle temporelle. Plutôt que de
modifier celui-ci directement, on utilise la même stratégie que pour la
1
θ=  61 (37) « méthode d’Euler modifiée » pour les EDO, qui est du 2nd ordre en temps.
1 + 4 α sin2

k δx
2
On fait donc la moyenne des seconds membres des deux algorithmes
qui montre que l’algorithme implicite est toujours stable quel que soit le précédents. On obtient l’algorithme de Crank–Nicolson :
choix des pas δx et δt.
un+1 D (uj+1 − 2uj + uj−1 ) + (unj+1 − 2unj + unj−1 )
n+1 n+1 n+1
− unj
" #
On obtient plus simplement le même résultat à partir des valeurs propres
=
j
,
de M(α)−1 = (11 − αL1D )−1 , inverses de 1 + 4α sin2 p 2Lπδx

. δt 2 (δx)2
(38)
Toutefois, cet algorithme étant basé sur une approximation du premier
qui est du deuxième ordre en δt comme en δx.
ordre en δt seulement, une description précise de la solution transitoire
requiert le choix d’un δt assez petit.

(UPMC) Thème 4 : Résolution numérique des EDP 2 et 3 mars 2017 51 / 60 (UPMC) Thème 4 : Résolution numérique des EDP 2 et 3 mars 2017 52 / 60
 Problèmes dynamiques Algorithme de Crank–Nicolson Problèmes dynamiques Algorithme de Crank–Nicolson

Schéma de progression de Crank–Nicolson Stabilité de l’algorithme de Crank–Nicolson

Le schéma de progression est désormais le suivant :
t L’étude de stabilité de von Neumann peut être reprise pour ce nouvel
algorithme et conduit à une décroissance géométrique des modes avec une
ujn+1 raison :
Figure 12 – Schéma (n+1) δt
1 − 2 α sin2
 
k δx
illustrant la progression n δt ujn
θ=
2
d’un pas de temps avec (39)
1 + 2 α sin2
 
k δx
l’algorithme de 2
Crank–Nicolson. x
ce qui permet de montrer que cet algorithme est lui aussi
j δx (nx+1) δx inconditionnellement stable (i.e. quel que soit le choix de δt).
0
(j −1) δx (j +1) δx
θCN (α) = θexplicite ( α2 ) × θimplicite ( α2 )
α α α α
− 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 53 / 60 (UPMC) Thème 4 : Résolution numérique des EDP 2 et 3 mars 2017 54 / 60

Problèmes dynamiques Algorithme de Crank–Nicolson Équation des ondes

Formulation matricielle Équation des ondes : différences finies

Les un+1
j s’obtiennent désormais en fonction des unj en résolvant le Évoquons brièvement la résolution de l’équation de d’Alembert :
système d’équations linéaires (pour j ∈ [1, nx ]) :
On discrétise le domaine comme pour l’équation de la chaleur ;
On se donne des conditions initiales arbitraires u(x, 0) et ∂u/∂t(x, 0) ;
α α α n α On se fixe pour simplifier des conditions aux limites de Dirichlet, ou
j−1 + (1 + α)uj
− un+1 j+1 = uj−1 + (1 − α)unj + unj+1
n+1
− un+1
2 2 2 2 de Neumann homogènes ;
(40) On essaye le schéma de différences finies « naturel » :
sauf pour j = 1 et j = nx où interviennent les CL.
∂2u 2
2∂ u
un+1 + un−1 − 2unj un + unj−1 − 2unj
= =c
j j 2 j+1
Ce système peut être écrit de façon matricielle sous la forme : c ⇒
∂t2 ∂x2 (δt)2 (δx)2
α −α αh n (42)
   
Un+1 = M Un + V + Vn+1
i
M (41)
2 2 2 qui est bien sûr du deuxième ordre en x et en t, et restaure une
certaine symétrie entre le temps et l’espace.
avec les mêmes notations que pour les deux autres algorithmes.

(UPMC) Thème 4 : Résolution numérique des EDP 2 et 3 mars 2017 55 / 60 (UPMC) Thème 4 : Résolution numérique des EDP 2 et 3 mars 2017 56 / 60
 Équation des ondes Équation des ondes

Algorithme explicite Critère CFL

δt
En posant α = c , on obtient les équations explicites à trois niveaux de
δx
temps :
Ce schéma « saute-mouton » donne une relation de dispersion (à 1D) :
un+1 =α 2
(unj+1 + unj−1 ) + 2(1 − α 2
)unj − un−1 . (43)
j j sin2 ω δt
= α2 sin2 k δx



2 2 (44)
soit, sous forme matricielle :
Si α 6 1, on a des solutions réelles (ie non divergentes). C’est le critère de
U n+1
= 211 + α L1D U − U
2
n n−1
+ α L1D V 2 n Courant–Friedrichs–Lewy prescrivant que α doit être « assez petit » pour
que la méthode soit stable. La valeur précise de la borne dépend de la
où V représente à nouveau les conditions aux limites sur les bords. méthode de discrétisation de l’équation, de la dimension de l’espace, et
éventuellement de la direction de propagation.
Ce schéma est appelé “leap-frog” en anglais, singulièrement traduit en De plus, la valeur α = 1 réalise la relation de dispersion exacte, et fait
« saute-mouton ». aussi disparaître unj de l’expression (43) de un+1
j .

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 57 / 60 (UPMC) Thème 4 : Résolution numérique des EDP 2 et 3 mars 2017 58 / 60

Équation des ondes Équation des ondes

Critère CFL Schéma de progression et critère CFL

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

t

L’idée générale est néanmoins toujours la même et peut être formulée des
3 façons équivalentes : Figure 13 – Progression
d’un pas de temps dans un+1
Le pas de temps δt doit être inférieur au temps mis par l’onde (à la l’algorithme explicite du (n + 1)δt
j

célérité c) pour traverser un pas d’espace δx. second ordre en temps et cδt
espace pour l’équation nδt
Le pas d’espace doit être supérieur à la distance parcourue par l’onde
durant δt. des ondes de d’Alembert
(n − 1)δt
un−1
La vitesse de « diffusion de l’information » δx/δt entre les noeuds de Critère CFL : cδt < δx 1/c j δt
la grille doit être supérieure à la célérité de l’onde. 0 L x
δ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 59 / 60 (UPMC) Thème 4 : Résolution numérique des EDP 2 et 3 mars 2017 60 / 60