Calcul Scientifique et Symbolique en Maths
Calcul Scientifique et Symbolique en Maths
Alain Yger
Institut de Mathématiques, Université Bordeaux 1, Talence 33405,
France
E-mail address: [Link]@[Link]
Version du 9 mai 2014.
Annexe C. TP3 : familiarisation avec MATLAB, travail sur les tableaux 119
Bibliographie 173
Index 175
CHAPITRE 1
Si ce type de langage script exécute les instructions bien plus lentement qu’un lan-
gage compilé (même remarque qu’à propos de Maple, mis à part que les choses
s’avèrent plus cruciales ici car les tableaux en jeu dans les calculs peuvent être
1. Il convient aussi ici de mentionner le logiciel libre Sage (téléchargeable sur le site
[Link] offrant toutes les potentialités du calcul symbolique implémenté
sous Maple, ce sous le langage informatique Python.
2. L’acronyme MATLAB vient d’ailleurs de là : MATLAB pour Matrix Laboratory.
1
2 1. CALCUL SYMBOLIQUE, CALCUL NUMÉRIQUE, ERREURS
Mais l’on a pris la précaution de déclarer 346 comme un nombre décimal (et non plus
comme un entier √ traité de fait comme un symbole), par exemple si l’on demande
l’évaluation de x en x = 346.0, Maple répond cette fois :
18.60107524
On peut ici s’entrainer avec d’autres fonctions Maple (par exemple l’évaluation
eval) en utilisant l’aide interactive disponible dans initiationMaple mentionnée
plus haut. Le logiciel Maple, conçu pour faire du calcul symbolique, c’est-à-dire de
la manipulation d’expresions mathématiques formelles, aura ici été détourné de sa
fonction
√ pour traiter ce calcul comme un calcul scientifique (la valeur fournie ici
pour 346 n’étant bien sûr qu’une valeur tronquée à certaines décimales après la
virgule).
Si l’on demande par contre à MATLAB de retourner le nombre pi en spécifiant format
long, il retournera ce nombre avec ses quinze premières décimales. Le nombre π
est stocké ainsi dans la machine et tout calcul numérique l’impliquant sera auto-
matiquement un calcul avec pertes . On reviendra dans la section suivante sur
la représentation des réels en virgule flottante et le codage machine des nombres ou
des symboles.
25 0
26 0
27 0
28 0
29 0
30 0
Le constat est clair : au delà de la 16-ème itération, la machine ne rend plus le
résultat escompté (qui bien sûr devrait être 1.5000). En fait, ce qui se passe ici
est que la taille des nombres impliqués a nécessité (au niveau de leur codage) un
espace de travail dont le volume dépasse le seuil fixé par la machine : il faut savoir
en effet que pour coder un nombre réel en double précision (ce que fait MATLAB
par défaut), la machine dispose de 64 bits (seulement 32 bits pour le codage en
simple précision). Nous verrons au paragraphe suivant comment ces 64 bits sont
organisés pour le codage des réels en virgule flottante. Mais d’ores et déjà,
nous nous rendons compte que la défaillance du calcul ici est manifestement liée à
des erreurs d’arrondi, la capacité de discernement de la machine ne s’avérant plus
suffisante lorsque les entrées sont des expressions faisant intervenir trop de digits,
comme des entiers positifs de trop grande taille (comme c’est le cas ici lorsque
l’exposant de 10 se met à dépasser 16). Nous reviendrons sur cette question dans
la sous-section 1.2.3.
1.2.2. Un second exemple. Le second exemple proposé, outre qu’il permet
une première prise en main de MATLAB, met en évidence un autre type de faille
du calcul scientifique. À la fois la machine et les mathématiques y sont cette fois
pour quelque chose (nous le verrons beeucoup plus loin dans ce cours lorsque nous
dégagerons la notion de conditionnement).
Rappelons tout d’abord quelques bases concernant la résolution des systèmes linéaires
de N équations à N inconnues
M ·X =B,
où M est une matrice N × N à coefficients réels et B un vecteur colonne à entrées
réelles ; si le rang de la matrice M est égal à N (ce qui signifie det M 6= 0), alors
l’application
X ∈ RN 7−→ M · X
est bijective et l’unique vecteur colonne
x1
..
X= .
xN
de l’équation M · X = B est donné par
X = M −1 · B ,
où
1
(1.1) M −1 = [cofacteurs (M )]t
det M
(At désignant la transposée d’une matrice A).
On verra au chapitre suivant pourquoi la méthode (algorithmique) du pivot s’avère
plus judicieuse pour résoudre un tel système linéaire que la résolution via le calcul de
6 1. CALCUL SYMBOLIQUE, CALCUL NUMÉRIQUE, ERREURS
la matrice M −1 (que nous utiliserons ici). Nous allons dans cette section justement
utiliser MATLAB qui, comme son nom l’indique, un logiciel s’articulant essentiellement
autour du calcul matriciel. On part d’une matrice M très simple que nous déclarons
sous MATLAB (ce calcul est une première initiation à cet environnement, initiation
que vous allez approfondir en TP) :
>> M= [10 7 8 7 ; 7 5 6 5 ; 8 6 10 9 ; 7 5 9 10] ;
M =
10 7 8 7
7 5 6 5
8 6 10 9
7 5 9 10
Le calcul du déterminant de M montre que det M = 1. On déclare un vecteur
colonne B par
>> B = [32 ; 23 ; 33 ; 31];
B =
32
23
33
31
La résolution immédiate du système donne :
>> M^(-1)*B
ans =
1.0000
1.0000
1.0000
1.0000
L’opération ∗ introduite ici correspond à la multiplication usuelle entre deux ma-
trices (ou tableaux) rectangulaires M et N multipliaples , c’est-à-dire telles que
le nombre de lignes de M (matrice de gauche) soit égal au nombre de colonnes de
N (matrice de droite) 6.
Si ce problème numérique correspond à un calcul numérique sollicité par le résultat
d’une expérience, il est vraisemblable que les entrées de l’appareil (représenté par
l’action de la matrice M ) sont connues non pas exactement, mais avec une marge
d’erreur. Il en est de même pour les coordonnées de la sortie B. Perturbons
donc légèrement les entrées de M et tentons à nouveau de résoudre le système. Les
nouvelles matrices M et B perturbées sont les suivantes :
Mperturb =
ans =
1.1693
0.6593
1.1329
0.9322
On voit que l’on est très loin de la solution
1
1
X= 1
1
obtenue lorsque l’entrée M n’est pas perturbée !
On peut envisager de perturber aussi B, par exemple en
>> Bperturb=[32.01 ; 22.99 ; 33.01 ; 30.997]
Bperturb =
32.0100
22.9900
33.0100
30.9970
f−1 · B
Le résultat du calcul de M e est encore pire (on y voit apparaitre une entrée
négative) !
>> Mperturb^(-1)*Bperturb
ans =
1.9353
-0.6112
8 1. CALCUL SYMBOLIQUE, CALCUL NUMÉRIQUE, ERREURS
1.4540
0.7420
Ce n’est pas ici la machine qui est en jeu, mais les mathématiques elles mêmes ! On
peut deviner le problème en réalisant que
>> M^(-1)
ans =
X
N
n= bj β j , bj ∈ {0, ..., β − 1}.
j=0
19 = 1 × 24 + 0 × 23 + 0 × 22 + 1 × 2 + 1 = [1 0 0 1 1].
L’encodage d’un réel en machine, dit en virgule flottante, se fait suivant le stan-
dard IEEE754 (1985, révisé en 2008), soit dans l’un des trois formats binaires
binary32, binary64, binary128, soit dans l’un des deux formats décimaux que
sont decimal64, decimal128.
1.2. REPRÉSENTATION DES NOMBRES EN MACHINE, TYPES D’ERREUR 9
Dans un des systèmes binaires (on prendra comme exemple binary64, dit système
en double précision binaire), un nombre réel x (ou plutôt une valeur approchée x
de ce nombre réel) est encodé sur 1 + 11 + 52 = 64 bits 7 :
– le premier bit est réservé pour le signe du nombre 8 ; on note s ∈ {0, 1} la
valeur de ce bit ;
– les 11 bits suivants servent à encoder l’exposant, i.e l’entier e ∈ Z défini
(lorsque x est non nul) par le fait que 2−e |x| ∈ [1/2, 1[, ce qui signifie que le
nombre 2−e |x| admet un unique développement binaire propre (c’est-à-dire
dont les coefficients ne sont pas tous égaux à 1)
b0 b1 bk
(1.2) 2−e−1 |x| = + 2 + · · · + k+1 + · · · , b0 = 1, bj ∈ {0, 1} ∀ j ∈ N∗ ;
2 2 2
la valeur codée avec ces 11 bits est donc un entier 0 ≤ v ≤ 211 − 1 = 2047
et l’exposant e = v − 1023 peut prendre toute valeur entière entre −1023
(v = 0) et 1024 (v = 2047) ;
– les 52 bits restants sont utilisés pour coder le mot [b1 · · · b52 ], mot de 52
lettres, chacune valant 0 ou 1 ; le mot [1 b1 · · · b52 ] est appelé mantisse 9 de
x.
Les nombres −0 et +0 sont donc encodés respectivement avec v = 0 (tous les bits
b1 , ..., b52 étant mis à 0) tandis que ±∞ et NaN ( Not A Number ) sont encodés
avec v = 2047 (respectivement lorsque tous les bits b1 , ..., b52 sont mis à zéro ou
non).
Exemple 1.1. Par exemple, pour π :
7074237752028440 252 + 2570638124657944 N
π' = = 2 1 +
251 251 252
le numérateur N = 2570638124657944 ≤ 251 s’exprimant ici comme un mot de 52
caractères (0 ou 1) en base 2. On a donc s = 0, e = 1, c’est-à-dire v = 1024, le mot
[b1 · · · b52 ] correspondant à l’écriture en base 2 du nombre
N = 2570638124657944 ≤ 251 .
L’encodage de π sera ainsi (si l’on respecte l’ordre des bits qui a été indiqué)
l’écriture en base 2 du nombre
s × 263 + v × 252 + N = 4614256656552045848
= [0100000000001001001000011111101101010100010001000010110100011000]
Si Nb est le nombre de bits impliqués dans le codage de la mantisse (Nb = 52 en
double précision, Nb = 23 en simple précision, Nb = 112 en quadruple précision),
l’erreur d’arrondi entre une vraie mantisse b et la mantisse codée b est donc
majorée en module par
|b − b| ≤ 2−Nb −1
7. En simple précision, soit dans binary32, le découpage est 1 + 8 + 23 = 32 bits, tandis qu’en
quadruple précision, soit dans binary128, le découpage est 1 + 15 + 112 = 128 bits.
8. On note à ce propos qu’une ambiguité existe pour x = 0 et que la machine nous force à
distinguer −0 et +0.
9. Le premier bit matérialisé en b0 = 1 dans (1.2) est dit bit caché ; mettre ce bit à 0 conduit
à la représentation des nombres dénormalisés ou sous-normaux (subnormal), bien utiles dans les
opérations pour éviter par exemple une division par zéro conduisant brutalement à ±∞.
10 1. CALCUL SYMBOLIQUE, CALCUL NUMÉRIQUE, ERREURS
(on arrondit au nombre codable avec Nb + 1 chiffres le plus proche, qu’il soit
plus petit ou plus grand que m, exactement comme on le fait en décimal). L’er-
reur relative commise lorsque l’on arrondit x = ±b × 2e+1 en x = ±b × 2e+1 =
±[1 b1 · · · bNb ] × 2e+1 est donc majorée par
|x − x| |b − b| |b − b|
= ≤ ≤ 2 × 2−Nb −1 = 2−Nb .
|x| b 1
2
Cette erreur relative 2−Nb (où Nb est le nombre de bits utilisés pour coder la
mantisse), est appelée erreur machine ; c’est elle qui sera responsable des erreurs
d’arrondi dans les calculs (voir la section suivante).
1.2.4. Arrondis dans les calculs entre flottants. On fixe ici une précision
(simple, double ou quadruple), donc un entier Nb (23, 52 ou 115 suivant que l’on
est en simple, double ou quadruple) et un entier Mv (28 − 1, 211 − 1, 215 − 1 suivant
que l’on est en simple, double ou quadruple précision). Si y est un nombre flottant,
i.e un réel du type y = (−1)s [1 b1 . . . bNb ] 2e avec s ∈ {0, 1} et e tel que 0 ≤ v ≤ Mv ,
on pose
ulp(y) := 2e−Nb .
Cette notation (anglosaxonne) vaut pour Unit in the Last Place .
Définition 1.1 (les cinq modes d’arrondi correct du standard IEEE 754). On
dit que y est arrondi au plus proche du nombre réel x
– avec arrondi pair (roundTiesToEven) si y est un nombre flottant tel que
ulp(y)
|y − x| ≤
2
avec la convention que si x est exactement la valeur médiane entre deux
flottants y1 et y2 , alors on privilégie celui dont la mantisse est paire ;
– avec arrondi away (roundTiesToAway) si y est un nombre flottant tel que
ulp(y)
|y − x| ≤
2
avec la convention que si x est exactement la valeur médiane entre deux
flottants y1 et y2 , alors on privilégie y2 si x > 0, y1 si x < 0.
On dit que y est un arrondi dirigé de x
– vers 0 (roundTowardZero) si |y − x| ≤ ulp(y) et |y| ≤ |x| ;
– vers −∞ (roundTowardNegative) si |y − x| ≤ ulp(y) et y ≤ x ;
– vers +∞ (roundTowardPositive) si |y − x| ≤ ulp(y) et y ≥ x.
Les opérations arithmétiques classiques que l’on introduit entre nombres flottants
sont l’addition, la soustraction, la multiplication et la division, que l’on peut implé-
menter en flottant comme :
(a, b) 7→ a ⊕ b := arrondi au plus près de a + b
(a, b) 7→ a b := arrondi au plus près de a − b
(a, b) 7→ a ⊗ b := arrondi au plus près de a × b
(a, b) 7→ a b := arrondi au plus près de a/b.
À la place du choix qui est fait ici, on peut choisir l’un des cinq modes d’arrondi
proposés dans la Définition 1.1.
1.3. LES BOUCLES LOGIQUES (IF/ELSE, FOR, WHILE) SUR QUELQUES EXEMPLES 11
1.3. Les boucles logiques (if/else, for, while) sur quelques exemples
Cette section visant à rappeler quelques rudiments de programmation pourra
aussi être mise à profit pour une famirialisation avec la rédaction de routines .m
sous MATLAB. On y retrouvera sous une forme algorithmique cette fois le bagage
arithmétique introduit en MISMI.
10. La preuve en 1998 par T. Hales de la√conjecture de Kepler (1611) stipulant que densité
de l’empilement cubique à faces centrées (π/ 18 ' .74) maximise la densité d’un empilement de
sphères égales en est un exemple instructif.
12 1. CALCUL SYMBOLIQUE, CALCUL NUMÉRIQUE, ERREURS
Définition 1.2. Si a et b sont deux entiers positifs non tous les deux nuls, on
appelle plus grand diviseur commun (en abrégé PGCD) de a et b le plus grand de
tous les nombres entiers positifs divisant à la fois a et b. Si b = 0, le PGCD de a
et b vaut donc a. On dit que a et b sont premiers entre eux si leur PGCD est égal
à 1. Le PPCM de deux entiers strictement positifs a et b est le plus petit entier
strictement positif multiple à la fois de a et de b et l’on a la relation
PPCM(a, b) × PGCD(a, b) = a × b , ∀ a, b ∈ N∗ .
Le calcul du PGCD de deux nombres entiers positifs a et b avec b non nul fait
apparaitre une démarche mathématique constructive (donc implémentable sur une
machine) que l’on appelle un algorithme. Ce terme vient du surnom Al-Khwarizmi
du mathématicien ouzbek Al Khuwarismi, 780-850, (dont le début du titre d’un des
ouvrages fournit d’ailleurs aussi le mot algèbre).
Notons (comme sous la syntaxe du logiciel de calcul MATLAB que nous utilisons ici)
[q,r] = div(a,b)
l’instruction qui calcule, étant donnés deux nombres entiers tels que b 6= 0, l’unique
couple (q, r) avec r ∈ {0, ..., b−1} tel que a = bq+r, donné par la division euclidienne
de a par b. On pourrait tout aussi bien utiliser les commandes de Maple. La fonction
div s’exprime ainsi :
function [q,r] = div (x,y) ;
q = floor (x./y) ;
r = x- q*y ;
Considérons la suite d’instructions qui conduit au calcul du PGCD de deux entiers
a et b ; le nombre à calculer est noté PGCD dans la suite d’instructions ci-dessous :
function PGCD=PGCD(a,b);
x=a ;
y=b ;
while y>0
[q,r] = div(x,y);
if r==0
PGCD = y;
y = 0 ;
else
[q1,r1] = div(y,r);
x = r;
PGCD = x ;
y=r1 ;
end
end
Si l’on traduit ceci en français, voici la lecture :
fonction PGCD=PGCD(a,b);
x=a;
y=b;
tant que y est non nul, faire
[q,r] = div(x,y);
1.3. LES BOUCLES LOGIQUES (IF/ELSE, FOR, WHILE) SUR QUELQUES EXEMPLES 13
si r=0
PGCD = y;
y=0;
sinon
[q1,r1]= div(y,r);
x=r;
PGCD = x;
y=r1;
fin
fin
Si maintenant, on lit ceci en langage mathématique et de manière exhaustive
(toutes les instructions sont listées et l’on n’effectue pas de ré-assignement des
variables), la démarche que traduit cette routine est :
a = bq + r
b = rq1 + r1
r = r1 q2 + r2
.. .
. = ..
rN −2 = qN rN −1 + rN
rN −1 = rN qN +1 + 0
(comme les restes successifs r, r1 , ... décroissent strictement et qu’il y a un nombre
fini d’entiers entre 0 et b, il vient forcément un moment où rN divise rN −1 , rN étant
le dernier reste obtenu non nul). Le PGCD de a et b est aussi celui de b et r, de r
et r1 , et ainsi, en cascade, de rN et 0 ; il vaut donc rN .
CONCLUSION : Le PGCD de a et b est donc égal au dernier reste non nul rN
dans ce très célèbre algorithme de division dit algorithme d’Euclide.
Par exemple, pour a = 13 et b = 4,
13 = 4 × 3 + 1
4 = 4×1+0
le dernier reste non nul étant ici r0 = 1. On a donc PGCD (a, b) = 1.
La routine MATLAB que l’on a mis en place pour implémenter la fonction PGCD inclut
une boucle de calcul
while y > 0
...
...
...
end
Tant que la variable y garde une valeur strictement positive, une certaine opération
est répétée en boucle. Cette variable y, ré-initialisée après chaque retour de boucle,
matérialise le dernier reste atteint dans la division euclidienne (tant que se reste
reste non nul). A l’intérieur même de cette boucle, figure un nœud de décision :
if r = 0
...
...
14 1. CALCUL SYMBOLIQUE, CALCUL NUMÉRIQUE, ERREURS
...
else
...
...
...
end
La première alternative conduit à l’arrêt du processus précisément lors de la boucle
en cours (puisque on assigne à y la valeur 0 sous cette alternative) tandis que la
seconde appelle la boucle suivante avec des variables x et y qui entre temps ont été
ré-initialisées.
Avec cet exemple de routine PGCD, on a un exemple concret tant de la commande
while que de l’alternative if/else .
Theorème 1.1. Soient a et b deux nombres entiers relatifs non tous les deux
nuls et d leur PGCD ; il existe au moins un couple (u0 , v0 ) ∈ Z2 tel que
au0 + bv0 = d
(une telle relation est appelée identité de Bézout lorsque d = 1).
11. En voici un exemple d’application parlant : c’est grâce à ce type de résultat que l’on peut
combiner deux prises de vue à des temps d’exposition premiers entre eux (disons par exemple
a et b secondes) d’un mobile se déplacant à vitesse constante pour réaliser un cliché instantané
net de ce mobile en mouvement (imaginez deux photos d’une avenue la nuit avec les trainées
lumineuses des phares de voitures, lorsque l’on veut précisément voir se fixer sur la pellicule le
flux de véhicules).
1.3. LES BOUCLES LOGIQUES (IF/ELSE, FOR, WHILE) SUR QUELQUES EXEMPLES 15
On écrit
d = rN −2 − qN rN −1
= rN −2 − qN (rN −3 − qN −1 rN −2 )
= −qN rN −3 + (1 + qN qN −1 )rN −2
.. ..
. = .
= u0 a + v0 b
Plus que l’énoncé du théorème 2.1 lui même, ce qui est très important (parce que
très utile) est qu’il s’agisse d’une assertion dont la démonstration est constructive,
c’est-à-dire s’articule sur un algorithme, dit algorithme d’Euclide étendu. Voici, en
quelques lignes de code, la fonction qui calcule, étant donnés deux entiers relatifs a
et b avec b 6= 0 à la fois le PGCD d de a et b, mais aussi une paire d’entiers relatifs
(u, v) telle que d = au + bv (il s’agit, on le remarquera, d’un algorithme inductif
qui s’auto-appelle) :
fonction [PGCD,u,v]=bezout(a,b);
x=a ;
y= abs(b) ;
[q,r]=div(x,y);
if r==0
PGCD = y;
u=0 ;
v=1;
else
[d,u1,v1]=bezout(y,r);
PGCD=d;
u=v1;
v=sign(b)*(u1- q*v1);
end
On retrouve ici une alternative
if
...
...
else
...
...
end
mais cette fois non plus une boucle
while ...
...
...
...
end
mais (dans un des volets de l’alternative) un ré-appel au programme (ici appelé be-
zout) comme sous-programme de lui-même. Il faut prendre garde ici que la présence
précisément de l’alternative interdit que la procédure ne boucle sur elle même. Cette
16 1. CALCUL SYMBOLIQUE, CALCUL NUMÉRIQUE, ERREURS
alternative devient à une certaine étape un test d’arrêt . La notion dégagée ici est
celle de récursivité en programmation.
jusqu’au degré 4 seulement (sous-section 1.4.2) constituera une transition avec les
méthodes numériques (Newton, sécante, fausse position, dichotomie) qui feront l’ob-
jet du chapitre à venir.
ans =
3.141592653589793
Ce que nous venons de remarquer dans la section précédente à propos des limites
de la représentation des nombres en virgule flottante lorsqu’il s’agit de calcul scien-
tifique (simple ou double précision pour le codage des nombres réels) ne vaut plus
lorsqu’il s’agit de représenter sous forme décimale un nombre rationnel du type
N/D où N et D sont deux entiers strictement positifs sous un logiciel de calcul
symbolique tel que MATLAB ou Mathematica. En effet, ce que peut faire le logiciel
dans ce cas, c’est effectuer l’algorithme de division euclidienne de N par D (celui
que vous connaissez depuis l’école primaire), ce autant de fois qu’on le veut, et
afficher donc un nombre arbitrairement grand de décimales du nombre N/D. On
obtient ici une valeur numérique de N/D avec une précision arbitraire (bien sûr, la
seule limite est la capacité de stockage mémoire de la machine).
Illustrons cela avec une formule d’analyse (que vous prouverez ultérieurement dans
le cours d’Analyse 3 au semestre 4) permettant d’obtenir le nombre π comme limite
1.4. LA FORCE ET LES LIMITES DU CALCUL SYMBOLIQUE 17
648084627718263852275691373137426993671519326356091089553028236126486890826\
6204946805311497184/2330778846653263981508241030629200773293192084180236783\
989076833198484591101069034917139195411721198713640327975684786306003115020\
03066966840474629472267552661106168111125
> v:=4:
for i from 1 to 200 do
v:=v + 4*(-1)^(i)/(2*i+1);
end do:
> v;
29409153714228755488352788592404817226666042608582540000644959155354541277336\
844643332400591621231400189173721343489670702108917127711991109422911095010\
3818394313354582815284/9346423175079588565848046532823095100905700257562749\
18 1. CALCUL SYMBOLIQUE, CALCUL NUMÉRIQUE, ERREURS
503796198101125923210315286830017728173601002006841697715182495993087072491\
2303229853703030326418379288617103573412561125
Les résultats sont ici explicités sous forme de fraction puisque les expressions (1.5)
donnant u = u100 et v = v100 sont des expressions rationnelles et que Maple se
comporte comme un logiciel de calcul symbolique. C’est l’algorithme de division
euclidienne que Maple met en jeu lorsque l’on invoque la commande evalf en
imposant le nombre de décimales (ici par exemple 50). On obtient ici :
> evalf(u,50);
3.1365926848388167504149697050776129667152913517315
> evalf(v,50);
3.1465677471829564012877876601898324180868624240507
Le nombre π est ici encadré par ces deux valeurs. Malheureusement, ce résultat est
très mauvais : il nous permet simplement d’affirmer que la première décimale de π
vaut 1. La raison pour laquelle nous ne pouvons avoir de résultat meilleur ici tient
au fait que la convergence de la suite
Xn
(−1)k
4
2k + 1 n≥0
k=0
vers sa limite, donc des deux suites adjacentes (un )n≥0 et (vn )n≥0 vers π, est ex-
cessivement lente. L’erreur entre la somme de tous les termes d’une série alternée
(c’est de cela qu’il s’agit ici car on ajoute, avec les (−1)k /(2k + 1), des termes dont
le signe alterne avec l’indice k) et la somme de ses n premiers termes est en effet
seulement majorée en module par le premier terme négligé au cran n, c’est-à-dire
le (n + 1)-ième. Ce qui donne ici
X
n
1 4 4
π−4 (−1)k ≤ = ∀ n ∈ N.
2k + 1 2(n + 1) + 1 2n + 3
k=0
On peut aisément montrer qu’en fait ici cette erreur est équivalente à 2/n lorsque
n tend vers l’infini. Il s’agit donc là d’une convergence très lente et on est bien loin
de la convergence exponentielle ! Il faudrait donc travailler avec 2000 termes pour
espérer par exemple une erreur en 10−3 ! On est donc loin du compte avec seulement
199 ou 200 termes comme ce que l’on prend ici.
Pourtant une formule algébrique aussi simple que
π
atan (1) = = 4 atan (1/5) − atan (1/239).
4
1.4. LA FORCE ET LES LIMITES DU CALCUL SYMBOLIQUE 19
> evalf(f,500);
3.141592653589793238462643383279502884197169399375105820974944592307816406286\
208998628034825342117067982148086513282306647093844609550582231725359408128\
481117450284102701938521105559644622948954930381964428810975665933446128475\
648233786783165271201909145648566923460348610454326648211411836337351028535\
076359081339898934439794046156228875441930324313039229345995246595044357372\
067926237220688440078983131937289229037299185437241143101068704481169242885\
4990477850854148747668313991452839287642153241703
> g:=4*(4/5-1/239):
for i from 1 to 200 do
g:=g + 4*((-1)^(i)/(2*i+1))*(4*(1/5)^(2*i+1)- (1/239)^(2*i+1));
end do:
> evalf(g,500);
3.141592653589793238462643383279502884197169399375105820974944592307816406286\
12. Cette formule (et d’autres similaires), basée ici sur l’identité algébrique (1.6), est très an-
cienne. On la doit au mathématicien anglais John Machin (1680-1751), à qui l’on doit le calcul
des cent premières décimales de π. Évidemment aujourd’hui, les outils informatiques, ne serait-ce
que Maple comme ici, permettent d’aller beaucoup plus loin.
20 1. CALCUL SYMBOLIQUE, CALCUL NUMÉRIQUE, ERREURS
208998628034825342117067982148086513282306647093844609550582231725359408128\
481117450284102701938521105559644622948954930381964428810975665933446128475\
648233786783165271201909145648566923460348610454326648213472484619116142621\
584823258499409657096915875845384795905299386889117132202327371985285038293\
873677009013743380924863132735294216568471255262677552078624814206855028421\
6586487825916492887319186809407951507093524812775
Le fait que f100 et g100 aient ici leurs 277 premières décimales égales justifie, puisque
π se trouve entre ces deux nombres, que ces 277 premières décimales soient en
fait les 277 premières décimales de π. En travaillant avec fn et gn avec n >>
100, on déterminerait les décimales de π bien au delà de ce seuil de 277. Toute la
force du calcul symbolique, mettant en jeu ici l’algorithme d’Euclide pour traiter la
division de deux nombres entiers afin de déterminer le développement décimal de
fractions telles f100 ou g100 , a été exploité ici pour obtenir π avec une telle précision.
Ceci n’a rien a voir avec l’encodage de π dans un logiciel de calcul scientifique tel
MATLAB car nous avons mis en œuvre ici (avec la division euclidienne des entiers)
une démarche algorithmique. Il convient de noter tout de fois que la démarche
consistant à remplacer l’approximation (1.4) (trop lente) par l’approximation (1.7)
est, elle, une démarche relevant du calcul scientifique : c’est ce que l’on appelle
l’accélération de convergence. On la retrouvera dans ce cours avec l’interpolation,
le calcul approché d’intégrales et la démarche proposée par le mathématicien néo-
zélandais Alexander Craig Aitken, 1895-1967, dont on reparlera plus loin.
1.4.2. Les commandes solve, fsolve sous Maple ; au delà du degré 5 ?
Les équations algébriques à coefficients complexes sont résolubles par radicaux tant
qu’elles sont de degré inférieur où égal à 4. On connait bien les cas des degrés 1 et
2. Le cas du degré 3 consiste à ramener la résolution de l’équation à celle de
X3 + P X + Q = 0
(ce sont les formules de Cardan 13). La résolution par radicaux des équations de
degré 4 a été proposée par le mathématicien italien Ludovico Ferrari (1522-1565) ;
elle consiste à se ramener au cas des équations du type
x4 + px2 + qx + r = 0
et à ramener la résolution de cette équation à celle (proposée par Cardan pour le
degré 3) de l’équation de degré 3 cette fois :
X 3 + 2pX 2 + (p2 − 4r)X − q 2 .
La commande solve sous un logiciel de calcul symbolique tel que Maple reprend,
lorsque le degré de l’équation à résoudre est inférieur ou égal à 4, donc que la
résolution par radicaux est possible, la démarche explicite proposée par Cardan
13. Le mathématicien italien de la Renaissance Gerolamo Cardano (plus connu sous le nom
de Cardan, 1501-1576), est aussi connu, avec son traité de la théorie des jeux, comme l’un des
précurseurs de la théorie de probabilités.
1.4. LA FORCE ET LES LIMITES DU CALCUL SYMBOLIQUE 21
(1/3)
/ (1/2)\ 1 1
\116 + 12 93 / + ------------------------- - -
(1/3) 3
/ (1/2)\
3 \116 + 12 93 /
/ (1/3) \
1 (1/2) | 1 / (1/2)\ 2 | 1
+ - I 3 |- - \116 + 12 93 / + ------------------------- , --
2 | 6 (1/3)| 12
| / (1/2)\ |
\ 3 \116 + 12 93 / /
(1/3)
/ (1/2)\ 1 1
\116 + 12 93 / + ------------------------- - -
(1/3) 3
/ (1/2)\
3 \116 + 12 93 /
/ (1/3) \
1 (1/2) | 1 / (1/2)\ 2 |
- - I 3 |- - \116 + 12 93 / + -------------------------|
2 | 6 (1/3)|
| / (1/2)\ |
\ 3 \116 + 12 93 / /
La racine réelle et les deux racines complexes conjuguées sont ici séparées par des
virgules. Pour une équation de degré 4, telle que x4 + 3x2 + 5x + 1 = 0, Maple
propose la résolution par radicaux de Ricatti :
> solve(x^4+3*x^2+5*x+1 = 0);
22 1. CALCUL SYMBOLIQUE, CALCUL NUMÉRIQUE, ERREURS
(1/3)
1 / (1/2)\ 22 1 1
-1, - - \244 + 12 1005 / + --------------------------- + -, --
6 (1/3) 3 12
/ (1/2)\
3 \244 + 12 1005 /
(1/3)
/ (1/2)\ 11 1
\244 + 12 1005 / - --------------------------- + -
(1/3) 3
/ (1/2)\
3 \244 + 12 1005 /
/ (1/3) \
1 (1/2) | 1 / (1/2)\ 22 |
+ - I 3 |- - \244 + 12 1005 / - ---------------------------|,
2 | 6 (1/3)|
| / (1/2)\ |
\ 3 \244 + 12 1005 / /
(1/3)
1 / (1/2)\ 11 1
-- \244 + 12 1005 / - --------------------------- + -
12 (1/3) 3
/ (1/2)\
3 \244 + 12 1005 /
/ (1/3) \
1 (1/2) | 1 / (1/2)\ 22 |
- - I 3 |- - \244 + 12 1005 / - ---------------------------|
2 | 6 (1/3)|
| / (1/2)\ |
\ 3 \244 + 12 1005 / /
Cependant, même avec des équations de degré 4, il peut arriver que le logiciel soit
incapable d’effectuer cette résolution par radicaux, pourtant en principe possible,
comme sur l’exemple de l’équation x4 +3x3 +5x+1 où la non-réponse est éloquente :
> solve(x^4+3*x^3+5*x+1 = 0);
1.4. LA FORCE ET LES LIMITES DU CALCUL SYMBOLIQUE 23
/ 4 3 \
RootOf\_Z + 3 _Z + 5 _Z + 1, index = 1/,
/ 4 3 \
RootOf\_Z + 3 _Z + 5 _Z + 1, index = 2/,
/ 4 3 \
RootOf\_Z + 3 _Z + 5 _Z + 1, index = 3/,
/ 4 3 \
RootOf\_Z + 3 _Z + 5 _Z + 1, index = 4/
La résolution par radicaux devient de toutes façon impossible au delà du degré (du
fait du résultat de Galois assurant la non simplicité du groupe An pour n ≥ 5,
comme vous le verrez plus tard) et ceci se reflète dans la réponse de Maple, logiciel
de calcul symbolique, donc formel, à la commande solve. Ainsi sur l’exemple de
l’équation x5 + x3 /2 + 1 = 0 :
/ 5 3 \ / 5 3 \
RootOf\2 _Z + _Z + 2, index = 3/, RootOf\2 _Z + _Z + 2, index = 4/,
/ 5 3 \
RootOf\2 _Z + _Z + 2, index = 5/
En revanche, la commande fsolve fournit pour l’exemple de x5 + x3 /2 + 1 = 0 la
solution réelle (unique) de cette équation de manière approchée, soit :
> fsolve(x^5+(1/2)*x^3+1 = 0, x);
-0.9098248906
Comme nous l’avons vu dans au chapitre 1, l’approche numérique est la seule en-
visageable pour envisager de résoudre des équations algébriques de degré supérieur
ou égal à 5. Ceci est a fortiori encore plus vrai lorsqu’il s’agit d’équations du type
{f (x) = 0}, où f est une fonction non algébrique (en général tabulée dans le noyau
du logiciel interprété), avec tous les problèmes d’évaluation et d’arrondi que cela
peut poser (cf. la sous-section 1.2.4 au chapitre 1). L’objectif de ce chapitre est de
présenter un certain nombre de méthodes itératives conduisant au calcul approché
des racines d’une équation {f (x) = 0} (Newton, fausse position , sécante, dicho-
tomie), de dégager ensuite les outils mathématiques impliqués dans leur validation
(formule de Taylor, pour l’essentiel), enfin de dégager la notion d’ordre.
x 7−→ f (x) = x5 + x3 /2 + 1
sur [−1, −1/2]. Nous savons qu’elle est croissante (strictement), que sa valeur en
−1 vaut −1/2 et que sa valeur en −1/2 vaut −1/32 − 1/16 + 1 > 0. Cette fonction
doit donc (d’après le théorème des valeurs intermédiaires) s’annuler en un unique
point de [−1, −1/2].
f(b)
a=x 0 x1 x2 x3 ξ b
T1
T0
y=f(x)
f(a)
25
26 2. LA RÉSOLUTION NUMÉRIQUE DES ÉQUATIONS NON LINÉAIRES
1. Revoir le cours de MIS 101 (pour l’énoncé, on pourra par exemple se reporter à [Y0], section
3.2) et celui d’Analyse 1 [Yan] (pour la preuve et une étude plus poussée).
2.1. PRÉSENTATION DES MÉTHODES 27
Comme autre exemple d’application, prenons une fonction tabulée sous MATLAB dont
la dérivée est aussi tabulée (f (x) = cos x, f 0 (x) = − sin(x)) et cherchons la valeur
approchée de π en remarquant que la fonction cos est strictement décroissante,
strictement concave sur ]0, 2] et que cos(0) = 1, cos(2) < 0. On peut calculer le
zéro ξ = π/2 de cos sur [0, 2] en partant de init = 2 et en utilisant la méthode
de Newton. On constate qu’au terme de 5 itérations, on obtient les 14 premières
décimales de π.
Il faut noter que, si l’algorithme itératif ainsi décrit converge sous les hypothèses
de stricte monotonie et de stricte concavité (ou convexité) pour la fonction f , sous
lesquelles nous nous sommes placés 2, il n’en est pas de même en général, même si
f 0 garde un signe fixe sur [a, b] ; des changements de convavité du graphe peuvent
nous faire sortir du cadre [a, b]. Même si d’ailleurs la fonction f est de classe C 2 et
strictement monotone (avec f 0 > 0 ou f 0 < 0) sur R tout entier, de tels changements
de concavité peuvent aisément induire la non convergence de l’algorithme lorsqu’il
est initié en un point arbitraire de R ; l’exemple de la fonction
5x − x3
x ∈ [−1, 1] 7→ ,
4
pour lequel la méthode, initiée à −1 ou 1, ne fait que rebondir entre ces deux
valeurs, en est un exemple (notons qu’il y a ici changement de concavité en x = 0).
Voir la figure 2.2.
(1,1)
3
y= (5x−x )/4
−1 0 1
(−1,−1)
Figure 2.2. La méthode de Newton prise en défaut
2. Même si, comme nous le verrons dans la section 2.2.2, ces hypothèses de stricte concavité
ou convexité peuvent être allégées.
3. Philosophe, mathématicien, physicien, esprit universel , l’anglais Isaac Newton (1642-
1727) fut, avec le mathématicien prussien Leibniz, l’un des pères incontestés du calcul différentiel
moderne.
2.1. PRÉSENTATION DES MÉTHODES 29
f(a)
y= f(x)
S1,2
S 0,1
x 0= a x3 x2 x1 = b
ξ
f(b)
soit
x1 − x0
x2 = x1 − f (x1 ) .
f (x1 ) − f (x0 )
Une fois ce calcul de x2 effectué, intervient un processus de choix ([if ... else
...]) : on calcule f (x2 ) et l’on opère comme suit :
– soit f (x2 ) est du signe de f (x1 ) (c’est le cas sur notre figure 2.3), auquel
cas le calcul de x3 se fait comme précédemment, mais cette fois à partir des
valeurs de x0 et x1 réactualisées, c’est-à-dire :
x0 ←− x0
x1 ←− x2
x0 ←− x2
x1 ←− x1
Ce sera d’ailleurs cette seconde formulation (2.3) que nous exploiterons dans la sous-
section 2.2.3 pour préciser la vitesse de convergence de la méthode et comparer les
performances de cet algorithme avec celui de Newton.
On note qu’un des avantages de la méthode de fausse position (par rapport à la
méthode de Newton), indépendamment des questions de rapidité de convergence
dont nous reparlerons, est que la mise en place du calcul itératif (2.2) ou (2.3) ne
nécessite que la connaissance des valeurs (éventuellement tabulées) de la fonction
f , pas celles de sa dérivée f 0 (comme c’était le cas avec l’itération (2.1)).
2.1. PRÉSENTATION DES MÉTHODES 31
f(b)
y=f(x)
x0 =a x3 ξ x4 x2 x1 = b
f(a)
4. La valeur xk est la fausse position du zéro ξ, fausse position que l’algorithme itératif
conduit ici s’emploie à corriger.
32 2. LA RÉSOLUTION NUMÉRIQUE DES ÉQUATIONS NON LINÉAIRES
n’était pas incorporé dans la routine, le résultat fourni par le logiciel (ici, par
exemple, MATLAB) serait NaN ( Not a Number ). Ici, dans le cas où f (xk ) =
f (xk−1 ) aux erreurs machine près, on peut décider du compromis ξ = (xk−1 +xk )/2.
Nous pouvons ici comparer cette routine avec celle régissant la méthode de fausse
position (que nous donnerons plus loin). Il faut noter toutefois que le fait de ne
plus chercher à piéger ξ entre xk−1 et xk (en préalable au calcul de xk+1 ) peut
fort bien nous faire sortir du cadre. Dans les bons cas (voir la figure 2.4), sous
des hypothèses de stricte concavité, la convergence de la suite (xk )k≥0 vers ξ se
présente comme une convergence alternée (xk et xk+1 étant de part et d’autre de
ξ). La convention consistant à poser ξ = (xk−1 +xk )/2 est bien dans ce cas justifiée.
La routine correspondant à l’algorithme de fausse position est similaire, mais intro-
duit en plus un processus de décision. Il convient donc de modifier ainsi la routine
précédente :
function xi=fausseposition(init1,init2,N);
x1=init1;
x2=init2;
for i=1:N
delta = f(x2) - f(x1);
if delta == 0
xi= (x1 + x2)./2;
else
y = x2 - (f(x2) .* (x2-x1))./delta ;
fy = f(y).* f(x1);
if fy < 0
x1 = x1;
x2 = y;
else
x1= x2;
x2= y;
end
end
end
xi=x2;
et initiée donc avec des données initiales x0 , ..., xp−1 , est dite méthode itérative à
p pas. On dit que cette méthode est d’ordre q ∈]0, +∞[ si et seulement si q est la
borne supérieure de l’ensemble des nombres réels strictement positifs κ tels que,
dès que la suite (xk )k≥0 converge vers une limite finie x∞ (zéro d’une fonction
f remplissant, au voisinage de x∞ , les conditions sous lesquelles la méthode est
implémentée), on a
|xk+1 − x∞ |
(2.5) lim sup < ∞.
k→+∞ |xk − x∞ |κ
Remarque 2.1. Souvent, on se contente d’une minoration κ de q, ce qui nous
permet de dire que l’ordre de la méthode est au moins égal à κ. On peut par
exemple choisir κ de manière à ce que la clause (2.5) soit remplie pour toute suite
(xk )k≥0 convergent vers une limite x∞ (zéro de f sous les conditions dans lesquelles
on travaille). Pour prouver que l’ordre est exactement égal à κ, et donc être
sûr que le q = κ ainsi choisi réalise la borne supérieure de l’ensemble de tous les
κ possibles (telle que (2.5) soit satisfaite pour toute suite (xk )k≥0 convergent vers
une limite finie x∞ , zéro réel d’une fonction f sous les conditions sous lesquelles est
envisagé le problème), il convient en plus d’exhiber une suite (xk )k≥0 générée par
l’algorithme itératif (2.4) à partir des p données initiales x0 , ..., xp−1 , convergent
vers une limite finie x∞ (qu’il convient donc de connaitre, ce qui complique ici le
problème, et qui soit zéro de la fonction f sous les conditions imposées), et telle
que, pour cette suite,
|xk+1 − x∞ |
lim inf = l ∈]0, ∞[.
k→+∞ |xk − x∞ |κ
Ceci n’est pas toujours si facile ! En revanche, un minorant de l’ordre q est donné
par tout réel strictement positif κ tel que, pour toute suite (xk )k≥0 générée par
l’algorithme et convergent vers une limite finie x∞ (qui soit un zéro réel de la
fonction f en jeu), il existe une constante C (dépendant a priori de la suite (xk )k≥0 ),
telle que
|xk+1 − x∞ | ≤ C|xk − x∞ |κ , ∀ k ≥ 0 .
2.2.2. La méthode de Newton est au moins d’ordre 2. Que la méthode
de Newton (pour la recherche des zéros réels d’une fonction f de classe C 2 dont la
dérivée ne s’annule pas au voisinage de tels zéros) soit d’ordre au moins égal à 2
repose sur la Proposition 2.1 suivante.
Proposition 2.1. Soit f une fonction de classe C 2 sur l’intervalle ]ξ −r, ξ +r[,
s’annulant en x = ξ, et telle que f 0 ne s’annule pas dans ]ξ − r, ξ + r[. Soit
|f 00 (u)|
γ= sup 0
.
x,u∈]ξ−r,ξ+r[ |f (x)|
Si γr < 2, la suite de Newton (xk )k∈N régie par (2.1), initiée en un point x0
de ]ξ − r, ξ + r[, converge vers ξ lorsque k tend vers l’infini, et l’on a même les
estimations
γ
|xk+1 − ξ| ≤ |xk − ξ|2
2
(2.6) γr 2k −1
|xk − ξ| ≤ |x0 − ξ| ∀ k ∈ N.
2
34 2. LA RÉSOLUTION NUMÉRIQUE DES ÉQUATIONS NON LINÉAIRES
on trouve que |ϕ| ≤ M = K/η sur [−, ]. Si l’on pose ρk = M |ek |, k ∈ N, on
constate que
ρk+1 ≤ ρk ρk−1 , k ≥ 1 .
En prenant le logarithme, il vient
log ρk+1 ≤ log ρk + log ρk−1 .
Si l’on introduit maintenant la suite récurrente (à deux pas) définie par les condi-
tions initiales
α0 = log ρ0 < 0, α1 = log ρ1 < 0
(quitte à supposer que x0 et x1 sont déjà dans ] − 1/M, 1/M [, on peut en effet
supposer ρ0 < 1 et ρ1 < 1) et par la relation inductive
αk+1 = αk + αk−1 ∀ k ≥ 1,
(celle qui donne les nombres de Fibonacci), on constate (par récurrence sur k) que
∀ k ∈ N , log ρk ≤ αk .
6
Mais l’on sait que
1 + √5 k 1 − √ 5 k
k k
αk = λξ + µη = λ +µ ,
2 2
où ξ et η sont les racines de X 2 − X − 1 = 0 ; la constante µ (calculée à partir des
conditions initiales α0 < 0 et α1 < 0) étant strictement négative (car tous les αk le
sont par récurrence), on a
1 + √5 k
αk ≤ −C
2
pour k assez grand pour un certain C > 0, donc, en passant aux exponentielles 7,
1 + √5 k
M |ek | ≤ (e−C ) 2 ,
(toujours pouvu que k soit assez grand). Ce raisonnement montre bien que que
l’ordre de la méthode de la sécante (ou de la méthode de fausse position basée sur
la même formule inductive à deux pas, écrite sous forme symétrique comme (2.3))
est au moins égal au nombre d’or.
6. Voir les cours de MISMI [Y0] et d’Analyse 1 [Yan] pour l’étude des suites récurrentes à
deux pas : uk+1 = auk−1 + buk .
7. On rappelle que log |ek | = log ρk − log M ≤ αk − log M pour tout k ∈ N.
2.3. LA MÉTHODE DE LA DICHOTOMIE 37
Son principe est très simple : on part encore d’un segment [a, b] sur laquelle f vérifie
les propriétés précédentes (f continue, strictement monotone, f (a) et f (b) de signes
opposés). On part de x0 = a, x1 = b et on calcule
x0 + x1
c= .
2
Si f (x0 )f (c) ≤ 0, on pose x2 = x0 et x3 = c et on continue ; si f (x0 )f (c) > 0, on
pose x2 = c et x3 = x1 et on poursuit ainsi, suivant le processus itératif décrit dans
la routine suivante :
function x=dichot(init1,init2,N);
x1=init1;
x2=init2;
for i=1:N
y=(x1+x2)/2;
u=f(x1)*f(y);
if u <= 0
x1=x1;
x2=y;
else
x1=y;
x2=x2;
end
end
x=x2
Notons que l’algorithme sur lequel se fonde la méthode de fausse position inclut
dans son synopsis une procédure de dichotomie.
Cette méthode est seulement linéaire, comme l’est, ce qui justifie l’épithète, le calcul
de y = (xk + xk−1 )/2 en fonction de xk et xk−1 . Elle peut être néanmoins exploitée
pour déterminer les nombres a et b à partir lesquels on s’autorisera à lancer une
méthode plus rapide (par exemple, celle de Newton ou celle de fausse position ou
de la sécante), offrant ainsi un premier dégrossissage de la situation avec une
première localisation grossière du zéro éventuel ξ entre x0 = a et x1 = b. Sur notre
exemple
x3
x 7−→ f (x) = x5 + + 1,
2
la méthode de dichotomie initiée avec x0 = −1 et x1 = −1/2 fournit par exemple :
>> dichot(-1,-.5,3)
-0.875000000000000
>> dichot(-1,-.5,5)
-0.906250000000000
>> dichot(-1,-.5,10)
-0.909667968750000
>> dichot(-1,-.5,20)
-0.909824848175049
>> dichot(-1,-.5,50)
-0.909824890637916
>> dichot(-1,-.5,100)
-0.909824890637916
38 2. LA RÉSOLUTION NUMÉRIQUE DES ÉQUATIONS NON LINÉAIRES
On constate que la vitesse de convergence ici est loin d’être comparable à celle
donnée par l’algorithme de Newton (ou par les algorithmes de fausse position ou
de la sécante). Nous avons en fait la Proposition suivante.
Proposition 2.4 (ordre de la méthode de dichotomie). La méthode de dicho-
tomie est une méthode linéaire à deux pas, d’ordre q = 1.
Démonstration. Comme la méthode repose sur la relation inductive
xk + xk−1
xk+1 =
2
(xk+1 est une fonction linéaire des deux entrées xk et xk−1 ), la méthode de dicho-
tomie (qui est aussi une méthode à deux pas) mérite son qualificatif de linéaire. On
vérifie de plus par récurrence que
1 1
xk+1 − ξ = (xk − ξ) + (xk−1 − ξ)
2 2
et donc, si ek := xk − ξ pour k ∈ N, que
ek + ek−1
ek+1 = .
2
Si ek+1 = O(eqk ) avec q > 1, on aurait
ek = −ek−1 + O(eqk ) ,
donc ek ' −ek−1 , ce qui imposerait q ≤ 1. L’ordre de la méthode de dichotomie ne
peut donc être strictement supérieur à 1. Il est de fait égal à 1 du fait de la linéarité
de la méthode.
CHAPITRE 3
et par le code suivant (sous MATLAB) ; ce code évalue le polynôme dont les coefficients
sont donnés par le vecteur ligne A sur un vecteur ligne d’entrées réelles ou complexes
XX :
function P=Horner(A,XX);
N=length(A);
NN=length(XX);
P=A(1)*ones(1,NN);
for i=2:N
P = XX.* P + A(i);
end
Applioqué à un scalaire X=x (NN=1), cet algorithme consommme N multipli-
cations et N additions (en fait N opérations du type x ∗ (·) + (·)) ; la méthode
naive consistant à calculer toutes les puissances xk , k = 1, ..., N nécessite aussi
N multiplications mais en plus de l’espace mémoire (ce que ne nécessite pas la
démarche algorithmique présentée ci dessus) ; il reste ensuite (toujours si l’on suit
cette méthode naive) N multiplications à effectuer pour évaluer
X
N
ak · xN −k = P (x) ;
k=0
cette méthode naive s’avère donc plus coûteuse (en temps et en espace mémoire).
L’algorithme décrit dans le synopsis présenté plus haut (synopsis synthétisant une
boucle, un test d’arrêt et un branchement) est dit algorithme de Hörner 1 . C’est
lui que, par exemple, les commandes telles que
>> P = [-3 4 1 0 2 -6];
>> y = polyval (P,x);
(sous MATLAB) exploitent pour le calcul du vecteur ligne de scalaires
(P (x0 ), ..., P (xM ))
si x0 , x1 , ..., xM sont M + 1 valeurs (ici nombres réels ou complexes car on travaille
avec un logiciel de calcul scientifique) spécifiées.
Sous Maple13, ce sont les commandes du type
> eval (P(x),x)
> eval (f(x,y),x)
qui réalisent l’évaluation d’un polynôme P (X) (dont les coefficients peuvent tout
aussi bien être des expression symboliques, par exemple des polynômes en une nou-
velle variable Y indépendante de X) sur un scalaire x ou une expression symbolique
(par exemple f (X, Y ), où f est une expression rationnelle en X et Y ).
1. C’est au mathématicien anglais William George Hörner (1786-1837) que l’on attribue cet
algorithme, même si probablement l’origine remonte bien au delà dans l’histoire des mathématiques
et du calcul.
3.1. QUELQUES GÉNÉRALITÉS EN PRISE AVEC L’ALGORITHMIQUE 41
coincide avec l’ensemble des polynômes multiples du dernier reste non nul ∆ dans
l’algorithme de division euclidienne de P1 par P2 (si deg P2 ≥ 0) ou de P2 par P1 (si
deg P1 ≥ 0). Le polynôme non nul ∆ de degré inférieur à max(deg P1 , deg P2 ) ainsi
défini est appelé PGCD (plus grand diviseur commun) de P1 et P2 ; il est défini à
la multiplication par un élément non nul de K près.
C’est en remontant les calculs conduits dans l’algorithme de division euclidienne (de
P1 par P2 ou de P2 par P1 ) que l’on réalise une identité polynômiale (dite identité
de Bézout 2), de la forme
Y
N
P (X) = a0 (X − λj ) ,
j=1
où λ1 , ..., λN sont des nombres complexes ; si λj est répété exactement mj fois dans
cette liste, on dit que l’entier strictement positif mj est la multiplicité de λj comme
zéro du polynôme P ; dire que λj ∈ C est zéro de multiplicité mj ∈ N∗ équivaut à
dire
P (λj ) = P 0 (λj ) = · · · P (mj −1) (λj ) = 0 ,
où P (k) désigne, pour k ∈ N, le k-ème polynôme dérivé de P .
avec X
w(k) = u(k1 )v(k2 ), k = 0, ..., N1 + N2 − 2.
k1 +k2 =k
car X N − 1 = (X − 1)(1 + X + · · · + X N −1 ).
Si les vecteurs-ligne [u(0), ..., u(N1 −1)], [v(0), ..., v(N2 −1)], [w(0), ..., w(N1 +N2 −2)]
sont complétés par des zéros en des vecteurs ligne U, V, W de longueur N , puis
3. Nous envisageons ici les choses du point de vue numérique (car nous ferons appel à des
nombres complexes tels les racines N -ièmes de l’unité dont on ne connait que des approximations).
Cependant, des idées similaires peuvent être introduites dans un cadre plus arithmétique, avec
la notion de transformée de Fourier arithmétique. Signalons aussi qu’un autre algorithme de
multiplication rapide (dans un contexte algébrique) a été proposé en 1960 par le mathématicien
russe Anatolii Karatsuba (1937 –).
4. Pour Discrete Fourier Transform .
44 3. POLYNÔMES, INTERPOLATION, ÉLIMINATION
1 h i
(3.5) t
W = dft (N ) · dft(N ) · t U . ∗ dft(N ) · t V
N
où .∗ désigne (comme par exemple dans les environnements MATLAB ou Scilab) la
multiplication des vecteurs colonne de longueur N entrée par entrée.
Lorsque N = 2p , avec p ∈ N∗ , on doit aux deux ingénieurs informaticiens américains
James William Cooley and John Wilder Tukey (autour de 1965) la construction d’un
algorithme permettant d’implémenter les opérations de multiplication matricielle
1
dft(2p ) · t U, resp. dft(2 c ,
p) · tW
2p
c sont deux vecteurs colonne de longueur 2p donnés, avec p2p−1 (au
où t U et t W
lieu de (2 ) = 22p ) multiplications (en fait seulement (p − 1)2p−1 si l’on prend
p 2
en compte que 2p−1 de ces multiplications sont des multiplications par −1, que
l’on peut donc considérer comme des additions au niveau de la complexité). Nous
énonçons ici ce résultat majeur, moteur de ce qui allait être la révolution numérique
des années 1968-1970.
(l’action de cette matrice, dite action papillon , est donc très simple) et à enchai-
ner, pour calculer une prise de fft(2k ), deux prises de fft(2k−1 ) avec 2k−1 prises
de dft(2). Il est commode de visualiser cette idée (sous-tendant la récursivité) grâce
au diagramme suivant (où WN := exp(−2iπ/N )).
3.1. QUELQUES GÉNÉRALITÉS EN PRISE AVEC L’ALGORITHMIQUE 45
x(0) X(0)
X(1)
x(2)
x(4) FFT X(2)
.... 2
k-1 .. .. .. ..
. . k. . k-1
x(2 -4 ) X(2 -2 )
k k-1
x(2 -2 ) X(2 -1 )
0 k-1
x(1) x W N X(2 )
1 k-1
x(3) x W N X(2 +1 )
FFT 2 k-1
0.8
0.6
0.4
0.2
0
−8 −6 −4 −2 0 2 4 6 8
16k
xk := −8 + , ..., k = 0, ..., 20 .
20
Associons à ces valeurs xk , k = 0, ..., N , des valeurs numériques (y0 , ..., yN ), réelles
ou complexes. Dans l’exemple que nous traitons ici, nous prendrons par exemple
>> y = 1./ (1+ x.^2);
ce qui signifie
y = (y0 , ..., yN )
et
1
yk = , k = 0, ..., N .
1 + x2k
Les points (xk , yk ) affichés ici avec des croix sont des points du graphe de la fonction
1
f : x ∈ [−8, 8] 7−→ .
1 + x2
Comme on le voit sur la figure 3.2, le graphe de la fonction affine par morceaux
interpolant ces points (en plein sur la figure) ne rend pas compte (par exemple sur
[−2, 2]) de la forme du graphe de f , en pointillés sur la figure, qui présente par
exemple une tangente horizontale au point (0, 1).
Le R (resp. C)- espace vectoriel des polynômes à coefficients réels (resp. complexes)
de degré au plus N est un R (resp. C)-espace vectoriel de dimension N + 1 dont une
base est constituée de l’ensemble des monômes {1, X, X 2 , ..., X N }. Ecrire qu’un tel
polynôme
X
N
P (X) = ak X N −k
k=0
prend des valeurs spécifiéees y0 , ..., yN (réelles ou complexes) aux points distincts
x0 , ..., xN revient à écrire les N + 1 contraintes
X
N
P (xk ) = ak xkj = yj , j = 0, ..., N .
k=0
où M (x0 , ..., xN ) est une matrice carrée dont on vérifie que le déterminant (dit de
Vandermonde 5) est égal à
−1
xN0 xN
0 ... x0 1
N −1
xN x1 ... x1 1
1
.. .. .. .. Y
. . . . = (xk − xl ) 6= 0 .
−1
xN
N −1 xN
N −1 ... xN −1 1 0≤k<l≤N
N −1
xN
N xN ... xN 1
Ceci prouve que le système (3.6) est un système de Cramer qui admet une unique
solution (a0 , ..., aN ). Le polynôme
PN (X ; Y ) = a0 X N + a1 X N −1 + · · · + aN
correspondant est l’unique polynôme de degré N prenant les valeurs prescrites des
coordonnées du vecteur ligne Y
Y = (y0 , ..., yN ) ,
(dans cet ordre) respectivement aux points x0 , ..., xN . Ce polynôme est dit polynôme
d’interpolation de Lagrange 6 attaché aux données
{x0 , y0 } , {x1 , y1 } , ... , {xN , yN } .
La matrice M (x0 , ..., xN ) est en général mal conditionnée et l’inverser n’est pas la
bonne stratégie !
Certes, on vérifie immédiatement, si l’on pose
P (X) = (X − x0 ) · · · (X − xN ) ,
que notre polynôme cherché s’écrit
X
N
Pk (X) − Pk (xk )
PN (X; Y ) = yk
Pk0 (xk )(X − xk )
k=0
X
N
Pk (X)
= 0 yk
P (xk )(X − xk )
k=0 k
Q
N
(X − xl )
X
N l=0
l6=k
(3.7) = yk .
Q
N
k=0 (xk − xl )
l=0
l6=k
Tout ce que nous avons dit dans cette section s’adapte au cas où x0 , ..., xN sont
N +1 nombres complexes distincts et y0 , ..., yN nombres complexes. Dans ce nouveau
contexte, avant de clôre cette section, il n’est pas inintéressant de remarquer que
la construction de polynômes d’interpolation de Lagrange fournit, lorsque K = C,
une résolution directe de l’identité de Bézout envisagée dans la sous-section 3.1.2.
On a en effet la
Proposition 3.2 (interpolation de Lagrange et identité de Bézout). Soient P
et Q deux polynômes à coefficients complexes, de degrés respectifs p > 0 et q > 0,
n’ayant tous les deux que des racines simples (λ1 , ..., λp pour P , µ1 , ..., µq pour Q) 7.
Si P et Q n’ont aucune racine commune et si A désigne le polynôme d’interpolation
de Lagrange (de degré exactement q − 1) interpolant les valeurs 1/P (µj ), avec j =
1, ..., q aux points µ1 , ..., µq et B le polynôme d’interpolation de Lagrange (de degré
exactement p − 1) interpolant les valeurs 1/Q(λj ), avec j = 1, ..., p, aux points
λ1 , ..., λp , on a l’identité de Bézout
1 = A(X)P (X) + B(X)Q(X) .
7. Cette condition de simplicité des zéros peut être levée, mais il faut disposer de la notion
de polynôme d’interpolation de Lagrange à des points éventuellement multiples, notion que nous
n’avons pas développé en cours, mais qui pourra être introduite par exemple en TP. En fait, on a
l’identité de Bézout dès que l’on pose (racines multiples ou non pour P ou Q)
X
q h Q(X) i
A(X) = ResY =µj
j=1
(X − Y )P (Y )Q(Y )
X
p h P (X) i
B(X) = ResY =λj
j=1
(X − Y )P (Y )Q(Y )
si le résidu en u = 0 d’une fraction rationnelle R(u) = N (u)/um D(u) (avec D(0) 6= 0) est
le coefficient de u−1 dans le développement de R obtenu après division suivant les puissances
croissantes de N par D. En effet, on constate que A(X)P (X) + B(X)Q(X) s’écrit aussi
X
q h P (X)Q(X) − P (Y )Q(Y )) i X
p h P (X)Q(X) − P (Y )Q(Y )) i
ResY =µj + ResY =λj ≡ 1.
j=1
(X − Y )P (Y )Q(Y ) j=1
(X − Y )P (Y )Q(Y )
50 3. POLYNÔMES, INTERPOLATION, ÉLIMINATION
où P1 , ..., Pm sont m polynômes en n variables sans zéros communs dans Cn 8. Dans
ce cas, le recours à des identités obtenues directement suivant le mécanisme suggéré
par l’interpolation de Lagrange peut s’avèrer dans bien des cas algorithmiquement
moins complexe que la (trop) coûteuse méthode d’élimination.
3.2.2. Différences divisées et méthode récursive d’Aitken. Utiliser la
formule (3.7) n’est pas non plus la meilleure stratégie pour calculer le polynôme d’in-
terpolation de Lagrange (comme d’ailleurs le fait de tenter de résoudre le système
(3.6), ce à cause des probèmes de conditionnement mentionnés auparavant). On
préfère utiliser la méthode (due à Newton) des différences divisées consistant à
chercher à exprimer P (· ; Y ) non dans la base usuelle
{1, X, ..., X N }
des R (resp. C)-polynômes de degré au plus N à coefficients réels (resp. complexes),
mais à l’exprimer dans la base
{1 , (X − x0 ) , (X − x0 )(X − x1 ) , ... , (X − x0 ) · · · (X − xN −1 )}
(on peut prendre N points parmi les N + 1 points xk proposés et choisir un ordre
arbitraire). Afin de spécifier la suite x0 , ..., xN choisie et les valeurs interpolées
y0 , ..., yN , nous noterons le polynôme d’interpolation de Lagrange PN (X; Y ) plus
précisément :
PN (X; Y ) = Lagrange[x0 , ..., xN ; y0 , ..., yN ](X),
Ce polynôme s’exprime alors comme
Lagrange[x0 , ..., xN ; y0 , ..., yN ](X) = y[x0 ] + y[x0 , x1 ] (X − x0 )
+ y[x0 , x1 , x2 ] (X − x0 )(X − x1 ) + · · ·
(3.8) −2 −1
Y
N Y
N
· · · + y[x0 , x1 , ..., xN −1 ] (X − xj ) + y[x0 , x1 , ..., xN ] (X − xj ),
j=0 j=0
où les nombres complexes y[x0 , ..., xn ], n = 0, ..., N , sont appelés différences divisées
des yn par les xn (dans l’ordre imposé n = 0, ..., N ) et sont à extraire d’un tableau
de nombres organisé suivant les règles
y[xn ] = xn
(3.9) y[x0 , ..., xn−1 ] − y[x1 , ..., xn ] , n = 0, ..., N.
y[x0 , ..., xn ] =
x0 − xn
Le calcul des différences divisées s’organise suivant un algorithme triangulaire aisé
à décrire sous forme d’un tableau à N + 2 colonnes, numérotées de −1 à N . La
colonne d’indice k (k = 0, ..., N ) de ce tableau représente la suite des nombres
(k)
un := y[xn , xn+1 , ..., xn+k ], n = 0, ..., N − k. Cette colonne d’indice N est donc
réduite à un élément (y[x0 , ..., xN ]) tandis que la colonne d’indice k = 0 est la
colonne des entrées yn = y[xn ], n = 0, ..., N . Il est commode de placer (pour
mémoire) en position k = −1 la colonne des nœuds x0 , ..., xN , colonne qui sera
rappelée à chaque étape de la progression vers la droite dans la construction du
8. C’est à David Hilbert que l’on doit dans ce cas le résultat impliquant l’existence de tels
polynômes A1 , ..., Am .
3.2. INTERPOLATION DE LAGRANGE ET DIFFÉRENCES DIVISÉES 51
se visualise en calcul numérique (sous la forme d’un tableau triangulaire, comme sur
la présentation algorithmique décrite plus haut) la table des différences divisées :
>> with(Student[NumericalAnalysis]);
>> P:=[[x0,y0],...,[xN,yN]];
>> PI:= PolynomialInterpolation(P,independentvar=’x’,method=newton):
>> DividedDifferenceTable(PI);
En calcul formel, on a par exemple (sous Maple) la syntaxe :
> with (CurveFitting);
> PolynomialInterpolation([x0, ... , xN], [y0, ...,yN], z, form = Newton);
Le lemme qui soutend pareille démarche est dû à A. Aitken 9 :
Lemme 3.2 (lemme d’Aitken). Soient x0 , ..., xN N + 1 nombres réels distincts
et y0 , ..., yN N + 1 nombres réels. Si Q désigne le polynôme d’interpolation de
Lagrange interpolant les valeurs y0 , ..., yN −1 respectivement aux points x0 , ..., xN −1 ,
R le polynôme d’interpolation de Lagrange interpolant y1 , ..., yN respectivement aux
points x1 , ..., xN , on a, si P désigne le polynôme d’interpolation de f aux points
x0 , ..., xN , la formule d’Aitken :
(X − x0 )R(X) − (X − xN )Q(X)
P (X) = .
xN − x0
Démonstration. On note
(X − x0 )R(X) − (X − xN )Q(X)
Pe(X) = .
xN − x0
Il suffit de remarquer que, pour k = 1, ..., N − 1,
(xk − x0 )R(xk ) − (xk − xN )Q(xk )
Pe(xk ) =
xN − x0
(xk − x0 )yk − (xk − xN )yk
= = yk
xN − x0
(puisqu’à la fois Q et R interpolent les valeurs yk aux points x1 , ..., xN ) et que
d’autre part
−(x0 − xN )Q(x0 )
Pe(x0 ) = = Q(x0 ) = y0
xN − x0
(xN − x0 )R(xN )
Pe(xN ) = = R(xN ) = yN
xN − x0
puisque Q interpole y0 au point x0 et que R interpole yN au point xN . Comme le
degré de Pe est exactement égal à N (puisque Q et R sont deux polynômes unitaires
de degré exactement N − 1), Pe coincide bien avec le polynôme d’interpolation de
Lagrange aux points x0 , ..., xN .
9. Mathématicien néo-zélandais, Alexander Craig Aitken, 1895-1967, est aussi connu comme
un calculateur mental prodige ; le lemme cité ici lui est certainement bien antérieur, mais non
sans relation avec les méthodes d’accélération de convergence qu’il développa.
3.2. INTERPOLATION DE LAGRANGE ET DIFFÉRENCES DIVISÉES 53
Plus généralement, voici la procédure d’Aitken, telle que l’on peut, par exemple,
l’implémenter de manière récursive sous MATLAB :
> with(CurveFitting);
[ArrayInterpolation, BSpline, BSplineCurve, Interactive, LeastSquares,
PolynomialInterpolation, RationalInterpolation, Spline, ThieleInterpolation]
> PolynomialInterpolation([[0, 6], [1, 4], [2, 3], [3, 7]], z);
2 3 3 2 7
- z - - z - - z + 6
3 2 6
> PolynomialInterpolation([0, 3, a, 8], [2, 5, b, 7], z, form = Lagrange);
(z - 3) (z - a) (z - 8) 5 z (z - a) (z - 8) b z (z - 3) (z - 8)
- ----------------------- + ------------------- + -------------------
12 a -45 + 15 a a (a - 3) (a - 8)
7 z (z - 3) (z - a)
+ -------------------
320 - 40 a
>> P = polyfit(x,y,N);
On suppose ici que x est distinct de tous les xj (dans le cas contraire, la formule
(3.10) est immédiate). Cette fonction admet m + 1 zéros distincts dans [a, b]. Le
théorème de Rolle appliqué m fois assure que que sa dérivée d’ordre m + 1 s’annule
en un point ξx de ]a, b[, ce qui donne le résultat voulu.
Pour rendre compte des problèmes pratiques que peut soulever l’interpolation de
Lagrange, reprenons l’exemple de la fonction
1 1 1 1
t 7→ = + .
1 + t2 2 t+i t−i
Notons que la dérivée d’ordre N + 1 de cette fonction est
(−1)N +1 (N + 1)! (t + i)N +2 + (t − i)N +2
t 7→ × .
2 (1 + t2 )N +2
56 3. POLYNÔMES, INTERPOLATION, ÉLIMINATION
50
1
−50
0.5
−100
−150
−200
0
−250
−300
−350 −0.5
−8 −6 −4 −2 0 2 4 6 8 −5 −4 −3 −2 −1 0 1 2 3 4 5
10. Carl Runge, mathématicien appliqué (mais aussi théorien) et physicien allemand (1856-
1927) le souligna.
11. Pour une suite discrète {z0 , ..., zN } de nombres complexes, la norme quadratique (ou eu-
P
clidienne) kzk = ( j |zj |2 )1/2 s’interprète en physique comme la racine carrée de l’énergie. Cette
norme dérive d’un produit scalaire sur CN .
3.3. INTERPOLATION POLYNOMIALE APPROCHÉE AU SENS DES MOINDRES CARRÉS57
sont linéairement indépendants, ce qui implique que la matrice constituée avec leurs
colonnes est de rang N + 1 (son déterminant est un déterminant de Vandermonde,
voir les cours d’Algèbre 1 et Algèbre 2). Une base de WM est constituée des res-
trictions à {x0 , ..., xN } des fonctions
Yk : x 7→ xk , k = 0, ..., M .
Ceci s’écrit
N
X X
M
yj − al xlj xkj = 0 , k = 0, ..., M ,
j=0 l=0
15
14
13
12
11
10
4
−0.5 0 0.5 1 1.5
13. Voir le cours d’Algèbre 1 : si P et Q sont premiers entre eux dans K[X] et si P divise un
produit du type B(X) × Q(X), alors P divise B.
14. James Sylvester (1814-1897), avocat, musicologue et mathématicien anglais, ami de Cayley
(vous connaissez depuis le cours d’Algèbre 2 le célèbre théorème de Cayley-Hamilton) ; Sylvester
fut l’un des pionniers de la théorie des déterminants ; c’est à lui d’ailleurs que l’on doit le qualificatif
matrice .
60 3. POLYNÔMES, INTERPOLATION, ÉLIMINATION
> Q := 2*y^2*x^3+y^2*x^4-y*x^3+x*y+1;
3.4. UNE SENSIBILISATION AU PRINCIPE DE L’ ÉLIMINATION ALGÉBRIQUE 61
2 3 2 4 3
Q = 2 y x + y x - y x + x y + 1
> evala(Resultant(3*x^2*y+y^4*x+5*y^2*x^3+8*x*y+2,
2*y^2*x^3+y^2*x^4-y*x^3+x*y+1, x));
20 18 17 16 15 14 13 12 11
y + 20 y + 4 y + 8 y + 340 y - 71 y + 142 y + 1925 y - 156 y
10 9 8 7 6 5
+ 847 y + 3494 y + 1178 y + 1280 y - 328 y + 64 y
> evala(Resultant(3*x^2*y+y^4*x+5*y^2*x^3+8*x*y+2,
2*y^2*x^3+y^2*x^4-y*x^3+x*y+1, y));
17 16 15 14 13 12 11 10
16 x + 120 x + 396 x + 774 x + 972 x + 739 x + 256 x - 29 x
9 8 7 6 5 2
- 34 x - 2 x + 46 x + 34 x - 6 x + x
X
M
Gj (X2 , ..., Xn ) = Aj,k (X1 , ..., Xn )Fk (X1 , ..., Xn ) , j = 1, ..., M 0 .
k=1
62 3. POLYNÔMES, INTERPOLATION, ÉLIMINATION
15. Il faut toutefois souligner que ce procédé devient très coûteux en temps et espace de calcul
dès que le nombre de variables X1 , ..., Xn augmente ! D’autres méthodes, inspirées de l’algorithme
d’Euclide (la plus importante est celle qui passe par la construction de bases de Gröbner), ont été
introduites depuis les années 1970 (à l’aube de la révolution numérique) pour gérer à moindre coût
les systèmes d’équations en beaucoup de paramètres que fait par exemple surgir la robotique.
3.4. UNE SENSIBILISATION AU PRINCIPE DE L’ ÉLIMINATION ALGÉBRIQUE 63
65
66 4. MÉTHODES ITÉRATIVES EN ALGÈBRE LINÉAIRE
Ainsi, si p ≥ 1,
X
k+p−1
kXk+p − Xk k ≤ kXl+1 − Xl k
l=k
κk
≤ κk (1 + κ + · · · + κp−1 )kX1 − X0 k ≤ kX1 − X0 k
1−κ
P
puisque la série géométrique k≥0 κk est convergente (de somme 1/(1 − κ)). La
suite (Xk )k≥0 est donc bien de Cauchy, donc converge vers un point Y , l’erreur
d’approximation de Y par Xk au cran k de l’itération étant majorée par
kX1 − X0 k
kXk − Y k ≤ κk .
1−κ
C’est une convergence très rapide (exponentielle puisque κ peut s’écrire κ = e−u
avec u > 0). Comme T est continue (car contractant, même strictement, les dis-
tances), on déduit de T (Xk−1 ) = Xk que T (Y ) = Y , donc que Y est le point fixe
de T .
2. On pourra aussi consulter avec profit (pour plus de détails) le texte rédigé par Michael
Eisermann [Eiser] dont je me suis inspiré pour cette section (consulter aussi [Eiser0] pour une
version plus piétonne ). On pourra aussi, pour des exemples plus élémentaires (au niveau de
l’enseignement secondaire) se référer au document ressources [TermS] mis en ligne par le Ministère
de l’Education Nationale à l’appui des nouveaux programmes de la spécialité en Terminale S, pages
15 à 26.
4.1. LE THÉORÈME DU POINT FIXE DANS KN (K = R OU C) 67
on a, pour cette norme, kGk ≤ 1. De fait, on a même kGk = 1 puisque 1 est valeur
propre de G.
Supposons un instant qu’aucune page du réseau ne soit une impasse, i.e. ne pointe
sur aucune autre page. Si un robot dépourvu de la moindre capacité de discernement
se déplace sur la toile (à partir d’un instant initial noté k = 0) et que l’on note µk,j
la probabilité que notre robot se trouve sur la page j au clic k, on a
X
N
(4.3) µk+1,i = P clic sur i le robot est en j µk,j ∀ i = 1, ..., N.
j=1
Comme le robot est idiot, il clique au hasard et la probabilité qu’il fasse son (k + 1)-
clic vers la page i depuis la page j (où il est arrivé après k clics) vaut 1/Lj = gji .
Les relations (4.3) se lisent matriciellement
[µk+1,1 , ..., µk+1,N ] = [µk,1 , ..., µk,N ] · G.
On peut certes corriger le fait que certaines pages puissent se révéler être des im-
passes en supposant que toute page pointe sur elle-même, ce que nous ferons. Mais
afin de gommer les effets qu’entraine ce palliatif , on suppose qu’avec une pro-
babilité 1 − κ (on prend couramment .15), le robot, au moment de décider où aller
au bout de k clics, décide d’aller avec la probabilité 1/N vers une page arbitraire
du réseau. Ceci revient à modifier la matrice stochastique G en posant
1−κ
Gκ = ones(N, N ) + κ G ,
N
où ones(N, N ) est la matrice N ×N dont les entrées sont toutes des 1. Les relations
matricielles deviennent
[µk+1,1 , ..., µk+1,N ] = [µk,1 , ..., µk,N ] · Gκ .
Une mesure d’équilibre [µ1 , ..., µN ] (au seuil de tolérance 1−κ) est par définition
une distribution de probabilité sur {1, ..., N } telle que
[µ1 , ..., µN ] = [µ1 , ..., µN ] · Gκ .
En fait, on voit que ceci est équivalent à dire que [µ1 , ..., µN ] (traité comme vecteur
ligne) est un point fixe de l’application affine Tκ de RN dans RN (les vecteurs étant
ici traités en ligne) définie par
1−κ
Tκ : [x1 , ..., xN ] 7→ ones(1, N ) + κ [x1 , ..., xN ] · G.
N
Comme le rayon spectral de G vaut 1, cette application Tκ est κ ' .85 < 1- contrac-
tante et le théorème du point fixe assure l’existence et l’unicité de la mesure
68 4. MÉTHODES ITÉRATIVES EN ALGÈBRE LINÉAIRE
d’équililibre (au seuil de tolérance 1 − κ), en même temps que la possibilité d’ap-
procher asymptotiquement (avec une erreur décroissant exponentiellement) cette
mesure d’équilibre en partant d’une distribution de probabilité arbitraire
(i indice de ligne, j indice de colonne) ; si par contre, on prend comme norme sur
KN la norme k k = k k1 , la norme kAk1 correspondante d’une matric A de taille
(N, N ) à coefficients dans K est
X
N
(4.4) kAk1 = k[ai,j ]k1 = sup |ai,j | .
j
i=1
C’est donc très différent ! Si l’on fait le choix de la norme euclidienne sur KN , la
norme d’une matrice réelle A de taille (N, N ) n’est d’ailleurs pas si simple à expri-
mer. Il faut introduire les valeurs propres (réelles positives) de t A · A (symétrique si
K = R, hermitienne si K = C) et prendre le maximum des racines carrées positives
de ces valeurs propres 6 :
kAk2 = sup{|λ| ; λ2 valeur propre de t A · A}.
4. Pour p ∈]0, 1[, l’inégalité triangulaire est en défaut, et ce n’est donc plus une norme.
5. Qu’il s’agisse d’une norme est ici moins immédiat ; l’inégalité triangulaire est une célèbre
inégalité due au mathématicien russe, géomètre tant des espaces que des nombres, Hermann
Minkowski (1864-1909).
6. Ces racines carrées positives sont appelées valeurs singulières de la matrice A, on y
reviendra.
70 4. MÉTHODES ITÉRATIVES EN ALGÈBRE LINÉAIRE
cependant en rien une obligation). Ce que nous avons dit pour les matrices carrées
vaut donc aussi pour les matrices rectangulaires.
Étant données deux matrices carrées et un choix de norme dans KN , on a toujours :
kA · Bk ≤ kAk × kBk
(ce pour le choix de la norme dérivé du choix de norme qui a été fait dans KN ,
disons par exemple ici la norme euclidienne, mais on peut très bien avoir fait un
autre choix).
4.2.2. Rayon spectral d’une matrice carrée. La notion de rayon spectral
d’une matrice carrée (réelle ou complexe) s’avère une première notion très impor-
tante du pont de vue de l’analyse numérique et du calcul matriciel appliquée. Cette
notion conditionne en général la justification de l’utilisation du théorème du point
fixe.
Soit A une matrice N × N à coefficients réels ou complexes.
Définition 4.1 (rayon spectral). On appelle rayon spectral de la matrice A le
maximum des modules des racines du polynôme caractéristique det[A − XIN ] dans
le corps algébriquement clos C.
La proposition suivante sera pour nous capitale :
Proposition 4.1. Si A est une matrice (N, N ) coefficients complexes de rayon
spectral ρ(A) et si k k est une norme arbitraire sur KN , (K = R ou C), induisant
donc une norme matricielle sur le K-espace des matrices (N, N ), on a toujours
l’inégalité :
(4.5) ρ(A) ≤ kAk .
D’autre part, pour tout > 0, il existe toujours une norme k k() convenable sur
Kn telle que l’on ait, pour la norme matricielle induite,
kAk() ≤ ρ(A) + .
Démonstration. Prouvons d’abord la première assertion. Prenons une valeur
propre (par exemple λ1 ) de module maximum et V un vecteur propre associé. On
a
kA · V k = kλ1 V k = |λ1 | × kV k ≤ kAk kV k
par définition de la norme de A :
kA · Xk
kAk = sup .
X∈K \{0} kXk
N
7. C’est-à-dire, telle que les vecteurs colonnes forment une base orthonormée de KN pour le
produit scalaire usuel, soit t U · U = U · t U = IN .
4.2. QUELQUES NOTIONS PRÉALABLES À L’ALGORITHMIQUE MATRICIELLE 71
avec la même valeur propre. Il suffit donc de prouver la seconde assertion pour une
matrice triangulaire supérieure T . Fixons C assez grand et introduisons la matrice
diagonale DC = diag[1, C, C 2 , ..., C N −1 ]. On remarque que l’on a
(
−1 ti,j C j−i si j ≤ i
(DC · T · DC )i,j =
0 si j > i.
Remarque 4.1 (le cas des matrices diagonalisables sur C). Dans le cas où A est
diagonalisable sur C, il existe une base (V1 , ..., VN ) de vecteurs propres dans CN . En
prenant comme norme d’un vecteur V ∈ KN ⊂ CN le maximum des modules des
coordonnées de V ∈ CN dans cette base (V1 , ..., VN ), un induit sur le K-espace des
matrices (N, N ) une norme matricielle telle que l’on ait exactement kAk = ρ(A).
Le point important concernant le rayon spectral d’une matrice carrée est qu’il peut
être calculé algorithmiquement, par une méthode itérative que nous décrivons ici.
Proposition 4.2 (algorithme itératif pour le calcul approché du rayon spec-
tral). Soit A une matrice à coefficients dans K (K = R ou K = C), diagonalisable 8
sur C, telle que les valeurs propres (distinctes ou confondues) aient des modules
s’organisant comme suit 9
|λ1 | > |λµ+1 | ≥ |λµ+2 | ≥ ... ≥ |λN | ≥ 0
(donc λ1 = λ2 = · · · = λµ ). Soit (V1 , ..., VN ) une base de CN constituée de vecteurs
propres pour A telle que V1 , ..., Vµ soit une base du sous espace propre associé λ1 .
Soit X0 un vecteur
X0 = x1 V1 + · · · + xN VN ,
où l’un au moins des xk , k = 1, ..., µ, est non nul 10. Sous ces hypothèses, l’algo-
rithme itératif initié à X0 et régi ensuite par
A · Xk
Xk+1 = , k≥0
kA · Xk k
(une norme sur KN ayant été arbitrairement choisie) est tel que
lim kA · Xk k = |λ1 |,
k→+∞
8. C’est le cas, avec une probabilité 1, pour une matrice dont les coefficients sont pris au
hasard (suivant une loi uniforme) dans K.
9. Il n’y a donc qu’une seule valeur propre de plus grand module, c’est le cas par exemple pour
une matrice dont toutes les entrées sont positives, ou une matrice dont toutes les valeurs propres
sont réelles (par exemple une matrice symétrique réelle ou une matrice hermitienne).
10. C’est encore le cas d’un X0 pris au hasard (de manière uniforme) dans KN , ce avec une
probabilité égale à 1.
72 4. MÉTHODES ITÉRATIVES EN ALGÈBRE LINÉAIRE
Remarque 4.2. Le choix d’une norme doit être fait au préalable ; sous MATLAB,
la norme que l’on choisit en priorité est la norme euclidienne (ou k k2 ), mais l’on
pourrait prendre aussi n’importe laquelle des normes k kp , p ∈ [1, ∞]. Voici (sous
MATLAB) la routine algorithmique exprimée dans l’énoncé (ici avec la norme eucli-
dienne) :
function r=rayonspectral1(A,X,k)
x=X;
for i=1:k
y=A*x;
x=y/norm(y);
end
r=norm(y);
Démonstration. Montrons d’abord par récurrence sur k que, pour tout entier
k ≥ 1,
Ak · X0
Xk = .
kAk · X0 k
Ceci est vrai pour k = 1 par définition de X1 et on a, pour k ≥ 1,
A · Xk Ak · X kAk+1 · X k −1 Ak+1 · X0
0 0
Xk+1 = =A· × = ,
kA · Xk k kA · X0 k
k kA · X0 k
k kAk+1 · X0 k
ce qui prouve le résultat au cran k + 1. Or
XN h Xµ XN λ k i
j
Ak · X0 = λkj xj Vj = λk1 xj Vj + xj Vj
j=1 j=1 j=µ+1
λ1
X
µ
= λk1 xj Vj + ~k
j=1
avec
k~k k < kX0 k(|λµ+1 |/|λ1 |)k .
On a donc
X
µ X
µ
avec
kX0 k
η̃k ≤ (|λµ+1 |/|λ1 |)k kAk .
P
µ
|λ1 | xj V j
j=1
On achève donc la preuve en combinant les encadrements (4.6) et (4.7). Comme
(|λµ+1 |/|λ1 |)k (qui gouverne la décroissance vers 0 de ηk et η̃k est en e−ρk avec
ρ > 0 (puisque |λµ+1 | < |λ1 |), on a bien une vitesse exponentielle de convergence
de kA · Xk k vers |λ1 |.
12. Algébriste et géomètre allemand Carl Gustav Jacobi (1804-1851) marqua l’essor des
mathématiques au XIX-ème siècle (déterminants, fonctions elliptiques, algèbre linéaire, problèmes
géométriques d’intersection,...).
4.3. LES ALGORITHMES ITÉRATIFS POUR LA RÉSOLUTION DE M · X = B 75
illustre (sur un exemple (3, 3)) l’approche préparatoire de Jacobi, tandis que la
décomposition
5 2 2 5 0 0 0 −2 −2
1 6 3 = 1 6 0 − 0 0 −3
3 4 −8 3 4 −8 0 0 0
Supposons |λ| ≥ 1. Soit i0 tel que |xi0 | = sup1≤i≤N |xi | > 0. De deux choses l’une.
4.3. LES ALGORITHMES ITÉRATIFS POUR LA RÉSOLUTION DE M · X = B 77
– Soit i0 = N , auquel cas la dernière relation dans (4.12), qui s’écrit aussi
X
N −1
mN,N xN = − mN,j xj ,
j=1
Dans tous les cas, quelque soit donc la valeur de i0 ∈ {1, ..., N }, on a mis en
contradiction le fait que
X
|mi,i | > |mi,j | ∀ i = 1, ..., N.
j6=i
L’hypothèse |λ| ≥ 1 est donc absurde, et l’on a bien prouvé ainsi (par l’absurde)
que ρ(Tinf [M ]−1 · F ) < 1.
>> M^(-1)*B
ans =
0.649375600384246
0.625452420067077
0.457215765882871
0.640405560134302
Implémentons plutôt ici l’algorithme de Jacobi :
M =D−E
en isolant les termes diagonaux de M dans la matrice D. Ici
>> Jacobi(M,B,X0,100)
ans =
0.649375600384245
0.625452420067076
0.457215765882871
0.640405560134301
>> Jacobi(M,B,X0,150)
ans =
0.649375600384246
0.625452420067077
0.457215765882871
0.640405560134302
4.3. LES ALGORITHMES ITÉRATIFS POUR LA RÉSOLUTION DE M · X = B 79
>> Jacobi(M,B,X0,200)
ans =
0.649375600384246
0.625452420067077
0.457215765882871
0.640405560134302
Si l’on part d’un autre vecteur X0 :
>> X00= [168 ; -754 ; -432 ; 1218]
>> Jacobi(M,B,X00,50)
ans =
0.649378833791283
0.625456199615299
0.457218432024114
0.640408621290100
>> Jacobi(M,B,X00,100)
ans =
0.649375600384397
0.625452420067254
0.457215765882996
0.640405560134445
>> Jacobi(M,B,X00,150)
ans =
0.649375600384246
0.625452420067077
0.457215765882871
0.640405560134302
>> Jacobi(M,B,X00,200)
ans =
0.649375600384246
0.625452420067077
0.457215765882871
0.640405560134302
On retrouve bien la convergence, ce vers manifestement la solution du système de
Cramer M · X = B.
(ici (e1 , ..., eN ) désigne la base canonique de RN ). Par conséquent la matrice Tinf [S]
est (facilement) inversible (comme matrice triangulaire inférieure à termes diago-
naux tous non nuls). Le fait que l’algorithme de Gauß-Seidel soient opérationnel
pour la résolution de systèmes de Cramer S · X = B lorsque S est une telle matrice
est une conséquence du résultat suivant (qu’il convient ensuite de combiner avec la
Proposition 4.3).
Proposition 4.5 (décomposition de Gauß-Seidel pour une matrice symétrique
réelle définie positive). Soit S une matrice symétrique réelle définie positive de
taille (N, N ) et S = Tinf [S] − F sa décomposition de Gauß-Seidel (4.10). On a
ρ(Tinf [S]−1 · F ) < 1.
Démonstration. Comme S est définie positive, tous les termes diagonaux
hS(ej ), ej i, j = 1, ..., n (e1 , ..., en désignant ici la base canonique de Rn ) sont stric-
tement positifs. La matrice symétrique réelle définie positive S s’écrit
S = Tinf [S] − F = (diag[S] + t T ) + T = (D + t T ) + T,
où T désigne la matrice triangulaire supérieure (avec diagonale nulle) −F . La ma-
trice Tinf [M ] = D + t T est une matrice triangulaire inférieure inversible (en vertu
de la remarque ci-dessus concernant la stricte positivité des coefficients de S situés
sur la diagonale) qui peut s’écrire
√ √ −1 √ √
D + tT = D · D) · (D + t T ) · ( D)−1 · D
√ √ √ √
= D · IN + ( D)−1 · t T · ( D)−1 · D.
Les valeurs propres de la matrice
A = (Tinf [S])−1 · F = −(D + t T )−1 · T
(4.14) √ √ √ −1 √
= −( D)−1 · IN + ( D)−1 · t T · ( D)−1 · ( D)−1 · T
√ √
sont les mêmes (au signe près) que celles de la matrice D · A · ( D)−1 , donc,
compte-tenu de (4.14), les mêmes (au signe près) que celles de la matrice
B = (IN + Te)−1 · t Te,
où √ √
Te := ( D)−1 · t T · ( D)−1 .
Si X est un vecteur propre (dans CN ) de norme 1, associé à une valeur propre
complexe λ de la matrice B, on a
(IN + Te)−1 · t Te · X = λ X,
soit
Te · X = λ(IN + Te) · X = λ X + λ Te · X.
t
telles que
0 ··· ··· 0
.. .. .. .. · t V .
M = U 1 · A . . . . 1
0 ··· ··· 0
La matrice symétrique réelle strictement positive A · A si K = R (resp. hermitienne
t
Supposons que u1 , ..., ur soient les r vecteurs colonnes de la matrice (r, r) U (ortho-
gonale ou unitaire suivant que K = R ou K = C) donnée dans (4.18), complétés par
des zéros en des vecteurs de Kp (ou directement les r premiers vecteurs colonnes
de la matrice (p, p) U orthogonale ou unitaire donnée, c.f. la Remarque 4.3, par la
commande svd sous MATLAB) et v1 , ..., vN les N vecteurs colonnes de la matrice V .
Si Y ∈ Kp , le vecteur
Xr
hY, uj i
X = L−1 app [Y ] = vj
j=1
σj
représente le vecteur X de KN de norme minimale parmi tous les vecteurs x ∈ KN
qui minimisent l’application
x ∈ KN 7−→ kL(x) − Y k22
On l’appelle pseudo-inverse de Y via l’application L. Lorsque L est surjective (c’est-
à-dire r = p ≤ N ), ce vecteur X = L−1app [Y ] figure donc la projection orthogonale
du vecteur nul sur le sous-espace affine de KN défini comme
EL (Y ) := {x ∈ KN ; L(x) = Y }.
16. Voir le cours d’Algèbre 2. Si K = R, il s’agit d’une matrice symétrique réelle et il convient
de remplacer partout matrice unitaire par matrice orthogonale .
84 4. MÉTHODES ITÉRATIVES EN ALGÈBRE LINÉAIRE
Cette démarche itérative, impliquant cette fois le concept d’orthogonalité, a été pro-
posée par le mathématicien polonais Stefan Kaczmarz (1895-1940). Pour prouver
la convergence de la suite (Xk )k≥0 initiée au vecteur nul vers la projection de ce
vecteur sur le sous-espace affine des vecteurs u ∈ KN tels que Lj (u) = Y pour
j = 1, ..., p, on pourra par exemple s’appuyer sur la figure 4.1 suivante.
W3
W
4
X 0,4 X3 0,3
,
X1 X W2
0,2
f
X1,2
X
1,3
0 = X0
W1 x
X
0,1
87
88 5. SCHÉMAS NUMÉRIQUES SIMPLES POUR LA RÉSOLUTION DES EDO
3. De fait, la fonction se prolonge en une fonction de classe C 1 à I0 tout entier, cette fonction
vérifiant d’ailleurs Y 0 (t) = F (t, Y (t)) pour tout t ∈ I0 .
4. On cherche, pour (t0 , y0,0 , ..., y0,N−1 ) = (t0 , Y0 ) ∈ U , les couples (I, y) tels que I soit un
intervalle de R avec (t, y(t), y 0 (t), ..., y (N−1) (t)) ∈ U pour tout t ∈ I,
(5.7) y (N) (t) = F (t, y(t), y 0 (t), ..., y (N−1) (t)) ∀t ∈ I
5.2. RÉSOLUTION NUMÉRIQUE DES EDO 89
où F est une fonction continue dans un ouvert U de RN+1 vérifiant la condition
(5.2) dans cet ouvert, se ramène à celle des équations différentielles Y 0 = F (t, Y ).
Il suffit en effet de remarquer que dire que (I, y) vérifie l’équation d’ordre N
(5.7) avec les conditions initiales (5.8) équivaut à dire que (I, Y ), où Y (t) =
(y(t), y 0 (t), ..., y (N−1) (t)) vérifie le système
Y00 (t) = Y1 (t)
Y10 (t) = Y2 (t)
..
.
(5.9) Yk0 (t) = Yk+1 (t)
..
.
0
YN−2 (t) = YN−1 (t)
0
YN−1 (t) = F (t, Y0 (t), ..., YN−1 (t)).
avec les conditions initiales Y (t0 ) = (y0,0 , ..., y0,N−1 ).
On se place dans le contexte présenté dans la section 5.1, où la condition de Lipschitz
forte (5.4) est supposée remplie par la fonction F . On se limitera aussi ici au cas
N = 1 (le cas général se traitant coordonnée-fonction par coordonnée-fonction ou
matriciellement). On notera donc F (t, y) = f (t, y) pour (t, y) ∈ I0 × R.
5.2.1. Schémas numériques explicites ou implicites : l’exemple d’Eu-
ler. Nous allons introduire ici un principe 5 basé sur la démarche suivante (lorsque
[t0 , t0 + T ] ⊂ I0 ) :
(1) on choisit un pas maximal h0 > 0 et une fonction continue
Φ[f ] : [t0 , t0 + T ] × R × [0, h0 ] −→ R.
(2) pour h = T /N ≤ h0 , on construit la suite récurrente (yh,k )k≥0 solution de
yh,k+1 − yh,k
(5.10) = Φ[f ](tk , yh,k , h), k = 0, ..., N − 1, (ici tk = t0 + k h)
h
et initiée à yh,0 = y0 , y0 étant donné dans R (condition initiale).
L’objectif visé est que, si le pas h est fixé assez petit (en tout cas inférieur à h0 ),
yh,k approche la valeur de la solution f de l’équation différentielle y 0 (t) = f (t, y(t))
(avec condition initiale y(t0 ) = y0 ) au point tk = t0 + k h (pour simplifier, on omet
dans la notation tk la dépendance implicite en h) du maillage
t0 < t + h < t + 2h < · · · < t0 + (N − 1) h < t0 + N h = t0 + T.
On se repose pour cela sur le fait que
yh,k+1 − yh,k
h
et que soient remplies les conditions initiales
(5.8) y(t0 ) = y0,0 , y 0 (t0 ) = y0,1 , ..., y (N−1) (t0 ) = y0,N−1 (conditions initiales),
peut être interprété comme la valeur approchée de t 7→ y 0 (t) soit au point médian de
[tk , tk+1 ] (calcul numérique de dérivée dit centré ), soit au point tk ou tk+1 (cal-
cul numérique de dérivée dit décentré ). En conséquence, il n’est plus question
ici, comme nous l’avions fait jusque là dans les méthodes itératives décrites aux cha-
pitres 2 et 3 (Newton, sécante, dichotomie, Jacobi, Gauß-Seidel, etc.) d’ écraser la
valeur yh,k au fur et à mesure du déroulement de l’algorithme. Ces valeurs yh,0 = y0 ,
yh,1 ,..., yh,k , ..., doivent au contraire être ici stockées en mémoire, ce en vu de l’af-
fichage final du graphe approché de la solution f sur le segment temporel
[t0 , t0 + T ], en l’occurence de celui de l’application discrétisée
k ∈ {0, ..., N } 7−→ yh,k ' f (t0 + kh).
Un tel schéma numérique est dit explicite car le calcul de yh,k+1 se fait à partir de
la connaissance de yh,k tant que k = 0, ..., N − 1.
On peut aussi envisager les schémas implicites où, dans l’étape 2 du processus
décrit, on remplace (5.10) par
yh,k − yh,k−1
(5.11) = Ψ[f ](tk , yh,k−1 , yh,k , h), k = 1, ..., N, (ici tk = t0 + k h),
h
où
Ψ[f ] : (t, y, ξ, h) ∈ [t0 , t0 + T ] × R × R × [0, h0 ] 7−→ Ψ(t, y, ξ, h) ∈ R
est une fonction continue (toujours en maintenant la condition initiale yh,0 = 0) ;
cette fois la détermination de yh,k à partir de yh,k−1 passe par la résolution de
l’équation implicite (d’inconnue ξ)
ξ − yh,k−1
= Φ[f ](tk , yh,k−1 , ξ, h).
h
Qu’il s’agisse de la méthode explicite (basée sur (5.10)) ou implicite (basée sur
(5.11)), la récurrence permettant de calculer de manière inductive les yh,k lorsque
le pas h ≤ h0 est fixé est une récurrence à un terme (yh,k+1 fonction de yh,k ). C’est
la raison pour laquelle on appelle ces méthodes méthodes à un pas.
et
(5.15)
Z b XM
bk+1 − ak+1
tk dt = = I[tk ; [a, b]] = Iapp [tk ; [a, b]] = λj xkj ∀ k = 0, ..., M.
a k + 1 j=0
se révèle exacte pour les fonctions polynomiales de degré 1 (elle devient fausse
par contre pour les fonctions polynomiales de degré 2).
– Le cas où M = 1 et où, toujours par souci de compromis, on prend x0 = a et
x1 = b. Dans ce cas, on trouve λ0 = λ1 = (b − a)/2 et la formule approchée
devient
b−a
(5.18) Iapp [f ; [a, b]] = (f (a) + f (b)).
2
Ce calcul approché correspond à la méthode des trapèzes. Le calcul approché
cesse aussi d’être exact pour les fonctions polynomiales de degré 2.
– Le cas M = 2 où, toujours par souci de compromis, on prend x0 = a,
x1 = (a + b)/2 et x2 = b. On trouve dans ce cas
b−a 2(b − a)
λ0 = λ2 = & λ1 =
6 3
(il est facile de résoudre le système de Cramer (5.15) dans ce cas) et la formule
approchée devient
b−a a + b
(5.19) Iapp [f ; [a, b]] = f (a) + 4 f + f (b) .
6 2
Cette méthode est dite méthode de Simpson 7. Un miracle se produit encore
ici : la formule (5.17) est encore valide (comme on le vérifie aisément) pour
les fonctions polynomiales de degré 3 (elle est fausse par contre pour les
fonctions polynomiales de degré 4).
Si on note (x − t)+ := sup(x − t, 0), l’erreur E[f ; [a, b]] commise entre I[f ; [a, b]] et
Iapp [f ; [a, b]] dans une formule de Newton-Cotes est celle que l’on commet avec la
fonction
Z
1 b
x ∈ [a, b] 7→ (x − t)p+ f (p+1) (t) dt.
p! a
Cette erreur s’exprime aussi (si l’on utilise le théorème de Fubini) comme
Z h i
1 b (p+1)
E[f ; [a, b]] = f (t) E x 7→ (x − t)p+ ; [a, b] dt
p! a
En particulier, en utilisant la formule de la moyenne 9, on trouve, si la fonction
h i
t ∈ [a, b] 7−→ E x 7→ (x − t)p+
8. Ces formules sont apparues à l’occasion du travail de relecture par le mathématicien anglais
Roger Cotes (1682-1716) des Principia d’Isaac Newton.
R R
9. Qui assure ab u(t) v(t) dt = u(ξ) ab v(t) dt pour un certain ξ ∈ [a, b] si u et v sont continues
sur [a, b] et v y garde un signe constant.
10. Ce sera le cas dans nos exemples, on le verra.
94 5. SCHÉMAS NUMÉRIQUES SIMPLES POUR LA RÉSOLUTION DES EDO
Le calcul de E[f : [a, b]] (dans tous les cas de figure pour une méthode de Newton-
Cotes exacte à l’ordre p) conduit (pour une fonction C ∞ sur [a, b]) à une estimation
d’erreur du type
(5.21) |E[f ; [a, b]]| = I[f ; [a, b]] − Iapp [f ; [a, b]] ≤ C[f ] × (b − a)p+2 .
dans laquelle figure déjà sous le premier intégrant la fonction x 7→ (x − t)p+ . Les cal-
culs peuvent être menés explicitement dans les trois premiers exemples mentionnés
(rectangles et trapèzes avec p = 1, Simpson avec p = 3) :
– pour la méthode des rectangles, on trouve, pour t ∈ [a, b],
(
h i (t−a)2
si t ≤ a+b
p
Erect x 7→ (x − t)+ ; [a, b] = (t−b) 2
2
2
2 si t ≥ a+b
2
2 si t ≥ a+b
2
M augmente, certes plus p augmente, mais aussi plus la méthode devient instable
car les coefficients uM,0 , ..., uM,M dans
Z b X
M
f (t) dt ' uM,k f (xM,k )
a k=0
explosent avec M , pour la même raison en fait que celle impliquant les effets
de bord dans l’interpolation de Lagrange 11 sous-jacente aux méthodes de Newton-
Cotes : on ne peut approcher sans risque une fonction dont les dérivées successives
s’amplifient près des extrémités de [a, b] par des fonctions polynomiales dont les
dérivées d’ordre assez grand sont automatiquement identiquement nulles ! L’idée
pratique pour pallier à pareille difficulté consiste à utiliser des méthodes composites.
On découpe [a, b] en M segments [ak , bk ] de même longueur (donc de longueur
h = (b − a)/M ) et, sur chacun de ces segments, on utilise une méthode de Newton-
Cotes d’ordre p. L’ordre de la méthode composite est alors p−1 car il faut multiplier
hp (contrôlant l’erreur sur chaque segment) par M = (b − a)/h, ce qui donne une
erreur en hp−1 .
Le calcul approché I(h) de l’intégrale par le procédé composite associé se présente
donc sous la forme
Z b
(5.25) I(h) = f (t) dt + αhp−1 + o(hp−1 ) .
a
Connaitre explicitement α n’est pas évident, mais si l’on refait le calcul approché
en divisant le pas par 2, on trouve (ceci en effet revient à remplacer h par h/2 dans
(5.25))
Z b
(5.26) I(h/2) = f (t) dt + α(h/2)p−1 + o(hp−1 )
a
(tous les o(hp−1 ) ont ici été notés de manière identique pour alléger les notations).
En soustrayant (5.26) à (5.25), il vient
I(h) − I(h/2)
α(h/2)p−1 ∼
2p−1 − 1
au voisinage de h = 0, ce qui implique donc que l’étude du comportement de
I(h) − I(h/2)
h 7−→ p−1
(2 − 1)(h/2)p−1
lorsque h tend vers 0 nous permet de déterminer explicitement α (au moins de
manière approchée) et d’être par là même capable de contrôler l’ordre de grandeur
(et le signe si α 6= 0) de l’erreur numérique
Z b
f (t)dt − I(h) .
a
Il y a aussi un autre moyen d’exploiter pareille idée, aux fins cette fois d’en tirer une
approximation plus efficace (c’est-à-dire convergeant plus rapidement) de l’intégrale
inconnue. Ce procédé à été initié assez récemment par L.F. Richardson 12 et est
aujourd’hui très utilisé dans le calcul approché des limites de suites ou des sommes
de séries numériques.
On multiplie cette fois la relation (5.26) par 2p−1 et l’on soustrait la relation obtenue
à la formule (5.25), ce qui donne
Z b
(1 − 2p−1 ) f (t) dt = I(h) − 2p−1 I(h/2) + o(hp−1 ) ,
a
ou encore
Z b
I(h) − 2p−1 I(h/2)
f (t) dt = + o(hp−1 ) .
a 1 − 2p−1
Si l’on pose
e
par I(h) p−1
est cette fois en o(h ), alors que l’erreur commise en la remplaçant par
I(h) était en O(hp−1 ), d’où un gain significatif dans la vitesse de convergence de
l’erreur vers 0.
Le procédé d’extrapolation (et d’accélération de convergence) de Richardson (ap-
pliqué ici aux fins de calcul numérique d’intégrales) peut même être itéré si l’on
sait a priori que
Z b
(5.28) f (t) dt = I(h) + α0 hp−1 + α1 hp + · · · + αk hp−1+k + o(hp−1+k ) .
a
En écrivant aussi
Z b
f (t) dt = I(h/2) + α0 (h/2)p−1 + α1 (h/2)p + · · · + αk (h/2)p−1+k + o(hp−1+k )
a
avec p̃ := p + 1, k̃ := k − 1 et
avec p̌ := p̃ + 1 = p + 2, ǩ := k̃ − 1 = k − 2 et
e − 2p̃−1 I(h/2)
I(h) e
ˇ
I(h) := .
1 − 2p̃−1
On est ainsi en position de recommencer, et ainsi de suite jusqu’à ce que l’on ait
épuisé le développement limité. L’opération peut ainsi être itérée k fois et conduit
à une approximation de l’intégrale
Z b
f (t) dt
a
avec une erreur en o(hp−1+k ), tous les termes de la partie principale développement
jusqu’à cet ordre ayant été aspirés . Ceci est d’autant plus important que l’on
sait, d’après une formule sommatoire dite formule d’Euler-MacLaurin 13, que le
terme d’erreur dans la méthode des trapèzes composite se présente précisément
sous la forme
Z b
(5.31) f (t) dt − I(h) = α0 h2 + α2 h4 + · · · + α2(p−1) h2(K−1) + O(h2K ) ,
a
ce qui permet de lui appliquer le processus d’extrapolation de Richardson de manière
itérative comme ci-dessus. Le procédé décrit précédemment et transposé à cet
exemple particulier est ce que l’on appelle aujourd’hui la méthode de Romberg 14.
C’est une méthode très utilisée du fait de son efficacité.
Remarque 5.1 (pourquoi parler d’ extrapolation ?). On peut naturellement
se poser la question suivante : pourquoi parler d’ extrapolation à propos de la
méthode de L. F. Richardson (ou de celle de W. Romberg qui s’en déduit) ? La
raison en est que la quantité
I(h) − 2p−1 I(h/2)
1 − 2p−1
13. Cette formule, établie par Leonhard Euler et le mathématicien écossais Colin MacLaurin
autour de 1735, stipule en effet que, si f est de classe C 2K sur le segment [0, M ],
Z M h f (0) M f (M ) i X b2j (2j−1)
X−1 K
f (t) dt = + f (j) + − f (M ) − f (2j−1) (0)
0 2 j=1
2 k=1
(2j)!
Z M
B2K (t − [t])
+ f (2K) (t) dt ,
0 (2K)!
où les b2 , ..., b2K sont des nombres rationnels indépendants de f (les nombres de Bernoulli d’indices
2, 4, ..., 2K) et B2K un certain polynôme (lui aussi indépendant de f ), dit polynôme de Bernoulli
d’indice 2K. Il suffit d’appliquer cette formule à la fonction
b−a
t 7−→ f (a + ht) , h =
M
(définie sur [0, M ]) pour en déduire le résultat (5.31) voulu. Vous démontrerez plus tard (voir
[Y2], section 5.2.4) la formule importante d’Euler-MacLaurin ; il s’agit d’une formule dans le
même esprit que la formule de Taylor avec reste intégral.
14. Du nom du mathématicien allemand Werner Romberg (1909-2003) qui l’introduisit dans
ses travaux en intégration numérique.
98 5. SCHÉMAS NUMÉRIQUES SIMPLES POUR LA RÉSOLUTION DES EDO
(tk = t0 + kh, k = 0, 1, 2, ...) de proche en proche suivant le principe décrit sur les
quatre étapes suivantes :
tk,1 := t0 + kh = tk
tk,2 := tk,3 = tk,1 + h/2
tk,4 := tk,2 + h/2 = tk,3 + h/2 = tk,1 + h = t0 + (k + 1)h
Ak,1 := f (tk,1 , yh,k )
(5.37) Ak,2 := f (tk,2 , yh,k + h Ak,1 /2) = f (tk,1 + h/2, yh,k + h Ak,1 /2)
Ak,3 := f (tk,3 , yh,k + h Ak,2 /2) = f (tk,1 + h/2, yh,k + h Ak,2 /2)
Ak,4 := f (tk,3 , yh,k + h Ak,3 )
h
yh,k+1 := yh,k + Ak,1 + 2Ak,2 + 2Ak,3 + Ak,4 .
6
Du point de vue algorithmique, ceci peut être codé ainsi sous MATLAB une fois que la
fonction de deux variables f : (t, y) 7→ f (t, y) a été déclarée ’inline’ au préalable
par
>> f = inline (’expression MATLAB en t et y’, ’t’,’y’);
Le but est ici de construire les valeurs approchées de la solution de l’EDO y’=f(t,y)
aux points t0 + (k-1)*T/N, k=1,...,N+1, lorsque la condition initiale est y(t0)
=y0 (t0 désigne ici l’instant initial et T la longueur du segment [t0,t0+T] sur lequel
se trouve échantillonnée la solution y) :
17. On doit cette méthode (initialement introduite pour la résolution de l’équation de la chaleur)
à la mathématicienne britannique Phillis Nicholson (1917-1968) et à son compatriote le physicien
John Craig (1916-2006).
18. Cette démarche a été introduite par Carl Runge, mathématicien et physicien allemand
(1856-1927) et développée numériquement par le mathématicien allemand Martin Kutta (1867-
1944), connu aussi pour ses travaux en aérodynamique.
5.2. RÉSOLUTION NUMÉRIQUE DES EDO 101
voir (5.23)) conduisent à une erreur contrôlée en h3−2 = h2 . Il en est de même pour
le schéma numérique correspondant à Euler modifié (reposant sur l’utilisation de la
formule des rectangles, qui est aussi une méthode d’ordre 3, voir (5.22)). Le schéma
de Runge-Kutta évoqué ci-dessus conduit, on l’a vu, à une erreur en h5−1 = h4 .
20. On se reportera à [Vial] pour une présentation plus détaillée (enrichie de l’aspect
numérique) dont je me suis ici beaucoup inspiré. On pourra aussi consulter avec profit le do-
cument ressources mis en ligne par le Ministère de l’Education Nationale à l’appui des nouveaux
programmes de Spécialité Mathématiques en Terminale S [TermS], pages 54 à 59.
21. Le mathématicien et statisticien autrichien Alfred James Lotka (1880-1949) et le
mathématicien et physicien italien Vito Volterra (1860-1940) l’introduisirent vers 1925, ouvrant
la voie à la dynamique des populations.
5.3. UN MODÈLE DE SYSTÈME AUTONOME : LE MODÈLE PROIE-PRÉDATEUR 103
u=uinit;
v=vinit;
for i=1:N;
u=u+pas*(-(b*c/d)*v -b*u*v);
v=v+pas*((a*d/b)*u +d*u*v);
U=[U,u];
V=[V,v];
end;
Voici maintenant le synopsis de la méthode d’Euler explicite conduisant au tracé des
trajectoires autour du point x = c/d, y = a/b correspondant au point d’équilibre
instable (0, 0) :
function [U,V] = proieinstable(uinit,vinit,pas,a,b,c,d,N);
% [U,V]=proieinstable(uinit,vinit,pas,a,b,c,d,N);
U=uinit;
V=vinit;
u=uinit;
v=vinit;
for i=1:N;
u=u+pas*(a*u -b*u*v);
v=v+pas*(-c*v +d*u*v);
U=[U,u];
V=[V,v];
end;
Ce synopsis correspond à la résolution du système approché :
uk+1 − uk
= uk (a − bvk )
h
(5.41) vk+1 − vk
= vk (−c + duk ),
h
k = 0, ..., N,
initiée en u0 =uinit et v0 =vinit. Pour ce même système, la méthode d’Euler
implicite conduirait à
uk − uk−1
= uk (a − bvk )
h
(5.42) vk − vk−1
= vk (−c + duk ),
h
k = 1, ..., N,
avec les mêmes conditions initiales u0 et v0 , l’expression de (uk , vk ) en termes
de (uk−1 , vk−1 ) à partir des relations (5.42) étant cette fois implicite et non plus
explicite. Que ce soit pour le modèle discret ou pour le modèle continu, on peut
dresser un tableau de variations pour prédire l’évolution du processus depuis une
position donnée. Le plan est ainsi partitionné en neuf regions dans lesquelles on
précise les sens de variation de x et y. Il y a attraction vers (0, 0) le long de l’axe
des y, répulsion le long de l’axe des x. Le calcul numérique (via Euler explicite) 22
fait apparaitre des cycles lorsque l’on initie le processus discret à (x0 , y0 ) avec x0 > 0
22. Avec les routines explicitées ci-dessus sous MATLAB. On suggère les valeurs numériques a = 3,
b = 1, c = 2, d = 1, le pas pas h étant pris égal à .05.
5.3. UN MODÈLE DE SYSTÈME AUTONOME : LE MODÈLE PROIE-PRÉDATEUR 105
et y0 > 0. Les cycles deviennent cependant des cycles épais si (x0 , y0 ) s’écarte
de plus en plus du point d’équilibre stable (c/d, a/b).
12
10
0
0 2 4 6 8 10 12
Ce que l’on observe ici numériquement peut, même au niveau L2, être justifié
théoriquement. Si t 7→ (x(t), y(t)) est solution de (5.39) et que x(t)y(t) reste non
nul, on a
x0 (t) y 0 (t)
(−c + dx(t)) = (a − by(t)).
x(t) y(t)
Ceci s’écrit encore
x0 (t) y 0 (t)
c +a = dx0 (t) + by 0 (t).
x(t) y(t)
En intégrant, on trouve
log |x(t)|c |y(t)|a = d x(t) + b y(t) + const.
On obtient ainsi une intégrale première du système différentiel (5.39). Si on fixe
x(t) = x 6= 0 et la constante const, il apparait qu’il existe un unique y ∈]0, a/b[ tel
que
log y a + c log x = dx + by + const =⇒ a log y − by = const + dx − c log x.
Cette remarque peut être utilisée comme point de départ pour montrer (en ex-
ploitant le tableau de variation) qu’une trajectoire initiée en un point (x0 , y0 ) avec
x0 > 0 et y0 > 0 est automatiquement un cycle.
Au voisinage du point (c/d, a/b) (xn = c/d + un , yn = a/b + vn , un , vn voisins de
(0, 0)), l’étude du système (5.39) linéarisé se ramène à celle de
bc
un+1 = un − τ vn
(5.43) d
ad
vn+1 = vn + τ un .
b
Ici, on se ramène à la résolution d’une équation aux différences à deux pas pour
trouver xn et yn explicitement. On retrouve ici l’algèbre linéaire et le calcul matriciel
(équation caractéristique). On fait une étude similaire au voisinage de l’autre point
d’équilibre (instable cette fois), (0, 0).
ANNEXE A
Exercice A.4 (le changement de base dans les systèmes de numération). Ou-
vrez une nouvelle feuille de travail vierge sous Maple. Apprenez à utiliser les onglets
T et [> pour donner un titre (en texte interte ) à votre feuille de travail (par
exemple : TP1-exo4). Utilisez (après avoir étudié dans le Help le sens de leur fonc-
tion) les routines
> Digits := m ;
> convert (N,binary) ;
> convert (evalf(r),binary) ;
> convert (f,binary) ;
lorsque N, r, f désignent respectivement un entier naturel, un rationnel strictement
positif, un nombre réel flottant strictement positif, pour calculer l’exposant et la
mantisse (voir le cours, section 1.2.3) des nombres suivants : N = 78679, r =
89765432/456843, f = π, lorsque le nombre Digits est fixé d’abord égal à 10,
puis à 20. Sauvez votre session (en lui donnant comme nom TP1-exo4) dans votre
répertoire TPMaple12. Ce sera un fichier .mw.
Exercice A.5 (la génération de fonctions). Ouvrez une nouvelle feuille de
travail vierge sous Maple. Donnez lui un titre (par exemple TP1-exo5). Une fonction
d’une variable réelle (par exemple g) se déclare sur le modèle par exemple de la
fonction x 7→ x2 :
> g := x -> x^2 ;
Vérifiez le avec cette fonction, puis déclarez les fonctions h : x 7→ sin x/(1 + x2 ),
h ◦ g et h ◦ h (exploiter pour cela les commandes @ et subs). Tracez les graphes de
ces fonctions sur [−1, 1], [−10, 10], puis sur R. Sauvez votre session (en lui donnant
comme nom TP1-exo5) dans votre répertoire TPMaple12.
Exercice A.6 (l’utilisation de la palette de gauche expressions ). Ouvrez
une nouvelle feuille de travail vierge sous MAPLE. Donnez lui un titre (par exemple
TP1-exo6).
(1) Calculez la limite en π/4 de la fonction cosinus, puis une approximation
du résultat (avec Digits:=10, puis Digits:=20).
(2) Calculez la limite en 0 de x 7→ (1 − cos x)/x2 .
(3) Tester ce que répond le logiciel lorsqu’on lui demande de calculer une
primitive (intégrale sans bornes précisées) pour
log x
, x 7→ e−x , x 7→ (x2 + 1) log x.
2
x 7→ arctan(x), x 7→
1+x
Dans quel cas (parmi ces quatre exemples) le résultat est-il réellement
concluant et exploitable du point de vue pratique ? Quel est dans ce cas
l’outil que le logiciel est parvenu à mettre en œuvre ?
(4) Calculez la somme des entiers de 1 à n et exploiter la commande factor
pour factoriser cette somme. Evaluer cette somme lorsque n = 1234 en
utilisant eval et subs.
(5) Calculez, pour n ∈ N,
Q
0<k pair ≤2n k
Π(n) = Q ,
0<k impair <2n k
110 A. TP1 : PRISE EN MAIN DES LOGICIELS MAPLE ET MATLAB
Sauvez votre session (en lui donnant comme nom TP1-exo6) dans votre répertoire
TPMaple12.
Exercice A.7 (manipulation de nombres complexes, formule de Machin). Cet
exercice est la mise en pratique d’un exemple qui sera exploité en cours (section
1.4.1 du polycopié). Mais il s’agit ici que vous conduisiez vous même les calculs sous
Maple. Ouvrez une nouvelle feuille de travail vierge sous Maple. Donnez lui un titre
(par exemple TP1-exo7).
(1) Calculez (avec Maple) module et argument des nombres complexes 1 + i
et 2 + 3i.
(2) Vérifiez (toujours sous Maple) la formule de John Machin (mathématicien
anglais, 1680-1751) :
(5 + i)4 = 2(239 + i)(1 + i). (∗)
Peut-on se fier par contre à Maple pour valider la formule
arctan(1) = 4 arctan(1/5) − arctan (1/239) (∗∗)
ou faut-il raisonner mathématiquement à partir de la formule (*) validée,
elle, grâce au logiciel ? Faites le si nécessaire dans ce cas.
(3) On admet ici que la formule (**) se lit aussi :
X (−1)k
n
π
= arctan (1) = lim 4 (1/5)2k+1 − (1/239)2k+1 (∗ ∗ ∗)
4 n→+∞ 2k + 1
k=0
car Z x
dt
arctan(x) = 2
∀x ∈ R
0 1+t
(en particulier pour x = 1/5 et x = 1/239) et que
1 X n
= lim (−1) k 2k
t ∀ t ∈] − 1, 1[
1 + t2 n→+∞
k=0
(en particulier pour tout t ∈ [0, 1/5] ⊃ [0, 1/239]). On pose, pour tout
n ∈ N,
X (−1)k
n
un := 4 4 (1/5)2k+1 − (1/239)2k+1 .
2k + 1
k=0
Montrez que les deux suites (u2n )n≥1 et (u2n−1 )n≥1 sont adjacentes.
(4) Calculez, en utilisant la palette de gauche Expressions sous Maple, les
deux fractions u199 et u200 , puis leurs écritures avec Digits=500. Com-
bien de décimales exactes de π peut-on déduire des résultats affichés ?
Justifiez clairement votre réponse en vous appuyant sur un raisonnement
mathématique (reposant sur la formule (∗ ∗ ∗) partiellement admise).
A. TP1 : PRISE EN MAIN DES LOGICIELS MAPLE ET MATLAB 111
Sauvez votre session (en lui donnant comme nom TP1-exo7) dans votre répertoire
TPMaple12.
Exercice A.8 (familiarisation avec les routines solve et fsolve). Ouvrez
une nouvelle feuille de travail vierge sous Maple. Donnez lui un titre (par exemple
TP1-exo8).
(1) Résolvez les équations :
x2 + x + 1 = 0, x3 + x + 1 = 0, x2 + 3x + 1 > 0, sin x = 0.
Expliquez pourquoi dans les deux premiers cas, on obtient bien toutes
les racines complexes, alors que ce n’est manifestement pas le cas dans le
dernier cas.
(2) Résolvez les inéquations (cette fois dans R) :
x2 + 3x + 1 > 0, x3 + x + 1 > 0.
(3) Dites que répond le logiciel lorsqu’on lui demande de résoudre
x5 − 2x4 + 3x2 − 3/2 = 0, x ∈ C,
soit l’instruction :
> solve (x^5 - 2*x^4 + 3*x^2 - 3/2 =0,x) ;
Est-ce vraiment une surprise ? Dites ce qui se passe par contre avec la
commande :
> solve(x^5 - 2*x^4 + 3*x^2 - 1.5 = 0,x);
Expliquez pourquoi la réaction du logiciel est dans ce cas différente. Que
se passe t-il sur cet exemple si l’on remplace la routine solve par fsolve ?
Préciser ce qui se passe si l’on utilise la commande avec option :
> fsolve (x^5 - 2*x^4 +3*x^2 - 3/2 = 0,x,complex);
Sauvez votre session (en lui donnant comme nom TP1-exo8) dans votre répertoire
TPMaple12.
Exercice A.9 (un (tout) premier contact avec MATLAB). Fermez l’environne-
ment Maple et lancez MATLAB avec la commande matlab depuis la fenêtre de votre
Terminal. Une fois MATLAB ouvert, placez vous, en utilisant l’onglet dans le bandeau
supérieur, dans votre répertoire TPMATLAB.
(1) Ouvrez (en utilisant l’onglet file) le fichier test1.m et analysez-en la
syntaxe. Que devrait être le résultat final ?
(2) Lancez test1 derrière le prompt de MATLAB. Qu’observez vous ? Sachant
que le logiciel code les réels en double précision (binary64), peut-on
prédire à partir de quel seuil d’itérations N le résultat de test1 n’est
plus fiable ? Calculez pour cela log10 (252 ) ; qu’observez vous ?
(3) Ouvrez (toujours en utilisant l’onglet file) le fichier PGCD.m. Essayez
d’analyser comment l’algorithme présenté ici traduit la démarche mathé-
matique conduisant au calcul du PGCD (souvenez vous de l’algorithme
d’Euclide présenté en MISMI). Quel sous-programme est appelé dans cette
procédure ? Testez la commande PGCD en le lançant derrière le prompt :
>> d=PGCD (a,b);
112 A. TP1 : PRISE EN MAIN DES LOGICIELS MAPLE ET MATLAB
(avec des valeurs des entiers naturels a et b > 0 fixées). Testez aussi le
sous-programme que cette procédure appelle.
(4) Vérifiez que la matrice
2 −4 −1 −3
A= 1 0 −5 7
−6 −8 0 7
se déclare sous MATLAB ainsi :
>> A = [ 2 - 4 - 1 -3 ; 1 0 - 5 7 ; -6 - 8 0 7] ;
>> A
Déclarez, si vous avez compris, une matrice B cette fois à 4 lignes et 3
colonnes (mettez des entrées réelles arbitraires). Que se passe -il lorsque
vous tentez les commandes :
>> C = A*B;
>> C = B*A;
>> BB = B’;
>> C = BB*A;
>> C = A*BB;
>> C = A.*BB;
>> C = BB.*A;
Dans chaque cas où vous n’avez pas reçu un message d’erreur, affichez BB
ou C en faisant
>> BB
>> C
(sans point virgule cette fois) et tentez d’analyser ce qui se passe. Les
bases de calcul matriciel acquises en S2 vous seront ici utiles.
Sauvez votre session (elle le sera automatiquement dans le répertoire TPMATLAB sous
lequel vous travaillez sous forme d’un fichier .mat) ainsi (vérifiez bien après que ce
fichier a été créé) :
>> save TPMATLAB-1
Exercice A.10 (préparation au TP suivant : Euclide étendu, boucles et récur-
sivité).
(1) Révisez dans votre cours de MISMI la démarche conduisant à la recherche
d’une solution de l’identité de Bézout : étant donnés deux entiers relatifs
a et b avec b 6= 0, trouver deux nombres entiers u et v tels que
au + bv = PGCD (|a|, |b|).
(2) Ouvrez sous MATLAB la routine bezout.m (en vous mettant dans le répertoire
TPMATLAB).
(3) Réfléchissez sur le synopsis de la procédure (ici récursive, on notera qu’elle
s’auto-appelle) écrite dans ce programme. Essayez de la relier à la démarche
mathématique vue en MISMI pour trouver une solution (u, v) à l’identité
de Bézout (et en même temps d’ailleurs le PGCD de |a| et |b|).
(4) Faites des tests pour valider cette routine sur des exemples.
ANNEXE B
ête éliminé). Dans les boucles do les éléments de syntaxe by ... et while [...]
peuvent, le cas échéant, être omis si leur présence n’est pas justifiée par l’algorithme.
On notera que l’instruction : au lieu de ; lors de l’instruction end évite les affichages
intermédiaires au fur et à mesure de l’exécution du code.
La réalisation de codes sous Maple12 sera toujours faite ici en se positionnant
derrière le prompt et en suivant sous les onglets du bandeau le chemin
Insert → Code Edit Region
L’exécution du code, une fois celui ci rédigé dans la fenêtre adéquate, se fait avec
le clic droit.
end do:
for i from 0 to 7 do
[i,i*i];
ifactor (i!);
‘-----‘;
end do;
X
N
k2
k=1
X
N
k3
k=1
Exercice B.4 (boucles do ... end do et if ... else ... end if enchainées
et conjecture de Syracuse). Cet exercice est centré sur une célèbre conjecture, dite
Conjecture de Syracuse, soulevée autour de 1930 par le mathématicien allemand
Lothar Collatz. Voir par exemple, pour une histoire de cette conjecture et l’ état
de l’art aujourd’hui, le site
[Link]
116 B. TP2 : ASSIGNATION/RÉASSIGNATION DE VARIABLES, BOUCLES (MAPLE)
Cette conjecture stipule que, pour tout entier N ∈ N∗ , la suite générée par récurrence
suivant
u(1) = N
(
u(k)
2 si u(k) est pair (soit type(u(k),even) = true)
u(k + 1) =
3u(k) + 1 si u(k) est impair (soit type(u(k),odd) = true)
∀k ≥ 1
finit par atteindre en un temps fini la valeur 1.
(1) Expliquez pourquoi, une fois cette valeur 1 atteinte pour k = k(N ), la
suite va nécessairement répéter indéfiniment le motif
1 4 2 1 4 2 1 4 2 1 4 2 1 4 2 1 ...
(2) Ouvrez une nouvelle feuille de travail, que vous intitulerez (en mode Texte)
Conjecture de Syracuse. En utilisant la syntaxe
u:= N :
for i from 1 to M while u > 1 do
if ... then
...
else
...
end if :
...
end do ;
réalisez un code affichant les valeurs prises par la suite (u(k))k≥1 (initiée
à la valeur N ) tant que k ≤ M et que u(k) reste différent de 1, ainsi que
la valeur du premier cran 1 ≤ k ≤ M (s’il existe) tel que u(k) = 1. Testez
ce code avec les valeurs de N = 127, N = 537, N = 1235 en prenant (ici
arbitrairement) M = 500 comme marge de sécurité . La conjecture
de Syracuse est-elle bien validée pour ces valeurs de N ? Que vaut dans
chacun de ces trois cas le cran k ≥ 1 tel que u(k) prenne la valeur 1 pour
la première fois ? Peut on affirmer ce cran est de plus en plus grand au fur
est à mesure que N augmente ?
(3) Modifiez le code que vous venez d’élaborer pour que soient affichées uni-
quement au terme de son exécution :
– le cran k (s’il existe entre 1 et M ) où u(k) prend la valeur 1 pour
la première fois ;
– la valeur maximale prise par l’entier u(k) avant d’atteindre cette
valeur 1.
Sauvez votre feuille de travail en l’enregistrant sous votre répertoire TPMaple12
comme TP2-Syracuse.
Exercice B.5 (Algorithmes d’Euclide et d’Euclide étendu). Ouvrez une nou-
velle feuille de travail sous Maple12.
(1) Apprenez à vous familiariser avec les commandes iquo et irem (fournis-
sant quotient et reste dans la division euclidienne de deux entiers) et quo
et rem (fournissant quotient et reste dans la division euclidienne de deux
polynômes à coefficients entiers ou rationnels).
B. TP2 : ASSIGNATION/RÉASSIGNATION DE VARIABLES, BOUCLES (MAPLE) 117
(2) Dans l’exercice 9 de la feuille 1 (question 3), vous avez observé la syntaxe
du code PGCD.m (et du code div.m que ce code appelle) fournissant le
calcul du PGCD de deux entiers sous MATLAB. On rappelle ces deux codes
ici, tels qu’ils sont donnés dans le polycopié de cours :
function PGCD=PGCD(a,b);
x=a ;
y=b ;
while y>0
[q,r] = div(x,y);
if r==0
PGCD = y;
y = 0 ;
else
[q1,r1] = div(y,r);
x = r;
PGCD = x ;
y=r1 ;
end
end
En vous inspirant de ces modèles, élaborez et testez :
– un code sous Maple12 (utilisant les instructions iquo et irem) per-
mettant, au terme de son exécution, de calculer le PGCD de deux
entiers naturels a et b, avec b > 0 ;
– un code sous Maple12 (utilisant les instructions quo et rem) permet-
tant, au terme de son exécution, de calculer le PGCD (dans Q[X])
de deux polynômes A et B à coefficients rationnels, avec B 6= 0.
(3) Dans l’exercice 10 du TP1, vous avez observé la syntaxe du code bezout.m
permettant de manière récursive (ce code s’auto-appelle) de calculer, étant
donnés deux entiers relatifs a et b avec b 6= 0, une solution (u, v) de
l’identité de Bézout :
PGCD(|a|, |b|) = au + bv. (∗)
On rappelle ici ce code, tel qu’il figure dans le polycopié de cours :
fonction [PGCD,u,v]=bezout(a,b);
x=a ;
y= abs(b) ;
[q,r]=div(x,y);
if r==0
PGCD = y;
u=0 ;
v=1;
118 B. TP2 : ASSIGNATION/RÉASSIGNATION DE VARIABLES, BOUCLES (MAPLE)
else
[d,u1,v1]=bezout(y,r);
PGCD=d;
u=v1;
v=sign(b)*(u1- q*v1);
end
En vous inspirant de ce modèle, élaborez et testez :
– un code sous Maple12 (utilisant cette fois les instructions iquo et
irem) permettant, au terme de son exécution, de calculer, étant
donné un couple d’entiers relatifs (a, b) avec b 6= 0), le PGCD de
|a| et |b| ainsi qu’un couple d’entiers (u, v) solution de l’identité de
Bézout (∗) ;
– un code sous Maple12 (utilisant les instructions quo et rem) permet-
tant, au terme de son exécution, de calculer le PGCD (dans Q[X])
de deux polynômes A et B à coefficients rationnels, avec B 6= 0,
ainsi qu’un couple de polynômes U et V de Q[X] tels que
A(X)U (X) + B(X)V (X) = PGCD(A, B)
(le PGCD étant ici considéré à un facteur multiplicatif λ ∈ Q∗
près).
Sauvez votre travail dans votre dossier TPMaple12 en l’enregistrant sous le nom
TP2-Euclide.
Exercice B.6 (familiarisation avec les commandes eval et subs). Soit l’ex-
pression algébrique en trois variables :
P (X, Y, Z) = 3XY + 5XY Z 2 − 2XY 3 Z.
(1) En utilisant la commande eval (on en examinera avec ?eval la syntaxe),
déclarez sous Maple12 la fonction de la variable t :
t 7−→ P (tα , tβ , tγ ),
où α, β, γ sont trois nombres rationnels strictement positifs auxquels des
valeurs seront ensuite assignées.
(2) Affichez le graphe de la fonction :
h i
t ∈ [0, 1] 7→ P (x, y, z) cos(xy) .
x=t1/5 ,y=t1/8 ,z=t1/4
1. Attention : mettre ; derrière une instruction vous évite de voir l’affichage (c’est exactement
le contraire de ce qui se passe avec Maple12 !) ; ne rien mettre derrière vous fait au contraire courir
le risque de voir s’afficher un résultat, ce qui n’est en règle générale pas souhaitable, surtout lorsque
les variables en jeu sont de grands tableaux.
119
120 C. TP3 : FAMILIARISATION AVEC MATLAB, TRAVAIL SUR LES TABLEAUX
>> A = rand(10,8);
>> A
>> Adouble = [A A];
>> Adouble
>> B = randn(12,9);
>> B
>> Bdouble = [B;B];
>> Bdouble
>> AA = A(1:2:10,1:2:8);
>> AA
>> BB = B(1:3:8,2:2:9);
>> ONES = ones (7,5);
>> ONES
>> ZEROS = zeros (10,7);
>> ZEROS
>> EYE = eye (10,5);
>> EYE
>> AAA = flipud (A);
>> AAA
>> BBB = fliplr (B);
>> BBB
Analysez l’opération qu’exécute chacune de ces instructions. Qu’exécute
en particulier l’instruction Mat1=Mat(1:pas1:M,1:pas2:N) si Mat est un
tableau numérique à M lignes et N colonnes, pas1 et pas2 étant des entiers
respectivement entre 1 et M et 1 et N ?
(2) Construisez une matrice Apad à 16 lignes et 16 colonnes en bordant la
matrice A générée à la question (1) de manière symétrique (des deux côtés)
par des blocs de zéros 2. Faites ensuite la même chose avec cette fois des
blocs de 1.
(3) Créez une matrice à 19 lignes et 15 colonnes dont les lignes et les colonnes
sont celles de la matrice A générée à la question (1), mais séparées cette
fois entre elles par des lignes et des colonnes de zéros.
(4) Quelle est la méthode la plus judicieuse (au niveau du temps d’exécution)
pour déclarer sous MATLAB les matrices :
0 0 3 0 1 0 1 1 3 1 1 1
5 0 0 0 0 8 5 1 1 1 1 8
0 2 0 0 1 0 V =
U = 1 2 1 1 1 0 ?
0 0 1 0 8 0 1 1 1 1 8 1
1 0 0 0 0 5 1 1 1 1 1 5
Déclarez ces deux matrices de la manière que vous jugerez la plus efficace
possible (au niveau du temps d’exécution).
Sauvez votre session dans votre répertoire TPMATLAB (celui dans lequel vous êtes
précisément en train de travailler) avec l’instruction save TP3exo1. Vérifiez que
2. Pareille opération, fort utile dans la pratique, soit pour compléter une matrice rectangulaire
en une matrice carrée pour des calculs ultérieurs, soit pour s’accorder une marge de sécurité ,
telle un cadre , autour d’un tableau ou d’une image, s’appelle le zeropadding.
122 C. TP3 : FAMILIARISATION AVEC MATLAB, TRAVAIL SUR LES TABLEAUX
vous avez bien ainsi créé un fichier [Link] dans votre répertoire TPMATLAB.
Faites ensuite l’instruction clear all pour rendre à nouveau vierge tout votre
Workspace. Vérifiez que c’est bien le cas. Que se passe t’il si vous donnez les ins-
tructions :
>> load TP3exo1
>> who
Faites clear all à nouveau.
Exercice C.2 (décomposition de Haar d’une image). Voici encore un exercice
prétexte à la manipulation de tableaux sous MATLAB. Vous y verrez aussi l’étroite
relation entre tableaux de nombres (donc matrices) et images. Téléchargez depuis le
site [Link] le fichier (de
données) imageTP3exo2. Une fois ceci fait, chargez ce fichier dans votre Works-
pace en tapant sous la Command Window les instructions :
>> load imageTP3exo2
>> I=imageTP3exo2;
Grâce à la commande size, vérifiez que la variable ainsi déclarée I est un tableau à
256 lignes et 256 colonnes 3. Testez les valeurs de quelques entrées I(i,j) pour voir
qu’il s’agit en fait d’une matrice dont les entrées sont des nombres entiers positifs
(mais considérés ici comme des nombres flottants en double précision, comme la
réponse 1 à l’instruction isfloat(I) vous le prouvera, vérifiez le).
(1) Visualisez la matrice I de deux manières :
– en utilisant la visualisation en 3D suivant mesh(I) ;
– en utilisant la visualisation en 2D suivant image(I), imagesc(I)
(voire aussi imshow(I,[]) si vous êtes au CREMI, cf. plus loin).
Quelle est, pour ce type de matrice I, de ces deux visualisations, celle qui
est la plus adaptée ? Quel est l’effet de la commande
>> II = imrotate(I,angle);
(où angle désigne la valeur en flottant d’un angle entre 0 et 360 degrés
exprimé en flottant) ? 4. On utilisera l’une des routines du type image
mentionnées plus haut pour visualiser l’image II ainsi qu’une instruction
convenable (faire help axis) pour que l’affichage respecte la taille (ici
carrée) de l’image. Quelle est la difficulté liée à la réalisation mathématique
d’une telle transformation de matrice (pensez au passage entre le repérage
cartésien dans une grille de pixels et le repérage polaire) ? Que se passe
t-il en particulier lorsqu’un pixel (i, j) tourne d’un certain angle ? Tombe
t’on encore sur un pixel de la grille cartésienne ?
(2) Ouvrez un fichier .m vierge en utilisant l’onglet File du bandeau. Vous allez
essayer dans ce fichier de rédiger un code, dont la première instruction
(déclaration des entrées et sorties de la fonction) sera
function [RR,DH,DV,DO] = Haar (ipt)
3. Notez que 256 est une puissance de 2, ce qui n’est, on le verra plus tard, nullement un
hasard en ce qui concerne les formats d’images jpeg.
4. Cette instruction imrotate bien utile n’est malheureusement pas présente dans la bi-
bliothèque de Scilab.
C. TP3 : FAMILIARISATION AVEC MATLAB, TRAVAIL SUR LES TABLEAUX 123
qui est censé calculer, étant donnée une matrice de flottants ipt de taille
(2N , 2N ), les quatre matrices, de taille (2N −1 , 2N −1 ), dont les entrées sont
respectivement les nombres
ipt(2i − 1, 2j − 1) + ipt(2i − 1, 2j) + ipt(2i, 2j − 1) + ipt(2i, 2j)
RR(i, j) =
2
ipt(2i − 1, 2j − 1) + ipt(2i − 1, 2j) − ipt(2i, 2j − 1) − ipt(2i, 2j)
DH(i, j) =
2
ipt(2i − 1, 2j − 1) − ipt(2i − 1, 2j) + ipt(2i, 2j − 1) − ipt(2i, 2j)
DV(i, j) =
2
ipt(2i − 1, 2j − 1) − ipt(2i − 1, 2j) − ipt(2i, 2j − 1) + ipt(2i, 2j)
DO(i, j) =
2
Testez votre code en l’exécutant avec comme matrice ipt la matrice I :
>> [RR, DH, DV, DO] = Haar(I);
Affichez avec image (ou mieux imagesc) les quatre matrices RR, DH, DV,
DO. Pour le rendu des couleurs et du graphisme, vous pouvez jouer avec la
commande colormap (exemples les instructions colormap hot,colormap
jet, colormap pink, colormap jet, etc., faire help colormap pour dé-
cider). Pour un meilleur rendu graphique, vous pouvez aussi faire appel à
l’instruction
>> imshow(MATRICE,[ ])
Cette commande est dans la bibliothèque du Toolbox Image Processing
(disponible sous MATLAB au CREMI). Justifiez la terminologie : RR =
Résumé, DH = Détails Horizontaux, DV = Détails Verticaux, DO = Détails
Obliques. Recommencez le traitement de la matrice I cette fois sur la ma-
trice RR. Qu’observez vous sur les quatre images de taille (64, 64) ainsi
obtenues ? Poursuivez si vous voulez pour obtenir quatre images de taille
(32,32) en décomposant le nouveau résumé. La décomposition que vous
êtes en train d’effectuer ici s’appelle décomposition de Haar : elle permet
d’ isoler les structures cohérentes des détails (ou accidents )
d’une image digitalisée. Sauvez votre fichier .m sous le nom Haar si vous
ne l’avez pas encore fait et fermez le.
(3) Ouvrez un nouveau fichier vierge .m depuis l’onglet File du bandeau.
Rédigez dans ce fichier un code (correspondant encore à une fonction)
dont la première instruction sera
function II = HaarInverseAux (I1,I2,I3,I4)
qui, étant donnée quatre matrices I1,I2,I3,I4 de taille (2N −1 , 2N −1 ),
fabrique une matrice II de taille (2N , 2N ) telle que
II(2i − 1, 2j − 1) = I1(i, j), II(2i − 1, 2j) = I2(i, j)
II(2i, 2j − 1) = I3(i, j), II(2i, 2j) = I4(i, j)
pour toute paire d’entiers (i, j) entre 1 et 2N −1 . Sauvez ce nouveau fichier
.m sous le nom HaarInverseAux. Vous pourrez être amenés à utiliser ce
code comme code auxiliaire par la suite.
(4) Ouvrez un nouveau fichier .m où vous réaliserez un code (toujours une
fonction) dont la première instruction sera cette fois
124 C. TP3 : FAMILIARISATION AVEC MATLAB, TRAVAIL SUR LES TABLEAUX
5. Ce devrait en principe être 0 pour la machine lorsque Nb<N ; si ce n’est pas le cas, ceci
signifie que l’algorithme n’a pu aboutir au terme du nombre garde-fou N d’itérations imparti a
priori.
ANNEXE D
> op(8,{5,6,7,n,8,n,9,n+1});
> op(9,{5,6,7,[n,8,n,9,n+1]});
> op(3..4,{5,6,7,[n,8,n,9,n+1]});
> op(4..6,{5,6,7,n,8,n,9,n+1});
> op(4..6,{5,6,7,[n,8,n,9,n+1]});
> op([5,6,7,n,{8,n,9,n+1}]);
> nops([5,6,7,n,{8,n,9,n+1}]);
Pour appliquer une fonction (disons g, préalablement déclarée) à une suite, une
liste, ou un ensemble (et donc calculer la suite, la liste ou l’ensemble constituées
des images par g des éléments de la suite, de la liste ou bien de l’ensemble initial),
on utilise l’instruction map. Exemple :
> g:= x-> x^3 - sin(x);
> map(g,S);
> map(g,L);
> map(g,E);
Recommencez après avoir assigné une valeur numérique à n. Entrainez vous.
Dans ce TP, vous allez apprendre à rédiger des programmes autonomes. Sup-
posons que l’on souhaite rédiger une procedure sous Maple12, c’est-à-dire créer
une fonction Maple qui, étant donnée un certain nombres d’arguments (ou encore
données en input, par exemple argument1,argument2,...,argumentk), ren-
Pourquoi est-il essentiel de faire figurer la boucle if n=0 then else ...
end if; dans le synopsis ? Déduisez de cela la construction d’une procé-
dure récursive
fonction1TP4 := proc(raison,init,n)
qui, étant données les trois arguments raison, init, n, génère la suite
arithmétique de premier terme init et de raison raison.
(2) Sur le modèle établi à la question 1, rédigez une procédure récursive
FibonacciRec:=proc(n)
permettant de calculer le n-ième nombre de Fibonacci. On rappelle que la
suite des nombres de Fibonacci est la suite (Fn )n≥0 telle que F0 = F1 = 1
et
Fn = Fn−2 + Fn−1 ∀ n ≥ 2.
Adaptez ce code pour l’insérer dans une procédure récursive
fonction2TP4 :=proc(a,b,init0,init1,n)
qui, étant donnés les cinq arguments a,b,init1,init2,n, génère la suite
(un )n≥0 initiée à u0 = init0 et u1 = init1 et régie ensuite par la
récurrence linéaire à deux pas :
un = a un−1 + b un−2 ∀ n ≥ 2.
(3) En utilisant la commande time (combinée avec la commande seq per-
mettant de générer des suites ou des listes), générez la liste L des temps
de calcul CPU pris par l’exécution de FibonacciRec (i), i=20..35. En
utilisant listplot (sous with (plots);), affichez (en fonction de l’indice
i) la liste obtenue en transformant L par l’application log. Que constatez
vous ?
(4) Pour chaque entier n supérieur ou égal à 2, soit cn le nombre d’appels
au programme FibonacciRec impliqués dans FibonacciRec(n). Montrez
que la suite (cn )n≥2 vérifie la relation cn = cn−1 + cn−2 + 1, puis déduisez
en que cn = 2Fn − 1. Lorsque n tend vers +∞, prouvez en utilisant cette
fois la relation matricielle
Fn 1 1 Fn−1
=
Fn−1 1 0 Fn−2
√
et votre cours d’Algèbre 3 que Fn ∼ ((1 + 5)/2)n . Comparez la√pente
visible sur le graphe obtenu précédemment avec le nombre d’or (1+ 5)/2.
Pouviez vous prévoir le résultat observé ainsi ?
(5) Modifiez la syntaxe de la procédure FibonacciRec en utilisant l’instruc-
tion additionnelle option remember, sous la forme
FibonacciRec := proc(n) option remember;
Testez là maintenant avec des valeurs de n supérieures à 35. Que constatez
vous ? Modifiez la procédure obtenue en une procédure FibonacciRecbis
de manière à ce que se trouvent affichées en sortie toutes les valeurs in-
termédiaires F0 , ..., Fn (au lieu de simplement la valeur finale Fn ). Testez
là pour n = 30 et comparez avec les temps de calculs impliqés dans l’uti-
lisation de FibonacciRec aux questions 3 et 4.
130 D. TP4 : PROCÉDURES SOUS MAPLE, REPRÉSENTATION GRAPHIQUE
Sauvez votre feuille de travail comme TP4exo1 dans votre répertoire TPMaple12
puis fermer cette feuille de travail.
Sauvez votre travail cette fois comme le fichier TP4exo3 (toujours dans le répertoire
TPMaple12) et fermez ensuite cette feuille de travail.
Exercice D.4 (une approximation de π par des suites adjacentes). Ouvrez
une nouvelle feuille de travail.
(1) Réalisez une procédure approxPi1TP4:=proc(n) (basée sur une boucle
DO) générant simultanément les deux suites (un )n≥0 et (vn )n≥0 de nombres
strictement positifs initiées respectivement par u0 = 1., v0 = 2. (attention
de déclarer ces entrées comme des nombres flottants !) et obéissant au
couple de relations inductives :
un + vn
un+1 =
√ 2
vn+1 = un+1 vn .
Lancez cette procédure sous par exemple Digits:=20 avec n=50. Que
constatez vous ? On admet que les deux suites (u√ n )n≥0 et (vn )n≥0 sont
adjacentes et ont pour limite commune le nombre 27/π. Modifiez votre
procédure en une procédure approxPi2TP4:=proc(n) qui fournisse √ l’ap-
proximation de π déduite de l’égalité approchée (un + vn )/2 ' 27/π.
(2) Modifiez la procédure approxPi2TP4 en une procédure
approxPi3TP4:=proc(epsilon,n)
qui, étant donné un seuil epsilon fixé, retourne
√ l’approximation de π
déduite de l’égalité approchée (un + vn )/2 ' 27/π dès que la condition
|un − vn |
≤
un + vn
est satisfaite (et affiche le nombre i≤n d’itérations de la boucle DO néces-
saires).
(3) Réalisez une procédure inductive approxPiRecTP4:=proc(n) générant si-
multanément les deux suites (un )n≥0 et (vn )n≥0 de la question 1, ce de
manière exacte cette fois (on prend u0 = 1 et v0 = 2, déclarés cette fois
comme rationnels et non plus comme flottants). Comparez, pour n=12
(attention surtout à ne pas prendre de valeur plus grande, car le temps de
calcul risque d’exploser !) le temps d’exécution de cette nouvelle procédure
approxPiRecTP4(12) avec celui de la procédure approxPi1TP4(12) (dans
laquelle u0 et u1 sont déclarés comme les entiers 1 et 2 et non plus comme
les flottants 1.0 et 2.0). Que constatez vous ?
Sauvez votre feuille de travail comme TP4exo4 (dans TPMaple12) et fermez cette
feuille.
Exercice D.5 (un modèle de courbe fractale : le flocon de Von Koch). Les
courbes planes qui ont la propriété d’autosimilarité, c’est-à-dire qui, regardées à la
loupe avec n’importe quel grossissement, reproduisent à toutes les échelles le même
motif, se prêtent (en ce qui concerne leur tracé) à la réalisation de procédures
récursives. Le flocon de Von Koch , introduit vers 1900 par le mathématicien
suédois Helge von Koch (1870-1924), en est une illustration ; regardez par curiosité
avant de commencer le site dédié :
[Link]
D. TP4 : PROCÉDURES SOUS MAPLE, REPRÉSENTATION GRAPHIQUE 133
(1) Soient a et b deux nombres complexes, définissant un segment (noté [a, b])
du plan complexe. On introduit les trois nombres complexes :
b−a 2a + b
u(a, b) := a + =
3 3
b−a 2b + a
v(a, b) := b − =
3 3
w(a, b) := u(a, b) + eiπ/3 (v(a, b) − u(a, b))
et l’on note Φ l’application qui au segment [a, b] (d’origine a et d’extrémité
b) associe la ligne brisée Φ([a, b]) joignant (dans cet ordre) les quatre points
a, u(a, b), w(a, b), v(a, b), b. Si L est une ligne brisée de sommets d’affixes
a1 , a2 , ..., aN , la ligne brisée Φ(L) désigne par extension la ligne brisée
obtenue en mettant bout à bout les diverses lignes brisées Φ([aj , aj+1 ])
pour j = 1, ..., N − 1. Exprimez en fonction de N (nombre de sommets de
la ligne brisée L) le nombre de sommets de la ligne brisée Φ(L). Réalisez
une procédure
VonKoch1TP4:=proc(L)
qui, étant donnée une liste L de nombres complexes, considérée comme
la liste [a1,...,aN] des affixes des sommets successifs d’une ligne brisée
L, renvoie (dans l’ordre, et de manière approchée) la liste des affixes des
sommets de la ligne brisée Φ(L).
(2) En utilisant la procédure VonKoch1TP4, réalisez une nouvelle procédure
VonKoch2TP4:=proc(n)
qui renvoie (dans l’ordre) la liste des affixes des sommets de la ligne brisée
Φ[n] ([0, 1]), où Φ[n] désigne l’application Φ itérée n fois. Vu que le nombre
de sommets de la ligne brisée Φ[n] ([0, 1]) croit vers l’infini comme 2 × 4n
(dites pourquoi), vous ne testerez cette procédure VonKoch2TP4 que sur
des valeurs de n entre 0 et 5. En utilisant les instructions
> with (plots);
> L:=VonKoch2TP4(n):
> complexplot(L,style = line, scaling = constrained);
affichez la ligne brisée obtenue en itérant n fois Φ à partir du segment
[0, 1]. Veillez à ne prendre comme valeurs de n que des valeurs entre 1 et 7
(notez que 2 × 47 = 32768 est déjà très grand !). Que se passe-t-il si vous
remplacez la dernière instruction par :
> complexplot(L,style=line);
La courbe ainsi obtenue est une approximation du flocon de Von Koch ; ce
flocon de Von Koch est un exemple de courbe fractale, continue mais ne
présentant de tangente en aucun point 2 ; c’est en fait la limite uniforme
des courbes ainsi tracées – en fonction de n – lorsque n tend vers l’infini).
2. Vous verrez plus tard dans votre cursus mathématique que pareille courbe fractale n’est
pas rectifiable , et que l’on ne peut donc pas parler de longueur du flocon de Von Koch.
La nature fourmille de tels modèles fractaux (en relation avec ce que l’on appelle le chaos
dynamique ) : pensez par exemple à la découpe d’une côte volcanique telle celle du Groendland,
à celle de la ligne d’arête des aiguilles de Chamonix ... Vous trouverez aussi beaucoup de modèles
(mathématiques cette fois) d’images fractales sur le web (courbes de Mandelbrojt, frontières de
domaines de Julia, etc.).
134 D. TP4 : PROCÉDURES SOUS MAPLE, REPRÉSENTATION GRAPHIQUE
3. Notéz que 2 × 48 = 131072, ce qui est vraiment très grand pour la taille d’une liste de
flottants sous Maple12 !
D. TP4 : PROCÉDURES SOUS MAPLE, REPRÉSENTATION GRAPHIQUE 135
qui, étant donnée une liste d’affixes L (figurant la liste des sommets or-
donnés d’une ligne brisée L) renvoie la liste (ordonnée) des affixes des
sommets de la ligne brisée Ψ[n] (L) dès l’instant où le pas de cette ligne
brisée devient strictement inférieur au seuil epsilon.
(3) En utilisant complexplot (après avoir déclaré l’environnement plots par
l’instruction with(plots);), représentez la liste L obtenue via
> L:=Levy2TP4(.03,[I,1,I]):
Sauvez votre feuille de travail comme TP4exo6 dans votre répertoire TPMaple12,
puis fermez cette feuille.
ANNEXE E
(formule (3.8) du polycopié, page 50), rédigez dans ce fichier une fonction :
function LXX = LagrangeNewton (X,Y,XX)
qui, étant donnée une liste de nombres complexes distincts X (de longueur
N+1) et une liste de nombres complexes Y de même longueur, évalue le
polynôme de Lagrange Z 7→ Lagrange[X;Y](Z) aux points de la liste
1. Ceci sera repris dans un TP ultérieur sous Maple12, cette fois dans le cadre du calcul
arithmétique sans pertes, sous l’angle du calcul symbolique et de la cryptograhie.
137
138 E. TP5 : LES FONCTIONS SOUS MATLAB ET L’INTERPOLATION
>> plot(XX,f(XX));
(utilisez la routine conj qui conjugue les entrées des vecteurs ou des ma-
trices). Validez ce code sur un exemple, par exemple :
>> X=rand(16,1) + i*rand(16,1);
>> Y=fft(X,16);
>> XX=InverseFFT(Y,16);
>> max(abs(X-XX))
Que constatez vous ? Pourquoi ne trouvez vous pas 0 comme vous vous
y attendriez ? Ne perdez pas de vue que vous êtes ici en train de faire
du calcul scientifique, et non du calcul symbolique. On reviendra dans
un TP ultérieur sous Maple12 cette fois sur la fft discrète, dans le cadre
(arithmétique) du calcul symbolique cette fois. Sauvez votre fichier (appelé
InverseFFT) comme un fichier .m dans votre répertoire TPMATLAB pour le
conserver dans vos archives, tout en sachant que la routine
>> X=ifft(X,N);
(déjà implémentée dans le noyau du logiciel) correspond de fait à la même
fonction.
(2) Ouvrez un fichier .m vierge. Rédigez dans ce fichier une procédure
function P1P2 = ProduitPolynomes (P1,P2)
qui, étant donnés deux polynômes P1 et P2 déclarés sous forme de listes
comme
>> P1 = [P1(1) ... P1(N1)];
>> P2 = [P2(1) ... P2(N2)];
142 E. TP5 : LES FONCTIONS SOUS MATLAB ET L’INTERPOLATION
lorsque
P1(X) = P1(1) X N1−1 + · · · + P1(N1-1) X + P1(N1)
P2(X) = P2(1) X N2−1 + · · · + P2(N1-1) X + P2(N2),
retourne (sous forme d’une liste de longueur N1+N2-1) les coefficients du
polynôme P1*P2, les monômes de ce polynôme étant rangés dans l’ordre
décroissant. La première étape de la procédure est de déterminer le pre-
mier entier p tel que 2p > N1+N2-1 (l’entier le plus proche d’un nombre
réel flottant x lorsque l’on va à droite de ce nombre est ceil(x) ; c’est
floor(x) si l’on va à gauche). Utilisez ensuite les routines fft(.,N) et
ifft(.,N) comme indiqué dans le cours, avec précisément N = 2p , ce
après avoir si nécessaire complété les listes P1 et P2 par des zéros. Sau-
vez le fichier ProduitPolynomes.m dans votre répertoire TPMATLAB après
l’avoir testé (et validé) sur des polynômes de bas degré (inférieur ou égal
à 10).
(3) La commande MATLAB permettant de déterminer l’écriture binaire d’un
entier M (bien sûr inférieur à 252 ) est
>> P = dec2bin(M);
mais, attention, la sortie P est ici au format char, pas au format numérique
double précision double. Ce format double est pourtant nécessaire ici car
les calculs de fft impliquent la matrice de nombres fottants
[WNkj ]0≤j,k≤N−1 .
Il faut donc convertir P au format double pour travailler numériquement ;
en faisant un petit test, vous verrez que le caractère 0 correspond dans
cette conversion à un certain entier positif P=double(’0’), tandis que
1 correspond (heureusement) à P+1=double(’1’), ce qu’il convient donc
de prendre en compte pour les conversions. À vous de trouver ces deux
nombres double(’0’) et double(’1’), car vous voulez, vous, 0 et 1
comme nombres. Sur un nouveau fichier vierge .m, rédigez une procédure
function M1M2=ProduitEntiers(M1,M2)
qui, étant donnés deux entiers naturels tous deux inférieurs à 226 (cela
vaut mieux, dites pourquoi !), utilise les routines dec2bin, fft(.,64),
ifft(.,64) pour calculer leur produit M1*M2. Pensez à associer à tout
entier naturel M (inférieur à 252 ) le polynôme
>> PolM = double(dec2bin(M))-double(’0’) ;
(pour avoir des coefficients numériquement égaux à 0 et 1 et non plus
cette fois à double(’0’) et double(’1’) comme le donne la conversion
numérique du format char au format double), puis ensuite à évaluer le
polynôme PolM1*PolM2 (calculé comme à la question 2) en x=2.
Exercice E.4 (L’algorithme de Cooley-Tukey). La routine sous MATLAB
>> Y=fft(X,N);
(lorsque N est une puissance de 2) correspond à l’implémentation de l’algorithme de
Cooley-Tukey (cf. le polycopié de cours pages 44-45). Vous pouvez aussi consulter
pour plus de détails les sites wikipedia
E. TP5 : LES FONCTIONS SOUS MATLAB ET L’INTERPOLATION 143
[Link]
[Link]
(le site en anglais est plus riche). Le but de cet exercice (prétexte, comme dans
le TP4, à votre familiarisation avec le principe de récursivité) est de réaliser une
routine réalisant la même opération, ce afin de comprendre le schéma récursif sur
lequel se fonde la procédure. Ouvrez un fichier .m vierge sur lequel vous allez rédiger
une procédure récursive :
function Y=CooleyTukey(X,p)
transformant, si p désigne un entier supérieur ou égal à 1 et X un vecteur colonne de
longueur 2p , le vecteur X en son image par la multiplication à gauche par la matrice
h i
W2jk
p
p−1
, où W2p := exp(−2iπ/2p ). (†)
0≤j,k≤2
qui sont, elles aussi, des fonctions oscillantes, donc de même nature que la fonction
à interpoler f. Vous verrez plus tard que c’est le principe de ce que l’on appelle
faire l’analyse de Fourier de l’information fournie par f. Soit N = 2p .
(1) Vérifiez, si p est un entier naturel non nul donné, que le système de 2p
équations à 2p inconnues uk , k = −2p−1 , ..., 2p−1 − 1,
h X−1
2p−1 i
uk e2iπkθ = f (j/2p ), j = 0, ..., 2p − 1 (∗∗)
θ=j/2p
k=−2p−1
>> U= InterpolTrigo1(f,p);
Testez votre routine sur un polynôme trigonométrique f1 (mais non 1-
périodique) déclaré en ligne
>> f5= inline(’4*sin(pi/2*x) + cos(5*x) - 2*sin(3*x)’,’x’);
suivant les instructions :
>> XX=0:1/2000:1;
>> PXXT=InterpolTrigo2(f5,p,XX);
>> plot(XX,f5(XX),’r’)
>> hold
>> plot(XX,PXXT)
Prenez pour cela des valeurs de p entre 6 (2p = 64) et 10 (2p = 1024).
Qu’observez vous lorsque p augmente concernant la qualité de l’approxi-
mation de f par son polynôme trigonométrique interpolant ? Recommen-
cez avec cette fois la fonction :
>> f6= inline(’4*sin(pi/2*x) + cos(5*x) - 2*sin(3*x)’,’x’);
Après avoir évalué f(1)-f(0), essayez de d’expliquer ce qui crée le phéno-
mène au bords de [0, 1] dans le cas du premier exemple, et ne semble plus
le créer (en tout cas de manière aussi nette) ici. Ce phénomène s’appelle,
lorsqu’il se produit phénomène de Gibbs ; c’est le pendant du phénomène
de Runge observé dans le cadre de l’interpolation polynomiale par le po-
lynôme de Lagrange (Exercice 1). Voir par exemple le site (succint) sur
wikipedia :
[Link]
Sauvez votre fichier InterpolTrigo2.m dans votre répertoire TPMATLAB.
Testez aussi votre programme avec une fonction polynomiale :
>> f7 = inline(’x.^7+3*x.^5-2*x+1’,’x’);
(c’est la fonction f1 de l’exercice 2). Est-il préférable dans ce dernier cas
d’utiliser l’interpolation par des polynômes trigonométriques plutôt que
l’interpolation de Lagrange par des polynômes ?
(3) Pour mieux vous convaincre du phénomène de Gibbs, déclarez la fonction
f8 définie par
>> f8=inline(’(1/2)*(sign ((x-1/3).*(2/3-x))+1)’,’x’);
Affichez le graphe de f8 :
>> XX=0:1/2000:1;
>> plot(XX,f8(XX));
Sur le même graphe, affichez (en utilisant la commande hold) avec diverses
couleurs les graphes successifs des divers PXXT :
>> PXXT = InterpolTrigo2(f8,p,XX);
>> hold
>> plot(XX,PXXT,’color’);
146 E. TP5 : LES FONCTIONS SOUS MATLAB ET L’INTERPOLATION
>> Q:=diff(P,X);
>> R:=resultant(P,Q,X);
>> degree(R);
calculez les discriminants ordinaires ∆3 , ∆4 , ∆5 , ∆6 , ainsi que le degré
total de ces trois polynômes. Essayez de deviner une formule donnant le
degré de ∆d , puis justifiez cette formule en vous souvenant que ∆d doit
se représenter comme un certain déterminant de Sylvester ; lequel ?
Sauvez votre feuille de travail comme [Link] dans votre répertoire TPMaple12,
puis fermez là.
Exercice E.7 (intersection de deux courbes planes). Une courbe plane définie
dans C × C (espace vectoriel de dimension 2 sur C) comme le lieu des zéros d’un
polynôme en deux variables de degré total égal à 3 est appelée cubique (tandis
qu’elle est appelée conique si elle est définie comme le lieu des zéros d’un polynôme
de degré total égal à 2). Ouvrez une nouvelle feuille de travail sous Maple12.
(1) Vérifiez que les sous-ensembles de C × C définis par
n o
Γ1 := (x, y) ∈ C × C ; xy(y − x) − y 2 + 6x − y − 6 = 0
n o
Γ2 := (x, y) ∈ C × C ; (y − 2x)2 (y − 3x) − (x + y)2 − x − y − 1 = 0
sont bien des cubiques et que le seul zéro commun des parties homogènes
de plus haut degré des polynômes
P:= X*Y*(Y-X) - Y^2 + 6*X - Y - 6 ;
Q:= (Y-2*X)^2*(Y-3*X) - (X+Y)^2 - X -Y -1 =0;
définissant ces courbes est l’origine (0, 0) de C2 . On admettra que dans
ce cas (c’est un important théorème de géométrie algébrique attribué à
Etienne Bézout) les deux cubiques Γ1 et Γ2 se coupent en exactement
9 = 3 × 3 (le produit des degrés totaux de leurs équations définissantes
P=0 et Q=0) points de C × C (les points d’intersection des deux courbes
pouvant être éventuellement multiples).
(2) Utilisez l’affichage graphique 3
>> with(plots);
>> implicitplot([P=0,Q=0],X=-10..10,Y=-10..10,
color=[red,blue],gridrefine=2);
pour tracer les intersections des deux courbes Γ1 et Γ2 avec le monde
réel R2 . Qu’observez vous ? Les deux courbes réelles se coupent-elles et
en combien de points visibles ? Se coupent t’elles transversalement ? Que
peut-on dire des autres points d’intersection des deux courbes Γ1 et Γ2
(puisque l’on sait a priori qu’il doit y en avoir 9 dans C2 ) ?
(3) En utilisant le principe de l’élimination (et donc la commande resultant)
– d’abord pour éliminer la variable X entre les équations P(X,Y)=0 et
Q(X,Y)=0,
– ensuite pour éliminer la variable Y entre ces deux mêmes équations,
3. Prenez soin de décortiquer au préalable le sens des options de la routine implicitplot
(affichant une courbe dans R2 donnée par son équation cartésienne) après avoir fait ?implicitplot
dans votre feuille de travail pour analyser toutes les potentialités de cette commande implicitplot.
TP6 : LE PRINCIPE DE L’ÉLIMINATION ALGÉBRIQUE SOUS MAPLE 149
puis la routine
>> zeros :=[fsolve(R,T,complex)];
(qui donne sous forme de liste les d zéros complexes, calculés de manière
approchée, d’un polynôme R ∈ C[T ] de degré d) déterminer, concernant
les points d’intersection (x, y) des deux courbes Γ1 et Γ2 dans C × C :
– 9 valeurs possibles pour x, parmi lesquelles on observe que se trouve
une et une seule valeur réelle x0 ;
– 9 valeurs possibles pour y, parmi lesquelles on observe que se trouve
une et une seule valeur réelle y0 .
Déduisez-en les valeurs approchées des coordonnées de l’unique point
de R2 appartenant aux deux cubiques Γ1 et Γ2 . Le résultat est il bien
conforme à ce que vous avez observé graphiquement à la question 2 ?
(4) Rédigez une procédure
testzeros = proc(epsilon,P,Q,LX,LY)
qui, étant donnés deux polynômes en deux variables P et Q et deux listes
de nombres complexes LX, LY (déclarés en flottants), renvoie p la liste des
couples (L[i],L[j]) en lequels l’évaluation de la fonction |P |2 + |Q|2
est strictement inférieure au seuil . En prenant convenable (essayez
10−k , k variant entre 3 et 8), exécutez ce code testzeros avec comme
autres entrées P, Q, et les deux listes LX, LY de 9 nombres complexes
construites à la question 3, pour obtenir la liste des valeurs approchées des
9 couples de nombres complexes correspondant aux positions des 9 points
d’intersection des cubiques Γ1 et Γ2 cette fois dans C × C. Retrouvez vous
bien le point réel observé à la question 2 et calculé de manière approché
à la question 3 ?
(5) On considère, à la place de la cubique Γ2 , la cubique Γ3 définie par
n o
Γ01 := (x, y) ∈ C × C ; xy(x − y) − x2 + 6y − x − 6 = 0 .
qui fournisse, étant donnés en input deux polynômes P et Q en une variable X (de
degrés respectifs p et q, à coefficients pouvant dépendre de manière polynomiale ou
rationnelle d’autres variables Y,Z...), fournisse en output le déterminant de la ma-
trice de Sylvester (de taille (p + q, p + q), voir (3.14) dans la section 3.4 du polycopié
de cours) de P et Q, considérés comme des polynômes en X. Ouvrez une nouvelle
feuille de travail sous Maple12.
(1) Chargez les bibliothèques suivantes :
>> with(PolynomialTools);
>> with(LinearAlgebra);
et étudiez ce que font les routines :
>> ?CoefficientList;
>> ?Determinant;
ainsi que la syntaxe régissant leur utilisation. Consultez aussi
>> ?Matrix;
>> ?MVassignment;
pour voir comment modifier une matrice S de taille D×D donnée composée
initialement de zéros (déclarée par S:=Matrix(D)).
(2) Rédigez une procédure
Sylvester := proc(P,Q,X)
fournissant, lorsque P et Q sont deux polynômes non nuls de degrés stric-
tement positifs en la variable X, la matrice de Sylvester (i.e la matrice
(3.14) des notes de cours).
(3) Rédigez une procédure
resultantTP := proc(P,Q,X)
fournissant, lorsque P et Q sont deux polynômes non nuls de degrés stric-
tement positif en la variable X, leur résultant. Recalculez les déterminants
ordinaires ∆3 et ∆4 de l’exercice 1 (question 3) avec cette fois votre rou-
tine resultantTP et comparez avec les résultats obtenus à l’exercice 1
(question 3) avec la routine (intégrée au logiciel) resultant.
(4) Étant donnés deux polynômes P et Q de degrés respectifs p > 0 et q > 0,
montrez que l’on ne change pas le déterminant de la matrice de Sylvester
(matrice (3.14) dans les notes de cours) en remplacant le dernier vecteur
colonne par le vecteur colonne :
q−1
X ∗ P(X)
..
.
X ∗ P(X)
P(X)
p−1 (∗)
X ∗ Q(X)
..
.
X ∗ Q(X)
Q(X)
TP6 : LE PRINCIPE DE L’ÉLIMINATION ALGÉBRIQUE SOUS MAPLE 151
(pensez pour cela à ajouter au dernier vecteur colonne la somme de tous les
vecteurs colonnes précédents multipliés par Xp+q−1 , Xp+q−2 , ...). Réalisez
une procédure
resultantTPmodif:=proc(P,Q,X,S,T);
qui, étant donnés deux polynômes P et Q en une variable X (de degrés stric-
tement positifs p et q, les coefficients de ces polynômes pouvant dépendre
de manière polynomiale ou rationnelle d’autres variables Y,Z...), calcule
le déterminant de la matrice Sylvester(P,Q,X) modifiée par le fait que
la dernière colonne y ait été remplacée par
q−1
X ∗S
..
.
X∗S
S
p−1 , (∗∗)
X ∗ T
..
.
X∗T
T
où S et T sont deux nouvelles variables. La sortie de cette procédure doit
donc être un polynôme en les variables X, S, T et les coefficients des po-
lynômes P et Q. Vous pourrez au préalable consulter l’aide
>>?Vector;
pour voir comment déclarer le vecteur colonne (**). Calculez, étant donnés
deux polynômes P et Q comme ci-dessus
>> R:= resultantTPmodif(P,Q,X,S,T);
>> resultant(P,Q,X) - simplify(P*diff(R,S)+Q*diff(R,T));
Expliquez d’où ce résultat provient. Si le résultant de P et Q par rapport à
X est non nul, pourquoi obtient-on ainsi une identité de Bézout 1=A*P+B*Q
avec deg A < q et deg B < q ?
Sauvez votre feuille de travail comme [Link] dans votre répertoire TPMaple12,
puis fermez la.
Exercice E.9 (Intégration symbolique de fractions rationnelles et élimination).
Les méthodes reposant sur l’élimination algébrique ainsi que sur la division eucli-
dienne (algorithme de Bézout étendu) jouent un rôle clef dans le noyau de logiciels
de calcul symbolique tels Maple12 lorsqu’il s’agit d’effectuer l’intégration symbo-
lique des fractions rationnelles. C’est ceci que cet exercice se propose d’illustrer.
Ouvrez une nouvelle feuille de travail sous Maple12.
(1) Soit Q un polynôme en une variable à coefficients entiers, de degré stric-
tement plus grand que 1, n’ayant aucune racine double dans C. Pourquoi
le résultant (par rapport à X) de Q et Q0 := dQ/dX est-il non nul ?
Pourquoi existe-t-il bien dans ce cas deux polynômes U et V à coefficients
rationnels, avec de plus deg(U ) ≤ deg Q − 2 et deg V ≤ deg Q − 1 tels que
l’on ait l’identité Q(X)U (X) + Q0 (X)V (X) ≡ 1 ?
(2) Regardez la syntaxe de la routine gcdex :
152 TP6 : LE PRINCIPE DE L’ÉLIMINATION ALGÉBRIQUE SOUS MAPLE
>> gcdex(P1,P2,X,’U’,’V’);
>> U;
>> V;
>> simplify(U*P1+V*P2);
Prenez un polynôme Q de degré 3 sans racine double et à coefficients
entiers, par exemple
>> Q:= X^3 - 2*X^2 + 3*X+5;
et calculez les polynômes U et V tels que U*Q+V*diff(Q,X)=1 donnés par
la routine gcdex. Calculez d’autre part
>> resultant(Q,diff(Q,X),X);
Qu’observez vous ? Est-ce une surprise et avez vous une explication (pen-
sez r̀egarder ici les notes de cours, section 3.4) ?
(3) Vérifiez, si P et Q sont deux polynômes (avec deg Q > 0 et Q sans racine
double dans C) que, si U et V sont tels que l’on ait l’identité de Bézout
1 = U Q + V Q0 , alors, pour tout entier n > 1,
Z X
dt 1 V (X)P (X)
P (t) = +
(Q(t))n 1 − n ((Q(X))n−1
Z X
1 d dt
+ P (t)U (t) + [P V ](t) .
n − 1 dt (Q(t))n−1
(4) En utilisant la procédure gcdex (algorithme d’Euclide étendu résolvant
l’identité de Bézout, voir la question 2) et la procédure int
>> F:=int(P,X);
appliquée aux polynômes (qui fournit une primitive d’un polynme en X
et qu’il vous serait aisé de construire vous-même), rédigez une procédure
récursive
HermiteItere:=proc(P,Q,X,n)
qui, étant donnés deux polynômes P, Q en la variable X (avec Q de degré
strictement positif et sans racine double) à coefficients entiers retourne
une liste [R,H], où R est un élément de Q(X) de dénominateur Qn−1 et H
un élément de Q(X) de degré strictement inférieur au degré de Q tels que
l’on ait : Z X Z X
P(t) H(t)
n
dt = R(X) + dt.
(Q(t)) Q(t)
(5) Soient H et Q deux polynômes à coefficients rationnels tels que 0 ≤ deg H <
deg Q = d et
>> RX := resultant(H-Y*diff(Q,X),Q,X);
Pourquoi RX est-il un polynôme en Y à coefficients rationnels ? Quel est
son degré ? Soient α1 , ..., αd les d racines (complexes, mais algébriques,
toutes simples par hypothèses) de Q. Montrez (en revenant à la définition
du résultant) que :
H(αl )
RX(z) = 0 ⇐⇒ ∃ l ∈ {1, ..., d} t.q. z = 0 .
Q (αl )
Montrez que, si z est une racine de RX, alors le polynôme
TP6 : LE PRINCIPE DE L’ÉLIMINATION ALGÉBRIQUE SOUS MAPLE 153
(la seconde routine est un peu compliquée, mais je n’ai pas trouvé mieux
pour regrouper les facteurs avec radicaux figurant au dénominateur de la
fraction rationnelle Derj initiale). Calculez aussi le résultant RXj et sa
factorisation avec les commandes
>>RXjbis:=Split(RXj,Y);
>>convert(RXjbis,radical);
En conclusion de ces exemples, on voit que la procédure de Rothstein-
Trager s’avère d’autant plus intéressante que le résultant RX a des zéros
multiples. Le nombre de zéros distincts de RX peut s’avérer alors nettement
plus petit que le degré du polynôme Q (qui est aussi le degré en Y du
résultant RX, voir la question 5). De plus, si par chance les racines distinctes
de RX sont rationnelles (comme dans l’exemple 3), le calcul de la liste des
PGCD Rz (pour z balayant la liste des zéros distincts de RX) ne fait pas
apparaitre de factorisation intempestive impliquant des radicaux. Si les
zéros distincts de RX sont calculables, mais complexes, on notera cependant
que le regroupement des logarithmes conjugués via la formule formelle
ln X-(a+I*b) − ln X-(a-I*b) = 2 ∗ I ∗ arctan((X-a)/b)
si a ∈ R et b ∈ R∗ n’est pas opéré ici. Il conviendrait de l’effectuer
pour avoir une forme simplifiée (et réelle) de l’expression de la primitive
obtenue, ce que nous ne ferons pas faute de temps. Vous pouvez ajouter
cette opération supplémentaire de concaténation à titre d’exercice.
Sauvez votre feuille de travail comme [Link] dans votre répertoire TPMaple12,
puis fermez la.
ANNEXE F
1. Cela ne donne pas a priori une majoration par epsilon de l’erreur commise entre le rayon
spectral et sa valeur approchée (au terme de Niter itérations), mais seulement un contrôle de
155
156 F. TP7 : LE THÉORÈME DU POINT FIXE EN ACTION SOUS MATLAB
>> X1 = [1;1;1;1];
>> X2 = [1;1;0;0];
>> X3 = [1;0;0;0];
>> X4 = rand(4,1);
Quelle valeur de Niter trouvez vous dans chacun de ces cinq cas ? Vérifiez
que la valeur de r obtenue alors est bien en accord avec le résultat fourni
à la question 1 par la commande eig(A) donnant les quatre valeurs
propres (dans ce cas réelles) de la matrice A. Recommencez avec cette
fois epsilon = 10−6 .
(3) On suppose que A est une matrice réelle et que le rayon spectral r(A) est
la valeur absolue d’une valeur propre réelle simple de la matrice A. Vérifiez
que c’est bien le cas pour la matrice A donnée à la question 1. Modifiez la
routine rayonspectral1 en une routine AppVPdom ( Approximation du
Vecteur Propre associé à la valeur propre dominante ) :
[VP,Niter] = AppVPdom(A,X,epsilon,k);
qui fournisse, avec les mêmes données que dans rayonspectral1 en input,
– le vecteur VPNiter+1 de la suite initiée à X1 = X/norm(X) et régie
ensuite par l’équation récurrente
A · VPk
VPk+1 = , k ≥ 1.
norm A · VPk
– le nombre d’itérations Niter effectué dans la boucle sachant que
cette boucle s’arrête dès que, pour la première fois :
min (norm(VP_(Niter) -VP_(Niter-1)),
norm (VP_(Niter)+ VP_(Niter-1))) <= epsilon
Testez cet algorithme avec la matrice A en prenant epsilon=eps,
k=200 et respectivement X = X1,X2,X3,X4 comme à la question 2.
Quelle valeur de Niter obtenez vous dans chacun de ces quatre cas ?
Recommencez avec cette fois epsilon = 10−6 . Calculez les vecteurs
propres (normalisés) de la matrice A en utilisant :
>> [V,D] = eig(A);
>> V
Comparez le vecteur propre Y=V(:,1) ainsi obtenu (correspondant
à la valeur propre de valeur absolue r(A)) et le vecteur VP obtenu
via
>> [VP,Niter] = AppVPdom(A,X,epsilon,k);
(4) Modifiez la routine AppVPdom construite à la question 3 en une nouvelle
routine
function [VP,Niter,errY] = AppVPdom1(A,X,Y,epsilon,k);
qui, en plus des sorties VP et Niter (comme pour la routine AppVPdom
construite à la question 3), fournit aussi la liste errY des erreurs succes-
sives 2
cette erreur à un facteur multiplicatif près, ce contrôle étant de l’ordre de epsilon/(1 − ρ), où ρ
désigne le rapport entre r(A) et le module de la première valeur propre de module strictement
inférieur à r(A). Voir pour cela la preuve de la Proposition 4.2 du cours.
2. Même remarque que précédemment à propos du contrôle d’erreur : cette tolérance epsilon
contrôle l’erreur entre VPNiter et un vrai vecteur propre normalisé (pour la valeur propre ± r(A))
en epsilon/(1 − ρ), où ρ a été défini dans la note 1 précédente.
F. TP7 : LE THÉORÈME DU POINT FIXE EN ACTION SOUS MATLAB 157
3. D’après l’étude faite en cours, cf. la preuve du Théorème 4.1 (du point fixe), l’erreur entre
LNiter et sa limite est alors majorée par epsilon/(1-kappa).
158 F. TP7 : LE THÉORÈME DU POINT FIXE EN ACTION SOUS MATLAB
norm(XX_(Niter)-XX_(Niter-1)) <=epsilon
4. D’après l’étude faite en cours, cf. la preuve du Théorème 4.1 (du point fixe), l’erreur entre
XXNiter et sa limite (à savoir la solution du système de Cramer que l’on tente d’approcher) est
alors majorée par epsilon/(1-r(A)), où A = D−1 ∗ E ou A = T−1 ∗ F suivant le cas (Jacobi ou
Gauß-Seidel).
F. TP7 : LE THÉORÈME DU POINT FIXE EN ACTION SOUS MATLAB 159
7 5 9 10
utilisée dans la section 1.2.2 du cours. Comparez avec le résultat donné
par cond(A). En utilisant les routines
>> help cond
>> cond(A,1)
>> cond(A,inf)
calculez le conditionnement de la matrice A relativement au choix de la
norme matricielle k k1 et de la norme matricielle k k∞ .
(2) Perturbez la matrice A en lui ajoutant la perturbation :
0 0 .1 .2
.08 .04 0 0
DeltaA := 0
−.02 −.11 0
−.01 −.01 0 −.02.
Calculez les solutions XX et XXX des deux systèmes de Cramer
A ∗ XX = B, (A + DeltaA) ∗ XXX = B
si B=[32;23;33;31], d’abord par la méthode directe :
>> XX = A^{-1} * B ;
>> XXX= (A+DeltaA)^(-1)*B ;
puis par les méthodes itératives :
>> [XX1,Niter] = GausssSeidel(A,B,zeros(4,1),10^(-8),k);
>> [XXX1,Niter] = GaussSeidel(A,B,zeros(4,1),10^(-8),k);
Calculez le coefficient d’amplification d’erreur relative :
>> (norm(XXX-XX)/norm(XXX))*(norm(A)/norm(DeltaA))
Que remarquez vous ? Quel est l’ordre de grandeur de ce coefficient d’am-
plification des erreurs relatives ?
(3) On perturbe maintenant le vecteur B (de la question 2) par le vecteur
DeltaB = [.01;-.01;.01;-.01]. Résoudre (de deux manières différentes,
comme à la question 2 (directement ou alors de manière itérative en uti-
lisant l’algorithme de Gauß-Seidel) le système de Cramer :
A ∗ XXXX = DeltaB.
Calculez encore le coefficient d’amplification de l’erreur relative :
>> (norm(XXXX-XX)/norm(XXXX))*(norm(B)/norm(DeltaB))
Que remarquez vous à nouveau ? Quel est l’ordre de grandeur de ce coef-
ficient d’amplification des erreurs relatives ?
Projets (à traiter par binôme ou trinôme)
161
162 PROJETS (À TRAITER PAR BINÔME OU TRINÔME)
à [Peyr] pour une présentation (contenant des routines sous Scilab) au traitement
numérique des images (voir aussi le document ressources du ministère [TermS]).
7. Le principe de l’interpolation trigonométrique et le phénomène de
Gibbs. Le point de départ de ce projet sera l’exercice 5 de la feuille de TP 5, que l’on
approndira avec d’autres exemples de fonctions ayant beaucoup de discontinuités
(construites par exemple sur le principe des fonctions fractales). Le logiciel utilisé
sera MATLAB, accompagné si nécessaire du Toolbox SignalProcessing.
8. Le théorème du point fixe et les moteurs de recherche sur le web
(Google). On présentera des modèles de graphes orientés et la recherche de dis-
tribution d’équilibre. On pourra s’aider des articles introductifs de M. Eisermann
[Eiser] ou [Eiser0], ou bien du document ressources établi par le Ministère de
l’Education Nationale pour les enseignants en Spécialité Mathématiques en Termi-
nale S [TermS]. Le logiciel utilisé sera de préférence MATLAB (plus adapté au calcul
scientifique), mais l’on pourra aussi utiliser Maple.
9. Les méthodes d’accélération de convergence à l’épreuve des approxi-
mations du nombre π. On présentera quelques méthodes conduisant à la re-
cherche de décimales du nombre π et faisant intervenir le principe de l’accélération
de convergence (formule de John Machin, exercice 4 du TP4, couplé avec le principe
de l’accélération de convergence de Richardson, etc.). Le logiciel exploité pour ce
projet sera Maple.
Annales : examen 2011-2012, session 1 (1h30)
Exercice I.
1. Soit A la matrice 7 × 7 suivante :
−6 1 2 −1 −1/2 1/4 1/2
1 5 1/3 −1/3 0 −1 1/2
0 −1 4 0 1/2 −1/3 2
A=
−1/2 1 1/3 −3 1/2 1/3 0 .
2 −1 −1/2 1/3 −5 1 0
−1/2 1 1/3 2 −1 7 1
2 −3 −1 1 0 −2 10
1
0
−2
AX =
−1
(∗)
0
3
1
−5 1 2 −1 −1/2 1/4 1/2
1 5 1/3 −1/3 0 −1 1/2
0 −1 4 0 1/2 −1/3 2
B=
−1/2 1 1/3 −3 1/2 1/3 0
?
2 −1 −1/2 1/3 −5 1 0
−1/2 1 1/3 2 −1 7 1
2 −3 −1 1 0 −2 10
165
166 ANNALES : EXAMEN 2011-2012, SESSION 2 (1H30)
Exercice 2.
1. Soit N un entier strictement positif donné. Donner
√ la définition du polynôme
d’interpolation de Lagrange de la fonction x 7→ 1/ x aux points xj = 1 + j/N ,
j = 0, ..., N .
2. Expliquer le principe de la méthode de Aitken conduisant au calcul récursif de
ce polynôme de Lagrange PN .
√
3. Donner une majoration de sup[1,2] |1/ x − PN (x)|.
Exercice 3.
1. Expliciter les méthodes d’Euler explicite, d’Euler implicite et d’Euler modifiée
pour la résolution de l’équation différentielle :
y 0 (t) = sin(y 2 (t)) + t2
sur l’intervalle [0, 1], avec la condition initiale y(0) = 1. On prendra comme pas
h = 1/N , où N est un entier strictement positif.
2. Quel est l’ordre de chacune de ces trois méthodes ?
Exercice 4.
1. Soient P et Q les polynômes en deux variables :
P (X, Y ) = X 3 Y 2 + 5X 2 Y − X(3Y 2 + 1) − Y 2 − 1
Q(X, Y ) = X 2 (Y 3 + 2Y + 4) − XY + 1.
Soit E l’ensemble des couples (x, y) de points de C2 tels que P (x, y) = Q(x, y) =
0. Déterminer des polynômes A et B d’une variable, à coefficients réels, dont on
précisera le degré, tels que :
(x, y) ∈ E =⇒ A(x) = B(y) = 0 .
2. Peut on résoudre par radicaux l’équation B(y) = 0 dans C ? Que retournerait
MAPLE si on lui soumettait l’instruction solve(B(y)=0) ? Expliciter une méthode
approchée d’ordre deux permettant de calculer asymptotiquement les racines réelles
de l’équation B(y) = 0.
Annales : examen 2012-2013, session 1 (1h30),
texte + corrigé
Le théorème de Rolle assure d’autre part que, puisque la suite (xn )n≥1
n’est pas stationnaire, alors nécessairement f(xn ) 6= f(xn+1 ) pour n assez
grand : si ce n’était pas le cas, on devrait avoir en effet, pour un tel n,
xn = xn+1 , car, sinon, il y aurait (par Rolle) un zéro de f0 entre xn et
xn+1 , ce qui est impossible du fait que la suite (xn )n≥1 est censée converger
vers ξ tel que f0 (ξ) 6= 0 et que f0 est continue en ce point ξ. Pour n assez
grand, on a donc
1
xn − xn−1 = f(xn−1 ) × 0 + o(1) = o(1).
f (ξ)
167
168 ANNALES : EXAMEN 2012-2013, SESSION 1 (1H30), TEXTE + CORRIGÉ
On en déduit que
f(ξ) = lim f(xn−1 ) = 0.
n→+∞
L’ordre de convergence √ de la méthode de la sécante est au moins égal au
nombre d’or (1 + 5)/2 d’après le cours (Proposition 2.3) et les TP. En
fait, il y a égalité.
(3) Rédiger, en complétant le code Maple suivant, une procédure récursive qui
calcule xN pour N ≥ 1, a,b étant des rationnels fixés (la fonction f ayant
été préalablement déclarée sous Maple) :
algorithme := proc(f,a,b,N) option remember;
local ... ;
if N=1 then
... ;
elif N=2 then
... ;
else
if ... then
... ;
else
... ;
... ;
... ;
end if;
end if;
end proc;
Expliquer pourquoi les deux boucles if ... end if sont essentielles. Que
se passerait-il en particulier si l’on omettait la première d’entre elles ?
Quels sont à la fois le sens et l’intérêt de l’instruction option remember ?
seconde boucle if ... end if est tout aussi essentielle, mais pour une autre rai-
son : afin d’éviter cette fois une intempestive division par zéro lors de l’implémentation
de la procédure (la procédure récursive générant la suite (xn )n≥1 contient d’ailleurs
dans sa définition même cette seconde boucle if ... end if). L’instruction option
remember (dans une procédure récursive) permet de conserver en mémoire les
résultats et donc de ne pas refaire un calcul déjà fait lors des appels successifs
à la récursivité ; l’intérêt de cette instruction réside dans l’allègement qui en résulte
concernant le temps d’exécution time du code (évalué en temps CPU).
Exercice 2.
(1) Rappeler comment déclarer sous l’environnement MATLAB la fonction
t2 × (cos(t)) × (ln(t2 + 1))
f : t ∈ R 7−→ ,
2 + t4 − (sin(t))2
de manière à cette fonction puisse aussi être directement évaluée aussi
bien sur un vecteur ligne t que sur un simple scalaire t.
On utilise pour cela la commande :
>> f =
inline (’(t.^2).*cos(t).*(log(t.^2+1))./(2+t.^4-(sin(t)).^2)’,’t’);
(2) Rédiger un code (sous MATLAB) :
function Im = trapezes (a,b,f,m) ;
qui, étant donnés deux réels a < b et un entier m ∈ N∗ , renvoie en
sortie, une fois exécuté, l’intégrale approchée Im de f sur le segment
[a, b], calculée suivant la méthode des trapèzes composite, avec comme
pas hm = (b − a)/M, où M = 2m . Que vaut (en vous reportant au cours) le
plus grand entier k ∈ N∗ tel que l’on puisse affirmer :
Z b b − a k
f(t) dt − Im = O
a 2m
lorque m croit vers l’infini ?
Il suffit d’utiliser sur chaque segment de la subdivision de [a, b] la formule
(5.18), puis d’ajouter : les extrémités a et b sont prises en compte avec un
facteur hm /2, tandis que les 2m autres nœuds de la subdivision sont pris
en compte avec un facteur hm . Cela donne donc le code MATLAB suivant :
function Im = trapezes (a,b,f,m) ;
h = (b-a)/(2^m) ;
t = a:h:b ;
Im = h *((f(t(1)) + f(t(2^m+1)))/2 + sum(f(t(2:2^m))));
L’ordre de la méthode des trapèzes vaut p = 3 (voir la formule (5.23)
dans la sous-section 5.2.2 du cours). Lorsque cette méthode des trapèzes
est utilisée de manière composite avec un pas h = (b − a)/M, l’erreur
Rb
absolue commise entre l’intégraleexacte a f(t) dt et sa version approchée
est majorée en M×O ((b−a)/M)3 = O ((b−a)/M)2 . Le plus grand entier
k possible ici pour que l’assertion exigée soit valide est k = 2.
170 ANNALES : EXAMEN 2012-2013, SESSION 1 (1H30), TEXTE + CORRIGÉ
end
Prendre epsilon=eps revient à décider que epsilon correspond à l’er-
reur machine (2−52 si l’on travaille en double précision). Ceci signifie que
l’exécution du code s’arrête (si toutefois on y parvient en moins de k
itérations) dès que la machine devient incapable de distinguer du point de
vue numérique les états au cran Niter et au cran Niter + 1 de la variable
XX.
Bibliographie
[Eiser0] M. Eisermann, l’algorithme Pagerank de Google, une promenade sur la toile, disponible
en ligne sur :
[Link]
[Eiser] M. Eisermann, Comment fonctionne Google ?
[Link]
[MathL2] J.M. Marco, P. Thieullen, J.P. Marco (ed.), Mathématiques L2, Pearson Education,
2007.
[MathAp] A. Yger & J.A. Weil (ed.), Mathématiques Appliquées L3, Pearson Education, Paris,
2009.
[Peyr] G. Peyré, le traitement numérique des images, disponible sur :
[Link]
[TermS] Document Ressources pour la classe de Terminale Générale et Technologique,
Mathématiques série S, enseignement de spécialité :
[Link]
/20/8/LyceeGT_ressources_SpeMath_Matrices_218208.pdf
[Vial] G. Vial, Le système proie-prédateur de Volterra-Lotka, Mars 2011, disponible en ligne sur :
[Link]
[Y0] A. Yger, Mathématiques de base (MIS 101, cours 2007-2008)
[Link]
[Yan] A. Yger, Analyse 1 (MIS201, cours 2012-2013) :
[Link]
[Y1] A. Yger, Calcul Symbolique et Scientifique, polycopié de l’UE MHT 304 (2010-2011) :
[Link]
[Y2] A. Yger, Algorithmique Numérique, polycopié de l’UE MHT 632 (2010-2011) :
[Link]
[Zim] P. Zimmerman Peut-on calculer sur ordinateur ?, Leçon de Mathématiques d’aujourd’hui,
Bordeaux, Octobre 2010 :
[Link]
173
Index
175
176 INDEX
sécante, 29
méthode de la, 31
ordre de la méthode de la, 35
synopsis de l’algorithme de la, 31
Schmidt
Ehrard, 57
Seidel
Philipp Ludwig von, 75
Simpson
méthode de, 92
Thomas, 92
singulière, valeur, 69
solve
commande MAPLE, 20
sous-normal
nombre, 9