Introduction au calcul scientifique
Introduction au calcul scientifique
On peut définir le CALCUL SCIENTIFIQUE comme la discipline qui permet de reproduire sur un ordinateur un phénomène
ou un processus décrit par un modèle mathématique.
A NALYSE MATHÉMATIQUE
P HÉNOMÈNE PHYSIQUE , ÉCO - M ODÈLE MATHÉMATIQUE Bien posé
NOMIQUE , BIOLOGIQUE . . . Mise en équations : Bien conditionné
Observation expérimentale équations différentielles, inté- Propriétés de la solution (stabi-
Modèle conceptuel grales, stochastiques. . . lité,. . .)
Solutions analytiques
C ALCULS
P ROGRAMMATION
Langage de programmation (C, A NALYSE NUMÉRIQUE
C++, Fortran, Java, Python, Mat- Méthodes de discrétisation
lab, Scilab, Octave. . .) Analyse des algorithmes (rapi-
Structure des données dité, précision, souplesse)
Implémentation de l’algorithme Estimation des erreurs
Optimisation
P OST P ROCESSING
Visualisation C ALCUL SCIENTIFIQUE
Analyse des résultats
L’ordinateur est aujourd’hui un outil incontournable pour simuler et modéliser des systèmes complexes, mais il faut
encore savoir exprimer nos problèmes (physiques, économiques, biologiques. . .) en langage formalisé des mathématiques
pures sous la forme d’équations mathématiques (différentielles, intégrales. . .). Nous sommes habitués à résoudre les pro-
blèmes de façon analytique, alors que l’ordinateur ne travaille que sur des suites de nombres. On verra qu’il existe souvent
plusieurs approches pour résoudre un même problème, ce qui conduit à des algorithmes différents. Un des objectifs de ce
cours est de fournir des bases rigoureuses pour développer quelques algorithmes utiles dans la résolution de problèmes en
mathématique, économie, physique. . .
Un algorithme, pour être utile, doit satisfaire un certain nombre de conditions. Il doit être :
rapide : le nombre d’opérations de calcul pour arriver au résultat escompté doit être aussi réduit que possible ;
précis : l’algorithme doit savoir contenir les effets des erreurs qui sont inhérentes à tout calcul numérique (ces erreurs
peuvent être dues à la modélisation, aux données, à la représentation sur ordinateur ou encore à la troncature) ;
souple : l’algorithme doit être facilement transposable à des problèmes différents.
Le choix et l’optimisation des algorithmes numériques mis en pratique sont absolument cruciaux tant pour les calculs
de type industriel souvent très répétitifs et devant donc pouvoir être exécutés en un temps très court, que pour les cal-
culs de référence pour lesquels la seule limite est la patience de celui qui les fait. Par exemple, en fluidodynamique, en
laissant tourner une station de travail pendant quelques jours, les numériciens résolvent des systèmes frisant le milliard
d’inconnues. L’expérience montre qu’entre une approche numérique standard et une approche soigneusement réfléchie
1
Introduction au calcul scientifique
et optimisée un gain de temps de calcul d’un facteur 100, voire davantage, est souvent observé. Il est clair qu’on peut pas-
ser ainsi, grâce à cet effort, d’un calcul totalement déraisonnable à un calcul parfaitement banal : tout l’enjeu de l’analyse
numériques est là ! C’est dire l’importance pour tous scientifique de bien connaître ces méthodes, leurs avantages et leurs
limites.
p
Exemple Calcul de A
Sur ordinateur, l’addition de deux entiers peut se faire de façon exacte mais non le calcul d’une racine carrée. On procède alors par
approximations successives jusqu’à converger vers la solution souhaitée. Il existe pour cela divers algorithmes. Le suivant est connu
depuis l’antiquité (mais ce n’est pas celui que les ordinateurs utilisent).
Soit A un nombre réel positif dont on cherche la racine carrée. Désignons par x 0 la première estimation de cette racine (généralement
le plus grand entier dont le carré est inférieur à A ; par exemple, si A = 178, alors x 0 = 13 car 132 = 169 < 178 et 142 = 196 > 178) et par
ε0 l’erreur associée : p
A = x 0 + ε0 .
Cherchons une approximation de ε0 . On a
A = (x 0 + ε0 )2 = x 02 + 2x 0 ε0 + ε20 .
Supposons que l’erreur soit petite face à x 0 , ce qui permet de négliger le terme en ε20 :
A ' x 02 + 2x 0 ε0 .
Remplaçons l’erreur ε0 par un ε00 , qui en est une approximation, de telle sorte que
A = x 02 + 2x 0 ε00 .
On en déduit que
A − x 02
ε00 =
2x 0
donc la quantité
1 A
µ ¶
x 1 ≡ x 0 + ε00 =+ x0
2 x0
constitue une meilleure approximation de la racine que x 0 (sous réserve que le développement soit convergent). De plus, rien ne nous
empêche de recommencer les calculs avec x 1 , puis x 2 , etc., jusqu’à ce que la précision de la machine ne permette plus de distinguer
le résultat final de la véritable solution. On peut donc définir une suite, qui à partir d’une estimation initiale x 0 devrait en principe
converger vers la solution recherchée. Cette suite est
1 A
µ ¶
x k+1 = + xk , x 0 > 0.
2 xk
A−x 2
3. Calculer l’erreur associée ε0k+1 = 2x k+1 .
k+1
4. Tant que l’erreur est supérieure à un seuil fixé, recommencer au point 2
Le tableau ci-dessous illustre quelques itérations de cet algorithme pour le cas où A = 5 :
k xk ε0k
0 2.0000000000 0.2360679775
1 2.2500000000 0.0139320225
2 2.2361111111 0.0000431336
3 2.2360679779 0.0000000004
4 2.2360679775 0.0000000000
On voit que l’algorithme converge très rapidement et permet donc d’estimer la racine carrée d’un nombre moyennant un nombre li-
mité d’opérations élémentaires (additions, soustractions, divisions, multiplications). Il reste encore à savoir si cet algorithme converge
toujours et à déterminer la rapidité de sa convergence. L’analyse numérique est une discipline proche des mathématiques appliquées,
qui a pour objectif de répondre à ces questions de façon rigoureuse.
Les erreurs
Le simple fait d’utiliser un ordinateur pour représenter des nombres réels induit des erreurs. Par conséquent, plutôt que
de tenter d’éliminer les erreurs, il vaut mieux chercher à contrôler leur effet. Généralement, on peut identifier plusieurs
niveaux d’erreur dans l’approximation et la résolution d’un problème physique.
2
Au niveau le plus élevé, on trouve l’erreur qui provient du fait qu’on a réduit la réalité physique à un modèle mathéma-
tique. De telles erreurs limitent l’application du modèle mathématique à certaines situations et ne sont pas dans le champ
du contrôle du Calcul Scientifique.
On ne peut généralement pas donner la solution explicite d’un modèle mathématique (qu’il soit exprimé par une in-
tégrale, une équation algébrique ou différentielle, un système linéaire ou non linéaire). La résolution par des algorithmes
numériques entraîne immanquablement l’introduction et la propagation d’erreurs d’arrondi. De plus, il est souvent né-
cessaire d’introduire d’autres erreurs liées au fait qu’un ordinateur ne peut effectuer que de manière approximative des
calculs impliquant un nombre infini d’opérations arithmétiques. Par exemple, le calcul de la somme d’une série ne pourra
être accompli qu’en procédant à une troncature convenable. On doit donc définir un problème numérique, dont la solution
diffère de la solution mathématique exacte d’une erreur, appelée erreur de troncature. La somme des erreurs d’arrondis et
de troncature constitue l’erreur de calcul. L’erreur de calcul absolue est la différence entre x, la solution exacte du modèle
mathématique, et x̃, la solution obtenue à la fin de la résolution numérique, tandis que (si x 6= 0) l’erreur de calcul relative
est définie par l’erreur de calcul absolue divisé par x. Le calcul numérique consiste généralement à approcher le modèle
mathématique en faisant intervenir un paramètre de discrétisation, que nous noterons h et que nous supposerons posi-
tif. Si, quand h tend vers 0, la solution du calcul numérique tend vers celle du modèle mathématique, nous dirons que le
calcul numérique est convergent. Si de plus, l’erreur (absolue ou relative) peut être majorée par une fonction de C h p où
C est indépendante de h et où p est un nombre positif, nous dirons que la méthode est convergente d’ordre p. Quand, en
plus d’un majorant, on dispose d’un minorant C 1 h p (C 1 étant une autre constante (≤ C ) indépendante de h et p), on peut
remplacer le symbole ≤ par '.
3
1. Résolution d’équations non linéaires
Recherche des solutions de l’équation non linéaire f (x) = 0 où f est une fonction donnée
Un des problèmes classiques en mathématiques appliquées est celui de la recherche des valeurs pour lesquelles une
fonction donnée s’annule. Dans certains cas bien particuliers, comme pour les fonctions x 7→ x + 1, x 7→ cos(2x) ou encore
x 7→ x 2 − 2x + 1, le problème est simple car il existe pour ces fonctions des formules qui donnent les zéros explicitement.
Toutefois, pour la plupart des fonctions f : R → R il n’est pas possible de résoudre l’équation f (x) = 0 explicitement et il faut
recourir à des méthodes numériques. Ainsi par exemple une brève étude de la fonction f (x) = cos(x) − x montre qu’elle
possède un zéro à proximité de 0.7 mais ce zéro ne s’exprime pas au moyen de fonctions usuelles et pour en obtenir une
valeur approchée il faut recourir à des méthodes numériques.
Plusieurs méthodes existent et elles différent pas leur vitesse de convergence et par leur robustesse. Lorsqu’il s’agit de
calculer les zéros d’une seule fonction, la vitesse de la méthode utilisée n’est souvent pas cruciale. Cependant, dans certains
applications il est nécessaire de calculer les zéros de plusieurs milliers de fonctions et la vitesse devient alors un élément
stratégique. Par exemple, si on veut représenter graphiquement l’ensemble de points (x, y) du plan pour lequel x 2 sin(y) +
e x+y − 7 = 0, on peut procéder comme suit : pour une valeur donnée de x, on cherche l’ensemble des valeurs de y pour
lesquelles x 2 sin(y)+e x+y −7 = 0, i.e. l’ensemble des zéros de la fonction f (y) = x 2 sin(y)+e x+y −7 = 0, et on représente tous
les couples obtenus. On choisit une nouvelle valeur pour x et on répète l’opération. Pour obtenir une courbe suffisamment
précise, il faut procéder à la recherche des zéros d’un très grand nombre de fonction et il est alors préférable de disposer
d’une méthode rapide.
Soit f : R → R une fonction continue donnée dont on veut évaluer numériquement un ou plusieurs zéros xb, c’est-à-dire
qu’on cherche tous les xb tels que f (b
x ) = 0. Les méthodes numériques pour approcher xb consistent à :
¬ localiser grossièrement le (ou les) zéro(s) de f en procédant à l’étude du graphe de f et/ou à des évaluations qui sont
souvent de type graphique ; on note x 0 cette solution grossière ;
construire, à partir de x 0 , une suite x 1 , x 2 , x 3 , . . . telle que limk→∞ x k = xb où f (b
x ) = 0. On dit alors que la méthode est
convergente.
x k+1 = G(x k ), k = 0, 1, 2, . . .
dans lequel on part d’une valeur donnée x 0 pour calculer x 1 , puis à l’aide de x 1 on calcul x 2 etc. La formule même est
dite formule de récurrence. Le procédé est appelé convergent si x k tend vers un nombre fini lorsque k tend vers +∞. Il
est bien évident qu’une méthode itérative n’est utile que s’il y a convergence vers les valeurs cherchées.
On peut parfaitement envisager des méthodes itératives multi-niveaux, comme par exemples les schémas à trois niveaux
dans lesquels on part de deux valeurs données x 0 et x 1 pour calculer x 2 , puis à l’aide de x 1 et x 2 on calcule x 3 etc.
Définition Ordre de convergence
Soit p un entier positif. On dit qu’une méthode (à deux niveaux) convergente est d’ordre p s’il existe une constante C
telle que
x − x k |p .
x − x k+1 | ≤ C |b
|b
ou encore
x k+1 − xb
lim = C.
k→∞ (x k − xb)p
Si p = 1 (et C < 1) on parle de convergence linéaire, si p = 2 on parle de convergence quadratique.
4
1. Résolution d’équations non linéaires
Pour localiser grossièrement le (ou les) zéro(s) de f on va d’abord étudier de la fonction f , puis on va essayer d’utiliser un
corollaire du théorème des valeurs intermédiaires et le théorème de la bijection afin de trouver un intervalle qui contient
un et un seul zéro.
Théorème des valeurs intermédiaires
Soit f une fonction continue sur un intervalle I = [a; b] de R. Alors f atteint toutes les valeurs intermédiaires entre f (a)
et f (b). Autrement dit :
? si f (a) ≤ f (b) alors pour tout d ∈ [ f (a), f (b)] il existe c ∈ [a; b] tel que f (c) = d ;
? si f (a) ≥ f (b) alors pour tout d ∈ [ f (b), f (a)] il existe c ∈ [a; b] tel que f (c) = d .
Ce théorème garantit juste l’existence d’un zéro. Pour l’unicité on essayera d’appliquer le théorème de la bijection dont
l’énoncé est rappelé ci-dessous.
Théorème de la bijection
Soit f une fonction continue et strictement monotone sur un intervalle I de R, alors f induit une bijection de I dans
f (I ). De plus, sa bijection réciproque est continue sur I , monotone sur I et de même sens de variation que f .
5
1. Résolution d’équations non linéaires
[a k ; b k ] en [a k ; c k ] et [c k ; b k ] où c k est
ak + bk
ck = .
2
? Dans la méthode de L AGRANGE, plutôt que de diviser l’intervalle [a k ; b k ] en deux intervalles de même longueur, on
découpe [a k ; b k ] en [a k ; c k ] et [c k ; b k ] où c k est l’abscisse du point d’intersection de la droite passant par (a k , f (a k ))
et (b k , f (b k )) et l’axe des abscisses, autrement dit c k est solution de l’équation
f (b k ) − f (a k )
(c − b k ) + f (b k ) = 0
bk − ak
qui est
bk − ak a k f (b k ) − b k f (a k )
ck = bk − f (b k ) = .
f (b k ) − f (a k ) f (b k ) − f (a k )
Dans les deux cas, pour l’itération suivante, on pose soit [a k+1 ; b k+1 ] = [a k ; c k ] soit [a k+1 ; b k+1 ] = [c k ; b k ] de sorte à ce que
f (a k+1 )· f (b k+1 ) < 0. La suite (c k )k∈N converge vers xb puisque la longueur de ces intervalles tend vers 0 quand k tend vers
+∞.
Soit ε l’erreur maximale qu’on peut commettre, les algorithmes s’écrivent alors comme suit :
D ICHOTOMIE : L AGRANGE :
Require: a, b > a, ε, f : [a, b] → R Require: a, b > a, ε, f : [a, b] → R
k ←0 k ←0
ak ← a ak ← a
bk ← b bk ← b
ak + bk bk − ak
xk ← xk ← ak − f (a k )
2 f (b k ) − f (a k )
while b k − a k > ε or | f (x k )| > ε do while b k − a k > ε or | f (x k )| > ε do
if f (a k ) f (x k ) < 0 then if f (a k ) f (x k ) < 0 then
a k+1 ← a k a k+1 ← a k
b k+1 ← x k b k+1 ← x k
else else
a k+1 ← x k a k+1 ← x k
b k+1 ← b k b k+1 ← b k
end if end if
a k+1 + b k+1 b k+1 − a k+1
x k+1 ← x k+1 ← a k+1 − f (a k+1 )
2 f (b k+1 ) − f (a k+1 )
k ← k +1 k ← k +1
end while end while
On n’est pas obligé de stoker tous les intervalles et les itérées, on peut gagner de la mémoire en les écrasant à chaque
étape :
D ICHOTOMIE : L AGRANGE :
Require: a, b > a, ε, f : [a, b] → R Require: a, b > a, ε, f : [a, b] → R
a +b b−a
x← x ←a− f (a)
2 f (b) − f (a)
while b − a > ε or | f (x)| > ε do while b − a > ε or | f (x)| > ε do
if f (a) f (x) < 0 then if f (a) f (x) < 0 then
b←x b←x
else else
a←x a←x
end if end if
a +b b−a
x← x ←a− f (a)
2 f (b) − f (a)
end while end while
Remarque
Avec la méthode de la dichotomie, les itérations s’achèvent à la m-ème étape quand |x m − xb| ≤ |I m | < ε, où ε est une
tolérance fixée et |I m | désigne la longueur de l’intervalle I m . Clairement I k = b−a
2k
, donc pour avoir une erreur |x m − xb| < ε,
6
1. Résolution d’équations non linéaires
Notons que cette inégalité est générale : elle ne dépend pas du choix de la fonction f .
1+T ¡
(1 + T )n − 1 − M .
¢
f (T ) = v
T
Étudions la fonction f :
? f (T ) > 0 pour tout T > 0,
? limT →0+ f (T ) = nv − M < (n − 1)v, limT →+∞ f (T ) = +∞,
? f 0 (T ) = v2 1 + (1 + T )n (T n − 1) > 0 pour tout T > 0 (comparer le graphe de −1/(1 + T )n et de nT − 1)
¡ ¢
T
En étudiant la fonction f on voit que, comme nv < M dès que n > 1, elle admet un unique zéro dans l’intervalle ]0, +∞[ (on peut
même prouver qu’elle admet un unique zéro dans l’intervalle ]0, M [).
Supposons que v = 1000 ¤ et qu’après 5 ans M est égal à 6000 ¤. En étudiant la fonction f on voit qu’elle admet un unique zéro dans
l’intervalle ]0.01, 0.1[. Si on applique la méthode de la dichotomie avec ε = 10−12 , après 37 itérations de la méthode de dichotomie, la
méthode converge vers 0.061402411536. On conclut ainsi que le taux d’intérêt T est approximativement égal à 6.14%. Ce résultat a été
obtenu en faisant appel à la fonction dichotomie définie à la page 24 comme suit :
Exemple
Soit f (x) = −39 − 43x + 39x 2 − 5x 3 . On cherche a estimer x ∈ [1; 5] tel que f (x) = 0.
y D ICHOTOMIE
f (3) = 1
f (2.5) = 0.3984375
1 2
f (2) = −0.1875 2.5 3 5 x
f (1) = −1
I 0 = [1; 5]
I 1 = [1; 3]
I 2 = [2; 3]
I 3 = [2; 2.5]
7
1. Résolution d’équations non linéaires
y L AGRANGE
f (1) = −1
I 0 = [1; 5]
I 1 = [1; 2.3̄]
I 2 = [2.11; 2.3̄]
La méthode de dichotomie est simple mais elle ne garantit pas une réduction monotone de l’erreur d’une itération à
l’autre : tout ce dont on est assuré, c’est que la longueur de l’intervalle de recherche est divisée par deux à chaque étape. Par
conséquent, si le seul critère d’arrêt est le contrôle de la longueur de I k , on risque de rejeter de bonnes approximations de
xb. En fait, cette méthode ne prend pas suffisamment en compte le comportement réel de f . Il est par exemple frappant que
la méthode ne converge pas en une seule itération quand f est linéaire (à moins que le zéro xb ne soit le milieu de l’intervalle
de recherche initial).
ce qui donne
x k − x k−1
x = xk − f (x k ).
f (x k ) − f (x k−1 )
8
1. Résolution d’équations non linéaires
On peut se demander comment exploiter cette procédure pour calculer les zéros d’une fonction donnée. Remarquons
qu’on peut voir ` comme un point fixe du cosinus, ou encore comme un zéro de la fonction f (x) = x − cos(x). La méthode
proposée fournit donc un moyen de calculer les zéros de f . Précisons ce principe : soit f : [a, b] → R la fonction dont on
cherche le zéro. Il est toujours possible de transformer le problème
(Pb-1) “chercher x tel que f (x) = 0”
en un problème équivalent (i.e. admettant les mêmes solutions)
(Pb-2) “chercher x tel que x − ϕ(x) = 0”.
Pour que les deux problèmes soient équivalent, la fonction auxiliaire ϕ : [a, b] → R doit être choisie de manière à ce que
ϕ(bx ) = xb si et seulement si f (b
x ) = 0 dans [a; b] (on dit alors que le problème (Pb-2) est consistant avec le problème (Pb-1)).
Clairement, il existe une infinité de manières pour opérer cette transformation. Par exemple, on peut poser ϕ(x) = x −
f (x) ou plus généralement ϕ(x) = x + γ f (x) avec γ ∈ R∗ quelconque. On peut même remplacer γ par une fonction de x
pour autant qu’elle ne s’annule pas.
Naturellement une telle suite n’est pas forcement convergente. Par contre, si elle converge, c’est-à-dire si la suite x k a une
limite que nous notons `, et si ϕ est continue, alors cette limite est nécessairement un point fixe de ϕ puisque
µ ¶
` = lim x k+1 = lim ϕ(x k ) = ϕ lim x k = ϕ(`).
k→∞ k→∞ k→∞
On utilise alors l’algorithme itératif suivant pour construire la suite (comme l’ordinateur ne peux pas construire une infinité
de termes, on calcul les premiers termes de la suite et on s’arrête dès que la différence entre deux éléments de la suite est
inférieure à une tolérance ε > 0 donnée) :
Require: x 0 , ε, ϕ : [a, b] → R
k ←0
x 1 ← x 0 + 2ε
while |x k+1 − x k | > ε do
x k+1 ← ϕ(x k )
k ← k +1
end while
On va maintenant s’intéresser à la convergence de la suite construite par une méthode de point fixe.
Théorème Convergence (globale) des itérations de point fixe
Considérons une fonction ϕ : [a; b] → R. On se donne x 0 ∈ [a; b] et on considère la suite x k+1 = ϕ(x k ) pour k ≥ 0. Si les
deux conditions suivantes sont satisfaites :
1. condition de stabilité : ϕ(x) ∈ [a, b] pour tout x ∈ [a, b]
9
1. Résolution d’équations non linéaires
y y
b b
x1
f
x3
x5
x7
x6
x4
x6
x5
x4 x2
x3
f
x2
x1
2. condition de contraction stricte : il existe K ∈ [0; 1[ tel que |ϕ(x) − ϕ(y)| ≤ K |x − y| pour tout x, y ∈ [a, b]
alors
? ϕ est continue,
? ϕ a un et un seul point fixe x b dans [a, b],
? la suite x k+1 = ϕ(x k ) converge vers x
b pour tout choix de x 0 dans [a, b].
Démonstration.
Continuité La condition de contraction stricte implique que ϕ est continue puisque, si on prend une suite (y k )k∈N ∈ [a, b]
qui converge vers un élément x de [a, b], alors nous avons |ϕ(x)−ϕ(y n )| ≤ K |x − y n | et par suite limk→∞ ϕ(y k ) = ϕ(x).
Existence Commençons par prouver l’existence d’un point fixe de ϕ. La fonction g (x) = ϕ(x) − x est continue dans [a, b]
et, grâce à la condition de stabilité, on a g (a) = ϕ(a) − a ≥ 0 et g (b) = ϕ(b) − b ≤ 0. En appliquant le théorème des
valeurs intermédiaires, on en déduit que g a au moins un zéro dans [a, b], i.e. ϕ a au moins un point fixe dans [a, b].
Unicité L’unicité du point fixe découle de la condition de contraction stricte. En effet, si on avait deux points fixes distincts
xb1 et xb2 , alors
|b x 1 ) − ϕ(b
x 1 − xb2 | = |ϕ(b x 2 )| ≤ K |b
x 1 − xb2 | < |b
x 1 − xb2 |
ce qui est impossible.
Convergence Prouvons à présent que la suite x k converge vers l’unique point fixe xb quand k tend vers +∞ pour toute
donnée initiale x 0 ∈ [a; b]. On a
0 ≤ |x k+1 − xb| = |ϕ(x k ) − ϕ(b
x )| ≤ K |x k − xb|
où K < 1 est la constante de contraction. En itérant k + 1 fois cette relation on obtient
Il est important de disposer d’un critère pratique assurant qu’une fonction ϕ est contractante stricte. Pour cela, rappelons
quelques définitions.
Théorème
Si ϕ : [a; b] → [a; b] est de classe C 1 ([a, b]) et si |ϕ0 (x)| < 1 pour tout x ∈ [a, b], alors la condition de contraction stricte est
satisfaite avec K = max|ϕ0 (x)|.
[a;b]
10
1. Résolution d’équations non linéaires
Démonstration. Considérons la fonction affine g qui transforme l’intervalle [0; 1] dans l’intervalle [x; y] :
g : [0; 1] → [x; y]
t 7→ x + t (x − y)
Alors
Z y Z 1 Z 1
ϕ(x) − ϕ(y) = ϕ0 (ζ) dζ = ϕ0 (g (t ))g 0 (t ) dt = (x − y) ϕ0 (x + t (x − y)) dt
x 0 0
et donc
Z 1
¯ 1¯ 0
¯ ¯ Z
¯ϕ(x) − ϕ(y)¯ = ¯(x − y) 0
ϕ ¯ϕ (x + t (x − y))¯ dt ≤ K ¯(x − y)¯ .
¯ ¯ ¯ ¯ ¯ ¯ ¯ ¯
¯ (x + t (x − y)) dt ¯ ≤ ¯(x − y)¯
¯
0 0
Le théorème de convergence globale assure la convergence, avec un ordre 1, de la suite (x k )k∈N vers le point fixe xb pour
tout choix d’une valeur initiale x 0 ∈ [a; b]. Mais en pratique, il est souvent difficile de déterminer à priori l’intervalle [a; b] ;
dans ce cas, le résultat de convergence locale suivant peut être utile.
Théorème d’O STROWSKI ou de convergence (locale) des itérations de point fixe
Soit [a; b] ∈ R et ϕ : [a; b] → [a; b] une application de classe C 1 ([a; b]). Soit xb ∈ [a; b] un point fixe de ϕ. On peut distinguer
trois cas :
¬ Soit |ϕ0 (b x )| < 1. Par continuité de ϕ0 il existe un intervalle [b x − ε; xb + ε] ⊂ [a; b] sur lequel |ϕ0 (b
x )| < K < 1, donc ϕ est
contractante stricte sur cet intervalle. On a nécessairement ϕ([b x −ε; xb +ε]) ⊂ [bx −ε; xb +ε] et par conséquent la suite
(x k )k∈N converge vers xb pour tout x 0 ∈ [bx − ε; xb + ε]. On dit que xb est un point fixe attractif. De plus,
? si 0 < ϕ0 (b x ) < 1 la suite converge de façon monotone, c’est-à-dire l’erreur x k − xb garde un signe constant quand
k varie ;
? si −1 < ϕ0 (b x ) < 0 la suite converge de façon oscillante, c’est-à-dire l’erreur x k − xb change de signe selon la parité
de k.
Si |ϕ0 (b x − ε; xb + ε] ⊂ [a; b] tel que la suite (x k )k∈N converge vers xb pour tout
x )| > 1, alors il n’existe aucun intervalle [b
x − ε; xb + ε] à l’exception du cas x 0 = xb. On dit que
x 0 ∈ [b xb est un point fixe répulsif. De plus,
? si ϕ0 (bx ) > 1 la suite diverge de façon monotone,
? si ϕ0 (bx ) < −1 la suite diverge en oscillant.
® Si |ϕ0 (b
x )| = 1 on ne peut en général tirer aucune conclusion : selon le problème considéré, il peut y avoir conver-
gence ou divergence.
Exemple
? La fonction ϕ(x) = cos(x) vérifie toutes les hypothèses du théorème d’O STROWSKI : elle est de classe C ∞ (R) et |ϕ0 (x b)| = | sin(xb)| '
0.67 < 1, donc il existe par continuité un intervalle [c, d ] qui contient x tel que |ϕ 0 (x )| < 1 pour x ∈ [c, d ].
p b p
b
? La fonction ϕ(x) = x 2 − 1 possède deux points fixes x b1 p= (1 + 5)/2 et xb2 = (1 − 5)/2 mais ne vérifie l’hypothèse du théorème
d’O STROWSKI pour aucun d’eux puisque |ϕ(xb1,2 )| = |(1 ± 5)/2| > 1. Les itérations de point fixe ne convergent pas.
? La fonction φ(x) = x − x 3 admet x b = 0 comme point fixe. On a φ0 (xb) = 1 et x k → xb pour tout x 0 ∈ [−1; 1] car
? si x 0 = ±1 alors x k = x
b pour tout k ≥ 1,
? si x 0 ∈]−1, 1[ alors on montre que x k ∈]−1, 1[ pour tout k ≥ 1. De plus, la suite est monotone décroissante si 0 < x 0 < 1, monotone
croissante si −1 < x 0 < 0 donc elle converge vers ` ∈ [−1; 1]. Les uniques candidats limites sont les solutions de l’équation ` = φ(`)
et par conséquente x k → 0.
? La fonction φ(x) = x + x 3 admet aussi x
b = 0 comme point fixe. À nouveau φ0 (xb) = 1 mais dans ce cas la suite diverge pour tout choix
de x 0 6= 0.
Le théorème d’O STROWSKI dit que, de manière générale, la méthode de point fixe ne converge pas pour des valeurs arbi-
traires de x 0 , mais seulement pour des valeurs suffisamment proches de xb, c’est-à-dire appartenant à un certain voisinage
de xb. Au premier abord, cette condition semble inutilisable : elle signifie en effet que pour calculer xb (qui est inconnu), on
devrait partir d’une valeur assez proche de xb. En pratique, on peut obtenir une valeur initiale x 0 en effectuant quelques
itérations de la méthode de dichotomie ou en examinant le graphe de f . Si x 0 est convenablement choisi alors la méthode
de point fixe converge.
Proposition Calcul de l’ordre de convergence d’une méthode de point fixe
Soit xb un point fixe d’une fonction ϕ ∈ C p+1 pour un entier p ≥ 1 dans un intervalle [a; b] contenant xb. Si ϕ(i ) (bx ) = 0 pour
1 ≤ i ≤ p et ϕ(p+1) (b
x ) 6= 0, alors la méthode de point fixe associée à la fonction d’itération ϕ est d’ordre p + 1 et
x k+1 − xb ϕ(p+1) (b
x)
lim = .
b)p+1
k→+∞ (x k − x (p + 1)!
11
1. Résolution d’équations non linéaires
x ) = xb et ϕ(i ) (b
où ξ est entre x et xb. Comme ϕ(b x ) = 0 pour 1 ≤ i ≤ p, cela se simplifie et on a
ϕ(p+1) (ξ)
ϕ(x) = xb + (x − xb)p+1 .
(p + 1)!
ϕ(p+1) (ξ)
x k+1 − xb = (x k − xb)p+1 .
(p + 1)!
Lorsque k → +∞, x k tend vers xb et donc ξ, qui se trouve entre x k et xb, tend vers xb aussi. Alors
ϕ(p+1) (b
x)
Pour un ordre p fixé, la convergence de la suite vers xb est d’autant plus rapide que (p+1)! est petit.
Pour choisir la fonction ϕ il est nécessaire de prendre en compte les informations données par les valeurs de f et, éven-
tuellement, par sa dérivée f 0 ou par une approximation convenable de celle-ci (si f est différentiable). Écrivons pour
cela le développement de TAYLOR de f en xb au premier ordre : f (b x − x) f 0 (ξ) où ξ est entre xb et x. Le problème
x ) = f (x)+(b
“chercher xb tel que f (b
x ) = 0”
devient alors
x − x) f 0 (ξ) = 0”.
“chercher xb tel que f (x) + (b
Cette équation conduit à la méthode itérative suivante :
“pour tout k ≥ 0, étant donné x k , déterminer x k+1 en résolvant l’équation f (x k )+(x k+1 − x k )q k = 0, où q k est
égal à f 0 (ξk ) (ou en est une approximation) avec ξk un point entre x k et x k+1 .”
La méthode qu’on vient de décrire revient à chercher l’intersection entre l’axe des x et la droite de pente q k passant par
le point (x k , f (x k )), ce qui s’écrit sous la forme d’une méthode de point fixe avec
f (x k )
x k+1 = ϕ(x k ) ≡ x k − , k ≥ 0.
qk
Considérons maintenant quatre choix particuliers de q k et donc de ϕ qui définissent des méthodes célèbres :
b−a b−a
Méthode de la Corde 1 : qk = =⇒ ϕ(x) = x − f (x)
f (b) − f (a) f (b) − f (a)
f (x)
Méthode de la Corde 2 : q k = f 0 (x 0 ) =⇒ ϕ(x) = x − 0
f (x 0 )
f (x)
Méthode de N EWTON : q k = f 0 (x k ) =⇒ ϕ(x) = x − 0
f (x)
Proposition
Si la méthode de la corde converge, elle converge à l’ordre 1 ; si la méthode de N EWTON converge, elle converge à l’ordre
2 si la racine est simple, à l’ordre 1 sinon.
12
1. Résolution d’équations non linéaires
Démonstration.
f (b)− f (a)
Méthodes de la Corde qk = (méthode 1) ou q k = f 0 (x 0 ) (méthode 2). Si f 0 (b
b−a x ) = 0 alors ϕ0 (b
x ) = 1 et on ne peut pas
0
assurer la convergence de la méthode. Autrement, la condition |ϕ (b x )| < 1 revient à demander que 0 < f 0 (b x )/q k < 2.
Ainsi la pente de la corde doit avoir le même signe que f 0 (bx ) et, pour la méthode 1, l’intervalle de recherche [a; b] doit
être tel quel
f (b) − f (a)
b−a <2 .
f 0 (b
x)
La méthode de la corde converge en une itération si f est affine, autrement elle converge linéairement, sauf dans le
f (b)− f (a)
cas (exceptionnel) où f 0 (b
x ) = b−a (méthode 1) ou f 0 (bx ) = f 0 (x 0 ) (méthode 2), i.e. ϕ0 (b
x ) = 0 (la convergence est
alors au moins quadratique).
Méthode de Newton Soit la méthode de N EWTON pour le calcul de xb zéro de f :
f (x)
ϕ(x) = x − .
f 0 (x)
? Si f 0 (b
x ) 6= 0 (i.e. si xb est racine simple), on trouve
f (x) (x − xb)h(x)
ϕ(x) = 1 − = 1− ,
f 0 (x) mh(x) + (x − xb)h 0 (x)
h(x) m(m − 1)h(x) + 2(x − xb)h 0 (x) + (x − xb)2 h 00 (x)
¡ ¢
1
ϕ0 (x) = ¢2 , ϕ0 (b
x) = 1 − .
m
¡
mh(x) + (x − xb)h 0 (x)
Si la valeur de m est connue a priori, on peut retrouver la convergence quadratique en modifiant la méthode de
N EWTON comme suit :
f (x)
ϕ(x) = x − m 0 .
f (x)
Attention
À noter que même si la méthode de N EWTON permet en général d’obtenir une convergence quadratique, un mauvais
choix de la valeur initiale peut provoquer la divergence de cette méthode (notamment si la courbe représentative de f
présente au point d’abscisse x 0 un tangente à peu près horizontale). D’où l’importance d’une étude préalable soignée de
la fonction f (cette étude est d’ailleurs nécessaire pour toute méthode de point fixe).
(
y = f 0 (x k )(x − x k ) + f (x k ),
y = 0.
On obtient
f (x k )
x = xk −
f 0 (x k )
13
1. Résolution d’équations non linéaires
y
f
y = f 0 (x 1 )(x − x 1 ) + f (x 1 )
xb
x2 x1 x0 x
y = f 0 (x 0 )(x − x 0 ) + f (x 0 )
Soit f : R → R une fonction continue et soit xb ∈ [a, b] un zéro de f . Cette fois-ci, pour calculer x k+1 on prend l’intersection
de l’axe des abscisses avec la droite passant par le point (x k , f (x k )) et parallèle à la droite passant par les points (a, f (a))
et (b, f (b)), i.e. on cherche x solution du système linéaire
f (b)− f (a)
(
y= b−a (x − x k ) + f (x k ),
y = 0,
ce qui donne
b−a
x = xk − f (x k ).
f (b) − f (a)
Il s’agit de la méthode de la corde 1. Cette méthode permet d’éviter qu’à chaque itération on ait à évaluer f 0 (x k ) car on
f (b)− f (a)
remplace f 0 (x k ) par b−a . Une variante de la méthode de la corde consiste à calculer x k+1 comme l’intersection entre
l’axe des abscisses et la droite passant par le point (x k , f (x k )) et parallèle à la droite tangente au graphe de f passant par
le point (x 0 , f (x 0 )), i.e. on cherche x solution du système linéaire
(
y = f 0 (x 0 )(x − x k ) + f (x k ),
y = 0,
ce qui donne
f (x k )
x = xk −
f 0 (x 0 )
Dans cette variante on remplace f 0 (x k ) par f 0 (x 0 ).
y
y f
f
y = f 0 (x 0 )(x − x 0 ) + f (x 0 )
f (b)− f (a)
y = b−a (x − x 0 ) + f (x 0 )
f (b)− f (a)
y = b−a (x − x 1 ) + f (x 1 )
xb
a x1 x2 b x1
x2 x0 x
xb x0 x y = f 0 (x 0 )(x − x 1 ) + f (x 1 )
Exemple
On se trouve en possession d’une calculatrice qui ne sait effectuer que les opérations addition, soustraction et multiplication. Lorsque
a > 0 est donné, on veut calculer sa valeur réciproque 1/a. Le problème peut être ramené à résoudre l’équation x = 1/a ce qui équivaut
à chercher le zéro de la fonction
f : R+
∗ →R
1
x 7→ − a
x
Selon la formule de N EWTON on a
x k+1 = (1 + a)x k + x k2 ,
14
1. Résolution d’équations non linéaires
une récurrence qui ne requiert pas de divisions. Pour a = 7 et partant de x 0 = 0.2 par exemple, on trouve x 1 = 0,12, x 2 = 0,139 2,x 3 =
0,142 763 520 0, x 4 = 0,142 857 081 5, etc. Cette suite converge vers 1/7 ' 0,142 857 142 857.
Exemple Comparaison des méthodes de N EWTON pour différentes formulation de la fonction initiale
Dans R∗
+ on veut résoudre l’équation
x = e 1/x . (1.1)
En transformant l’équation donnée de différentes manières, on arrive à différentes formules de récurrence :
1. L’équation (1.1) équivaut à chercher le zéro de la fonction
f : R∗
+ →R
x 7→ x − e 1/x
f (x ) x − e 1/xk x − e 1/xk
x k+1 = x k − 0 k = x k − k 1/x = x k − x k2 k2 .
f (x k ) 1+ e
k x + e 1/xk k
x k2
g : R∗
+ →R
y 7→ 1 − ye y
et x k = 1/y k .
3. L’équation (1.1) est encore équivalente à chercher le zéro de la fonction
h : R∗
+ →R
x 7→ 1 − x ln(x)
La représentation graphique de f montre qu’il n’existe qu’une seule racine. Comme f (1.7) f (1.9) < 0, elle se trouve dans l’intervalle
[1.7; 1.9]. En partant de x 0 = 1.8 on trouve les suites suivantes :
Critères d’arrêt
Supposons que (x n )n∈N soit une suite qui converge vers xb zéro de la fonction f . Nous avons le choix entre deux types de
critères d’arrêt pour interrompre le processus itératif d’approximation de xb : ceux basés sur le résidu et ceux basés sur
l’incrément. Nous désignerons par ε une tolérance fixée pour le calcul approché de xb et par e n = xb − x n l’erreur absolue.
Nous supposerons de plus f continûment différentiable dans un voisinage de la racine.
Contrôle du résidu : les itérations s’achèvent dès que | f (x n )| < ε. Il y a des situations pour lesquelles ce test
s’avère trop restrictif ou, au contraire, trop optimiste.
x )| ' 1 alors |e n | ' ε : le test donne donc une indication satisfaisante de l’erreur ;
? si | f 0 (b
x )| ¿ 1, le test n’est pas bien adapté car |e n | peut être assez grand par rapport à ε (voir la figure 1.2 à droite) ;
? si | f 0 (b
? si enfin | f 0 (b x )| À 1 alors |e n | ¿ ε et le test est trop restrictif (voir la figure 1.2 à gauche).
Contrôle de l’incrément : les itérations s’achèvent dès que |x n+1 − x n | < ε. Soit (x n )n∈N la suite produite par la
méthode de point fixe x n+1 = ϕ(x n ). Comme xb = ϕ(b
x ) et x n+1 = ϕ(x n ), si on développe au premier ordre on sait
qu’il existe ξn ∈ I xb,xn tel que
15
1. Résolution d’équations non linéaires
y y
f (x k )
f
xk xk f (x k )
xb xb
ek x ek x
F IGURE 1.2.: Deux situations pour lesquelles le résidu e k = x k − xb est un mauvais estimateur d’erreur : | f 0 (x)| À 1 (à gauche),
| f 0 (x)| ¿ 1 (à droite), pour x dans un voisinage de xb.
on en déduit que
x n+1 − x n
en = .
1 − ϕ0 (ξn )
Par conséquent, ce critère fournit un estimateur d’erreur satisfaisant si ϕ0 (x) ' 0 dans un voisinage de xb. C’est le
cas notamment des méthodes d’ordre 2, dont la méthode de N EWTON. Cette estimation devient d’autant moins
bonne que ϕ0 s’approche de 1.
Notons d’ailleurs que si la méthode de point fixe converge avec K < 1 et si on considère le critère d’arrêt |x n+1 −
x n | < ε alors
ε
|e n | = |x n − xb| ≤ ≤ 2ε.
1−K
En effet, il suffit de considérer les inégalités suivantes :
16