1/18

Chapitre 1 : Introduction à L’Analyse Numérique

...
 2/18

1 Définition

L’analyse numérique est la conception et l’étude d’algorithmes pour

obtenir des solutions à des ensembles d’équations issus de modèles
issus de la physique, de la biologie, de la finance ...

2 Motivations :
Recherche et développement : études expérimentales coûteuses
Les modèles considérés sont composés d’ensemble d’équations dont
on ne sait pas déterminer de solutions explicites
Proposer une solution approchée, calculée à l’aide de l’ordinateur.
3 Développer des algorithmes efficaces
Convergence et stabilité de la méthode numérique
Coût algorithmique
 3/18

Plan du cours
1 Introduction à l’analyse numérique
2 Interpolation et approximation
3 Intégration numérique
4 Résolution de systèmes linéaires
5 Equations non linéaires
6 Equations différentielles ordinaires
7 Equations aux dérivées partielles
 4/18

Quelques exemples
Mouvement du pendule
Equation

θ l  g
θ00 (t ) + l sin(θ(t )) = 0,


θ(0) = θ0 , θ0 (0) = θ1



mg
Solution exacte ?

Solution approchée pour θ << 1

g

θ (t ) + l θ(t ) =q0,  q

 00
 q 
g g
θ(t ) = θ0 cos l
θ



l t + g 1 sin lt

Utiliser une approximation numérique de la solution !

 5/18

Quelques exemples

Calculer les racines du polynôme p (x ) = ax 2 + bx + c

 6/18

Quelques exemples
Temps de calcul pour l’inversion d’une matrice
 7/18

1 Représentation des réels sur l’ordinateur

2 Conditionnement, stabilité et complexité

 8/18

Représentation des réels en base b

Tous nombre réel x peut être représenté sous la forme

x = ±m b ±e , avec b ≥ 2.

avec la mantisse m :

m = m1 b −1 + m2 b −2 + ..., et mi ∈ {0, 1, ..., b − 1}

et l’exposant e :

e = e0 b 0 + e1 b 1 + ...es −1 b s −1 , avec s∈N

Représentation unique ?
Ordinateur : mémoire finie !
La mantisse est tronquée au bout de r chiffres,
Valeur maximale pour s
 9/18

Le standard IEE (754-1985)

Float double précision : utilisation de 8 octets ( 64 bits ), b = 2

Signe: 1 bit e : 11 bits m : 52 bits

Pour la mantisse : utilisation de 52 bits

m = 2−1 + m2 2−2 + m3 2−3 + ... + m53 2−53 ,

Pour l’exposant : 11 bits : e ∈ [−1022, 1025] avec

e = c0 20 + c1 21 + ... + c10 210 − 1022, avec ci ∈ {0, 1}

Les valeurs e = −1022 et e = 1025 sont réservées à la

représentation de 0 et de Inf
On appelle F l’ensemble fini de ces nombres
 10/18

Exercices

Exercice
Calculer xmax le nombre le plus grand de F et xposmin , le nombre positif le
plus petit de F

Exercice
Proposer deux algorithmes pour déterminer respectivement une
approximation numérique de xmax et xposmin
 11/18

Erreur d’arrondi

Seuls les éléments de F sont autorisés

Un nombre réel x ± m 2e , m = 0.1m2 m3 ... pour xposmin < |x | < xmax
est arrondi de la manière suivante

0.1m2 ...m53 2e

 si m54 = 0,
rd (x ) = sign(x ) 
(0.1m2 ...m53 + 2−53 )2e
 si m54 = 1,

Exercice
Calculer le nombre le plus petit tel que rd (x + 1) > 1 et proposer un
algorithme pour retrouver numériquement ce nombre
 12/18

Erreur arrondi

Si xposmin < |x | < xmaxx , alors

L’erreur d’arrondi absolue est
1 −53 e
|x − rd (x )| ≤ 2 2 .
2
L’erreur d’arrondi relative est
 
 x − rd (x )  1 2−53 2e
 ≤ ≤ 2−53 ' 10−16
x 2 |m|2e

 13/18

1 Représentation des réels sur l’ordinateur

2 Conditionnement, stabilité et complexité

 14/18

Conditionnement d’un problème numérique

Definition
Le conditionnement représente la sensibilité du résultat par rapport à de
petites variations des données ...
On dit qu’un problème est
- bien conditionné si une petite variation des données entraîne une petite
variation du résultat
- mal conditionné si une petite variation des données entraîne une grande
variation du résultat
Exemple : soit f : R → R
Développement de Taylor
f (x + δx ) = f (x ) + f 0 (x ) ∗ δx + o (δx )
δf := f (x + δx ) − f (x ) ' f 0 (x )δx
Conditionnement
δf x xf 0 (x )
k = '
f δx f (x )
 15/18

Conditionnement des opérations arithmétiques

élémentaires

Soit f ∈ C2 (Rn , R). Quel est le conditionnement de f ?

D’après le développement de Taylor de f ,
n
X ∂f
δf = f (x + dx ) − f (x ) = δx + 0(|δx |2 ).
i =1
∂xi i

Le conditionnement est donc déterminé par les nombres

∂f xi
ki = .
∂xi f (x )

Exercice
Calculer le conditionnement de l’addition et de la multiplication
 16/18

Stabilité d’un algorithme

Definition
La stabilité d’un algorithme se réfère à la propagation des erreurs au cours
des étapes du calcul, à la capacité de l’algorithme de ne pas trop amplifier
d’éventuels écarts, à la précision des résultats obtenus.

Exemple : l’algorithme qui calcule les racines du polynôme

p (x ) = ax 2 + bx + c basé sur les formules
√ √
−b + b 2 − 4ac −b − b 2 − 4ac
x1 = , x2 =
2a 2a

s’avère instable en pratique : exemple avec a = 10−20 , b = 1, c = 1

Exercice
Proposer un autre algorithme plus stable
 17/18

Complexité algorithmique

Definition (Coût algorithmique)

Nombre d’opérations élémentaires effectuées par l’algorithme (+, *, sqrt,
puissance ...)
Pn
Exemple : Evalutation du polynôme p (x ) = i =0 ai x i :
Pn n ( n + 1)
Nombre de multiplications : i =0 i= 2
Nombre d’additions : n
Coût algorithmique → O (n2 )

Definition




 O (nk ) → coût polynomial

O (a n ) → coût exponentiel





O (n!) → coût factoriel

 18/18

Complexité algorithmique

Pn
Exemple : Evalutation du polynôme p (x ) = i =0 ai x i :
Schéma de Horner

p (x ) = a0 + x (a1 + x (a2 + x (...an )))

Coût algorithmique ?
 1/41

Chapitre 2 : Interpolation et approximation

...
 2/41
A partir d’un ensemble de points (xi , yi )

i =0:n , on recherche dans un
espace de fonctions,

Interpolation : une fonction s qui interpole les noeuds (xi , yi ),

s (xi ) = yi , ∀i = 0 : n
Approximation au sens des moindres carrés : une fonction s qui
minimise l’énergie,
n
X
|yi − s (xi )|2
i =0
 3/41
Exemples d’interpolation :

fig: Interpolation polynomiale, linéaire par morceaux et spline

Exemples d’approximation au sens des moindres carrés :

fig: Régression linéaire, quadratique

 4/41

1 Interpolation polynomiale

2 Interpolation polynomiale par morceaux

3 Approximation au sens des moindres carrés

 5/41

Interpolation polynomiale :

Définition
On note Pn l’ensemble des polynômes réels de degré n
n o
Pn = p (x ) ; p (x ) = a0 + a1 x + a2 x 2 + ... + an x n

avec ai ∈ R

On recherche un polynôme pn ∈ Pn tel que pour tout i = 0 : n,

pn (xi ) = yi .

Théorème
Il existe un unique polynôme pn ∈ Pn qui interpole les noeuds (xi , yi )

i =0:n .
 6/41

Polynôme d’interpolation de Vandermonde :

On souhaite trouver un polynôme de degré 2 qui interpole les noeuds

(−1, 2), (0, 3), (1, 6) .


On cherche le polynôme p2 sous la forme p2 (x ) = a0 + a1 x + a2 x 2 tel que

p2 (−1) = 2, p2 (0) = 3 et p2 (1) = 6. Alors
  



 p2 (−1) = 2 


 a0 − a1 + a2 = 2 


 a0 = 3
  
p2 (0) = 3 =⇒  = 3 =⇒  a1 = 2
  

  a0 

 
 

p2 (1) = 6
 a0 + a1 + a2 = 6
 a2 = 1


La solution du problème est donc p2 (x ) = 3 + 2x + x 2

 7/41

Polynôme d’interpolation de Vandermonde :

Plus généralement,
n Pno est un espace vectoriel dont la base canonique
s’écrit 1, X , X , ..., X n et
2

n
X
pn (x ) = ak x k .
k =0

Alors
 



 p n ( x0 ) = y0



 a0 + a1 x0 + a2 x02 + ... + an x0n = y0
a0 + a1 x1 + a2 x1 + ... + an x1
  2 n
pn (x1 ) = y1 = y1


 


 
.. =⇒ ..
. .


 




 


 

pn (xn ) = yn
 a0 + a1 xn + a2 x 2 + ... + an x n


= yn
n n
 8/41

Polynôme d’interpolation de Vandermonde :





p (x0 ) = y0  1 x0 x02 · · ·

x0n

  a0
 
  y0


  1 x x 2 · · ·
p (x1 ) = y1 x1n a1 y1

     

 1 1
.. =⇒  . ..  =  .. 
     



 . .
 .  
  .   . 


1 xn xn2 · · · xnn
    

p (xn ) = yn
 an yn
Les coefficients ai du polynôme d’interpolation s’obtiennent donc en
résolvant le système linéaire suivant
Ba = y ,
où (B )i ,j = xji −1 , (a )i = ai et (y )i = yi .

Propriété
Le déterminant de la matrice B vérifie
Y
det (B ) = (xj − xi ).
0≤i <j ≤n
 9/41

Polynôme d’interpolation de Vandermonde :

Le polynôme d’interpolation pn s’obtient donc par la résolution du système

linéaire Ba = y ! Mais en pratique
B est une matrice mal conditionnée !
Coût de la résolution du système linéaire en O (n3 ) !

Idée : Exprimer le polynôme pn dans une autre base de Pn et pour

laquelle, la décomposition de pn est explicite !
 10/41

Polynôme d’interpolation de Lagrange

On souhaite trouver le polynôme d’interpolation associé aux noeuds
(−1, 2), (0, 3), (1, 6) .

On introduit alors les trois polynômes L0 , L1 et L2 définis par
1 1
L0 (x ) = x (x − 1), L1 (x ) = 1 − x 2 , et L2 (x ) = x (1 + x ),
2 2
qui vérifient




 L0 (−1) = 1, L0 (0) = 0, L0 (1) = 0,
L1 (−1) = 0, L1 (0) = 1, L1 (1) = 0,





L2 (−1) = 0, L2 (0) = 0, L2 (1) = 1,



Exercice
Montrer que la famille {L0 , L1 , L2 } est libre et expliciter la décomposition
de p2 dans cette base
 11/41

Polynôme d’interpolation de Lagrange

n o
Plus généralement, on introduit une base de Pn notée Lnk , où les
k = 0 :n
polynômes Lnk sont définis par la propriété suivante :

1 si k = i

Lnk (xi ) = 

0 sinon


Exercice
Montrer que le polynôme Lnk s’explicite sous la forme
Qn
i =0,i ,k (x − xi )
Lnk (x ) = Qn .
i =0,i ,k (xk − xi )

Exercice
Montrer que le polynôme d’interpolation pn au noeuds (xi , yi )

i = 0 :n vérifie
n
X
pn = yi Lni
i =0
 12/41

fig: Interpolation polynomiale de la fonction f (x ) = 1

1+x 2
sur [−1, 1] avec
respectivement n = 6, n = 14 et n = 24

fig: Interpolation polynomiale de la fonction f (x ) = 1

1+x 2
sur [−5, 5] avec
respectivement n = 6, n = 14 et n = 24
 13/41

Erreur d’interpolation

On suppose que f : [a , b ] → R est une fonction suffisamment régulière.

On vient d’expliquer comment obtenir le polynôme d’interpolation pn au

noeuds (xi , f (xi )) i =0:n .


On souhaite maintenant donner une estimation de l’erreur E (x ) commise

entre pn et f :
E (x ) = f (x ) − pn (x ).
 14/41

Erreur d’interpolation

Théorème
Soit f ∈ C n+1 [a , b ] avec a = x0 < x1 < x2 ... < xn = b. Alors ∀x ∈ [a , b ], il
existe ξx ∈ [a , b ] tel que
n
f (n+1) (ξx ) Y
f (x ) − p (x ) = (x − xi )
(n + 1)! i =0

Démonstration : Appliquer le théorème de Rolle à la fonction φ

définie par

(f (x ) − pn (x )) ni=0 (t − xi )
Q
φ(t ) = f (t ) − P (t ) − Qn
i =0 (x − xi )
 15/41

Erreur d’interpolation

Si les noeuds sont équirépartis sur [a , b ] alors

 
Y n  n!
 (x − xi ) ≤ h n+1 ,
 i =0  4

où h = (b − a )/n, et
!n+1
b −a 1
 
|E | ≤ sup f (n+1) (x ) .
x ∈[a ,b ] n 4(n + 1)

A noter, rien n’assure que E → 0 lorsque n → ∞ !

Exercice
Montrer le résultat précédent
 16/41

Noeuds d’interpolation de tchebychev

Idée : Optimiser la position des noeuds xi ∈ [a , b ] afin de minimiser
 
Y n 
max  (x − xi )
x ∈[a ,b ] 
i =0

Solution :

a+b b −a
xi = + x̂i ,
2 2
où
π 2i + 1
!
x̂i = cos
2 n+1 x0 x1 x2 x3 x4 x5 x6 x7 x 8 x 9 x 10

On obtient alors l’estimation

1 (b − a )n+1  
sup f (n+1) (x ) ,
 
max |E (x )| ≤
(n + 1)! 22n+1 x ∈[a ,b ]

et limn→∞ E = 0.
 17/41

fig: Interpolation polynomiale de la fonction f (x ) = 1+1x 2 sur [−5, 5] avec

respectivement n = 6, n = 14 et n = 24 et noeuds equirépartis

fig: Interpolation polynomiale de la fonction f (x ) = 1+1x 2 sur [−5, 5] avec

respectivement n = 6, n = 14 et n = 24 et noeuds d’interpolation de tchebychev
 18/41

1 Interpolation polynomiale

2 Interpolation polynomiale par morceaux

3 Approximation au sens des moindres carrés

 19/41

Interpolation polynomiale par morceaux

Motivation :
L’interpolation polynomiale a tendance à produire des oscillations
lorsque le degré des polynômes est trop élevé : c’est le phénomène
de Runge ...
-> Utiliser des polynômes de degré plus petit sur des sous-intervalles
de [a , b ]
 20/41

Interpolation polynomiale par morceaux

Noeuds d’interpolation a = x0 < x1 < x2 < ... < xn = b
On note Ii = [xi −1 , xi ] et hi = |Ii | = |xi − xi −1 | pour tout i = 1 : n
On recherche une interpolation de f dans les espaces Shk ,r définis par

Skh ,r ([a , b ]) = s ∈ C r ([a , b ]) ; s|Ii ∈ Pk (Ii ), ∀i = 1 : n



p3
p2
p4
f(x 2)
f(x 3) p1
f(x 4)

f(x 1)

f(x 0)

I1 I2 I3 I4
x0 x1 x2 x3 x4
 21/41

Interpolation linéaire par morceaux

Espace Skh ,r ([a , b ]) avec k = 1 et r = 0
Espace vectoriel de dimension n + 1
2n coefficients associés aux n polynômes (sur chaque intervalle Ii )
n − 1 contraintes pour la continuité de s : pi (xi ) = pi +1 (xi ).
Solution explicite
! !
x − xi x − x i −1
s|iI (x ) = yi −1 + yi
xi −1 − xi x i − x i −1
 22/41

Exercice
Avec φi définie par
 x −x




i −1
xi −xi −1 si x ∈ [xi −1 , xi ]
 x i + 1 −x
φi (x ) =  x ∈ [xi , xi +1 ]

 xi +1 −xi si


0

sinon

Montrer que {φi }i =0:n forme une base de S1h,0 ([a , b ]) et que

n
X
s= yi φi .
i =0

Exercice
Soit f ∈ C 2 ([a , b ]), montrer que
  1
max |hi |2 max |f 00 (x )|

max s (x ) − f (x ) ≤
 
x ∈[a ,b ] 8 i =1:n x ∈[a ,b ]
 23/41

Interpolation Spline

Motivation : On souhaite plus de régularité...

En particulier, on recherche parmi toutes les fonctions C 2 ([a , b ]) qui
interpolent f aux noeuds {xi }i =0:n , celle qui minimise l’énergie élastique
Z b  00 2
J (s ) = s (x ) dx .
a

Definition
La fonction s définie sur [a , b ] est une spline cubique (par rapport à la
subdivision a = x0 < x1 < ... < xn ) si s ∈ S3,2 ([a , b ]). Si de plus
s 00 (a ) = s 00 (b ) = 0, alors s est une spline cubique naturelle.
 24/41

Interpolation Spline

L’espace vectoriel S3,2 ([a , b ]) est de dimension n + 3

n polynômes de degré 3 → 4n coefficients
contrainte continuité → n − 1 conditions
contrainte dérivée continue → n − 1 conditions
contrainte dérivée seconde continue → n − 1 conditions
Il faut donc imposer 2 conditions supplémentaires pour espérer avoir
l’unicité de l’interpolation spline

Théorème
Il existe une spline cubique qui interpole les noeuds (xi , yi ) i =0:n . Si de

plus les coefficients s 00 (a ) et s 00 (b ) sont fixés, alors cette spline est
unique.
 25/41

Une expression explicite de la spline interpolante s’obtient en résolvant un

système linéaire de taille n + 1. Pour simplifier les calculs, nous
supposons que la répartition des noeuds est uniforme, c’est à dire

hi = |xi − xi −1 | = h = (b − a )/n

La restriction de s à Ii est un polynôme de degré 3 noté pi :

s|Ii = pi ∈ P 3 (Ii )

p3
p2
p4
f(x 2)
f(x 3) p1
f(x 4)

f(x 1)

f(x 0)

I1 I2 I3 I4
x0 x1 x2 x3 x4
 26/41

Sur chaque intervalle Ii , la spline vérifie

− x i −1 ) 3 − xi ) 3 h2
!
(x (x (x − xi −1 )
s|Ii (x ) = pi (x ) = yi00 − yi00−1 + yi − yi00
6h 6h 6 h
h 2 00 (x − xi )
!
− y i −1 − y
6 i −1 h

où les coefficients yi00 s’obtiennent par la résolution du système linéaire

suivant
  s 00 (a ) i
  00 
1 0 0 0 ··· 0   y0
 
  
  6 y2 −y1 − y1 −y0 /h
h
 1 4 1 0 ··· 0   y 00 
  1 h h 
.. .. .. ..   .. ..
   
. . . .   . .
  
 0 0   
  .  =  ..
.. .. ..

  ..

. . . .

0 ··· 0
   
   00   h y −y
  6 n n−1 − yn−1 −yn−2 /h
i
0 ··· 0 1 4 1   y
 
   n −1 h h

yn00
 
0 ··· ··· ··· 0 1 s 00 (b )

 27/41
Continuité de la dérivée seconde de s
La dérivée seconde du polynôme pi00 est un polynôme de degré 1, qu’on
exprime sous la forme

( x − xi ) ( x − x i −1 )
pi00 (x ) = yi00−1 + yi00 ,
(xi −1 − xi ) (xi − xi −1 )

où les coefficients yi00 = s 00 (xi ) sont pour l’instant indéterminés à

l’exception de y000 = s 00 (a ) et yn00 = s 00 (b ).

En intégrant 2 fois, on obtient

(x − xi −1 )3 (x − xi )3
pi (x ) = yi00 − yi00−1 + qi (x ),
6h 6h
où qi (x ) ∈ P1 (Ii ).
Dans la suite, on recherche qi sous la forme

(x − xi −1 ) (x − xi )
qi (x ) = αi − βi .
h h
 28/41

Conditions d’interpolation s (xi ) = yi

On rappelle que

(x − xi −1 )3 (x − xi )3 (x − xi −1 ) (x − xi )
pi (x ) = yi00 − yi00−1 + αi − βi .
6h 6h h h
Alors
2
 
pi (xi )

 = yi αi

 = yi − h6 yi00
=⇒ 2

pi (xi −1 ) = yi −1
 βi

 = yi −1 − h6 yi00−1
Et

− xi −1 )3 − xi )3 h2
!
(x (x (x − xi −1 )
pi (x ) = yi00 − yi00−1 + yi − yi00
6h 6h 6 h
h 2 00 (x − xi )
!
− yi −1 − y
6 i −1 h
 29/41

Continuité de la dérivée première de s

Avec
 2 2
 

 p 0 (x ) = y 00 (x −xi −1 ) − y 00 (x −xi ) + (yi −yi −1 ) + h (y 00 − y 00 )
2h i −1 2h
 h 6 i −1
 i

 i i

 0 00 (x −xi )2 00 (x −xi +1 )2 (y i + 1 −y i ) h 00 00
( ) = + + ( )


 p
 i +1 x y i +1 2h
− y i 2h h 6 iy − yi +1

Et   

 p 0 (x ) = y 00 h + (yi −yi −1 ) + h (y 00 − y 00 )
i i 2
 i  h 6 i −1 i



 0 00 h (y i + 1 −y i ) h 00 00 )
pi +1 (xi ) = −yi 2 + + (



h 6 i y − y i +1

L’égalité pi0 (xi ) = pi0+1 (xi ) implique que

 yi+1 −yi yi −yi −1 
−
= 6  h h
 
yi00−1 + 4yi00 + yi00+1 
h
 30/41
En conclusion,
les coefficients yi00 s’obtiennent donc en résolvant le système linéaire
suivant
  s 00 (a ) i
  00 
1 0 0 0 ··· 0   y0
 
  
  6 y2 −y1 − y1 −y0 /h
h
 1 4 1 0 ··· 0   y 00 
  1 h h 
.. .. .. ..   .. ..
   
. . . .   . .
  
 0 0   
  .  =  ..
.. .. ..

  ..

. . . .

0 ··· 0
   
   00   h y −y
  6 n n−1 − yn−1 −yn−2 /h
i
0 ··· 0 1 4 1   y
 
   n −1 h h

yn00
 
0 ··· ··· ··· 0 1 s 00 (b )


Et sur chaque intervalle Ii , la spline vérifie

− x i −1 ) 3 − xi ) 3 h2
!
(x (x ( x − x i −1 )
s|Ii (x ) = yi00 − yi00−1 + yi − yi00
6h 6h 6 h
h 2 00 (x − xi )
!
− y i −1 − y
6 i −1 h
 31/41

Exemple d’interpolation spline 2d

Théorème (Erreur d’interpolation)

Soit f ∈ C 4 [a , b ] avec a = x0 < x1 < x2 ... < xn = b. Alors
  1  
max f (x ) − s (x ) ≤ h 4 max f (4) (x )
x ∈[a ,b ] 2 x ∈[a ,b ]
 32/41

Théorème (Minimisation d’énergie)

Soit f ∈ C 2 [a , b ] telle que f (xi ) = yi . Alors la spline cubique naturelle s
vérifie Z b Z b
(s 00 (x ))2 dx ≤ (f 00 (x ))2 dx
a a

Démonstration : Montrer que

Z b 
f 00 (x )s 00 (x ) − s 00 (x )2 dx = 0,
a

puis en déduire le résultat.

 33/41

1 Interpolation polynomiale

2 Interpolation polynomiale par morceaux

3 Approximation au sens des moindres carrés

 34/41

Approximation au sens des moindres carrés

Motivation :
Lorsque les données sont bruitées, l’interpolation de f aux noeuds
(xi , yi ) i =0:n n’est pas efficace. Il est alors préférable de rechercher une

fonction dans un espace vectoriel plus petit qui minimise une erreur L 2 par
rapport aux attaches au données.
 35/41

Approximation au sens des moindres carrés

Plus précisément, l’idée consiste à introduire un ensemble de fonctions

linéairement indépendantes
n o
g0 , g1 , ...gp

et parmi l’ensemble des combinaisons linéaires des gk , on recherche la

Pp
fonction g = k =0 ak gk qui minimise l’énergie
 p
2
n 
X X 
φ(a0 , a1 ..., ap ) = y −
 i ak gk (xi ) .
i =0 k =0

 36/41

Régression linéaire
Pour la régression linéaire, on utilise g0 = 1 et g1 = x et on minimise
n 
y − (a + a x )2 .
X
φ(a0 , a1 ) = i 0 1 i
i =0

Le minimum de φ est alors atteint si ∇φ = 0 :


∂a0 φ(a0 , a1 ) = −2 i =0 (yi − (a0 + a1 xi )) = 0
 Pn

∂a1 φ(a0 , a1 ) = −2 Pn xi (yi − (a0 + a1 xi )) = 0,


i =0

où encore
P  n
X
n
(n + 1)a0 + i =0 xi a1 = yi
i =0
 n  n
X  P  X
 xi  a0 + ni=0 xi2 a1 = xi yi
i =0 i =0
 37/41

Régression linéaire

Les coefficients a0 et a1 s’identifient donc à

 P 
1 n Pn 2 Pn Pn
a0 =

 D  i =0 yi i =0 xi − i = 0 xi i = 0 xi yi

i = 0 yi ,

a1 = 1 Pn Pn Pn

D (n + 1) i =0 xi yi − i = 0 xi

Pn 2 Pn 2
où D = (n + 1) i =0 xi −( i =0 xi ) .
 38/41

Définition du gradient d’une fonctionnelle

Soit φ une fonctionnelle de Rp → R.

La différentielle de φ en x dans la direction y est noté D φ(x )(y )

φ(x + ty ) − φ(x )
D φ(x )(y ) = lim .
t →0 t
Le gradient de φ en x par rapport à un produit scalaire < ., . > donné
est alors défini par

< ∇φ(x ), y >= D φ(x )(y ), ∀y ∈ Rp .

Pp
Dans le cas du produit scalaire canonique < x , y >= xy,
i =0 i i

(∇φ(x ))i = ∂xi φ(x ).

 39/41

Définition du gradient d’une fonctionnelle

Exercice
Calculer le gradient des fonctionnelles suivantes




 φ1 (x ) =< x , b >
φ2 (x ) =< Ax , b >





φ3 (x ) = 1 < Ax , Ax >


2
 40/41

Moindres carrés, cas général

On introduit p + 1 fonctions {gi }i =0:p et on recherche le minimum de

l’énergie
 p
2
n 
X X 
φ(a0 , a1 ..., ap ) = y −
 i ak gk (xi ) .
i =0 k =0

La fonctionnelle φ s’exprime sous la forme

φ(a ) =< y − Ba , y − Ba >,

avec

a0  g0 (x0 ) g1 (x0 ) · · · gp (x0 )  y0 

     
  
 a1   g (x ) g (x ) · · · gp (x1 )   y 
 0 1 1 1  1 
..  , B =  .. .. .. ..  , y =  ..  ,
  
a = 

 .   . . . .   . 
    
ap g0 (xn ) g1 (xn ) · · · gp (xn ) yn
 41/41

Moindres carrés, cas général

Le gradient de φ s’identifie à
 
∇φ(a ) = 2 B t Ba − B t y .

Les coefficients a sont donc solution du système linéaire suivant

B t Ba = B t y .

Exercice
nEcrire la omatrice B dans le cas d’approximation polynomiale :
gk = x k .
k =0:p

Exercice
Retrouver l’expression de a0 et a1 dans le cas d’une régression linéaire.
 1/17

Chapitre 3 : Intégration numérique

...
 2/17

Objectif : Proposer des méthodes numériques pour calculer des

intégrales :
Z b
I(f ) = f (x )dx .
a

Motivations : Exemple, évaluer des fonctions qui ne s’expriment pas à

l’aide de fonctions usuelles
Z t Z t
2
φ(t ) = e −πx dt , ln(t ) = 1/xdx .
0 1

Point de vue : Introduire une discrétisation de l’intervalle [a , b ],

a = x0 < x1 < ... < xn = b et regarder des approximations de I(f ) de la
forme
n
X
In (f ) = αi f (xi ).
i =0

Il s’agit de bien choisir les coefficients αi pour obtenir une bonne

approximation de I(f ) !
 3/17

Première exemple : méthode des rectangles

a x1 x2 x3 x4 a x1 x2 x3 x4 a x1 x2 x3 x4
b b b

Quelle est la meilleure approximation de l’intégrale I(f ) ?

 (b −a ) Pn−1



 In (f ) = n i =0 f (xi ),
 (b −a ) Pn
In (f ) = f (x )

n  i =1 i P


 (b −a )

f (a )/2 + in=−11 f (xi ) + f (b )/2

In (f ) =

n
 4/17

1 Formules de Newton-Cotes

2 Formules de Newton-Cotes composite

 5/17

Principe des formules de Newton-Cotes : intégrer des

interpolations polynomiales de f .

L’interpolation polynomiale de Lagrange s’écrit

n
X
f (x ) ' pn (x ) = f (xi )Lni (x ),
k =0

Une approximation de I(f ) s’obtient avec

Z b n Z
X b !
In (f ) = pn (x )dx = Lni (x )dx f (xi ),
a i =0 a
 6/17

Les coefficients αi s’identifient aux intégrales

Z b Z b Y ( x − xj )
αi = Lni (x )dx = dx
a a j ,i
(xi − xj )

Dans le cas d’une discrétisation uniforme, on a de plus

b n Z n Y n
(x − xj )
Z Y (t − j )
αi = dx = h dt ,
a j =0,j ,i
(xi − xj ) 0 j =0,j ,i i − j

avec h = (b − a )/n.
 7/17

Formules du trapèze
Avec n = 1 et x0 = a, x1 = b, la formule d’approximation I1 s’identifie à
b −a
I1 ( f ) = [f (a ) + f (b )]
2

Exercice
Montrer cette formule

Exercice
Calculer une approximation de ln(2)

Théorème
Soit f ∈ C 2 ([a , b ]), alors il existe η ∈ [a , b ] tel que

(b − a )3
I(f ) − I1 (f ) = − f 00 (η).
12
 8/17

La démonstration de ce théorème repose essentiellement sur la formule

d’erreur d’interpolation polynomiale et le théorème de la moyenne :

Théorème (Théorème de la moyenne)

Soit g une fonction continue et h une fonction de signe constant sur [a , b ].
Alors il existe η ∈ [a , b ] tel que
Z b Z b
g (x )h (x )dx = g (η) h (x )dx .
a a

Exercice
Démontrer le théorème de la moyenne.
 9/17
Démonstration :
On rappelle que pour tout x ∈ [a , b ], il existe ξx ∈ [a , b ] tel que
1 00
f (x ) − p1 (x ) = f (ξx )(x − a )(x − b ).
2
Alors Z b
1
I(f ) − I1 (f ) = f 00 (ξx )(x − a )(x − b )dx .
2 a
On applique le théorème de la moyenne avec g (x ) = f 00 (ξx ), et
h (x ) = (x − a )(x − b ) : on en déduit qu’il existe η ∈ [a , b ] tel que
Z b
1
I(f ) − I1 (f ) = f 00 (η) (x − a )(x − b )dx .
2 a

Au final, avec Z b
1
(x − a )(x − b )dx = − (b − a )3 ,
a 6
et
(b − a )3
I(f ) − I1 (f ) = − f 00 (η)
12
 10/17

Formule de Simpson
Avec n = 2 et x0 = a, x1 = (a + b )/2, x2 = b, la formule d’approximation
I2 s’identifie à
" ! #
b −a a+b
I2 ( f ) = f (a ) + 4f + f (b )
6 2

Exercice
Montrer cette formule

Exercice
Calculer une approximation de ln(2)

Théorème
Soit f ∈ C 4 ([a , b ]), alors il existe η ∈ [a , b ] tel que

1 (4)
I(f ) − I2 (f ) = −(b − a )5 f (η).
2880
 11/17

Intuition de la preuve :
1) Erreur d’interpolation polynomiale

1 (3)
f (x ) − p2 (x ) = f (ξx )(x − a )(x − (a + b )/2)(x − b )
3!
2) Théorème de la moyenne sur [a , (a + b )/2] et [(a + b )/2, b ] :
!4 
1 b −a 
I(x ) − I2 (f ) = f (3) (η1 ) − f (3) (η2 )
4 2

avec η1 ∈ [a , (a + b )/2] et η2 ∈ [(a + b )/2, b ]

3) Théorème de la moyenne

f (3) (η1 ) − f (3) (η2 ) = (η2 − η1 )f (4) (η)

 12/17

1 Formules de Newton-Cotes

2 Formules de Newton-Cotes composite

 13/17

Formules de Newton-Cotes composites

Motivation : Augmenter la précision des intégrales sans augmenter le

degré des polynômes (phénomène de Runge) !.

=⇒ Utiliser des approximations polynomiales par morceaux.

On utilise une discrétisation de l’intervalle [a , b ], a = x0 < x1 < ... < xn = b

pour découper l’intégrale I(f )
Z b X Z xi !
f (x )dx = f (x )dx
a i =1:n xi −1

puis on applique les méthodes précédentes sur chacune des nouvelles

intégrales.
 14/17

Formule des trapèzes composites

Avec hi = |xi − xi −1 |, on a
Z xi
hi
f (x )dx ' (f (xi −1 + f (xi )) .
xi −1 2
et on pose
n
X hi
I1c (f , n) = (f (xi −1 ) + f (xi )) .
2
i =1
Dans le cas uniforme
 n −1

h  X 
I1c (f , h ) = f (a ) + 2 f (xi ) + f (b )
2 
i =1

Théorème
Si f ∈ C 2 ([a , b ]) alors

  h2
I(f ) − I1c (f , h ) ≤ (b − a ) max |f 00 (x )|
12 x ∈[a ,b ]
 15/17

Formule Simpson composites

Avec hi = |xi − xi −1 |, on a
Z xi
hi x i −1 + x i
   
f (x )dx ' f (xi −1 + 4f + f ( xi ) .
xi −1 6 2
et on pose
n
X hi
 x
i −1 + xi  
I2c (f , h ) = f (xi −1 ) + 4f + f ( xi ) .
6 2
i =1
Dans le cas uniforme
 n −1 n

h  X X x + x  
i − 1 i
I2c (f , n) = f (a ) + 2 f (x i ) + 4 f + f ( b )  .
6  2 
i =1 i =1

Théorème
Si f ∈ C 4 ([a , b ]) alors

  h4
I(f ) − I2c (f , h ) ≤ (b − a ) max |f (4) (x )|
2880 x ∈[a ,b ]
 16/17

Estimation d’erreur a posteriori

Motivation : On souhaite déterminer une estimation numérique de
l’erreur commise entre I(f ) et Ic2 (f , h ) pour un pas de discrétisation h.

Cas de la Formule Simpson composites :

On suppose que l’erreur entre I(f ) et Ic2 (f , h ) est de la forme

E (f , h ) = I(f ) − Ic2 (f , h ) ' ωf h 4

pour h suffisamment petit

En remarquant que
!4
h
E (f , h /2) = I(f ) − Ic2 (f , h /2) ' ωf ,
2
une approximation de E (f , h ) s’obtient avec
Ic2 (f , h /2) − Ic2 (f , h )
E (f , h ) '
1 − (1/2)4
 17/17

Estimation d’erreur a posteriori

Algorithme (Simpson composites)

Donnée : Précision  > 0,
E = 1, h = 1, I = Ic2 (f , h )
While E ≥ 
Ip = I,
h = h /2,
I = Ic2 (f , h )
E = (I − Ip )/(1 − 1/24 )
End
 1/49

Chapitre 4 : Résolution de systèmes linéaires

...
 2/49

Ce chapitre s’intéresse à différentes méthodes numériques pour la

résolution de systèmes linéaires

Ax = b .

Motivations : La plupart des méthodes numériques conduisent à la

résolution d’un système linéaire qui représente en général, la tache la plus
coûteuse d’un point de vue algorithmique !

2 types de méthodes :
Méthodes directes : factorisation de la matrice A et résolution exacte
du système linéaire.
Méthodes itératives : résolution approchée du système linéaire

Coût algorithmique : (à l’exception de matrices particulières) O (n3 )

 3/49

Cas simples

Matrices diagonales

x1 = b1 /a1,1
    
 a1,1 0   x1   b1  

.. ..  =  ..
  

xi = bi /ai ,i
    
. .   .  =⇒ 

   
    
0 an,n xn
 
bn xn = bn /an,n



Coût algorithmique -> O (n)

Matrices triangulaires

 a1,1 a1,2 · · · a1,n x1 b1

     
 0 a
     
 xn = bann
a2,n x2 b2

2,2
     

 bn−1 −an−1,n xn
n −1 =
x
 .. .. .. ..  =  ..  =⇒ 
     
 . . Pan−1,n−1
. .   . 
  
b − n a x

xk = k j=k +1 k ,j j
  


   
0 ··· 0 an,n xn bn ak ,k

Coût algorithmique -> O (n2 )

 4/49

1 Méthodes directes
Décomposition LU
Décomposition Cholesky

2 Conditionnement et erreur pour la résolution de systèmes linéaires

3 Méthodes itératives classiques

Algorithme du point fixe
Méthode du gradient
 5/49

Décomposition LU : idée générale

La décomposition LU consiste à factoriser A sous la forme

PA = LU ,

avec
 
 1 0 ··· 0   u1,1 u1,2 · · ·

u1,n

..
  
.
 0 u u2,n
2,2

 l2,1 1 
et U =  . .. ..
 
L =  .. ..  ..

. .

. .

0
  
  
ln,1 ··· l n −1 , n 1
 0 ··· 0 un,n

et P, une matrice de permutation (en particulier P t = P −1 ).

Motivation : La résolution du système Ax = b s’effectue en 2 étapes :


Ly = P t b

LUx = P t b =⇒ 

Ux = y

 6/49

Décomposition LU sans pivot

Hypothèse : Supposons que toutes les sous matrices principales de A
sont inversibles :
 
 a1,1 · · · a1,k 
 . ..
avec Ak =  ..

det (Ak ) , 0,  , ∀k = 1 : n,

.
 
ak ,1 · · · uk ,k

alors, il existe une telle décomposition LU avec P = Id , c’est à dire

A = LU .

Exemple en dimension 3 :
 
 a1,1 a1,2 a1,3  
a1,1 = det (A1 ) , 0


A =  a2,1 a2,2 a2,3  et
 
  
a1,1 a2,2 − a1,2 a2,1 = det (A2 ) , 0

a3,1 a3,2 u3,3
 7/49

Avec  
 1 0 0   a2,1
q2,1 =

G1 =  −q2,1 1 0  ,
  a11
avec


q3,1 =
 a3,1

−q3,1 0 1

a11

on a
 (1)
  

 a2,2 = a2,2 − a1,2 q2,1
 a1,1 a1,2 a1,3  

 (1)
a2,3 = a2,3 − a1,3 q2,1

(1) (1)

 ,
 
A (1)

= G1 A =  0 a2,2 a2,3 avec  (1)
 (1) (1)
 

 a3,2 = a3,2 − a1,2 q3,1
0 a3,2 a3,3 

a (1)



3,3
= a3,3 − a1,2 q3,1

On recommence avec A (1) :

  (1)
 1 0 0  a3,2
G2 =  0 1 0  , avec q3,2 =
 
(1)

0 −q3,2 1 a22
 8/49

 
 a1,1 a1,2 a1,3 
(1) (1) (2) (1) (1)
 , a3,3 = a3,3 − q3,2 a2,3 ,
 
G2 G1 A = A (2) =  0 a2,2 a2,3 où
(2)
 
0 0 a3,3

Au final  
 a1,1 a1,2 a1,3 
−1  (1) (1)
 = L U ,
 
A = (G2 G1 )  0 a2,2 a2,3
(2)
 
0 0 a3,3

où L = (G2 G1 )−1 = G1−1 G2−1 et donc

    
 1 0 0   1 0 0   1 0 0 
L =  q2,1 1 0   0 1 0  =  q2,1 1 0
    
  
q3,1 0 1 0 q3,2 1 q3,1 q3,2 1
 9/49

Algorithme (Décomposition LU sans Pivot )

For k = 1 : n − 1
For i = k + 1 : n
(k − 1 ) (k −1)
qi ,k = ai ,k /ak ,k
For j = k + 1 : n
(k ) (k −1) (k −1)
ai ,j = ai ,j − qi ,k ak ,j
End
End
End

Exercice
Calculer le coût de cet algorithme
 10/49

Exemple avec
   
 1 4 7   18 
A =  2 5 8  , et b =  24 
   
  
3 6 11 32
On pose q2,1 = a2,1 /a1,1 = 2 et q3,1 = a3,1 /a1,1 = 3
   
 1 4 7 18   1 4 7 18 
 .

8 24  −→  2 −3 −6 −12
 2 5
   
 
3 6 11 32 3 −6 −10 −22

(1) (1)
On pose q3,2 = a3,2 /a2,2 = (−6)/(−3) = 2 .
   
 1 4 7 18   1 4 7 18 
 2 −3 −6 −12  −→  2 −3 −6 −12  .
   
3 −6 −10 −22 3 2 2 2
 11/49

Au final,
    
 1 4 7   1 0 0   1 4 7 

A =  2 5 8  =  2 1 0   0 −3 −6 
   

3 6 11 3 2 1 0 0 2

La solution du système s’obtient en résolvant le système

        
 1 4 7   x1   18 
      x1   3 
 0 −3 −6   x2  =  −12  =⇒  x2  =  2 
        
0 0 2 x3 2 x3 1

En effet,     
 1 0 0   18   18 
 2 1 0   −12  =  24 
    
3 2 1 2 32
 12/49

Algorithme LU avec pivot

Motivation : La décomposition LU est instable lorsque les coefficients

(i )
ai +1,i +1 deviennent trop petits et ne s’applique plus si l’hypothèse des
sous-matrices principales inversibles n’est pas vérifiée.

En pratique, il est plus efficace d’utiliser une décomposition LU avec pivot

PA = LU .
 13/49

Exemple avec
   
 0 2 2   4 
A =  4 2 4  , et b =  10 
   
   
2 4 4 10
On commence par choisir le pivot, le plus grand élément de la première
colonne puis on procède comme pour l’algorithme LU sans pivot
     
 0 2 2 4   0 1 0   4 2 4 10 
 4 2 4 10  −→ P1 =  1 0 0  −→  0 2 2 4 
     
2 4 4 10 0 0 1 2 4 4 10

Puis
     
 4 2 4 10   1 0 0   4 2 4 10 
 0 2 2 4  −→ P2 =  0 0 1  −→  1/2 3 2 5 
     
1/2 3 2 5 0 1 0 0 2 2 4
 14/49

   
 4 2 4 10   4 2 4 10 
 1/2 3 2 5  −→  1/2 3 2 5  .
   
0 2 2 4 0 2/3 2/3 2/3

Alors, 



 x3 = 1

x2 = (5 − 2)/3 = 1





x1 = (10 − 4 − 2)/4 = 1


et     
 1 0 0   0 1 0   0 1 0 
P = P2 P1 =  0 0 1   1 0 0  =  0 0 1  ,
    
    
0 1 0 0 0 1 1 0 0
   
 1 0 0   4 2
4 
L =  1/2  ,

1 0  et U =  0 3 2
  
3 /2 1 0 0 2/3
  
0
 15/49

Exercice
Résoudre le système linéaire suivant en utilisant la décomposition LU
avec pivot
   
 3 1 6   2 
A =  2 1 3  , et b =  4  .
   
   
1 1 1 7
 16/49

1 Méthodes directes
Décomposition LU
Décomposition Cholesky

2 Conditionnement et erreur pour la résolution de systèmes linéaires

3 Méthodes itératives classiques

Algorithme du point fixe
Méthode du gradient
 17/49

Décomposition de Cholesky
La décomposition de Cholesky consiste à factoriser A sous la forme
A = LL t , où  
 l1,1 0 ··· 0 
 .. 
 l2,1 l2,2 . 
L =  .. ..

. .

0 


ln,1 ··· l n −1 , n ln,n


Motivation : Comme pour la décomposition LU, la résolution du système

Ax = b s’effectue en 2 étapes :

Ly = b

t
,

LL x = b =⇒ 
L t x = y


L’avantage par rapport à la factorisation LU est

de n’avoir qu’une seule matrice à stocker
de présenter un algorithme 2 fois plus rapide !
 18/49

Definition
Une matrice A est définie positive si

< x , Ax >≥ 0 et < x , Ax >= 0 =⇒ x = 0.

Un exemple de telles matrices sont les matrices à diagonale strictement

dominante !
Théorème
Si A est une matrice symétrique définie positive, alors il existe une unique
matrice L triangulaire inférieure avec coefficients diagonaux positifs telle
que
A = LL t .
 19/49

Exemple en dimension 2 :
!
a c
Soit A = une matrice définie positive.
c b
En particulier, l’hypothèse de positivité implique que

a =< e1 , Ae1 > > 0


ab − c 2 = det (A ) > 0



Avec
l12,1
! !
l1,1 0 l1,1 l2,1
L= et LL t = .
l2,1 l2,2 l2,1 l1,1 l12,1 + l22,2
On en déduit que  √



 l1,1 = a
c

2 , 1 = √a
l



 q
c2


l2,2 = b −

a
 20/49
Cas général

 l
 1,1 0 ··· 0
 l
  1,1 l2,1 ··· ln,1
  a
  1,1 a1,2 ··· a1,n


 ..   ..   .. 
 l2,1 l2,2 .   0 l2,2 .   a1,2 a2,2 .
     
 =  .

 .   .
.. .. ..

 .   .   .
 . .   . .   . .

0 ln−1,n an−1,n 
ln,1 ··· ln−1,n ln,n 0 ··· 0 ln,n a1,n ··· an−1,n an,n

La première ligne montre que a1,i = l11 li ,1 et

 √
l1,1 = a1,1


li ,1 = a1,i , pour i > 1


l1,1

La première colonne de L est maintenant connue.

La deuxième ligne montre que a2,i = `2,1 `i ,1 + `2,2 `2,i pour i ≥ 2 et
 q
`2,2 = a2,2 − `2,1


 2

2,1 `i ,1
`i ,2 = a2,i −` , pour i > 2



`2,2

Les deux premières colonnes de L sont maintenant connues.

 21/49

Cas général

 l
 1,1 0 ··· 0
 l
  1,1 l2,1 ··· ln,1
  a
  1,1 a1,2 ··· a1,n


 ..   ..   .. 
 l2,1 l2,2 .   0 l2,2 .   a1,2 a2,2 .
     
 =  .

 .   .
.. .. ..

 .   .   .
 . .   . .   . .

0 ln−1,n an−1,n 
ln,1 ··· ln−1,n ln,n 0 ··· 0 ln,n a1,n ··· an−1,n an,n

A la k ieme ligne : les k − 1 premières colonnes de L sont supposées

connues, alors
k
X
ak ,i = lk , j li , j , pour i ≥ k
j =1

et  q
`k ,k = ak ,k − j =1 `k2 ,l

 Pk −1

 
`i ,k = ak ,i − kj =−11 `k ,j `i ,j /`k ,k


 P

Les k premières colonnes de L sont connues.

 22/49

Algorithme (Décomposition de Cholesky)

For k = 1 : q
n
`k ,k = ak ,k − kj =−11 `k2,l
P

For i = k + 1 : n 
`i ,k = ak ,i − `k ,j `i ,j /`k ,k
Pk − 1
j =1

End
End

Exercice
Calculer le coût algorithmique de la décomposition de Cholesky
 23/49

Exercice
Appliquer l’algorithme de Cholesky pour factoriser la matrice
 
 4 2 1 
A =  2 4 2 
 
 
1 2 4
 24/49

1 Méthodes directes
Décomposition LU
Décomposition Cholesky

2 Conditionnement et erreur pour la résolution de systèmes linéaires

3 Méthodes itératives classiques

Algorithme du point fixe
Méthode du gradient
 25/49

Conditionnement d’une matrice

Motivation : Il s’agit ici de comprendre comment une erreur δA sur la

matrice A et une erreur δb sur le vecteur b influencent la solution x du
système linéaire
Ax = b .
En particulier, en notant x + δx la solution de

(A + δA ) (x + δx ) = b + δb ,

on souhaite contrôler l’erreur relative

kδx k
kx k
 26/49

Norme vectorielle

Definition
k.k est une norme sur Rn si ∀(x , y ) ∈ (Rn )2 et α ∈ R,
Positivité
kx k ≥ 0 et kx k = 0 ⇔ x = 0
Linéarité
kαx k = |α|kx k
Inégalité triangulaire
kx + y k ≤ kx k + ky k

Exemple  qP
n 2
k x k = i = 1 xi




 2

kx k1 = ni=1 |xi |
 P





kx k∞ = maxi =1:n {|xi |}

 27/49

Norme matricielle induites

Definition
Soit k.k une norme sur Rn . On note k|.k| la norme induite sur l’espace
vectoriel des matrices carrées de taille (n × n) définie par

k|A k| = max kAx k

kx k≤1

Exercice
Montrer que k|.k| est bien une norme sur l’espace des matrices carrées
(n × n )

Exemples :
√
k|A k|2 = λmax






k|A k|1 = maxj ni=1 |aij |
 P



k|A k|∞ = maxi Pn |aij |,


j =1

où λmax est la plus grande valeur propre de A t A .

 28/49

Conditionnement

Definition
Soit A ∈ Rn×n une matrice inversible. Le conditionnement de A par rapport
à une norme donnnée k.k est défini par

cond (A ) = k|A k| k|A −1 k|

Exercice
Calculer le conditionnement par rapport à la norme k.k∞ de la matrice

1.2969 0.8648
!
A=
0.2161 0.1441
 29/49

Théorème
Soit x et x + δx les solutions des systèmes linéaires

Ax = b


A (x + δx ) = b + δb



Alors
kδx k kδb k
≤ cond (A )
kx k kb k

Démonstration : Il suffit de remarquer que

1 k|A k|
Ax = b =⇒ kb k ≤ k|A k| kx k =⇒ ≤
kx k kb k
et
A δx = δb =⇒ kδx k ≤ k|A −1 k| kδb k
 30/49

Théorème
Soit x et x + δx les solutions des systèmes linéaires

Ax = b


(A + δA ) (x + δx ) = b + δb



Alors !
kδx k cond (A ) kδb k kδA k
≤ +
kx k (1 − k|A k| k|δA k|) kb k kA k
 31/49
Démonstration : Avec δx = δxb + δxA , où
A δxb = δb ,
il ressort que
A δxb = δb =⇒ kδxb k ≤ k|A −1 k| kδb k
et
(A + δA ) (x + δxb + δxA ) = b + δb =⇒ δxA = −A −1 δA (x + δx )
=⇒ kδxA k ≤ k|A −1 k| k|δA k| (kx k + kδx k)
On en déduit que
kδx k ≤ kδxA k + kδxb k ≤ k|A −1 k| kδb k + k|A −1 k| k|δA k| (kx k + kδx k)
Enfin avec
1 k|A k|
Ax = b =⇒ kb k ≤ k|A k| kx k =⇒ ≤
kx k kb k
on en déduit que
1  
kδx k ≤ k|A −1 k| kδb k + k|A −1 k| k|δA k| kx k
1− k|A −1 k| k|δA k|
 32/49

1 Méthodes directes
Décomposition LU
Décomposition Cholesky

2 Conditionnement et erreur pour la résolution de systèmes linéaires

3 Méthodes itératives classiques

Algorithme du point fixe
Méthode du gradient
 33/49

Méthode itératives classiques

Motivations :
Le coût de la résolution des méthodes directes est en O (n3 ) !
Dans de nombreuses situations, les matrices A sont creuses (peu
des coefficients non nuls) et la décomposition LU perd cette propriété
→ problème de place mémoire !
 34/49

1 Méthodes directes
Décomposition LU
Décomposition Cholesky

2 Conditionnement et erreur pour la résolution de systèmes linéaires

3 Méthodes itératives classiques

Algorithme du point fixe
Méthode du gradient
 35/49

Algorithme du point fixe

Le principe est de rechercher la solution x ∗ du système linéaire

Ax = b ,

comme le point fixe d’une application φ : Rn → Rn :

φ(x ∗ ) = x ∗ .

L’algorithme du point fixe construit alors la suite xk définie par


x0 ∈ Rn


xk +1 = φ(xk )



en espérant qu’avec de bonnes propriétés sur φ,

lim xk = x ∗
k →∞
 36/49

Exemple en dimension 1

Avec des fonctions φ de la forme φ(x ) = ax − b.

φ (x)
f(x) f(x)

φ (x)

x* x 2 x1 x0 x* x 0 x 1 x 2 x3

une condition nécessaire et suffisante pour que xk → x ∗ est

|a | < 1
 37/49

Contraction et existence de point fixe

Definition
Une application φ : Rn → Rn est une contraction pour la norme k.k s’il
existe λ < 1 tel que
kφ(x ) − φ(y )k ≤ λkx − y k

Théorème
Si φ est une contraction, alors φ admet un unique point fixe noté x ∗ .

Démonstration :
Existence : admise (suite de Cauchy)
Unicité : si x1 = φ(x1 ) et x2 = φ(x2 ) alors

kx1 − x2 k = kφ(x1 ) − φ(x2 )k

≤ λkx1 − x2 k
< kx1 − x2 k si x1 , x2
 38/49

Convergence de l’algorithme du point fixe

Théorème
Soit φ : Rn → Rn une contraction et x0 ∈ Rn . Alors la suite xk définie

x0 ∈ Rn


xk +1 = φ(xk ) ∀k > 0



converge vers l’unique point fixe de x ∗ et

1 Estimation a posteriori

1
k x ∗ − xk k ≤ kxk − xk +1 k
1−λ
2 Estimation a priori

λk
kx ∗ − xk k ≤ kx1 − x0 k
1−λ
 39/49

Démonstration :
1 Il suffit de remarquer que

kx ∗ − xk k ≤ kx ∗ − xk +1 + xk +1 − xk k ≤ kx ∗ − xk +1 k + kxk +1 − xk k
≤ kφ(x ∗ ) − φ(xk )k + kxk +1 − xk k
≤ λkx ∗ − xk k + kxk +1 − xk k

2 Avec l’inégalité (1), on a

1 1
kx ∗ − xk k ≤ k xk + 1 − xk k ≤ kφ(xk ) − φ(xk −1 )k
1−λ 1−λ
λ λk
≤ k x k − x k −1 k ≤ kx1 − x0 k
1−λ 1−λ
 40/49

Algorithme avec une précision 

Algorithme (Méthode du Point fixe)

Donnée : x0 ∈ Rn ,  > 0
x1 = φ(x0 ), k = 1
While kxk − xk −1 k ≤ c
xk +1 = φ(xk )
End

Exercice
Comment choisir c pour que l’algorithme s’arrête avec une précision de 
sur x ∗ ?
 41/49

Application pour la résolution de système linéaire

On s’intéresse à des applications φ de la forme

φ(x ) = M −1 (Nx + b ) .

où A = M − N.
Exercice
Montrer que si x ∗ est un point fixe de φ alors

Ax ∗ = b

L’application φ est une contraction si

k|M −1 N k| < 1.

On remarque de plus que xk +1 = φ(xk ) est solution de système linéaire

Mxk +1 = Nxk + b .

Idée : choisir M proche de A et facile à inverser !

 42/49

Méthode de Jacobi
La méthode de Jacobi utilise
M = D, N = −L − U .

Exercice
Montrer que  
n
1  X 
(xk +1 )i = b −
 i a ( x )
i ,j k j 

aii 
j =1,j ,i

Théorème
Si A est à diagonale strictement dominante, alors la méthode de Jacobi
converge

Démonstration : Avec M −1 N = −D −1 (L + U ), il ressort que

  
X  a 
i ,j
kM −1 N k∞ = maxj    < 1
aii
i ,j
 43/49

Méthode de Gauss-Seidel

La méthode de Gauss-Seidel utilise

M = D + L, N = −U .

Exercice
Montrer que
i −1 n
1
 X X 
(xk +1 )i = bi − aij (xk +1 )j − aij (xk )j (1)
aii
j =1 j =i +1

Théorème
Si A est symétrique définie positive, alors la méthode de Gauss-Seidel
converge

Démonstration : Admise ( regarder la norme k|M −1 N k|2 ...)

 44/49

1 Méthodes directes
Décomposition LU
Décomposition Cholesky

2 Conditionnement et erreur pour la résolution de systèmes linéaires

3 Méthodes itératives classiques

Algorithme du point fixe
Méthode du gradient
 45/49

Méthode du gradient

Hypothèse : A symétrique définie positive.

Le principe est de regarder la solution x ∗ du système linéaire

Ax = b

comme le minimum de l’énergie

1
J (x ) = < Ax , x > − < b , x > .
2
On construit alors une suite {xk } qui minimise J de la forme

x0 ∈ Rn


xk +1 = xk + αk dk

 ∀k ≥ 1,

où dk est une direction de descente et αk un pas de descente.

46/49
 47/49
Choix de dk :
La direction dk est l’opposé du gradient de J :

dk = −∇J (xk ) = −(Axk − b ) = b − Axk

Choix de αk :
On choisit αk comme le pas qui minimise J dans la direction dk :

J̃ (α) = J (xk + αdk ).

Avec

J̃ 0 (α) = < ∇J (xk + αdk ), dk >

= < A (xk + αdk ) − b , dk >
= α < Adk , dk > − < dk , dk >,

on en déduit que
< dk , dk >
J̃ 0 (αk ) = 0 ⇐⇒ αk =
< Adk , dk >
 48/49
Algorithme (Méthode du gradient)
Donnée : x0 ∈ Rn ,  > 0
While kAxk − b k ≤ k|A −1 k|
dk = b − Axk
αk =< dk , dk > / < Adk , dk >
xk +1 = xk + αk dk
End
Exemple avec A = (1, 0; 0, 5), b = (0, 0)0 et x0 = (1, 0.2)0
 49/49
Algorithme (Méthode du gradient conjugué)
Donnée : x0 ∈ Rn ,  > 0
While kAxk − b k ≤ k|A −1 k|
rk = b − Axk
dk = rk − i <k <<ddii ,,Ark>
P
d
Adi > i
αk =< dk , rk > / < Adk , dk >
xk +1 = xk + αk dk
End
Exemple avec A = (1, 0; 0, 5), b = (0, 0)0 et x0 = (1, 0.2)0
 1/27

Chapitre 5 : Résolution d’équations non-linéaires

...
 2/27
Ce chapitre s’intéresse aux méthodes numériques pour la résolution
d’équations non linéaires de la forme
f1 (x1 , . . . , xn )
 
 .. 
f (x ) =  .  = 0
fn (x1 , . . . , xn )
 

Exemples :
√
Calcul de la racine carré a : Premier algorithme, Babyloniens
Méthode des Moindres carrés non linéaires : Trouver la fonction de la
forme
g (α, β, x ) = αx β
qui minimise l’énergie
n 
β 2
X 
φ(α, β) = yi − αxi
i =1

La solution s’obtient en résolvant

∇φ = 0.
 3/27

1 Méthode de la Dichotomie

2 Méthode du point fixe

3 Méthode de Newton
 4/27

Méthode de la Dichotomie

Hypothèses
Soit f : [a , b ] → R, continue telle que f (a )f (b ) < 0

Principe
On regarde la fonction f au milieu
de [a,b], c = (a + b )/2,

 a c x* b
Si f (a )f (c ) < 0 =⇒ ∃x ∗ ∈ [a , c ]


Si f (a )f (c ) > 0 =⇒ ∃x ∗ ∈ [c , b ]


 5/27

Algorithme (Méthode de la dichotomie)

Donnée : aa = a, b0 = b, x0 = (a + b )/2 et  > 0
For k = 0 : N
If f (xk )f (ak ) < 0
ak +1 = ak
bk +1 = xk
else
ak +1 = xk
bk +1 = bk
end
xk +1 = (ak +1 + bk +1 )/2
End
 6/27

Exercice
Avec Ik = bk − ak , montrer que

b −a
Ik = .
2k

En déduire que le nombre d’itérations N nécessaire à l’algorithme pour

obtenir une précision de  sur la solution vérifie

ln((b − a )/)
N ≥
ln(2)

Exercice
√
Application pour l’estimation de 2 à une précision de 10−2 avec
f (x ) = x 2 − 2.
 7/27

Limite de l’approche de la Dichotomie

Cas pathologiques !
Vitesse de convergence linéaire.
Extension en dimension supérieure ?
 8/27

1 Méthode de la Dichotomie

2 Méthode du point fixe

3 Méthode de Newton
 9/27

Méthode du point fixe

Principe :
Transformer le problème f (x ) = 0 en un problème équivalent de la forme
g (x ) = x.
√
Exemple pour le calcul d’une racine carré a :
Cette racine peut s’obtenir en résolvant f (x ) = x 2 − a = 0, ou encore en
regardant les points fixes des fonctions g suivantes




 g1 (x ) = a /x

g2 (x ) = x 2 + x − a




g3 (x ) = (1 − 1/m)x + a /(mx ),



Question : Comment choisir une fonction g optimale ?

 10/27

Méthode du point fixe

Algorithme (Méthode du Point fixe)

Donnée : x0 ∈ Rn ,  > 0
for k = 0 : N
xk +1 = g (xk )
End

φ (x)
f(x) f(x)

φ (x)

x* x 2 x1 x0 x* x 0 x 1 x 2 x3
 11/27

Théorème (Version globale en dimension 1)

Soit g : [a , b ] → [a , b ] de régularité C 1 ([a , b ]) telle que

|g 0 (x )| ≤ λ < 1, ∀x ∈ [a , b ].

Alors g possède un unique point fixe x ∗ ∈ [a , b ] et xk → x ∗ lorsque

k → ∞. Plus précisément,
1 Estimation a posteriori

1
k x ∗ − xk k ≤ kxk − xk +1 k
1−λ
2 Estimation a priori

λk
kx ∗ − xk k ≤ kx1 − x0 k
1−λ
 12/27

Démonstration :
Existence : Avec h (x ) = g (x ) − x, remarquer que h (a )h (b ) ≤ 0
Unicité : Avec x1 = g (x1 ) et x2 = g (x2 ), remarquer que
  Z x2
|x1 − x2 | = g (x1 ) − g (x2 ) ≤ |g 0 (x )|dx ≤ λ|x1 − x2 |
x1

Vitesse de convergence : voir cours résolution systèmes linéaires.

Exercice
Montrer que le nombre d’itérations N permettant à l’algorithme d’obtenir
un précision de l’ordre de  sur la solution x ∗ du problème vérifie
 
(1−λ)
ln kx1 −x0 k
N ≥
ln(λ)
 13/27

Théorème (Version locale en dimension 1)

Soit x ∗ un point fixe d’une fonction g de régularité C 1 au voisinage de x ∗
telle que
|g 0 (x ∗ )| < 1.
Alors il existe un réel α > 0 tel que ∀x0 ∈ [x ∗ − α, x ∗ + α], l’algorithme du
point fixe converge vers le point fixe x ∗

Démonstration :
Appliquer le théorème global sur l’intervalle Iα = [x ∗ − α, x ∗ + α] avec α
suffisamment petit.
 14/27

Cas de plusieurs points fixes

Definition
Un point fixe x ∗ d’une fonction g est dit

attractif

 si |g10 (x ∗ )| < 1
|g10 (x ∗ )| > 1

répulsif
 si

Definition
Le bassin d’attraction d’un point fixe x ∗ est l’ensemble des valeurs initiales
x0 pour lesquelles xk tend vers x ∗ lorsque k tend vers l’infini

Exercice
Etudier le bassin d’attraction des points fixes de

g (x ) = x + (x 2 − 1)x .
 15/27

Vitesse de convergence
Le taux de convergence est donné par |g 0 (x ∗ )|.




 si |g 0 (x ∗ )| > 1 =⇒ l’algorithme diverge
si |g (x )| < 1 et g (x ) , 0 =⇒

 0 ∗ 0 ∗

 convergence linéaire

si |g 0 (x ∗ )| = 0


=⇒ convergence quadratique

Exercice (Cas où f (x ) = x 2 − a)
Etudier la convergence de l’algorithme du point fixe dans le cas des
fonctions suivantes




 g1 (x ) = a /x

g2 (x ) = x 2 + x − a




g3 (x ) = (1 − 1/m)x + a /(mx )



√
Déterminer une approximation de 2.
 16/27

Méthode du point fixe en dimension supérieure

Théorème (Version globale)

On introduit
B (x , R ) = y ∈ R n ; ky − x k ≤ R ,


Soit g : B (x , R ) → B (x , R ) est une contraction, i.e. ∃λ < 1 tel que

kg (x ) − g (y )k ≤ λkx − y k, ∀(x , y ) ∈ B (x , R ).

Alors g admet un point fixe x ∗ ∈ B (x , R ) et xk → x ∗

Théorème (Version locale)

Soit g : Rn → Rn de régularité C 1 au voisinage d’un point fixe x ∗ de g. Si
kDg (x ∗ )k < 1, alors il existe un réel α > 0 tel que ∀x0 ∈ B (x ∗ , α),
l’algorithme du point fixe converge vers x ∗ .
 17/27

Exercice
Déterminer les solutions de

x12 + x22 − 2
!
F (X ) = = 0.
x12 − x22

Montrer que ces solutions s’obtiennent aussi comme les points fixes
de l’application

 (1 − 1/m)x + 2−x22

 

1 mx1
G (X ) =  
 (1 − 1/m)x + x12 
2 mx2

Etudier la convergence de l’algorithme du point fixe dans le cas où

m = 4.
 18/27

1 Méthode de la Dichotomie

2 Méthode du point fixe

3 Méthode de Newton
 19/27

Méthode de Newton : cas dimension 1

Principe
Approcher la fonction f par sa tangente en (x0 , f (x0 )) :

1) p (x ) = f 0 (x0 )(x − x0 ) + f (x0 )

2) Résoudre p (x ) = 0 :

f (x0 )
=⇒ x1 = x0 −
f 0 (x0 )
x*
x2 x1 x0

Algorithme (Méthode du Point fixe)

f (x0 )
Donnée : x0 ∈ Rn ,  > 0, x1 = x0 − f 0 (x0 )
While kxk +1 − xk k ≤ /2,
f (xk )
xk +1 = xk − f 0 (xk )
End
 20/27

Méthode de Newton : cas dimension 1

Remarque
Choisir x0 suffisamment proche de
x∗

x0 x1 x2
x*

Remarque
Il s’agit d’une méthode de point fixe avec

f (x )
g (x ) = x −
f 0 (x )

En effet, si f 0 (x ∗ ) , 0 alors ”g (x ) = x ⇐⇒ f (x ) = 0”.

 21/27

Théorème (Cas de zéro simple)

Soit f ∈ C 2 ([a , b ]) et x ∗ ∈ [a , b ] tel que

f ( x ∗ ) = 0


f 0 (x ∗ ) , 0,



Alors il existe α > 0 tel que ∀x0 ∈ Iα = [x ∗ − α, x ∗ + α], l’algorithme de

Newton converge vers x ∗ .
Sa vitesse de convergence est de plus quadratique :

M
|xk + 1 − x ∗ | ≤ |xk − x ∗ |2 , ∀k ∈ N
2
supx ∈Iα |f 00 (x )|

où M = .
infx ∈Iα |f 0 (x )|

 22/27
Démonstration :
Idée : appliquer le théorème de convergence du point fixe local à la
f (x )
fonction g (x ) = x − f 0 (x ) .
Régularité de g : f étant C 2 et f 0 (x ∗ ) , 0, il existe un voisinage
Iα = [x ∗ − α, x ∗ + α] tel que g ∈ C 1 (Iα ).
Valeur de g 0 (x ∗ ) : en remarquant que

f 0 (x )2 − f (x )f 00 (x ) f (x )f 00 (x )
g 0 (x ) = 1 − = ,
f 0 (x )2 f 0 (x )2

on en déduit que |g 0 (x ∗ )| = 0.
Vitesse de convergence : La formule de Taylor montre que

1
0 = f (x ∗ ) = f (xk )−f 0 (xk )(x ∗ −xk )+ f 00 (ξ)(x ∗ −xk )2 , avec ξ ∈ [xk , x ∗ ]
2
et
f (xk ) 1 f 00 (ξ) ∗
xk +1 − x ∗ = (xk − x ∗ ) − = (x − xk )2
f 0 (xk ) 2 f 0 (xk )
 23/27

Exercice
2
√ le cas de la fonction f (x ) = x − 2 et
Ecrire l’algorithme de Newton dans
calculer une approximation de 2.

Exercice
Ecrire l’algorithme de Newton dans le cas de la fonction f (x ) = exp (x ) − 1.
Vérifier que la vitesse de convergence au voisinage de zéro est bien
quadratique.

Exercice
Même exercice pour la fonction f (x ) = sin(x )2 . La vitesse de convergence
au voisinage de zéro est-elle quadratique ? Pourquoi ?
 24/27

Méthode de Newton : Cas de racines multiples

Hypothèses : f (x ) = (x − x ∗ )m h (x ) avec m > 1 et h suffisamment
régulière.

Exercice
Montrer que l’algorithme de Newton converge mais que sa vitesse de
convergence est seulement linéaire.

Exercice
En appliquant l’algorithme de Newton à la fonction f̃ (x ) = f (x )1/m , montrer
que la variante de l’algorithme de Newton

f (xk )
xk +1 = xk − m ,
f 0 (xk )

admet cette fois-ci une vitesse de convergence quadratique. Application

dans le cas de la fonction f (x ) = sin(x )2 .
 25/27

Algorithme de Newton, cas dimension supérieure

Hypothèse : Soit f : Rn → Rn de régularité C 2 .
Principe : Le développement de Taylor montre que

f (x ) = f (x0 ) + Df (x0 )(x − x0 ) + O (kx − x0 k2 )

où  ∂ f1 ∂ f1
...

 ∂x ∂ xn

 1
Df (x ) =  ..
. .. 

. 
∂ fn ∂fn 

...

∂x1 ∂ xn

Algorithme (Algorithme de Newton)

Donnée : x0 ∈ Rn ,  > 0, x1 = x0 − Df −1 (x0 )f (x0 )
While kxk +1 − xk k ≤ /2,
xk +1 = xk − Df −1 (xk )f (xk )
End
 26/27

Théorème
Soit f : Rd → Rd de régularité C 2 et x ∗ un point fixe de f tel que Df (x ∗ ) soit
inversible. Alors il existe un réel α > 0 tel que ∀x0 ∈ B (x ∗ , α), l’algorithme
de Newton converge vers x ∗ . Sa vitesse de convergence est de plus
quadratique.

Exercice
Déterminer les solutions de

x12 + x22 − 2
!
F (X ) = = 0,
x12 − x22

puis calculer 2 itérations de l’algorithme de Newton en partant de

!
2
X0 =
2
 27/27

Exemple : Méthode des Moindres carrés non linéaires

On cherche la fonction de la forme

g (α, β, x ) = αx β

qui minimise l’énergie

n 
β 2
X 
φ(α, β) = yi − αxi .
i =1

Exercice
Montrer que les coefficients optimaux (α, β) vérifient
β 2 β
α
 Pn Pn 
F (α, β) =  i =1 (xi ) − i =1 yi xi  = 0,
 
β β
α2 ni=1 ln(xi )(xi )2 − α ni=1 yi xi ln(xi )
P P

Expliciter l’algorithme de Newton pour déterminer les coefficients

(α, β).
 1/33

Chapitre 6 : Méthodes numériques pour les équations

différentielles ordinaires (EDO)

...
 2/33
Ce chapitre s’intéresse aux méthodes numériques pour la résolution
d’équations différentielles ordinaires (EDO)

y 0 (t ) = f (t , y (t ))


(1)
y (t0 ) = y0 ∈ Rn ,



où f : R × Rn → Rn est supposée suffisamment régulière.

Exemple: Le pendule
θ l

ml θ00 = −α`θ0 − mgsin(θ)


θ(0) = θ0 , et θ0 (0) = θ0 ,


0
mg

Le problème peut se réécrire sous la forme (1) avec

θ θ0
! ! !
y2
y= , f (t , y ) = , et y (0) = ,
θ0 − mα y2 − g` sin(y1 ) θ00
 3/33

Existence + unicité de solution

Soit 
y 0 (t ) = f (t , y (t )), ∀t ∈ [a , b ]


(1)
y (t0 ) = y0 ∈ Rn , avec t0 ∈ [a , b ]



Théorème
Hypothèses : On suppose que f : [a , b ] × Rn → Rn est continue,
localement Lipschitzienne par rapport à la deuxième variable, i.e,
∀(u, v ) ∈ (R2 )2 et ∀t ∈ [a , b ], il existe un réel L > 0 tel que

kf (t , u) − f (t , v )k ≤ L ku − v k.

Alors (1) admet une unique solution C 1 sur Iα = [t0 − α, t0 + α].

 4/33

Exercice
Exemple avec  p
y 0 (t ) = |y (t )|,

 ∀t ∈ R

y (0) = 0

√
La fonction f (t , u) = |u| n’est pas lipschitzienne par rapport à u ! Montrer
que cette équation admet plusieurs solutions !

Exercice
Exemple avec 
y 0 (t ) = y (t )2 , ∀t ∈ R


y (0) = y0 > 0



La fonction f associée est bien continue et localement lipschitz. Montrer

que ce problème admet pour solution

1
y (t ) = , pour t < t0 + y0−1 .
y0−1 + t0 − t
 5/33

Discrétisation de l’équation

On suppose que le problème


y 0 (t ) = f (t , y (t )), ∀t ∈ [a , b ]


(1)
y (t0 ) = y0 ∈ Rn , avec t0 ∈ [a , b ]



admet une solution y suffisamment régulière sur [t0 , t0 + T ].

Discrétisation de l’intervalle [t0 , t0 + T ] avec N points :
T
Avec h = N , on pose tn = t0 + hn pour tout n = 0 : N.
Approximation numérique de y aux noeuds tn

yn ' y (tn ).

Schéma numérique de la forme

yn+1 = yn + h φ(tn , yn ).
 6/33

Convergence d’une méthode numérique

Définition (Consistance et ordre d’une méthode numérique)

Avec
y (tn+1 ) − y (tn )
τn + 1 ( h ) = − φ(tn , y (tn )),
h
La méthode numérique est dite consistante si limh →0 τn+1 (h ) = 0.
Elle est de plus d’ordre p si

kτn+1 (h )k ≤ Ch p .

Définition (Stabilité)
La stabilité d’un schéma numérique consiste à contrôler, pour un pas de
temps δt fixé, l’erreur
en = y (tn ) − yn ,
au cours des itérations.
 7/33

1 Méthodes d’Euler explicite-implicite

2 Méthodes d’ordre supérieur

 8/33

Méthode d’Euler explicite

Idée :
Z t1 Z t1
y (t1 ) − y (t0 ) = 0
y (t )dt = f (t , y (t )dt ' hf (t0 , y (t0 ))
t0 t0

Algorithme (Euler explicite)

Donnée : y0 ∈ Rn ,
for k = 0 : N − 1
yn+1 = yn + h f (tn , yn )
End
 9/33

Euler explicite : Consistance et ordre de la méthode

On remarque que

1
τn+1 (h ) = [y (tn+1 ) − y (tn ) − hf (tn , y (tn ))]
h
1
= [y (tn + h ) − y (tn ) − hy 0 (tn )]
h
1 2 1 00
= h y (η) avec η ∈ [tn , tn+1 ],
h 2
Avec M = supt ∈[t0 ,t0 +T ] |y 00 (t )| , on en déduit que


Théorème (Consistance méthode Euler explicite)

M
h.
|τn+1 (h )| ≤
2
La méthode d’Euler explicite est donc consistante et d’ordre 1.
 10/33

Euler explicite : Stabilité

Avec en = y (tn ) − yn , on montre que

en+1 = y (tn+1 ) − yn+1

Z tn +1
= y (tn ) + f (t , y (t ))dt − (yn + hf (tn , yn ))
tn
Z tn +1
= en + f (t , y (t ))dt − hf (tn , yn )
tn
"Z tn +1 #
= en + f (t , y (t ))dt − hf (tn , y (tn )) + h [f (tn , y (tn )) − f (tn , yn )]
tn
= en + h τn+1 (h ) + h [f (tn , y (tn )) − f (tn , yn )]

On rappelle que f est L -localement Lipschitzienne et donc

h2
|en+1 | ≤ (1 + Lh )|en | + M.
2
 11/33

Euler explicite : Stabilité

Au final,

h2
|en+1 | ≤ (1 + Lh )|en | + M
2
h2
≤ (1 + Lh )2 |en−1 | + M (1 + (1 + h L ))
2
n −1
h2 X
≤ (1 + Lh )n |e0 | + M (1 + h L )k
2
k =0
h 1
≤ exp (Ltn )|e0 | + M [exp (tn L ) − 1]
2 L
On en déduit le théorème suivant
Théorème (Stabilité algorithme Euler explicite)
Mh
ky (tn ) − yn k ≤ exp (LT )ky (t0 ) − y0 k + T exp (LT )
2
 12/33

Exemple: Le pendule
θ l

ml θ00 = −α`θ0 − mgsin(θ)


θ(0) = θ0 , et θ0 (0) = θ0 ,


0
mg

θ
!
Avec y = , cette équation se réécrit sous la forme
θ0
!
y2
y (t ) = f (t , y ) =
0
,
− mα y2 − g` sin(y1 )

Exercice
Écrire l’algorithme d’Euler explicite dans le cas particulier de l’équation du
pendule.
 13/33

Le pendule sans frottement : simulation numérique avec

la méthode d’Euler explicite
 14/33

Méthode d’Euler implicite

Idée :
Z t1 Z t1
y (t1 ) − y (t0 ) = 0
y (t )dt = f (t , y (t )dt ' hf (t1 , y (t1 ))
t0 t0

Algorithme (Euler explicite)

Donnée : y0 ∈ Rn ,
for k = 0 : N − 1
Déterminer yn+1 , solution de

F (y ) = y − h f (tn+1 , y ) − yn = 0

End

Exercice
Montrer que le schéma d’Euler explicite est bien consistant d’ordre 1
 15/33

Équation du pendule

Exemple: Le pendule
θ l

ml θ00 = −α`θ0 − mgsin(θ)


θ(0) = θ0 , et θ0 (0) = θ0 ,


0
mg

θ
!
Avec y = , cette équation se réécrit sous la forme
θ0
!
y2
y (t ) = f (t , y ) =
0
,
− mα y2 − g` sin(y1 )

Exercice
Écrire l’algorithme d’Euler implicite dans le cas particulier de l’équation du
pendule. On utilisera l’algorithme de Newton pour évaluer yn+1 en fonction
de yn .
 16/33

Le pendule sans frottement : simulation numérique avec

la méthode d’Euler implicite
 17/33

Comparaison Euler Explicite-implicite sur un système

raide

On s’intéresse au problème suivant


y 0 (t ) + Ly (t ) = 0,

 ∀t ∈ R+
y (0) = y0 ∈ Rn



Calculer la solution exacte de ce problème

Déterminer l’approximation numérique yn par la méthode d’Euler
Explicite.
Déterminer l’approximation numérique yn par la méthode d’Euler
implicite.
 18/33

1 Méthodes d’Euler explicite-implicite

2 Méthodes d’ordre supérieur

 19/33

Schémas numériques d’ordre supérieur à 1 pas

On suppose que le problème

y 0 (t ) = f (t , y (t )), ∀t ∈ [a , b ]


(1)
y (t0 ) = y0 ∈ Rn , avec t0 ∈ [a , b ]



admet une solution y suffisamment régulière sur [t0 , t0 + T ].

Discrétisation de l’intervalle [t0 , t0 + T ] avec N points :
T
Avec h = N , on pose tn = t0 + hn pour tout n = 0 : N.
Approximation numérique de y aux noeuds tn

yn ' y (tn ).

Schéma numérique de la forme

yn+1 = yn + h φ(tn , yn ).

Objectif : Trouver des schémas numériques d’ordre supérieur à 1.

 20/33

Méthode de Taylor
Idée :
h2
y (tn+1 ) = y (tn + h ) = y (tn ) + y 0 (tn )h + y 00 (tn ) + O (h 3 )
2
d h2
= y (tn ) + f (tn , y (tn ))h + [f (t , y (t ))]|t =tn + O (h 3 ),
dt 2
où
d
[f (t , y (t ))] = ∂t f (t , y (t )) + ∂y f (t , y (t ))y 0 (t )
dt
= ∂t f (t , y (t )) + ∂y f (t , y (t ))f (t , y (t ))
On en déduit que
h2
y (tn+1 ) = y (tn ) + f (tn , y (tn ))h + ∂t f (tn , y (tn ))
2
h2
+ ∂y f (tn , y (tn ))f (tn , y (tn )) + O (h 3 )
2
 21/33

Méthode de Taylor

Algorithme
Donnée : y0 ∈ Rn ,
For k = 0 : N − 1
h2
yn+1 = yn + f (tn , yn ))h + 2
[∂t f (tn , yn ) + ∂y f (tn , yn )f (tn , yn )]
End

Consistance : Ordre 2 : |τn+1 | ≤ Ch 2

Extension dimension supérieure :
Utiliser les dérivées partielles d’ordre supérieur

∂2tt f , ∂2ty f , et ∂2yy f ...

Difficultés : Pas de connaissance systématique de ces dérivées

partielles !
 22/33

Exemple: Le pendule
θ l

ml θ00 = −α`θ0 − mgsin(θ)


θ(0) = θ0 , et θ0 (0) = θ0 ,


0
mg

θ
!
Avec y = , cette équation se réécrit sous la forme
θ0
!
y2
y (t ) = f (t , y ) =
0
,
− mα y2 − g` sin(y1 )

Exercice
Écrire l’algorithme de Taylor d’ordre 2 dans le cas particulier de l’équation
du pendule.
 23/33

Méthode de Runge Kutta d’ordre 2

Idée : On souhaite déterminer des coefficients (a1 , a2 , a3 , a4 ) tels que

y (tn+1 ) = y (tn ) + a1 h f (tn , y (tn ))

+a2 h f (tn + a3 h , y (tn ) + a4 h ) + O (h 3 )

En remarquant que

f (tn + a3 h , y (tn ) + a4 h ) = f (tn , y (tn )) + a3 h ∂t f (tn , y (tn )) + a4 h ∂y f (tn , y (tn ))

On en déduit que

A = y (tn ) + a1 h f (tn , y (tn )) + a2 h f (tn + a3 h , y (tn ) + a4 h )

= y (tn ) + (a1 + a2 )h f (tn , y (tn )) + a2 a3 h 2 ∂t f (tn , y (tn ))
+a2 a4 h 2 ∂y f (tn , y (tn )) + O (h 3 )
 24/33

On rappelle que le développement de Taylor montre

h2
y (tn+1 ) = y (tn ) + f (tn , y (tn ))h + ∂t f (tn , y (tn ))
2
h2
+ ∂y f (tn , y (tn ))f (tn , y (tn )) + O (h 3 ).
2

Les coefficients (a1 , a2 , a3 , a4 ) doivent donc satisfaire les conditions

suivantes : 
a1 + a2 = 1



= 1/2



 a2 a3

= f (tn , y (tn ))/2

a2 a4

 25/33

Méthode d’Euler modifiée (RK2)

Choix des coefficients (a1 , a2 , a3 , a4 ) :
a1 = a2 = 12 ,
a3 = 1, et a4 = f (tn , y (tn ))

Algorithme (Runge Kutta 2)

Donnée : y0 ∈ Rn ,
for k = 0 : N − 1
k1 = yn + hf (tn , yn )
yn+1 = yn + h2 (f (tn , yn ) + f (tn + h , k1 ))
End
Lien avec la méthode des trapèzes
Z tn +1
y (tn+1 ) − y (tn ) = f (t , y (t ))dt
tn
h
= [f (tn , y (tn )) + f (tn+1 , y (tn+1 ))] + O (h 2 )
2
h
= [f (tn , y (tn )) + f (tn+1 , yn + hf (tn , y (tn ))))] + O (h 2 )
2
 26/33

Méthode du point du milieu

Choix des coefficients (a1 , a2 , a3 , a4 ) :
a1 = 0, a2 = 1, a3 = 1/2, et a4 = f (tn , y (tn ))/2

Algorithme (Méthode du point du milieu)

Donnée : y0 ∈ Rn ,
for k = 0 : N − 1
k1 = yn + h /2 f (tn , yn )
yn+1 = yn + hf (tn + h /2, k1 )
End

Lien avec la méthode des rectangles

Z tn +1
y (tn+1 ) − y (tn ) = f (t , y (t ))dt = hf (tn + h /2, y (tn + h /2)) + O (h 2 )
tn
!
h h
= h f tn + , yn + f (tn , y (tn ))) + O (h 2 )
2 2
 27/33

Le pendule sans frottement : simulation numérique avec

la méthode RK2
 28/33

Le pendule sans frottement : Exemple simulation

numérique avec la méthode RK2
 29/33

Méthode Runge Kutta d’ordre 4 (RK4)

Même idée avec des développements de Taylor d’ordre supérieur :

Algorithme (Méthode RK4)

Donnée : y0 ∈ Rn ,
for k = 0 : N − 1
k1 = f (tn , yn )
k2 = f (tn + h /2, yn + h /2k1 )
k3 = f (tn + h /2, yn + h /2k2 )
k3 = f (tn + h , yn + hk3 )
yn+1 = yn + h6 (k1 + 2k2 + 2k3 + k4 )
End

Méthode d’ordre 4, très utilisée en pratique ...

Voir le lien avec l’intégration de Simpson !
 30/33

Méthodes d’ordre supérieur multi-pas : formule

D’Adams-Bashforth
Idée : Utiliser les estimations de yn , yn−1 , yn−2 ...yn−k pour obtenir une
approximation de yn+1 .

Exemple avec k=1

Interpolation linéaire de y 0 (t ) = f (t , y (t )) au noeuds tn et tn−1 . Avec
fn = f (tn , yn ), alors
t − tn t − tn−1
y 0 (t ) = f (t , y (t )) ' fn−1 + fn ,
h h
et finalement,
Z tn +1 Z tn +1
y (tn−1 ) − y (tn ) = y 0 (t )dt = f (t , y (t ))dt
tn tn
Z tn +1
t − tn t − tn−1
' fn−1 + fn dt
tn h h
 31/33

Méthode d’ordre supérieur multi-pas : formule

D’Adams-Bashforth

Exercice
Montrer que
h
y (tn−1 ) − y (tn ) ' [3fn − fn−1 ]
2

Algorithme (Méthode d’Adams-Bashforth d’ordre 2)

Donnée : y0 ∈ Rn ,
for k = 0 : N − 1
fn = f (tn , yn )
yn+1 = yn + h2 (3fn − fn−1 )
End

Extension ordre supérieur // explicite-implicite !

 32/33

Utilisation d’un pas adaptatif

Motivation : L’utilisation d’un pas de temps h uniforme est très

coûteux dans certain cas, et il devient plus efficace d’utiliser un pas de
temps hn adaptatif, choisi plus petit dans les zones où y évolue fortement !

Rappel :
Schéma numérique : yn+1 = yn + hn φ(tn , yn )
p p +1
Consistance : τn (hn ) ' ωn hn + O (hn )
P
Consistance : maxn |en | ≤ C (T ) n hn |τn (hn )|
 33/33

Estimation de ωn à l’instant tn
Avec
ynh+1 = yn + h φ(tn , yn )
/2 /2 h /2 /2
ynh+ 1
= ynh+1/2
+ φ(tn , ynh+1/2
), où ynh+ 1/2
= yn + h φ(tn , yn )
2
alors
/2
ynh+ 1
− ynh+1
ωn ' .
h p +1 (1 − 2−p )
Déterminer hn afin de satisfaire maxn |en | ≤ 
Il suffit d’imposer
C (T )|τn | ≤ /T .
En effet X X 
max |en | ≤ C (T ) hn |τn (hn )| ≤ hn ≤
n
n
T
Alors !1/p
 
|τn | ≤ ⇐⇒ hn .
TC (T ) TC (T )|ωn |
 1/25

Chapitre 7 : Méthodes numériques pour les équations

aux dérivées partielles (EDP)

...
 2/25
Ce chapitre s’intéresse aux méthodes des différences finies pour la
résolution d’équations aux dérivées partielles. On regardera en particulier
le cas des équations suivantes :

Equation de Poisson :

−∆u(x ) = f (x )

Equation de la chaleur :

∂t u(x , t ) = ∆u

Equation de transport :

∂t u(x , t ) + v .∇u = 0

Equation des ondes :

∂2tt u(x , t ) = c 2 ∆u
 3/25

Méthodes des différences finies

Approximation des dérivées d’ordre 1 :
Schéma centré :
u(x + h ) − u(x − h )
u0 (x ) = + O (h 2 )
2h
Schéma décentré à gauche :
u(x ) − u(x − h )
u0 (x ) = + O (h )
h
Schéma décentré à droite :
u(x + h ) − u(x )
u0 (x ) = + O (h )
h
Approximation des dérivées d’ordre 2 :
u(x + h ) − 2u(x ) + u(x + h )
u00 (x ) = 2
+ O (h 2 )
h

Exercice
Montrer les approximations précédentes en explicitant les constantes
dans les O (h ) et O (h 2 )
 4/25

1 Exemple d’un problème aux limites : l’équation de Poisson

2 Exemple d’un problème d’évolution : l’équation de la chaleur

3 Exemple d’un problème d’évolution : l’équation de transport

4 Exemple d’un problème d’évolution : l’équation des ondes

 5/25

Equation de Poisson
On s’intéresse à la résolution numérique de l’équation

−u00 (x ) = f (x ) pour x ∈ [0, 1]


(P1 )
u(0) = g0 , u(1) = g1



Discrétisation de l’intervalle [0, 1] avec N + 2 points :

On pose xi = i h pour i ∈ [0, N + 1] et où h = N + 1
1.
Approximation numérique de u et f aux noeuds xi

ui ' u(xi ),

 u0 = g0 et uN +1 = g1

fi = f (xi )


Approximation de l’EDP :

1
− (ui −1 − 2ui + ui +1 ) = fi , ∀i = 1 : N
h2
 6/25

Approximation de l’équation de Poisson

Avec

 2 −1 0 · · · 0 
 
 f1 + g0 /h 2
   
 u1 .. 
 −1 . . . .. ..
 
. .

 u2 f2 . 
   
 . ..
  
1 
uh =  ..  , bh =   , Ah = 2  0 . . . .. .. 
. . 0 
  
. h 
 .. .. .. ..
    
 uN −1 fN −1
 . . . . −1 
   
fN + g1/h 2
   
uN
0 ...

0 −1 2

La solution uh s’obtient en résolvant le système linéaire

Ah uh = bh .

Propriété
La matrice Ah est symétrique définie positive.
z11 +(z1 −z2 )2 +···+(zN −1 −zN )2 +zN2
Preuve Montrer que < Z , Ah Z >= h2
 7/25

Propriétés de la matrice Ah
Exercice
Montrer que les vecteurs {v n }n=1:N dont les composantes sont définies par
vin = sin(nπxi ), sont des vecteurs propres de Ah associés respectivement
aux valeurs propres
2
λnh = 2 (1 − cos (nπh )) .
h
En déduire que pour N suffisamment grand, les valeurs propres λn de Ah
sont comprises entre
4
λ1h ' π2 et λN
h ' 2,
h
et que
4 1
kAh k2 ' 2 et kAh−1 k2 ' 2 .
h π
On admet aussi par la suite que
1
kAh−1 k∞ ≤ .
8
 8/25

Consistance de la méthode des différences finies

Soit u la solution du problème (P1 ) et uh le vecteur défini par (uh )i = u(xi ).
Définition
Le schéma numérique est dit consistant si

Πh (u) = Ah uh − bh → 0 quand h → 0.

En particulier, le schéma est d’ordre p pour une norme k.k donnée si

kΠh (u)k ≤ Ch p .

Propriété
Supposons que la solution u de (P1 ) soit de régularité C 4 ([0, 1]). Alors le
schéma est consistant d’ordre 2 pour la norme k.k∞ , i.e.

h2 n o
kΠh (u)k∞ ≤ sup |u(4) (x )|
12 x ∈[0,1]
 9/25

Convergence de la méthode des différences finies

Démonstration : utiliser un développement de Taylor de u :

u(x + h ) − 2u(x ) + u(x − h ) h 2 (4)

u00 (x ) = +2 u (ηx ), avec ηx ∈ (x −h , x +h ).
h2 4!

Théorème
Supposons que la solution u de (P1 ) soit de régularité C 4 ([0, 1]). Alors, le
schéma est convergent d’ordre 2 pour la norme k.k∞ , i.e.

h2 n o
kuh − uh k∞ ≤ sup |u(4) (x )|
96 x ∈[0,1]

Démonstration :
1
Utiliser la propriété précédente et kAh−1 k∞ ≤ 8
 10/25

1 Exemple d’un problème aux limites : l’équation de Poisson

2 Exemple d’un problème d’évolution : l’équation de la chaleur

3 Exemple d’un problème d’évolution : l’équation de transport

4 Exemple d’un problème d’évolution : l’équation des ondes

 11/25

Equation de la chaleur

On s’intéresse à la résolution numérique de l’équation de la chaleur





 ∂t u(x , t ) = ∆u(x , t ) + f (x , t ), pour (x , t ) ∈ [0, T ] × [0, 1]
u(x , 0) = u0 , pour x ∈ [0, 1]

(P2 )




u(0, t ) = u(1, t ) = 0, pour t ∈ [0, T ]



Discrétisation de l’intervalle [0, 1] avec N + 2 points :

1
On pose xi = i h pour i ∈ [0 : N + 1] et où h = N + 1.
Discrétisation de l’intervalle [0, T ] avec M + 1 points :
On pose tn = n δt pour n ∈ [0 : M ] et où δt = M1 .
Approximation numérique de u et f aux noeuds (xi , tn )

uin ' u(xi , tn ), fin = f (xi , tn ), u0n = uM

n
+1 = 0, et ui = u0 (xi )
0
 12/25

Méthode d’Euler explicite

Approximation de l’EDP :

uin+1 − uin uin+1 − 2uin + uin−1

= + fin ,
δt h2
∀(i , n) ∈ [1 : N ] × [0, M − 1].
Schéma numérique :

uhn+1 = (Id − δt Ah )uhn + δt fhn ,

où
uhn = (u1n , u2n , · · · , uN ) , et fhn = (f1n , f2n , · · · , fNn )t ,
n t
 13/25

Méthode d’Euler explicite

n n
On note uh le vecteur défini par (uh )i = u(xi , tn ) où u est la solution du
problème (P2 ).

Propriété
Supposons que la solution u de (P2 ) soit de régularité
C 2 ([0, T ], C 4 ([0, 1])). Alors le schéma est consistant d’ordre 2 en espace
et d’ordre 1 en temps pour la norme k.k∞ , i.e.
 
kΠnh (u)k∞ ≤ C δt + h 2 ,

où
Πnh (u) = (uhn+1 − uhn+1 )/δt + δt Ah unh − fhn

Démonstration : Utiliser un développement de Taylor.

 14/25

Propriété
Le schéma d’Euler explicite est L 2 -stable sous la condition

1 2
δt ≤ h .
2

Démonstration : Regarder les valeurs propres de la matrice (Id − δt Ah ).

 15/25

Méthode d’Euler implicite

Approximation de l’EDP :

uin+1 − uin uin++11 − 2uin+1 + uin−+11

= + fin+1
δt h2
∀(i , n) ∈ [1 : N ] × [0, M − 1].
Schéma numérique :
Le vecteur uhn+1 est donc solution du système linéaire

(Id + δt Ah )uhn+1 = uhn + δt fhn+1

Exercice
On suppose que la solution u de (P2 ) est suffisamment régulière. Montrer
que le schéma d’Euler implicite est consistant d’ordre 2 en espace et
d’ordre 1 en temps et qu’il est inconditionnellement L 2 -stable
 16/25

Schéma de Crank-Nicolson
Approximation de l’EDP :

 n+1
uin+1 − uin  ui +1 − 2uin+1 + uin−+11


=  + f n +1  /2
δt h2 i 
 n
 ui +1 − 2uin + uin−1

+ fi  /2
n

+ 
h2
Schéma numérique :
Le vecteur uhn+1 est donc solution du système linéaire
 
(Id + δt Ah /2)uhn+1 = (Id − δt Ah /2)uhn + δt fhn + fhn+1 /2

Exercice
On suppose que la solution u de (P2 ) est suffisamment régulière. Montrer
que le schéma d’Euler implicite est consistant d’ordre 2 en espace et en
temps et qu’il est inconditionnellement L 2 -stable
 17/25

1 Exemple d’un problème aux limites : l’équation de Poisson

2 Exemple d’un problème d’évolution : l’équation de la chaleur

3 Exemple d’un problème d’évolution : l’équation de transport

4 Exemple d’un problème d’évolution : l’équation des ondes

 18/25

Equation de transport

On s’intéresse à la résolution numérique de l’équation de transport


∂t u(x , t ) + c ∂x u(x , t ) = 0


(P3 )
u(x , 0) = u0



qui a pour solution

u(x , t ) = u0 (x − ct ).
Cette équation vérifie en particulier des principes de stabilité qu’on
souhaiterait conserver numériquement !
Stabilité L ∞ :
c0 ≤ u0 ≤ c1 =⇒ c0 ≤ u(x , t ) ≤ c1
Stabilité L 2
ku(x , t )k2 ≤ ku0 k2
 19/25

Schéma explicite différences finies centré

Schéma numérique :

uin+1 − uin uin+1 − uin−1

+c = 0 avec u0i = u0 (xi )
δt 2h

Ce schéma est inconditionnellement instable et s’avère inutilisable. Il

ne vérifie aucun des principes de stabilité précédents !

Exercice
Avec u0 (x ) = e ipx , calculer la solution exacte de (P3) et la solution
numérique uhn . Remarques ?
 20/25

Schéma explicite décentré en amont

Schéma numérique ( cas c > 0) :
uin+1 − uin uin − uin−1
+c = 0 avec u0i = u0 (xi )
δt h

Propriété
Ce schéma est L ∞ -stable sous la condition de Courant-Friedrichs-Levy
(CFL)
c δt ≤ h .

Démonstration : Montrer que sous la CFL, uin+1 est une combinaison

n o
convexe des uin .
i ∈Z

Exercice
Si u0 ∈ C 2 (R), montrer que le schéma numérique est consistant à l’ordre 1
en temps et en espace!
 21/25

1 Exemple d’un problème aux limites : l’équation de Poisson

2 Exemple d’un problème d’évolution : l’équation de la chaleur

3 Exemple d’un problème d’évolution : l’équation de transport

4 Exemple d’un problème d’évolution : l’équation des ondes

 22/25

Equation des ondes

On s’intéresse à la résolution numérique de l’équation des ondes





 ∂tt u(x , t ) = c 2 ∆u(x , t )
u(x , 0) = u0 (x )

(P4 )




∂t u(x , 0) = u0 (x )


0

Formule de d’Alembert
Z x +ct
1 1
u(x , t ) = [u0 (x − ct ) + u0 (x + ct )] + u00 (y )dy
2 2c x −ct

Conservation de l’énergie
Z 
1 
E(t ) = (∂t u(x , t ))2 + c 2 (∂x u(x , t ))2 dxdt
2
 23/25

Approche directe explicite

Approximation de l’EDP :

uin+1 − 2uin + uin−1 un − 2uin + uin−1


2 i +1 ui0 = u0 (xi )

,

=c et
δ2t h2 u1 = u0 (xi ) + δt u0 (xi )


i 0

Schéma numérique :
 
uhn+1 = 2Id − c 2 δ2t Ah uhn − uhn−1 .

Propriété
On suppose que la solution u de (P4 ) est suffisamment régulière. Alors, le
schéma est consistant à l’ordre 2 en temps et en espace. Ce schéma est
de plus stable sous la condition CFL

h
c δt ≤ .
2
 24/25

Approche Euler implicite

Idée : Réécrire le problème (P4 ) sous la forme d’un système d’EDP

d’ordre 1 :

u (x , t )
! !
0 Id
∂t Y (x , t ) = Y (x , t ), avec Y (t ) =
c 2∆ 0 ∂t u(x , t )

Euler implicite : Yhn+1 solution de

!
Id −δt Id
Bh Yhn+1 = Yhn avec Bh =
+δt c 2 Ah Id

Propriété
Ce schéma est d’ordre 1 en temps et d’ordre 2 en espace . Il est de plus
inconditionnellement stable
 25/25

Schéma de Crank-Nicolson

Exercice
Montrer que l’approche Crank-Nicolson conduit au schéma consistant
d’ordre 2,
Ch Yhn+1 = Dh Yhn
avec

−δt Id /2 +δt Id /2
! !
Id Id
Ch = et Dh =
+δt c 2 Ah /2 Id −δt c 2 Ah /2 Id