Mat 249Bernard.Parisse@ujf-grenoble.fr2009/10 |
Table des matières
1 Présentation du module
Dans ce module, on introduira moins de notions nouvelles que
dans d’autres modules de mathématiques, par contre on insistera
sur le calcul effectif, si possible efficace, et sur le controle
de la précision des résultats,
ceci explique la part importante consacrée aux TP
(18h de cours en 12 séances, 18h de TD en 12
séances et 24h de TP). On présentera par
exemple des méthodes de calcul des fonctions usuelles (racine
carrée, trigonométriques, ...), il s’agira
non seulement de savoir calculer une valeur numérique, mais aussi
de pouvoir majorer l’écart entre la valeur trouvée et
la valeur exacte, en utilisant des théorèmes du cours.
Les calculs se feront dans la mesure du possible sur ordinateur
ou sur calculatrices.
Les thèmes abordés seront :
-
calcul exact et approché, représentation des données
- Suites récurrentes, méthode du point fixe, de Newton
- Interpolation polynômiale (évaluation, interpolation de Lagrange)
- Séries de Taylor et approximation des fonctions usuelles
- Arithmétique des polynômes (PGCD, Bézout, factorisation,
décomposition en éléments simples)
- Intégration
- Algèbre linéaire
L’évaluation se fait sur :
-
1/4 : un DS à mi-semestre
- 1/4 : certains compte-rendus de TP (à rédiger seul ou en binome),
- 1/2 : l’examen final
Les calculatrices et les netbooks de taille d’écran plus petits
que 10 pouces sont autorisées au DS et à l’examen final (prêt possible
de netbooks pour le semestre).
2 Représentation des nombres et autres données, calcul exact/approché
Résumé:
Types de base : entier machine, entier long, flottant machine et
multiprécision (Base 2, base 10, BCD).
Types composés : complexes, polynomes (représentation dense/creuse),
symboles, listes (vecteurs, matrices), expressions, fonctions.
Erreur relative, erreur absolue, erreur d’arrondi, +/-, */%
Algorithme, complexité, exemple puissance modulaire, algorithme de Horner.
Les principaux ensembles de nombres en mathématiques sont les entiers
positifs ℕ et relatifs ℤ, les rationnels ℚ, les réels ℝ
et les complexes ℂ. Sur ordinateur, on peut représenter ces
nombres de manière exacte dans certains cas, approchée dans d’autres.
2.1 Entiers courts et longs
Proposition 1
Division euclidienne de deux entiers : si a et b sont
deux entiers, a ≥ 0, b>0, il existe un unique couple (q,r) tel que
Preuve : On prend pour q le plus grand entier tel que a−bq>0.
La division euclidienne permet d’écrire un nombre entier, en utilisant
une base b et des caractères pour représenter les entiers
entre 0 et b−1. Nous écrivons les nombres entiers en base b=10
avec comme caractères les chiffres de 0 à 9.
Les ordinateurs utilisent des circuits binaires pour stocker
les informations, il est donc naturel d’y travailler en base 2
en utilisant comme caractères 0 et 1 ou en base 16 en utilisant
comme caractères les chiffres de 0 à 9 et les lettres de A à F.
En général, pour trouver l’écriture d’un nombre en base b (par
exemple b=2),
on effectue des divisions euclidienne successives par b du nombre puis
de ses quotients successifs jusqu’à ce que le quotient soit 0 et on
accolle les restes obtenus (premier reste à droite, dernier
reste à gauche).
Inversement, pour retrouver un entier d à partir
de son écriture dn...d0, on traduit les divisions euclidiennes
successives en
d | = | ( ... ((dn b +dn−1)b + dn−2)...+d1)b+d0 |
| = | dn bn + dn−1 bn−1 + ... + d0 |
|
Par exemple, vingt-cinq s’écrit en base 16 19
car 25 divisé
par 16 donne quotient 1, reste 9.
En base 2, on trouverait 00011001
car 25=24+23+1.
On peut effectuer les opérations arithmétiques de base
(+,-,*, division) directement en base 2 (ou 16). Par exemple
la table de l’addition est 0+0=0, 0+1=1+0=1 et 1+1=0 je retiens 1,
donc :
01001111
+ 01101011
----------
10111010
Exercice :
comment passe-t-on simplement de la représentation d’un nombre en
base 2 à un nombre en base 16 et réciproquement ?
Les microprocesseurs peuvent effectuer directement
les opérations arithmétiques de base sur les entiers “machine”
(déclinés en plusieurs variantes selon la taille et la
possibilité d’avoir un signe). Noter que la division de deux
entiers a et b n’a pas la même signification que la division de deux
réels, comme elle ne tomberait pas forcément juste,
on calcule le quotient et le reste de la division euclidienne.
Ces entiers machines permettent
de représenter de manière exacte des petits entiers relatifs
par exemple un entier machine signé sur 4 octets est
compris entre [−231,231−1]. Selon le microprocesseur
les 4 octets représentant l’entier sont stockés
par adresse mémoire décroissante ou croissante (big ou little endian).
Sur certains systèmes
(dits BCD), on écrit les entiers en base 10, chaque chiffre
occupant 4 bits (qui normalement sert à stocker un chiffre en base 16).
Les microprocesseurs correspondants ont un flag leur permettant d’effectuer
les opérations sur des nombres vu en représentation BCD (base 10) ou
hexadécimale (base 16).
Ces entiers machines permettent de faire très rapidement
du calcul exact sur les
entiers, mais à condition qu’il n’y ait pas de dépassement
de capacité, par exemple pour des entiers 32 bits, 231+231
renverra 0. Ils sont très utilisés en calcul formel pour
les algorithmes dits modulaires (on travaille modulo un entier assez
petit). Pour travailler avec des entiers plus grands, on doit
utiliser
des entiers de taille plus grande, mais il faut alors programmer les
opérations de base et décider d’un mécanisme de stockage, par
exemple en représentant un entier par une zone mémoire
commencant par la taille et suivie par l’écriture à l’aide
d’entiers machines de l’entier (en base 232). Bien entendu,
plus les entiers sont grands, plus les opérations seront longues,
par exemple l’addition de deux entiers longs de taille N nécessite
un temps proportionnel à N, leur multiplication par l’algorithme
élémentaire nécéssite un temps proportionnel à N2
(mais il existe des algorithmes plus efficaces, par exemple
Karatsuba ou FFT, cf. Knuth, The Art of Computer Programming
ou la documentation de la librairie GMP).
2.2 Les rationnels.
On sait donc représenter les entiers, pour les
rationnels, il suffit
de les représenter comme un couple d’entiers correspondant à leur
écriture sous forme de fraction irréductible avec un
dénominateur positif.
Proposition 2
L’algorithme d’Euclide
permet de calculer le PGCD (plus grand commun diviseur) de 2 entiers,
écrit ici en syntaxe Xcas :
pgcd(x,y):={
local r;
while (y!=0){
r:=irem(x,y); // reste de x par y
x:=y; // PGCD(x,y)=PGCD(y,r) donc on decale
y:=r;
}
return x; // c'est le resultat car PGCD(x,0)=x
}
Preuve :
on utilise le fait qu’un nombre divise a et b si et seulement si
il divise r=a−bq et b. Le PGCD de a et b est donc le PGCD
de b et du reste de la division euclidienne de a par b. Comme
le reste est en valeur absolue plus petite que |b|, la taille
des variables x,y,r
décroit à chaque itération. Arrive un
moment où le reste est nul, le PGCD est alors l’entier par lequel
on a divisé.
Il existe des variantes de cet algorithme un peu plus efficaces
lorsque les nombres sont représentés en base 2 (PGCD binaire,
voir par exemple A. Cohen).
On utilise cet algorithme et la division euclidienne pour
simplifier une fraction d’entiers par le PGCD du numérateur et
du dénominateur pour l’écrire sous forme irréductible.
Les calculs
sont maintenant exacts et sans limitation de capacité (ou presque,
la taille des entiers longs est bornée parce que
la taille du champ mémoire fixant la longueur de stockage est
bornée) mais souvent trop lents pour les calculs numériques
usuels (par exemple pour calculer la valeur approchée de cosinus
23 degrés 27 minutes).
On utilise alors un autre type dont
les calculs de base sont gérés par le microprocesseur (ou son
coprocesseur arithmétique).
2.3 Les réels
On se ramène d’abord au cas des réels positifs, en machine
on garde traditionnellement un bit pour stocker le signe du réel
à représenter.
2.3.1 Virgule fixe et flottante.
La première idée qui vient naturellement serait d’utiliser
un entier et de déplacer la virgule
d’un nombre fixe de position, ce qui revient à mulitplier
par une puissance (négative) de la base. Par exemple en base 10 avec un
décalage de 4, 1234.5678
serait représenté par 12345678
et 1.2345678
par
12345
(on passe de l’entier au réel par multiplication
par 10−4. L’inconvénient d’une telle représentation est
qu’on ne peut pas représenter des réels grands ou petits,
comme par exemple le nombre d’Avogadro, la constante de Planck, etc.
D’où l’idée de ne pas fixer la position de la virgule, on parle
alors de représentation à virgule flottante ou de nombre flottant : on
représente un nombre par deux entier, l’un appelé mantisse
reprend les chiffres significatifs du réel sans virgule, l’autre
l’exposant, donne la position de la virgule. Attention, le séparateur
est un point et non une virgule dans la grande
majorité des logiciels scientifiques.
On sépare
traditionnellement la mantisse de l’exposant par la lettre e
.
Par exemple 1234.5678
peut être représenté
par 12345678e-8
(mantisse 12345678
, exposant -8)
mais aussi par 1234567800e-10
.
Naturellement, sur un ordinateur, il y a des limites pour les entiers
représentant la mantisse m et l’exposant e.
Si on écrit les nombres en base
b, la mantisse m s’écrira avec un nombre n fixé de chiffres (ou
de bits en base 2), donc m ∈ [0,bn[. Soit un réel x représenté
par
Si m∈ [0,bn−1[, alors on peut aussi écrire x=m′ be−1 avec
m′=mb ∈ [0,bn[, quelle écriture faut-il choisir?
Intuitivement, on sent qu’il vaut mieux prendre m′ le plus grand
possible, car cela augmente le nombre de chiffres significatifs (alors
que des 0 au début de m ne sont pas significatifs).
Ceci est confirmé par le calcul de l’erreur d’arrondi pour
représenter un réel. En effet, si x est un réel non nul, il ne
s’écrit pas forcément sous la forme mbe, on doit l’arrondir,
par exemple au plus proche réel de la forme mbe. La distance
de x à ce réel est inférieure ou égale à la moitié
de la distance entre deux flottants consécutifs,
mbe et (m+1)be, donc l’erreur d’arrondi
est inférieure ou égale à be/2. Si on divise par x ≥ mbe,
on obtient une erreur relative d’arrondi majorée par 1/(2m).
On a donc intérêt à prendre m le plus grand possible pour
minimiser cette erreur. Quitte à mulitplier par b, on peut
toujours se ramener (sauf exceptions, cf. ci-dessous),
à m ∈ [bn−1,bn[, on a alors
une erreur d’arrondi relative majorée par
On appelle flottant normalisé un flottant tel que m ∈
[bn−1,bn[. Pour écrire un réel sous forme de flottant
normalisé, on écrit le réel en base b, et on déplace
la virgule pour avoir exactement n chiffres non nuls avant la
virgule et on arrondit (par exemple au plus proche).
L’exposant est égal au décalage effectué.
Notez qu’en base 2, un flottant normalisé commence forcément
par 1, ce qui permet d’économiser un bit dans le stockage.
Ainsi, l’erreur d’arrondi commise lorsqu’on
représente un réel (connu exactement) par un double normalisé
est une erreur relative inférieure à
de 2−53 (b=2 et n=52+1 pour les doubles).
Exemples :
Il existe une exception à la possibilité de normaliser les flottants,
lorsqu’on atteint la limite inférieure de l’exposant e.
Soit en effet em le plus petit exposant des flottants normalisés
et considérons les flottants x=bem(1+1/b) et y=bem. Ces
flottants sont distincts,
mais leur différence n’est plus représentable par un flottant normalisé.
Comme on ne souhaite pas représenter x−y par 0,
(le test x==y renvoie faux), on introduit les flottants
dénormalisés, il s’agit de
flottants dont l’exposant est l’exposant minimal représentable sur
machine et dont la mantisse appartient à [0,bn−1[. Par exemple
0 est représenté par un flottant dénormalisé de mantisse 0
(en fait 0 a deux reprśentation, une de signe positif et une de
signe négatif).
Enfin, on utilise traditionnellement une valeur de l’exposant pour
représenter les nombres plus grands que le plus grand réel
reprśentable sur machine (traditionnellement appelé plus ou
moins infini)
et les erreurs (par exemple 0./0. ou racine carrée d’un nombre
réel négatif, traditionnellement appelé NaN, Not a Number).
Exercice : quels sont les nombres réels représentables exactement
en base 10 mais pas en base 2 ?
Si on écrit 1/10 en base 2 avec 53 bits de précision, puis que
l’on arrondit avec 64 bits de précision, ou si on écrit 1/10 en
base 2 avec 64 bits de précision, obtient-on la même chose ?
Les ordinateurs reprśentent générallement les flottants en base 2
(cf. la section suivante pour
plus de précisions), mais cette représentation n’est pas utilisée
habituellement par les humains, qui préfèrent compter
en base 10. Les ordinateurs effectuent donc la conversion dans
les routines d’entrée-sortie. Le format standard utilisé
pour saisir ou afficher un nombre flottant dans un logiciel
scientifique est composé d’un nombre à virgule
flottante utilisant le point comme séparateur décimal (et
non la virgule) suivi si nécessaire de la lettre e
puis de l’exposant,
par exemple 1.23e-5
ou 0.0000123
. Dans les
logiciels de calcul formel, pour distinguer un entiers
représentés par un entier
d’un entier représenté par un flottant on écrit
l’entier suivi de .0
par exemple 23.0
.
Remarque :
Les microprocesseurs ayant un mode BCD peuvent avoir un format
de représentation des flottants en base 10, les nombres décimaux
comme par exemple 0.3 peuvent être représentés exactement.
Certains logiciels, notamment maple, utilisent par défaut des
flottants logiciels en base 10 sur des microprocesseurs sans mode BCD,
ce qui entraine une baisse de
rapidité importante pour les calculs numériques (on peut
partiellement améliorer les performances en utilisant evalhf
en maple).
2.3.2 Les flottants au format double
Cette section développe les notions de la section précédente
pour les flottants machine, utilisables dans les langage de
programmation usuels, elle peut être omise en première lecture.
La représentation d’un double
en mémoire se compose de 3 parties : le bit
de signe s=± 1 sur 1 bit,
la mantisse M ∈ [0,252[ sur 52 bits,
et l’exposant e ∈ [0, 211[ sur 11 bits. Pour les nombres
“normaux”, l’exposant est en fait compris entre 1 et 211−2,
le nombre représenté est le rationnel
Pour écrire un nombre sous cette forme, il faut d’abord chercher par
quel multiple de 2 il faut le diviser pour obtenir un réel r dans
[1,2[, ce qui permet de déterminer l’exposant e. Ensuite on
écrit la représentation en base 2 de r−1 ∈ [0,1[.
Exemples :
-
-2
Signe négatif. Il faut diviser sa valeur absolue
2 par 21 pour être entre 1 et 2 dont
e+1−210=1, l’exposant est e=210. On a alors r=1, r−1=0.
Représentation
1 10000000000 00000000...0000
- 1.5=3/2
Signe positif, compris entre 1 et 2 dont l’exposant vérifie
e+1−210=0 soit
e=210−1=29+28+27+26+25+24+23+22+21+20.
On a r−1=1/2=2−1. D’où la représentation
0 01111111111 10000000...0000
- 6.4=32/5
Positif. Il faut le diviser par 22 pour avoir 8/5 ∈ [1,2[
donc e+1−210=2 soit e=210+1. Ensuite r=3/5 qu’il faut
écrire en base 2 (cf. section précédente),
on écrit donc les 52 premiers éléments du développement
avec une règle d’arrondi du dernier bit au nombre le plus proche.
Ici le bit suivant le dernier 1001
est un 1
, on arrondit
donc à 1010
. D’où la représentation
0 1000000001 100110011001...10011010
On observe que la représentation en base 2 de 6.4 a du être
arrondie (car elle est infinie en base 2) bien qu’elle soit exacte
(finie) en base 10.
Seuls les entiers et les rationnels dont le dénominateur est une puissance
de 2 peuvent être représentés exactement.
Ceci entraine des résultats qui peuvent surprendre
comme par exemple le fait que
0.3 - 3*0.1
n’est pas nul.
Des représentations spéciales (avec e=0 ou e=211−1)
ont été introduites
pour représenter ± ∞ (pour les flottants plus grands
en valeur absolue que le plus grand flottant représentable), et pour
représenter les
nombres non nuls plus petits que le plus petit flottant représentable
de la manière exposée ci-dessus (on parle de flottants dénormalisés),
ainsi que le nombre NaN (Not a Number) lorsqu’une opération a un résultat
indéfini (par exemple 0/0).
Remarque : Sur les processeurs compatibles avec les i386,
le coprocesseur arithmétique i387 gère en interne des flottants
avec 80 bits dont 64 bits de mantisse. Sur les architectures 64 bits
(x86 ou AMD), le jeu d’instruction SSE permet de travailler avec
des flottants de 128 bits. Le compilateur gcc permet d’utiliser
ces flottants longs avec le type long double
ou
les types __float80
et __float128
en utilisant
un drapeau de compilation du type -msse
2.3.3 Opérations sur les flottants
Les opérations arithmétiques de base sur les flottants
se font de la manière suivante :
-
addition et soustraction : on détecte s’il faut additionner
ou soustraire en valeur absolue en analysant les signes,
on détermine l’exposant le plus grand et on décale la partie mantisse
du flottant dont l’exposant est le plus petit pour se ramener à additionner
deux entiers (partie mantisses correspondant au même exposant),
on décale à nouveau la partie mantisse en modifiant l’exposant
après l’opération pour normaliser le flottant
- multiplication : on additionne les exposants et on multiplie
les parties mantisses (vus comme des entiers), on tronque et
on ajuste l’exposant si nécessaire
- division : on soustrait les exposants et on divise les parties
mantisses (division “à virgule”), on tronque et on
ajuste l’exposant si nécessaire
La représentation des nombres réels par des doubles présente
des avantages, les opérations arithmétiques
sont faites au plus vite par le microprocesseur
(quoique sur les microprocesseurs plus anciens, par exemple Saturn des
calculatrices HP, Z80 des calculatrices TI8x, la multiplication
et la division n’est pas une opération de base du microprocesseur,
elle doit être codée à partir des opérations d’addition, soustraction,
décalage de bits, et/ou logique. A contrario,
les coprocesseurs arithmétiques (intégrés sur les microprocesseurs
de PC) proposent même
le calcul des fonctions usuelles (trigonométriques, racine carrée, log et exp)
sur le type double) et utilisent des formats de représentation interne
ayant plus de 64 bits pour les doubles, ce qui permet de limiter
dans certains cas les erreurs d’arrondi).
Par contre, des erreurs vont être introduites,
on parle de calcul approché par opposition au calcul exact sur les
rationnels. En effet, la représentation doit d’abord arrondir
tout réel qui n’est pas un rationnel dont le dénominateur
est une puissance de 2. Ensuite chaque opération va entrainer
une propagation de ces erreurs et va y ajouter une erreur d’arrondi
sur le résultat.
Enfin, l’utilisation du type double peut provoquer un dépassement
de capacité (par exemple 100!*100!
).
Pour diminuer ces erreurs et les risques de dépassement de
capacité, il existe des types flottants multiple précision,
qui permettent de travailler avec un nombre fixé à l’avance
de décimales. Les calculs sont plus longs mais les erreurs
plus faibles. Attention, il s’agit toujours de calcul approché!
De plus, pour des quantités dont la valeur est déterminée
de manière expérimentale, la source principale de propagation
d’erreurs est la précision des quantités initiales, il ne sert
souvent à rien d’utiliser des types flottants multiprécision car les
erreurs dus à la représentation (double) sont négligeables devant
les erreurs de mesure.
2.3.5 Erreur absolue, relative et propagation des erreurs.
On a vu précédemment que pour représenter un réel, on devait
l’arrondir, ce qui introduit une erreur même si le réel est
connu exactement (par exemple 1/10).
Voyons comment se propagent les erreurs
dans les opérations arithmétiques
de base : on distingue l’addition, la multiplication
et l’inversion. La soustraction se ramène à l’addition car
le calcul de l’opposé n’introduit aucune erreur nouvelle.
Pour l’addition, si |x −x0| ≤ ε0 et si |y−y0| ≤ ε1
alors par l’inégalité triangulaire (|a+b|≤ |a|+|b|), on a :
|(x+y)−(x0+y0)| ≤ |x−x0| + | y−y0 | ≤
ε0 + ε1 |
on dit que les erreurs absolues s’additionnent.
Définition 1
L’erreur absolue
est définie comme un majorant de la valeur absolue
de la différence entre le nombre réel et
son représentant double :
Mais comme il faut représenter x0+y0
en machine, on doit ajouter une erreur d’arrondi, qui est
proportionnelle à la valeur absolue de x0+y0 d’où la notion
d’erreur relative :
Définition 2
L’erreur relative est égale à l’erreur absolue divisée par
la valeur absolue du nombre
Remarquons au passage que les erreurs de mesure expérimentales sont
pratiquement toujours des erreurs relatives.
Donc lorsqu’on effectue une addition (ou une soustraction) de deux réels
sur machine, on doit additionner les deux erreurs absolues sur les opérandes
et ajouter une erreur d’arrondi (relative de 2−53, à titre
d’exercice, on pourra vérifier que cette erreur
d’arrondi est majorée par l’erreur absolue de la somme
x+y dès l’instant où
x et y ont eux-même une erreur d’arrondi).
Lorsqu’on effectue une multiplication de deux nombres x,y dont les
représentants x0,y0 sont non nuls, on a
⎪
⎪
⎪
⎪ | | | ⎪
⎪
⎪
⎪ | =
| ⎪
⎪
⎪
⎪ | | | | −1) | ⎪
⎪
⎪
⎪ |
= | ⎪
⎪
⎪
⎪ |
( | | −1)( | | −1)+( | | −1)+( | | −1)
| ⎪
⎪
⎪
⎪ |
l’erreur relative est donc la somme des erreurs relatives et du produit
des erreurs relatives (on peut souvent négliger le produit devant la somme). Il
faut aussi y ajouter une erreur relative d’arrondi de 2−53 sur x0 y0.
On observe que la multiplication est une opération posant moins
de problèmes que l’addition, car on manipule toujours des erreurs
relatives, par exemple si l’erreur relative sur deux doubles
x et y non nuls
est de 2−53, alors l’erreur relative sur xy sera de
2−53 + 2−53 + 2−106 + 2−53 ≈ 3 × 2−53 |
Lorsque l’erreur relative sur les données est grande devant 2−53,
l’erreur relative d’arrondi final est négligeable, on peut alors dire que
les erreurs relatives s’additionnent pour un produit (c’est aussi vrai
pour un quotient: exercice!).
Par contre, si on additionne deux nombres dont le représentant de
la somme est proche de 0, la somme des erreurs absolues peut
devenir non négligeable par rapport à la somme des représentants, entrainant
une erreur relative très grande. Par exemple si x est représenté
par x0=1+2−52 avec
une erreur d’arrondi de 2−53 et
y par y0=−1 avec la même erreur d’arrondi,
l’addition de x et y renvoie 2−52 avec une erreur
absolue de 2 * 2−53 (ici il n’y a pas d’arrondi lorsqu’on fait la somme).
C’est une erreur relative de 1 (qui domine largement
l’erreur d’arrondi) ce qui signifie que dans la mantisse, seul le
premier bit sur les 52 a un sens, la perte de précision est très grande.
Une autre conséquence importante est que l’addition de réels sur machine
n’est pas une opération associative,
par exemple
(2.0−53+2.0−53)+1.0 → 1+2−52 |
alors que
Si on a plusieurs termes
à additionner, il faut commencer par additionner entre eux
les termes les plus petits, pour que les petits termes ne soient
pas absorbés un à un dans les erreurs d’arrondi (les petits ruisseaux
font les grands fleuves).
Exercice : pour calculer la valeur
numérique d’une dérivée de fonction, il vaut mieux
calculer (f(x+h)−f(x−h))/(2h) que (f(x+h)−f(x))/h. Attention
à ne pas prendre h trop petit, sinon x+h=x.
Remarquons néanmoins que les erreurs calculées ici sont des majorations
des erreurs réelles (ou si on préfère l’erreur obtenue dans le pire
des cas), statistiquement les erreurs sur les résultats sont moindres.
Il est d’ailleurs souvent trop difficile de calculer la majoration
rigoureuse de l’erreur pour des calculs complexes.
Lorsqu’on doute de la précision d’un calcul, un test peu couteux consiste
à refaire ce calcul en utilisant des flottants en précision plus
grande et tester si le résultat varie en fonction du nombre de chiffres
significatifs utilisés. On peut aussi faire varier légèrement
les données et observer la sensibilité du résultat.
Si on veut travailler en toute rigueur sans
pour autant calculer les erreurs à priori, il faut utiliser un logiciel
utilisant des intervalles pour représenter les réels (par exemple
la bibliothèque C MPFI).
2.4 Types composés.
Après les nombres réels, on passe aux
nombres complexes :
on utilise un couple (partie réelle, imaginaire) de fractions (exacts)
ou de flottants et les règles habituelles sur les complexes.
Après les nombres, l’objet le plus utilisé dans les systèmes
de calcul formel est probablement le polynôme,
toute simplification
d’une expression se ramène à un moment donné à mettre sous
forme irréductible une fraction de polynômes. Les principales
représentations possibles sont :
-
les polynômes à 1 variable, représentation dense, on stocke
la liste des coefficients du polynôme par ordre croissant ou décroissant
- les polynômes à 1 variable, représentation creuse, on stocke
des paires coefficients, degré pour les coefficients non nuls
- les polynômes à plusieurs variables, représenté récursivement
de manière dense ou creuse (i.e. P(x1,...,xn) vu comme polynôme
en xn à coefficients polynômes dépendant des variables x1,...,xn−1),
ce sont des cas particuliers des 2 cas précédents
- les polynômes à plusieurs variables distribués, on stocke
des monômes, qui sont des paires coefficient, liste d’entiers,
la liste représentant les exposant des variables dans le monôme.
- la représentation symbolique (par exemple x y2−5x+y3)
beaucoup plus difficile à manipuler directement
Algorithmes de base sur les polynômes :
l’évaluation en un point (Horner, cf. TD/TP), la multiplication et
division euclidienne et le PGCD (même algorithme que pour les
entiers mais avec la division euclidienne des polynômes, il existe
des algorithmes plus efficaces, cf. le chapitre sur les polynômes)
Lien avec la représentation en base z (TD).
Les polynômes peuvent servir à représenter des nombres non
rationnels de manière exacte, par exemple les nombres algébriques
(qui sont solutions d’une équation polynomiale).
Les symboles ou noms de variable désignent par exemple
le nom d’une inconnue dans un polynôme, ils sont représentés
par une chaine de caractére et peuvent être affectés à une valeur
pendant une session (la valeur dépend d’un contexte d’exécution
et le remplacement du symbole par sa valeur affectée est appelé
évaluation).
Les expressions, par exemple sin(x)+2*x^2
,
elles peuvent être représentées par des arbres. L’évaluation
d’une expression consiste à remplacer les symboles de l’expression
par leur valeur, puis à effectuer les opérations en tenant compte
de la substitution. Il est parfois souhaitable de ne pas
effectuer certaines opérations de substitution, on empêche
l’évaluation, explicitement (''
) ou implicitement
(par exemple l’affectation n’évalue pas le symbole qu’on va affecter).
Les fonctions ne doivent pas être confondues avec les
expressions, elles associent à leurs arguments une expression. Par
exemple sin
est une fonction, alors que sin(x)
est une expression.
Les conteneurs contiennent plusieurs objets et permettent d’associer
à un indice un objet. Il en existe de plusieurs types, par exemple
les listes et les séquences
dont l’indice est un
entier compris entre 1 (ou 0) et la taille (-1), les tables dont
l’indice est plus général, et les tableaux
(utilisés pour les vecteurs, matrices)
qui sont essentiellement des listes ou des listes de listes de même taille.
Les séquences sont des listes d’objets ordonnés “non récursifs”
(ils ne peuvent contenir des séquences), alors que les listes
peuvent contenir des listes, sinon il n’y a pas de différences.
Dans les logiciels de calcul formel, la plupart du temps
les séquences se notent en indiquant les éléments séparés par
des virgules. Les listes s’en distinguent par les délimiteurs []
.
Il faut prendre garde au fait qu’en général affecter par exemple
l[1]:=3;
à une variable libre l
crée une table et
non une liste.
Remarque: certains logiciels accédent à certains types de
conteneurs uniquement par référence (par exemple maple pour les
vecteurs et matrices), dans ce dernier cas une seule copie des objets
du conteneur existe si on copie de la manière habituelle
une variable contenant un vecteur ou
une matrice dans une autre variable,
la modification d’un élément du conteneur modifie
alors toutes les copies pointant sur ce conteneur. Cette méthode
est plus efficace mais peut être surprenante.
3 Suites itératives et applications
Résumé:
Théorème du point fixe, méthode de Newton,convexité.
Exemple: calcul de valeur approchée de racines carrées,
Résolution d’équations.
3.1 Le point fixe
Soit f une fonction continue sur un intervalle I=[a,b] de ℝ, et
à valeurs dans I (attention à bien choisir I pour que l’image
de I par f reste dans I).
On s’intéresse à la suite
Supposons que un converge vers une limite l ∈ I lorsque
n → +∞, alors la limite doit vérifier
puisque f est continue. On dit que l est un point fixe de f.
Ceci amène à l’idée d’utiliser ces suites pour résoudre numériquement
l’équation f(x)=x. Nous allons donner un théorème permettant
d’assurer que la suite (1) converge, et que la limite
est l’unique solution de f(l)=l sur I.
Définition 3
On dit que f est contractante de rapport k<1 sur I si
∀ x,y ∈ I, |f(y)−f(x)| ≤ k |y−x| |
En pratique, les fonctions f que l’on considèrera seront continument
dérivables, donc d’après le théorème des accroissements finis
f(y)−f(x)=f′(θ) (y−x), θ ∈ [x,y] |
ainsi pour vérifier que f est contractante, on étudie la valeur absolue
de f′ sur I, il suffit de montrer que cette valeur absolue
est strictement inférieure à un réel k<1 pour conclure (il faut
donc chercher le maximum de |f′| sur I. Attention, il s’agit du
maximum de |f′| et pas du maximum de f′, ce qui revient à chercher
le maximum de f′ et de −f′).
On a alors le
Théorème 1 (du point fixe)
si f est contractante de I=[a,b]
dans I de rapport k
alors la suite (1) converge vers l’unique
solution de f(l)=l dans I. On a de plus les encadrements :
|un−l| ≤ kn |b−a|,
|un −l | ≤ | |
(2) |
Démonstration : Tout d’abord si f est contractante, on montre à partir
de la définition de la continuité que f est continue.
Soit g(x)=f(x)−x, alors g est continue, positive en a et négative
en b, il existe donc l∈[a,b] tel que g(l)=0 (théorème des
valeurs intermédiaires). Soit un une suite définie par
(1). On a alors pour tout n
|un+1−l|=|f(un)−f(l)| ≤ k |un−l| |
Donc par une récurrence évidente :
ce qui entraine d’ailleurs que |un−l| ≤ kn |a−b|.
Comme k ∈ [0,1[ , la suite géométrique kn converge vers 0
lorsque n tend vers l’infini, donc un tend vers l.
Notons que l est unique car si l′ est une autre solution
alors |l−l′|=|f(l)−f(l′)| ≤ k|l−l′| donc (1−k)|l−l′| ≤ 0,
or 1−k>0 et |l−l′| ≥ 0 donc |l−l′| doit être nul.
Passons à la preuve de la majoration (2) qui est importante
en pratique car elle donne un test d’arrêt de calcul des
termes de la suite récurrente, on écrit pour m>0 :
un−l= un − un+1 + un+1 − un+2 + ... + un+m−1− un+m
+ um−l |
puis on majore avec l’inégalité triangulaire
|un−l| ≤ | | |un+j−un+j+1| + |um−l| |
puis on applique le fait que f est contractante de rapport k
|un−l| ≤ | | kj |un−un+1| + |um−l| |
soit
|un−l| ≤ | | |un−un+1| + |um−l| |
On fait alors tendre m vers l’infini d’où le résultat.
Exemples :
Cherchons une valeur approchée de √2 par cette méthode.
Il faut d’abord trouver une fonction f dont √2 est un point
fixe, par exemple
On vérifie que f(√2)=√2), puis que f′=−1/(x+1)2
donc f décroit. On va voir si les hypothèses du théorème du point
fixe s’appliquent sur par exemple [1,2]. Comme f est décroissante
f([1,2])=[f(2),f(1)]=[4/3,3/2] qui est bien inclus dans [1,2] .
De plus f′ est comprise entre −1/(1+1)2=−1/4 et −1/(2+1)2=−1/9 donc
|f′|<1/4, f est contractante de rapport 1/4. On peut donc
itérer la suite à partir par exemple de u0=1 et on va converger
vers √2 (en s’en rapprochant à chaque cran d’un rapport
inférieur à 1/4).
Considérons l’équation en x
x− e sin(x) =t, e ∈ [0,1[ |
c’est l’équation du temps utilisée en astronomie pour trouver la
position d’une planète sur son orbite elliptique (e étant l’excentricité
de l’ellipse).
Il n’y a pas de formule exacte permettant de calculer x en fonction de t.
Si on a une valeur numérique pour t, on peut trouver une valeur
numérique approchée de x par la méthode du point fixe, en réécrivant
l’équation sous la forme
On observe que f envoie ℝ dans [t−e,t+e] donc on peut prendre
I=[t−e,t+e], de plus |f′|≤ e <1, f est contractante
de rapport e ∈ [0,1[, le théorème s’applique, il suffit de
prendre une valeur initiale dans [t−e,t+e] et d’itérer la suite
jusqu’à obtenir la précision désirée. Par exemple si on veut une
valeur approchée de x à 10−6 près, il suffira que la différence
entre deux termes successifs de la suite un vérifie
on aura alors bien :
Cette méthode n’est pas toujours optimale, car la vitesse de convergence
vers la limite l est dite “linéaire”, c’est-à-dire
que le temps de calcul pour
avoir n décimales est proportionnel à n (ou
encore il faut effectuer un nombre d’itérations
proportionnel à n, chaque itération faisant
gagner en précision de l’ordre du rapport k de contractance).
En effet, supposons que f′ est continue en l et que 0<L=|f′(l)|<1 .
Il existe alors un intervalle I=[l−η,l+η] tel que
Le théorème des accroissements finis donne alors
|un+1 − l | = |f(un)−f(l)| = |f′(θ)| |un−l|,
θ ∈ [un,l] |
Si u0 ∈ I, alors θ ∈ I donc |u1−l| ≤ |u0−l| et
u1 ∈ I, par récurrence on a pour tout n, un ∈ I
| |un−l| ≤ |un+1 − l| ≤ | | |un−l| |
on a donc par récurrence
⎛
⎜
⎜
⎝ | | ⎞
⎟
⎟
⎠ | | |u0−l| ≤ |un−l| ≤
| ⎛
⎜
⎜
⎝ | | | ⎞
⎟
⎟
⎠ | | |u0−l|
|
Donc pour avoir |un−l| ≤ є il suffit que
⎛
⎜
⎜
⎝ | | | ⎞
⎟
⎟
⎠ | | |u0−l| ≤ є ⇒
n ≥ | | |
et il faut que
⎛
⎜
⎜
⎝ | | ⎞
⎟
⎟
⎠ | | |u0−l| ≤ є
⇒
n ≥ | |
Si f est suffisamment régulière,
il existe une méthode plus rapide lorsqu’on est proche de la racine ou lorsque
la fonction a des propriétés de convexité, c’est la méthode de Newton.
Et même si Newton n’est pas applicable, une simple dichotomie
peut être plus efficace si la constante de contractance est
supérieure à 1/2 (y compris prés de la solution de f(x)=x).
3.2 La méthode de Newton.
La méthode de Newton est une méthode de résolution de l’équation
f(x)=0, attention à la différence avec le théorème du point fixe
qui permet de résoudre numériquement f(x)=x.
Si x0 est proche de la racine r
on peut faire un développement de Taylor à l’ordre 1 de la
fonction f en x0 :
f(x)=f(x0)+(x−x0)f′(x0)+O((x−x0)2) |
Pour trouver une valeur approchée de r, on ne garde que la partie
linéaire du développement, on résout :
f(r)=0 ≈ f(x0) + (r−x0) f′(x0) |
donc (si f′(x0)≠ 0) :
Graphiquement, cela revient à tracer la tangente à la courbe représentative
de f et à chercher où elle coupe l’axe des x.
On considère donc la suite récurrente définie par une valeur u0
proche de la racine et par la relation :
Il y a deux théorèmes importants, l’un d’eux prouve que
si u0 est “assez proche” de r alors la suite un converge vers r,
malheureusement il est difficile de savoir en pratique si on est
“assez proche” de u0 pour que ce théorème s’applique. Le second
théorème donne un critère pratique facile à vérifier qui assure
la convergence, il utilise les propriétés de convexité de la fonction.
Théorème 2
Soit f une fonction de classe C2 (2 fois continument dérivable)
sur un intervalle fermé I. Soit r une racine simple de f
située à l’intérieur de I
(telle que f(r)=0 et f′(r)≠ 0). Alors il existe ε>0
tel que la suite définie par
converge vers r.Si on a |f′′| ≤ M et |1/f′| ≤ m sur un intervalle
[r−η,r+η] contenu dans I, alors on peut prendre tout réel
ε>0 tel que ε < 2/(mM) et ε ≤ η.
Démonstration : on a
un+1−r = un − r − | | =
| (un−r)f′(un)−f(un) |
|
f′(un) |
| |
En appliquant un développement de Taylor de f en un à l’ordre 2,
on obtient pour un réel θ
situé entre r et un :
0 = f(r)=f(un)+(r−un) f′(un) + (r−un)2 | | |
donc :
(un−r)f′(un)−f(un)= (un−r)2 | | |
d’où :
On commence par choisir un intervalle [r−ε,r+ε]
contenant strictement r et tel que |f′′|<M et |1/f′|<m
sur [r−ε,r+ε] (c’est toujours possible car
f′′ et 1/f′ sont continues au voisinage de r puisque f′(r)≠ 0).
Si un est dans cet intervalle, alors θ aussi donc
|un+1−r| ≤ |un−r|2 | | ≤ | |
|un−r|, |
On a |un−r| ≤ ε, on diminue si nécessaire
ε pour avoir ε < 2/(Mm), on a alors :
|un+1−r| ≤ k |un−r|, k= | | <1 |
donc d’une part un+1 est encore dans l’intervalle
[r−ε,r+ε]
ce qui permettra de refaire le même raisonnement au rang
suivant, et d’autre part
on a une convergence au moins géométrique vers r.
En fait la convergence est bien meilleure
lorsqu’on est proche de r grace au carré dans |un−r|2,
plus précisément, on montre par récurrence que
|un−r| ≤ |u0 − r|2n | ⎛
⎜
⎜
⎝ | | | ⎞
⎟
⎟
⎠ | |
il faut donc un nombre d’itérations proportionnel à ln(n)
pour atteindre une précision donnée.
Remarque : ce théorème se généralise sur ℂ et même sur ℝn.
Exemple : pour calculer √2, on écrit l’équation x2−2=0
qui a √2 comme racine simple sur I=[1/2,2],
on obtient la suite récurrente
Si on prend η=1/2, on a f′=2x et f′′=2
donc on peut prendre M=2 et m=1 car |1/f′|≤ 1 sur
[√2−1/2,√2+1/2]. On a 2/(mM)=1, on peut donc
prendre ε=1/2, la suite convergera pour tout
u0 ∈ [√2−1/2,√2+1/2].
Plus générallement, on peut calculer une racine k-ième d’un réel a
en résolvant f(x)=xk−a par la méthode de Newton.
L’inconvénient de ce théorème
est qu’il est difficile de savoir si la valeur de départ qu’on
a choisie se trouve suffisamment près d’une racine pour que
la suite converge. Pour illustrer le phénomène,
on peut par exemple colorer les points du plan
complexe en n+1 couleurs selon que la suite définie par la méthode
de Newton converge vers l’une des n racines d’un polynôme de degré
n fixé au bout de par exemple 50 itérations (la n+1-ième couleur
servant aux origines de suite qui ne semblent pas converger).
Passons maintenant à un critère très utile en pratique :
Définition 4 (convexité)
Une fonction f continument dérivable sur un intervalle I de ℝ
est dite convexe si son graphe est au-dessus de la tangente en tout point
de I.
Il existe un critère simple permettant de savoir si
une fonction de classe C2 est convexe :
Théorème 3
Si f est C2 et f′′ ≥ 0 sur I alors f est convexe.
Démonstration :
L’équation de la tangente au graphe en x0 est
Soit
g(x)=f(x)−(f(x0)+f′(x0)(x−x0)) |
on a :
g(x0)=0, g′(x)=f′(x)−f′(x0), g′(x0)=0,
g′′=f′′ ≥ 0 |
donc g′ est croissante, comme g′(x0)=0, g′ est négative
pour x<x0 et positive pour x>x0, donc g est décroissante
pour x<x0 et croissante pour x>x0. On conclut alors que
g ≥ 0 puisque g(x0)=0. Donc f est bien au-dessus
de sa tangente.
On arrive au deuxième théorème sur la méthode de Newton
Théorème 4
Si f(r)=0, f′(r)>0 et si f′′ ≥ 0 sur [r,b] alors
pour tout u0 ∈ [r,b] la suite de la méthode de Newton
est définie, décroissante, minorée par r et converge vers
r. De plus
Démonstration :
On a f′′ ≥ 0 donc si f′(r)>0 alors f′>0
sur [r,b], f est donc strictement croissante sur [r,b]
on en déduit que f>0 sur ]r,b] donc un+1 ≤ un.
Comme la courbe représentative de f est au-dessus de la tangente,
on a un+1 ≥ r (car un+1 est l’abscisse du point
d’intersection de la tangente avec l’axe des x).
La suite un est donc décroissante minorée par r, donc convergente
vers une limite l≥ r. À la limite, on a
donc l=r car f>0 sur ]r,b].
Comme (un) est décroissante, on a bien 0 ≤ un −r,
pour montrer l’autre inégalité, on applique le théorème
des accroissements finis, il existe θ ∈ [r,un] tel que
comme f(r)=0, on a
et la deuxième inégalité du théorème
en découle parce que f′ est croissante.
Variantes :
Il existe des variantes, par exemple si f′(r)<0 et f′′ ≥ 0
sur [a,r]. Si f′′ ≤ 0, on considère g=−f.
Application :
On peut calculer la valeur approchée de la
racine k-ième d’un réel a>0 en appliquant ce deuxième
théorème. En effet si
a>0, alors xk−a est 2 fois continument dérivable et
de dérivée première kxk−1 et
seconde k(k−1)xk−2 strictement positives sur ℝ+∗ (car k ≥ 2).
Il suffit donc de prendre une valeur de départ u0 plus grande que
la racine k-ième, par exemple 1+a/k (en effet
(1+a/k)k ≥ 1+k a/k=1+a).
En appliquant l’inégalité du théorème, on a :
Pour avoir une valeur approchée de (a)1/k à ε près,
on peut donc choisir comme test d’arrêt
Par exemple pour √2, le test d’arrêt serait
un2−2 ≤ 2 ε.
4 Développement de Taylor, séries entières, fonctions usuelles
Résumé: Séries entières. Calcul des fonctions
transcendantes usuelles.
Soit f une fonction indéfiniment dérivable sur un intervalle I de ℝ
et x0 ∈ I. On peut alors effectuer le développement de Taylor de f
en x0 à l’ordre n
Tn(f)(x)= f(x0) + (x−x0) f′(x0) + ... +
(x−x0)n | | |
et se demander si Tn(f) converge lorsque n tend vers
l’infini, si la limite est égale à f(x) et si on peut facilement
majorer la différence entre f(x) et Tn(f)(x). Si c’est le
cas, on pourra utiliser Tn(f)(x) comme valeur approchée de f(x).
On peut parfois répondre à ces questions simultanément en regardant le
développement de Taylor de f avec reste : il existe θ compris
entre x0 et x tel que
Rn(x) := f(x)− Tn(f)(x) = (x−x0)n+1 | | |
C’est le cas pour la fonction exponentielle que nous allons
détailler, ainsi que les fonctions sinus et cosinus.
4.1 La fonction exponentielle
Soit f(x)=exp(x) et x0=0, la dérivée n-ième de f
est exp(x), donc Rn(x)=exp(θ)xn+1/(n+1)! avec θ
compris entre 0 et x, ainsi si x est positif
|Rn(x)| ≤ ex xn+1/(n+1)! et si x est négatif,
|Rn(x)| ≤ xn+1/(n+1)!. Dans les deux cas, la limite de Rn
est 0 lorsque n tend vers l’infini, car pour n ≥ 2x, on a
on a donc pour tout x réel
Comment en déduire une valeur approchée de ex? Il suffira d’arrêter
la sommation lorsque R:=xn+1/(n+1)! si x<0 ou lorsque
R:=ex xn+1/(n+1)! si x>0 est inférieur à
l’erreur absolue souhaitée, le plus tôt étant le mieux pour des
raisons d’efficacité et pour éviter l’accumulation d’erreurs
d’arrondi.
Si on veut connaitre ex à une erreur relative ε donnée
(par exemple ε=2−53 pour stocker le résultat dans un double)
il suffit que R/ex < ε, donc si x est positif, il suffit
que xn+1/(n+1)!<ε, on peut donc arrêter la sommation
lorsque le terme suivant est plus petit que ε.
On observe que plus x est grand, plus n devra
être grand pour réaliser le test d’arrêt, ce qui est facheux
pour le temps de calcul.
De plus, le résultat final peut être petit alors que les termes
intermédiaires calculés dans la somme peuvent être grands, ce qui
provoque une perte de précision relative, par exemple si on
veut calculer e−10 ou plus générallement l’exponentielle
d’un nombre négatif de grande valeur absolue.
Exercice : combien de termes faut-il calculer dans le développement
de l’exponentielle de -10 pour que le reste soit plus petit
que 2−53 ? Quel est la valeur du plus grand terme rencontré dans
la suite ? Quelle est la perte de précision relative occasionné
par cette méthode de calcul ?
On peut utiliser les propriétés de la fonction exponentielle
pour éviter ce problème. Pour les nombres négatifs, on peut
utiliser l’équation e−x=1/ex (ne change pas l’erreur relative).
Pour les grands réels, on peut utiliser e2x=(ex)2
(multiplie par 2 l’erreur relative).
On peut aussi, si on connait une valeur approchée
de ln(2), effectuer la division euclidienne de x par ln(2)
avec reste symétrique :
x = a ln(2) + r, a ∈ ℤ, |r| ≤ | | |
puis si r est positif, on somme la série de T(f)(r), si r
est négatif, on calcule T(f)(−r) et on inverse, on applique alors :
Il faut toutefois noter que ln(2) n’étant pas connu exactement,
on commet une erreur d’arrondi absolu sur r d’ordre a η,
où η est l’erreur relative sur ln(2),
il faut donc ajouter une erreur d’arrondi relative de x/ln(2) η
qui peut devenir grande si x est grand. Puis il faut ajouter
la somme des erreurs d’arrondi due au calcul de er, que l’on
peut minimiser en utilisant la méthode de Horner pour évaluer
Tn(f)(r) (car elle commence par sommer les termes de plus haut degré
qui sont justement les plus petits termes de la somme).
Les coprocesseurs arithmétiques qui implémentent la fonction exponentielle
ont un format de représentation interne des double avec une mantisse
plus grande que celle des double (par exemple 64 bits au lieu de 53),
et une table contenant des constantes dont ln(2) avec cette précision,
le calcul de ex par cette méthode entraine donc seulement une erreur
relative d’arrondi au plus proche sur le résultat converti en double
(donc de 2−53).
Notons que en général x lui-même
a déjà été arrondi ou n’est connu qu’avec une précision relative.
Or si x>0 est connu avec une erreur relative de ε
(donc une erreur absolue de ε |x|, alors
donc on ne peut pas espérer mieux qu’une erreur relative de
eε |x|−1 sur l’exponentielle de x. Si ε x est petit
cette erreur relative (impossible à éviter, quel que soit
l’algorithme utilisé pour calculer l’exponentielle)
est d’ordre ε |x|. Si ε x est
grand alors l’erreur relative devient de l’ordre de 1, et la valeur
de l’exponentielle calculée peut être très éloignée de la valeur
réelle! Notons que pour les double, il y aura dans ce cas débordement
soit vers l’infini soit vers 0
(par exemple si x est supérieur à 709, l’exponentielle renvoie infini).
Exercice : refaire les mêmes calculs pour les fonction sinus ou cosinus.
On utilise par exemple sin(x+π)=−sin(x), sin(−x)=−sin(x),
sin(x)=cos(π/2−x) pour se ramener au calcul de sin(x)
ou de cos(x) sur [0,π/4].
sin(x)= | | (−1)n | | ,
cos(x)= | | (−1)n | | |
Cette méthode a toutefois ces limites, car il peut devenir impraticable
de calculer la dérivée n-ième d’une fonction (par exemple avec tan(x)),
et encore plus de la majorer. D’où l’intérêt de développer une théorie
des fonctions qui sont égales à leur développement de Taylor à l’infini
d’une part, et d’avoir d’autres méthodes pour majorer le reste, nous
présentons ici le cas des séries alternées.
4.2 Séries entières.
Les séries de type prendre la limite lorsque n tend vers
l’infini du développement de Taylor en x=0 sont de la forme
On peut s’intéresser plus générallement à ∑n=0∞an xn
lorsque an est un complexe quelconque,
c’est ce qu’on appelle une série entière, on peut aussi les voir comme
des polynômes généralisés.
S’il existe un point x0 tel que
|an x0n| est borné (ce sera le cas en particulier
si la série converge en x0), alors
|an xn| = |an x0n| | | | |n ≤
M | | | |n
|
la série converge donc en x si |x|<|x0| et on
peut majorer le reste de la série au rang n par
la vitesse de convergence est donc du même type que pour le
théorème du point fixe (le nombre de termes à calculer
pour trouver une valeur approchée avec k décimales
dépend linéairement k, les constantes sont d’autant
plus grandes que |x| est grand).
Théorème 5
S’il existe un rang n0, un réel M>0 et un complexe x0 tels que
pour n>n0, on ait :
alors la série converge pour |x|<|x0|
et pour n≥ n0, on a :
On en déduit qu’il existe un réel
positif R≥ 0 éventuellement égal à +∞ tel que
la série converge (la limite de la somme jusqu’à l’infini existe)
lorsque |x|<R et n’existe pas lorsque
|x|>R, ce réel est appelé rayon de convergence de la série.
Par exemple ce rayon vaut +∞ pour l’exponentielle, le sinus
ou le cosinus. Il
est égal à 1 pour la série géométrique ∑xn (car elle diverge
si |x|>1 et converge si |x|<1).
On ne peut pas dire ce qui se passe génériquement lorsqu’on
est à la limite, c’est-à-dire lorsque |x|=R (si R≠
+∞). Mais cela n’a en fait pas trop d’importance en pratique
car même si la série converge, elle converge souvent trop lentement
pour donner de bonnes approximations. En fait, la vitesse de
convergence d’une série entière de rayon R≠ +∞ est
en gros la même que celle d’une série géométrique de raison |x|/R.
Lorsque 2 séries ont un rayon de convergence non nul, alors on
peut effectuer leur somme, leur produit comme des polynômes et la
série somme/produit a un rayon de convergence au moins égal au plus
petit des 2 rayons de convergence des arguments. On peut inverser une série
entière non nulle en 0 en appliquant
et on obtient une série entière de rayon de
convergence non nul. On peut aussi composer deux séries entières
g et f en g∘ f (avec les règles de calcul de composition
des polynômes) si f(0)=0. On peut enfin dériver
et intégrer une série entière terme à terme dans son rayon de convergence.
On dit qu’une fonction est développable en série entière en 0 si
elle est égale à son développement de Taylor en 0 sommé jusqu’en l’infini
dans un disque de centre 0 et de rayon non nul. Les fonctions
exponentielle, sinus, cosinus sont donc développables en série entière en 0.
La fonction tangente également car le dénominateur cosinus est non nul en 0,
mais son rayon de convergence n’est pas l’infini et le calcul des an
est assez complexe.
La fonction (1+x)α est développable en séries entières
pour tout α ∈ ℝ avec un rayon de convergence 1 (ou l’infini
pour α entier positif).
(1+x)α= 1 + α x + | | x2 +
... + | | xn + ... |
Pour α=−1, c’est la série
géométrique de raison −x, en effet si |x|<1 :
En intégrant par rapport à x, on obtient que ln(1+x) est développable
en série entière en 0 de rayon de convergence 1 et
On peut calculer de manière analogue le développement en série entière
de arctan(x) en iintégrant celui de 1/(1+x2), de même pour arccos(x)
et arcsin(x) en intégrant celui de (1−x2)−1/2.
On peut donc calculer ln, arctan, ... par ces formules,
mais il faut répondre
à la question où arrête-t-on la somme pour obtenir une précision donnée?
Dans le cas de ln(1+x),
on pourrait répondre comme avec l’exponentielle en majorant la dérivée
n+1-ième, mais ce n’est plus faisable pour arctan, arcsin, arccos.
On va donner un autre critère qui ne nécessite pas
de calculer cette dérivée mais utilise l’alternance des signes
dans la somme.
4.3 Série alternée
Théorème 6
Soit Sn= ∑k=0n (−1)k uk la somme jusqu’au rang n d’une
série de réels tels que la suite des uk décroit à partir
d’un rang n0 et tend vers 0
lorsque k→ +∞. Alors Sn converge vers une limite
S. Si n≥ n0, la limite est comprise entre deux sommes
partielles succesives Sn et Sn+1 et le reste est majoré par
la valeur absolue du premier terme non sommé :
Démonstration :
on montre que les suites vn=S2n et wn=S2n+1 sont
adjacentes. On a
vn+1−vn= S2n+2−S2n= (−1)2n+2 u2n+2 + (−1)2n+1
u2n+1 = u2n+2−u2n+1 ≤ 0 |
donc vn est décroissante, de même wn est croissante,
et vn−wn=u2n+1 est positif et tend vers 0. On en déduit que
vn et wn convergent vers la même limite S telle que
vn>S>wn et les inégalités du théorème s’en déduisent.
Remarque
lorsqu’on utilise une suite alternée pour
trouver une valeur approchée, il faut que un tende assez
vite vers 0, sinon il y aura perte de précision sur la mantisse
lorsqu’on effectuera u2n−u2n+1. On sommera aussi les termes
par ordre décroissant pour diminuer les erreurs d’arrondi.
4.4 La fonction logarithme
Si nous voulons calculer ln(1+x) pour x ∈ [0,1[ avec une
précision ε, il suffit de calculer
pour n tel que la valeur absolue du terme suivant soit plus petit
que ε :
en effet, les signes sont alternés et la suite xk+1/k+1
décroit vers 0.
Si la suite décroit lentement vers 0, cette méthode est
mauvaise numériquement et en temps de calcul
car il y a presque compensation entre
termes successifs donc perte de précision sur la mantisse et il
y a beaucoup de termes à calculer. C’est le cas pour le logarithme,
si x est voisin de 1, il faut calculer n termes pour avoir une
précision en 1/n, par exemple 1 million de termes pour avoir
une précision de 1e−6 (sans tenir compte des erreurs d’arrondi).
Si x est proche de 1/2 il faut de l’ordre
de −ln(ε)/ln(2) termes ce qui est mieux, mais encore
relativement grand (par exemple
50 termes environ pour une précision en 1e−16,
13 termes pour 1e−4).
On a donc intérêt à se
ramener si possible à calculer la fonction en un x où la
convergence est plus rapide (donc |x| le plus petit possible).
Par exemple pour le calcul de ln(1+x) on peut :
-
utiliser la racine carrée
on observe que :
il faut toutefois faire attention à la perte de précision sur
X par rapport à x lorsque x est petit.
- utiliser l’inverse
ln(1+x)=−ln(1/(1+x))=−ln(1 + | | ) |
lorsque x est proche de 1, −x/(1+x) est proche de −x/2, on
a presque divisé par 2. Attention toutefois, on se retrouve alors
avec une série non alternée, mais on peut utiliser
(3) pour majorer le reste dans ce cas.
- trouver une valeur approchée y0 de ln(1+x) à une précision
faible, par exemple 1e−4, et utiliser la méthode de Newton pour
améliorer la précision. Soit en effet y=ln(1+x), alors
ey=1+x, on pose f(y)=ey−(1+x), on utilise la suite itérative
Comme y0 est proche à 1e−4 de y, on peut espérer avoir
une valeur approchée de y à 1e−16 en 2 itérations. Notez
que y est proche de 0, on est dans un domaine où le calcul
de ey est rapide et précis et de plus la méthode de Newton
“corrige” les erreurs intermédiaires.
Nous sommes donc en mesure de calculer précisément le logarithme
ln(1+x) pour disons |x|<1/2. Pour calculer ln sur ℝ+,
on se ramène à [1,2] en utilisant l’écriture
mantisse-exposant, puis si x∈[3/2,2] on peut en prendre la racine
carrée pour se retrouver dans l’intervalle souhaité.
On peut aussi effectuer une division par √2.
Remarquons que si x est connu à une erreur relative ε
près, comme
ln(x(1 ± ε))=ln(x) + ln(1 ± ε) |
ln(x) est connu à une erreur absolue de
|ln(1 ± ε)| ≈ ε. Si ln(x) est proche
de 0, on a une grande perte de précision relative.
Finalement, nous savons calculer ln et exp sous réserve
d’avoir dans une table la valeur de ln(2). Pour calculer
ln(2) précisément, on peut utiliser
ln(2)=−ln(1/2)=−ln(1−1/2) |
et le développement en série calculé en mode exact avec des
fractions à un ordre suffisant, on majore le reste en utilisant
que le terme général de la série ln(1+x) est borné par
M=1 en x=1, donc d’après (3) :
(on peut même obtenir 1/(n2n) car on a besoin de M uniquement
pour les termes d’ordre plus grand que n, on peut donc prendre M=1/n).
Par exemple, pour avoir ln(2) avec une mantisse de 80 bits,
on effectue une fois pour toutes avec un logiciel
de calcul formel :
a:=sum((1/2)^k/k,k=1..80)
|
puis la division en base 2 avec 81 bits de précision
iquo(numer(a)*2^81,denom(a))
Exercice : pour les fonctions trigonométriques, il faut une
méthode de calcul de π. On peut par exemple faire le calcul
de 16 arctan(1/5)−4arctan(1/239) en utilisant le développement
de la fonction arctan à un ordre suffisant.
4.5 Autres applications
On peut calculer certaines intégrales de la même manière,
par exemple
mais aussi des fonctions définies par des intégrales (cas de nombreuses
fonctions spéciales).
4.5.1 Exemple : la fonction d’erreur (error fonction, erf)
Cette fonction est définie à une constante multiplicative près par :
On peut développer en séries entières l’intégrand (rayon
de convergence +∞), puis intégrer terme
à terme, on obtient
Ce développement converge très rapidement pour |x|≤ 1.
Par contre, pour |x| grand, il faut calculer beaucoup de termes
avant que le reste soit suffisamment petit pour être négligeable,
et certains termes intermédiaires sont grands, ce qui provoque
une perte de précision qui peut rendre le résultat calculé
complètement faux. Contrairement à la fonction exponentielle,
il n’y a pas de possibilité de réduire l’argument à une plage
où la série converge vite. Il faut donc
-
soit utiliser des flottants multiprécision, avec une précision
augmentée de la quantité nécessaire pour avoir un résultat fiable
- soit, pour les grandes valeurs de x, utiliser un développement
asymptotique (en puissances de 1/x) de
ainsi que
Le développement asymptotique
s’obtient par exemple en changeant de variable
u=t2 et en effectuant des intégrations par parties répétées
en intégrant e−u et en dérivant u−1/2 et ses dérivées
successives. Ce type de développement asymptotique a la propriété
inverse du développement en 0: les termes successifs commencent
par décroitre avant de croitre et de tendre vers l’infini. Il faut
donc arrêter le développement à un rang donné (dépendant de x)
et il est impossible d’obtenir une précision meilleure pour cette
valeur de x par un développement asymptotique (on parle parfois de
développement des astronomes).
Exercice : donner une valeur approchée de f(1) à 1e−16 près.
Combien de termes faut-il calculer dans la somme pour trouver
une valeur approchée de f(7) à 1e−16 près ? Comparer la valeur
de f(7) et la valeur absolue du plus grand terme de la série,
quelle est la perte de précision relative si on effectue les
calculs en virgule flottante ? Combien de chiffres significatifs
faut-il utiliser pour assurer une précision finale de 16 chiffres
en base 10 ? Calculer le développement asymptotique en l’infini
et déterminer un encadrement de f(7) par ce développement. Combien
de termes faut-il calculer pour déterminer f(10) à 1e−16 près
par le développement asymptotique et par le développement en séries ?
Quelle est la meilleure méthode pour calculer f(10) ?
4.5.2 Recherche de solutions d’équations différentielles
On peut aussi appliquer les techniques ci-dessus pour calculer
des solutions de certaines équations différentielles dont les
solutions ne s’expriment pas à l’aide des fonctions usuelles,
on remplace dans l’équation la fonction inconnue par son
développement en séries et on cherche une relation de récurrence
entre an+1 et an. Si on arrive à montrer par exemple
qu’il y a une solution ayant un développement alternée,
ou plus générallement,
si on a une majoration |an+1/an|<C, alors le reste de la
série entière est majoré par |anxn|/(1−|Cx|) lorsque
|x|<1/C, on peut alors calculer des valeurs approchées
de la fonction solution à la précision souhaitée en utilisant
le développement en séries entières.
4.5.3 Exemple : fonctions de Bessel d’ordre entier
Soit m un entier positif fixé, on considère l’équation
différentielle
x2 y′′ + x y′ + (x2−m2)y=0 |
dont on cherche une solution série entière
y=∑k=0∞ak xk . En remplacant dans l’équation, si
x est dans le rayon de convergence de la série (rayon supposé
non nul), on obtient
| k(k−1)ak xk + | | k ak xk
+ | | (x2−m2) ak xk =0
|
soit encore
0 | = | |
| = | −m2 a0 + (1−m2)a1 x + | | [(k2−m2) ak +ak−2]xk |
|
|
Par exemple, prenons le cas m=0. On a alors a0 quelconque, a1
nul et pour k≥ 2
Donc tous les a d’indice impair sont nuls. Les pairs sont non nuls
si a0≠ 0, et ils sont de signe alterné.
Soit x fixé, on observe que pour 2k > |x|,
|a2k x2k| < |a2k−2 x2k−2| |
donc la série ∑k=0∞ak xk est alternée à partir
du rang partie entière de |x| plus un. Donc elle converge pour
tout x (le rayon de convergence de y est +∞)
et le reste de la somme jusqu’à l’ordre 2n est
inférieur en valeur absolue à :
Par exemple, pour avoir une valeur approchée à 1e−10 près de
y(x) pour a0=1 et |x|≤ 1, on calcule y=∑k=02n ak xk ,
on s’arrête au rang n tel que
|a2n+2 x2n+2| ≤ |a2n+2| ≤ 10−10 |
On remarque que :
donc n=7 convient.
Pour m ≠ 0, on peut faire un raisonnement analogue (les
calculs sont un peu plus compliqués).
On a ainsi trouvé une solution y0 de l’équation
différentielle de départ dont on peut facilement calculer
une valeur approchée (aussi facilement que par exemple la fonction sinus
pour |x| ≤ 1),
on peut alors trouver toutes les solutions de l’équation
différentielle (en posant y=y0 z et en cherchant z).
Exercice : faire de même pour les solutions de
y′′−xy=0 (fonctions de Airy).
4.6 Développements asymptotiques et séries divergentes
Un développement asymptotique est une généralisation d’un développement de Taylor, par exemple
lorsque le point de développement est en l’infini. De nombreuses fonctions ayant
une limite en l’infini admettent un développement asymptotique en l’infini, mais ces
développements sont souvent des séries qui semblent commencer par converger
mais sont divergentes. Ce type de développement s’avère néanmoins très utile lorsqu’on
n’a pas besoin d’une trop grande précision sur la valeur de la fonction.
Nous allons illustrer ce type de développement sur un exemple, la fonction
exponentielle intégrale, définie à une constante près par
On peut montrer que l’intégrale existe bien, car l’intégrand est positif et inférieur à e−t
(qui admet −e−t comme primitive, cette primitive ayant une limite en +∞).
Pour trouver le développement asymptotique de f en +∞, on effectue des intégrations
par parties répétées, en intégrant l’exponentielle et en dérivant la fraction rationnelle
f(x) | = | |
| = | |
| = | |
| = | |
| = | ... |
| = | e−x | ⎛
⎜
⎜
⎝ | | − | | + | | + ... + | | ⎞
⎟
⎟
⎠ |
− | ∫ | | | | dt |
|
| = | S(x) + R(x) |
|
où
S(x)=e−x
| ⎛
⎜
⎜
⎝ | | − | | + | | + ... + | | ⎞
⎟
⎟
⎠ | ,
R(x)=− | ∫ | | | | dt
(4) |
Le développement en séries est divergent puisque pour x>0 fixé et n tendant vers l’infini
mais si x est grand, au début la série semble converger, de manière très rapide :
On peut utiliser S(x) comme valeur approchée de f(x) pour x grand si on sait majorer
R(x) par un nombre suffisamment petit. On a
On retrouve une majoration du type de celle des séries alternées,
l’erreur relative est inférieure
à la valeur absolue du dernier terme sommé divisé par e−x/x.
Pour x fixé assez grand, il
faut donc trouver un rang n, s’il en existe un,
tel que (n+1)!/xn+1<є où
є est la précision relative que l’on s’est fixée.
Par exemple, si x≥ 100, n=11 convient pour
є=12!/10012=5e−16 (à peu près
la précision relative d’un “double”).
Ceci permet d’avoir une approximation de la fonction avec une bonne
précision et peu de calculs, mais contrairement aux séries entières,
il n’est pas possible d’améliorer cette précision de manière
arbitraire en poussant le développement plus loin, il y a une
précision maximale possible (qui dépend de x).
Ce type de développement asymptotique peut être effectué pour
d’autres fonctions du même type, par exemple
Digression: calcul approché de la constante d’Euler γ
On peut montrer que
existe (par exemple en cherchant un équivalent de un+1−un qui vaut
−1/2n2)
et on définit γ comme sa limite. Malheureusement, la convergence
est très lente et cette définition n’est pas applicable pour obtenir la valeur
de γ avec une très grande précision.
Il y a un lien entre γ et la fonction exponentielle intégrale, plus précisément
lorsque x→ 0, f(x) admet une singularité en −ln(x),
plus précisément f(x)+ln(x)
admet un développement en séries (de rayon de convergence +∞), car :
Que vaut la constante du membre de droite :
Il se trouve que C=−γ (voir plus bas une démonstration condensée) et donc :
Pour obtenir une valeur approchée de γ, il suffit donc de prendre un x assez grand
pour pouvoir calculer f(x) par son développement asymptotique à la précision requise,
puis de calculer l’intégrale du membre de droite par le développement en séries en x=0
(en utilisant une précision intermédiaire plus grande puisque ce développement en séries
va sembler diverger au début avant de converger pour n suffisamment grand).
Par exemple, on pose x=13, on calcule f(13) par (4)
avec n=13 (qui correspond au moment où le terme général
de la série est minimum puisque le rapport de deux termes successifs
est en n/x)
et une erreur absolue inférieure à e−13 13!/1314=4e−12
f(13) ≈ exp(-13)*sum((-1)^n*n!/13.^(n+1),n=0..13)
puis on remplace dans (6), avec
dont on obtient une valeur approchée,
en faisant la somme jusqu’au rang 49 (pour lequel
le terme général est de l’ordre de 1e-12
),
le reste de cette somme R50 est positif et est inférieur à
(-1)^
50*13.^
51/51/51!) qui est de l’ordre de
8e-12
evalf(sum((-1)^n*13^(n+1)/(n+1)/(n+1)!,n=0..49))
La somme argument de evalf
étant exacte, il n’y a pas de problèmes de perte de précision,
on peut aussi faire les calculs intermédiaires en arithmétique approchée,
on doit alors prendre 4 chiffres significatifs de plus
pour tenir compte de la valeur du plus grand terme
sommé dans la série, terme que l’on détermine par exemple par
seq(13.^
(n+1)/(n+1)/(n+1)!,n=0..20)
ce terme vaut 13^11/11/11!
soit 4000 environ)
Digits:=16; sum((-1)^n*13.^(n+1)/(n+1)/(n+1)!,n=0..49)
On obtient finalement comme valeur approchée de γ
-exp(-13)*sum((-1)^n*n!/13.^(n+1),n=0..13)-ln(13)+
sum((-1)^n*13^(n+1)/(n+1)/(n+1)!,n=0..49)
soit 0.577215664897
avec une erreur inférieure à 1.2e-11
.
Bien entendu, cette méthode est surtout intéressante si on veut calculer
un grand nombre de décimales de la constante d’Euler, sinon
on peut par exemple appliquer la méthode d’accélération de Richardson à
la suite convergente (5) qui définit γ
ou d’autres méthodes d’accélération (en transformant par
exemple la série en série alternée). On calcule alors
de deux manières différentes f(x) pour x plus grand (déterminé
par la précision qu’on peut obtenir par le développement
aymptotique de f).
On peut calculer π de la même manière avec le développement
en séries et asymptotique
de la fonction sinus intégral (on remplace exponentielle par sinus dans
la définition de f) et l’égalité (dont un schéma de preuve est aussi
donné plus bas)
Calcul de C (et preuve de (7)):
Pour cela on effectue une intégration par parties, cette fois en intégrant 1/t
et en dérivant l’exponentielle (moins 1 dans la première intégrale).
C | = | |
| = | [(e−t−1)ln(t)]01 + | ∫ | | ln(t) e−t dt + [e−t ln(t)]1+∞
+ | ∫ | | ln(t) e−t dt |
|
| = | |
|
Pour calculer cette intégrale, on utilise
l’égalité (qui se démontre par récurrence en faisant une
intégration par parties) :
On va à nouveau intégrer par parties,
on intègre un facteur multiplicatif 1
et on dérive l’intégrand, on simplifie, puis
on intègre t et on dérive l’autre terme, puis t2/2, etc.
C | = | [te−t ln(t)]0+∞ − | ∫ | | t e−t( | | −ln(t)) dt |
|
| = | 0 − | ∫ | | e−t dt + | ∫ | | t e−t ln(t) dt |
|
| = | −1 + [ | | e−t ln(t)]0+∞
− | ∫ | | | | e−t( | | −ln(t)) dt |
|
| = | −1 − | ∫ | | | | e−t + | ∫ | | | | e−t ln(t) dt |
|
| = | |
| = | ... |
| = | −1 − | | − ... − | | + | ∫ | | | | e−t ln(t) dt |
|
| = | −1 − | | − ... − | | + ln(n) + In |
|
|
où
In= | ∫ | | | | e−t (ln(t)−ln(n)) dt |
Pour déterminer In on fait le changement de variables t=nu
Or en faisant le même changement de variables t=nu :
n!= | ∫ | | tn e−t dt = nn+1 | ∫ | | en(ln(u)−u) du
|
Donc
Lorsque n tend vers l’infini, on peut montrer que In → 0, en effet les intégrales
sont équivalentes à leur valeur sur un petit intervalle autour de u=1, point où l’argument
de l’exponentielle est maximal,
et comme l’intégrand du numérateur a une amplitude ln(u) qui s’annule en u=1,
il devient négligeable devant le dénominateur. Finalement on a bien C=−γ.
On peut remarquer qu’en faisant le même calcul que C
mais en remplacant e−t par e−α t pour ℜ(α)>0, donne
limIn=−ln(α) (car le point critique où la dérivée
de la phase s’annule est alors 1/α). Ceci peut aussi se vérifier
pour α réel en faisant le changement de variables α t=u
∫ | | (e−α t−1) | | dt + | ∫ | | e−α t | | dt
= −γ −ln(α) |
En faisant tendre α vers −i, −ln(α)
tend vers ln(i)=iπ/2 et on obtient
∫ | | (eit−1) | | dt + | ∫ | | ei t | | dt
= −γ + i | | |
dont la partie imaginaire nous donne (7), et la
partie réelle une autre identité sur γ faisant intervenir
la fonction cosinus intégral.
5 Polynômes : arithmétique, factorisation,
interpolation
5.1 Arithmétique des polynomes: Bézout et applications
On considère les polynômes à une variable à coefficients
dans ℝ ou ℂ ou ℚ. Les algorithmes de base déjà
évoqués sont l’évaluation en un point (méthode de Horner),
l’addition, la soustraction, la multiplication et la division
euclidienne de A par B ≠ 0 :
A = B Q + R, deg(R)< deg(B)
|
A l’aide de la division euclidienne, on peut calculer le PGCD de
deux polynômes par l’algorithme d’Euclide. Nous allons présenter
l’algorithme d’Euclide étendu (ou de Bézout)
Théorème 7
Étant donnés 2 polynômes A et B, il existe deux polynômes U,
V tels que
AU + B V =
pgcd (A,B) , deg(U) < deg(B),
deg(V) < deg(A)
|
Algorithme :
On construit en fait 3 suites (Un), (Vn) et (Rn) telles que :
-
on initialise U0=1, V0=0, R0=A et U1=0, V1=1, R1=B
- on calcule les indices n+2 en fonction de n et n+1
en effectuant la division euclidienne de Rn par Rn+1
Rn = Qn Rn+1 + Rn+2,
Un+2= Un − Qn Un+1,
Vn+2 = Vn − Qn Vn+1
|
- on s’arrête au dernier reste non nul
Exemple :
A=x3−1, B=x2+1, les rangs 0 et 1 sont donnés ci-dessus.
Au rang 2, Q0 est le quotient euclidien de A par B (fonction
quo
) donc x, d’où
Puis on divise x2+1 par −x−1, quotient −x+1, donc
U3=x−1, V3=1+x(−x+1)=1+x−x2, R3=2 |
Preuve de l’algorithme :
On montre facilement par récurrence que la relation AUn + B Vn =
Rn est conservée. Comme Rn est la suite des restes, le dernier
reste non nul est bien le pgcd de A et B.
D’autre part, examinons les degrés des Vk. Supposons que
deg(A) ≥ deg(B) (sinon on échange A et B).
Au rang n=0, V0=0 donc V2=−Q0 V1, aux rangs suivants
le degré de Qn est non nul (car le degré de Rn+1 est
strictement inférieur au degré de Rn)
On montre donc par récurrence que la suite des degrés de Vn
est croissante et que :
deg (Vn+2) = deg(Qn)+deg(Vn+1) |
Comme deg(Qn)=deg(Rn)-deg(Rn+1), on en déduit que
deg(Vn+2)+deg(Rn+1)
= deg(Vn+1)+deg(Rn) = ... =
deg(V1)+deg(R0 ) = deg(A)
|
Donc si n+2 est le rang du dernier reste non nul, Vn+2=V et
degV=degA-degRn+1 est donc strictement
inférieur au degré de A (car Rn+1,
l’avant-dernier reste non nul,
est de degré plus grand ou égal à 1).
On en déduit enfin que le degré de U est strictement
inférieur au degré de B, car AU=R−BV, le degré
de BV est strictement inférieur à celui de B plus
celui de A.
L’identité de Bézout permet de résoudre plus générallement
une équation du type
où A,B,C sont trois polynômes donnés, à condition que
C soit divisible par le pgcd de A et B.
L’ensemble des solutions s’obtient à partir d’une solution
particulière U,V de Bézout, notons c=C/gcd(A,B),
on a alors
A(cU) + B(cV)=c gcd(A,B) = C |
et l’ensemble des solutions est donné par
u=cU−PB, v=cV+PA où P est un polynôme
quelconque.
Si le degré de C est plus petit que le degré de A
plus le degré de B, il existe une solution
“priviligiée”, on prend pour u le reste de la division
euclidienne de cU par B, v est alors le reste
de la division euclidienne de cV par A pour des
raisons de degré.
Exemple : si on veut résoudre
on multiplie U=x−1 et V=1+x−x2 par x2 ce qui donne une solution
l’ensemble des solutions est de la forme
et la solution priviligiée (de degrés minimaux) est
−x+1=rem(x2(x−1),x2+1),
x2−x+1=rem(x2(1+x−x2),x3−1) |
L’identité de Bézout intervient dans de nombreux problèmes
en particulier la décomposition en éléments simples d’une fraction
rationnelle. Si le dénominateur D d’une fraction se factorise
en produit de 2 facteurs D=AB premiers entre eux, alors
il existe deux polynômes u et v tels que N=Au+Bv, donc
Si de plus N/D est une fraction propre (degré de N plus petit
que celui de D), alors u/B et v/A sont encore des fractions
propres (en calculant le reste de la division euclidienne pour u et
v comme expliqué ci-dessus).
Par exemple :
| =
| (−x+1)(x3−1)+(x2−x+1)(x2+1) |
|
(x3−1)(x2+1) |
|
= | | + | |
Les applications sont diverses, citons
-
le calcul de primitive de fraction rationnelles (et tout
ce qui s’y ramène), par exemple
Puis on fait apparaitre la dérivée du dénominateur au numérateur
pour éliminer les x, 2x=(x2+1)′
pour faire le calcul complet, il faut aussi décomposer
la fraction restante (exercice!)
- Le calcul de la fonction exponentielle (à nouveau).
Au lieu d’utiliser T
le développement de Taylor en 0 par exemple à l’ordre
10, on cherche une fraction rationnelle N/D ayant le même
développement de Taylor que l’exponentielle en 0 avec
degré de N et de D majorés par 5. Pour trouver N et D
on multiplie la condition N/D=T+O(x11) par D ce qui donne
on applique l’algorithme de Bézout aux polynômes x11 et T
en s’arrêtant prématurément, lorsque le reste est de degré 5,
on montre alors que le reste est N et le coefficient de Bézout
de T est D. On peut alors montrer que l’approximation
est un peu meilleure, et nécessite moins d’opérations
(il y a une certaine symétrie entre les termes de N et D).
- le calcul de transformée de Laplace inverse de fractions
rationnelles, l’idée est la même, sauf qu’on remplace l’intégrale
par la transformée de Laplace inverse (et les formules donnant
la transformée inverse de 1/(x−p), 1/(x2+p2), p/(x2+p2)
respectivement exp(px), sin(xp)/p, cos(px))
(calcul non exigible à l’examen)
- le calcul du terme d’ordre n du développement de Taylor
en 0 d’une fraction rationnelle. On décompose, et on se ramène à des
séries dont le terme général est connu, comme (a+x)−n.
Par exemple pour connaitre le développement de 1/(x2−3x+2),
on factorise le dénominateur 1/((x−1)(x−2)), on décompose
et on développe, le terme d’ordre n est donc 1−(1/2)n+1.
Il faut néanmoins savoir factoriser un polynôme, ce dont nous
parlerons dans la section suivante.
Exercice :
Calculer l’intégrale
en utilisant l’identité de Bézout pour décomposer la fraction rationnelle.
Trouver à l’aide de cette décomposition le terme d’ordre
n du développement de Taylor de la fraction à intégrer, vérifier
avec un logiciel de calcul formel que les termes d’ordre 0 à 3 sont
corrects.
Una autre application est l’élimination dans les systèmes
polynomiaux, par exemple considérons le système de 2 équations
à 2 inconnues (intersection d’une ellipse et d’un cercle) :
x2+y2−9=0, x2+2y2−2xy−7=0 |
En calculant les coefficients de Bézout des 2 polynômes en x
x2+y2−9 et x2+2y2−2xy−7 et en multipliant au besoin par
le PPCM (plus grand commun multiple) des dénominateurs, on obtient
à droite de l’équation de Bézout un polynôme ne dépendant que
de y et qui s’annule aussi aux solutions du système, on peut alors
résoudre en y (en factorisant) puis en x. Ici par exemple
ce polynome est 5y4−32y2+4.
Cette méthode se systématise, le polynome obtenu par élimination
d’une variable est appelé résultant.
5.2 Factorisation des polynômes
Soit P un polynôme de degré non nul. Factoriser P n’a pas
une signification unique, tout dépend d’une part si on veut
une factorisation exacte ou approchée, et d’autre part quels seront
les types des coefficients de la factorisation (complexes, réels,
entiers).
5.2.1 Multiplicité des racines.
On dit que r est une racine de multiplicité k de P si
P(x)=(x−r)k Q et Q(r)≠ 0.
En faisant le développement de Taylor de P en r à l’ordre
degré de P, on voit que cela équivaut à :
P(r)=P′(r)=...=P[k−1](r)=0, P[k](r) ≠ 0 |
En particulier si P(r)=0, on peut factoriser P par X−r.
On peut donc détecter les racines de multiplicité supérieure
à 1 en cherchant un facteur commun à P et P′, en effet x−r
divisera P et P′.
Théorème 8
Si P et P′ sont premiers entre eux (pgcd = 1), alors les
racines de P sont simples (de multiplicité 1).
Il existe un algorithme (dû à Yun) qui permet d’écrire un
polynome quelconque comme produit de polynômes dont les racines
sont simples en effectuant uniquement des calculs de PGCD de
polynomes.
yun(P):= {
local W,Y,G,res;
W:=P;
Y:=diff(W,x);
res:=NULL;
while(true){
if (Y==0) {
return res[1..size(res)-1],W;
};
G:=gcd(Y,W);
res:=res,G;
W:=normal(W/G);
Y:=normal(Y/G-diff(W,x));
};
}
L’instruction squarefree
ou équivalente
de votre logiciel de calcul formel effectue cette décomposition.
5.2.2 Factorisation dans ℂ.
Reste maintenant à trouver des racines!
On a le :
Théorème 9 (d’Alembert)
Soit P un polynome de degré non nul, alors P admet au moins
une racine complexe.
On peut alors factoriser P par X−r si r est la racine, et
recommencer avec le quotient, d’où le corollaire.
Théorème 10
Soit P un polynome de degré n non nul, alors P admet n
racines complexes (comptées avec multiplicité) x1,...,xn,
on a donc :
où an est le coefficient dominant de P.
Démonstration du théorème de d’Alembert :
On va montrer que le minimum de la valeur absolue de P est atteint
en un nombre complexe puis que ce minimum est forcément nul.
Soit
P(x)=an xn + ... + a0, an ≠ 0 |
Lorsque |x| tend vers l’infini, |P(x)| tend vers l’infini, en
effet
P(x)= an xn ( 1 + | | | | + ...
+ | | | | ) ≈|x|→ ∞ an xn |
plus précisément il existe R>0 tel que si
|x|>R alors |P(x)|>|an| |x|n/2. Quitte à augmenter R on
peut donc supposer que |P(x)|>|P(0)| si |x|>R, donc il existe
un complexe x0 qui réalise le minimum de |P| sur ℂ (ce
minimum est en fait le minimum pour |x|≤ R). On va montrer par
l’absurde que ce minimum est nul (donc que x0 est la racine
cherchée). Supposons donc que P(x0) ≠ 0. On fait le
développement de Taylor de P en x0 à l’ordre n=degré de P, donc
le développement n’a pas de reste :
P(x) − P(x0) = (x−x0) P′(x0) + .. + (x−x0)n
| |
Comme P n’est pas constant, l’un des termes du membre de droite est
non nul, soit k l’indice du premier terme non nul, on a alors :
P(x) = P(x0) + (x−x0)k
| | + o((x−x0)k)
|
Comme P(x0) ≠ 0, on peut le factoriser en :
P(x) = P(x0)( 1 + (x−x0)k | | +
o((x−x0)k)
|
on pose alors x=x0+t w où w est une racine k-ième (cela
existe dans ℂ) de
on a alors :
lorsque t est positif, suffisamment petit, on a 0 < 1−tk+o(tk+1
< 1, donc |P(x)|<|P(x0)|, ce qui est absurde (x0 réalisant
le minimum de P sur ℂ).
Remarque :
Si on développe la relation (8), on obtient
des relations entre les coefficients du polynome et les racines, par
exemple :
an−1=an ∑j=1n (−xj), ...,
a0 = an | | (−xj),
|
5.2.3 Calcul approché des racines complexes simples
La section précédente nous a montré qu’on pouvait
se ramener à la recherche de racines simples, ce qui
donne envie d’essayer la méthode de Newton. On a malheureusement
rarement la possibilité de pouvoir démontrer qu’à partir d’une valeur
initiale donnée, la méthode de Newton converge,
parce que les racines peuvent être complexes, et même si elles
sont réelles, on n’a pas forcément de résultat sur la convexité
du polynôme (cf. cependant une application des suites de
Sturm dans la section suivante qui permet de connaitre le signe
de P′′ sur un intervalle sans le factoriser).
On effectue donc souvent des itérations de Newton, en partant de
0.0, en espérant s’approcher suffisamment d’une racine pour que
le théorème de convergence théorique s’applique. On se fixe
un nombre maximal d’itérations, si on le dépasse on prend alors
une valeur initiale aléatoire complexe et on recommence.
Une fois une racine déterminée, on l’élimine en calculant
le quotient euclidien Q de P par X−r (par l’algorithme de Horner),
puis on calcule les racines du quotient Q (qui sont des racines de P).
Un problème pratique apparait alors, c’est que r n’est pas exact
donc le quotient Q non plus, au fur et à mesure du calcul des
racines de P, on perd de plus en plus de précision.
Il existe une amélioration simple, si r′ est une racine
approchée de Q, alors elle est racine approchée de P
et on a toutes les chances qu’elle soit suffisamment proche
d’une racine de P pour que le théorème s’applique, on
effectue alors 1 ou 2 itérations de Newton avec r′ mais pour P
(et non Q) afin d’améliorer sa précision comme racine de P.
5.2.4 Factorisation dans ℝ, localisation des racines
Pour factoriser un polynôme à coefficients réels, on commence
par le factoriser dans ℂ. On observe ensuite que si r est
une racine complexe non réelle de P, alors son conjugué
l’est aussi (il suffit de prendre le conjugue de la relation P(r)=0)
et avec la même multiplicité (les dérivées successives
de P étant aussi à coefficients réels). On regroupe alors
les facteurs correspondant à des racines complexes conjuguées :
(X−r)(X− | | ) = X2 − (r+ | | )X+r | |
= X2 − 2 ℜ(r) X + |r|2 |
Finalement, on a le :
Théorème 11
La factorisation d’un polynôme à coefficients réels
sur ℝ donne un produit de facteurs de degré 1 (correspondant
à des racines réelles) et de degré 2 (correspondant à
des paires de racines complexes conjuguées)
Il existe un algorithme utilisant l’algorithme de calcul du PGCD
de P et P′ qui permet de déterminer le nombre de
racines réelles d’un polynôme P sans racine multiple
sur ℝ ou dans un intervalle de R.
Théorème 12
On définit la suite de polynômes A0=P, A1=P′, ..., Ak,0
en prenant l’opposé du reste de la division euclidienne des deux
précdents :
Ai = Ai+1 Qi+2 − Ai+2
(9) |
Soit Ak, le dernier reste non nul, c’est un polynôme constant puisque
P n’a pas de racine multiple.
On définit s(a) comme étant le nombre de changements de signes
de la suite Ai(a) en ignorant les 0.
Alors le nombre de racines réelles de A0=P sur l’intervalle
]a,b] est égal à s(a)−s(b).
Exemple :
Quel est le nombre de racines réelles de P=x3+x+1
sur [−2,2]? sur [0,2]?
On a donc
A0=x3+x+1, A1=P′=3x2+1,
A2=−rem(A0,A1,x)=− | | x−1,
A3=− | | |
En x=−2 on obtient la suite −9,13,1/3,−31/4 (2 changements de
signe), en x=2 on obtient la suite 11,13,−7/3,−31/4 (1 changement
de signe), il y a donc 1 racine réelle entre -2 et 2. En x=0
on obtient la suite 1,1,−1,−31/4 (1 changement de signe) donc la
racine réelle est entre -2 et 0.
Preuve
On considére la suite des signes en un point : elle ne peut contenir
deux 0 successifs (sinon toute la suite vaudrait 0 en ce point en appliquant
(9), or Ak est constant non nul). Elle ne peut pas
non plus contenir ...,+,0,+,... ni ...,-,0,-,...
à cause de la convention de signe
sur les restes de (9). Donc si b est une racine
de Ai pour 0<i<k, alors en b on a soit ...,-,0,+,... soit
...,+,0,-,... . Regardons le premier cas (le deuxième cas se traite
de manière analogue), pour x proche de b, on va avoir
...,-,-,+,... ou ...,-,+,+,... dans les 2 cas la contribution au
nombre de changements de signe est constant (égal à 1).
Comme Ak est constant, seules les racines de A0=P
sont susceptibles de faire varier s. Comme A1=P′, le sens de
variations de A0 au voisinage d’une racine de A0 est déterminé
par le signe de A1, donc lorsque x augmente en traversant une
racine r de P, il y a deux possibilités soit P est croissant
et on passe de -,+,... à +,+,..., soit P est décroissant et
on passe +,-,... à -,-,.... Dans les deux cas, on diminue
s d’une unité.
Application :
Si il n’existe pas de racines réelles dans un intervalle donné,
alors le polynôme garde un signe constant sur cet intervalle,
que l’on peut déterminer en calculant la valeur du polynôme
en un point de cet intervalle. On peut ainsi établir
dans certains cas que la méthode de Newton pour trouver une racine
d’un polynôme convergera.
Par exemple pour le polynôme P=3x5−10x3+30x2−x−45, on a
P′′=60(x3−x+1), est positif sur ℝ+ (exercice : calculer
la suite de Sturm correspondante pour le vérifier). On vérifie
que P(1)<0 et P′(1)>0 donc il existe une racine r>1 telle
que P′(r)>0, toute valeur de départ de Newton supérieure
à r assure la convergence.
Remarque :
On peut aussi déterminer les racines réelles d’un polynôme
à coefficients rationnels
en faisant uniquement des calculs exacts par dichotomie.
Cette méthode de localisation des racines réelles se
généralise d’ailleurs au cas complexe. On peut ainsi déterminer les
racines complexes d’un polynôme à coefficients complexes
rationnels de manière déterministe à la précision voulue
(cf. Eisermann).
5.2.5 Factorisation exacte
Soit P un polynôme à coefficients entiers. Lorsqu’on demande
à un logiciel de calcul formel de factoriser P, par défaut
il ne calcule pas les racines complexes approchées, mais renvoie
une factorisation exacte, sous forme de produit de facteurs
à coefficients entiers. Les degrés des facteurs peuvent
être plus grand que 2. Par exemple x4+x+1 ne peut pas
être factorisé en produit de polynômes à coefficients entiers
(bien qu’il ait 2 facteurs de degré 2 dans ℝ et 4 de degré 1
dans ℂ).
Commencons par une méthode simple de calcul des racines rationnelles
de P (les racines rationnelles correspondent à des facteurs entiers de
degré 1 de la forme qX−p de P).
Soit x=p/q une
racine rationnelle écrite sous forme de fraction
irréductible de P=an Xn+...+a0, on a alors
0 = P( | | ) = an | |
+ an−1 | | + ... + a0
= | an pn + an−1 pn−1q+...+a1 p qn−1+ a0 qn |
|
qn |
|
Donc :
p(an pn−1 + an−1 pn−1q+...+a1 qn−1) = − a0 qn |
et p divise donc a0 qn. Comme p/q est irréductible,
cela entraine que p divise a0. De même q divise an.
Il suffit donc de tester quelles sont les racines de P parmi
toutes les fractions irréductibles
de la forme un diviseur de a0 sur un diviseur de an (attention
à ne pas oublier les diviseurs négatifs!).
Exemple: racines rationnelles de 2x2+3x+1=0. On a p divise 1
donc vaut 1 ou -1, q divise 2 donc vaut 1 ou 2. On teste
donc 1, -1, 1/2, -1/2. On obtient ici la factorisation complète
du polynome (les racines sont -1 et -1/2)
Remarques :
-
Pour un polynome aléatoire, on ne trouvera aucune racine
rationnelle.
- Cette méthode n’est pas très efficace, car factoriser
un entier peut être long, le nombre de tests peut être très
grand (si an et a0 ont beaucoup de facteurs), les logiciels
de calcul formel utilisent des méthodes appelées p-adiques pour trouver
les racines rationnelles d’un polynome (on calcule d’abord les racines
de P modulo p puis modulo pk pour k assez grand).
On pourrait aussi penser à calculer les racines complexes approchées
et voir si en multipliant par an on est proche d’un entier,
on testerait alors le rationnel correspondant.
Pour déterminer les facteurs à coefficients entiers de plus
grand degré, il n’existe pas de méthode aussi simple.
On peut calculer des valeurs approchées des racines complexes
et essayer de créer des paquets de racines complexes, puis tester
si an∏r ∈ paquet (X−r) est à coefficients entier
(aux erreurs d’arrondi près). Par exemple si on calcule les
racines complexes approchées de x6+2x3−x2+1, on pourra composer
un facteur de degré 3 à coefficients entiers en rassemblant
les racines de x3+x+1.
Les logiciels de calcul formel utilisent des algorithmes modulaires
et p-adiques (consistant à factoriser le polynome modulo p).
5.3 Approximation polynomiale
Étant donné la facilité de manipulation qu’apportent les
polynomes, on peut chercher à approcher une fonction par un
polynôme. La méthode la plus naturelle consiste à chercher
un polynôme de degré le plus petit possible
égal à la fonction en certains points x0,...,xn
et à trouver une majoration de la différence entre la fonction
et le polynôme.
Le polynome interpolateur de Lagrange répond à cette question.
Soit donc x0,...,xn des réels distincts et y0,...,yn
les valeurs de la fonction à approcher en ces points (on posera
yj=f(xj) pour approcher la fonction f). On cherche
donc P tel que P(xj)=yi pour j ∈ [0,n].
Commencons par voir s’il y a beaucoup de solutions. Soit P et Q
deux solutions distinctes du problème, alors P−Q est non nul
et va s’annuler en x0, ...,xn donc possède n+1 racines donc
est de degré n+1 au moins. Réciproquement, si on ajoute
à P un multiple du polynome A=∏j=0n (X−xj), on obtient
une autre solution. Toutes les solutions se déduisent donc
d’une solution particulière en y ajoutant un polynome de degré
au moins n+1 multiple de A.
Nous allons maintenant construire
une solution particulière de degré au plus n.
Si n=0, on prend P=x0 constant. On procède ensuite par
récurrence. Pour construire le polynôme correspondant
à x0,...,xn+1 on part du polynoôme Pn correspondant à
x0,...,xn et on lui ajoute un multiple réel de A
Ainsi on a toujours Pn+1(xj)=yj pour j=0,..n, on calcule
maintenant αn+1 pour que Pn+1(xn+1)=yn+1.
En remplacant avec l’expression de Pn+1 ci-dessus, on obtient
Pn(xn+1)+ αn+1 | | (xn+1−xj) = yn+1 |
Comme tous les xj sont distincts, il existe une solution unique :
On a donc prouvé le :
Théorème 13
Soit n+1 réels distincts x0,...,xn et n+1
réels quelconques y0,...,yn.
Il existe un unique polynôme P de degré inférieur ou égal à
n, appelé polynome de Lagrange, tel que :
Exemple : déterminons le polynome de degré inférieur ou égal
à 2 tel que P(0)=1, P(1)=2, P(2)=1. On commence par P0=1.
Puis on pose P1=P0+ α1X=1+ α1X.
Comme P(1)=2=1+ α1 on en tire α1=1
donc P1=1+X. Puis on pose P2=P1+ α2X(X−1), on a
P2(2)=3+2 α2=1
donc α2=−1, finalement P2=1+X−X(X−1).
Reste à estimer l’écart entre une fonction et son polynome
interpolateur, on a le :
Théorème 14
Soit f une fonction n+1 fois dérivable sur un intervalle I=[a,b]
de ℝ, x0,...,xn des réels distincts de I.
Soit P le polynome de Lagrange donné par les xj et yj=f(xj).
Pour tout réel x ∈ I,
il existe un réel ξx ∈ [a,b] (qui dépend de x) tel
que :
Ainsi l’erreur commise dépend d’une majoration de la taille
de la dérivée n+1-ième sur l’intervalle, mais aussi
de la disposition des points xj par rapport à x. Par exemple
si les points xj sont équidistribués, le terme
|∏j=0n(x−xj)| sera plus grand près du bord de I qu’au
centre de I.
Preuve du théorème : Si x est l’un des xj l’égalité est vraie. Soit
on considère maintenant la fonction :
g(t)=f(t)−P(t) − C | | (t−xj) |
elle s’annule en xj pour j variant de 0 à n ainsi qu’en x
suite au choix de la constante C, donc g s’annule au moins n+2 fois
sur l’intervalle contenant les xj et x, donc g′ s’annule au moins
n+1 fois sur ce même intervalle, donc g′′ s’annule au moins
n fois, etc. et finalement g[n+1] s’annule une fois
au moins sur cet intervalle. Or
g[n+1] = f[n+1] − C (n+1)! |
car P est de degré inférieur ou égal à n
et ∏j=0n(x−xj) − xn+1 est de degré
inférieur ou égal à n. Donc il existe bien un réel ξx dans
l’intervalle contenant les xj et x tel que
Calcul efficace du polynôme de Lagrange.
Avec la méthode de calcul précédent, on remarque que le
polynôme de Lagrange peut s’écrire à la Horner sous la forme :
P(x) | = | α0 + α1 (x−x0) + ... + αn
(x−x0)...(x−xn−1) |
| = | α0 + (x−x0)( α1 + (x−x1)(α2 + ... +
(x−xn−2)(αn−1+(x−xn−1) αn)...)) |
|
ce qui permet de le calculer rapidement une fois les αi
connus.
On observe que
On va voir que les αk peuvent aussi se mettre sous forme
d’une différence.
On définit les différences divisées d’ordre n par récurrence
f[xi]=f(xi), f[xi,...,xk+i+1]=
| f[xi+1,...,xk+i+1]−f[xi,...,xk+i] |
|
xk+i+1−xi |
| |
On va montrer que αk=f[x0,...,xk].
C’est vrai au rang 0, il suffit donc de le montrer au rang k+1 en
l’admettant au rang k. Pour cela on observe qu’on peut construire
le polynôme d’interpolation en x0,...,xk+1 à partir des polynômes
d’interpolation Pk en x0,...,xk et Qk en x1,...,xk+1
par la formule :
Pk+1(x)= | (xk+1−x)Pk + (x−x0)Qk |
|
xk+1−x0 |
|
en effet on vérifie que Pk+1(xi)=f(xi) pour i∈ [1,k] car
Pk(xi)=f(xi)=Qk(xi),
et pour i=0 et i=k+1, on a aussi Pk+1(x0)=f(x0) et
Pk+1(xk+1)=f(xk+1).
Or αk+1 est le coefficient dominant de Pk+1 donc
c’est la différence du coefficient dominant de Qk et de Pk
divisée par xk+1−x0, c’est-à-dire la définition de
f[x0,...,xk+1] en fonction de f[x1,...,xk+1] et f[x0,...,xk].
Exemple : on reprend P(0)=1, P(1)=2, P(2)=1. On a
xi | f[xi] | f[xi,xi+1] | f[x0,x1,x2] |
0 | | | |
| | | |
1 | 2 | | |
| | (1−2)/(2−1)=−1 | |
2 | 1 | | |
|
donc P(x)=1+(x−0)(1+(x−1)(−1))=1+x(2−x).
On peut naturellement utiliser l’ordre que l’on souhaite pour les
xi, en observant que le coefficient dominant de P ne dépend pas de
cet ordre, on en déduit que f[x0,...,xk] est indépendant de
l’ordre des xi, on peut donc à partir du tableau ci-dessus
écrire P par exemple avec l’ordre 2,1,0, sous la forme
P(x)=1+(x−2)(−1+(x−1)(−1))=1+(x−2)(−x) |
6 Intégration numérique
Les fractions rationnelles admettent une primitive que l’on calcule
en décomposant la fraction avec Bézout comme expliqué précédemment.
Mais elles font figure d’exceptions,
la plupart des fonctions n’admettent pas de primitives qui s’expriment
à l’aide des fonctions usuelles. Pour calculer une intégrale,on
revient donc à la définition d’aire sous la courbe, aire que
l’on approche, en utilisant par exemple un polynome de Lagrange.
Le principe est donc le suivant : on découpe l’intervalle d’intégration
en subdivisions [a,b]=[a,a+h] + [a+h,a+2h]+...[a+(n−1)h,a+nh=b, où
h=(b−a)/n est le pas de la subdivision, et
sur chaque subdivision, on approche l’aire sous la courbe.
6.1 Les rectangles et les trapèzes
Sur une subdivision [α,β], on approche la fonction par un segment.
Pour les rectangles, il s’agit d’une horizontale : on peut prendre
f(α), f(β) (rectangle à droite et gauche)
ou f((α+β)/2) (point milieu), pour les trapèzes on utilise
le segment reliant [α,f(α)] à [β,f(β)].
Exemple : calcul de la valeur approchée de ∫01 t3 dt
(on en connait la valeur exacte 1/4=0.25) par ces méthodes en subdivisant
[0,1] en 10 subdivisions (pas h=1/10), donc α=j/10 et β=(j+1)/10
pour j variant de 0 à 9.
Pour les rectangles à gauche, on obtient sur une subdivision
f(α)=(j/10)3 que l’on multiplie par la longueur de la subdivision
soit h=1/10 :
Pour les rectangles à droite, on obtient
Pour le point milieu f((α+β)/2)=f((j/10+(j+1)/10)/2)=f(j/10+1/20)
Enfin pour les trapèzes, l’aire du trapèze délimité par l’axe des x,
les verticales y=α, y=β et les points sur ces verticales
d’ordonnées respectives f(α) et f(β) vaut
donc
| | | | ⎛
⎜
⎜
⎝ | ( | | )3 +( | | )3
| ⎞
⎟
⎟
⎠ | =
| | = 0.2525 |
Dans la somme des trapèzes, on voit que chaque terme apparait deux fois
sauf le premier et le dernier.
Plus générallement, les formules sont donc les suivantes :
| rectangle gauche | = | | (11) |
rectangle droit | = | | (12) |
point milieu | = | | (13) |
trapezes | = |
h | ⎛
⎜
⎜
⎝ | | + | | f(a+jh) | ⎞
⎟
⎟
⎠ |
|
| (14) |
|
où h=(b−a)/n est le pas de la subdivision, n le nombre de subdivisions.
On observe sur l’exemple que le point milieu et les trapèzes donnent
une bien meilleure précision que les rectangles. Plus généralement,
la précision de l’approximation n’est pas la même selon le choix
de méthode.
Ainsi pour les rectangles à gauche (le résultat est le même
à droite), si f est continument dérivable, de dérivée majorée
par une constante M1 sur [a,b], en faisant un
développement de Taylor de f en α, on obtient
| | ∫ | | f(t) dt − | ∫ | | f(α) dt |
= | | ∫ | | f′(θt)(t−α) dt |
≤ M1 | ∫ | | (t−α) dt = M1 | |
Ainsi dans l’exemple, on a M1=3, l’erreur est donc majorée par 0.015
sur une subdivision, donc par 0.15 sur les 10 subdivisions.
Pour le point milieu, on fait le développement en (α+β)/2 à l’ordre
2, en supposant que f est deux fois continument dérivable :
Dans l’exemple, on a M2=6, donc l’erreur sur une subdivision est
majorée par 0.25e−3, donc sur 10 subdivisions par 0.25e−2=0.0025.
Pour les trapèzes, la fonction g dont le graphe est le segment reliant
[α,f(α)] à [β,f(β)] est f(α)+(t−α)/(β−α)f(β),
c’est en fait un polynome de Lagrange, si f est deux fois continument
dérivable, on peut donc majorer la différence
entre f et g en utilisant (10), on intègre la valeur
absolue ce qui donne
| | ∫ | | f(t) dt − | ∫ | | g(t) dt |
≤ | ∫ | | | | | (x−α)(x−β)|
≤ M2 | | |
où M2 est un majorant de |f′′| sur [a,b].
Lorsqu’on calcule l’intégrale sur [a,b] par une de ces méthodes,
on fait la somme sur n=(b−a)/h subdivisions de longueur
β−α=h, on obtient
donc une majoration de l’erreur commise sur l’intégrale :
-
pour les rectangles à droite ou gauche nM1h2/2=M1 h (b−a)/2
- pour le point milieu M2 h2 (b−a)/24
- pour les trapèzes M2h2 (b−a)/12.
Lorsque h tend vers 0, l’erreur tend vers 0, mais pas à la même vitesse,
plus rapidement pour les trapèzes et le point milieu
que pour les rectangles. Plus on
approche précisément la fonction sur une subdivision, plus la
puissance de h va être grande, plus la convergence sera rapide
lorsque h sera petit, avec toutefois une contrainte fixée par la
valeur de Mk, borne sur la dérivée k-ième de f (plus
k est grand, plus Mk est grand en général). Nous allons voir dans la suite
comment se comporte cette puissance de h en fonction de la facon
dont on approche f.
6.2 Ordre d’une méthode
On appelle méthode d’intégration l’écriture d’une approximation
de l’intégrale sur une subdivision sous la forme
∫ | | f(t) dt ≈ I(f)= | | wj f(yj) |
où les yj sont dans l’intervalle [α,β], par exemple équirépartis
sur [α,β]. On utilise aussi la définition :
∫ | | f(t) dt ≈ I(f)=
(β−α) | | wj f(yj) |
On prend toujours ∑j wj=β−α (ou ∑j wj=1) pour
que la méthode donne le résultat exact si la fonction est constante.
On dit qu’une méthode d’intégration est d’ordre n si il y a égalité
ci-dessus pour tous les polynômes de degré inférieur ou égal à n
et non égalité pour un polynôme de degré n+1.
Par exemple, les rectangles à droite et gauche sont d’ordre 0,
le point milieu et les trapèzes sont d’ordre 1. Plus
générallement, si on approche f par son polynôme
d’interpolation de Lagrange en n+1 points (donc par un polynôme
de degré inférieur ou égal à n), on obtient une méthode
d’intégration d’ordre au moins n.
Si une méthode est d’ordre n avec des wj≥ 0 et
si f est n+1 fois continument dérivable,
alors sur une subdivision, on a :
| | ∫ | | f−I(f)| ≤ Mn+1 | |
( | | +1)
(15) |
En effet, on fait le développement de Taylor de f par exemple
en α à l’ordre n
f(t) | = | |
Tn(f) | = | f(α)+(t−α)f′(α)+...+
| | f[n](α) |
|
|
Donc
| | ∫ | | f− | ∫ | | Tn(f)|
≤ | ∫ | | | | |f[n+1](θt)|
≤ | ⎡
⎢
⎢
⎣ | Mn+1 | | | ⎤
⎥
⎥
⎦ | |
De plus,
|I(f) −I(Tn(f))| =|I | ⎛
⎜
⎜
⎝ | f[n+1](θt)
| | | ⎞
⎟
⎟
⎠ | | |
| ≤ | |
| ≤ | |
|
Donc comme la méthode est exacte pour Tn(f), on en déduit que
| = | | | ∫ | | f− | ∫ | | Tn(f)+I(Tn(f))− I(f)| |
|
| ≤ | | | ∫ | | f− | ∫ | | Tn(f)|+|I(Tn(f))− I(f)| |
|
| ≤ | |
|
Si les wj≥ 0, alors ∑j=1k |wj|=∑j=1k wj=β−α
et on obtient finalement (15)
On remarque qu’on peut améliorer la valeur de la constante
en faisant tous les développement de Taylor
en (α+β)/2 au lieu de α,
Après sommation sur les n subdivisions, on obtient que :
Théorème 15
Pour une méthode d’ordre n à coefficients positifs et une fonction f
n+1 fois continument dérivable
| | ∫ | | f−I(f)| ≤ Mn+1 | | (b−a)
( | | +1) |
On observe que cette majoration a la bonne puissance de h sur
les exemples déja traités, mais pas forcément le meilleur coefficient
possible, parce que nous avons traité le cas général d’une
méthode d’ordre n.
6.3 Simpson
Il s’agit de la méthode obtenue en approchant la fonction
sur la subdivision [α,β] par son polynome de Lagrange
aux points α,(α+β)/2,β. On calcule l’intégrale
par exemple avec un logiciel de calcul formel, avec Xcas :
factor(int(lagrange([a,(a+b)/2,b],[fa,fm,fb]),x=a..b))
qui donne la formule sur une subdivision
I(f) = | | (f(α)+4f( | | ) + f(β)) |
et sur [a,b] :
I(f) = | | | ⎛
⎜
⎜
⎝ | f(a)+f(b)+ 4 | | f(a+jh+ | | )
+ 2 | | f(a+jh) | ⎞
⎟
⎟
⎠ |
(16) |
Si on intègre t3 sur [0,1] en 1 subdivision par cette méthode,
on obtient
c’est-à-dire le résultat exact, ceci est aussi vérifié pour f polynome
de degré inférieur ou égal à 2 puisque l’approximation de Lagrange
de f est alors égale à f. On en déduit que la méthode de Simpson
est d’ordre 3 (pas plus car la méthode
de Simpson appliquée à l’intégrale de
t4 sur [0,1] n’est pas exacte). On peut même améliorer (cf. par exemple
Demailly) la constante générale de la section précédente pour la majoration
de l’erreur en :
Cette méthode nécessite 2n+1 évaluations de f (le calcul
de f est un point étant presque toujours
l’opération la plus couteuse en temps d’une
méthode de quadrature), au lieu de n pour les rectangles
et le point milieu et n+1 pour les trapèzes. Mais on a une majoration
en h4 au lieu de h2 donc le “rapport qualité-prix” de la méthode
de Simpson est meilleur, on l’utilise donc plutot que les
méthodes précédentes sauf si f n’a pas la régularité
suffisante (ou si M4 est trop grand).
6.4 Newton-Cotes
On peut généraliser l’idée précédente, découper la subdivision
[α,β] en n parts égales et utiliser le polynôme d’interpolation
en ces n+1 points x0=α, x1, ..., xn=β.
Ce sont les méthodes de Newton-Cotes,
qui sont d’ordre n au moins. Comme le polynôme d’interpolation
dépend linéairement des ordonnées, cette méthode est bien
de la forme :
De plus les wj sont universels (ils ne dépendent pas de
la subdivision), parce qu’on peut faire
le changement de variables x=α+t(β−α) dans l’intégrale
et le polynôme d’interpolation et donc se ramener à [0,1].
Exemple : on prend le polynôme d’interpolation en 5 points
équidistribués sur une subdivision [a,b] (méthode de Boole).
Pour calculer les
wj, on se ramène à [0,1], puis on tape
int(lagrange(seq(j/4,j,0,4),[f0,f1,f2,f3,f4]),x=0..1)
et on lit les coefficients de f0
à f4
qui sont les w0 à w4: 7/90, 32/90, 12/90, 32/90, 7/90.
La méthode est d’ordre au moins 4 par construction, mais on vérifie
qu’elle est en fait d’ordre 5 (exercice), la majoration de l’erreur
d’une méthode d’ordre 5 est
| | ∫ | | f −I(f)| ≤ | | (1+ | | ) h6 (b−a) |
elle peut être améliorée pour cette méthode précise en
En pratique, on ne les utilise pas très souvent, car d’une part
pour n≥ 8, les wj ne sont pas tous positifs, et
d’autre part, parce que la constante Mn devient trop
grande. On préfère utiliser la méthode de Simpson en utilisant
un pas plus petit.
Il existe aussi d’autres méthodes, par exemple les quadratures de Gauss
(on choisit d’interpoler en utilisant des points non équirépartis
tels que l’ordre de la méthode soit le plus grand possible)
ou la méthode de Romberg qui est une méthode d’accélération
de convergence basée sur la méthode des trapèzes (on prend
la méthode des trapèzes en 1 subdivision de [a,b], puis
2, puis 22, ..., et on élimine les puissances de h
du reste ∫f−I(f) en utilisant un théorème d’Euler-Mac Laurin
qui montre que le développement asymptotique de
l’erreur en fonction de h ne contient que des puissances paires
de h). De plus, on peut être amené à faire varier le pas h
en fonction de la plus ou moins grande régularité de la fonction.
6.5 En résumé
Intégration sur [a,b], h pas d’une subdivision, Mk majorant
de la dérivée k-ième de la fonction sur [a,b]
| formule | Lagrange degré | ordre | erreur |
rectangles | (11), (12) | 0 | 0 | M1 h (b−a)/2 |
point milieu | (13) | 0 | 1 | M2 h2 (b−a)/24 |
trapèzes | (14) | 1 | 1 | M2 h2 (b−a)/12 |
Simpson | (16) | 2 | 3 | M4 h4 (b−a)/2880 |
7 Algèbre linéaire
7.1 Le pivot de Gauss
Cet algorithme permet de créer des zéros en effectuant des manipulations
réversibles sur les lignes d’une matrice. Ces lignes peuvent
représenter les coefficients d’un système linéaire, on obtient
alors un système linéaire équivalent, ou les coordonnées
d’un système de vecteur, on obtient alors les coordonnées
d’un système de vecteur engendrant le même sous-espace vectoriel. On
peut également représenter 2 matrices A et B reliés par une relation
Ax=B, cette relation reste alors vraie au cours et donc
après la réduction.
7.1.1 L’algorithme
L’algorithme est le suivant:
-
on initialise c=1 et l=1,
c désigne le numéro de colonne c à réduire, et l le
numéro de ligne à partir duquel on cherche un “pivot”
(au début l et c valent donc 1, en général
les 2 augmentent de 1 à chaque itération)
- Si c ou l est plus grand que le nombre de colonnes ou
de lignes on arrête.
- Si la colonne c n’a
que des coefficients nuls à partir de la ligne l, on incrémente
le numéro de colonne c de 1 et on passe à l’étape 2. Sinon,
on cherche la ligne dont le coefficient est en valeur absolue
le plus grand possible (en calcul approché) ou le plus simple
possible (en calcul exact), on échange cette ligne avec la ligne l.
Puis on effectue pour toutes les lignes sauf l ou pour toutes
les lignes à partir de l+1 (selon qu’il s’agit d’une réduction
de Gauss complète ou d’une réduction de Gauss sous-diagonale) la
manipulation réversible
On incrémente c et l de 1 et on passe à l’étape 2.
7.1.2 Efficacité de l’algorithme
Si la matrice possède L lignes et C colonnes,
le nombre maximal d’opérations pour réduire une ligne est C divisions,
C multiplications, C soustractions, donc 3C opérations
arithmétiques de base. Il y a L−1 lignes à réduire à chaque
étape et min(L,C) étapes à effectuer, on en déduit que le
nombre maximal d’opérations pour réduire une matrice est
3LCmin(L,C). Pour une matrice carrée de taille n, cela fait
3n3 opérations.
7.1.3 Erreurs d’arrondis du pivot de Gauss
Comme |ajc| ≤ |alc|, une étape de réduction multiplie
au plus l’erreur absolue des coefficients par 2. Donc la
réduction complète d’une matrice peut multiplier au pire l’erreur
absolue sur les coefficients par 2n (où n est le nombre
d’étapes de réduction, inférieur au plus petit du nombre de lignes
et de colonnes). Ceci signifie qu’avec la précision d’un double,
on peut au pire perdre toute précision pour des matrices pas
si grandes que ça (n=52). Heureusement, il semble qu’en pratique,
l’erreur absolue
ne soit que très rarement multipliée par un facteur supérieur à 10.
Par contre, si on ne prend pas la précaution de choisir le pivot
de norme maximale dans la colonne, les erreurs d’arrondis se
comportent de manière bien moins bonnes, cf. l’exemple suivant.
Exemple
Soit à résoudre le système linéaire
є x + 1.0 y = 1.0 , x + 2.0 y = 3.0 |
avec є =2−54 (pour une machine utilisant des doubles pour
les calculs en flottant,
plus générallement on choisira є tel que (1.0+3є)−1.0
soit indistinguable de 0.0).
Si on résoud le système exactement,
on obtient x=1/(1−2є) (environ 1)
et y=(1−3є)/(1−2є) (environ 1).
Supposons que l’on n’utilise pas la stratégie du pivot partiel,
on prend alors comme pivot є, donc on effectue la
manipulation de ligne L2 ← L2 − 1/є L1 ce qui
donne comme 2ème équation (2.0−1.0/є)y=3.0−1.0/є.
Comme les calculs sont numériques, et à cause des erreurs
d’arrondis, cette 2ème équation sera remplacée par
(−1.0/є)y=−1.0/є d’où y=1.0, qui sera remplacé
dans la 1ère équation, donnant є x = 1.0−1.0y=0.0 donc
x=0.0.
Inversement, si on utilise la stratégie du pivot partiel, alors
on doit échanger les 2 équations L2′=L1 et L1′=L2 puis on effectue
L2 ← L2′ − є L1′, ce qui donne
(1.0−2.0є) y = 1.0 − 3.0 є , remplacée en raison
des erreurs d’arrondi par 1.0*y=1.0 donc y=1.0, puis on remplace
y dans L1′ ce qui donne x=3.0−2.0y=1.0.
On observe dans les deux cas que la valeur de y est proche de la
valeur exacte, mais la valeur de x dans le premier cas est
grossièrement eloignée de la valeur correcte.
On peut aussi s’intéresser à la sensibilité de la solution d’un
système linéaire à des variations de son second membre. Le traitement
du sujet à ce niveau est un peu difficile, cela fait intervenir
le nombre de conditionnement de la matrice A du système (qui
est essentiellement la valeur absolue du rapport de la valeur propre
la plus grande sur la valeur propre la plus petite), plus
ce nombre est grand, plus la solution variera (donc plus on
perd en précision).
7.2 Applications de Gauss
7.2.1 Base d’un sous-espace
On réduit la matrice des vecteurs écrits en ligne, puis
on prend les lignes non nulles, dont les vecteurs
forment une base du sous-espace vectoriel engendré par
les vecteurs du départ.
Exemple : base du sous-espac engendré par (1,2,3), (4,5,6), (7,8,9).
On réduit la matrice, la 3ème ligne est nulle
donc on ne garde que les 2 premières lignes (1,0,−1), (0,1,2)
(remarque: une réduction sous-diagonale aurait suffi).
7.2.2 Déterminant
On réduit la matrice (carrée) en notant le nombre
d’inversions de ligne, et on fait le produit des coefficients
diagonaux, on change le signe si le nombre d’inversions de lignes
est impair.
7.2.3 Réduction sous forme échelonnée (rref)
On réduit la matrice puis on divise chaque ligne par son premier
coefficient non nul. Si la matrice représentait un système
linéaire inversible on obtient la matrice identité sur les colonnes
sauf la dernière et la solution en lisant la dernière colonne.
La relation conservée est en effet Ax=b où x est la solution
de l’équation, et à la fin de la réduction A=I.
Par exemple pour résoudre le système
⎧
⎪
⎨
⎪
⎩ | ccc
x + 2y + 3z | = | 5 |
4x + 5y + 6z | = | 0 |
7x + 8y | = | 1 |
|
|
on réduit sous forme échelonnée la matrice
[[1,2,3,5],[4,5,6,0],[7,8,0,1]]
, ce qui donne
[[1,0,0,-9],[0,1,0,8],[0,0,1,-2/3]]
, d’où la solution
x=−9, y=8, z=−2/3.
On accolle la matrice identité à droite de la matrice à inverser. On effectue
la réduction sous forme échelonnée, on doit obtenir à droite
l’identité si la matrice est inversible, on a alors à gauche
la matrice inverse.
La relation conservée est en effet A x=B où x est l’inverse
de la matrice de départ, et en fin de réduction A=I.
Par exemple, pour calculer l’inverse
de [[1,2,3],[4,5,6],[7,8,0]]
, on réduit avec rref
[[1,2,3,1,0,0],[4,5,6,0,1,0],[7,8,0,0,0,1]]
.
On réduit la matrice sous forme échelonnée. Puis on introduit des lignes
de 0 afin que les 1 en tête de ligne se trouvent sur la diagonale
de la matrice. On enlève ou on rajoute des lignes de 0 à la fin
pour obtenir une matrice carrée. Une base du noyau est alors
formée en prenant chaque colonne correspondant à un 0 sur la
diagonale, en remplaçant ce 0 par -1. On vérifie qu’on obtient bien 0
en faisant le produit de ce vecteur par la matrice réduite. De
plus les vecteurs créés sont clairement linéairement indépendants
(puisqu’ils sont échelonnés), et il y en a le bon nombre (théorème
noyau-image).
Exemple : calcul du noyau de [[1,2,3,4],[1,2,7,12]]
, on
réduit la matrice avec rref, ce qui donne
[[1,2,0,-2],[0,0,1,2]]
, on ajoute une ligne de 0 entre
ces 2 lignes pour mettre le 1 de la 2ème ligne sur la diagonale
ce qui donne [[1,2,0,-2],[0,0,0,0],[0,0,1,2]]
, puis
on ajoute une ligne de 0 à la fin pour rendre la matrice carrée.
On obtient ainsi le système équivalent de matrice
[[1,2,0,-2],[0,0,0,0],[0,0,1,2],[0,0,0,0]]
.
La 2ème colonne donne le premier vecteur de la base du noyau,
(2,−1,0,0), la 4ème colonne donne le deuxième vecteur
(−2,0,2,−1), on vérifie aisément que ces 2 vecteurs
forment une famille libre du noyau, donc une base car
la dimension du noyau est égale à 2 (dimension de l’espace de départ
moins le rang de la matrice, c’est-à-dire le nombre de lignes
non nulles de la matrice réduite).
7.2.6 La méthode de factorisation LU
Nous ne la développons pas à ce niveau, elle permet d’écrire
une matrice A comme produit de deux matrices triangulaire
inférieures et supérieures, ce qui permet de ramener la résolution
de système à la résolution de deux systèmes triangulaires.
7.3 Réduction exacte des endomorphismes
On calcule le polynome caractéristique ou le polynome minimal,
on le factorise, et on calcule ensuite le noyau de A−λ I
pour les λ racines. Il existe des méthodes évitant
le calcul de noyau, méthode de Fadeev-Laguerre-Souriau que
nous ne présentons pas ici.
7.3.1 Polynome caractéristique
On peut le calculer en développant le déterminant detA−λ I, mais
il est plus efficace de le calculer par interpolation.
Soit A une matrice carrée de taille n,
on sait que son polynome caractéristique
est un polynome de degré n, il suffit
de connaitre sa valeur en n+1 points distincts, on calcule
donc n+1 déterminants detA−λ I en remplaçant λ par
sa valeur (il y a plus de déterminants à calculer mais
ce sont des déterminants sans paramètre λ donc beaucoup
plus simple à calculer), ce qui permet de reconstruire le polynome
caractéristique par interpolation de Lagrange.
Exercice : pour [[1,-1],[2,4]]
, calculer
det(A−λ I) en λ=0,1,2 puis le polynome d’interpolation,
vérifier que c’est bien le polynome caractéristique.
Il faut effectuer n+1 calculs de déterminants, ce qui nécessite
O(n4) opérations. Il existe des méthodes plus efficaces, par
exemple le calcul du polynome minimal probabiliste présenté
plus bas (O(n3) opérations).
7.3.2 Polynome minimal
Définition 5
Le polynome minimal d’une matrice A est un polynôme M de degré minimal
tel que M(A)=0 et de coefficient dominant égal à 1.
Un tel polynome divise tous les polynomes
tels que P(A)=0, il divise le polynome caractéristique
de A et il a les mêmes racines que le polynome caractéristique.
Preuve:
D’abord M divise tous les polynomes
tels que P(A)=0, car si R désigne le reste de la division
de P par M alors R(A)=(P−QM)(A)=P(A)−Q(A)M(A)=0, donc
R est nul car son degré est plus petit que celui de M.
En particulier le polynome minimal divise le polynome caractéristique C,
car C(A)=0 (on peut montrer que C(A)=0 en faisant
le produit de la matrice A−λ I par sa comatrice,
on obtient le déterminant fois l’identité, soit C(λ)I. Comme
C(λ)I−C(A) peut se factoriser par λ I −A en appliquant
(17) à chaque monome de C, on
en déduit que C(A) se factorise par λ I −A, donc C(A)=0
en regardant les termes de plus haut degré de ces polynomes en
λ à coefficients matriciels).
Montrons enfin que les racines du polynome caractéristique sont racines
du polynome minimal. En effet soit λ une racine
du polynome caractéristique alors A−λ I n’est pas inversible.
Or M(A)−M(λ)I se factorise par A−λ I car
Ak−λk I=(A−λ I) | | λk−1−j Aj
(17) |
donc M(A)−M(λ)I ne peut pas être inversible. Comme M(A)=0
on en déduit que M(λ)I n’est pas inversible donc M(λ)=0,
λ est une racine de M.
Donc si le polynome caractéristique n’a pas de racines multiples,
il est égal au polynome minimal.
Pour calculer M, on peut chercher une relation de degré minimal entre
les puissances de A, en les voyant comme des vecteurs à n2
composantes (ce qui revient à aplatir en un long vecteur
tous les coefficients de la matrice). Cela revient à calculer
le noyau de l’application linéaire dont les colonnes sont les coefficients
des puissances de A (de 0 à n), en gardant le premier vecteur
obtenu par l’algorithme calcul du noyau ci-dessus.
Cette méthode est toutefois couteuse, car il faut réduire une matrice
ayant n2 lignes et n+1 colonnes. Il existe une autre méthode
moins couteuse et qui marche presque toujours. Elle consiste à
calculer le polynome minimal de A par rapport à un vecteur v
c’est-à-dire le polynome de degré minimal (et coefficient dominant
1) tel que Mv(A)v=0. Comme M(A)=0, on a M(A)v=0, donc
Mv divise M qui divise le polynome caractéristique. Si par
chance, on trouve que Mv est de degré n, alors Mv sera
égal à M et au polynome caractéristique. On fait donc le calcul
du noyau de l’application linéaire dont les colonnes sont les Ajv
pour j variant de 0 à n. Si l’on trouve un espace de dimension 1,
alors Mv est de degré n et on a simultanément le polynome
minimal et caractéristique avec le polynome correspondant à
ce vecteur du noyau. Si le degré n’est pas n, on peut essayer
un ou quelques autres vecteurs, et faire le PPCM des polynomes
minimaux obtenus. Si on obtient un polynome de degré n
on conclut, sinon on peut tester si ce polynome évalué en A
est nul, ce sera alors le polynome minimal.
Exemple, on reprend la matrice [[1,-1],[2,4]]
,
et comme vecteur aléatoire v=(1,0), on a Av=(1,−1)
et A(Av)=(−1,−5). On calcule donc le noyau de la matrice
[[1,1,-1],[0,-1,-5]]
(on écrit en colonnes v, Av, A2v),
on trouve que (−6,5,−1) engendre le noyau, donc le polynome
minimal relatif au vecteur v est (au signe près) −6+5x−x2.
Comme il est de degré maximal 2, c’est le polynome minimal
et caractéristique.
7.4 Réduction approchée des endomorphismes
On pourrait trouver des valeurs propres approchées d’une matrice
en calculant le polynome caractéristique ou minimal puis en le
factorisant numériquement. Mais cette méthode n’est pas idéale
relativement aux erreurs d’arrondis (calcul du polynome caractéristiaue,
de ses racines, et nouvelle approximation en calculant le noyau
de A−λ I), lorsqu’on veut calculer quelques valeurs propres
on préfère utiliser des méthodes itératives directement sur A
ce qui évite la propagation des erreurs d’arrondi.
7.4.1 Méthode de la puissance
Elle permet de déterminer la plus grande valeur propre en valeur absolue
d’une matrice diagonalisable lorsque celle-ci est unique.
Supposons en effet que les valeurs propres de A soient
x1,...,xn avec |x1| ≤ |x2| ≤ ... ≤ |xn−1| < |xn|
et soient e1,...,en une base de vecteurs propres correspondants.
On choisit un vecteur aléatoire v et on calcule la suite
vn=Avn−1=An v . Si v a pour coordonnées V1,...,Vn)
dans la base propre, alors
vn = | | Vj xjn ej
= xnn wn, wn=∑Vj | ⎛
⎜
⎜
⎝ | | ⎞
⎟
⎟
⎠ | | ej
|
L’hypothèse que xn est l’unique valeur propre
de module maximal entraine alors
que limn → +∞ wn = Vn en puisque la suite
géométrique de raison xj/xn converge vers 0.
Autrement dit, si Vn≠ 0 (ce qui a une probabilité 1 d’être
vrai pour un vecteur aléatoire),
vn est équivalent à Vn xnn en. Lorsque n
est grand, vn est presque colinéaire au vecteur propre
en (que l’on peut prendre comme vn divisé par sa norme),
ce que l’on détecte en testant si vn+1 et vn
sont presques colinéaires, et de plus le facteur de colinéarité
entre vn+1 et vn est presque xn, la valeur propre
de module maximal.
Exercice : tester la convergence de Anv vers l’espace propre
associé à λ=3 pour la matrice [[1,-1],[2,4]]
et le vecteur (1,0).
Lorsqu’on applique cette méthode a une matrice réelle, il peut arriver
quíl y ait deux valeurs propres conjuguées de module maximal. Le même
type de raisonnement montre que pour n grand,
vn+2 est presque colinéaire à l’espace engendré par vn et vn+1,
la relation vn+2+ x vn+1 + x2 vn=0 permet de calculer les
valeurs propres.
La convergence est de type série géométrique, on gagne le même nombre
de décimales à chaque itération.
7.4.2 Itérations inverses
La méthode précédente permet de calculer la valeur propre
de module maximal d’une matrice. Pour trouver une valeur
propre proche d’une quantité donnée x, on peut appliquer
la méthode précédente à la matrice (A−xI)−1. En effet,
les valeurs propres de cette matrice sont les (xi−x)−1 dont
la norme est maximale lorsqu’on se rapproche de xi.
7.4.3 Elimination des valeurs propres trouvées
Si la matrice A est symétrique, et
si en est un vecteur propre normé écrit en colonne, on peut considérer
la matrice B=A−xn en ent qui possède les mêmes valeurs
propres et mêmes vecteurs propres que A avec même multiplicité,
sauf xn qui est remplacé par 0.
En effet les espaces propres de A sont orthogonaux
entre eux, donc
Ben=xnen −xn en ent en = 0,
Bek = xk ek − xn en ent ek = xk ek |
On peut donc calculer
la 2ème valeur propre (en valeur absolue), l’éliminer
et ainsi de suite.
Si la matrice A n’est pas symétrique, on peut utiliser une technique
analogue lorsque 0 n’est pas valeur propre de A (on peut s’y
ramener en ajoutant à A un multiple de l’identité).
En effet on peut construire un vecteur propre de B pour une
valeur propre xk ≠ 0 à partir d’un vecteur propre
de B, en cherchant y tel que tel que
On obtient pour le membre de gauche :
Bek − yB en =Bek= (A−xn en ent)ek
= xk ek − xn en.ek en |
et pour le membre de droite
d’où l’équation
Néanmoins cette méthode n’est pas stable, en particulier si la valeur
propre ek est proche de 0, car les vecteurs propres se rapprochent
alors tous de en.
8 Quelques références
-
Analyse numérique et équations différentielles, Demailly J.-P.,
Presses Universitaires de Grenoble, 1996
- The Art of Computer Programming, Vol. 2: Seminumerical
algorithms, Knuth D., Addison-Wesley, 1998
- Mathématiques concrètes, illustrées par la TI-92 et la TI-89.
Lemberg H, et Ferrard J.-M., Springer, 1998
- Maths et Maple, J.M. Ferrard, Dunod, 1998
- Handbook of Mathematical Functions, Abramowitz and Stegun,
disponible en ligne sur
http://www.math.sfu.ca/~cbm/aands/toc.htm
- Arithmétique flottante,
Rapport de l’INRIA de V. Lefèvre et P. Zimmermann,
téléchargeable sur
http://www.inria.fr/rrrt/rr-5105.html
- Modern Computer Arithmetic, R. P. Brent, P. Zimmermann,
téléchargeable sur
http://www.loria.fr/~zimmerma/mca/pub226.html
- Matrix computations, Golub and Loa, Hopkins University Press, 1989
- Gantmacher
Index
-
atan, 4.2
- Bezout, 5.1
- bit, 2.3.2
- complexe, 2.4
- contractante, 3.1
- convexe, 3.2
- cos, 4.1
- determinant, 7.2.2
- diagonalisation, 7.3
- division euclidienne, 2.1
- double, 2.3.2
- erreur, 2.3.5, 7.1.3
- exp, 4.1
- exposant, 2.3.2
- expression, 2.4
- factorisation, 5.2.1, 5.2.2, 5.2.3, 5.2.4, 5.2.5
- flottant, 2.3.2
- fonction, 2.4
- Gauss, 7.1
- integration, 6
- interpolation, 5.3
- inverse, 7.2.4
- iterations inverses, 7.4.2
- ker, 7.2.5
- LU, 7.2.6
- lagrange, 5.3, 5.3
- liste, 2.4
- ln, 4.4
- mantisse, 2.3.2
| - matrice, 2.4
- multiplicite, 5.2.1
- Newton, 3.2, 3.2
- Newton-Cotes, 6.4
- noyau, 7.2.5
- ordre, 6.2
- pivot, 7.1
- point fixe, 3.1
- point milieu, 6.1
- polynome, 2.4
- polynome caracteristique, 7.3.1
- polynome minimal, 7.3.2
- puissance, 7.4.1
- quadrature, 6
- racine, 5.2.1, 5.2.3
- racines rationnelles, 5.2.5
- rationnel, 2.2
- rectangle, 6.1
- reduction, 7.2.3
- rref, 7.2.3
- Simpson, 6.3
- Sturm, 5.2.4
- sequence, 2.4
- serie alternee, 4.3
- serie entiere, 4.2
- sin, 4.1
- squarefree, 5.2.1
- symbole, 2.4
- Taylor, 4
- trapeze, 6.1
- vecteur, 2.4
|
A La moyenne arithmético-géométrique.
A.1 Définition et convergence
Soient a et b deux réels positifs,
on définit les 2 suites
u0=a, v0=b, un+1= | | , vn+1= | √ | |
(18) |
On va montrer que ces 2 suites sont adjacentes et convergent donc vers
une limite commune notée M(a,b) et il se trouve que la convergence
est très rapide, en raison de l’identité :
un+1−vn+1= | | ( | √ | | − | √ | | )2
= | | (un−vn)2
(19) |
la convergence est quadratique.
On suppose dans la suite que a≥ b sans changer la généralité puisque échanger a et b
ne change pas la valeur de un et vn pour n>0. On a alors un ≥ vn
(d’après (19) pour n>0) et un+1 ≤ un car
et vn+1=√unvn ≥ √vnvn=vn. Donc (un) est décroissante
minorée (par v0), (vn) est croissante majorée (par u0), ces 2 suites sont
convergentes et comme un+1=un+vn/2, elles convergent vers la même limite
l qui dépend de a et b et que l’on note M(a,b).
On remarque aussi que M(a,b)=bM(a/b,1)=aM(1,b/a).
Précisons maintenant la vitesse de convergence lorsque a ≥ b>0.
On va commencer par estimer le nombre
d’itérations nécessaires pour que un et vn soient du même ordre de grandeur.
Pour cela, on utilise la majoration
ln(un+1)−ln(vn+1) ≤ ln(un)−ln(vn+1) = | |
(ln(un)−ln(vn)) |
donc
ln | | = ln(un)−ln(vn) ≤
| | (ln(a)−ln(b)) = | | ln | | |
Donc si n ≥ ln( ln(a/b)/m)/ln(2) alors
lnun/vn ≤ m (par exemple, on peut prendre m=0.1 pour
avoir un/vn ∈ [1,e0.1]). Le nombre minimum d’itérations n0 est proportionnel
au log du log du rapport a/b.
Ensuite on est ramené à étudier la convergence de la suite arithmético-géométrique
de premiers termes a=un0 et b=vn0 et même en tenant compte
de M(a,b)=aM(1,b/a) à a=1 et b=vn/un donc 0≤ a−b ≤ 1−e−0.1.
Alors l’équation (19) entraine
puis (par récurrence)
Donc comme M(a,b) est compris entre vn et un, l’erreur relative sur la limite
commune est inférieure à une précision donnée є
au bout d’un nombre d’itérations proportionnel au ln(ln(1/є)).
Typiquement dans la suite, on souhaitera calculer M(1,b) avec b de l’ordre
de 2−n en déterminant n chiffres significatifs,
il faudra alors O(ln(n)) itérations pour se ramener à M(1,b) avec b∈ [e−0.1,1]
puis O(ln(n)) itérations pour avoir la limite avec n chiffres significatifs.
Le cas complexe
On suppose maintenant que a, b ∈ ℂ avec ℜ(a)>0, ℜ(b)>0. On va voir que
la suite arithmético-géométrique converge encore.
Étude de l’argument
On voit aisément (par récurrence)
que ℜ(un)>0 ; de plus ℜ(vn) > 0 car par définition de la racine carrée
ℜ(vn)≥ 0 et est de plus non nul car le produit de deux complexes d’arguments dans
]−π/2,π/2[ ne peut pas être un réel négatif.
On en déduit que
arg(un+1)=arg(un+vn) se trouve dans l’intervalle de bornes
arg(un) et arg(vn) et que arg(vn+1)=1/2(arg(un)+arg(vn))
donc
| arg(un+1−arg(vn+1) | ≤ | | |arg(un)−arg(vn)| |
Après n itérations, on a
Après quelques itérations, un et vn seront donc presque alignés.
Faisons 4 itérations.
On peut factoriser par exemple vn et on
est ramené à l’étude de la suite de termes initiaux a=un/vn d’argument
arg(un)−arg(vn) petit
(inférieur en valeur absolue à π/16) et b=1. On suppose donc dans la suite que
Étude du module
On a :
| = | | ⎛
⎜
⎜
⎜
⎜
⎜
⎜
⎜
⎝ | | | + | | ⎞
⎟
⎟
⎟
⎟
⎟
⎟
⎟
⎠ |
Posons un/vn=ρn eiθn, on a :
| = | | ⎪
⎪
⎪
⎪
⎪
⎪ | √ | | eiθn/2
+ | | e−iθn/2 | ⎪
⎪
⎪
⎪
⎪
⎪ |
|
| = | | | ⎪
⎪
⎪
⎪
⎪
⎪ | ( | √ | | + | | )cos | |
+ i ( | √ | | − | | )sin | | | ⎪
⎪
⎪
⎪
⎪
⎪ |
|
| = | |
| = | |
|
Si ρ désigne le max de ρn et 1/ρn, on a alors la majoration
donc en prenant les logarithmes
lnρn+1 ≤ | | lnρ= | | |lnρn|
(20) |
On rappelle qu’on a la majoration
qui va nous donner la minoration de ρn+1
en prenant les log et en minorant ln(1−x) par −2x
lnρn+1 ≥ | | (−|lnρn|+ln(1 − | | ))
≥ − | | (|lnρn|+ | | ) |
Finalement avec (20)
On en déduit
|lnρn| ≤ | | lnρ0 + | | + ... +
| | + | |
= | | lnρ0 + | | |
La convergence du ln(un/vn) vers 0 est donc géométrique, donc un et vn convergent
quadratiquement.
A.2 Lien avec les intégrales elliptiques
Le calcul de la limite commune des suites un et vn en fonction
de a et b n’est pas trivial
au premier abord. Il est relié aux intégrales elliptiques, plus
précisément on peut construire une intégrale dépendant
de deux paramètres a et b et qui est invariante par
la transformation un,vn → un+1,vn+1 (18)
On a en effet
On pose alors
où t → u est une bijection croissante de t∈]0,+∞[ vers
u ∈ ]−∞,+∞[, donc
| = | ∫ | | | dt/2(1+ab/t2) |
|
| |
(( | | )2+1/4(t−ab/t)2)(ab+1/4(t−ab/t)2) |
|
|
|
|
|
| = | |
|
On note au passage que I est définie si a,b ∈ ℂ vérifient ℜ(a)>0, ℜ(b)>0,
on peut montrer que la relation ci-dessus s’étend (par holomorphie).
Lorsque a=b=l (par exemple lorsqu’on est à la limite),
le calcul de I(l,l) est explicite
donc
On peut transformer I(a,b) en posant t=bu
Puis en posant u=tan(x) (du=(1+u2) dx)
I(a,b)= | | | ∫ | |
| | |
1+tan(x)2 |
|
1+(b/a)2tan(x)2 |
|
| dx |
et enfin en posant tan2(x)=sin(x)2/1−sin(x)2
Si on définit pour m<1
alors on peut calculer K en fonction de I, en posant
m=1−b2/a2 soit b2/a2=1−m
d’où l’on déduit la valeur de l’intégrale elliptique en fonction
de la moyenne arithmético-géométrique :
Dans l’autre sens, pour x et y positifs
et finalement
A.3 Application : calcul efficace du logarithme.
On peut utiliser la moyenne arithmético-géométrique pour
calculer le logarithme efficacement, pour cela on cherche le développement
asymptotique de K(m) lorsque m tend vers 1. Plus précisément,
on va poser 1−m=k2 avec k ∈ ]0,1], donc
en posant y=π/2−x, et
la singularité de l’intégrale pour k proche
de 0 apparait lorsque y est proche de 0.
Si on effectue un développement de Taylor en y=0, on trouve
sin(y)2+k2 cos(y)2 = k2 + (1−k2) y2 + O(y4) |
Il est donc naturel de comparer K(m) à l’intégrale
qui se calcule en faisant par exemple le changement de variables
ou directement avec Xcas,
supposons(k>0 && k<1);
J:=int(1/sqrt(k^2+(1-k^2)*y^2),y,0,pi/2)
qui donne après réécriture :
J= | |
| ⎛
⎜
⎜
⎜
⎜
⎝ | ln | ⎛
⎜
⎜
⎝ | | ⎞
⎟
⎟
⎠ | +
ln | ⎛
⎜
⎜
⎜
⎜
⎝ | | | ⎛
⎜
⎜
⎜
⎜
⎝ | | |
+ | √ | | | ⎞
⎟
⎟
⎟
⎟
⎠ | ⎞
⎟
⎟
⎟
⎟
⎠ | ⎞
⎟
⎟
⎟
⎟
⎠ |
(22) |
et on peut calculer le développement asymptotique de J en 0
series(J,k=0,5,1)
qui renvoie :
J =ln | ⎛
⎜
⎜
⎝ | | ⎞
⎟
⎟
⎠ | +O( | ⎛
⎜
⎜
⎝ | | ⎞
⎟
⎟
⎠ | | ) |
on peut alors préciser ce développement par
series(J+ln(k)-ln(pi),k=0,5,1)
qui renvoie (après simplifications et où la notation Õ peut contenir des logarithmes)
⎛
⎜
⎜
⎝ | | + | | ⎞
⎟
⎟
⎠ | k2 +
Õ(k4) |
donc
J=−ln(k)+ln(π)+ | ⎛
⎜
⎜
⎝ | | +
| | ⎞
⎟
⎟
⎠ | k2 + Õ(k4)
(23) |
Examinons maintenant K−J, il n’a plus de singularité en y=0, et il admet une limite
lorsque k → 0, obtenue en remplacant k par 0
(K−J)|k=0 = | ∫ | | | ⎛
⎜
⎜
⎝ | | − | | ⎞
⎟
⎟
⎠ |
dy = | ⎡
⎢
⎢
⎣ | ln | ⎛
⎜
⎜
⎝ | tan | ⎛
⎜
⎜
⎝ | | ⎞
⎟
⎟
⎠ | ⎞
⎟
⎟
⎠ | − ln(y) | ⎤
⎥
⎥
⎦ | | =
ln( | | ) |
D’où pour K
Kk → 0 = ln | ⎛
⎜
⎜
⎝ | | ⎞
⎟
⎟
⎠ | + O( | ⎛
⎜
⎜
⎝ | | ⎞
⎟
⎟
⎠ | | ) |
Pour préciser la partie du développement de K en puissances de k, nous allons
majorer K−J−ln(4/π), puis J−ln(π/k).
Posons
A=sin(y)2+k2 cos(y)2, B=y2+(1−y2)k2 |
Majoration de K−J−ln(4/π)
L’intégrand de la différence K−J−ln(4/π) est
Soit
K−J−ln( | | )= | ∫ | |
| (y−sin(y))[(1−k2)y sin(y)(y+sin(y))− | √ | | ( | √ | | + | √ | | )] |
|
|
|
|
(27) |
On décompose l’intégrale en 2 parties [0,k] et [k,π/2].
Sur [0,k] on utilise (25), on majore chaque terme séparément
et on minore A et B par
A=k2+(1−k2)sin(y)2 ≥ k2, B=k2+(1−k2)y2 ≥ k2 |
Donc
| ≤ | |
| ≤ | |
| ≤ | |
+ ln(sin( | | )) −ln( | | ) − ln(cos( | | )) |
|
| ≤ | |
| ≤ | |
| ≤ | |
|
Sur [k,π/2], on utilise (27)
et on minore A et B par
A=sin(y)2+k2 cos(y)2 ≥ sin(y)2, B=y2+(1−y2)k2 ≥ y2 |
on obtient
| | ∫ | | | ≤ | ∫ | |
| (y−sin(y))|C| |
|
y sin(y) (y+sin(y)) |
| ,
|
où :
C | = | (1−k2)y sin(y)(y+sin(y))−A | √ | | +B | √ | | |
|
| = | −A( | √ | | −y)−B( | √ | | −sin(y))
−Ay−Bsin(y) + (1−k2)y sin(y)(y+sin(y)) |
|
| = | −A( | √ | | −y)−B( | √ | | −sin(y)) − k2(y+sin(y)) |
|
|
Donc
|C| | ≤ | A( | √ | | −y)+B( | √ | | −sin(y)) + k2(y+sin(y)) |
|
| ≤ | |
| ≤ | |
|
et
| | ∫ | | | ≤
| ∫ | |
| (y−sin(y))k2( | | + | | + (y+sin(y))) |
|
|
y sin(y) (y+sin(y)) |
|
On peut majorer y−sin(y) ≤ y3/6, donc
On majore enfin A et B par 1,
Le premier morceau se calcule par intégration par parties
| = | | | ⎛
⎜
⎜
⎜
⎜
⎝ | [− | | ]kπ/2
+ | ∫ | | | |
| ⎞
⎟
⎟
⎟
⎟
⎠ |
|
| = | | | ⎛
⎜
⎜
⎝ | | + [ln(sin(y))] | | | ⎞
⎟
⎟
⎠ |
|
| = | | | ⎛
⎜
⎜
⎝ | | −ln(sin(k)) | ⎞
⎟
⎟
⎠ |
|
| ≤ | |
|
Le deuxième morceau se majore en minorant sin(y)≥ (2y)/π
Finalement
|K−J−ln( | | )|
≤ k2 | ⎛
⎜
⎜
⎝ | − | | ln(k) + | | + | | + | | +
| | | ⎞
⎟
⎟
⎠ |
où J est donné en (22).
Majoration de J−ln(π/k)
On a
|J − ln | ⎛
⎜
⎜
⎝ | | ⎞
⎟
⎟
⎠ | |
= | ⎪
⎪
⎪
⎪
⎪
⎪ | ( | | −1) ln | ⎛
⎜
⎜
⎝ | | ⎞
⎟
⎟
⎠ |
+ | |
ln | ⎛
⎜
⎜
⎜
⎜
⎝ | | | ⎛
⎜
⎜
⎜
⎜
⎝ | | |
+ | √ | | | ⎞
⎟
⎟
⎟
⎟
⎠ | ⎞
⎟
⎟
⎟
⎟
⎠ | ⎪
⎪
⎪
⎪
⎪
⎪ |
et on va majorer la valeur absolue de chaque terme de la somme.
Pour k≤ 1/2, on a
Pour le second terme, on majore le facteur 1/√1−k2 par 2/√3,
l’argument du logarithme est inférieur à 1 et supérieur à
| (1 − | | +1− | | )
= 1 − k2 ( 1− | | ) > 1−k2
|
donc le logarithme en valeur absolue est inférieur à
donc, pour k≤ 1/2,
|J−ln | ⎛
⎜
⎜
⎝ | | ⎞
⎟
⎟
⎠ | | ≤
| | ln | ⎛
⎜
⎜
⎝ | | ⎞
⎟
⎟
⎠ |
+ k2 | |
|
Finalement, pour k<1/2
|K−ln | ⎛
⎜
⎜
⎝ | | ⎞
⎟
⎟
⎠ | |
≤ k2 | ⎛
⎜
⎜
⎜
⎜
⎝ | | + | |
+ | | + | |
− ( | | + | | ) ln(k) | ⎞
⎟
⎟
⎟
⎟
⎠ |
(28) |
que l’on peut réécrire
| | | −ln | ⎛
⎜
⎜
⎝ | | ⎞
⎟
⎟
⎠ | |
≤ k2(3.8−0.8ln(k))
(29) |
La formule (29)
permet de calculer le logarithme d’un réel positif
avec (presque) n bits
lorsque k ≤ 2−n/2 (ce à quoi on peut toujours se ramener
en calculant le logarithme d’une puissance 2m-ième de x ou
le logarithme de 2mx, en calculant au préalable ln(2)).
Par exemple, prenons k=2−27, on trouve (en 8 itérations)
M(1,2−27)=M1=0.0781441403763.
On a, avec une erreur inférieure à 19 × 2−54=1.1× 10−15
On peut donc déduire une valeur approchée de π si on connait
la valeur approchée de ln(2) et réciproquement.
Si on veut calculer les deux simultanément, comme les relations entre ln
et π seront des équations homogènes, on est obligé
d’introduire une autre relation. Par exemple pour calculer une
valeur approchée de π on calcule la différence
ln(229+1)−ln(229) dont on connait le développement au premier
ordre, et on applique la formule de la moyenne arithmético-géométrique.
Il faut faire attention à la perte de précision lorsqu’on fait
la différence des deux logarithmes qui sont très proches, ainsi
on va perdre une trentaine de bits, il faut grosso modo calculer les
moyennes arithmético-géométrique avec
2 fois plus de chiffres significatifs.
L’intérêt de cet algorithme apparait lorsqu’on veut calculer
le logarithme avec beaucoup de précision, en raison de la
convergence quadratique de la moyenne arithmético-géométrique
(qui est nettement meilleure que la convergence linéaire
pour les développements en série, ou logarithmiquement
meilleure pour l’exponentielle), par contre elle n’est pas
performante si on ne veut qu’une dizaine de chiffres significatifs.
On peut alors calculer les autres
fonctions transcendantes usuelles, telle l’exponentielle,
à partir du logarithme, ou les fonctions trigonométriques
inverses (en utilisant des complexes) et directes.
On trouvera dans Brent-Zimmermann quelques considérations permettant
d’améliorer les constantes dans les temps de calcul par rapport
à cette méthode (cela nécessite d’introduire des fonctions
spéciales θ) et d’autres formules pour calculer π.
Ce document a été traduit de LATEX par HEVEA