M´ethodes num´eriques de base

Le but de ce chapitre est de fournir, de mani`ere succincte, les principes

de base des m´ethodes num´eriques les plus utilis´ees. Mˆeme si d´esormais le
premier contact avec les m´ethodes num´eriques ne se fait plus au travers
d’une programmation de ces m´ethodes dans les langages de bas-niveau, mais
plutoˆt par l’utilisation d’outils tels que Matlab, Scilab, Maple, Mathematica...
permettant la manipulation de concepts math´ematiques ´evolu´es, il reste
indispensable de connaˆıtre les fondements des principales m´ethodes pour
les utiliser de fac¸on pertinente.

1 R´esolution des ´equations du type f(x) = 0

Nous nous restreignons, par souci de simplicit´e, `a la recherche de
racines r´eelles, z´eros de fonctions r´eelles continues.
L’existence et une premi`ere localisation des solutions utilisent le th
´eor`eme des valeurs interm´ediaires.

Th´eor`eme 1.1 (Th´eor`eme des valeurs interm´ediaires) : Si f est une

fonction continue sur [a,b], et si f(a)f(b) ≤ 0, alors il existe au moins un point c ∈
[a,b] tel que f(c) = 0.

Si de plus f est strictement monotone sur [a,b], la racine est unique dans
[a,b].

Page 1 of 44
 14

12

10

ï2

ï4

ï6
ï 1.5 ï1 ï 0.5 0 0.5 1 1.5

Figure III.1 – Fonction pr´esentant 3 racines.

1.1 M´ethode de dichotomie

La m´ethode de dichotomie est un proc´ed´e syst´ematique de raffinement
de la localisation d’une racine. Le mot dichotomie (dicho = deux, tomie =
coupe) exprime clairement le principe de la m´ethode.
Soit [a,b] un intervalle initial de localisation de la racine cherch´ee s.
Supposons que l’on ait f(a)f(b) < 0, l’algorithme de la dichotomie s’´ecrit :

Tant que( abs( b-a ) > epsilon) ! test d’arret m=(a+b)/2

si (f(a) * f(m)) < )
b=m ! s est dans [a,m]
sinon
a=m ! s est dans [m,b]
Fin si
Fin tant que
Cet algorithme r´eduit a` chaque pas l’amplitude de la localisation d’un
facteur 2. L’erreur est donc r´eduite d’un facteur 2 a` chaque it´eration. En 20
it´erations, par exemple l’erreur sera 10 −6 fois l’erreur initiale. Cette m´ethode
est relativement lente. Par contre elle converge dans tous les cas ou` la
fonction change de signe au voisinage de la racine (ce qui exclut les racines de

Page 2 of 44
multiplicit´es paires). C’est une m´ethode que nous qualifierons de m´ethode
tout-terrain, lente mais quasiment infaillible.

1.2 M´ethodes de point-fixe

Les m´ethodes de point-fixe permettent de construire des algorithmes
plus rapides que la dichotomie (parfois) mais surtout des algorithmes qui se
g´en´eralisent simplement au cas de probl`emes en dimension sup´erieure `a
un. On ram`ene l’´equation f(x) = 0 a` une ´equation ´equivalente de forme
point-fixe
x = g(x)
Ceci nous permettra d’obtenir simplement une m´ethode it´erative de la

forme x0 donn´e : initialisation

(III.1)
xn+1 = g(xn)

Si cette it´eration converge, elle converge vers le point-fixe de g, donc de

Figure III.2 – Forme f(x) = 0 et forme point-fixe ´equivalente d’une

´equation.

mani`ere ´equivalente vers le z´ero recherch´e de f. La condition de

convergence essentielle est une condition de contraction sur la fonction g.

Page 3 of 44
D´ef. 1.1 (Application contractante) : On dit qu’une application d´efinie de
[a,b] dans [a,b] est contractante, ou que c’est une contraction, s’il existe un
nombre 0 ≤ k < 1 tel que, pour tout couple de points distincts (x1,x2) de [a,b], on
ait :
|g(x1) − g(x2)| ≤ k|x1 − x2|

k est le facteur de contraction. Il donne la vitesse de convergence de l’it´eration.

Dans le cas ou` g est d´erivable, la condition de contraction se ram`ene `a la

condition suivante sur la d´eriv´ee : |g(x)| ≤ k < 1 ∀x ∈ [a,b]

Remarque 1.1 : Les notions de contraction et le th´eor`eme de convergence des

it´erations associ´e peuvent s’´ecrire et se d´emontrer dans le cadre g´en´eral
des espaces vectoriels norm´es (espaces de Banach). Cette possibilit´e de g´en
´eralisation tr`es large est un des int´erˆets principaux des m´ethodes de point-
fixe (voir la d´emonstration g´en´erale du th´eor`eme ??).

1.3 Vitesse de convergence et ordre d’une m´ethode it

´erative
Nous nous pla¸cons dans le cas d’une it´eration xn+1 = g(xn) convergente et
nous supposons la fonction g suffisamment d´erivable. L’application de la
formule de Taylor au voisinage de la racine s donne :

M´ethodes d’ordre un
Si g(s) = 0 , la limite du rapport des ´ecarts est :

L’´ecart au pas n+1 est donc du mˆeme ordre que l’´ecart au pas n. Le facteur
de r´eduction d’erreur est asymptotiquement donn´e par |g(s)|. Plus petite
sera la valeur de |g(s)|, plus vite se fera la convergence.

Page 4 of 44
M´ethodes d’ordre deux
Si g(s) = 0, l’erreur au pas n + 1 est un infiniment petit d’ordre ≥ 2 par rapport
`a l’erreur au pas n. On obtient en effet :

La convergence est dite quadratique. La r´eduction du nombre des it´erations

est spectaculaire d`es l’ordre 2. A partir d’une erreur de 0` .1, on obtient
10−8 en trois it´erations.
On peut essayer, pour chaque probl`eme particulier, de construire une it
´eration de point-fixe convergente. Il est ´evidemment plus int´eressant
d’utiliser des m´ethodes g´en´erales applicables pour toute ´equation f(x) = 0.
Voici une famille de m´ethodes classiques tr`es utiles dans la pratique.

1.4 M´ethode de Newton et Quasi-Newton

On obtient ´evidemment une forme point-fixe ´equivalente `a f(x) = 0 en
consid´erant la famille x = x − λ(x)f(x), avec λ d´efinie et non nulle sur un
intervalle [a,b] contenant la racine s. Parmi tous les choix possibles pour λ, le

choix qui conduit `a la convergence la plus rapide est d´eriv´ee de

f). La m´ethode obtenue ainsi est la m´ethode de Newton. Il est facile de
v´erifier que l’it´eration

(III.2)

est d’ordre deux si elle converge. Evidemment la m´ethode de Newton n’est´

plus efficace si f s’annulle, donc dans le cas de racines multiples. Il existe des r
´esultats de convergence globale de Newton. Ils supposent des hypoth`eses
de monotonie et de concavit´e de signe constant sur f (pas de point
d’inflexion). La m´ethode de Newton est tr`es classique. Elle s’interpr`ete g
´eom´etriquement comme m´ethode de la tangente.
La m´ethode de Newton n´ecessite le calcul des d´eriv´ees f(xn). C’est un
inconv´enient dans la pratique ou` l’on ne dispose pas toujours d’expression
analytique pour la fonction f.

Page 5 of 44
 Figure III.3 – Interpr´etation g´eom´etrique de la m´ethode de Newton.

Une solution simple est fournie par la m´ethode de la s´ecante ou de

fausseposition dans laquelle on remplace le calcul de f(xn) par l’approximation

Ce qui donne

(III.3)
Les proc´edures de r´esolution d’´equations de type f(x) = 0 que l’on
trouve dans les outils g´en´eriques d’algorithmique (Matlab par exemple),
combinent en g´en´eral une premi`ere approche de la racine par quelques pas
de dichotomie, suivis, pour la convergence fine, par une m´ethode rapide afin
d’obtenir une grande pr´ecision en peu d’it´erations. L’algorithme propos´e
par Matlab est duˆ a` Dekker (1969). La m´ethode rapide utilis´ee est une
interpolation quadratique inverse, assurant l’ordre deux de convergence
(comme Newton), sans calcul de la d´eriv´ee de f.
1.5 M´ethode de Newton et Quasi-Newton pour les syst`emes
La m´ethode de Newton se g´en´eralise en dimension sup´erieure a` un. On
peut en effet montrer que le choix du n + 1e it´er´e xn+1 est tel que

Page 6 of 44
 f(xn+1) = O(|xn+1 − xn|2)

En effet un d´eveloppement simple donne :

f(xn+1) = f(xn) + (xn+1 − xn)f(xn) + O(|xn+1 − xn|2) On voit donc

que le choix

assure bien f(xn+1) = O(|xn+1 − xn|2)

On retrouve ainsi l’ordre 2 de la m´ethode de Newton.

Dans le cas d’un syst`eme de N ´equations non-lin´eaires, on peut ´ecrire

F(Xn+1) = F(Xn) + F (Xn)(Xn+1 − Xn) + O(Xn+1 − Xn2)

la mˆeme id´ee conduit au choix suivant pour l’it´eration vectorielle :

Xn+1 = Xn − {F (Xn)}−1F(Xn) (III.4)

ou` F (Xn) d´esigne la matrice jacobienne de coefficients ).

Pratiquement le n+1e it´er´e est calcul´e a` chaque it´eration par r´esolution de
syst`emes lin´eaires

F (Xn)[Xn+1 − Xn] = −F(Xn) (III.5)

La matrice F (Xn) doit ˆetre assembl´ee et recalcul´ee et le syst`eme doit ˆetre r

´esolu a` chaque it´eration. Ceci rend la m´ethode de Newton tr`es couˆteuse
en temps de calcul.
Pour ´eviter ces inconv´enients, on utilise des m´ethodes dites de Quasi-
Newton dans lesquelles sont propos´ees des approximations de l’inverse de la
matrice Jacobienne. Une m´ethode classique et efficace de ce type est la m
´ethode BFGS (Broyden-Fletcher-Goldfarb-Shanno).

Page 7 of 44
1.5.1 Application `a l’optimisation
Dans un contexte d’optimisation, la recherche du minimum d’une
fonctionnelle J peut ˆetre ramen´ee a` la recherche du point qui annule son
gradient ∇J(x), x d´esignant la param´etrisation du probl`eme (voir chapitre
??). On pourra alors utiliser les m´ethodes de type Newton ou Quasi-Newton
(nous renvoyons `a la litt´erature pour les m´ethodes d’optimisation en
dimension un du type de la m´ethode de la section dor´ee qui ne n´ecessite
pas le calcul des d´eriv´ees). Ces m´ethodes demandent le calcul exact ou
approch´e de la matrice des d´eriv´ees partielles secondes ou matrice
Hessienne. La m´ethode BFGS permet de construire directement une
approximation de l’inverse de la
Hessienne H, en d´emarrant de la matrice identit´e (H0 = Id) et en appliquant
l’it´eration suivante :

Hp)T )γp
Hp+1 = Hp +
(III.6)

ou` p indique l’it´eration d’optimisation, δxp = xp − xp−1 la variation du param`etre

et γp = ∇J (xp+1)−∇J (xp). Voir chapitre ?? pour des d´eveloppements sur
l’optimisation.

2 Interpolation
Une collection de valeurs not´ees yi ´etant donn´ees pour un ensemble
d’abscisses xi, pour i = 0 `a n, l’interpolation est le proc´ed´e qui consiste a` d
´eterminer une fonction, d’une famille choisie a priori, prenant les valeurs
donn´ees aux abscisses correspondantes. Le choix de fonctions polynomiales
est le plus classique. Dans ce cas, le polynˆome d’interpolation est le
polynoˆme de degr´e minimal passant par les n + 1 points donn´es. Ce
polynoˆme est unique et il est de degr´e inf´erieur ou ´egal `a n. Si l’on exprime
le polynoˆme recherch´e dans la base canonique sous la forme

P(x) = a0 + a1x + a2x2 + ... + anxn

Page 8 of 44
on doit r´esoudre un syst`eme lin´eaire de n + 1 ´equations `a n + 1 inconnues.
Sous cette forme le calcul est donc couˆteux. La solution de ce syst`eme est
tr`es sensible aux erreurs d’arrondis.
1

0.8

0.6

0.4

0.2

ï 0.2

ï 0.4

ï 0.6

ï 0.8

ï1
ï5 ï4 ï3 ï2 ï1 0 1 2 3 4 5

Figure III.4 – Interpolation polynomiale de degr´e 10 pour une fonction

sinuso¨ıde.

2.1 Polynoˆmes de Lagrange

Une solution simple, ´el´egante et ´economique de ce probl`eme est
fournie par l’utilisation de la base des polynoˆmes de Lagrange.
On consid`ere les n + 1 polynˆomes Li de degr´e ≤ n qui v´erifient, pour tout i et
j compris entre 0 et n, les ´egalit´es :
Li(xi) = 1
(III.7)
Li(xj) = 0

Les polynˆomes Li sont d´etermin´es de fa¸con unique par les n + 1 ´equations

ci-dessus. Il est facile de montrer qu’ils forment une base de l’espace des
polynˆomes de degr´e inf´erieur ou ´egal a` n et qu’ils s’´ecrivent :

(III.8)

Page 9 of 44
Exprim´e dans cette nouvelle base, le polynoˆme d’interpolation s’´ecrit

) (III.9)
La relation ci-dessus, facile a` v´erifier, explique l’int´erˆet de la base de
Lagrange. Les coefficients du polynoˆme d’interpolation cherch´e sont, dans
cette base, tout simplement les valeurs yi donn´ees. Exprim´e autrement, le
changement de base, de la base canonique a` la base de Lagrange, a transform
´e le syst`eme a` r´esoudre en un syst`eme a` matrice identit´e.

2.2 Limites de l’interpolation polynomiale

L’interpolation polynomiale est la base de nombreuses techniques num
´eriques, en particulier les techniques d’int´egration approch´ee. Elle se g´en
´eralise de fa¸con naturelle aux cas de dimension sup´erieure `a un.
Cependant elle a des limites :
– th´eoriques : on n’est pas assur´e de la convergence du polynoˆme
d’inter-polation vers la fonction interpol´ee lorsque l’on fait tendre le
nombre de points d’interpolation (et donc le degr´e du polynoˆme) vers
l’infini (voir le ph´enom`ene de Runge pour la fonction );
– num´eriques : mˆeme dans le cas ou` la convergence th´eorique est as-
sur´ee, les instabilit´es de calcul provenant de l’accumulation des
erreurs d’arrondis, auxquelles le proc´ed´e d’interpolation polynomiale
est particuli`erement sensible, limite l’usage de cette technique d`es
que le nombre de points d’interpolation d´epasse la dizaine;

Page 10 of 44
 2

1.5

0.5

ï 0.5
ï5 ï4 ï3 ï2 ï1 0 1 2 3 4 5

Figure III.5 – Divergence de l’interpolation polynomiale pour la fonction

. Ph´enom`ene de Runge. En pointill´es : la fonction, en traits pleins :
le polynˆome d’interpolation de degr´e 10 construit sur 11 points r
´eguli`erement espac´es.

– pratiques : remarquons que dans de nombreux cas, les valeurs donn

´eesr´esultent d’exp´eriences ou de calculs pr´ealables. Ces valeurs sont
donc approximatives. Le probl`eme r´eel n’est alors plus un probl`eme
d’interpolation, mais plutoˆt un probl`eme de meilleure approximation
pour lequel les m´ethodes de moindres carr´es, pr´esent´ees plus bas,
sont mieux adapt´ees.

2.3 Interpolation par des splines

Pour ´eviter l’inconv´enient, signal´e plus haut, de l’augmentation du degr
´e du polynoˆme et de l’instabilit´e qui en r´esulte, lorsque le nombre de
points est grand, tout en restant dans un proc´ed´e d’interpolation, on
subdivise l’ensemble des points donn´es en plusieurs sous-ensembles. On r
´ealise les interpolations sur ces petits sous-ensembles, ce qui permet de se
limiter `a des polynˆomes de bas degr´e. Les fonctions polynomiales par
morceaux obtenues sont `a la base des ´el´ements finis de Lagrange.
Les interpolations ci-dessus produisent des fonctions globalement
continues mais non continuˆment d´erivables.

Page 11 of 44
 Les splines cubiques d’interpolation sont des fonctions cubiques par mor-

Figure III.6 – Une fonction affine par morceaux.

Figure III.7 – Une fonction polynomiale de degr´e deux par morceaux.

ceaux, globalement C2. On obtient leur expression analytique, segment par

segment, en imposant les conditions suivantes aux points xi d’interpolation
s(xi) = yi donn´e pour i = 0,..n, s et s continues

Les inconnues du probl`eme sont alors les d´eriv´ees secondes Ci de la spline

1

0.9

0.8

0.7

0.6

0.5

0.4

0.3

0.2

0.1

0
ï5 ï4 ï3 ï2 ï1 0 1 2 3 4 5

Figure III.8 – Spline cubique d’interpolation pour .

Page 12 of 44
aux points xi. On suppose la d´eriv´ee seconde de la spline affine par
intervalles. On int`egre deux fois en prenant en compte les conditions de
continuit´e de la d´eriv´ee et les valeurs donn´ees yi aux points xi. On en d
´eduit les expressions suivantes de la spline sur chaque intervalle [xi,xi+1] :

avec hi = xi+1 − xi, et ou` les Ci sont solutions du syst`eme tridiagonal :

3
pour i = 1,..n − 1, compl´et´e, en g´en´eral, par .

Voici par exemple (Figure III.8) la spline cubique d’interpolation de la

fonction sur 10 intervalles. On observe la stabilit´e de cette

interpolation par contraste avec le r´esultat obtenu (Figure III.5) par
interpolation polynomiale.

3 Approximation au sens des moindres carr´es

L’instabilit´e du proc´ed´e d’interpolation polynomiale lorsque le nombre
de points augmente, d’une part, l’incertitude des r´esultats de mesure, d’autre
part, conduisent a` pr´ef´erer a` l’interpolation des m´ethodes
d’approximation. Ainsi il est clair que l’exp´erimentateur qui rel`evera 100
points quasiment align´es sera plus int´eress´e par la droite passant “ au
mieux ” par ces 100 points plutoˆt que par le polynoˆme de degr´e 99 r
´ealisant l’interpolation exacte.
La plus c´el`ebre et la plus utile des m´ethodes d’approximation est la m
´ethode des moindres carr´es. La formalisation de l’id´ee intuitive d’une

Page 13 of 44
droite repr´esentant “au mieux” un nuage de points au sens des moindres
carr´es se fait de la mani`ere suivante.

3.1 Droite des moindres carr´es

Soient N valeurs y1,y2,...yi,...yN donn´ees aux N abscisses x1,x2,...xi,...xN.
Le polynˆome P de degr´e un : P(x) = a0+a1x (repr´esent´e par une droite) qui

Figure III.9 – Droite des moindres carr´es

r´ealise la meilleure approximation au sens des moindres carr´es des valeurs

yi donn´ees aux points xi est celui qui minimise la somme des carr´es des
´ecarts entre les yi et les P(xi), soit

(III.11)

S apparaˆıt comme le carr´e de la norme euclidienne du vecteur de

composantes

ydu vecteur le plus proche du vecteur i − (a0 + a1xi). La minimisation de S

Ys’interpr`ete donc comme la recherche∈ IRN de composantesU, de

composantesyi, dans le

Page 14 of 44
sous-espace de dimension deux engendr´e par les vecteurs toutes ´egales `a 1,
et X, de composantes xi. Comme la norme utilis´ee est la norme euclidienne, le
vecteur le plus proche est le projet´e orthogonal. On obtient ses composantes
a0 et a1 en ´ecrivant les relations d’orthogonalit´e :
N




 (Y − a00U − a11X|U) = ii=1=1[yii − (a00+ a11xii)] 1 = 0i

(III.12)

 (Y − a U − a X|X) = [y − (a + a x )] x = 00 1

Ceci conduit au syst`eme dit des ´equations normales pour a et a , coefficients

de la droite des moindres carr´es (ou de r´egression) cherch´ee.
N N



Ni =1 xii  a0  i yi 

N
iN =
N=1 (III.13)

x x2  a1 xiyi

i=1 i=1 i=1

Page 15 of 44
3.2 G´en´eralisation : polynˆome des moindres carr´es
Il est facile de g´en´eraliser le calcul pr´ec´edent au cas de la recherche du

polynoˆme de degr´eau sens des moindres carr´es des≤ m, avec m << Nyi. Ce

polynoˆme minimise, qui r´ealise la meilleure approximation

(III.14)
On obtient les relations d’orthogonalit´e :

(III.15)
ou` l’on a not´e Xm le vecteur de composantes xmi . Les valeurs des coefficients
ak du polynˆome des moindres carr´es s’en d´eduisent apr`es r´esolution du
syst`eme lin´eaire d´eduit de (III.15).
Remarque 3.1 : Le syst`eme des moindres carr´es ci-dessus est mal conditionn
´e (il est de plus en plus sensible aux erreurs d’arrondis `a mesure que m
augmente). On se limite habituellement `a des polynˆomes de degr´e peu ´elev´e.
La r´esolution pratique des probl`emes de moindres carr´es se fait par des
algorithmes sp´ecifiques d’´elimination pour des syst`emes rectangulaires surd
´etermin´es. Ces algorithmes utilisent la factorisation QR (technique de
Householder).
On peut envisager, comme dans le cas de l’interpolation, la recherche
d’une meilleure approximation par des splines en d´ecoupant l’ensemble des
points en sous-ensembles plus petits. Cette id´ee de meilleure approximation
au sens des moindres carr´es par des splines est `a la base de nombreuses
techniques d’approximation et de repr´esentation de donn´ees.
4 Int´egration num´erique
La plupart des formules d’int´egration num´erique proviennent de m
´ethodes d’interpolation polynomiales. Le calcul de l’int´egrale d’une fonction
est approch´e par le calcul exact de l’int´egrale d’un polynoˆme interpolant

Page 16 of 44
cette fonction en certains points xk pour k = 1,..p. On obtient ainsi une forme g
´en ´erale pour les quadratures num´eriques
b
p

f(x)dx
= Akf(xk) (III.16)

a
k=1

Les xk sont alors les points d’int´egration et les coefficients Ak les poids de la
formule de quadrature.
Le couˆt d’une formule est mesur´e par le nombre de calculs de f n´ecessaires,
donc par le nombre de points d’int´egration. Le crit`ere d’efficacit´e choisi est
le degr´e maximal de l’espace des polynoˆmes pour lequel la formule est
exacte. La pr´ecision d’une formule d’int´egration num´erique est mesur´ee
par son ordre.
D´ef. 4.1 : On dit qu’une formule est d’ordre k si k est le plus grand entier pour
lequel elle donne l’int´egrale exacte de tout polynˆome de degr´e ≤ k. On montre
que l’erreur d’int´egration est alors, pour toute fonction suffisamment r
´eguli`ere, un infiniment petit d’ordre k du pas d’int´egration h.

4.1 Formules des rectangles

Ce sont les formules a` un point d’int´egration qui proviennent
d’interpolation par des polynoˆmes constants.

) (III.17)
Le couˆt d’une telle formule `a un point est celui d’une ´evaluation de f. Les
choix classiques sont α = a, α = b et le choix donnant la meilleure formule dite
formule du point-milieu (exacte pour les fonctions affines) est .

4.2 Formule des trap`ezes

On consid`ere cette fois l’interpolation par un polynoˆme de degr´e un
construit sur les points a et b. On obtient la formule classique des trap`ezes :

Page 17 of 44
 )] (III.18) Cette formule, `a deux points,
est clairement exacte pour les polynˆomes de degr´e ≤ 1.

4.3 Formule de Simpson

En utilisant l’interpolation sur les trois points a, , on obtient la
formule de Simpson, que l’on v´erifie ˆetre exacte pour les polynˆomes de
degr´e

)] (III.19)
Cette formule `a trois points n´ecessite donc trois ´evaluations de f par
segment (voir cependant, paragraphe 4.5 dans le cas d’un segment global
subdivis´e en sous-intervalles, les formules composites de calcul d’int´egrales
par les m´ethodes des trap`ezes et de Simpson).

4.4 Formules de Gauss

Dans les formules pr´ec´edentes le choix des points d’int´egration ´etait fix
´e (extr´emit´es et/ou milieu des intervalles d’int´egration). Dans les formules
de type Gauss, les points d’int´egration sont choisis de mani`ere a` obtenir la
pr´ecision la plus ´elev´ee.
La Formule de Gauss Legendre a` 2 points, est exacte pour les polynˆomes
de degr´e ≤ 3 :

. (III.20)

avec et

Page 18 of 44
4.5 Formules composites, maillages et m´ethodes
adaptatives
Toutes les formules pr´esent´ees ci-dessus sont des formules de base,
utilisables sur de petits ´el´ements en dimension un, deux ou trois. Pour faire
un calcul r´eel, il faut pr´ealablement d´ecouper le domaine d’int´egration
global en un ensemble de petits sous-domaines ´el´ementaires. Voici, pour
fixer les

Figure III.10 – Choix optimal de points d’int´egration par Matlab pour la

fonction .

id´ees, le calcul global, de l’int´egrale d’une fonction f, par la m´ethodes des

trap`ezes, sur un intervalle [a,b] d´ecoup´e uniform´ement en N sous-
intervalles de longueur h (le pas d’int´egration) :

et le mˆeme calcul par Simpson

ou` P = N/2.

Page 19 of 44
Cependant, un d´ecoupage g´eom´etrique uniforme en sous-intervalles ´egaux
n’est pas optimal. L’op´eration de discr´etisation g´eom´etrique ou “maillage”
que l’on retrouvera dans le contexte des m´ethodes de diff´erences, d’´el
´ements ou de volumes finis est d´eterminante pour la pr´ecision du r´esultat.
Il faut mettre plus de points d’int´egration l`a ou` la fonction pr´esente des
variations relatives rapides. Dans les outils g´en´eriques comme Matlab, on
trouve des proc´edures de quadrature adaptatives qui choisissent
automatiquement les points d’int´egration (ou la taille des mailles) selon les
variations de la fonction a` int´egrer. Typiquement, ces m´ethodes utilisent
des indicateurs d’erreurs locaux bas´es sur l’´ecart entre deux calculs effectu
´es avec deux m´ethodes de base diff´erentes, trap`ezes et Simpson par
exemple (voir Fig III.10).

Page 20 of 44
5 R´esolution des ´equations diff´erentielles
Nous limiterons l’expos´e au cas simple d’´equations diff´erentielles
ordinaires de la forme

 yTrouver la fonction(x) = f(x,y(x)) y :∀x ∈x [→a,by](x) telle que :

(III.21)

y(a) = y0
Les probl`emes diff´erentiels de ce type sont appel´es probl`emes de Cauchy
ou probl`emes a` valeurs initiales. Si f est continue et si elle v´erifie une
condition de Lipschitz par rapport a` la deuxi`eme variable (Il existe L > 0 tel
que ∀x ∈ [a,b] et ∀y1, y2 on ait : |f(x,yy pour toute valeur initiale. On dit qu’il1) −
f(x,y2)| ≤ L|y1 − y2| ), alors le
probl`eme admet une solution unique est “bien pos´e”. On a vu, lors du
chapitre pr´ec´edent, que certains probl`emes diff´erentiels bien pos´es du
point de vue th´eorique pouvaient s’av´erer impossibles a` r´esoudre num
´eriquement car instables. Num´eriquement, il faudra en effet consid´erer
l’ensemble des solutions voisines de la solution exacte cherch´ee, solutions
voisines correspondant a` de petites perturbations des conditions initiales. Si
ces solutions ne s’´ecartent pas trop de la solution de r´ef´erence exacte, on
aura un probl`eme stable et on pourra construire des approximations num
´eriques convenables.
En g´en´eralisant l’´ecriture de l’´equation (III.21) au cas d’une fonction
inconnue vectorielle, on pourra traiter des syst`emes diff´erentiels et des
´equations ou syst`emes d’ordre sup´erieur a` un par des extensions
naturelles des techniques que nous pr´esentons ci-dessous dans le cas de la
dimension un par souci de simplicit´e.

Page 21 of 44
5.1 Principe g´en´eral des m´ethodes num´eriques
La solution exacte d’une ´equation diff´erentielle du type (III.21) est une
fonction continue. Les ordinateurs ne peuvent fournir qu’un nombre fini de r
´esultats num´eriques. Tout commence donc par un choix pr´ealable d’un
nombre fini de points xi sur [a,b]. Ceci s’appelle une discr´etisation ou un
maillage du domaine g´eom´etrique (ici le segment [a,b]). On limitera le calcul
au calcul approch´e de la solution en ces points.
Le choix des points xi est ´evidemment crucial pour la qualit´e de la solution
num´erique obtenue. Le maillage doit permettre de repr´esenter de fa¸con pr
´ecise la solution. Comme cette solution est au d´epart inconnue, on proc`ede
par des techniques d’adaptation de maillage a posteriori. On calcule une
premi`ere solution sur un premier maillage. On d´eduit de ce premier calcul
les zones de forte variation de la solution. On raffine le maillage dans ces
zones.
Encore une fois dans un souci de simplicit´e, nous pr´esenterons ici les m
´ethodes num´eriques dans le cas de maillage a` pas uniformes. Le probl`eme
de l’adaptation de maillage sera trait´e dans un cadre g´en´eral chapitre ??.
On peut construire des sch´emas d’int´egration d’´equations diff´erentielles
de diverses mani`eres. Par exemple :
– les sch´emas d’int´egration `a un pas peuvent ˆetre obtenus en
utilisantdes formules d’int´egration num´erique approch´ee. En
introduisant des sous-pas, on obtient les sch´emas de Runge et Kutta
qui sont les plus pratiques et les plus utilis´es;
– les sch´emas multi-pas peuvent ˆetre construits par d´eveloppement
deTaylor.
Deux crit`eres principaux gouvernent le choix d’une m´ethode :
– l’ordre, qui mesure la pr´ecision du sch´ema ou erreur de troncature
faiteen remplac¸ant l’´equation diff´erentielle exacte par son
approximation. L’ordre de la m´ethode donnera, a` convergence, l’ordre
de l’erreur commise en fonction du pas de discr´etisation.
– la stabilit´e, qui concerne le comportement de la solution approch
´eediscr`ete et la propagation des erreurs d’arrondis dans le cas d’un
calcul r´eel pour un pas de discr´etisation fix´e. Le sch´ema est stable si
la solution discr`ete reste born´ee quel que soit le nombre de pas de
discr´etisation.

Page 22 of 44
Remarque 5.1 : Un crit`ere suppl´ementaire et important de choix des sch´emas
concerne la facilit´e de mise en œuvre, notamment lors d’un red´emarrage des
calculs. Imaginons un calcul instationnaire ne pouvant faire l’objet d’un calcul
complet, soit en raison de la mod´elisation (comme en m´et´eorologie par exemple,
ou` de nouvelles conditions initiales et aux limites doivent ˆetre assimil´ees chaque
jour par le mod`ele), soit en raison de l’impl´ementation et de la dur´ee du calcul.
Par exemple, sur un calculateur parall`ele, le temps d’attente est li´e au temps de
calcul et `a la quantit´e de m´emoire requise. Dans le premier cas, comme dans le
second, on aura recours aux approches `a un pas de type Runge et Kutta pour
assurer une pr´ecision ´elev´ee, plutˆot qu’aux sch´emas multi-pas. En effet, quelle
que soit la pr´ecision demand´ee, les m´ethodes de Runge et Kutta ne n´ecessitent
que le stockage d’un seul ´etat pour un red´emarrage des calculs.
5.2 M´ethodes `a un pas
Dans ces m´ethodes la valeur approch´ee yn+1 de la fonction inconnue y
pour l’abscisse xn+1 est calcul´ee en fonction des seules valeurs de l’abscisse pr
´ec´edente xn, de l’approximation yn et du pas de discr´etisation h. Si yn+1
s’obtient par une formule explicite de la forme

yn+1 = yn + Φ(xn,yn,h)

on dit que la m´ethode est explicite.

Si par contre yn+1 est donn´ee par une relation de la forme g´en´erale

yn+1 = yn + Φ(xn,yn,yn+1,h)
on ne pourra l’obtenir que par la r´esolution d’une ´equation. On dit que la m
´ethode est implicite.
La fonction Φ d´efinit la m´ethode utilis´ee.
Ces sch´emas sont obtenus, par exemple, en int´egrant l’´equation diff
´erentielle et en utilisant des formules d’int´egration num´erique pour le
second membre. L’ordre du sch´ema sera ´egal au degr´e du polynˆome pour
lequel l’int´egration est exacte + 1.

A titre d’exemple, on obtient les sch´emas suivants :`

Page 23 of 44
 – A l’ordre 1, avec une formule exacte pour les polynoˆmes constants
parmorceaux, on obtient le sch´ema d’Euler explicite :

yn+1 − yn = hf(xn,yn)

– Toujours a` l’ordre 1, mais en utilisant le point d’arriv´ee, on obtient

lesch´ema d’Euler implicite :

yn+1 − yn = hf(xn+1,yn+1)

– A l’ordre 2, en utilisant la m´ethode des trap`ezes, on obtient le sch

´ema des trap`ezes ou de Crank-Nicolson :

En g´en´eral un sch´ema explicite a une complexit´e plus faible mais impose

des conditions de stabilit´e plus restrictives (voir plus loin).
y ✻
yn +1
yn ✶
✯ y(x )
❃

✒
y0
✲
x x x n +1 x
0 n

Figure III.11 – Interpr´etation de la m´ethode d’Euler comme m´ethode de la

tangente. En chaque point xn, utiliser un d´eveloppement de Taylor a` l’ordre
un revient `a remplacer la courbe solution par sa tangente.

Page 24 of 44
5.3 Interpr´etations de la m´ethode d’Euler explicite
La m´ethode d’Euler est le prototype le plus simple des m´ethodes num
´eriques de r ´esolution des ´equations diff´erentielles. Pour l’´equation
y(x) = f(x,y(x))
∀x ∈ [a,b] y(a)
= y0
elle s’´ecrit :
y0 donn
´e
(III.22)
yn+1 = yn + hf(xn,yn)
C’est la plus simple des m´ethodes explicites a` un pas.
1. La m´ethode d’Euler provient du d´eveloppement de Taylor d’ordre
unde y au voisinage de xn. On peut montrer que, lorsque cette m´ethode
converge, l’erreur est un infiniment petit d’ordre un en h.
2. G´eom´etriquement, elle revient a` remplacer localement en chaque
pointxn, la courbe solution par sa tangente. On voit donc qu’au cours du
processus num´erique, on va passer, `a chaque pas d’une courbe
solution a` une courbe voisine correspondant a` une condition initiale l
´eg`erement diff´erente. Si le probl`eme est stable, on pourra obtenir la
convergence. Par contre, les r´esultats peuvent vite devenir
catastrophiques dans un cas instable (exemple ?? du chapitre 1).
3. Enfin, en utilisant l’´equivalence entre un probl`eme diff´erentiel et une
´equation int´egrale, on peut interpr´eter, on l’a vu plus haut, la m
´ethode d’Euler comme le r´esultat de l’application de la formule des
rectangles bas´ee sur le point xn

5.4 M´ethodes de Runge et Kutta

Dans les m´ethodes de r´esolution des probl`emes a` valeurs initiales, le
processus de calcul est un processus fini. On avance de N pas, `a partir du

Page 25 of 44
temps initial jusqu’au temps final et on s’arrˆete. Chaque valeur est donc
calcul´ee une fois pour toutes. Sauf technique plus complexe d’adaptation de
maillages, il n’y a pas de r´eit´eration pour am´eliorer le r´esultat. Il faudra
donc utiliser des m´ethodes suffisamment pr´ecises. Ceci explique le recours
a` des m´ethodes d’ordre ´elev´e. Les m´ethodes de Runge et Kutta sont les g
´en´eralisations de la m´ethode d’Euler `a des ordres sup´erieurs a` un. Elles
s’obtiennent a` partir de formules d’int´egration num´eriques plus pr´ecises
que la formule des rectangles.
Consid´erons tout d’abord l’utilisation de la formule des trap`ezes. Elle
conduit a` la m´ethode

y0 donn´e
h
(III.23)
yn+1 = yn + [f(xn,yn) + f(xn+1,yn+1)] 2
Cette m´ethode est une m´ethode implicite. Le calcul de la nouvelle valeur
yn+1 n´ecessite la r´esolution d’une ´equation. Si l’on veut obtenir une m´ethode
explicite du mˆeme ordre, on peut proc´eder de la mani`ere
 suivante : y0 donn´e



hf (III.24)

Ceci peut s’interpr´eter comme une it´eration de point-fixe (limit´ee ici `a un

pas) pour r´esoudre l’´equation du sch´ema implicite des trap`ezes (voir plus
bas ??). On obtient ainsi la m´ethode de Runge et Kutta d’ordre 2 : RK2.
De mˆeme l’utilisation de la formule d’int´egration de Simpson est a` la base
de la formule de Runge et Kutta d’ordre 4 : RK4. C’est l’une des formules les
plus utilis´ees, elle s’´ecrit :

 y0 donn´e, puis pour n ≥ 0

Page 26 of 44
 (III.25)
Pour r´eduire la complexit´e en stockage de ce sch´ema (pas de stockage des
coefficients ki), on peut utiliser le sch´ema de Runge-Kutta sans stockage
suivant : yn+1 = yn + θphf(yn+1),p = 1...q, θp ∈]0,1]

ou` l’on passe de n `a n + 1 apr`es q sous-it´erations, sans stocker les valeurs

interm´ediaires. Les coefficients θp doivent ˆetre cal´es pour r´ealiser une int
´egration exacte `a un ordre donn´e pour une fonction f donn´ee (ce qui est
une limitation de cette technique).
Consid´erons l’´equation y(x) = λy(x), et prenons q = 2, le sch´ema s’´ecrit alors
(on introduit pour la compr´ehension les vi interm´ediaires) :

 v0 = yn

 v1 = v0 + hθ1λv0

v2 = v0 + hθ2λv1 yn+1 = v2
on a donc :
yn+1 = yn + hθ2λyn + h2θ1θ2λ2yn
Or y(x) = λ2y(x) d’ou` par identification : .

Page 27 of 44
5.5 Application aux syst`emes diff´erentiels
300

250

200

150

100

50
0 2 4 6 8 10 12 14 16 18 20

600

500

400

300

200

100

0
0 2 4 6 8 10 12 14 16 18 20

Figure III.12 – Syst`eme proie-pr´edateur : les proies sont repr´esent´ees en

trait continu, les pr´edateurs en pointill´es. En haut le cas b = 0.01 : proies en
milieu fini, en bas le cas b = 0 : pas de limite au d´eveloppement de la
population des proies, hormis la pr´esence de pr´edateurs.

On peut g´en´eraliser simplement l’application des m´ethodes de Runge et

Kutta aux syst`emes. Pour un syst`eme de deux ´equations coupl´ees de la
forme

donn´es


Page 28 of 44
 )= (x,y(x),z(x)) (III.26)
z (x) = g(x,y(x),z(x)) on obtient l’algorithme suivant pour Runge

et Kutta d’ordre 4 :

 y0,z0 donn´es, puis pour n ≥ 0

Voici trois applications classiques

1. Le syst`eme proies-pr´edateurs :

y1 = (a − by1 − cy2 )y1



y21 = (−α + γy21 )y2

Page 29 of 44
 (III.28) y (0) = 300 y (0) = 150 a = 2,b = 0.01 (ou b = 0), c =
0.01,α = 1,γ = 0.01

dont voici le programme en langage Matlab utilisant le sch´ema RK2.

x=0:0.1:20;
h=0.1; a= 2;
b=0.01; c=0.01;
alpha=1;
gamma=0.01;
y(1)=300; z(1)=
100;

for i=1:200
k1=h*(a -b*y(i)-c*z(i))*y(i); l1=h*(-alpha
+gamma*y(i))*z(i); k2=h*(a-b*(y(i)+k1) -c*(z(i)+l1))*(y(i)
+k1); l2=h*(-alpha +gamma*(y(i)+k1))*(z(i)+l1);
y(i+1)=y(i)+(k1+k2)/2; z(i+1)=z(i)+(l1+l2)/2;
end hold on
plot(x,y)
plot(x,z,’--’) hold
off

2. L’´equation du pendule amorti :

 (

2
θ
t) + αθ(t) + k sin(θ(t)) = 0 sur [0,4π]

3. Le syst`eme dynamique de Lorentz : (III.29)

x˙ = σ(y − x), y˙ = ρx − y − xz, z˙ = xy − βz (III.30)

Page 30 of 44
dont les conditions initiales sont pr´ecis´ees dans le programme ci-
dessous :

! Etude du systeme dynamique de Lorentz

x1=x0=-10.; y1=y0=20.;
z1=z0=-5.;
epsilon=1.e-2 ! perturbation 1 pourcent
sigma=10.*(1.+epsilon) ro=28.*(1.+epsilon)
beta=2.6667*(1.+epsilon)

Page 31 of 44
Figure III.13 – Oscillation du pendule amorti : les angles sont repr´esent´es en
pointill´es, leurs d´eriv´ees en trait continu. En bas le diagramme de phase,
angles / vitesses angulaires

irkmax=4 ! Runge et Kutta a 4 pas dt=1.e-3 ! pas de

temps
do kt=1,10000 ! boucle en temps
do irk=1,irkmax ! boucle Runge et Kutta sans stockage
alpha=1./(irkmax+1-irk);
x1=x0+dt*alpha*(sigma*(y1-x1));
y1=y0+dt*alpha*(ro*x1-y1-x1*z1);
z1=z0+dt*alpha*(x1*y1-beta*z1);
enddo
x0=x1; y0=y1; z0=z1; enddo

Ci-dessus le pas de temps a´et´e fix´e a priori, a` la suite d’essais num

´eriques. Nous pouvons, par une analyse de stabilit´e, proposer un
crit`ere pour son choix. Cependant, l’analyse de stabilit´e s’av`ere
souvent difficile pour les syst`emes coupl´es. Dans ce cas, on effectue
l’analyse pour chaque ´equation, en laissant invariantes les autres
variables. Le pas de temps sera alors le minimum des pas de temps
produits par ces analyses. Par exemple, dans le cas du syst`eme de
Lorentz on obtient :

5.6 M´ethodes `a pas multiples

Les m´ethodes de Runge et Kutta sont des m´ethodes `a un pas. Pour
obtenir des m´ethodes d’ordre de pr´ecision ´elev´e, on peut aussi augmenter
le nombre de pas. Dans ce cas, la solution au point n+1 sera fonction des
solutions aux p pas pr´ec´edents avec p > 1. La construction de ces sch´emas
peut se faire en utilisant la d´erivation num´erique et la pr´ecision du sch´ema
sera donn´ee par la pr´ecision de cette d´erivation. Ainsi, pour un sch´ema `a
l’ordre 2, on peut combiner les expressions ci-dessous :

Page 32 of 44
pour aboutir au sch´ema “saute-mouton” (leap-frog en anglais) :

) (III.31)
Un autre sch´ema a` deux pas largement utilis´e est le sch´ema implicite
d’ordre deux de Gear :

) (III.32)
On trouvera dans les ouvrages sp´ecialis´es un grand nombre de formules
multipas d’ordre ´elev´e (voir cependant la remarque 5.1 sur les avantages
pratiques des m´ethodes a` un pas de type Runge-Kutta). Remarquons
toutefois, qu’un sch´ema d’ordre ´elev´e n’a d’int´erˆet que si l’on est assur´e
de la r´egularit´e suffisante de la solution. Pour une solution discontinue ou
peu r´eguli`ere des m´ethodes d’ordre bas sont mieux adapt´ees. D’autre part,
il y a, en g´en´eral, sauf ´evidemment dans le cas de sch´emas implicites,
antagonisme entre ordre
´elev´e et stabilit´e.

5.7 Stabilit´e
Par opposition a` l’ordre de pr´ecision qui utilise la solution du probl`eme
continu, le concept de stabilit´e est bas´e sur la solution discr`ete. Il rend
compte du comportement r´eel de la solution approch´ee pour une valeur
pratique, donc non nulle, du pas h.
Lors d’un calcul r´eel, les erreurs d’arrondis, in´evitables, s’accumulent. Ceci
est particuli`erement ´evident dans le processus de r´esolution d’une
´equation diff´erentielle ou` l’on progresse pas a` pas `a partir d’une valeur
initiale. Il existe diverses conditions de stabilit´e. Tout d’abord la solution
num´erique doit rester born´ee. Cette exigence minimale de stabilit´e peut se
r´ev´eler insuffisante dans la pratique, la borne obtenue ´etant souvent une
exponentielle de la dur´ee qui donc croˆıt infiniment lorsque celle-ci
augmente.

Page 33 of 44
On introduit alors des crit`eres de stabilit´e plus exigeants afin que la solution
num´erique reproduise le comportement physique de la solution exacte. Par
exemple, pour des probl`emes dissipatifs, on imposera des conditions de
stabilit´e permettant d’obtenir une solution de norme d´ecroissante. Le
concept de A-stabilit´e, sous sa forme la plus simple, est bas´e sur l’analyse du
comportement, selon les valeurs du pas h, des solutions num´eriques de
l’´equation mod`ele

y(t) = −ay(t) avec y(0) donn´e et a r´eel > 0 (III.33)

dont la solution exacte est y(t) = y(0)e−at.

Etudions le cas du sch´ema d’Euler explicite. On obtient´ yn+1 = yn − ahyn,

soit
yn = (1 − ah)ny0
at
Si l’expression exacte e− est toujours d´ecroissante en temps et positive, ce
n’est pas le cas de l’expression approch´ee (1 − ah)n, qui selon les valeurs du
pas h > 0, peut tendre vers l’infini, vers z´ero ou prendre alternativement les
valeurs 1 et −1. Pour que la solution approch´ee reproduise le comportement
de la solution exacte, donc reste positive et d´ecroissante, il faut imposer une

condition de stabilit´e sur le pas h. On doit avoir . Pour des syst`emes

diff´erentiels de type X(t) + AX(t) = F (A est suppos´ee sym´etrique d´efinie
positive, donc `a valeurs propres r´eelles positives), la condition de stabilit´e

imposera, pour toute valeur propre λi de la matrice A, les conditions : .

Donc . Si certaines valeurs propres sont grandes, ceci imposera

des pas h tr`es petits, et donc des difficult´es pour le calcul des solutions sur
de longues dur´ees. On dit, dans ce cas, que le syst`eme diff´erentiel est raide
(stiff en anglais).
Etudions maintenant le cas d’un sch´ema implicite, le sch´ema d’Euler impli-´
cite : yn+1 = yn + f(xn+1,yn+1) Son application `a (III.33) donne :

Page 34 of 44
Dans ce cas, quel que soit h > 0, la solution num´erique est born´ee, positive et
d´ecroissante au cours du temps. On dit que le sch´ema est
inconditionnellement stable.

Remarque 5.2 : On retrouve ici les approximations stables et instables de la

fonction exponentielle pr´esent´ees au chapitre pr´ec´edent.
6 M´ethodes de r´esolution des syst`emes lin
´eaires
Les syst`emes lin´eaires sur-d´etermin´es (plus grand nombre
d’´equations que d’inconnues) sont r´esolus, en g´en´eral, en utilisant les
techniques de moindres carr´es. Nous nous limitons donc ici au cas des
syt`emes carr´es comportant autant d’´equations que d’inconnues.
On consid`ere le syst`eme suivant :

qui s’´ecrit matriciellement

AX = B
Soit φ l’application lin´eaire repr´esent´ee par la matrice A, si X repr´esente x
et B, b, AX = B correspond `a φ(x) = b.

6.1 Existence et unicit´e de la solution

Th´eor`eme 6.1 (Syst`emes de Cramer) : Une condition n´ecessaire et
suffisante pour qu’un syst`eme lin´eaire de n ´equations `a n inconnues admette
une solution unique pour tout second membre B ∈ IRn est de mani`ere
´equivalente que
– det(A) = 0
– Ker(φ) = {0}
– rg(φ) = dim(Im(φ)) = n
– les vecteurs lignes ou les vecteurs colonnes de A sont ind´ependants

Page 35 of 44
 – la seule solution de AX = 0 est X = 0

6.2 M´ethodes directes

Les m´ethodes directes de r´esolution des syst`emes lin´eaires sont des m
´ethodes dans lesquelles la solution est obtenue de fa¸con exacte en un nombre
fini d’op´erations. De fa¸con exacte s’entend, sur un ordinateur, aux erreurs
d’arrondis machine pr`es. Le prototype de m´ethode directe est la m´ethode du
pivot de Gauss. Cette m´ethode permet de ramener la r´esolution d’un syst`eme
g´en´eral a` la r´esolution d’un syst`eme triangulaire sup´erieur. La
triangularisation de Gauss consiste a` annuler, par ´etape, colonne par colonne,
les coefficients sous-diagonaux de la matrice par combinaison des lignes avec
une ligne de r´ef´erence ou ligne “pivot”. A la premi`ere ´etape, la ligne pivot` est
la premi`ere ligne, puis la ke a` l’´etape k et ainsi de suite jusqu’a` l’´etape n − 1.
Le coefficient diagonal Ak,k de la ligne de r´ef´erence s’appelle le pivot.
L’´elimination se fait selon

pour k = 1...N −1, i = k+1...N et j = k+1...N, sans oublier de combiner de la

mˆeme fac¸on les composantes du second-membre. La rencontre de pivot nul
peut n´ecessiter la permutation de lignes du syst`eme. Cependant pour
certaines classes de matrices, en particulier les matrices sym´etriques d
´efinies positives, on est assur´e de pouvoir triangulariser le syst`eme par
Gauss sans permutation. La m´ethode du pivot ´equivaut alors a` une
factorisation de type

A = LU

de la matrice A. L est une matrice triangulaire inf´erieure a` diagonale unit´e

et U une matrice triangulaire sup´erieure. Une fois obtenu le syst`eme
triangulaire sup´erieur ´equivalent au syst`eme initial, sa r´esolution se fait
ensuite explicitement par un processus de remont´ee. On commence par
calculer la derni`ere composante du vecteur inconnu en utilisant la derni`ere
´equation et on remonte ´equation par ´equation en d´eterminant les
composantes correspondantes.
Dans le cas d’une matrice A sym´etrique, on peut utiliser la sym´etrie pour
obtenir une factorisation de Crout :

Page 36 of 44
 A = LDLT

avec D diagonale.
Dans le cas d’une matrice A sym´etrique d´efinie positive, la m´ethode de
Choleski conduit a` une factorisation :

A = LLT

On trouve la matrice L, qui cette fois n’est plus `a diagonale unit´e, par un
algorithme d’identification de coefficients.
De

Aij = LikLjk pour j≤i

k=1,i

on d´eduit pour tout i : et pour tout j < i

Les m´ethodes directes pr´esentent l’avantage de fournir la solution exacte

(aux erreurs d’arrondis machine pr`es) en un nombre fini d’op´erations (de
l’ordre de pour Gauss). La m´ethode du pivot de Gauss s’applique `a tout
syst`eme inversible. Par contre les m´ethodes directes ont un couˆt important
en stockage m´emoire (bien que des techniques de stockage minimal associ
´ees a` des algorithmes de num´erotation optimale des inconnues permettent
de le r´eduire sensiblement). Ceci rend leur application pratiquement
impossible, en l’´etat actuel de la technologie, pour de gros syst`emes `a plus

Page 37 of 44
de 105 inconnues, et donc en particulier pour la r´esolution de probl`emes
industriels en dimension 3 d’espace.
Nous renvoyons `a la litt´erature pour plus de pr´ecisions sur les m´ethodes
directes (voir Lascaux-Th´eodor, Ciarlet, Lucquin-Pironneau).

6.3 M´ethode de Gauss-Seidel ou de relaxation

Soit encore
AX = B
le syst`eme lin´eaire `a r´esoudre. La m´ethode de Gauss-Seidel correspond
aussi a` la d´ecomposition de la matrice A sous la forme

A=D−E−F

avec D matrice diagonale constitu´ee des ´el´ements diagonaux aii de A,

−E, matrice triangulaire inf´erieure stricte, constitu´ee des ´el´ements
strictements sous-diagonaux de A : aij pour i > j, et

−F, matrice triangulaire sup´erieure stricte, constitu´ee des ´el´ements

strictements sur-diagonaux de A : aij pour i < j.

Mais on d´efinit cette fois la m´ethode de Gauss-Seidel comme la m´ethode it

´erative :
  X
(0)
donn´e
(III.34)

X(k+1) = (D − E)−1F X(k) + (D − E)−1B

La matrice d’it´eration de Gauss-Seidel est donc

L = (D − E)−1F

et la condition de convergence s’exprime par

ρ(L) < 1

Page 38 of 44
On observe que la m´ethode de Gauss-Seidel correspond a` l’´ecriture ligne
par ligne suivante :

(III.35)
Elle utilise les informations les plus r´ecentes sur les composantes des it´er
´es. La m´ethode de Gauss-Seidel converge en g´en´eral plus vite que la m
´ethode de Jacobi. Par contre, Gauss-Seidel n’est pas adapt´ee a` la
vectorisation.
Les m´ethodes de relaxation sont des extrapolations de la m´ethode de
Gauss-Seidel d´ependant d’un param`etre ω. Ce param`etre est d´etermin´e
pour obtenir une convergence plus rapide. Les m´ethodes de
relaxation s’´ecrivent : donn´e
(III.36)
X = (D − ωE)− (ωF + (1 − ω)D) X + ω(D − ωE)− B
+1) 1 (k) 1

La m´ethode de relaxation s’´ecrit ligne par ligne comme une extrapolation de

la composante obtenue par Gauss-Seidel. On a :
Xik+1(relax) = ω Xik+1(GS) + (1 ω
k
)X (relax)
i

La m´ethode de Gauss-Seidel correspond au cas = 1. Dans le cas de matrices A

sym´etriques d´efinies positives, les m´ethodes de relaxation convergent pour
0 < ω < 2.
6.4 M´ethodes de descente - M´ethode du gradient
Les m´ethodes de descente sont des m´ethodes it´eratives qui utilisent
l’´equivalence entre les probl`emes suivants (voir thEor` Eme¨ ??) :

  Trouver
X ∈ IRN tel que
(III.37)
 AX = B

et

Page 39 of 44
 
 Trouver X ∈ IRN tel que  (III.38)

1
J(X) = 2(AX,X) − (B,X) soit minimal

et qui sont donc limit´ees au cas des syst`emes lin´eaires dont la matrice A est
sym´etrique d´efinie positive. Elles se g´en´eralisent, par contre, au cas nonlin
´eaire de la minimisation de fonctionnelles strictement convexes quelconques
(voir chapitre ??).
Les m´ethodes de descente sont bas´ees sur le calcul de la solution comme
limite d’une suite minimisante de la forme quadratique J. Cette suite est
construite comme une suite r´ecurrente :
  X
(0)
donn´e
(III.39)
 X(k+1) = X(k) − α d (k)
k

aveccoeffidcient d´etermin´e de mani`ere `a minimiser(k) ∈ IRN vecteur

donnant la direction de descente a` l’´etapeJ dans la directionkdet(k)α. ∈
IR k

J(X(k) − αk d(k)) ≤ J(X(k) − αd(k)) ∀α ∈ IR

Apr`es d´eveloppement, on obtient αk comme valeur annulant la d´eriv´ee de J

par rapport a` α, soit :

avec
On peut montrer que la m´ethode de Gauss-Seidel est la m´ethode de
descente correspondant aux choix des axes de coordonn´ees comme
directions successives de descente. Dans le cas de la m´ethode du gradient,
on choisit comme direction de descente la direction du vecteur gradient de J

Page 40 of 44
au point X(k). Cette direction est la direction de variation maximale de J. On dit
encore direction de plus profonde descente, d’ou` le nom :“steepest descent
method” , employ´e dans certains ouvrages en anglais.
L’it´eration de la m´ethode du gradient s’´ecrit :

  X
(0)
donn´e
(III.40)

 X(k+1) = X(k) − αk G(k)

avec

6.4.1 Vitesse de convergence. Conditionnement.

A partir de` G(k+1) = AX(k+1) − B et de X(k+1) = X(k) − αk G(k), on obtient la r

´ecurrence suivante sur les gradients :

G(k+1) = G(k) − αk AG(k)

On en d´eduit G(k+1)2 = G(k)2 − 2αk(AG(k),G(k)) + αk2AG(k)2

avec

on obtient :

Et en utilisant les valeurs propres et les vecteurs propres de la matrice A

suppos´ee d´efinie positive, on d´eduit :

Page 41 of 44
D´ef. 6.1 : Pour une matrice sym´etrique r´eelle le nombre

rapport des plus grandes et plus petites valeurs propres de A est appel´e
nombre de conditionnement de la matrice A.

Plus il est proche de 1, plus vite la m´ethode du gradient converge. Les

techniques de pr´econditionnement ont, entre autres, pour but de rapprocher
les valeurs propres extrˆemes afin d’acc´elerer la convergence des m´ethodes
it´eratives.

Remarque 6.1 : Le nombre de conditionnement est, en particulier, ´egal `a 1

pour la matrice identit´e. D’ailleurs, dans ce cas, l’expression `a minimiser d
´ecrirait, en dimension deux, un ensemble de cercles concentriques. Les
gradients ayant la direction d’un rayon, on obtiendrait la solution en une it
´eration. Cette remarque est `a la base de l’id´ee de la m´ethode du gradient
conjugu´e.

6.5 M´ethode du gradient conjugu´e

Dans le cas d’une matrice sym´etrique d´efinie positive quelconque A, les
iso-J sont des hyper-ellipses. Elles redeviennent des cercles pour le produit
scalaire d´efinie par A :

X,Y −→ (X,Y )A = (AX,Y )

D’ou` l’id´ee de remplacer les directions de descente dans la m´ethode du

gradient, par des directions conjugu´ees, c’est `a dire orthogonales au sens du
produit scalaire (.,.)A.
On obtient alors l’algorithme :

 X0 ∈ IRN donn´e

k = 0,1...

Page 42 of 44
 (III.41)
Nous renvoyons a` la litt´erature pour une ´etude exhaustive de la m´ethode
du gradient conjugu´e. Retenons que la propri´et´e de A-orthogonalit´e
entraˆıne une convergence th´eorique en N it´erations. Ce qui en ferait une m
´ethode directe. Cependant, elle serait alors plus ch`ere que Choleski. De plus,
les erreurs d’arrondis font que les propri´et´es d’orthogonalit´e ne sont pas v
´erifi´ees exactement. Il faut donc consid´erer le gradient conjugu´e comme
une m´ethode it´erative. On montre que sa vitesse de convergence d´epend
´egalement du conditionnement de la matrice A. On est conduit a` pr
´econditionner A pour obtenir des performances int´eressantes. Parmi les pr
´econditionnements classiques, on citera le pr´econditionnement SSOR (O.
Axelsson) et par Choleski incomplet (Meijerink et Van der Vorst).
Dans le cas ou` la matrice du syst`eme n’est pas sym´etrique d´efinie
positive, il est plus difficile d’obtenir des m´ethodes it´eratives performantes.
Mentionnons la m´ethode GMRES (Saad et Schultz), dont il existe ´egalement
une version pour la r´esolution de syst`emes non-lin´eaires.

6.6 Application des m´ethodes de gradient au cas nonlin

´eaire
Dans le cas plus g´en´eral de minimisation d’une fonctionnelle J
strictement convexe non n´ecessairement quadratique, le gradient ∇J est non-
lin´eaire. Cependant les m´ethodes de gradient peuvent s’appliquer (on peut
mˆeme les appliquer sans garantie de convergence, car dans tous les cas on r
´eduira J, voir chapitre ??). La d´etermination du pas optimal αk a` l’it´eration
k se fait alors par une m´ethode de recherche de l’argument αk qui minimise la
fonction J(Xk −α∇J(Xk)). On est ramen´e `a un probl`eme en dimension un pour

Page 43 of 44
lequel diverses techniques existent, en particulier les algorithmes de
recherche de type section dor´ee.

L’extension de la m´ethode du gradient conjugu´e au cas non-lin´eaire n

´ecessite, de plus, la d´efinition des directions de descente “conjugu´ees”.
L’algorithme le plus efficace est celui de Polak-Ribi`ere. Les directions de
descente successives sont donn´ees par
dk+1 = ∇J(Xk+1) + βk+1 dk

avec cette fois

Page 44 of 44