Qu’est-ce qu’une image ?

Plusieurs réponses possibles

Traitement et analyse d’images numériques
Réponse "informatique" → un tableau de valeurs, une matrice de
pixels (picture element)
Partie 2: Dérivées, Opérateurs, Discrétisation
Image en niveau de gris, valeurs entre 0 (noir) et 255 (blanc)

Pierre Maurel
←→
Visages, IRISA/INRIA

[email protected]

Image encouleur, triplet de valeurs entre 0 et 255 (RGB : Red/Green/Blue)

(0, 0, 0) (128, 0, 0) (255, 0, 0)
 
 
←→   (0, 128, 0) (128, 128, 128) (255, 128, 0) 

 
https://perso.univ-rennes1.fr/pierre.maurel/IMA/ (0, 255, 0) (128, 255, 0) (255, 255, 255)

1 / 76 3 / 76

Qu’est-ce qu’une image ? Une image en tant que fonction : exemples

Plusieurs réponses possibles
f (x, y) = y f (x, y) = x + y f (x, y) = sin(y) f (x, y) = sin(x ∗ y)
A B C D
Réponse "mathématique" → une fonction , f : R × R → R
f associe un niveau de gris f (x, y) à une position (x, y)

1 2 3 4

A→4 B→3 C→1 D→2

En pratique : pour une image réelle, on

Autre visualisation de la fonction f (en fausses couleurs) n’a pas de formule explicite pour la −→ f =?
Ici, on a, par exemple (pour une image de taille 512 × 512) fonction f correspondante.
3 → f (320, 350) = 130 1 → f (1, 1) = 50 2 → f (400, 60) = 230

4 / 76 5 / 76
Manipulations basiques sur les images Autres possibilités

Images vectorielles
f (x, y) −→ composées d’entités mathématiques (cercles, lignes , etc)
chaque entité est décrite par une formule mathématique :

..............................................................................
g1 (x, y) = 255 − f (x, y) g2 (x, y) = f (y, x) g3 (x, y) = f (x, y)/3 g4 (x, y) = f (x, y) ∗ 2

A B C D

g1 → D g2 → A g3 → C g4 → B
Avantages : zoom infini, définir une image avec peu d’informations

manipulations de fonctions ... dérivées = ? Inconvénients : représentation de formes simples, pas d’images réalistes.

6 / 76 7 / 76

Autres possibilités Dérivées partielles : une image

Images 3D
On parle de voxels plutôt que de pixels.
En 1D : dérivée forte ⇒ forte variation de la fonction :

f (x) f 0 (x)

En 2D (ou plus) : la dérivée partielle dans une direction indique la variation de l’image dans cette
direction

∂f ∂f
| ∂x | | ∂y |
Images 2D+t
variante du 3D, il s’agit des vidéos

8 / 76 10 / 76
Dérivées partielles : une image Opérateurs usuels : Gradient

∂f
f ∂y
Le gradient en un point d’une image est un opérateur très important et
est utilisé dans de nombreux cas

f est un champ de scalaires, le gradient de f est un champ de vecteurs :

−−→
 
∂f ∂f
grad f = ∇f = , ...,
∂x1 ∂xn

f (10, 10) est un scalaire (∈ R) mais (∇f )(10, 10) est un vecteur (∈ R2 )

f (x, y) = x y2 + ex ⇒ ∇f (x, y) = y2 + ex , 2 x y


Exemple :
| {z } | {z }
∈R ∈R×R

11 / 76 12 / 76

Opérateurs usuels : Gradient Opérateurs usuels : Gradient

Exemple
f (x, y) = sin(y) −−→
le vecteur grad f (x, y) (ou ∇f (x, y)) indique la direction et l’intensité de la
∇f (x, y) = (0, cos(y)) plus grande variation du champ scalaire f autour du point (x, y)

interprétation : direction de "plus grande pente" et valeur de cette pente.

f ∇f

13 / 76 14 / 76
Opérateurs usuels : Gradient Opérateurs usuels : Gradient

les zones de fort gradient sont donc des zones de fortes variations les zones de fort gradient sont donc des zones de fortes variations

Norme duq gradient Norme duq gradient

−−→ ∂f 2 ∂f 2 −−→ ∂f 2 ∂f 2
f kgrad f k = ∂x + ∂y f kgrad f k = ∂x + ∂y

15 / 76 16 / 76

Opérateurs usuels : Gradient Opérateurs usuels : Gradient

Exemple de champ de gradients Exemple de champ de gradients

−−→ −−→
f grad f f grad f

17 / 76 18 / 76
Opérateurs usuels : Gradient Opérateurs usuels : Laplacien

Exemple de champ de gradients f est un champ de scalaires, le laplacien de f est également un champ
de scalaires :
∂2f ∂2f
∆ f = 2 + ... + 2
∂x1 ∂xn

Exemple : f (x, y) = x y2 + ex ⇒ ∆f (x, y) = ex + 2 x

| {z } | {z }
∈R ∈R

−−→ f ∆f
f grad f
Nombreuses applications physiques (eq de Poisson ∆φ = f )

19 / 76 20 / 76

Opérateurs usuels : Divergence Opérateurs usuels : Rotationnel

Divergence d’un champ w = (w1 , ..., wn ) (vecteurs −→ scalaires) :
 ∂    w = (w1 , ..., wn ) est un champ de vecteurs et son rotationnel est aussi un
∂x w1 champ de vecteurs :
∂w1 ∂wn 1
div w = + ... + =  ... · ...  = ∇ · w rot w = ∇ ∧ w
∂x1 ∂xn ∂
wn      
∂xn ∂/∂x u ∂w/∂y − ∂v/∂z
Cas 3D : rot w=∂/∂y∧v =∂u/∂z − ∂w/∂x
w(x, y) = y2 + ex , 2 x y ⇒ div w(x, y) = ex + 2 x


Exemple : ∂/∂z w ∂v/∂x − ∂u/∂y
| {z } | {z }
∈R
∈R2 Résumé :
interprétation physique : variation infinitésimale du volume autour d’un point grad : scalaire 7→ vecteur

div : vecteur 7→ scalaire

rot : vecteur 7→ vecteur
∆(laplacien) : scalaire 7→ scalaire

Remarque : Tous ces opérateurs sont linéaires.

−−→ −−→ −−→
grad(αf + βg) = α grad f + β grad g

21 / 76 22 / 76
Opérateurs usuels EDP : Équation aux Dérivées Partielles

Quelques formules Définition

−−→ Une EDP est une équation dont les inconnues sont des fonctions vérifiant
∆ = div(grad)
certaines conditions concernant leurs dérivées partielles.
−−→
rot(grad) = ∇ ∧ ∇ = 0
div(rot) = 0 =⇒ Omniprésence en sciences physiques :
−−→
rot(rot) = ∇ ∧ ∇w = grad(div) − ∆ Mécanique des fluides (équations de Navier-Stokes)
−−→ ~ ~ −−→ −−→
grad(A · B) = (~A · grad)~B + ~A ∧ rot(~B) + (~B · grad)~A + ~B ∧ rot(~A) Électromagnétisme (équations de Maxwell)
−−→ −−→ −−→
grad(fg) = f grad(g) + ggrad(f ) Réactions chimiques
−−→ Prévisions météorologiques
div(f ~A) = f div(~A) + grad(f ) · ~A
−−→ ...
rot(f ~A) = f rot(~A) + grad(f ) ∧ ~A
−−→ −−→
∆(fg) = f ∆g + 2grad(f ) · grad(g) + g∆f =⇒ Outil très puissant de représentation.

23 / 76 24 / 76

Équation de la chaleur
Équation de la chaleur sur une image
décrit le phénomène physique de conduction thermique (Fourier, 1811)
évolution de la température sans contraintes extérieures
T(x, y, t) le champ de température sur un domaine Ω :
∂T(x, y, t)
∀(x, y) ∈ Ω, = ∆T(x, y, t) et T(x, y, 0) = T0
∂t
2 2
(attention le laplacien est uniquement en "espace", sur x et y : ∆ T = ∂ T + ∂ T .)
∂x2 ∂y2

On parle de diffusion isotrope : pas d’orientation préférentielle.

25 / 76 26 / 76
Résolution de l’équation de la chaleur Résolution de l’équation de la chaleur

Cas 1D : diffusion isotrope du signal u0 sur R : On a donc :

∂u(x, t) ∂ 2 u(x, t) ∂u(x, t) ∂ 2 u(x, t)
∂u(x, t) ∂ 2 u(x, t) = ⇒ TF( ) = TF( )
= avec u(x, 0) = u0 (x), ∀x ∂t ∂x2 ∂t ∂x2
∂t ∂x2 ∂TF(u)
⇒ = −ω 2 TF(u)
∂t
2
On utilise la transformée de Fourier : (attention les constantes peuvent différer selon les définitions) ⇒ TF(u)(ω, t) = f (ω) e−ω t
(f est quelconque).
Z +∞ Z +∞
1
TF(f )(ω) = f (x)e−i ω x dx et TF−1 (f̂ )(x) = f̂ (ω)ei ω x dω Condition initiale : u(x, 0) = u0 (x) , donc TF(u)(ω, 0) = TF(u0 )(ω)
−∞ 2π −∞

on prend donc f (ω) = TF(u0 )(ω)

Ce qui permet de se "débarrasser" des dérivées spatiales
Z +∞ Z +∞ La transformée de Fourier transforme la multiplication en convolution
∂u(x, t) ∂u(x, t) −i ω x ∂ ∂TF(u)
TF( )(ω) = e dx = u(x, t)e−i ω x dx = 2
 2

∂t −∞ ∂t ∂t −∞ ∂t TF(u) = TF(u0 ) e−ω t ⇒ u = TF−1 TF(u0 ) e−ω t
 2

∂ 2 u(x, t) ⇒ u = u0 ∗ TF−1 e−ω t
TF( )(ω) = − ω 2 TF(u)(ω) (par intégration par parties et en supposant u à support borné)
∂x2
La transformée de Fourier d’une gaussienne est une gaussienne.

27 / 76 28 / 76

Résolution de l’équation de la chaleur Équation de la chaleur sur une image

Retour au 2D+t
On montre que la solution est : u(x, t) = g(x, t) ∗ u0 (x)
1 x2 √ ∂u(x, y, t)
avec g(x, t) = √ e− 4t → Gaussienne d’écart type σ = 2t = ∆u(x, y, t) et u(x, y, 0) = u0 (x, y)
2 πt ∂t

g est la fonction de Green associée à l’équation de la chaleur.

on peut montrer que la solution est donnée par

u(x, y, t) = G (x, y, σ(t)) ∗ u0 (x, y)

Solution de l’équation de la chaleur 1D
Si la donnée initiale u0 est suffisamment régulière, la solution explicite de où G (x, y, σ(t)) est une gaussienne dont l’écart-type σ est
l’équation proportionnel à t et ∗ est l’opérateur de convolution.
∂u(x, t) ∂ 2 u(x, t)
= avec u(x, 0) = u0 (x), ∀x
∂t ∂x2
est donnée par :
u(x, t) = G√2 t (x, t) ∗ u0 (x)

29 / 76 30 / 76
Équation de la chaleur sur une image Interprétation de l’équation de la chaleur par la
diffusion
(0 < t1 < t2 < t3 < t4 < t5 )

Diffusion isotrope de l’image I0

∂I
= ∆I = div(∇I)
∂t

On diffuse selon toutes les directions de ∇I

p p
En particulier on diffuse autant au niveau des contours qu’à l’intérieur
u(x, y, 0) = u0 (x, y) u(x, y, t1 ) = G(x, y, 2 t1 ) ∗ u0 (x, y) u(x, y, t2 ) = G(x, y, 2 t2 ) ∗ u0 (x, y)
−→ trop flou

Diffusion non-linéaire (« anisotrope »)

∂I  
= div f (x, y)∇I
∂t

là où f est petit : divergence faible, peu de variations temporelles

p p p
u(x, y, t3 ) = G(x, y, 2 t3 ) ∗ u0 (x, y) u(x, y, t4 ) = G(x, y, 2 t4 ) ∗ u0 (x, y) u(x, y, t5 ) = G(x, y, 2 t5 ) ∗ u0 (x, y)

là où f est grand : divergence grande, diffusion importante

31 / 76 32 / 76

Exemple de diffusion non linéaire : Perona-Malik Exemple de diffusion non linéaire : Perona-Malik
Principe : on veut préserver les contours et lisser les régions homogènes de Lissage des zones à faible gradient (réduction du bruit) → c(0) = 1
l’image
Atténuation de la diffusion lorsque le gradient est important (préservation
Perona-Malik, Scale-space and edge detection using anisotropic diffusion IEEE des singularités et contours) → lim c(x) = 0
Transactions on Pattern Analysis and Machine Intelligence, 1990 x→∞

article cité plus de 10000 fois ! (Google Scholar)

 2
1 u
la fonction f va dépendre de la norme du gradient : f (x, y) = c(k∇I(x, y)k)
exemples : c(u) = 2 c(u) = exp − 2
1 + λu 2 λ

Cette diffusion est en fait toujours isotrope. En effet c ne dépend pas de la direction du gradient mais
uniquement de sa norme.

33 / 76 34 / 76
Exemple de diffusion non linéaire : Perona-Malik Exemple de diffusion non linéaire : Perona-Malik

diffusion isotrope diffusion Perona-Malik diffusion isotrope diffusion Perona-Malik

35 / 76 36 / 76

Exemple de diffusion non linéaire : Perona-Malik Du continu au discret

théorie : fonctions "continues" (définies en tout point)

en fait : peu d’applications vraiment convaincantes en images MAIS
cas 1-D + temps : v : (x, t) ∈ R × R → v(x, t) ∈ R
beaucoup d’utilisations pour texturer des modèles (2d ou 3D)
cas 2-D + temps : v : (x, y, t) ∈ R2 × R → v(x, y, t) ∈ R

mais en pratique : fonctions discrétisées, définies sur une grille finie k ∆x

et pour certains instants n ∆t

continu discret

37 / 76 39 / 76
Du continu au discret Continu / discret

Problème : l’EDP que l’on cherche à résoudre n’est définie qu’en

continu, cf définition des dérivées partielles :

∂f f (a1 , ..., ai−1 , ai + h , ai+1 , ..., an ) − f (a1 , ..., an )

(a) = lim
∂xi h→0 h

−→ comment approximer ces dérivées ?

plusieurs possibilités : éléments finis, méthodes spectrales, différences

finies

On va donc chercher une solution à une nouvelle EDP discrète →

solution approchée de l’EDP originale (continue)

40 / 76 41 / 76

Continu / discret "Construction" de la solution discrète

“Avantage” : en discret, on n’a pas forcément besoin de trouver la

solution exacte (analytique) de l’EDP.

Pour l’équation de la chaleur (en continu), on a vu qu’on pouvait trouver

la solution analytique :

∂2u

∂u ? ?
= 2

∂t ∂x =⇒ u(x, t) = Gσ(t) ∗ u0 (x) u(x, 0) = u0 (x) u(x, ∆t) = F(u(x, 0)) u(x, 2 ∆t) = F(u(x, ∆t))
u(x, 0) = u0 (x)


Mais c’est très souvent impossible

En discret on peut espérer construire une solution à partir de la

condition initiale (fonction ou image initiale)

En particulier, s’il y a une variable temporelle (f (x1 , . . . , xn , t)), on peut ? ? ?

partir de la condition initiale et “construire” la solution à l’instant (n + 1) ∆t u(x, 3 ∆t) = F(u(x, 2 ∆t)) u(x, 4 ∆t) = F(u(x, 3 ∆t)) u(x, 5 ∆t) = F(u(x, 4 ∆t))
en fonction de la solution à l’instant n ∆t

42 / 76 43 / 76
Discrétisation des dérivées Conditions aux bords
En continu, f : R × R → R : on défini
∂f f (x + h, y) − f (x, y) I est une image définie sur Ω
= lim
∂x h→0 h Ω est borné
comment définir les dérivées de I sur les bords de Ω (noté Fr(Ω)) ?
En discret, la fonction/image est définie sur une grille (pixels) : Même problème pour la convolution (filtrage)

u(i ∆x, j ∆y) → on prolonge de manière arbitraire l’image

ne pas traiter le bord : nouvelle image plus petite que l’originale
Comment définir la dérivée partielle de u par rapport à x ? compléter l’extérieur de l’image originale avec des 0 (padding)
compléter l’extérieur de l’image originale en symétrisant ses valeurs aux
bords (mirroring)
u(i ∆x + ∆x, j ∆y) − u(i ∆x, j ∆y) u(i ∆x − ∆x, j ∆y) − u(i ∆x, j ∆y)
? ... ? compléter l’extérieur de l’image originale par extrapolation de ses valeurs
∆x −∆x aux bords
u(i ∆x + ∆x, j ∆y) − u(i ∆x − ∆x, j ∆y) compléter l’extérieur de l’image originale par périodisation
... ?
2 ∆x

44 / 76 45 / 76

Conditions aux bords Différences finies

∂I
Notation “petit o” : “négligeable devant”
Exple : par réflexion ( ∂N = 0 sur Fr(Ω))
f (x)
f (x) = x→a
o (g(x)) , lorsque −→ 0
g(x) x→a
Exemples : xn+1 = o (xn ) , 28 xn+1 + xn+3 = o (xn )
x→0 x→0

Formule de Taylor-Young :
f 00 (x) 2 f (n) (x) n
f (x + h) = f (x) + f 0 (x) h + h + ... + h + o (hn )
2 n! h→0

Formule de Taylor-Young pour une fonction à plusieurs variables :

∂u 1 ∂nu
u(x + h, y) = u(x, y) + (x, y) h + . . . + (x, y) hn + o (hn )
∂x n! ∂xn h→0

46 / 76 47 / 76
Différences finies Différences finies
Soit v : R × R → R une image Filtres associés
∂v
 
on a v(x + ∆x, y) = v(x, y) + ∆ x (x, y) + o (∆x) Gx ∗ I
et Gy = GTx
   
∂x ∇I = avec Gx = −1 1 ou −0.5 0 0.5
∆x→0 Gy ∗ I
v(x + ∆x, y) − v(x, y) ∂v 
0 1 0

donc = (x, y) + o(1)
∆x ∂x ∆I = D ∗ I, D = 1 −4 1
0 1 0
Pour une image, ∆x et ∆y n’ont pas de sens particulier et sont
généralement pris égaux à 1 pour simplifier les formules
∂u
Approximations classiques pour (x, y) :
∂x
Schéma "avant" : v(x + 1, y) − v(x, y)

Schéma "arrière" : v(x, y) − v(x − 1, y)

v(x + 1, y) − v(x − 1, y) ∂I ∂I
Schéma "centré" : ∂x → ∂y ↓ I
2
48 / 76 49 / 76

Différences finies Sensibilité au bruit

∂I ∂I
∂x → ∂y ↓ I

∂2I ∂2I
∆I = ∂x2 + ∂y2 I

50 / 76 51 / 76
Sensibilité au bruit Différentiation et lissage
pour augmenter la robustesse au bruit :

Lisser l’image avant différentiation (par filtrage linéaire)

∂I
:= Gx ∗ (H ∗ I)
∂x

Si H filtre moyenneur 3 × 3 → Filtre de Prewitt

∂I ∂Gσ ∗ I ∂Gσ
Si H filtre gaussien → := = ∗I
∂x ∂x ∂x

∂Gσ
←→
∂x

σ = 1.5

52 / 76 53 / 76

Différentiation et lissage Différentiation et lissage

Sans lissage pour comparaison

σ = 1.5

54 / 76 55 / 76
Différentiation et lissage Différentiation et lissage
Laplacien de Gaussienne (LoG) Laplacien de Gaussienne (LoG)
∂2I ∂ 2 Gσ ∗ I ∂ 2 Gσ ∂2I ∂ 2 Gσ ∗ I ∂ 2 Gσ
:= = ∗I := = ∗I
∂x2 ∂x2 ∂x2 ∂y2 ∂y2 ∂y2 σ = 1.5

∂2I ∂2I
∆I = 2
+ 2 := (∆Gσ ) ∗ I
∂x ∂y

∆Gσ ←→

56 / 76 57 / 76

2
∂v ∂ v
Différentiation et lissage Exemple : équation de la chaleur ∂t = α ∂x 2
Retour aux EDP
Sans lissage pour comparaison
on utilise le développement de Taylor de v(x, t) : R × R → R et on note
vnk = v(k ∆x, n ∆t) :

∂v
(premier ordre) v(k∆x, (n + 1)∆t) = vn+1
k = vnk + ∆t + o(∆t)
∂t
∂v ∆x2 ∂ 2 v
(second ordre) v((k + 1)∆x, n∆t) = vnk+1 = vnk + ∆x + + o(∆x2 )
∂x 2 ∂x2
∂v ∆x2 ∂ 2 v
(second ordre) v((k − 1)∆x, n∆t) = vnk−1 = vnk − ∆x + + o(∆x2 )
∂x 2 ∂x2
on en déduit des approximations des dérivées partielles :

∂v vn+1 − vnk
(k∆x, n∆t) = k + o(1)
∂t ∆t
∂2v vn − 2vnk + vnk−1
2
(k∆x, n∆t) = k+1 + o(1)
∂x ∆x2

58 / 76 59 / 76
 2
∂v ∂ v
Exemple : équation de la chaleur ∂t = α ∂x 2 Analyse numérique
Retour aux EDP

Comment choisir ∆x, ∆t, . . . ?

∂v ∂2v vkn+1 − vnk vn − 2vnk + vnk−1
Donc −α 2 =0 ⇒ − k+1 + o(1) = 0
∂t ∂x ∆t ∆x2
DEMO MATLAB

et finalement une approximation "raisonnable" de l’EDP semble être :

Précautions lors de la discrétisation numérique
∆t
un+1
k = (1 − 2r)unk + r(unk+1 + unk−1 ) avec r=α 2 1 Notion de consistance
∆x
2 Notion de stabilité

60 / 76 61 / 76

Précautions lors de la discrétisation numérique Précautions lors de la discrétisation numérique

EDP, Notation continue : v : R × R+ → R
(
L(v)(x, t) = g(x, t) t > 0, x ∈ R Définition : schéma convergent
(L est l’opérateur définissant l’EDP )
v(x, 0) = f (x) x∈R On dit qu’il y a convergence si, quand ∆x et ∆t → 0 :
(1)
(k ∆x, n ∆t) converge vers (x, t) =⇒ unk converge vers v(t, x)
∂v ∂2v
Exemple : pour l’équation de la chaleur L(v) = − α 2 et g = 0
∂t ∂x
→ "si les pas de discrétisation deviennent petits, alors la solution discrète
Notation discrète : unk ∈ R, pour n ∈ N (temps) et k ∈ Z (espace) tend bien vers la solution continue théorique"
(
L∆ (unk ) = g(k ∆x, n ∆t) En général, la convergence est dure à montrer mais :
(2)
u0k = f (k ∆x)
Théorème de Lax
Un schéma consistant pour un problème linéaire bien posé converge si et
Problématique : choisir L∆ (à partir de L) pour que la solution discrète seulement si il est stable.
unk de (2) soit une "bonne approximation" de la solution continue v(x, t)
de (1).
62 / 76 63 / 76
Consistance Consistance : équation de la chaleur

Consistance – Erreur de troncature Équation de la chaleur : on a montré que

On appelle "erreur de troncature" (ou "erreur de discrétisation") l’expression
∂v ∂2v vn+1 − vnk vn − 2vnk + vnk−1
−α 2 = k − α k+1 + o(1)
R∆ (φ) = L∆ (φnk ) − L(φ) ∂t ∂x ∆t ∆x2

Le schéma est dit consistant si, pour toute fonction φ, l’erreur de troncature avec vnk = v(k ∆x, n ∆t).
tend vers 0 lorsque ∆ tend vers 0 (i.e. lorsque tous les pas de discrétisation
tendent vers 0) !
∂2v vn+1 − vnk vn − 2vnk + vnk−1
 
∂v
⇒ −α 2 − k
− α k+1 −−−−−→ 0
→ "Lorsque tous les pas de discrétisation tendent vers 0, les ∂t ∂x ∆t ∆x2 ∆x,∆t→0
approximations discrètes des dérivées présentes dans L∆ tendent bien
vers les dérivées continues de L"

"Les dérivées sont bien approximées" Donc, le schéma de discrétisation obtenu précédemment est bien
∂v ∂2v
→ se montre facilement en utilisant la formule de Taylor. consistant avec l’équation de la chaleur = α 2.
∂t ∂x

64 / 76 65 / 76

Consistance : exercices
y(x + h) − y(x − h)
est-il une bonne approximation de y0 (x) ?
h
NON : 2y0 (x)

y(x + 2h) − 2y(x + h) + y(x)

approxime ?
h2 Stabilité d’un schéma numérique
00
y (x)

Déterminez une approximation consistante d’ordre 2 de y0 (x) faisant

intervenir une combinaison linéaire de y(x) , y(x − h) et y(x − 2h).

y(x − 2h) − 4y(x − h) + 3y(x)

y0 (x) = + o(h)
2h

66 / 76 67 / 76
Stabilité d’un schéma numérique Stabilité d’un schéma numérique

Définition
Stabilité du schéma un+1
k = (1 − 2r)unk + r(unk+1 + unk−1 ) en considérant la
On dit qu’un processus de calcul séquentiel (ou itératif) est « stable » si les n n
norme ku kL ∞ = sup |uk |.
erreurs d’arrondis ne s’amplifient pas lors de la progression des calculs. k
si (1 − 2r) ≥ 0, un+1 ≤ (1 − 2r)kun kL ∞ + rkun kL ∞ + rkun kL ∞
kun+1 k ≤ Kku0 k k
kun+1 kL ∞ = sup un+1
k ≤ kun kL ∞ , ∀n ∈ N
k
normes usuelles :
X sX kun+1 kL ∞ ≤ kun kL ∞ ≤ ku0 kL ∞
kvkL 1 = |vk | kvkL 2 = |vk |2 kvkL ∞ = sup |vk |
k
k k
Si r ≤ 12 , le schéma est stable.
Théorème de Lax → "si on a bien choisi nos discrétisations, soit un reste
borné et on a convergence, soit un explose"

68 / 76 69 / 76

Stabilité d’un schéma numérique Stabilité d’un schéma numérique : méthode de Fourier

En général : plus compliqué. Une méthode possible, passer en Fourier

Identité de Parseval : pour prouver la stabilité on peut raisonner en Fourier.
En continu En discret
Z +∞ +∞ on veut montrer que kun+1 kL 2 ≤ Kku0 kL 2
Transformée de 1 1 X
v(t, ω) = √ e−iωx v(t, x)dx ubk (ξ) = √ e−ikξ uk
Fourier
b
2π −∞ 2π k=−∞
Z +∞ Z +π or : kukL 2 = kb
ukL 2 (← Parseval)
Transformée de 1 1
v(t, x) = √ eiωxb
v(t, ω)dω uk = √ eikξ b
u(ξ)dξ
Fourier inverse 2π −∞ 2π −π
Remarques : •i désigne le complexe tel que i2 = −1 • les constantes (ici = √1 ) peuvent varier suivant les définitions.
u0 kL 2
un+1 kL 2 ≤ Kkb
on peut donc montrer que kb
2π

Identité de Parseval kvkL 2 (R) = kb

vkL 2 (R) kukL 2 = kb
ukL 2 (−π,π)
Remarque :
+∞
X Z +π
2
kukL 2 = |uk | et kûkL 2 = |û(ξ)|2 dξ
k=−∞ −π
F IGURE: Rappel sur les transformées de Fourier (continue et discrète)

70 / 76 71 / 76
Stabilité d’un schéma numérique Stabilité d’un schéma numérique : méthode de Fourier
Calculs préliminaires :
+∞ +∞ Pour étudier la stabilité d’un schéma numérique de la forme (avec f fonction linéaire)
1 X −ikξ 1 X −ikξ n
Si f : k → unk , on a donc bf (ξ) = √ e f (k) = √ e uk
2π k=−∞ 2π k=−∞ un+1 = f ( . . . , unk−1 , unk , unk+1 , . . . )
k
si on pose g : k → unk−1 , on a :
+∞ +∞
1 X −ikξ n 1 X
g(ξ) = √ e uk−1 = √ e−i(m+1)ξ unm = e−i ξ bf (ξ) n+1 n bn d n
Prendre la transf. de Fourier → ud = f ( . . . , ud
k−1 , uk , uk+1 , . . . )
b 1
2π k=−∞ 2π m=−∞ k

n n i p ξ bn n+1
Propriété : Transf. de Fourier : on se “débarrasse” des dérivées spatiales. 2 Se "débarrasser" des ud
k+p : k+p = e
ud uk ⇒ ud
k (ξ) = ρ(ξ)ubnk (ξ)

∂2v
d
en continu (résolution équation de la chaleur) 2 (t, ω) = −ω 2bv(t, ω) 3 Par récurrence :
∂x
n+1 n+1
∂2u ud
k (ξ) = ρ(ξ)n+1 ub0 (ξ)
k ⇒ |ud
k (ξ)| = |ρ(ξ)|n+1 |ub0k (ξ)|
en discret : 2 (k ∆x, n ∆t) ≈ unk+1 − 2unk + unk−1 unk+1 − 2unk + unk−1 (si ∆x = 1,
∂x | {z }
A
comme svt en image)
4 Il faut donc avoir |ρ(ξ)| ≤ 1, ∀ξ ⇒ Condition sur ∆t , ∆x, . . .

5 Conclusion avec l’identité de Parseval

A(ξ)
b = ei ξ ubnk (ξ) − 2 ubnk (ξ) + e−i ξ ubnk (ξ) = 2(cos(ξ) − 1) ubn (ξ)
72 / 76 73 / 76

Exemple : équation de la chaleur Exemple : équation de transport

EDP de transport/propagation : on cherche une fonction
∂v ∂2v
EDP : =α 2
∂t ∂x

 ∂v + c ∂v = 0 x ∈ R, t ∈ R+
v : R × R+ → R
On a vu une discrétisation consistante possible : vérifiant ∂t ∂x
(x, t) 7→ v(x, t)
v(x, 0) = v0 (x) x∈R

∆t
un+1
k = (1 − 2r)unk + r(unk+1 + unk−1 ) avec r=α 2 où c est une constante réelle.
∆x

 n+1 Solution : v(x, t) = v0 (x − ct)

n+1 (ξ) = 1 + 2r(cos(ξ) − 1) ubn (ξ) = 1 + 2r(cos(ξ) − 1)
ud ub0 (ξ)
et on peut montrer que
Stable si, et seulement si, |1 + 2r(cos(ξ) − 1)| ≤ 1, pour tout ξ
⇐⇒ −1 ≤ 1 + 2r(cos(ξ) − 1) ≤ 1, pour tout ξ Schéma centré → toujours instable
Si c > 0, l’information est propagée de gauche à droite → schéma arrière
⇐⇒ −1 ≤ r(cos(ξ) − 1) ≤ 0, pour tout ξ
Si c < 0, l’information est propagée de droite à gauche → schéma avant.
1
⇐⇒ 0 ≤ r ≤ , pour tout ξ ∆t
1 − cos(ξ) la condition de stabilité de ces 2 schémas peut s’écrire : : |c| ≤ 1.
∆x
1
⇐⇒ 0 ≤ r ≤
2 DEMO MATLAB
DEMO MATLAB
74 / 76 75 / 76
Exercice : équation de la chaleur (schéma implicite)
un schéma est dit explicite si on peut écrire un+1
k en fonction de uni pour
certains i
exemple : un+1
k = (1 − 2r)unk + r(unk+1 + unk−1 )
facile à implémenter

un schéma est dit implicite si ce n’est pas le cas

Le schéma suivant est un schéma implicite consistant avec l’équation de la
chaleur
un+1 − unk un+1 − 2un+1 n+1
+ uk+1
k
= k−1 k
2
∆t ∆x
Étudier sa stabilité
+∞
1 X
Rappel : ubk (ξ) = √ e−imξ um , (où i désigne le complexe i2 = −1)
2π m=−∞
On obtient b un+1 (ξ) = 1+2r(1−cos(ξ))
1
un (ξ) or on a 1 − cos(ξ) ≥ 0 donc
b
1
1+2r(1−cos(ξ))
≤ 1 et on obtient un schéma inconditionnellement stable (mais
difficile à implémenter car schéma implicite)
76 / 76