Mat 406

N.B. : la version HTML de ce document est interactive, vous pouvez exécuter certains champs d’entrée donnés en exemple en les modifiant éventuellement, il suffit de cliquer sur le bouton exe. Il est recommandé d’utiliser Firefox ou un navigateur compatible supportant MathML.

Table des matières

• 1  Présentation du module
• 2  Représentation des nombres et autres données, calcul exact/approché
  • 2.1  Représentation des entiers
  • 2.2  Les réels
    • 2.2.1  Virgule fixe et flottante.
    • 2.2.2  Les flottants au format double
    • 2.2.3  Opérations sur les flottants
    • 2.2.4  Erreurs
    • 2.2.5  Erreur absolue, relative, arrondi propagation des erreurs.
  • 2.3  L’arithmétique d’intervalle.
  • 2.4  Types composés.
    • 2.4.1  Les rationnels.
    • 2.4.2  Les complexes
    • 2.4.3  Les polynômes
    • 2.4.4  Calcul symbolique
    • 2.4.5  Listes, séquences, tables
• 3  Suites itératives et applications
  • 3.1  Rappel : suite récurrente
  • 3.2  Le point fixe
  • 3.3  La méthode de Newton.
• 4  Développement de Taylor, séries entières, fonctions usuelles
  • 4.1  La fonction exponentielle
  • 4.2  Séries entières.
  • 4.3  Série alternée
  • 4.4  La fonction logarithme
  • 4.5  Autres applications
    • 4.5.1  Exemple : la fonction d’erreur (error fonction, erf)
    • 4.5.2  Recherche de solutions d’équations différentielles
    • 4.5.3  Exemple : fonctions de Bessel d’ordre entier
  • 4.6  Développements asymptotiques et séries divergentes
• 5  Polynômes : arithmétique, factorisation, interpolation
  • 5.1  Arithmétique des polynomes: Bézout et applications
  • 5.2  Factorisation des polynômes
    • 5.2.1  Multiplicité des racines.
    • 5.2.2  Factorisation dans $\mathbb{C}$.
    • 5.2.3  Calcul approché des racines complexes simples
    • 5.2.4  Localisation d’une racine complexe près d’une racine approchée
    • 5.2.5  Localisation des racines réelles : Sturm
    • 5.2.6  Localisation des racines réelles : règle des signes de Descartes
    • 5.2.7  Factorisation exacte
  • 5.3  Approximation polynomiale
• 6  Intégration numérique
  • 6.1  Les rectangles et les trapèzes
  • 6.2  Ordre d’une méthode
  • 6.3  Simpson
  • 6.4  Newton-Cotes
  • 6.5  En résumé
• 7  Algèbre linéaire
  • 7.1  Le pivot de Gauss
    • 7.1.1  L’algorithme
    • 7.1.2  Efficacité de l’algorithme
    • 7.1.3  Erreurs d’arrondis du pivot de Gauss
  • 7.2  Applications de Gauss
    • 7.2.1  Base d’un sous-espace
    • 7.2.2  Déterminant
    • 7.2.3  Réduction sous forme échelonnée (rref)
    • 7.2.4  Inverse
    • 7.2.5  Noyau
    • 7.2.6  La méthode de factorisation $LU$
  • 7.3  Réduction exacte des endomorphismes
    • 7.3.1  Polynome caractéristique
    • 7.3.2  Polynome minimal
  • 7.4  Réduction approchée des endomorphismes
    • 7.4.1  Méthode de la puissance
    • 7.4.2  Itérations inverses
    • 7.4.3  Elimination des valeurs propres trouvées
• 8  Guide rapide KhiCAS sur calculatrices
  • 8.1  Installation
    • 8.1.1  Installation Casio Graph 90
    • 8.1.2  Installation TI Nspire CX/CX II
    • 8.1.3  Installation Numworks N0110 non verrouillée.
  • 8.2  Le shell KhiCAS
  • 8.3  Programmation
  • 8.4  Compatibilité avec Xcas
• 9  Quelques références
• A  La moyenne arithmético-géométrique.
  • A.1  Définition et convergence
  • A.2  Lien avec les intégrales elliptiques
  • A.3  Application : calcul efficace du logarithme.

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 19.5h de TP en 13 séances). 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 :

1. Calcul exact et approché, représentation des données
2. Suites récurrentes, méthode du point fixe, de Newton.
3. Séries de Taylor et approximation des fonctions usuelles
4. Arithmétique des polynômes (PGCD, Bézout, factorisation, décomposition en éléments simples)
5. Interpolation polynomiale
6. Intégration approchée.
7. Algèbre linéaire.

L’ordre des deux premiers thèmes sera inversé en cours.

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 sont autorisées au DS et à l’examen final (prêt possible 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 $\mathbb{N}$ et relatifs $\mathbb{Z}$, les rationnels $\mathbb{Q}$, les réels $\mathbb{R}$ et les complexes $\mathbb{C}$. Sur ordinateur, on peut représenter ces nombres de manière exacte dans certains cas, approchée dans d’autres.

2.1  Représentation des entiers

Preuve : On prend pour $q$ le plus grand entier tel que $a-bq \geq 0$.
Exemple : iquorem(23,7) exe

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 $d_n...d_0$, on traduit les divisions euclidiennes successives en $$\begin{matrix} d &=&( ... ((d_n b +d_{n-1})b + d_{n-2})...+d_1)b+d_0\\ &=& d_n b^n + d_{n-1} b^{n-1} + ... + d_0 \end{matrix}$$ Par exemple, vingt-cinq s’écrit en base 16 0x19 car 25 divisé par 16 donne quotient 1, reste 9
convert(25,base,16) exe

En base 2, on trouverait 0b11001 car $25=2^4+2^3+1$.
convert(25,base,2) exe

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 $[-2^{31},2^{31}-1]$.

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, $2^{30}+2^{30}+2^{30}+2^{30}$ renverra 0. Ils sont utilisables avec tous les langages de programmation traditionnels.

Les logiciels de calcul formel et certains logiciels de programmation permettent de travailler avec des entiers de taille beaucoup plus grande, ainsi qu’avec des rationnels, permettant de faire du calcul exact, mais on paie cette exactitude par un temps de calcul plus long, de plus pas mal de méthodes numériques ne gagnent rien à faire des calculs intermédiaires exacts. Néanmoins, l’utilisation d’un logiciel de calcul formel permettra dans certains cas d’illustrer certains phénomènes dus au calcul approché.

2.2  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.2.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 \in [0,b^n[$. Soit un réel $x$ représenté par $$x=mb^e, \quad m \in [0,b^n[$$ Si $m\in [0,b^{n-1}[$, alors on peut aussi écrire $x=m' b^{e-1}$ avec $m'=mb \in [0,b^n[$, 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 $mb^e$, on doit l’arrondir, par exemple au plus proche réel de la forme $mb^e$. La distance de $x$ à ce réel est inférieure ou égale à la moitié de la distance entre deux flottants consécutifs, $mb^e$ et $(m+1)b^e$, donc l’erreur d’arrondi est inférieure ou égale à $b^e/2$. Si on divise par $x \geq mb^e$, 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 \in [b^{n-1},b^n[$, on a alors une erreur d’arrondi relative majorée par $$\frac{1}{2b^{n-1}}$$

On appelle flottant normalisé un flottant tel que $m \in [b^{n-1},b^n[$. 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 :

• en base 10 avec $n=6$, pour représenter $\pi=3,14159265...$, on doit décaler la virgule de 5 positions, on obtient 314159.265... on arrondit à $314159$ donc on obtient 314159e-5.
• en base 2 avec $n=10$, pour représenter trois cinquièmes ($3/5$ en base 10, noté $11/101$ en base 2), on pose la division en base 2 de 11 par 101, ce qui donne

   11        | 101
   110       ---------
  -101       | 0.1001
  ----       |
    010      |
     100     |
     1000    |
    - 101    |
    -----    |
      011    |
  

  on retrouve le nombre de départ donc le développement est périodique et vaut 0.1001 1001 1001 .... On décale le point de 10 positions, on arrondit, donc trois cinquièmes est représenté par la mantisse 1001100110 et l’exposant -10. On observe aussi sur cet exemple que $3/5$ dont l’écriture en base 10 0.6 est exacte, n’a pas d’écriture exacte en base 2 (de même que 1/3 n’a pas d’écriture exacte en base 10).

Il existe une exception à la possibilité de normaliser les flottants, lorsqu’on atteint la limite inférieure de l’exposant $e$. Soit en effet $e_m$ le plus petit exposant des flottants normalisés et considérons les flottants $x=b^{e_m}(1+1/b)$ et $y=b^{e_m}$. 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, (puisque 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,b^{n-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éralement 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.2.2  Les flottants au format double

Cette section développe les notions de la section précédente pour les flottants machine selon la norme IEEE-754, 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=\pm 1$ sur 1 bit, la mantisse $M \in [0,2^{52}[$ sur 52 bits, et l’exposant $e \in [0, 2^{11}[$ sur 11 bits. Pour les nombres “normaux”, l’exposant est en fait compris entre 1 et $2^{11}-2$, le nombre représenté est le rationnel $$(1+\frac{M}{2^{52}}) 2^{e+1-2^{10}}$$ 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 \in [0,1[$. Exemples :

• -2
  Signe négatif. Il faut diviser sa valeur absolue 2 par $2^1$ pour être entre 1 et 2 dont $e+1-2^{10}=1$, l’exposant est $e=2^{10}$. 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-2^{10}=0$ soit $e=2^{10}-1=2^9+2^8+2^7+2^6+2^5+2^4+2^3+2^2+2^1+2^0$. 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 $2^{2}$ pour avoir $8/5 \in [1,2[$ donc $e+1-2^{10}=2$ soit $e=2^{10}+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.5 - 5*0.1 n’est pas nul.

Des représentations spéciales (avec $e=0$ ou $e=2^{11}-1$) ont été introduites pour représenter $\pm \infty$ (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

Voici un programme C++ affichant la représentation interne des flottants

2.2.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 arrondit 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

2.2.4  Erreurs

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. 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 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 et une plage d’exposants plus grande. 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. Dans ce cas, il est pertinent lorsqu’on évalue $f(x)$ avec $x$ mal connu de calculer aussi $f'(x)$, en effet : $$f(x(1+h))= f(x)+xh f'(x) + O(h^2)$$ l’erreur relative sur $f(x)$ est donc au premier ordre multipliée par $$|\frac{xf'(x)}{f(x)}|$$ Par exemple, l’erreur relative sur $e^x$ est au premier ordre l’erreur relative sur $x$ multipliée par $|x|$.

2.2.5  Erreur absolue, relative, arrondi 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 -x_0| \leq \varepsilon_0$ et si $|y-y_0| \leq \varepsilon_1$ alors par l’inégalité triangulaire ($|a+b|\leq |a|+|b|$), on a : $$|(x+y)-(x_0+y_0)| \leq |x-x_0| + | y-y_0 | \leq \varepsilon_0 + \varepsilon_1$$ on dit que les erreurs absolues s’additionnent.

Mais comme il faut représenter $x_0+y_0$ en machine, on doit ajouter une erreur d’arrondi, qui est proportionnelle à la valeur absolue de $x_0+y_0$ d’où la notion d’erreur relative :

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 $x_0,y_0$ sont non nuls, on a $$\left| \frac{xy-x_0 y_0}{x_0 y_0} \right| = \left| \frac{x}{x_0} \frac{y}{y_0} -1 \right| = \left| (\frac{x}{x_0}-1)(\frac{y}{y_0} -1)+(\frac{x}{x_0}-1)+(\frac{y}{y_0} -1) \right|$$ 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 $x_0 y_0$.

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} \approx 3 \times 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 $x_0=1+2^{-52}$ avec une erreur d’arrondi de $2^{-53}$ et $y$ par $y_0=-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 \rightarrow 1+2^{-52}$$ alors que $$2.0^{-53}+(2.0^{-53}+1.0) \rightarrow 1$$ Dans Xcas, il n’y a que 48 bits de mantisse :

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$ car le terme d’erreur est en $O(h^2)$ et non en $O(h)$. Attention toutefois à ne pas prendre $h$ trop petit, sinon $x+h=x$ en flottants et même si $x+h \neq x$, l’erreur absolue sur $f(x+h)-f(x-h)$ est (au moins) d’ordre $\varepsilon |f(x)|$, donc l’erreur relative est d’ordre $\varepsilon/h |f(x)|$. Par exemple pour h=1e-8 le reste est en $O(h^2)$ donc de l’ordre des erreurs d’arrondi mais l’erreur relative sur $f(x+h)-f(x-h)$ est d’ordre $\epsilon/h$ largement supérieure (en flottants double-précision). On choisira plutôt $h$ tel que $\epsilon/h$ soit proche de $h^2$, donc de l’ordre de 1e-5, qui fournira une valeur approchée avec une erreur relative de l’ordre de 1e-10. Exemple : calcul de la dérivée numérique de $\exp(\sin(x))$ en $x=1$

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, par exemple si on effectue $n$ calculs susceptibles de provoquer des erreurs indépendantes suivant une même loi d’espérance nulle, la moyenne des erreurs divisée par l’écart-type de la loi tend vers une loi normale centrée réduite. De manière plus déterministe, on a l’inégalité de Bienaymé-Tchebyshev $$P(|X|>\alpha) \leq \frac{n\sigma^2}{\alpha^2}$$ où $X$ est la variable aléatoire somme des $n$ erreurs, $\alpha$ l’erreur et $n\sigma^2$ la variance de la somme $n$ erreurs supposées indépendantes, cette probabilité tend vers 0 pour $n$ grand si $\alpha$ est d’ordre $n$, et ne tend pas vers 0 si $\alpha$ est de l’ordre de $\sqrt{n}$. Exemple : somme de $n=400$ nombres répartis sur $[-1,1]$ selon la loi uniforme (représentant des erreurs), on divise par $\sqrt{n}=$20, on effectue plusieurs tirages (par exemple 500) on trace l’histogramme et on compare avec la loi normale de moyenne nulle (l’espérance de la somme) et d’écart-type celui de la loi uniforme. m:=ranm(400,500,-1..1):;gl_x=-2..2;histogram(sum(m)/20/stddev(uniform,-1,1),-1,0.1); plot(normald(0,1),-2..2) exe

Attention, si on effectue la somme de $n$ réels $\sum_j x_j$, les erreurs d’arrondis ne satisfont pas à ces hypothèses. En effet, l’erreur d’arrondi à chaque opération est une erreur relative, l’erreur absolue correspondante est $\epsilon |x_1+x_2|$ puis $\epsilon |x_1+x_2+x_3|$ puis ... $\epsilon |x_1+x_2+...+x_n|$, que l’on peut majorer par $$\epsilon ((n-1)|x_1|+(n-2)|x_2|+...+|x_n||)$$ La majoration de l’erreur d’arrondi dépend donc de l’ordre des termes, on a intérêt à sommer en commençant par les termes les plus petits en valeur absolue. Voici des programmes C++ illustrant cela

• uniquement avec des flottants. À compiler avec la commande c++ erreur2.cc -o erreur2, ne pas utiliser d’option d’optimisation sinon les flottants intermédiaires seront probablement calculé en précision plus grande (64 bits de mantisse)
• en comparant avec des flottants multi-précision (nécessite d’avoir Giac installé sous Unix/Linux, se compile avec c++ erreur.cc -lgiac -lgmp -o erreur).

Mais on peut faire mieux, il est possible de corriger les erreurs d’arrondi dans une somme avec le programme suivant pour une liste (on peut bien sur adapter à la somme d’une expression dépendant d’une variable entière sans stocker de liste) :

Somme(l):={
  local x,s,c;
  s:=0.0;
  c:=0.0;
  pour x in l faire
    c += (x-((s+x)-s));
    s += x;
  fpour;
  print(c);
  return s+c;
}:;

exeonload
En effet, $c$ devrait valoir 0 sans erreurs d’arrondis, avec les erreurs d’arrondis, on a le premier calcul $s+x$ qui donnera une erreur opposée à celui du calcul de $s$ à la ligne suivante, le 2ième calcul effectué $(s+x)-s$ donne une erreur absolue en $\epsilon |x|$ au pire (car c’est une erreur relative par rapport à $(s+x)-s$), et la 3ième erreur d’arrondi est négligeable (puisque la somme vaut 0). On a donc une erreur absolue sur $s+c$ qui est au premier ordre au pire en $O(\epsilon \sum|x_i|)$, bien meilleure que la majoration $\epsilon ((n-1)|x_1|+(n-2)|x_2|+...+|x_n||)$ calculée précédemment.

Par exemple Digits:=14; n:=25000; l:=seq(1/j,j,1,n):; Somme(l); exe
à comparer avec S:=sum(1/j,j,1,n):; evalf(S); sum(1. /j,j,1,n); exe
(le calcul de $S$ est fait en exact, celui de sum(1. /j,j,1,n) est approché sans correction).

En conclusion, il est souvent très difficile de calculer une majoration rigoureuse de l’erreur pour des calculs (sauf les plus simples), et cette majoration est en général bien trop pessimiste. 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, ou 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 (section suivante)

2.3  L’arithmétique d’intervalle.

Certains systèmes de calcul formel peuvent manipuler directement des intervalles réels, par exemple par l’intermédiaire de la bibliothèque C MPFI. Les opérations arithmétiques sur des intervalles renvoient alors le meilleur intervalle possible contenant toutes les valeurs possibles lorsque les opérandes parcourent leurs intervalles respectifs. Exemple en Xcas (version 1.1.1 et ultérieures) : [-1..2]*[-1..2] renvoie [-2..4]. Attention ici on parcourt toutes les valeurs possibles de $xy, \ x \in [-1,2], y \in [-1,2]$. Ce qui est différent du carré d’un intervalle ou plus généralement de l’évaluation d’un polynôme en un intervalle, horner(x^2,[-1..2]) renvoie ainsi [0..4].

Les fonctions disponibles sont souvent moins riches qu’en arithmétique flottante, le calcul d’une fonction non monotone sur un intervalle peut s’avérer délicat, alors que si la fonction est monotone, il suffit de calculer l’image des deux bornes de l’intervalle. Pour les polynômes, Xcas décompose les coefficients en deux parties $P=P_+-P_-$ en fonction du signe, puis utilise la monotonie de $P_+$ et $P_-$ sur $\mathbb{R}^+$ et $\mathbb{R}^-$ respectivement.

L’arithmétique d’intervalle dans $\mathbb{C}$ est beaucoup plus difficile à mettre en oeuvre puisqu’il n’y a plus d’ordre ni de monotonie, on doit alors s’en remettre à des estimations sur les parties réelles et imaginaires qui ne tiendront pas compte du phénomène ci-dessus sur la différence entre $xy, \ x \in [-1,2], y \in [-1,2]$ et $x^2, \ x \in [-1,2]$.

2.4  Types composés.

2.4.1  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.

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 souvent gérés par le microprocesseur (ou son coprocesseur arithmétique).

2.4.2  Les complexes

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.

2.4.3  Les polynômes

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(x_1,...,x_n)$ vu comme polynôme en $x_n$ à coefficients polynômes dépendant des variables $x_1,...,x_{n-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 y^2-5x+y^3$) 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).

2.4.4  Calcul symbolique

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.

2.4.5  Listes, séquences, tables

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  Rappel : suite récurrente

Une suite récurrente à un cran est définie par une valeur initiale $u_0$ et une relation de récurrence que l’on va supposer ne pas dépendre explicitement de $n$, donc de la forme $$u_{n+1}=f(u_n)$$ On peut représenter graphiquement une suite de ce type par un graphe “en toile d’araignée”, par exemple pour les cinq premiers termes de $f(x)=\frac{x+2}{x+1}$ et $u_0=3$,

gl_x=-0.5..2.5;plotseq((x+2)/(x+1),[0,0,2],5)

exeonload
La suite peut converger ou pas, par exemple avec $f(x)=\tan(x)-0.1$ et $f(x)=\cos(x)$ plotseq(tan(x)-0.1,[0,-1.3,1.3],7); exe
plotseq(cos(x),[0,-2,2],5); exe

3.2  Le point fixe

Soit $f$ une fonction continue sur un intervalle $I=[a,b]$ de $\mathbb{R}$, 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 $$u_{n+1}=f(u_n), \quad u_0 \in I \qquad (1)$$ Supposons que $u_n$ converge vers une limite $l \in I$ lorsque $n \rightarrow +\infty$, alors la limite doit vérifier $$f(l)=l$$ 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$.

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'(\theta) (y-x), \quad \theta \in [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

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\in[a,b]$ tel que $g(l)=0$ (théorème des valeurs intermédiaires). Soit $u_n$ une suite définie par (1). On a alors pour tout $n$ $$|u_{n+1}-l|=|f(u_n)-f(l)| \leq k |u_n-l|$$ Donc par une récurrence évidente : $$|u_n-l| \leq k^n |u_0-l|$$ ce qui entraine d’ailleurs que $|u_n-l| \leq k^n |a-b|$. Comme $k \in [0,1[$, la suite géométrique $k^n$ converge vers 0 lorsque $n$ tend vers l’infini, donc $u_n$ tend vers $l$. Notons que $l$ est unique car si $l'$ est une autre solution alors $|l-l'|=|f(l)-f(l')| \leq k|l-l'|$ donc $(1-k)|l-l'| \leq 0$, or $1-k>0$ et $|l-l'| \geq 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$ : $$u_n-l= u_n - u_{n+1} + u_{n+1} - u_{n+2} + ... + u_{n+m-1}- u_{n+m} + u_m-l$$ puis on majore avec l’inégalité triangulaire $$|u_n-l| \leq \sum_{j=0}^{m-1} |u_{n+j}-u_{n+j+1}| + |u_m-l|$$ puis on applique le fait que $f$ est contractante de rapport $k$ $$|u_n-l| \leq \sum_{j=0}^{m-1} k^j |u_{n}-u_{n+1}| + |u_m-l|$$ soit $$|u_n-l| \leq \frac{1-k^m}{1-k} |u_{n}-u_{n+1}| + |u_m-l|$$ On fait alors tendre $m$ vers l’infini d’où le résultat.

Exemple 1 :
Cherchons une valeur approchée de $\sqrt{2}$ par cette méthode. Il faut d’abord trouver une fonction $f$ dont $\sqrt{2}$ est un point fixe, par exemple $$f(x)=\frac{x+2}{x+1}$$ On vérifie que $f(\sqrt{2})=\sqrt{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 $u_0=1$ et on va converger vers $\sqrt{2}$ (en s’en rapprochant à chaque cran d’un rapport inférieur à $1/4$).
f:=(x+2)/(x+1); f1:=factor(f'); x0:=fMin(f1,x=1..2); f1(x=x0); exe

u:=1.0; pour j de 1 jusque 10 faire u:=f(x=u); print(j,u); fpour; u; exe

Exemple 2 :
Considérons l’équation en $x$ $$x- e \sin(x) =t, \quad e \in [0,1[$$ c’est l’équation du temps utilisée en astronomie pour trouver la position (plus précisément l’angle avec le grand axe de l’ellipse) d’une planète à l’instant $t$ sur son orbite elliptique ($e$ étant l’excentricité de l’ellipse).

gl_x=-3.5..3.5;gl_y:=-2.2..2.2; 
E:=ellipse(-1,1,2,color=red); cercle(0,2); theta:=2; 
line(0,slope=tan(theta)); angle(0,5,5*exp(i*theta),"x");
D:=droite(x=2*cos(theta)); M:=inter(D,E,i,legend="M(t)");
E:=ellipse(-1,1,2,color=red); cercle(0,2); theta:=2; 
line(0,slope=tan(theta)); angle(0,5,5*exp(i*theta),"x");
D:=droite(x=2*cos(theta)); M:=inter(D,E,i,legend="M(t)",color=red)

exeonload
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 $$f(x)=t+e\sin(x) = x$$ On observe que $f$ envoie $\mathbb{R}$ dans $[t-e,t+e]$ donc on peut prendre $I=[t-e,t+e]$, de plus $|f'|\leq e <1$, $f$ est contractante de rapport $e \in [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 $u_n$ vérifie $$|u_{n+1}-u_n| \leq 10^{-6} (1-e)$$ on aura alors bien : $$|u_n-x| \leq \frac{|u_{n+1}-u_n|}{1-e} \leq 10^{-6}$$ Par exemple, pour $t=\pi/2$ et $E=0.0167$ (la Terre), au bout de 3 itérations
E:=0.0167; t:=pi/2.; f:=t+E*sin(x); u0:=t; u1:=f(x=u0); u2:=f(x=u1); u3:=f(x=u2); abs(u3-u2) exe

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-\eta,l+\eta]$ tel que $$x \in I \Rightarrow \frac{L}{2} \leq |f'(x)| \leq \frac{1+L}{2}$$ Le théorème des accroissements finis donne alors $$|u_{n+1} - l | = |f(u_n)-f(l)| = |f'(\theta)| |u_n-l|, \quad \theta \in [u_n,l]$$ Si $u_0 \in I$, alors $\theta \in I$ donc $|u_1-l| \leq |u_0-l|$ et $u_1 \in I$, par récurrence on a pour tout $n$, $u_n \in I$ $$\frac{L}{2} |u_n-l| \leq |u_{n+1} - l| \leq \frac{1+L}{2} |u_n-l|$$ on a donc par récurrence $$\left(\frac{L}{2}\right)^n|u_0-l| \leq |u_n-l| \leq \left( \frac{1+L}{2} \right)^n|u_0-l|$$ Donc pour avoir $|u_n-l| \leq \epsilon$ il suffit que $$\left( \frac{1+L}{2} \right)^n|u_0-l| \leq \epsilon \Rightarrow n \geq \frac{\ln(\frac{\epsilon}{|u_0-l|})}{\ln( \frac{1+L}{2}) }$$ et il faut que $$\left(\frac{L}{2}\right)^n |u_0-l| \leq \epsilon \Rightarrow n \geq \frac{\ln(\frac{\epsilon}{|u_0-l|})}{\ln( \frac{L}{2}) }$$

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$).

def dicho(f,a,b,eps):
    # local c,niter;
    if f(a)*f(b)>=0: 
        return "erreur: f(a)*f(b)>=0"
    while b-a>eps:
        c=(a+b)/2.0
        if f(a)*f(c)>0:
            a=c
        else:
            b=c
    return c

exeonload
f(x):=cos(x)-x; dicho(f,0.0,1.0,1e-8); exe

Cette méthode a toutefois l’avantage de se généraliser en dimension supérieure à 1, contrairement à la dichotomie.

3.3  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 $x_0$ est proche de la racine $r$ on peut faire un développement de Taylor à l’ordre 1 de la fonction $f$ en $x_0$ : $$f(x)=f(x_0)+(x-x_0)f'(x_0)+O((x-x_0)^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 \approx f(x_0) + (r-x_0) f'(x_0)$$ donc (si $f'(x_0)\neq 0$) : $$r \approx x_0 -\frac{f(x_0)}{f'(x_0)}$$ Graphiquement, cela revient à tracer la tangente à la courbe représentative de $f$ et à chercher où elle coupe l’axe des $x$.

gl_x=-3..3;
f:=x^2-2; 
G:=plot(f,x=-3..3);
x0:=1.9; M:=element(G,x0);
T:=tangent(G,x0,legend="",color=red);
N:=inter_unique(T,droite(y=0));
legend(-3,"xN="+round(abscisse(N),3))

exeonload
On considère donc la suite récurrente définie par une valeur $u_0$ proche de la racine et par la relation : $$u_{n+1} = u_n -\frac{f(u_n)}{f'(u_n)}$$ Remarque : C’est un cas particulier de point fixe pour l’équation $g(x)=x$ avec $$g(x)=x-\frac{f(x)}{f'(x)}, \quad g'(x)=-\frac{f(x) f'{'}(x)}{f'(x)^2}$$ et donc $g'$ s’annule en un point où $f$ est nulle. On espère ainsi avoir une convergence très rapide près d’une racine.

Il y a deux théorèmes importants, l’un d’eux prouve que si $u_0$ est “assez proche” de $r$ alors la suite $u_n$ converge vers $r$, malheureusement il est difficile de savoir en pratique si on est “assez proche” de $u_0$ 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.

Démonstration : on a $$u_{n+1}-r = u_n - r - \frac{f(u_n)}{f'(u_n)} = \frac{(u_n-r)f'(u_n)-f(u_n)}{f'(u_n)}$$ En appliquant un développement de Taylor de $f$ en $u_n$ à l’ordre 2, on obtient pour un réel $\theta$ situé entre $r$ et $u_n$ : $$0 = f(r)=f(u_n)+(r-u_n) f'(u_n) + (r-u_n)^2 \frac{f'{'}(\theta)}{2}$$ donc : $$(u_n-r)f'(u_n)-f(u_n)= (u_n-r)^2 \frac{f'{'}(\theta)}{2}$$ d’où : $$|u_{n+1}-r| \leq |u_n-r|^2 \frac{1}{|f'(u_n)|} \frac{|f'{'}(\theta)|}{2}$$ On commence par choisir un intervalle $[r-\varepsilon,r+\varepsilon]$ contenant strictement $r$ et tel que $|f'{'}|<M$ et $|1/f'|<m$ sur $[r-\varepsilon,r+\varepsilon]$ (c’est toujours possible car $f'{'}$ et $1/f'$ sont continues au voisinage de $r$ puisque $f'(r)\neq 0$). Si $u_n$ est dans cet intervalle, alors $\theta$ aussi donc $$|u_{n+1}-r| \leq |u_n-r|^2 \frac{Mm}{2} \leq \frac{|u_n-r|Mm}{2} |u_n-r|,$$ On a $|u_n-r| \leq \varepsilon$, on diminue si nécessaire $\varepsilon$ pour avoir $\varepsilon < 2/(Mm)$, on a alors : $$|u_{n+1}-r| \leq k |u_n-r|, \quad k=\frac{\varepsilon Mm}{2}<1$$ donc d’une part $u_{n+1}$ est encore dans l’intervalle $[r-\varepsilon,r+\varepsilon]$ 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 $|u_n-r|^2$, plus précisément, on montre par récurrence que $$|u_n-r| \leq |u_0 - r|^{2^n} \left( \frac{Mm}{2} \right)^{2^n-1}$$ 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 $\mathbb{C}$ et même sur $\mathbb{R}^n$.

Exemple : pour calculer $\sqrt{2}$, on écrit l’équation $x^2-2=0$ qui a $\sqrt{2}$ comme racine simple sur $I=[1/2,2]$, on obtient la suite récurrente $$u_{n+1} = u_n - \frac{u_n^2-2}{2u_n}$$ Si on prend $\eta=1/2$, on a $f'=2x$ et $f'{'}=2$ donc on peut prendre $M=2$ et $m=1$ car $|1/f'|\leq 1$ sur $[\sqrt{2}-1/2,\sqrt{2}+1/2]$. On a $2/(mM)=1$, on peut donc prendre $\varepsilon=1/2$, la suite convergera pour tout $u_0 \in [\sqrt{2}-1/2,\sqrt{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)=x^k-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). Session Xcas

Passons maintenant à un critère très utile en pratique :

Il existe un critère simple permettant de savoir si une fonction de classe $C^2$ est convexe :

Démonstration :
L’équation de la tangente au graphe en $x_0$ est $$y=f(x_0)+f'(x_0)(x-x_0)$$ Soit $$g(x)=f(x)-(f(x_0)+f'(x_0)(x-x_0))$$ on a : $$g(x_0)=0, \quad g'(x)=f'(x)-f'(x_0), \quad g'(x_0)=0, \quad g'{'}=f'{'} \geq 0$$ donc $g'$ est croissante, comme $g'(x_0)=0$, $g'$ est négative pour $x<x_0$ et positive pour $x>x_0$, donc $g$ est décroissante pour $x<x_0$ et croissante pour $x>x_0$. On conclut alors que $g \geq 0$ puisque $g(x_0)=0$. Donc $f$ est bien au-dessus de sa tangente.

On arrive au deuxième théorème sur la méthode de Newton

Démonstration :
On a $f'{'} \geq 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 $u_{n+1} \leq u_n$. Comme la courbe représentative de $f$ est au-dessus de la tangente, on a $u_{n+1} \geq r$ (car $u_{n+1}$ est l’abscisse du point d’intersection de la tangente avec l’axe des $x$). La suite $u_n$ est donc décroissante minorée par $r$, donc convergente vers une limite $l\geq r$. À la limite, on a $$l=l-\frac{f(l)}{f'(l)} \Rightarrow f(l)=0$$ donc $l=r$ car $f>0$ sur $]r,b]$.

Comme $(u_n)$ est décroissante, on a bien $0 \leq u_n -r$, pour montrer l’autre inégalité, on applique le théorème des accroissements finis, il existe $\theta \in [r,u_n]$ tel que $$f(u_n)-f(r)=(u_n-r)f'(\theta)$$ comme $f(r)=0$, on a $$u_n-r = \frac{f(u_n)}{f'(\theta)}$$ 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'{'} \geq 0$ sur $[a,r]$. Si $f'{'} \leq 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 $x^k-a$ est 2 fois continument dérivable et de dérivée première $kx^{k-1}$ et seconde $k(k-1)x^{k-2}$ strictement positives sur $\R^{+\ast}$ (car $k \geq 2$). Il suffit donc de prendre une valeur de départ $u_0$ plus grande que la racine $k$-ième, par exemple $1+a/k$ (en effet $(1+a/k)^k \geq 1+k a/k=1+a$). En appliquant l’inégalité du théorème, on a : $$0 \leq u_n - \sqrt[k]{a} \leq \frac{u_n^k - a}{k\sqrt[k]{a}^{k-1} } \leq \frac{u_n^k-a}{ka} \sqrt[k]{a} \leq \frac{u_n^k-a}{ka} (1+\frac{a}{k})$$ Pour avoir une valeur approchée de $\sqrt[k]{a}$ à $\varepsilon$ près, on peut donc choisir comme test d’arrêt $$u_n^k -a \leq \frac{ka}{1+\frac{a}{k}} \varepsilon$$ Par exemple pour $\sqrt{2}$, le test d’arrêt serait $u_n^2-2 \leq 2 \varepsilon$.

Remarque : Lorsque $u_n$ est proche d’une racine, le calcul approché de $f(u_n)$ peut poser des problèmes à cause de compensations entre deux valeurs proches que l’on soustrait. Par exemple pour $f(x)=x^2-2$, si $x$ est proche de $\sqrt{2}$, $x^2$ sera proche de 2, et le calcul de la différence va perdre énormément de précision relative.

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 $\mathbb{R}$ et $x_0 \in I$. On peut alors effectuer le développement de Taylor de $f$ en $x_0$ à l’ordre $n$ $$T_n(f)(x)= f(x_0) + (x-x_0) f'(x_0) + ... + (x-x_0)^n \frac{f^{[n]}(x_0)}{n!}$$ et se demander si $T_n(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 $T_n(f)(x)$. Si c’est le cas, on pourra utiliser $T_n(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 

La preuve se fait en considérant la fonction $$g(t) := T_{n,x}(f)(t)+C(x-t)^{n+1}-f(t), \quad T_{n,x}(f)(t)= f(x) + (t-x) f'(x) + ... + (t-x)^n \frac{f^{[n]}(x)}{n!}$$ qui vérifie $g(x)=0$ ainsi que toutes ses dérivées jusqu’à l’ordre $n$. On choisit $C$ pour que $g(x_0)=0$ et on déduit que la dérivée s’annule en $\theta_0$ entre $x_0$ et $x$, puis la dérivée seconde s’annule en $\theta_1$ entre $\theta_0$ et $x$ donc entre $x_0$ et $x$, etc. ce qui donne la valeur de $C$ pour la dérivée $n+1$-ième de $g$.

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 $x_0=0$, la dérivée $n$-ième de $f$ est $\exp(x)$, donc $R_n(x)=\exp(\theta)x^{n+1}/(n+1)!$ avec $\theta$ compris entre 0 et $x$, ainsi si $x$ est positif $|R_n(x)| \leq e^x x^{n+1}/(n+1)!$ et si $x$ est négatif, $|R_n(x)| \leq x^{n+1}/(n+1)!$. Dans les deux cas, la limite de $R_n$ est 0 lorsque $n$ tend vers l’infini, car pour $n \geq 2x$, on a $$\frac{x^{n+1}}{(n+1)!} = \frac{x^n}{n!} \frac{x}{n+1}\leq \frac{1}{2}\frac{x^n}{n!}$$ on a donc pour tout $x$ réel $$e^x = \lim_{n \rightarrow +\infty} T_n(f)(x) = \lim_{n \rightarrow +\infty} \sum_{k=0}^n \frac{x^k}{k!} = \sum_{k=0}^\infty \frac{x^k}{k!}$$