Mat 249

Table des matières

• 1  Présentation du module
• 2  Représentation des nombres et autres données, calcul exact/approché
  • 2.1  Entiers courts et longs
  • 2.2  Les rationnels.
  • 2.3  Les réels
    • 2.3.1  Virgule fixe et flottante.
    • 2.3.2  Les flottants au format double
    • 2.3.3  Opérations sur les flottants
    • 2.3.4  Erreurs
    • 2.3.5  Erreur absolue, relative et propagation des erreurs.
  • 2.4  Types composés.
• 3  Suites itératives et applications
  • 3.1  Le point fixe
  • 3.2  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 ℂ.
    • 5.2.3  Calcul approché des racines complexes simples
    • 5.2.4  Factorisation dans ℝ, localisation des racines
    • 5.2.5  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  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 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 :

1. calcul exact et approché, représentation des données
2. Suites récurrentes, méthode du point fixe, de Newton
3. Interpolation polynômiale (évaluation, interpolation de Lagrange)
4. Séries de Taylor et approximation des fonctions usuelles
5. Arithmétique des polynômes (PGCD, Bézout, factorisation, décomposition en éléments simples)
6. Intégration
7. 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

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 d_n...d₀, on traduit les divisions euclidiennes successives en

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=2⁴+2³+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 [−2³¹,2³¹−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, 2³¹+2³¹ 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 2³²). 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 à N² (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.

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⁻⁴. 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,bⁿ[. Soit un réel x représenté par

Si m∈ [0,bⁿ⁻¹[, alors on peut aussi écrire x=m′ b^e−1 avec m′=mb ∈ [0,bⁿ[, 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 ≥ 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 ∈ [bⁿ⁻¹,bⁿ[, on a alors une erreur d’arrondi relative majorée par

On appelle flottant normalisé un flottant tel que m ∈ [bⁿ⁻¹,bⁿ[. 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⁻⁵³ (b=2 et n=52+1 pour les doubles).

Exemples :

• en base 10 avec n=6, pour représenter π=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, (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ⁿ⁻¹[. 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,2⁵²[ sur 52 bits, et l’exposant e ∈ [0, 2¹¹[ sur 11 bits. Pour les nombres “normaux”, l’exposant est en fait compris entre 1 et 2¹¹−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 2¹ pour être entre 1 et 2 dont e+1−2¹⁰=1, l’exposant est e=2¹⁰. 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¹⁰=0 soit e=2¹⁰−1=2⁹+2⁸+2⁷+2⁶+2⁵+2⁴+2³+2²+2¹+2⁰. On a r−1=1/2=2⁻¹. D’où la représentation
  0 01111111111 10000000...0000
• 6.4=32/5
  Positif. Il faut le diviser par 2² pour avoir 8/5 ∈ [1,2[ donc e+1−2¹⁰=2 soit e=2¹⁰+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=2¹¹−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

2.3.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 (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 −x₀| ≤ ε₀ et si |y−y₀| ≤ ε₁ alors par l’inégalité triangulaire (|a+b|≤ |a|+|b|), on a :

on dit que les erreurs absolues s’additionnent.

Mais comme il faut représenter x₀+y₀ en machine, on doit ajouter une erreur d’arrondi, qui est proportionnelle à la valeur absolue de x₀+y₀ 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⁻⁵³, à 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₀,y₀ sont non nuls, on a

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⁻⁵³ sur x₀ y₀.

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⁻⁵³, alors l’erreur relative sur xy sera de

Lorsque l’erreur relative sur les données est grande devant 2⁻⁵³, 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₀=1+2⁻⁵² avec une erreur d’arrondi de 2⁻⁵³ et y par y₀=−1 avec la même erreur d’arrondi, l’addition de x et y renvoie 2⁻⁵² avec une erreur absolue de 2 * 2⁻⁵³ (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

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(x₁,...,x_n) vu comme polynôme en x_n à coefficients polynômes dépendant des variables x₁,...,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²−5x+y³) 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 u_n 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.

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

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∈[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

Donc par une récurrence évidente :

ce qui entraine d’ailleurs que |u_n−l| ≤ kⁿ |a−b|. Comme k ∈ [0,1[ , la suite géométrique kⁿ 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′)| ≤ 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 :

puis on majore avec l’inégalité triangulaire

puis on applique le fait que f est contractante de rapport k

soit

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)² 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)²=−1/4 et −1/(2+1)²=−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₀=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

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⁻⁶ près, il suffira que la différence entre deux termes successifs de la suite u_n 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

Si u₀ ∈ I, alors θ ∈ I donc |u₁−l| ≤ |u₀−l| et u₁ ∈ I, par récurrence on a pour tout n, u_n ∈ I

on a donc par récurrence

Donc pour avoir |u_n−l| ≤ є il suffit que

et il faut que

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 x₀ est proche de la racine r on peut faire un développement de Taylor à l’ordre 1 de la fonction f en x₀ :

Pour trouver une valeur approchée de r, on ne garde que la partie linéaire du développement, on résout :

donc (si f′(x₀)≠ 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 u₀ proche de la racine et par la relation :

Il y a deux théorèmes importants, l’un d’eux prouve que si u₀ 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₀ 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

En appliquant un développement de Taylor de f en u_n à l’ordre 2, on obtient pour un réel θ situé entre r et u_n :

donc :

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 u_n est dans cet intervalle, alors θ aussi donc

On a |u_n−r| ≤ ε, on diminue si nécessaire ε pour avoir ε < 2/(Mm), on a alors :

donc d’une part u_n+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 |u_n−r|², plus précisément, on montre par récurrence que

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 ℝⁿ.

Exemple : pour calculer √2, on écrit l’équation x²−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 u₀ ∈ [√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)=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).

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² est convexe :

Démonstration :
L’équation de la tangente au graphe en x₀ est

Soit

on a :

donc g′ est croissante, comme g′(x₀)=0, g′ est négative pour x<x₀ et positive pour x>x₀, donc g est décroissante pour x<x₀ et croissante pour x>x₀. On conclut alors que g ≥ 0 puisque g(x₀)=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′′ ≥ 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 ≤ u_n. Comme la courbe représentative de f est au-dessus de la tangente, on a u_n+1 ≥ 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≥ r. À la limite, on a

donc l=r car f>0 sur ]r,b].

Comme (u_n) est décroissante, on a bien 0 ≤ u_n −r, pour montrer l’autre inégalité, on applique le théorème des accroissements finis, il existe θ ∈ [r,u_n] 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 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 ℝ^+∗ (car k ≥ 2). Il suffit donc de prendre une valeur de départ u₀ 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 u_n²−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 x₀ ∈ I. On peut alors effectuer le développement de Taylor de f en x₀ à l’ordre 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 : il existe θ compris entre x₀ et x tel que

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, la dérivée n-ième de f est exp(x), donc R_n(x)=exp(θ)xⁿ⁺¹/(n+1)! avec θ compris entre 0 et x, ainsi si x est positif |R_n(x)| ≤ e^x xⁿ⁺¹/(n+1)! et si x est négatif, |R_n(x)| ≤ xⁿ⁺¹/(n+1)!. Dans les deux cas, la limite de R_n 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 e^x? Il suffira d’arrêter la sommation lorsque R:=xⁿ⁺¹/(n+1)! si x<0 ou lorsque R:=e^x xⁿ⁺¹/(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 e^x à une erreur relative ε donnée (par exemple ε=2⁻⁵³ pour stocker le résultat dans un double) il suffit que R/e^x < ε, donc si x est positif, il suffit que xⁿ⁺¹/(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⁻¹⁰ 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⁻⁵³ ? 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/e^x (ne change pas l’erreur relative). Pour les grands réels, on peut utiliser e^2x=(e^x)² (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 :

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 e^r, que l’on peut minimiser en utilisant la méthode de Horner pour évaluer T_n(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 e^x 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⁻⁵³).

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

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^∞a_n xⁿ lorsque a_n 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 x₀ tel que |a_n x₀ⁿ| est borné (ce sera le cas en particulier si la série converge en x₀), alors

la série converge donc en x si |x|<|x₀| 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).

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 ∑xⁿ (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 a_n 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).

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+x²), de même pour arccos(x) et arcsin(x) en intégrant celui de (1−x²)^−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

Démonstration :
on montre que les suites v_n=S_2n et w_n=S_2n+1 sont adjacentes. On a

donc v_n est décroissante, de même w_n est croissante, et v_n−w_n=u_2n+1 est positif et tend vers 0. On en déduit que v_n et w_n convergent vers la même limite S telle que v_n>S>w_n 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 u_n tende assez vite vers 0, sinon il y aura perte de précision sur la mantisse lorsqu’on effectuera u_2n−u_2n+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 x^k+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

  ln(1+x)= 2 ln(

  √

  1+x

  )

  on observe que :

  X=

  √

  1+x

  −1 =

  x 1+ √ 1+x

  ≤

  x 2

  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 +

  −x 1+x

  )

  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 y₀ 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 e^y=1+x, on pose f(y)=e^y−(1+x), on utilise la suite itérative

  yn+1 = yn −

  eyn−(1+x) eyn

  Comme y₀ 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 e^y 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) 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

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/(n2ⁿ) 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

  ∫

  +∞ x

  e−t2  dt

  ainsi que

  ∫

  +∞ 0

  e−t2  dt =

  √ π  2

  Le développement asymptotique s’obtient par exemple en changeant de variable u=t² 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 a_n+1 et a_n. 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 |a_n+1/a_n|<C, alors le reste de la série entière est majoré par |a_nxⁿ|/(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

dont on cherche une solution série entière y=∑_k=0^∞a_k x^k . En remplacant dans l’équation, si x est dans le rayon de convergence de la série (rayon supposé non nul), on obtient

soit encore

Par exemple, prenons le cas m=0. On a alors a₀ quelconque, a₁ nul et pour k≥ 2

Donc tous les a d’indice impair sont nuls. Les pairs sont non nuls si a₀≠ 0, et ils sont de signe alterné. Soit x fixé, on observe que pour 2k > |x|,

donc la série ∑_k=0^∞a_k x^k 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 a₀=1 et |x|≤ 1, on calcule y=∑_k=0²ⁿ a_k x^k , on s’arrête au rang n tel que

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 y₀ 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=y₀ 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

où

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)!/xⁿ⁺¹<є où є est la précision relative que l’on s’est fixée. Par exemple, si x≥ 100, n=11 convient pour є=12!/100¹²=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 u_n+1−u_n qui vaut −1/2n²) 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¹⁴=4e−12

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 R₅₀ est positif et est inférieur à (-1)^50*13.^51/51/51!) qui est de l’ordre de 8e-12

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

ce terme vaut 13^11/11/11! soit 4000 environ)

On obtient finalement comme valeur approchée de γ

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

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 t²/2, etc.

où

Pour déterminer I_n on fait le changement de variables t=nu

Or en faisant le même changement de variables t=nu :

Donc

Lorsque n tend vers l’infini, on peut montrer que I_n → 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 limI_n=−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

En faisant tendre α vers −i, −ln(α) tend vers ln(i)=iπ/2 et on obtient

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

Algorithme :
On construit en fait 3 suites (U_n), (V_n) et (R_n) telles que :

• on initialise U₀=1, V₀=0, R₀=A et U₁=0, V₁=1, R₁=B
• on calcule les indices n+2 en fonction de n et n+1 en effectuant la division euclidienne de R_n par R_n+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=x³−1, B=x²+1, les rangs 0 et 1 sont donnés ci-dessus. Au rang 2, Q₀ est le quotient euclidien de A par B (fonction quo) donc x, d’où

Puis on divise x²+1 par −x−1, quotient −x+1, donc

Preuve de l’algorithme :
On montre facilement par récurrence que la relation AU_n + B V_n = R_n est conservée. Comme R_n 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 V_k. Supposons que deg(A) ≥ deg(B) (sinon on échange A et B). Au rang n=0, V₀=0 donc V₂=−Q₀ V₁, aux rangs suivants le degré de Q_n est non nul (car le degré de R_n+1 est strictement inférieur au degré de R_n) On montre donc par récurrence que la suite des degrés de V_n est croissante et que :

Comme deg(Q_n)=deg(R_n)-deg(R_n+1), on en déduit que

Donc si n+2 est le rang du dernier reste non nul, V_n+2=V et degV=degA-degR_n+1 est donc strictement inférieur au degré de A (car R_n+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

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−x² par x² ce qui donne une solution

l’ensemble des solutions est de la forme

et la solution priviligiée (de degrés minimaux) est

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 :