Analyse Numerique
Analyse Numerique
Analyse Numerique
...
2/18
1 Définition
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 ?
Quelques exemples
Quelques exemples
Temps de calcul pour l’inversion d’une matrice
7/18
x = ±m b ±e , avec b ≥ 2.
avec la mantisse m :
et l’exposant e :
Représentation unique ?
Ordinateur : mémoire finie !
La mantisse est tronquée au bout de r chiffres,
Valeur maximale pour s
9/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
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
∂f xi
ki = .
∂xi f (x )
Exercice
Calculer le conditionnement de l’addition et de la multiplication
16/18
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.
Exercice
Proposer un autre algorithme plus stable
17/18
Complexité algorithmique
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
Coût algorithmique ?
1/41
...
2/41
A partir d’un ensemble de points (xi , yi )
i =0:n , on recherche dans un
espace de fonctions,
1 Interpolation polynomiale
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
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
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
Propriété
Le déterminant de la matrice B vérifie
Y
det (B ) = (xj − xi ).
0≤i <j ≤n
9/41
Exercice
Montrer que la famille {L0 , L1 , L2 } est libre et expliciter la décomposition
de p2 dans cette base
11/41
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
Erreur d’interpolation
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
(f (x ) − pn (x )) ni=0 (t − xi )
Q
φ(t ) = f (t ) − P (t ) − Qn
i =0 (x − xi )
15/41
Erreur d’interpolation
où h = (b − a )/n, et
!n+1
b −a 1
|E | ≤ sup f (n+1) (x ) .
x ∈[a ,b ] n 4(n + 1)
Exercice
Montrer le résultat précédent
16/41
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
et limn→∞ E = 0.
17/41
1 Interpolation polynomiale
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
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
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
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
hi = |xi − xi −1 | = h = (b − a )/n
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
− 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
( x − xi ) ( x − x i −1 )
pi00 (x ) = yi00−1 + yi00 ,
(xi −1 − xi ) (xi − xi −1 )
(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
(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
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
− 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
1 Interpolation polynomiale
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
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
Pn 2 Pn 2
où D = (n + 1) i =0 xi −( i =0 xi ) .
38/41
φ(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
Exercice
Calculer le gradient des fonctionnelles suivantes
φ1 (x ) =< x , b >
φ2 (x ) =< Ax , b >
φ3 (x ) = 1 < Ax , Ax >
2
40/41
avec
Le gradient de φ s’identifie à
∇φ(a ) = 2 B t Ba − B t y .
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
...
2/17
a x1 x2 x3 x4 a x1 x2 x3 x4 a x1 x2 x3 x4
b b b
1 Formules de Newton-Cotes
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
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
1 Formules de Newton-Cotes
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
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
...
2/49
Ax = b .
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
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
1 Méthodes directes
Décomposition LU
Décomposition Cholesky
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
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
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
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
En effet,
1 0 0 18 18
2 1 0 −12 = 24
3 2 1 2 32
12/49
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
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
Definition
Une matrice A est définie positive si
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
2,1 `i ,1
`i ,2 = a2,i −` , pour i > 2
`2,2
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
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
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
(A + δA ) (x + δx ) = b + δb ,
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
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
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
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
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
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
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
Ax = b ,
φ(x ∗ ) = x ∗ .
lim xk = x ∗
k →∞
36/49
Exemple en dimension 1
φ (x)
f(x) f(x)
φ (x)
x* x 2 x1 x0 x* x 0 x 1 x 2 x3
|a | < 1
37/49
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
Théorème
Soit φ : Rn → Rn une contraction et x0 ∈ Rn . Alors la suite xk définie
x0 ∈ Rn
xk +1 = φ(xk ) ∀k > 0
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
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
Exercice
Comment choisir c pour que l’algorithme s’arrête avec une précision de
sur x ∗ ?
41/49
φ(x ) = M −1 (Nx + b ) .
où A = M − N.
Exercice
Montrer que si x ∗ est un point fixe de φ alors
Ax ∗ = b
k|M −1 N k| < 1.
Mxk +1 = Nxk + b .
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
Méthode de Gauss-Seidel
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
1 Méthodes directes
Décomposition LU
Décomposition Cholesky
Méthode du gradient
Ax = b
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,
Choix de αk :
On choisit αk comme le pas qui minimise J dans la direction dk :
Avec
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
...
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
1 Méthode de la Dichotomie
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
Exercice
Avec Ik = bk − ak , montrer que
b −a
Ik = .
2k
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
1 Méthode de la Dichotomie
3 Méthode de Newton
9/27
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 ),
φ (x)
f(x) f(x)
φ (x)
x* x 2 x1 x0 x* x 0 x 1 x 2 x3
11/27
|g 0 (x )| ≤ λ < 1, ∀x ∈ [a , b ].
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
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
Démonstration :
Appliquer le théorème global sur l’intervalle Iα = [x ∗ − α, x ∗ + α] avec α
suffisamment petit.
14/27
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
kg (x ) − g (y )k ≤ λkx − y k, ∀(x , y ) ∈ B (x , R ).
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 Méthode de la Dichotomie
3 Méthode de Newton
19/27
f (x0 )
=⇒ x1 = x0 −
f 0 (x0 )
x*
x2 x1 x0
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 )
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
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 )
où ∂ f1 ∂ f1
...
∂x ∂ xn
1
Df (x ) = ..
. ..
.
∂ fn ∂fn
...
∂x1 ∂ xn
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
g (α, β, x ) = αx β
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
...
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 ,
Exemple: Le pendule
θ l
ml θ00 = −α`θ0 − mgsin(θ)
θ(0) = θ0 , et θ0 (0) = θ0 ,
0
mg
θ θ0
! ! !
y2
y= , f (t , y ) = , et y (0) = ,
θ0 − mα y2 − g` sin(y1 ) θ00
3/33
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.
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
1
y (t ) = , pour t < t0 + y0−1 .
y0−1 + t0 − t
5/33
Discrétisation de l’équation
yn ' y (tn ).
yn+1 = yn + h φ(tn , yn ).
6/33
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
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
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
h2
|en+1 | ≤ (1 + Lh )|en | + M.
2
11/33
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
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
yn ' y (tn ).
yn+1 = yn + h φ(tn , yn ).
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
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
En remarquant que
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
Exercice
Montrer que
h
y (tn−1 ) − y (tn ) ' [3fn − fn−1 ]
2
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
...
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
∂2tt u(x , t ) = c 2 ∆u
3/25
Exercice
Montrer les approximations précédentes en explicitant les constantes
dans les O (h ) et O (h 2 )
4/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
Approximation de l’EDP :
1
− (ui −1 − 2ui + ui +1 ) = fi , ∀i = 1 : N
h2
6/25
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
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
Πh (u) = Ah uh − bh → 0 quand h → 0.
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
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
Equation de la chaleur
Approximation de l’EDP :
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
Propriété
Le schéma d’Euler explicite est L 2 -stable sous la condition
1 2
δt ≤ h .
2
Approximation de l’EDP :
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
Equation de transport
Schéma numérique :
Exercice
Avec u0 (x ) = e ipx , calculer la solution exacte de (P3) et la solution
numérique uhn . Remarques ?
20/25
Propriété
Ce schéma est L ∞ -stable sous la condition de Courant-Friedrichs-Levy
(CFL)
c δt ≤ h .
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
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
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
u (x , t )
! !
0 Id
∂t Y (x , t ) = Y (x , t ), avec Y (t ) =
c 2∆ 0 ∂t u(x , t )
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