Chapitre Premier

Résolution de système linéaire.

En générale, les systèmes linéaires sont issues de problèmes physiques complexes (équations
aux dérivées partielles ou équations différentielles) et la dimension de la matrice est souvent très
grande. Par conséquent, la résolution peut avoir un coup de calcul plus ou moins important suivant
la nature/structure de la matrice. Pour cette raison, il existe différentes méthodes numériques dans
la littérature. Nous allons nous intéresser à la factorisation LU d’une matrice A afin de résoudre un
système Ax = f . Une telle factorisation est toujours possible quand la matrice A est symétrique
définie positive.
La factorisation LU consiste à écrire une matrice non- singulière A comme le produit de deux
matrices L et U. L est une matrice triangulaire inférieure et U est une matrice triangulaire supé-
rieure. Il s’agit de décomposer une matrice A deux matrices de type L et U telles que A = LU .

1.1 Matrices triangulaires.

Une matrice carrée A = (aij )1≤i,j≤n est dite triangulaire inférieure si aij = 0 pour tout i < j.
Une triangulaire inférieure généralement sous la forme :
a11 0 0 ··· 0
 
a
 21 a22 0 ··· 0  
 a31 a32 a33 · · · 0
 

 
A= ,
 
 . .. .. .. .. 
 .. . . . . 
 
0 
 

an1 an2 an3 . . . ann
Une matrice carrée A = (aij )1≤i,j≤n est dite triangulaire supérieure si aij = 0 pour tout j < i.
Une triangulaire supérieure généralement sous la forme :
a11 a12 a13 · · · a1n
 

22 a32 · · ·
 0 a a2n 
 
 0 0 a · · · a
 
 33 3n 

A= ,
 
 . . . . .
 .. . . . .

 . . . .  
 
 
0 0 ... 0 ann

1.1.1 Inverse d’une matrice : Méthode de GAUSS

Une opération élémentaire sur une matrice consiste à permuter deux lignes ou deux colonnes ou
remplacer une ligne (resp. colonne) par une combinaison linéaire de plusieurs lignes (resp. colonnes).

1
 1.1. MATRICES TRIANGULAIRES. 2
Une matrice élémentaire est une matrice obtenue en faisant une seule opération élémentaire sur
les lignes de la matrice unitaire In .

Proposition 1.1. Soit A une matrice de type n × m et E une matrice élémentaire de In . Alors
la matrice EA est égale à la matrice obtenue en effectuant les mêmes opérations sur A que sur In
pour obtenir E.

Proposition 1.2. Soit A une matrice carrée d’ordre n inversible. Il existe une famille de matrice
élémentaire E1 , E2 , · · · , Ek telles que :
(Ek · · · E2 E1 )A = In .
Comme chaque matrice élémentaire est inversible, on a :
A = E1−1 E2−1 · · · Ek−1 et A−1 = Ek−1 · · · E2−1 E1−1

La détermination de l’inverse d’une matrice en utilisant la proposition précédente est appelée

méthode de GAUSS.

Processus
Etape 1 : On construit une matrice de type n × (2n) sous la forme : (A|In )
Etape 2 : On effectue les opérations élémentaires sur la nouvelle matrice de taille n × 2n jusqu’à
obtenir une matrice de type n × 2n de la forme (In |B)
Etape 3 : On déduit l’inverse de A qui est A−1 = B.
 
1 2 1
Exemple 1.3. Déterminons l’inverse de la matrice A = 0 3 −1 ,
 

1 0 0
On pose
 
1 2 1 1 0 0
A=  0 3 -1 0 1 0 


1 0 0 0 0 1

1.1.2 Factorisation ou décomposition LU.

La factorisation ou la décomposition LU d’une matrice A de type n × m consiste à écrire la
matrice A sous la forme A = LU . Dans cette écriture la matrice triangulaire inférieure L est de type
n × n et ayant rien que des 1 sur sa diagonale principale. U est une matrice triangulaire supérieure
de type n × m.

Propriété 1.4. Soit une matrice A de type n×m. On suppose qu’il existe des matrices élémentaires
E1 , E2 , · · · , Ek telles que (Ek Ek−1 · · · E1 )A soit une matrice triangulaire supérieure. On effectue les
mêmes opérations correspondantes aux matrices élémentaires sur la matrice In dans le même ordre.
On note B cette nouvelle matrice. On a évidemment B = Ek Ek−1 · · · E1 .
La matrice B est inversible et est égale à la matrice inverse de L. On a : L = B −1
Alors la décomposition
 LU de A est donnée par :

 U = (Ek Ek−1 · · · E1 )A
A = LU avec  B = (Ek Ek−1 · · · E1 )

L = B −1
 
1 2 1
Exemple 1.5. Déterminons la factorisation LU de la matrice A = 0 3 −1 ,
 

1 0 0
1.2. APPLICATIONS AUX EDP. 3
1.1.3 Résolution des systèmes linéaires.
Soit un système d’équations linéaires défini par AX = b où A est une matrice de type n × m
sur R et b un vecteur de Rm .
Si
( A possède une décomposition LU alors le système linéaire peut se ramener sous la forme :
UX = Y
LY = b
L étant inversible, on détermine d’abord Y = L−1 b puis on résout le système U X = Y .
Les matrices triangulaires L et U auraient pu être inversées aisément en utilisant l’élimination
de Gauss-Jordan. Mais si l’on compte simplement le nombre d’opérations que cela représente pour
un système à n équations AX = b, on trouvera que la complexité algorithmique du calcul des
matrices inverses est supérieure, de sorte que si l’on veut résoudre ce système pour divers b, il est
plus intéressant de réaliser la décomposition LU une fois pour toute et d’effectuer les substitutions
de descente-remontée pour les différents b plutôt que d’utiliser l’élimination de Gauss-Jordan à de
multiples reprises.

Exercice 1.6. Utiliser la décomposition LU pour résoudre chacun des systèmes suivants :
  3x − 7y − 2z + 2t = a
 x + 2y + z = a  x + 2y − 3z = a


   −3x + 5y + z

= b
(S1 ) :  3y − z = b (S2 ) : −3x − 4y + 13z = b (S3 ) :
6x − 4y − 5t = c
2x + y − 5z
 

x = c 
= c 

−9x + 5y − 5z + 12t = d


1.2 Applications aux EDP.

Soit y : [a; b] −→ R une fonction de classe C 1 sur [a ; b]. Alors pour tout t ∈ [a; b] on a :
y(t + h) − y(t) y(t) − y(t − h) y(t + h) − y(t − h)
y ′ (t) = lim = lim+ = lim−
h,→0 h h,→0 h h,→0 2h
Soit n ∈ N∗ . On appelle discrétisation régulière de [a, b] à n pas ou n + 1 points l’ensemble des
points a + nh, pour n ∈ {0, 1, · · · , } où le pas h est donné par h = b−a n
Soit {ti }1≤i≤n une dicrétisation de [a ; b] et Dyi une approximation de y ′ (ti ). On appelle
1. différence finie progressive, l’approximation définie par :
Dynp = y(tn +h)−y(t
h
n)
= y(tn+1h)−y(tn )
2. différence finie rétrograde ou récursive, l’approximation l’approximation définie par :
n −h)
Dynr = y(t)−y(t
h
= y(tn )−y(t
h
n−1 )

3. différence finie centrée, l’approximation l’approximation définie par :

n −h)
Dync = y(t+h)−y(t
2h
= y(tn+1 )−y(t
h
n−1 )

Pour calculer l’erreur commise par ces approximations, on rappelle le développement de Taylor-
Lagrange.

Propriété 1.7. Soit y : [a; b] −→ R une fonction de classe C n+1 sur [a ; b] ( n + 1 fois dérivable
et la dérivée (n+1)ième de f est continue sur [a ; b]). Alors :
pour tout (x; y) ∈ [a; b]2 , il existe ϵ appartenant à l’intervalle ouvert défini par x et y tels que :
n
1 (k) 1
f (x)(x − y)k + f (n+1) (ϵ)(x − y)n+1
X
f (x) = f (y) + (1.1)
k=1 k! (n + 1)!

∀x ∈ [a; b], ∀h ∈ R∗ vérifiant (x + h) ∈ [a; b], alors il existe ϵ ∈] min{x, x + h}, max{x, x + h}[
tel quel
1.2. APPLICATIONS AUX EDP. 4

n
1 (k) 1
f (x)hk + f (n+1) (ϵ)hn+1
X
f (x + h) = f (y) + (1.2)
k=1 k! (n + 1)!

Le développement de Taylor-Lagrange s’écrit aussi sous la forme

n
1 (k)
f (x)hk + O(hn+1 )
X
f (x + h) = f (y) + (1.3)
k=1 k!

où l’expression O(hn+1 ) se lit grand O de hn+1 . C’est-à-dire qu’il existe M > 0 et C > 0 tels
1
que (n+1)! f (n+1) (ϵ)hn+1 ≤ C|hn+1 | ∀h ∈] − M ; M [

1.2.1 Problème de Cauchy

Soit f une application définie par
f : [a ; b] ×Rm −→ Rm
(t,y) 7−→ f(t,y)
y: [a ; b] −→ Rm
Le problème de Cauchy consiste à trouver une fonction dérivable telle
t 7−→ y(t)
que

(
y ′ (t) = f (t, y) ∀t ∈ [a; b]
(1.4)
y(a) = y0

Les données principales de ce problème sont : la fonction f et la condition au bord (y(a) = y0 )

Dans de nombreux cas, la variable t représente le temps et les composantes du vecteur y, une
famille de paramètres décrivant l’état d’un système matériel donné. L’équation différentielle (1.4)
traduit physiquement la loi d’évolution du système considéré.
En général, il est impossible de trouver analytiquement des solutions à ces problèmes. Il faut
alors construire des méthodes numériques pour obtenir des solutions approchées. Toutefois, il serait
bon de vérifier l’existence et l’unicité d’une solution.
Considérons le problème suivant :

1
(
y ′ (t) = (y(t)) 3 ∀t ∈ [0; +∞[
(1.5)
y(0) = 0
1
Ce problème est un problème de Cauchy en dimension 1 dont f (t, y) = y 3 . Le problème possède
3 solutions qui sont :q q
8t3 3
y(t) = 0 ; y(t) = 27
et y(t) = − 8t27
Donc le problème ne possède pas une solution unique.
Le théorème suivant assure l’existence et l’unicité sous certaines conditions en dimension 1.

Théorème 1.8. Soit le problème de Cauchy donné par la définition par la relation (1.4) On suppose
que la fonction f est continue sur un ouvert U de R × R et quelle est localement lipschitzienne en y.
C’est-à-dire pour tout couple (t, y) appartenant à U , il existe un voisinage W de t et un voisinage
V de y, il existe L > 0 tels que :
1.2. APPLICATIONS AUX EDP. 5

∀s ∈ W, ∀(u, v) ∈ V 2 , ∥f (s, u) − f (s, v)∥ ≤ L∥u − v∥

(1.6)

Sous ces hypothèses le problème de Cauchy (1.4) admet une unique solution.
∂f (t,y)
Propriété 1.9. Si ∂y
est continue et bornée, alors f satisfait la condition de Lipschitz.

Exercice 1.10. Pour chacune des E.D.O. suivantes écrire le problème de Cauchy associé
′′
+ αy ′ (t) + βcos(y(t)) = sin(t) ∀t ∈]0; +2π]


 y (t)
y(0) = 0
 ′

y (0) = 1
′′ L ′ R1

 LCy (t) + ( R2 + R1 C)y (t) + ( R2 + 1)y(t) = e ∀t ∈]0; +100]


 ′
y(0) = 0

y (0) = 0
′′ 2 ′

 y (t) = µ(1 − y (t))y (t) + y(t) ∀t ∈]0; +10]

y(0) = 1
 ′

y (0) = 1
 (3)

 y (t) − cos(t)y (2) (t) + 2sin(t)y (1) (t) − y(t) = 0 ∀t ∈]0; +2π]
 y(0) = u0

′
 y (0) = v0


y””(0) = w0


Résolution numérique du problème de Cauchy.

On veut résoudre numériquement le problème de Cauchy en utilisant les différences finies pro-
gressive, rétrograde ou centrée.
On pose t0 , · · · , tn une discrétisation régulière de [a; b]. On a ti = t0 + ih avec h = b−a n
L’équation y ′ (t) = f (t, y(t)) est discrétisée par y ′ (ti ) = f (ti , y(ti )), i ∈ [|0, n − 1|]
En utilisant la formule des différences finies progressive, on a : ∀i ∈ [|0, n|], ϵi ∈ [ti ; ti+1 ],

h2
y(ti+1 ) − y(ti ) = hf (ti , y(ti )) + y”(ϵi ) (1.7)
2
On note yi une approximation de y(ti ) pour i ∈ [|0; n − 1|].
La méthode d’Euler progressive est donnée par le schéma

(
yi+1 = yi + hf (ti , yi ) i ∈ [|0; n − 1|]
(1.8)
y0 = y(t0 )

Ce schéma est explicite, car il permet le calcul direct de yi+1 en fonction de yi .

De la même manière, la méthode d’Euler régressive est donnée par le schéma :

(
yi+1 = yi + hf (ti+1 , yi+1 ) i ∈ [|0; n − 1|]
(1.9)
y0 = y(t0 )

Ce schéma est implicite, car yi est définit implicitement en fonction de yi+1 . Il faut donc résoudre
à chaque pas de temps une équation non-linéaire en utilisant des méthodes de point fixe par exemple.
1.3. NOTIONS DE CONDITIONNEMENT ET STABILITÉ 6
Exercice 1.11. On veut résoudre
( numériquement le problème (P) suivant qui consiste à trouver
y ′ (t) = cos(t) + 1 t ∈ [0; 4π]
ytelleque :to (1.10)dont la solution exacte est
y0 = 0
y(t) = sin(t) + t.
1. Expliquer en détail comment utiliser le schéma d’Euler progressif pour résoudre le problème
(P) en précisant entre autres les données, les inconnues, les dimensions des variables, lien
entre yi et la fonction y.
(
−u′′ (xi ) = f (xi )
2. Montrer que le problème initial écrit au point xi comme suit peut
u(0) = 0
(
Au = f
être approché par le système par le système linéaire suivant : (S) : avec
u(0) = 0
u = (u1 , u2 , · · · , un−i )T , f = (f1 , f2 , · · · , fn−i )T et A est une matrice à déterminer.
3. Résoudre (S) par la décomposition LU.

Exercice 1.12. Soit f (x) = x + 1.

−u′′ (x) = f (x) ∀x ∈]0; 1[



On se propose de résoudre l’équation différentielle suivante : (E) : u(0) = 0


u(1) = 0
1. On suppose que que toute solution u de l’équation (E) est de classe C 4 sur [0 ; 1]. Montrer
en utilisant la formule de Taylor que :

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

u′′ (x) = + h2 ε(x), ∀ x ∈]0, 1[, h , 0
h2
tel que x + h, x − h ∈]0; 1[ où ε(x) est une fonction bornée sur [0; 1].
2. On subdivise le domaine [0 ;1] en 4 domaines [xi ; xi+1 ] de façon régulière avec xi = ih et
h = 14 . On pose fi = f (xi ) et ui = u(xi ) pour i = 0, 1, 2, 3.
(a) Montrer que l’on a l’approximation
ui−1 − 2ui + ui+1
u′′ (x) =
h2

−u′′ (xi ) = f (xi )




(b) Montrer que le problème initial écrit au point xi comme suit u(0) = 0


u(1) = 0

 Au = f

peut être approché par le système par le système linéaire suivant : (S) : u(0) = 0


u(1) = 0
T T
avec u = (u1 , u2 , · · · , un−i ) , f = (f1 , f2 , · · · , fn−i ) et A est une matrice à déterminer.
(c) Résoudre (S) par la décomposition LU.

1.3 Notions de conditionnement et stabilité

Ces deux notions, toujours présentes en analyse numérique, sont relatives à la propagation plus
ou moins importante des erreurs d’arrondi dans un calcul donné. Nous les étudions ici pour le calcul
d’une fonction f x ∈ R 7→ f (x).
1.3. NOTIONS DE CONDITIONNEMENT ET STABILITÉ 7
1.3.1 Conditionnement
Le conditionnement décrit la sensibilité de la valeur d’une fonction à une petite variation de son
(x∗ ) ∗
argument, c’est-à-dire : f (x)−f
f (x)
en fonction de x−x
x
lorsque xx∗ est petit.
f (x)−f (x∗ )
f (x) xf ′ (x)
Pour une fonction suffisamment régulière, on a évidemment : x−x∗ ≃ f (x)
D’où on tire :
x

Définition 1.13. On appelle conditionnement d’une fonction numérique f de classe C 1 en un

point x, le nombre

xf ′ (x)
cond(f )x = (1.11)
f (x)
√
Exemple 1.14. Pour f (x) = x on a cond(f)x = 12 unbonconditionnement, puisquel′ erreurrelativesurf seraau
x
Pour f (x) = x − a on a cond(f )x = x−a ce qui est un mauvais conditionnement surtout si x est
voisin de a.

1.3.2 Stabilité
La stabilité décrit la sensibilité d’un algorithme numérique
√pour le calcul
√ d’une fonction fq(x).
x
Le conditionnement de cette fonction est égal à : f (x) = x + 1 − x est cond(f )x = 12 x+1
Cette dernière expression étant proche de 1 2pour xgrand.Donc, sixestgrand, leconditionnementdef estbon.Cepend
√ √
(12345) = 12345 + 1 − 12345 = 111.113 − 111.108 = 0.500000.10−2 .
Or un calcul précis donne : f (12345) = 0.4500032.... × 10−2 . On a donc une erreur de 10% ce
qui est important et peu en accord avec le bon conditionnement de f. Ceci est dû à l’algorithme
utilisé dans ce calcul que l’on peut expliciter comme suit :
x0 = 12345 ;
x1 = x0 + 1 ;
√
x2 = x0 ;
√
x3 = x1 ;
x4 = x3 − x2

Il y a quatre fonctions à intervenir et, à priori, même si le conditionnement de f est bon, il se

peut que le conditionnement d’une ou plusieurs fonctions utilisées dans l’algorithme soit supérieur
à celui de f. C’est ce qui se produit ici pour la fonction x3 7→ x4 = x3 − x2 (x3 étant supposé fixé)
dont le conditionnement est grand lorsque x2 est voisin de x3 comme on l’a vu précédemment.
En conclusion, le choix d’un bon algorithme numérique est essentiel. Par exemple, cidessus, un
meilleur algorithme
√ est
√obtenu en utilisant :
1 √
f (x) = x + 1 − x = √x+1+ x
On obtient toujours avec 6 chiffres significatifs :
f (12345) = √12345+1+1 √
12345
= 0.450002.10−2
ce qui donne une erreur relative de 0, 0003%. Dans la même optique, l’exemple suivant peut
aussi être médité.
N
X xk
Le calcul de e−2 en utilisant le développement de Taylor ex = = SN pour N assez grand.
i=0 k!
1.3. NOTIONS DE CONDITIONNEMENT ET STABILITÉ 8
N SN N SN N SN
2 -11,0... 19 1629,87... 36 -0,001432...
3 61,0... 20 -996,45... 37 0,000472...
4 -227,0... 21 579,34... 38 -0,0001454...
5 637,0... 22 -321,11... 39 0,000049726...
6 -1436,6... 23 170,04... 40 -0,000010319...
7 2710,6... 24 -86,20... 41 0,000007694...
8 -4398,88... 25 41,91... 42 0,000002422...
9 6265,34... 26 -19,58... 43 0,000003928...
10 -7953,62... 27 8,80... 44 0,000003508...
11 9109,137... 28 -3,8130... 45 0, 000003623...
12 -9504,78... 29 1,5937... 46 0,000003592...
13 9109,13... 30 -0,6435... 47 0,000003600...
14 -8072,94... 31 0,2513... 48 0,000003598...
15 6654,55... 32 -0,0950... 49 0,000003599...
16 -5127,44... 33 0,0348... 50 0,000003598...
17 3709,05... 34 -0,01238... 51
18 -2528,47... 35 0,004283.. 52
La valeur de e− 12 est en fait de 0,0000061442.... La formule e−x = e1x donne lieu à un calcul
plus stable, même si ex est calculé comme ci-dessus.
Afin de limiter la propagation des erreurs d’arrondi, il faut essayer d’anticiper en utilisant
des algorithmes dont la stabilité est optimisée par un choix d’opérations intermédiaires à bon
conditionnement.
Les phénomènes soulignés dans cette partie ne pouvant être, en tout état de cause, complètement
éliminés, on peut essayer d’évaluer l’erreur totale à laquelle un algorithme est susceptible de donner
lieu : en faisant un calcul en double précision et en confrontant le résultat au même calcul fait en
simple précision.
 Chapitre Deux

Résolution d’équations non linéaires

L’objectif de ce chapitre est d’étudier les algorithmes les plus fréquemment utilisés pour résoudre
des équations non linéaires du type f (x) = 0 où x peut-être une variable réelle, complexe ou
vectorielle et f est une fonction.

2.1 Cas des fonctions d’une variable

2.1.1 Préliminaires : séparation des zéros
La séparation des zéros et une phase d’analyse mathématique minimale pour séparer les zéros,
c’est-à-dire déterminer des intervalles [ai , bi ] dans lesquels l’équation considérée a une solution et
une seule.
Ce travail devient relativement facile lorsque la fonction est continue et dérivable. Dans ce cas
si elle est strictement monotone sur un intervalle [ai , bi ] et que f (ai ) × f (bi ) ≤ 0 l’équation f (x) = 0
admet une seule solution sur l’intervalle [ai , bi ].

Exemple 2.1. On veut résoudre l’équation f (x) = x − 0.2sinx − 0.5 = 0. Il est évident que f
est dérivable et f ′ (x) = 1 − 0.2cos(x) > 0 pour tout x ∈ R. f est donc continue et strictement
monotone sur tout intervalle de R. f (0) = −0.5 < 0, f (π) = π − 0.5 > 0 alors l’équation f (x) = 0
admet une unique solution appartenant à [o, π].

La deuxième étape de la résolution est mettre la main sur la solution ou de trouver une valeur
approchée à cette solution. Les algorithme suivants vont nous permettre de gérer cette partie qui
est en réalité la plus importante.

2.1.2 Méthode de dichotomie (ou bissection)

Reprenons l’exemple (2.1) Nous avons vu que si f (x) = x − 0.2sin(x) − 0.5, f (0) < 0, f (π) > 0 ;
d’où un zéro dans ]0, π[. Si on calcule la valeur de f en ( π2 ), on constate que f( π2 )= π2 − 0.7 > 0. Donc
ce zéro est en fait dans [0, π2 ]. On peut alors calculer f( π4 ) = 0.14 > 0, qui montre qu’il est dans
[0, π4 ]. Il est clair qu’une répétition de ce procédé donne un encadrement de plus en plus précis du
zéro cherché et fournit donc un algorithme de calcul de ce zéro.

9
2.1. CAS DES FONCTIONS D’UNE VARIABLE 10
Algorithme
Soit une fonction f : [a, b] −→ R continue telle que f (a)f (b) < 0
1. Pour n = 0, · · · , N Faire :
an +bn
(a) m = 2
;
(b) Si f (a)f (b) < 0 Faire :
i. an+1 = an ;
ii. bn+1 = m ;
(c) Sinon Faire :
i. an+1 = m ;
ii. bn+1 = bn ;
(d) Retourner en (b).
On a : an+1 − bn+1 = 12 (an − bn ), soit par récurrence an − bn = 21n (a0 − b0 ). Il en résulte que
cette méthode est toujours convergente puisque an − bn ) tend vers 0 quand n tend vers l’infini. On
peut choisir le temps d’arrêt N pour que : 21N (a0 − b0 ) < ε où ε est la précision choisie On obtient
alors un encadrement de la solution cherchée à la précision voulue. Bien que cette situation soit
très alléchante, on verra qu’en fait cette méthode converge plutôt lentement comparée à celles qui
suivent.

2.1.3 Méthode de la sécante

Soit f admettant un zéro dans l’intervalle [x−1 , x0 ]. Pour obtenir une première approximation
de ce zéro, l’idée est de remplacer f par son interpolé linéaire sur [x1, x0], soit par : P (x) =
f (x0 ) + (x − x0 ) f (xx00)−f (x−1 )
−x−1
l’unique fonction linéaire dont les valeurs coïncident avec celles de f
en x−1 , x0 . L’approximation x1 est alors obtenue en résolvant l’équation P (x) = 0 : Soit x1 =
−x−1
x0 − f (x0 ) f (xx00)−f (x−1 )
Géométriquement ceci revient à remplacer la courbe d’équation y = f (x) par la droite d’équation
y = f (x0 ) + (x − x0 ) f (xx00)−f (x−1 )
−x−1
et xn+1 est l’intersection de cette droite avec la droite des abscisses.
Comme le montre le dessin, xn+1 semble plus voisin du zéro cherché que xn−1 ou xn . Pour
trouver une meilleure approximation, il suffit de répéter le procédé à l’aide des points (xn , xn+1 ).
En continuant, on obtient le zéro recherché. L’Algorithme se présente comme suit :
x0
" et x1 étant donnés,
∀ n = 0, 1, · · · , N
−xn−1
xn+1 = xn − f (xn ) f (xxnn)−f (xn−1 )

Critère d’arrêt

Cet algorithme n’est évidemment pas complet tant qu’on n’a pas précisé un critère d’arrêt.
Nous verrons plus loin que, généralement, xn converge vers la solution x∞ cherchée ce qui signifie
que pour n grand, xn est voisin de x∞ . Sans plus d’informations, il est difficile de savoir ce que
signifie numériquement "n grand", et c’est là une difficulté fréquente en analyse numérique. Un
critère d’arrêt souvent utilisé consiste à choisir a priori une tolérance ε et à terminer l’algorithme
lorsque |xn − xn+1 | ≤ ε
2.1. CAS DES FONCTIONS D’UNE VARIABLE 11
2.1.4 Méthode de Newton
Ici, au lieu d’assimiler la courbe y = f (x) à une sécante, on l’assimile à la tangente en un point
(xn , f (xn )) soit, la droite d’équation : y = f (xn ) + f ′ (xn )(x − xn ).
Son intersection avec l’axe des abscisses fournit une approximation de la solution x∞ . A nouveau,
on renouvelle
 le procédé jusqu’à obtenir une approximation suffisante, d’où l’algorithme :
x0 ∈ R
 ∀ n = 0, 1, · · · , N


xn+1 = xn − ff′(x n)
(xn )
Cet algorithme est initié à l’aide d’un seul point. Il peut être vu comme une modification de
−xn−1
l’algorithme précédent où le quotient f (xxnn)−f (xn−1 )
a été remplacé par f ′ (x1 n ) Nous verrons plus loin
qu’il donne lieu à une convergence beaucoup plus rapide.

2.1.5 Méthode du point fixe.

Elle consiste à d’abord remplacer l’équation f (x) = 0 par une équation g(x) = x ayant mêmes
solutions. On est ainsi ramené à la recherche des points fixes de l’application g. Le remplacement
de f (x) = 0 par g(x) = x est toujours possible en posant, par exemple, g(x) = f (x) + x. Ce n’est
évidemment pas forcément le meilleur choix.
Comme√exemple p, pour résoudre l’équation x2 − x − 2 = 0, on peut considérer g(x) = x2 − 2
ou g(x) = x + 2 ou g(x) = 1 + x2
L’algorithme
 se présente alors sous la forme :
x0 ∈ R
 ∀ n = 0, 1, 2, · · · ,


xn+1 = g(xn )

Propriété 2.2. Soit g une fonction

( définie de I = [a, b] à valeurs dans I, continue sur I et x0 ∈ R.
x0 ∈ R
On définie la suite (xn )n par
xn+1 = g(xn ), ∀n ≥ 0
Si xn converge vers x alors g(x) = x

2.1.6 Étude de convergence.

Méthode du point fixe
Commençons par traiter le cas du point fixe qui est fondamental d’un point de vue théorique.

Propriété 2.3. Soit une fonction g définie de I = [a, b] à valeurs dans I et dérivable telle qu’il
existe K ∈ [0, 1[ tels que |g ′ (x)| ≤ K ∀ x ∈ [a, b]. (
x0 ∈ R
Alors, pour tout x0 ∈ [a, b], la suite définie par converge vers
xn+1 = g(xn ), ∀n ≥ 0
l’unique point fixe de g.

Preuve 2.4. Notons d’abord que la suite est définie pour tout n puisque ∀ x ∈ [a, b], g(x) ∈ [a, b].
La fonction g a un point fixe et un seul dans l’intervalle [a, b], car la fonction h(x) = g(x) − x
vérifie :
— h est continue sur [a, b]
— h(x) = g(x) − 1 ≤ k − 1 < 0, donc h est strictement monotone.
— h(a) = g(a) − a ≥ 0 puisque g(a) ∈ [a, b]
— h(b) = g(b) − b ≥ 0 puisque g(b) ∈ [a, b].
2.1. CAS DES FONCTIONS D’UNE VARIABLE 12
Soit x ce point fixe, on a

xn+1 − x = g(xn ) − g(x) = (xn − x)g ′ (εn ) (2.1)

où εn ∈]a, b[ d’après le théorème des accroissements finis. Ainsi, en utilisant le fait g ′ soit bornée,
on obtient

|xn+1 − x| ≤ K|xn − x| (2.2)

et par récurrence, on a
|xn − x| ≤ K n |x0 − x| ∀ n ∈ N. Puisque 0 ≤ K < 1, ceci prouve que n−→∞
lim |xn − x| = 0

Quand une suite xn converge vers x selon la relation (2.2), on dit que la convergence est au
moins linéaire. Plus généralement, on définit
Définition 2.5. Soit xn une suite convergeant vers x. S’il existe un nombre p et une constante
C , 0 tels que

|xn+1 − x|
lim =C
n−→∞ |xn − x|p
(2.3)

on dit que la convergence est d’ordre p ; C est la constante d’erreur asymptotique.

S’il existe une constante C et un réel positif p tel que
|xn+1 − x| ≤ C|xn − x|p
pour tout n assez grand, on dit que la convergence est au moins d’ordre p.
1
Remarque 2.6. Nous pouvons dire que la convergence vers 0 de an − bn (an − bn ≤ 2n
(a0 − b0 ))
dans l’algorithme de dichotomie est linéaire.
Remarque 2.7. Il est clair qu’une convergence est d’autant plus rapide que son ordre est grand.
En effet, si |xn − x| ∈]0, 1[, |xn − x|2 est encore plus petit que |xn − x|.

Méthode de Newton
Partons de la méthode du point fixe. Supposons que la fonction g est suffisamment régulière.
Alors il existe εn entre xn et x tel que
1
xn+1 − x = g(xn ) − g(x) = g ′ (x)(xn − x) + g”(εn )(xn − x)2 (2.4)
2
Pour n grand, on déduit que

xn+1 − x ≃ g ′ (x)(xn − x) (2.5)

et la vitesse de convergence sera d’autant plus grande que g(x) est plus petit. Le cas le plus favorable
est celui où g(x) = 0. Si M est majorant de g”(x) sur [a, b], on a alors :

M
|xn+1 − x| ≤ |xn − x|2 (2.6)
2
2.2. SYSTÈME D’ÉQUATIONS NON LINÉAIRES 13
Alors la convergence est au moins d’ordre 2. Revenons maintenant sur l’algorithme de Newton que
nous pouvons interpréter comme l’algorithme de point fixe sur la fonction g(x) = x − ff′(x) (x)
La dérivée de g est donnée par :
′
g ′ (x) = 1 − ff ′ (x)
(x)
+ f (x)f ”(x)
f ′2 (x)
= f (x)f ”(x)
f ′2 (x)
Si f ′ (x) , 0 et comme f (x) = 0 on a g ′ (x) = 0
Ceci montre que la méthode de Newton converge de façon quadratique si elle converge.
Pour s’assurer de sa convergence, on peut essayer d’appliquer à g la propriété 2.3. Il est clair
que l’hypothèse de la propriété 2.3 est vérifiée sur un voisinage suffisamment petit de x(puisque
g(x) = 0 et on suppose g continue). L’hypothèse plus difficile à vérifier est que g envoie un voisinage
[a, b] de x dans lui-même(g([a, b]) ⊂ [a, b]). On a en fait le résultat suivant :

Théorème 2.8. Soit f deux fois continûment différentiable sur un intervalle ouvert de centre x
vérifiant : f (x) = 0 et f ′ (x) , 0.
Alors, si x0 est choisi assez près de x, la méthode de Newton :
xn+1 = x − ff′(x n)
(xn )
,n≥0
produit une suite xn convergeant au moins quadratiquement vers x.

2.2 Système d’équations non linéaires

Reprenons les méthodes étudiées dans le cas des fonctions numériques d’une variable réelle. Il
est clair que la méthode de dichotomie n’a aucun sens en dimension n ≥ 2. En revanche, les trois
autres méthodes peuvent s’adapter au cas d’une fonction F Rn 7→ Rn .
Passons rapidement sur les méthodes de point fixe qui nécessitent que le système d’équations à
résoudre soit écrit sous la forme



 x1 = g(x1 , x2 , · · · , xn )
x2 = g(x1 , x2 , · · · , xn )



(S) :  .. .. .. (2.7)



. . .
xn = g(x1 , x2 , · · · , xn )


Il est assez rare que la méthode de point fixe écrite à partir du système (2.7) converge, les
conditions de convergence étant plus difficiles à réaliser qu’en dimension un. Il faut en effet que pour
une certaine norme matricielle, la norme de la matrice Jacobienne de F soit strictement inférieure à
1 au point recherché. On peut bien sûr envisager des méthodes d’accélération de convergence comme
dans les cas précédents. Nous faisons le choix de présenter plutôt ici les méthodes généralisant la
méthode de Newton et celle de la sécante.

2.2.1 La méthode de Newton-Raphson

Considérons un système d’équations s’écrivant F(X) = 0 ou encore, de façon développée



x1 = g(x1 , x2 , · · · , xn )
x2 = g(x1 , x2 , · · · , xn )



(S) : .. .. .. (2.8)




. . .
xn = g(x1 , x2 , · · · , xn )


La généralisation naturelle de la formule de Newton unidimensionnelle

2.2. SYSTÈME D’ÉQUATIONS NON LINÉAIRES 14

xn+1 = xn − (f ′ (xn ))−1 f (xn ) (2.9)

fait intervenir la matrice Jacobienne JF de F en xn . La méthode de Newton-Raphson s’écrit

donc formellement

Xn+1 = Xn − [JF (Xn )]−1 F (Xn ). (2.10)

Le second membre de (2.10) est parfaitement défini car on effectue le produit d’une matrice
n × n par un vecteur de Rn . Dans la pratique, on ne calcule pas explicitement l’inverse de la
matrice Jacobienne, ce qui s’avèrerait trop coûteux, et on préfère écrire l’algorithme sous la forme
suivante :


X0 ∈ Rn ;
∀ n = 0, 1, · · · , N




 Tester : l’arrêt puis faire (2.11)

 Résoud : [JF (Xn )]Yn = F (Xn )
Xn+1 = Xn + Yn
Remarque 2.9. 1. Encore plus qu’en dimension 1, le choix de l’initialisation est crucial et le
risque de divergence, si on ne démarre pas à proximité de la solution cherchée, est grand.
2. La convergence est là aussi d’ordre 2, donc très rapide (quand il y a convergence !)
3. Dans le cas des systèmes de grande dimension, la méthode de Newton-Raphson s’avère assez
coûteuse puisqu’il faut évaluer à chaque itération n2 + n fonctions (les n2 dérivées partielles
de la matrice Jacobienne, plus les n fonctions coordonnées). De plus, il y a un système
linéaire n × n (dont la matrice est en général pleine) à résoudre. Pour pallier le premier
inconvénient, il est fréquent qu’on ne recalcule pas la matrice Jacobienne à chaque itération,
mais seulement de temps en temps. Bien sûr, à ce moment-là, la convergence n’est plus
quadratique, mais en temps de calcul il arrive qu’on y gagne beaucoup, en particulier quand
les dérivées de F sont très coûteuses à calculer.

2.2.2 La méthode de Broyden

La méthode de Broyden s’inspire directement de la remarque ci-dessus en poussant même la
logique un peu plus loin. Puisqu’il est coûteux (voire même dans certains cas impossible) de calculer
la matrice Jacobienne de F, on va se contenter d’en déterminer une valeur approchée à chaque
itération Bn . De plus la suite de matrices Bn se calculera simplement par récurrence.
On peut considérer que la méthode de la sécante, en dimension un, est basée un peu sur le
même principe. En effet elle revient à approcher
f ′ (xn ) par f (xxnn)−f (xn−1 )
−xn−1
En dimension n, on va simplement imposer à la suite de matrices Bn de vérifier la même relation :

Bn (Xn − Xn−1 ) = F (Xn ) − F (Xn−1 ). (2.12)

Bien sûr, la relation (2.12) ne définit pas la matrice Bn . Elle n’impose que sa valeur dans une
direction. Broyden a fait le choix simple de considérer qu’on pouvait passer de Bn à Bn+1 en
rajoutant simplement une matrice de rang 1. Ainsi, sa formule d’actualisation s’écrit :
2.2. SYSTÈME D’ÉQUATIONS NON LINÉAIRES 15

δFn − Bn δXn )(δXn )T

Bn−1 = Bn + ). (2.13)
δ∥Xn ∥2

où on a noté δXn = Xn+1 − X − n et δFn = F (Xn+1 ) − F (Xn ). On vérifie immédiatement que

la suite de matrice Bn définie par (2.13) satisfait bien (2.12). La méthode de Broyden s’écrit donc :

X0 , B0 sont donnés;


 ∀ n = 0, 1, · · · , N
Tester : l’arrêt puis faire


 (2.14)

 Résoud : Bn δn = −F (Xn )
T
Bn−1 = Bn + δFn −B∥X n δXn )(Xn )
n∥
2 )

Remarque 2.10. 1. On peut prendre comme matrice initiale B0 = Id, il faut évidemment
attendre un certain temps avant d’avoir une approximation raisonnable de la matrice Jaco-
bienne.
2. On peut démontrer qu’en général, comme pour la méthode de la sécante, la convergence est
superlinéaire. La suite de matrices Bn ne converge pas forcément vers la matrice Jacobienne
de F.

Exercice 2.11. On doit résoudre l’équation x2 − 2 = 0 par la méthode de point fixe.

1. Justifier que chacune des fonctions suivantes peut-être utilisée pour le calcule des solutions.
g1 : x 7→ x2
g2 : x 7→ 12 (x + x2 )
g3 : x 7→ 2x − x2
g4 : x 7→ 13 (2x + x2 )

2. Écrire l’algorithme du point fixe en utilisant la fonction g2 .

Exercice 2.12. Pour les fonctions suivantes, trouver un algorithme de point fixe (autre que la
méthode de Newton !) qui converge vers le plus petit zéro positif :
g1 : x 7→ x − x − 1 = 0
g1 : x 7→ x − tan(x) = 0
g1 : x 7→ e−x − cos(x) = 0
 Chapitre Trois

Interpolation et approximation
polynomiale.

Dans ce chapitre, on dispose d’une fonction f, connue par exemple uniquement par ses valeurs
en certains points, et on cherche à remplacer ou à approcher f par une fonction plus simple, le plus
souvent par un polynôme. Nous verrons dans ce contexte, l’interpolation qui consiste à rechercher
un polynôme qui passe exactement par les points donnés, l’interpolation par morceaux, en parti-
culier les fonctions splines où l’expression du polynôme est différente sur chaque sous-intervalle,
l’approximation uniforme ou l’approximation au sens des moindres carrés ou on cherche à appro-
cher au mieux la fonction, soit pour la norme uniforme (ou norme infinie), soit pour la norme
euclidienne. Cette dernière approche conduira naturellement aux fonctions trigonométriques et aux
séries de Fourier, et on présentera la Transformée de Fourier Rapide ou F.F.T.

3.1 Interpolation de Lagrange

3.1.1 Le polynôme d’interpolation
Soit f : [a, b] −→ R connue en n + 1 points distincts x0 , x1 , · · · , xn de l’intervalle [a, b]. Il s’agit
de construire un polynôme P de degré inférieur ou égal à n tel que

∀ i = 0, 1, · · · , n, P (xi ) = f (xi ). (3.1)

Théorème 3.1. Il existe un et un seul polynôme de degré inférieur ou égal à n solution de (3.1).
Le polynôme s’écrit

n
X
Pn (x) = f (xi )Li (x) (3.2)
i=0

avec

(x − xk )
Li (x) = Πnk=0,k,i (3.3)
(xi − xk )
Le polynôme Pn est appelé polynôme d’interpolation de Lagrange de la fonction f aux points
x0 , x1 , · · · , xn .
Les polynômes Li (x) sont appelés polynômes de base de Lagrange associés à ces points.
Pour n = 1 on a
(x−x1 ) (x−x0 ) (x−x1 ) (x−x0 )
L0 (x) = (x 0 −x1 )
, L1 (x) = (x 1 −x0 )
et P1 (x) = f (x0 ) (x 0 −x1 )
+ f (x1 ) (x 1 −x0 )

16
3.1. INTERPOLATION DE LAGRANGE 17
P1 peut s’écrit comme :
f (x1 ) − f (x0 )
P1 (x) = f (x0 ) + (x − x0 ) (3.4)
(x1 − x0 )

Remarque 3.2. L’écriture (3.2) du polynôme d’interpolation est intéressante au point de vue
théorique, mais peu du point de vue numérique : elle a un caractère peu algorithmique. De plus, son
évaluation requiert trop d’opérations élémentaires. On lui préfère la formule de Newton qui consiste
à écrire le polynôme d’interpolation Pn aux points x0 , x1 , · · · , xn sous une forme généralisant (3.4)
soit

Pn (x) = an0 + an1 (x − x0 ) + · · · + ann (x − x0 )(x − x1 ) · · · (x − xn−1 ) (3.5)

Notons d’abord que tout polynôme de degré inférieur ou égal à n peut se mettre sous cette
forme dès que les xi sont tous distincts. L’avantage de cette écriture est que la partie tronquée
an0 + an1 (x − x0 ) + · · · + ann (x − x0 )(x − x1 ) · · · (x − xn−2 )
n’est pas autre chose que Pn−1 (x) : en effet, il s’agit d’un polynôme de degré inférieur ou égal
à n − 1 tel que :
pour tout i = 0, 1, · · · , n − 1 la valeur en xi soit égale à Pn (xi ) = f (xi ).
Donc, connaissant Pn−1 , il suffit de calculer ann pour connaître Pn (ann est aussi le coefficient de
xn dans Pn écrit sous forme usuelle). Les ain sont donnés par la formule de Newton :
Théorème 3.3. (Formule d’interpolation de Newton)
Le polynôme d’interpolation de Lagrange de la fonction f aux points distincts x0 , x1 , · · · , xn est
donné par :

n
i−1
X
Pn (x) = f [x0 , x1 , · · · , xi ]Πk=0 (x − xk )
i=0

(3.6)

où f [.] désigne les différences divisées de f définies par :

i) f [xi ] = f (xi ) ∀ i = 0, 1, · · · , n (3.7)

f [x0 , x1 , · · · , xk−1 ] − f [x1 , · · · , xk ]
ii)f [x0 , x1 , · · · , xk ] = (3.8)
xk − x0
Propriétés des différences divisées
— f [x0 , x1 , ..., xn ] est invariant par des permutations sur les xi : en effet, permuter les xi ne
change pas le polynôme d’interpolation et donc ne change pas le coefficient du terme de plus
haut degré.
— si f = Q est un polynôme ( de degré q
0 si q < n
Q[x0 , x1 , ..., xn ] =
coefficient du terme de plus haut degré de Q si q = n.
En effet, si q ≤ n, l’interpolé de Q n’est autre que Q lui-même