ANALYSE NUMERIQUE

(Cours et Td)

Licence 3

M. AKPEMADO Komi

1
 Chapitre 1 Résolution de systèmes linéaires

Généralités

Introduction

Résoudre un système linéaire c’est de déterminer la ou les solutions de l’équation

algébrique Ax=b (1) A est une matrice carrée d’ordre n à valeurs dans C et b un
vecteur appelé second membre. On dispose des célèbres formules de Cramer.
Malheureusement ces formules sont totalement inadaptées pour calculer
efficacement la solution des systèmes linéaires. En effet leur cout en temps de calcul
sur un ordinateur est prohibitif. Il faut calculer (n+1) déterminants, et si on utilise la
méthode de développement par ligne ou par colonne, chaque déterminant
demande plus de n ! multiplications. Au total cette méthode de Cramer nécessite
(n+1) ! multiplications. Cela est inenvisageable. Pour n=50, si les calculs s’effectuent
sur un ordinateur de 109 opérations par seconde, le temps de calcul est de l’ordre de
4.8 x 1049 années. Même si on utilise une meilleure méthode pour calculer les
déterminants, la méthode de Cramer n’est pas compétitive par rapport aux méthodes
que nous allons étudier.

NB : Actu l’ordi le plus puissant au monde est le Supercalculateur Frontier. Un milliard

de milliard d’opérations par seconde.

Normes et rayon spectral:

Définition norme vectorielle :

x  x , R n → R+ vérifiant :

- x  0 et x = 0  x = 0
-  .x =  . x
- x+ y  x + y

Exemple de norme vectorielle :

- Norme 1 : X 1
=  xi
i

= x
2
- Norme 2 : X 2 i = X, X
i

- Norme  : X 
= Sup xi
i

Proposition : En dimension finie, toutes ces normes sont équivalentes ; en pratique,

on choisit celle qui nous arrange.

2
Définition norme matricielle :

C'est une norme dans l'espace vectorielle Mn(R), avec en plus A  B  A . B .

Norme matricielle induite :

A partir d'une norme vectorielle, on construit une norme matricielle induite :

 AX 
A = Sup =
 Sup ( AX ) .
X 0
 X  X =1 ou X 1

Notons de plus que I = 1 .

En particulier pour la norme 2, on a la relation : I 2

=  (A t A) avec  le rayon
spectral.

Propriétés :

- AX  A . X et AB  A . B
- Cas de A symétrique : A 2
=  ( A) = max i
i

- Dans le cas général, A

2
=  max , maximum des valeurs singulières  i = i avec
i les valeurs propres (positives) de la matrice symétrique A t A . (?)
- Soit A une matrice carré quelconque n x n. Pour toutes normes de matrice induite,
on a  ( A)  A .

Conditionnement d’une matrice inversible :

Considérons un système linéaire AX=b. L'expérience de Wilson met en évidence une
instabilité des calculs dans certains cas. Si l'on modifie sensiblement les paramètres de
la matrice A, ou du second membre b, on peut trouver des résultats complètement
différends !

- Cond ( A) = A . A −1
- Cond(A) = Cond( A)
- Cond ( A)  1 avec Cond(Id ) = 1
- Soit une matrice Q orthogonale : QX 2 = X 2 et Cond(Q) = 1
- Les bons conditionnements sont les petits conditionnements (au mieux 1 pour les
matrices orthogonales).
 
- Pour A symétrique : Cond2 ( A) = max ; dans la cas général : Cond2 ( A) = max
 min  min

3
Théorèmes sur le conditionnement :

X b
Théorème : AX = B et A( X + X ) = b + b . On a  Cond ( A) .
X b
X A
Théorème : AX = B et ( A + A)( X + X ) = b . On a  Cond ( A) .
X + X A

Remarque : Le conditionnement traduit l'amplification des grandeurs relatives.

Inverse numérique d’une matrice :

Soit A une matrice inversible et A-1 son inverse théorique.

M est un inverse de A du point de vue numérique si :

M − A−1
- −1
est petit, ce qui correspond à M=A-1
A
M −1 − A
- est petit, ce qui correspond à A=M-1
A
- AM − Id et MA − Id sont petits, ce qui correspond à
AM − Id = 0 et MA - Id = 0
On s'adapte au calcul que l'on veut faire. Notons que l'on résout rarement un système
linéaire en calculant l'inverse X = A−1b

A) Systèmes linéaires Ax=b : méthodes directes

0. Méthode de Cramer
Cette méthode nécessite le calcul d'un déterminant, dont la complexité est en n!. La
compléxité de la formule de Cramer est en n 2 n!. Par conséquent, cette méthode n'est
pas utilisable en informatique.

4
 1. Méthode du pivot de Gauss
 a  a
À la kème étape : pour i  k , Ligne i −  i ,k  Ligne k . On pose li ,k = i ,k .
a  a k ,k
 k ,k 

n3
L'algorithme à une complexité en .
3

Le problème des coefficients diagonaux nuls : il faut permuter les lignes lorsque
apparaît un zéro sur la diagonale.

Stratégie du pivot maximum partiel :

On permute la kème ligne et la ième ligne avec i tel que ai ,k = Sup al ,k . Pour des raisons
l k

de stabilité numérique, on divise par le terme de la colonne k, qui a la plus grande

valeur absolue.

Stratégie du pivot avec seuil :

On effectue une permutation uniquement des a k ,k   . Cependant cette méthode
peut être mise en défaut. Il vaut mieux utiliser la méthode du pivot maximum partiel,
qui n'est pas vraiment plus coûteuse en calcul.

Calcul de l’inverse
Une idée pour résoudre le sytème Ax=b serait de calculer l’inverse de A. En effet, on a
x = A −1b .

Méthode de Shermann – Morison

Soit A une matrice n x n, réelle inversible, dont on connaît l’inverse A-1 ; soit B une
perturbation, on cherche à calculer ( A + B ) .
−1

( )
On démontre que si B = XY t , alors ( A + B)−1 = I + c.A−1 B A−1 = A−1 + c.A−1 BA−1 , avec
−1
c= .
1 + Y t A −1 X

On suppose maintenant que la perturbation ne modifie qu’un seul en (i, j ) :

B =  .Ei , j =  . ei  e j .
t

 
X Y

 k , i j ,l

Notons A−1 = ( k , j )1k ,l n . Alors on trouve (A +  .Ei , j )
−1

k ,l
=  k ,l − 
1 +  . j ,i
.

Elimination de Gauss
On applique la stratégie du pivot partiel maximum, au calcul de l’inverse.

On considère une matrice A inversible, de dimension n

5
On résout Ax = el pour 1  l  n , ce qui revient à écrire A. A −1 = I . On procède par
élimination de Gauss, ce qui donne une matrice triangulaire supérieure (complexité
en n3). Ensuite, on réalise n remontées, une par vecteur el (complexité en n 2).

En résumé, on a Ax = el et U .x = M .el avec U = MA une matrice triangulaire

supérieure.

Remarques :

- Choix du pivot : pour des raisons de stabilité numérique, on divise toujours par le
terme de la colonne k qui a la plus grande valeur absolue (Cf. Pivot partiel
maximum).
- Recherche du pivot : si la plus grande valeur absolue est inférieure à la précision
machine , alors on arrête en disant que la matrice n’est pas inversible à  près
(test d’inversibilité informatique). Cependant, elle peut très bien être inversible du
point de vue mathématique !

Factorisation A=LU
Soit A une matrice régulière, de dimension n.

Lien avec la méthode de Gauss

On applique la méthode de Gauss sans pivoter. U est une matrice triangulaire
supérieure réelle, telle que  J k A = U .
k =1, n

1 
 
 1 
On donne la définition des matrices J k =   k +1,k 1  avec  = − al ,k pour
  l ,k
a k ,k
   
  n ,k 1 

lk.

J
−1
On démontre que l’inverse des Jk existe et que L = k est une matrice
k =1, n

triangulaire inférieure, avec des 1 sur la diagonale. En fait, on a

1 
 
 1 
−1
Jk =  −  k +1,k 1 .
 
   
 −  n ,k 1 


Calcul algorithmique des coefficients de L et de U

On procède par identification en utilisant les propriétés de L (triangulaire inférieure
avec une diagonale de 1) et de U (triangulaire supérieure). On cherche un
algorithme par ligne donnant les li, j pour 1  j  i − 1 et les ui , j pour 1  j  n .

6
 Pour i := 1 à n faire
Début
Pour j := 1 à i-1 faire
 j −1

li , j :=  ai , j −  li ,k .u k , j  a j , j ;
 k =1 
/* avec l i ,i = 1 */

Pour j := i à n faire
i −1
u i , j := a i , j −  l i ,k .u k , j ;
k =1
Fin ;

Remarque : L et U écrasent la matrice A.

Cas particulier des matrices à « structure bande »

  
 
   
   
 
    
   

   
 
  

- La largeur de la bande est W = P+Q+1

- On ne stocke dans une structure de donnée appropriée que la bande, soit un
O(W.N), ce qui est intéressant si W << N.

Propriété fondamentale :
La bande se conserve par factorisation LU. On parle de structure de données statique.
En revanche, si on utilise des techniques de pivotage au cours d’une factorisation LU
sur une matrice bande, alors la structure bande n’est pas conservée.
On utilise la méthode uniquement dans le cas, où à priori on ne pivote pas. La
condition suffisante la plus connue est « A matrice à diagonale strictement
dominante », c’est à dire i  1, n ai ,i  a
j =1, n
i, j .
j i

7
Cas des matrices symétriques définies positives (SDP)
Cf. factorisation Cholesky et Crout.

Factorisation PA=LU

P est la matrice des permutations.

L est une matrice triangulaire inférieure dont la diagonale se compose de 1.
U est une matrice triangulaire supérieure.

1 
 
 
Matrice de permutation élémentaire i  j : 


 1 
 0 1  → ligne i
 
 1 
Pi , j =   
 
 1 
 1 0  → ligne j
 
 1 
 
  
 1 


Le produit PA permute les lignes i et j de la matrice A. Le produit AP permute les

colonnes i et j de la matrice A.

8
 2) Factorisation A=BBt : Cholesky

Cette factorisation s'applique pour des matrices A symétriques et définies positives.

Théorème : Soit A une matrice SDP de dimension n. Il existe une et une seule matrice
B triangulaire inférieure, telle que A=BBt.

3) Factorisation A=LDLt : Crout

Théorème : Soit A une matrice SDP de dimension n. Il existe une et une seule matrice
L triangulaire inférieure et D diagonale positive, telles que A=LDLt.

Remarque : Crout a le même comportement numérique que Cholesky. Cependant, il

est plus utilisé car généralisable au cas complexe.

Soit A une matrice symétrique définie positive (quelconque ?).

L est une matrice triangulaire inférieure de dimension n dont la diagonale se

compose de 1.
D est une matrice diagonale de dimension n.

L est en fait la matrice B pour laquelle les termes d'une colonne sont divisés par les
coefficients bii. D est la matrice dont les coefficients diagonaux sont les bii².

Algorithme par ligne pour l’obtention de la factorisation

 1 00  d1 
   
On écrit L =   0  et D =   .
l   d n 
 i, j 1 
j −1
Il s'agit de calculer les li,j et les di. On a ai , j =  lik d k l jk + lij d j .
k =0

On établit un algorithme par ligne : Pour la 1ère ligne, on écrit d1 = a1,1 . Supposons
maintenant le calcul effectué jusqu'à la ligne i-1 ; le calcul de la ligne i est donné par
j −1

 (l d k l jk )
( )
ik j −1
li , j = ai , j − k =1
pour j  i − 1 et par d i = ai ,i −  lik d k . Ainsi une CNS pour
2

dj k =1

cette factorisation est que les di ne soient pas nuls.

On donne l’algorithme correspondant :

9
 Pour i de 1 à n faire
Début
Pour j de 1 à i-1 faire
j −1
li , j := ai , j −  (lik d k l jk ) / d j ;
k =1

/* l i ,i := 1 */

( )
j −1
d i := ai ,i −  lik d k ;
2

k =1
Fin ;

Remarque : L et D écrase la partie triangulaire inférieure de A.

On veut maintenant apporter une amélioration à cet algorithme en diminuant les

nombres de calcul ; on économise les produits l i , k .d k , que l’on stocke pour un i fixé
dans le vecteur C[k].

On donne l’algorithme modifié :

Pour i de 1 à n faire
Début
Pour j de 1 à i-1 faire
Début
j −1
C[ j] := ai , j −  l jk .C[k ] ;
k =1
C[ j ]
li , j := ;
dj
Fin ;
j −1
d i := ai ,i −  lik .C[k ] ;
k =1
Fin ;

10
Algorithme par colonne pour l’obtention de la factorisation
Principe : On construit L et D, en procédant par décomposition successives. On réalise
 v1
t

 1d 
 d 1
le découpage suivant : A = A1 =   , avec v1 un vecteur colonne de
 v1 B1 
 d1
 
dimension n-1, et B1 une matrice carré de dim n-1.

 
 1 O   d1 O
, A =   t
On pose L1 = 
v1 .v1
 2   , avec B1 = B1 − .
 v1 O 
B1 
d1
I n −1  
 d1 
Ainsi, on a : A = A1 = L1 A2 L1 . Puis, en réitérant n-1 fois, il vient :
t

1 
 
  
 
A = A1 = (L1  Ln −1 )An (L1  Ln −1 ) . Le produit des matrices Lk =  1
t
   
L D Lt   
 vk 
 1
 dk 
 1 
 
 1 
donnent la matrice finale L =   .
 
 v1 v2
 1
 d1 d2 

Ce qui correspond en effet à une factorisation de Crout, car la matrice D = (d i ) est

bien diagonale et la matrice L est bien triangulaire inférieure, par construction.

On donne l’algorithme de calcul de L et de D par colonne, qui procède en écrasant la

partie triangulaire inférieure de A, (les 1 de la diagonale de L ne sont pas stockés, car
on préfère y stocker D) :

Début
Pour k de 1 à n-1 faire
Pour j de k+1 à n faire
Début
d := a j ,k a k ,k ;

11
 Pour i de j à n faire ai , j := ai, j − ai ,k .d ;
a j , k := d ;
Fin ;
Fin ;

La structure donnée d’implémentation serait un vecteur à une dimension stockant

consécutivement les colonnes de la partie triangulaire inférieure de A.

Remarque fondamentale : ai , j est modifié par a i , k et a j , k .

Utilisation de la factorisation pour résoudre Ax=b

Il s’agit de résoudre Ax = b en utilisant la factorisation de Cholesky - Crout:
ALDLt X = b . On procède en trois étapes.

n2
- Descente : Posons y = DLt x et résolvons Ly = b , dont la complexité est en
2
(matrice triangulaire). On écrit un algorithme utilisant un accès par ligne dans la
matrice L.

y1 := b1 ;
i −1
Pour i de 2 à n faire yi = bi −  li , j y j ;
j =1

- Résolution du système diagonal : Posons z = Lt x et résolvons Dz = y (par n

divisions).

yi
Pour i de 1 à n faire z i = ;
di
n2
- Montée : Résolvons finalement Lt x = z , dont la complexité est en .
2
On écrit un algorithme utilisant un accès par ligne dans Lt, ce qui revient à un
accès par colonne dans L.
xn := yn ;
n
Pour i de n-1 à 1 faire xi := zi − l xj ;
t
i, j
j = i +1

On veut exploiter la symétrie. On ne stocke que L et D ; et on ne stocke pas les 1 de la

diagonale de L, qui sont pris en compte directement dans le calcul.

12
Considérons le fait que l’on stocke les données en mémoire de façon mono-
dimensionnelle. On va donc être obligé de choisir une organisation ligne par ligne,
ou colonne par colonne des données.
Si l’on choisit une organisation ligne par ligne, on s’aperçoit que l’algorithme de
descente est bien adapté à la structure de données par ligne, tandis que l’algorithme
de montée, qui pratique un accès aux données par colonne ne l’est pas. C’est
pourquoi, il faut écrire un algorithme de montée accédant aux données par ligne et
non plus par colonne.

Pour i de 1 à n faire xi := zi ;
Pour i de n à 2 par pas de -1 faire
Pour j de i-1 à 1 par pas de -1 faire
x j := x j − li , j xi ;

B) Systèmes linéaires Ax=b : méthodes itératives

Principe

Méthode itératives classiques :

On cherche à obtenir une suite de vecteur ( X n )n convergente vers X telle que AX = b .
En posant A = M − N où M est une matrice inversible, on propose la formule de
récurrence MX n +1 = NX n + b ou encore X n+1 = M −1 NX n + M −1b . Si X est solution du
problème, il vient que (X − X n ) = (M −1 N )  (X − X 0 ). B = (M −1 N ) est appelée la
n

matrice d'itération.
Ainsi pour que X n → X , il est nécessaire que (M −1 N ) → 0 .
n

Remarque : Le cadre normale de convergence des méthodes itératives classiques est :

« A matrice symétrique définie positive ».

Vitesse de convergence :
Théorème : (B n )n → 0 ssi  (B )  1 .

Vitesse de convergence : Considérons la matrice d'itération B symétrique, alors

B 2 =  (B ) et l'on obtient le résultat suivant : x k − x 2   (B ) . x0 − x 2 .
k

En conclusion, le rayon spectral de la matrice d'itération mesure la vitesse de

convergence. La convergence est d'autant plus rapide que le rayon spectral est plus
petit. On admettra que ce résultat ce généralise aux matrices quelconques.

13
 1) Méthode Jacobi
 −F
 
Notation : A =  D 
− E 
 
On pose A = M − N avec M = D et N = E + F selon le schéma ci-dessus. Ce qui donne la
formule de récurrence suivante : X n+1 = D −1 (E + F )X n + D −1b .

Soit xk = xki( ) i =1, n

le vecteur itéré, on donne l’algorithme donnant xk +1 en fonction de
xk .

Pour i de 1 à n faire
− 1  i −1 n 
x ki +1 =  a i , j x k + a x kj − bi 
j
i, j
ai ,i  j =1 j =i +1 

Théorème : Soit A une matrice symétrique, à diagonale strictement dominante, c'est

à dire a i ,i   a i , j alors la méthode est convergente car  (D −1 (E + F ))  1.
i j

A priori, cette méthode converge pour une matrice A SDP.

2) Méthode Gauss-Seidel
Boucle directe

On reprend le schéma de matrice précédent et on pose M = D − E et N = F . La

formule de récurrence s'écrit : (D − E )X n +1 = FX n + b .
i −1 n
On calcule x i
k +1 , l'ième composante du vecteur X k +1 : ai ,i x i
k +1 = − a i , j x j
k +1 − a i, j xkj + bi .
j =1 j =i +1

On propose l’algorithme suivant :

Pour i de 1 à n faire
1  i −1 n 
x ki +1 = −  ai , j x k +1 −  ai , j x k + bi 
j j

ai ,i  j =1 j =i +1 
On procède par écrasement du vecteur X k :

14
 t
 
 x j , 1  j  i −1 x ki +1 xk , i + 1  j  n  .
j
 k +
1
  
 déjà calculés calcul en cours 

Boucle rétrograde

Cette fois, on pose M = D − F et N = E .

On propose l’algorithme suivant :

Pour i de n à 1 faire
1  i −1 n 
x ki +1 = −  ai , j x k −  ai , j x k +1 + bi 
j j

a i ,i  j =1 j =i +1 

Conclusion
Théorème : Soit A une matrice symétrique définie positive, alors la méthode est
(
convergente car  (D − E ) F  1 .
−1
)

Comparaison : Si A est SDP les méthodes de Jacobi et de Gauss-Seidel convergent, et

Gauss-Seidel converge plus vite que Jacobi car
 (J ) =  (L)2 avec J = D −1 (E + F ) et L1 = (D − E )−1 F ou (D − F )−1 E .

3) Méthode de relaxation
D 1− w 
Pour cette méthode, on pose M = − E et N = M − A =   D + F , avec w  0 . On
w  w 
−1
D   1 − w  
donne la matrice d’itération : Lw =  − E    D + F  . Ainsi pour w = 1, on se
w   w  
ramène à la méthode de Gauss-Seidel.

Théorème : Soit A une matrice SDP. Il y a convergence si et seulement si  (Lw )  1 . On

démontre que ceci est vrai si et seulement si w  0,2 .

Remarque :

- w = 1 : Gauss-Seidel
- w  1 : Sous-relaxation
- w  1 : Sur-relaxation

15
Boucle directe
i −1
a i ,i 1  n
On veut calculer x ki +1 : x ki +1 = − ai , j x kj+1 − a
x kj +  − 1ai ,i x ki + bi . Comme dans
i, j
w j =1 j =i +1 w 
Gauss-Seidel, on note que l'algorithme procède par écrasement.

Pour i de 1 à n faire
w  i −1 n  w
x ki +1 = (1 − w)x ki −  ai , j x k +1 +
j
a i, j x kj  + bi
ai ,i  j =1 j =i +1  a i ,i

Choix de w :Soit A une matrice symétrique définie positive, et tridiagonale par blocs :
 O
 
 
A= .
 
O 
 
On cherche un w optimal, c'est à dire qui rend  (M −1 N ) minimum. On établit dans le

cas des matrices tridiagonales : woptimal =

2
( )
et  Lwoptimal = woptimal − 1 .
1 + 1 −  (L1 )
2

Note : On utilise aussi cette méthode de façon heuristique pour des matrices qui ne
sont pas tridiagonale.

Dichotomie : Si l’on cherche woptimal par dichotomie, l’expérience montre qu’il faut
obligatoirement procédé par valeurs supérieures (la fonction w   (Lw ) possède
pour w  woptimal une pente de 1).

Évaluation du rayon spectral de L :

 x k +1 = L1 .x k + c
On écrit les itérés de Gauss-Seidel :  d’où xk +1 − xk = L1 .(xk − xk −1 ) . On
 x k = L1 .x k −1 + c
x k +1 − x k
utilise la méthode de la puissance itérée et on déduit :  (L1 ) = lim . En
k → x − x
k k −1

pratique, on fait l’approximation pour des itérés dont l’écart a diminué de manière
significative.

Schéma d’implémentation :
- Soit x0 donné

16
- On réalise p itérations avec Gauss-Seidel, on obtient la suite des itérées jusqu’à xp
à partir de laquelle on évalue de manière approchée  (L1 ) et woptimal .
- Puis on applique la méthode de relaxation avec w = woptimal .
- On s’arrête lorsque l’on a atteint la précision désirée.

17
 Chapitre 2 Résolution numérique de problèmes
différentiels

Nombres de problèmes physiques aboutissent à des équations différentielles ou à un

système d'équations différentielles. La plupart des équations différentielles ne sont
pas solubles de manière analytique. Il faut donc les résoudre de manière numérique.

Généralités

Problème de Cauchy
Soit f : x0 , x0 + a R M → R M , continue, telle que (x, y )  f (x, y ). Etant donnée une
condition initiale y 0 , on cherche y : x0 , x0 + a → R M de classe C1 telle que
 y (x ) = f (x, y (x ))
 , problème de Cauchy dimension M et d'ordre 1.
 y (x0 ) = y 0

Le problème de Cauchy de dimension M et d'ordre P :

 (P )
(
 y (x ) = g x, y, y ,, y )
( P −1)
. On se ramène au problème de Cauchy


 0
y ( x ) = y 0 , y ( x 0 ) = y 
0 ,  , y ( P −1)
( x 0 ) = y ( P −1)
0

 y 
 
 y 
d'ordre 1 en posant z =  , de dimension PM .
 
 ( P −1) 
y 
 

Théorème (K - lipschitzienne)
Si il existe K tel que x, y1 , y 2 , f (x, y1 ) − f (x, y 2 )  K y1 − y 2 , alors le problème de
Cauchy d'ordre 1 possède une solution unique.

Méthode de résolution numérique

On se place dans le cadre de ce théorème, et l'on cherche à calculer une solution. On
divise l'intervalle x0 , x0 + a  en N de telle sorte que x N = x0 + a ; xn = x0 + nh avec le

. On voudrait que les y n soient proches de y ( x n ) . La formule générale

a
pas h =
N
employée pour la résolution numérique des équations différentielles est
 y n +1 = y n + h( x n , y n , h )
 .
 y0

18
Stabilité
 z n +1 = z n + h(x n , z n , h ) + h n
On introduit une perturbation,  , et l'on observe le
z0 = y0
comportement de la méthode. Posons  = sup  n . La méthode est stable si et
n

seulement si A / max z n − y n  A .
n

Consistance
La méthode est consistante si et seulement si
y ( x n +1 ) − y ( x n )
max −  ( x n , y ( x n ), h ) ⎯h⎯
⎯→ 0 . La méthode est consistante d'ordre p si
→0
n h
y ( x n +1 ) − y ( x n )
et seulement si A / max −  ( x n , y ( x n ), h )  Ah p .
n h

Convergence
La méthode est convergente si et seulement si max y n − y(xn ) ⎯h⎯
⎯→ 0 . Il y a
→0 n

convergence d'ordre p si et seulement si A / max y n − y (x n )  Ah p .

n

Si la méthode est stable et consistante [d'ordre p], alors la méthode est convergente
[d'ordre p].

Théorème stabilité
Si il existe K tel que x, y1 , y 2 , (x, y1 , h ) − (x, y 2 , h )  K y1 − y 2 , c'est-à-dire si  est K -
lipschitzienne, alors la méthode est stable.

Théorème consistance
- Si (x, y,0) = f (x, y ), on a consistance d'ordre 1.
   1  f 
- Si de plus,  ( x, y,0 ) =  ( x, y ) , alors on a consistance d'ordre 2.
 h  2  x 
  2  1  2 f 
- 
Si de plus,  2  ( x , y ,0 ) =  2 (x, y ) , alors on a consistance d'ordre 3.
  h  3  x 
- Etc. …

1) Méthode de Euler
(x, y, h) = f (x, y ) , ordre 1, stabilité acquise ( et f lipschitzienne). Intuitivement, on
y n+1 − y n
traduit la formule du taux d'accroissement f (xn , y n ) = . La formule sera donc
h
y n +1 = y n + hf (xn , y n ) .

19
  y' = y
Exemple : On considère l’E.D. 
 y ( 0) = 1

On la résout numériquement sur [0 ;5]

 yn +1 = yn + hyn

 y0 = 1

h = 0,1

n 1 10 20 30 40 50
xn 0,1 1 2 3 4 5
yn 1,1 2,59374 6,7275 17,4494 45,25926 117,39085
y(xn) 1,10517 2,71828 7,38906 20,08554 54,59815 148,41316
en 0,00517 0,12454 0,66156 2,63614 9,33889 31,02231

h = 0,01

n 1 100 200 300 400 500

xn 0,01 1 2 3 4 5
yn 1,01 2,70481383 7,31601785 19,7884663 53,5241172 144,772772
y(xn) 1,01005 2,71828 7,38906 20,08554 54,59815 148,41316
en 5E-05 0,01346617 0,07304215 0,29707374 1,07403279 3,64038757

h = 0,001

n 1 1000 2000 3000 4000 5000

xn 0,001 1 2 3 4 5
yn 1,001 2,71692393 7,38167565 20,0554512 54,4891355 148,042836
y(xn) 1,0010005 2,71828 7,38906 20,08554 54,59815 148,41316
en 5E-07 0,00135607 0,00738435 0,03008876 0,10901455 0,37032384

2) Méthode du développement de Taylor

1  f  1 h 2  f 
(x, y, h ) = f (x, y ) + h (x, y ) +  (x, y ) +  , consistante de l'ordre que l'on
2  x  3 2!  x 
veut, stabilité acquise. Mais en pratique les calculs de dérivées rendent la méthode
impraticable !

20
 3) Méthode de Runge - Kutta

Ordre 2
On cherche (x, y, h) = f (x, y ) + f (x + h, y + hf (x, y )) . L'ordre 1 impose  +  = 1 et
l'ordre 2 impose  =  et  =  = 1 / 2 .

- Si  = 1 , alors  = 0 et =  = 1 / 2 , d'où la formule

 
y n +1 = y n + hf  x n + , y n + f ( x n , y n ) . Soit le prédicteur y n +1 / 2 , la méthode s'écrit :
h h
 2 2 

 y n +1 / 2 = y n + f (x n , y n )
h
 2 La méthode est d'ordre 2, et correspond à une
 y n +1 = y n + hf (x n + h / 2, y n +1 / 2 )

amélioration de Euler : le taux d'accroissement est calculé pour le point milieu
(xn + h / 2, y n+1 / 2 ) .

Méthodes de la tangente améliorée, d’Euler améliorée, de Heun

L’idée est de chercher une fonction F telle que

F ( x, y, h) = a1 f ( x, y) + a2 f ( x + p1h, y + p2hf ( x, y))
Pour que la méthode soit convergente d’ordre 2 , on doit avoir :
 F ( x, y,0) = a1 f ( x, y ) + a2 f ( x, y ) = f ( x, y )

 F ( x, y,0) = a p df ( x, y ) + a p df ( x, y ) f ( x, y ) = 1 f (1) ( x, y ) = 1  df ( x, y ) + df ( x, y ) f ( x, y ) 
 h 2  dx 
2 1 2 2
 dx dy 2 dy 
a1 + a2 = 1

 1
 a2 p1 = a2 p2 =
2
Si on pose a2 =  , on obtient : F ( x, y, h) = (1 −  ) f ( x, y ) + f ( x +
h h
,y+ f ( x, y ))
2 2
On a donc une méthode convergente d’ordre 2.

h h
Si  = 1, on a la méthode de la tangente améliorée : F ( x, y, h) = f ( x + , y + f ( x, y ))
2 2

1
Si  = , on a la méthode d’Euler améliorée :
2
1 1
F ( x, y, h) = f ( x, y ) + f ( x + h, y + hf ( x, y ))
2 2

21
 3 1 3 2h 2h
Si  = , on a la méthode de Heun : F ( x, y, h) = f ( x, y ) + f ( x + , y + f ( x, y ))
4 4 4 3 3

1 3 2 2h
yn +1 = yn + h f ( xn , yn ) + h f ( xn + h, yn + f ( xn , yn )) 
4 4 3 3

 yn 0 = yn
 2h
 yn1 = yn + f ( xn , yn )
 3
 y = y + 1 h f ( x , y ) + 3 h f ( x + 2 h, y )
 n +1 n
4
n n
4
n
3
n1

 y' = y 2 + x
Exemple : On considère l’E.D.  On la résout numériquement sur [0 ; 1]
 y ( 0) = 0

 yn +1 = yn + hF ( xn , yn , h)

 y0 = 0

  
2

( ) h 
( ) 1 h
( )
2
h
F ( x, y, h) = (1 −  ) y 2 + x +   x + x + y 2   = y 2 + x + h + y 2 + x 
2
+y+ y2 + x
 2  2   2  4


Tangente améliorée h alpha

0,1 1
n 0 1 2 3 4 5
xn 0 0,1 0,2 0,3 0,4 0,5
yn 0 0,005 0,02001 0,04510018 0,08046261 0,1264784

n 5 6 7 8 9 10
xn 0,5 0,6 0,7 0,8 0,9 1
yn 0,1264784 0,18379726 0,2534407 0,33694681 0,43658688 0,5557065

22
Euler amélioré h alpha
0,1 0,5
n 0 1 2 3 4 5
xn 0 0,1 0,2 0,3 0,4 0,5
yn 0 0,005 0,02002376 0,04520432 0,08087513 0,12767911

n 5 6 7 8 9 10
xn 0,5 0,6 0,7 0,8 0,9 1
yn 0,12767911 0,18670939 0,25971218 0,34940505 0,46001149 0,59821583

Heun h alpha
0,1 0,75
n 0 1 2 3 4 5
xn 0 0,1 0,2 0,3 0,4 0,5
yn 0 0,005 0,02001424 0,04513564 0,08061261 0,12693526

n 5 6 7 8 9 10
xn 0,5 0,6 0,7 0,8 0,9 1
yn 0,12693526 0,18494123 0,25595927 0,34202165 0,4461994 0,5731619

Méthode de Runge Kutta d’Ordre 4

 ( x, y , h ) = (P1 + 2 P2 + 2 P3 + P4 )
1
-
6
- P1 = f (x, y )
- P2 = f (x + h / 2, y + h / 2.P1 )
- P3 = f (x + h / 2, y + h / 2.P2 )
- P4 = f (x + h, y + h.P3 )

Choix de h
On se donne une tolérance  > 0.

- On choisit arbitrairement h et on applique la méthode.

- On recommence avec h / 2 . Si le résultat obtenu pour un pas divisé par 2, est à 
près celui obtenu pour un pas de h alors le pas est bon, sinon on réitère. Attention
h
il faut comparer xn = x0 + nh et x 2 n = x 0 + 2n. .
2

23
On a également un système non linéaire dont on déduit

1 1
1 = ;  2 = ; 3 = 1
2 2
1 1 1 1 1
a10 = ; a20 = 0 ; a21 = ; a30 = 0 ; a31 = 0 ; a32 = 1 ; a40 = ; a 41 = ; a42 = ;
2 2 6 3 3
1
a43 =
6

 yn 0 = yn
 1
 yn1 = yn + h f ( xn , yn 0 )
 2
 1 h
 yn 2 = yn + h f ( xn + , yn1 )
 2 2
 yn 3 = yn + h f ( xn + h, yn 2 )
 1 1 h 1 h 1
 yn +1 = yn + 6 h f ( xn , yn ) + 3 h f ( xn + 2 , yn1 ) + 3 h f ( xn + 2 , yn 2 ) + 6 h f ( xn + h, yn 3 )

h 2 h3 h 4
Si on reprend notre exemple, on trouve yn +1 = (1 + h + + + ) yn d’où yn =
2 6 24
(1,10517083334)n

h = 0,1

n 1 10 20 30 40 50
xn 0,1 1 2 3 4 5
yn 1,10517083 2,71827974 7,38904477 20,0854907 54,5979826 148,41259
y(xn) 1,10517 2,71828 7,38906 20,08554 54,59815 148,41316
en 8,3334E-07 2,557E-07 1,5232E-05 4,9277E-05 0,00016741 0,00056985

h=1

1 2,70833333
n 0,1 1 2 3 4 5
xn 0,1 1 2 3 4 5
yn 1,10476577 2,70833333 7,33506944 19,8658131 53,8032438 145,717119
y(xn) 1,10517 2,71828 7,38906 20,08554 54,59815 148,41316
en 0,00040423 0,00994667 0,05399056 0,21972692 0,79490625 2,6960415

24
  y' = y 2 + x
Exemple : On considère l’E.D.  On la résout numériquement sur [0 ; 1]
 y ( 0) = 0

 yn 0 = yn
 1
 yn1 = yn + h f ( xn , yn 0 )
 2
 1 h
 yn 2 = yn + h f ( xn + , yn1 )
 2 2
y
 n3 = y n + h f ( xn + h , yn 2 )
 1 1 h 1 h 1
 yn +1 = yn + 6 h f ( xn , yn ) + 3 h f ( xn + 2 , yn1 ) + 3 h f ( xn + 2 , yn 2 ) + 6 h f ( xn + h, yn 3 )

RK4 h
0,1
n 0 1 2 3 4 5
xn 0 0,1 0,2 0,3 0,4 0,5
yn 0 0,00500188 0,02001485 0,04510414 0,08045029 0,12641868
yn1 0 0,01000313 0,03003488 0,06020586 0,1007739 0,15221776
yn2 0,0025 0,01250688 0,03255996 0,06278537 0,10345805 0,15507719
yn3 0,01000063 0,0155177 0,03212236 0,06000875 0,09956567 0,15146544

n 5 6 7 8 9 10
xn 0,5 0,6 0,7 0,8 0,9 1
yn 0,12641868 0,18364137 0,25312301 0,33638955 0,43571188 0,55446093
yn1 0,15221776 0,21532758 0,29132657 0,38204744 0,49020412 0,61983227
yn2 0,15507719 0,21845967 0,29486657 0,38618756 0,49522688 0,62617053
yn3 0,15146544 0,21677797 0,29712994 0,39494259 0,51380803 0,65911597

25
 Td Analyse Numérique

Exercice 1

1) Résoudre le système suivant à l’aide de la méthode de Gauss

4x -2y +3z = 9
2x - y +z = 3
x + 3y -z = 4

2) Résoudre le système suivant à l’aide de la méthode de Gauss-Jordan

x + y +z + w = -2

x + 2y +3z + 4w = -10

2x + 3y +4z + w = 4

3x + 4y +z + 2w = -10

Exercice 2

1 1 1 1
Soit la matrice A= (1 5 5 5)
1 5 14 14
1 5 14 14

1. Démontrer que A est symétrique et définie positive

2. Faites la décomposition de Cholesky de la matrice A

Exercice 3

1) Définir conditionnement d’un système

2) Considérons un système linéaire Ax=b.

Si l'on modifie sensiblement les paramètres de la matrice A, ou du second membre b,

on peut trouver des résultats complètement différents !

En altérant légèrement b en b+db et en notant x+dx la solution du nouveau système,

démontrer que :

26
‖𝑑𝑥‖ ‖𝑑𝑏‖
‖𝑥‖
≤ 𝜌(𝐴) ‖𝑏‖
où 𝜌(𝐴) est la mesure du conditionnement de A.

Exercice 4

6 1 1 𝑥 12
Soit le système linéaire AX=b avec A=(2 4 0) X=( ) b=( 0 )
𝑦
1 2 6 𝑧 6

1. résoudre le système AX=b par la méthode de Gauss

2
2. Approcher la solution par Jacobi avec 3 itérations à partir de X0=(2)
2
2
3. Approcher la solution par Gauss-Seidel avec 3 itérations à partir de X0=(2)
2
Exercice 5

La vitesse d’un corps en chute libre sous l’action de la pesanteur dans l’atmosphère
est solution de l’équation différentielle :
mv ′ (t) = −kv 2 (t) + mg ; 0 ≤ t ≤ 20s
{ avec v en m /s ; m=70kg ; g=9.8 ; k=0.27kg/m.
v(0) = 0

Estimer cette vitesse à l’aide de la méthode d’Euler avec le pas h=1.

Exercice 6

y′(x) = 2y(x)
On considère l’équation différentielle : {
y(0) = 1

1) Déterminer la solution exacte de cette équation différentielle.

2) On pose h=0.001

En utilisant l’Algorithme d’Euler, estimer y aux points xi=ih ;

i= 0…4

3) Refaire avec la méthode d’Euler Amélioré

Exercice 7

27
 y′(x) = y(x)𝑒 2𝑥
On considère l’équation différentielle : {
y(0) = 1

On pose h=0.1

1) Déterminer la solution exacte de cette équation différentielle.

2) En utilisant l’Algorithme d’Euler, estimer y aux points xi=ih ;

i= 0…4

3) En utilisant l’algorithme de Heun, estimer y aux points xi=ih ;

i= 0…4

28