Interpolation polynomiale

Objectifs du chapitre

À la fin de ce chapitre, vous serez capable de:

énoncer le problème d'interpolation, le distinguer de l'approximation au sens des moindres carrés, et démontrer qu'il possède une solution polynomiale de degré $\leq n$ et une seule à travers $n + 1$ nœuds distincts;
expliquer pourquoi la matrice de Vandermonde résout le problème en théorie et le ruine en pratique, en estimant la croissance de son conditionnement $κ_{\infty}$ ;
construire le polynôme d'interpolation sous la forme de Lagrange et sous la forme de Newton, et dire laquelle choisir selon que les nœuds sont figés ou que l'on en ajoute;
remplir un tableau de différences divisées, l'utiliser pour évaluer $p_{n}$ par le schéma de Horner, et relier $f [x_{0}, \dots, x_{k}]$ à une dérivée $k$ -ième;
démontrer la formule de l'erreur $f (x) - p_{n} (x) = \frac{f ^{(n + 1)} ( ξ )}{( n + 1 )!} \prod_{i} (x - x_{i})$ et lire ce que chacun de ses trois facteurs apporte;
reconnaître le phénomène de Runge, savoir qu'augmenter le degré peut dégrader l'approximation, et choisir entre les nœuds de Tchebychev et l'interpolation par morceaux pour y échapper.

Le problème d'interpolation

Les cinq premiers chapitres ont traité des problèmes où la fonction était donnée par une formule que l'on pouvait évaluer autant de fois que nécessaire. Celui-ci part de la situation inverse: on ne dispose que d'un tableau de valeurs. Un capteur a enregistré une température toutes les dix minutes; une table thermodynamique donne l'enthalpie de la vapeur d'eau tous les dix degrés; un calcul coûteux de mécanique des fluides a livré une traînée pour six angles d'incidence. Et l'on a besoin d'une valeur entre deux mesures, ou d'une dérivée, ou d'une intégrale.

Deux remarques sur cette définition, qui déterminent tout le chapitre.

Le nombre de conditions égale le nombre de degrés de liberté. Un polynôme de $P_{n}$ s'écrit $a_{0} + a_{1} x + \dots + a_{n} x^{n}$ : il a coefficients, et (6.1) impose conditions. C'est un système carré, et la section suivante montre qu'il est toujours inversible. Avec moins de nœuds que de coefficients, la solution ne serait pas unique; avec davantage, elle n'existerait en général pas.

Interpoler n'est pas approcher. L'interpolation exige le passage exact par chaque point. C'est le bon choix quand les $y_{i}$ sont fiables — une table de valeurs mathématiques, le résultat d'un calcul déterministe. Ce n'est pas le bon choix quand les $y_{i}$ sont des mesures bruitées: y passer exactement, c'est reproduire le bruit, et le polynôme s'en trouve plus faux entre les points qu'une droite qui ne passerait par aucun. La démarche inverse — se donner une famille plus petite que le nombre de données et minimiser un écart plutôt que l'annuler — est l'approximation au sens des moindres carrés; elle fait, avec les splines, l'objet du chapitre 7, et le présent chapitre ne l'aborde pas. Retenez la frontière: si vos données sont exactes, interpolez; si elles sont bruitées, approchez.

Question 6.1

Vous disposez de 40 mesures de température, entachées d'une incertitude de $\pm 0, 3$ °C. Un collègue propose de les interpoler par un polynôme de degré 39. Quelle objection est la plus pertinente?

Existence et unicité

Démonstration. Écrivons le polynôme cherché dans la base canonique, $p_{n} (x) = a_{0} + a_{1} x + \dots + a_{n} x^{n}$ . Les conditions (6.2) forment le système linéaire , où

V = 11 ⋮ 1 x_{0} x_{1} ⋮ x_{n} \dots \dots \dots x_{0}^{n} x_{1}^{n} ⋮ x_{n}^{n} (6.3)

est la matrice de Vandermonde des nœuds. Son déterminant vaut

det V = 0 \leq i < j \leq n \prod (x_{j} - x_{i}), (6.4)

produit non nul dès que les nœuds sont deux à deux distincts. La matrice est donc inversible, et le système admet une solution et une seule: existence et unicité sont acquises d'un coup.

Cette première démonstration a l'inconvénient d'invoquer (6.4), que nous n'avons pas établie. En voici une seconde, indépendante, qui n'utilise que le théorème fondamental de l'algèbre et qui donne l'unicité en trois lignes.

Unicité. Supposons que $p$ et $q$ appartiennent tous deux à $P_{n}$ et vérifient (6.2). Posons $d = p - q$ . Alors $d \in P_{n}$ , donc est de degré au plus , et

d (x_{i}) = p (x_{i}) - q (x_{i}) = y_{i} - y_{i} = 0, i = 0, \dots, n .

Le polynôme $d$ possède donc $n + 1$ racines distinctes. Or un polynôme non nul de degré au plus $n$ en a au plus $n$ . Donc $d$ est le polynôme nul, c'est-à-dire $p = q$ .

Existence. L'unicité que l'on vient d'obtenir dit que l'application linéaire $Φ : P_{n} \to R^{n + 1}$ , $p \mapsto (p (x_{0}), \dots, p (x_{n}))$ , est . Comme , le théorème du rang (Algèbre linéaire) la rend : tout vecteur est atteint, donc le polynôme cherché existe. La section suivante en donnera de surcroît une formule explicite, ce qui fournit une troisième preuve d'existence, entièrement constructive.

Chiffrons ce dernier point, car il est le plus important et le moins souvent montré. Prenons les nœuds équidistants $x_{i} = i / n$ de $[0, 1]$ , construisons la matrice (6.3) et calculons $κ_{\infty} (V) = ∥ V ∥_{\infty} ∥ V^{- 1} ∥_{\infty}$ en arithmétique (fractions rationnelles), pour éviter que la mesure du conditionnement soit elle-même victime du conditionnement.

$n$	$κ_{\infty} (V)$	$κ_{\infty} \cdot ε_{mach}$	erreur observée sur $a$

La dernière colonne est une mesure, pas une prédiction: on a pris le polynôme $p (x) = 1 + x + x^{2} + \dots + x^{n}$ , dont tous les coefficients valent exactement $1$ , on a formé $y_{i} = p (x_{i})$ , puis résolu par l'élimination de Gauss avec pivot partiel en double précision, et l'on reporte . À il ne reste que trois chiffres corrects sur seize; à , le coefficient le plus faux vaut au lieu de , c'est-à-dire qu'il . La majoration du chapitre 4, erreur relative , prédit correctement l'ordre de grandeur — elle est pessimiste d'un facteur 10 à 40, ce qui est le comportement habituel d'une borne du pire cas.

Le choix de l'intervalle n'y change rien d'essentiel. Sur $[- 1, 1]$ avec $x_{i} = - 1 + 2 i / n$ , le conditionnement est bien plus doux mais croît de la même façon: il vaut exactement $6$ pour $n = 2$ , $18$ pour $n = 3$ , pour , puis , et atteint à — un facteur d'environ par degré supplémentaire. Aucune normalisation ne rend la base canonique utilisable en grand degré: ses éléments deviennent presque colinéaires sur un intervalle borné, et c'est cela, et non la faute de l'élimination de Gauss, que le conditionnement mesure.

Question 6.2

Sur les nœuds équidistants de $[0, 1]$ avec $n = 15$ , le conditionnement de la matrice de Vandermonde vaut $1, 25 \cdot 1 0^{13}$ . Combien de chiffres significatifs corrects, au plus, pouvez-vous espérer sur les coefficients calculés en double précision? Donnez un entier.

La forme de Lagrange

La première écriture utile renverse le point de vue: au lieu de chercher les coefficients dans la base canonique, on fabrique une base de $P_{n}$ adaptée aux nœuds, dans laquelle les coefficients sont les données elles-mêmes.

Deux propriétés, et elles suffisent à tout.

Première propriété — le symbole de Kronecker. Pour tous $i, k$ ,

ℓ_{i} (x_{k}) = δ_{ik} = {10 si i = k, si i \neq = k . (6.6)

La vérification est immédiate: si $k \neq = i$ , le facteur d'indice $j = k$ du numérateur vaut $x_{k} - x_{k} = 0$ et annule tout le produit; si , chaque facteur vaut .

Seconde propriété — la partition de l'unité. Pour tout $x$ ,

i = 0 \sum n ℓ_{i} (x) = 1. (6.7)

La démonstration est un joli usage de l'unicité: la somme de gauche est un polynôme de degré au plus $n$ qui, d'après (6.6), vaut $1$ en chacun des $n + 1$ nœuds; le polynôme constant $1$ aussi; le théorème 6.1 les identifie. Autrement dit, interpoler la fonction constante $1$ redonne la constante $1$ — une propriété que l'on aimerait trouver évidente et qu'il vaut mieux avoir démontrée, car elle sert de test de recette à toute implémentation.

Démonstration. Le membre de droite de (6.8) est une combinaison linéaire de polynômes de degré $n$ , donc appartient à $P_{n}$ . En $x = x_{k}$ , la propriété (6.6) annule tous les termes sauf celui d'indice $i = k$ , qui vaut $y_{k} \cdot 1$ :

p_{n} (x_{k}) = i = 0 \sum n y_{i} ℓ_{i} (x_{k}) = i = 0 \sum n y_{i} δ_{ik} = y_{k} .

Les conditions (6.2) sont donc satisfaites, ce qui prouve l'existence; l'unicité a été établie au théorème 6.1.

Pour la base: la famille compte $n + 1$ éléments dans un espace de dimension $n + 1$ , il suffit donc de vérifier qu'elle est libre. Si $\sum_{i} c_{i} ℓ_{i} = 0$ comme polynôme, alors en évaluant en on obtient pour chaque . Les coordonnées d'un polynôme dans cette base sont d'ailleurs lisibles sans aucun calcul: ce sont ses valeurs aux nœuds, .

Le coût, et son défaut

Comptons les opérations. Évaluer $p_{n}$ en un point $t$ par (6.8) demande, pour chaque $i$ , un produit de $n$ fractions: $2 n$ soustractions, $n$ divisions et $n - 1$ multiplications, soit environ $4 n$ opérations. Multiplié par les $n + 1$ valeurs de , cela fait environ

4 n (n + 1) \approx 4 n^{2} op \overset{e}{ˊ} rations par point d’ \overset{e}{ˊ} valuation. (6.9)

C'est un coût en $O (n^{2})$ , sans aucun calcul préparatoire: on peut évaluer $p_{n}$ sans jamais former ses coefficients. Voici l'implémentation directe, en Python:

def ell(i, xs, t):
    """i-eme polynome de Lagrange, evalue en t."""
    p = 1.0
    for k in range(len(xs)):
        if k != i:
            p *= (t - xs[k]) / (xs[i] - xs[k])
    return p


def interpole_lagrange(xs, ys, t):
    return sum(ys[i] * ell(i, xs, t) for i in range(

Onze lignes, deux boucles imbriquées: le coût $4 n^{2}$ de (6.9) se lit directement sur le code.

Le défaut est ailleurs. Les dénominateurs $\prod_{j \neq = i} (x_{i} - x_{j})$ ne dépendent que des nœuds; on peut donc les calculer une fois pour toutes, et chaque évaluation ultérieure ne coûte plus que $O (n)$ — c'est l'idée de la , standard aujourd'hui pour l'interpolation en grand degré. Mais si l'on , tous les changent: chacun gagne un facteur, chaque dénominateur doit être multiplié à nouveau, et il n'existe aucune façon de réutiliser pour obtenir . Le travail est entièrement à refaire. C'est précisément ce que la forme de Newton corrige.

Question 6.3

Écrivez ell(i, xs, t), le $i$ -ième polynôme caractéristique de Lagrange évalué en $t$ , puis interpole_lagrange(xs, ys, t). Le programme affiche $ℓ_{1}$ en $x_{0} = 0$ , puis $ℓ_{1}$ en $x_{1} = 0, 5$ , puis la somme des trois $ℓ_{i}$ en $t = 0, 25$ .

lagrange.py

La forme de Newton et les différences divisées

L'idée de Newton est de choisir une base triangulaire pour les nœuds: le premier élément ne voit que $x_{0}$ , le deuxième $x_{0}$ et $x_{1}$ , et ainsi de suite.

Cherchons $p_{n}$ sous la forme $p_{n} (x) = \sum_{k = 0}^{n} c_{k} ω_{k} (x)$ et écrivons les conditions d'interpolation dans l'ordre. En , tous les sauf s'annulent, donc . En , seuls et survivent, donc et . Le système est : chaque nouvelle condition détermine un nouveau coefficient sans toucher aux précédents. Cette structure est toute la valeur de la forme de Newton, et elle mérite un nom pour les coefficients.

Démonstration. Notons $q_{k}$ le polynôme d'interpolation de $f$ aux nœuds $x_{0}, \dots, x_{k}$ , de sorte que $q_{n} = p_{n}$ . Le polynôme est de degré au plus et s'annule en — car et y prennent tous deux la valeur . Il est donc divisible par :

q_{k} (x) - q_{k - 1} (x) = γ_{k} ω_{k} (x)

pour un scalaire $γ_{k}$ , qui est le coefficient de $x^{k}$ dans $q_{k}$ puisque $q_{k - 1}$ est de degré au plus et unitaire. En sommant ces égalités pour et en ajoutant , on obtient

p_{n} (x) = f (x_{0}) + k = 1 \sum n γ_{k} ω_{k} (x),

qui est bien de la forme (6.12). Il reste à identifier $γ_{k}$ avec $f [x_{0}, \dots, x_{k}]$ .

Montrons-le par récurrence sur $k$ , en établissant du même coup la récurrence (6.11). Pour $k = 0$ , $γ_{0} = f (x_{0}) = f [x_{0}]$ . Supposons le résultat vrai à l'ordre pour toute famille de nœuds. Soient le polynôme interpolant en et celui l'interpolant en . On vérifie directement que

q_{k} (x) = \frac{( x - x _{0} ) s ( x ) - ( x - x _{k} ) r ( x )}{x _{k} - x _{0}} (6.13)

est de degré au plus $k$ et interpole $f$ aux $k + 1$ nœuds: en $x_{0}$ le premier terme s'annule et il reste $r (x_{0}) = f (x_{0})$ ; en il reste ; et pour , , de sorte que le numérateur vaut . Par unicité, (6.13) . Identifions maintenant les coefficients de des deux membres: à gauche , à droite , où et sont les coefficients dominants de et de , égaux par hypothèse de récurrence à et . C'est exactement (6.11), donc .

La symétrie est alors immédiate: $f [x_{0}, \dots, x_{k}]$ est le coefficient dominant du polynôme d'interpolation aux nœuds ${x_{0}, \dots, x_{k}}$ , et cet ne retient pas l'ordre.

Le tableau triangulaire

Les différences divisées se calculent en remplissant un tableau, colonne par colonne, chaque entrée étant la différence des deux entrées voisines de la colonne précédente, divisée par l'écart des nœuds extrêmes concernés. Les coefficients de (6.12) sont les entrées de la diagonale supérieure, c'est-à-dire la première de chaque colonne.

Le coût, et l'évaluation par Horner

Construire le tableau demande, pour la colonne $k$ , $n - k + 1$ entrées, chacune coûtant deux soustractions et une division. Au total

3 k = 1 \sum n (n - k + 1) = \frac{3 n ( n + 1 )}{2} \approx \frac{3}{2} n^{2} op \overset{e}{ˊ} rations, une fois pour toutes. (6.14)

Une fois les coefficients connus, on évalue $p_{n}$ par le schéma de Horner généralisé, qui factorise (6.12) de l'intérieur vers l'extérieur:

p_{n} (x) = c_{0} + (x - x_{0}) (c_{1} + (x - x_{1}) (c_{2} + \dots + (x - x_{n - 1}) c_{n})), (6.15)

soit $n$ multiplications et $2 n$ additions, c'est-à-dire $3 n$ opérations par point — linéaire en $n$ , contre $4 n^{2}$ pour la forme de Lagrange évaluée naïvement. Pour tracer une courbe en 600 points avec $n = 15$ , cela fait $45$ opérations par point au lieu de $960$ .

def differences_divisees(xs, ys):
    """Coefficients c[k] = f[x_0, ..., x_k] de la forme de Newton."""
    c = list(ys)
    for k in range(1, len(xs)):
        for i in range(len(xs) - 1, k - 1, -1):
            c[i] = (c[i] - c[i - 1]) / (xs[i] -

La boucle interne descend délibérément de la fin vers le début: elle écrase c[i] après avoir lu c[i-1], ce qui permet de tout faire sur place, dans un seul tableau de $n + 1$ nombres, sans jamais matérialiser le triangle. Le triangle reste le bon objet pour comprendre; il n'est pas le bon objet pour programmer.

Question 6.4

Vous avez déjà le polynôme d'interpolation de degré $n$ sous forme de Newton et vous recevez un nœud supplémentaire $x_{n + 1}$ . Remettez dans l'ordre les étapes qui donnent $p_{n + 1}$ au moindre coût.

Glissez les éléments pour les mettre dans le bon ordre

Relever la dernière entrée obtenue,

f [x_{0}, \dots, x_{n + 1}]

, comme nouveau coefficient

Calculer de proche en proche les différences divisées d'ordre 1 à

n + 1

de cette nouvelle ligne

Évaluer le nouveau polynôme par le schéma de Horner généralisé

Ajouter la ligne

(x_{n + 1}, f (x_{n + 1}))

au bas du tableau

Ajouter à

p_{n}

le terme

f [x_{0}, \dots, x_{n + 1}] (x - x_{0}) \dots (x - x_{n})

, sans toucher aux coefficients précédents

Les différences divisées sont des dérivées déguisées

La notation $f [x_{0}, x_{1}]$ pour un taux d'accroissement n'est pas un hasard: quand les nœuds se rapprochent, cette quantité tend vers $f^{'} (x_{0})$ . Le résultat général est le suivant.

Démonstration. Posons $g = f - p_{k}$ , où $p_{k}$ est le polynôme interpolant $f$ aux $k + 1$ nœuds. Alors et s'annule en , soit zéros distincts. Le théorème de Rolle (Analyse I, chapitre 8) appliqué sur chacun des intervalles consécutifs donne zéros distincts de ; appliqué à nouveau, zéros de ; et après applications, zéro de , situé dans l'intervalle ouvert délimité par le plus petit et le plus grand des nœuds. Or est de degré au plus de coefficient dominant d'après le théorème 6.3, donc , une constante. De on tire

f^{(k)} (ξ) = p_{k}^{(k)} (ξ) = k! f [x_{0}, \dots, x_{k}],

ce qui est (6.16). $□$

Deux conséquences. D'abord, en faisant tendre tous les nœuds vers un même point $x$ , on obtient la différence divisée confluente $f [x, x, \dots, x] = f^{(k)} (x) / k!$ , ce qui prolonge la définition (6.11) au cas de nœuds répétés — et c'est exactement le mécanisme de l'interpolation d'Hermite de la fin du chapitre. Ensuite, le tableau de différences divisées est un estimateur de dérivées: vérifions-le sur l'exemple 6.2. La différence d'ordre 2 vaut $- 0, 379444$ , donc $f^{''} (ξ) = 2 \times (- 0, 379444) = - 0, 758888$ pour un ; comme , une résolution numérique donne . De même l'ordre 3, , correspond à avec . Les deux existent bien et sont bien à l'intérieur.

Question 6.5

Écrivez differences_divisees(xs, ys), qui renvoie la liste des coefficients de Newton, et evalue_newton(c, xs, t), qui évalue le polynôme par le schéma de Horner généralisé. Le programme affiche les quatre coefficients de l'exemple 6.2, à dix décimales.

newton.py

La formule de l'erreur

Jusqu'ici, rien n'a dit à quelle distance $p_{n}$ se tient de $f$ entre les nœuds. C'est l'objet du théorème suivant, le plus important du chapitre: il structure tout ce qui suit, et il sert encore aux chapitres 7, 8 et 10.

Démonstration. Si $x$ est l'un des nœuds, les deux membres de (6.17) sont nuls et n'importe quel $ξ_{x}$ convient. Fixons donc $x \in [a, b]$ distinct de tous les nœuds; alors $ω (x) \neq = 0$ , et l'on peut poser la constante

K = \frac{f ( x ) - p _{n} ( x )}{ω ( x )},

qui dépend du $x$ fixé mais plus de rien d'autre. Introduisons la fonction auxiliaire de la variable $t$ :

g (t) = f (t) - p_{n} (t) - K ω (t), t \in [a, b] . (6.19)

C'est le cœur de la démonstration, et il vaut la peine de dire pourquoi elle marche: $g$ est construite pour s'annuler une fois de plus que $f - p_{n}$ , et c'est cette annulation supplémentaire qui va, par Rolle, faire descendre l'information jusqu'à la dérivée d'ordre $n + 1$ .

Régularité. $f \in C^{n + 1}$ par hypothèse, $p_{n}$ et $ω$ sont des polynômes, donc $g \in C^{n + 1} ([a, b])$ .

Zéros de $g$ . En chaque nœud $x_{i}$ : $f (x_{i}) - p_{n} (x_{i}) = 0$ par interpolation, et par définition de , donc . En : par le choix de ,

g (x) = f (x) - p_{n} (x) - \frac{f ( x ) - p _{n} ( x )}{ω ( x )} ω (x) = 0.

Comme $x$ est distinct des $n + 1$ nœuds, $g$ possède au moins $n + 2$ zéros distincts dans $[a, b]$ .

Descente de Rolle. Rangeons ces $n + 2$ zéros par ordre croissant. Sur chacun des $n + 1$ intervalles qu'ils délimitent, $g$ est continue, dérivable, et prend la même valeur $0$ aux deux extrémités: le théorème de Rolle (Analyse I, chapitre 8) y fournit un point où $g^{'}$ s'annule. Ces $n + 1$ points sont deux à deux distincts puisqu'ils vivent dans des intervalles d'intérieurs disjoints. Donc possède au moins zéros distincts. Le même raisonnement appliqué à donne au moins zéros distincts de , puis pour , et ainsi de suite. Après applications, possède , que nous notons ; il appartient à l'intervalle ouvert délimité par le plus petit et le plus grand des zéros de départ, donc à .

Calcul de $g^{(n + 1)}$ . Dérivons (6.19) exactement $n + 1$ fois. Le polynôme $p_{n}$ est de degré au plus $n$ , donc $p_{n}^{(n + 1)} \equiv 0$ . Le polynôme est , c'est-à-dire , donc . Il reste

g^{(n + 1)} (t) = f^{(n + 1)} (t) - K (n + 1)! .

Conclusion. En écrivant $g^{(n + 1)} (ξ_{x}) = 0$ , on obtient $K = f^{(n + 1)} (ξ_{x}) / (n + 1)!$ , et en reportant dans la définition de ,

f (x) - p_{n} (x) = K ω (x) = \frac{f ^{(n + 1)} ( ξ _{x} )}{( n + 1 )!} ω (x),

ce qui est (6.17). La majoration (6.18) s'en déduit en majorant séparément les deux facteurs sur $[a, b]$ . $□$

Ce que chacun des trois facteurs apporte

La formule (6.17) est un produit de trois choses de natures très différentes, et savoir laquelle on maîtrise est ce qui permet d'agir.

Le facteur $f^{(n + 1)} (ξ_{x})$ vient de la fonction, et on ne le choisit pas. Il mesure la «raideur» de $f$ à l'ordre $n + 1$ . Pour un polynôme de degré $\leq n$ il est nul, et l'interpolation est exacte — c'est le contrôle que nous avons utilisé dans les exercices de programmation. Pour $sin$ ou , il est borné par quel que soit , et l'interpolation converge magnifiquement. Pour la fonction de Runge, il plus vite que la factorielle, et c'est toute l'histoire de la section suivante. Le point , lui, est inconnu et dépend de : (6.17) est une , pas une formule de calcul; en pratique on passe toujours par la majoration (6.18).

Le facteur $1/ (n + 1)!$ est le seul qui pousse vers le bas, et il est très puissant — quand il gagne. La factorielle écrase tout: $1/16! = 4, 8 \cdot 1 0^{- 14}$ . C'est lui qui fait croire, à première lecture, qu'augmenter le degré améliore toujours l'approximation. Il ne gagne que si $M_{n + 1}$ croît moins vite que , ce qui est vrai pour les fonctions entières d'ordre modeste et faux pour toute fonction ayant une singularité complexe proche de l'intervalle.

Le facteur $ω (x) = \prod_{i} (x - x_{i})$ est le seul que nous contrôlions. Il ne dépend que du placement des nœuds, pas de $f$ . Deux observations le concernant gouvernent la pratique. D'une part, il s'annule aux nœuds: l'erreur est nulle là et maximale entre eux, ce qui donne à l'erreur d'interpolation son allure caractéristique en festons. D'autre part, pour des nœuds équidistants, $∣ ω ∣$ est : à sur , son maximum vaut et il est atteint dans les intervalles extrêmes (en ), alors que dans l'intervalle central il ne dépasse pas , soit . Minimiser par un meilleur choix de nœuds est le sujet de la section sur Tchebychev, et c'est le seul levier dont nous disposons.

Question 6.6

Vous tabulez $sin x$ sur $[0, π /2]$ avec un pas $h$ constant et vous interpolez linéairement entre deux entrées voisines. Quelle est la plus grande valeur de $h$ qui garantisse une erreur inférieure ou égale à $1 0^{- 6}$ ?

Le phénomène de Runge

Nous arrivons au fait le plus contre-intuitif de ce chapitre. Toutes les analyses précédentes suggèrent qu'augmenter le nombre de nœuds améliore l'approximation; le facteur $1/ (n + 1)!$ y pousse fortement. C'est faux en général, et le contre-exemple est devenu si célèbre qu'il porte le nom de son auteur, Carl Runge, qui le publia en 1901.

Interpolons-la sur des nœuds équidistants $x_{i} = - 1 + 2 i / n$ , pour $n = 5$ , $10$ , $15$ , et mesurons l'erreur maximale sur une grille de $2001$ points de $[- 1, 1]$ .

$n$	nœuds équidistants	nœuds de Tchebychev
5	0,433	0,556
10	1,92	0,109
15	2,11	0,083

Lisez d'abord la colonne du milieu. L'erreur maximale, déjà de $0, 433$ pour $n = 5$ , croît avec le degré: $1, 92$ à $n = 10$ , $2, 11$ à $n = 15$ . Elle dépasse donc largement l'amplitude de $f$ elle-même, qui vaut $1$ . Et cela continue: à l'erreur maximale atteint , et à elle vaut . Sur des nœuds équidistants, la suite vers ; elle diverge, et elle diverge vite.

Pourquoi, au regard de la formule de l'erreur

La formule (6.17) explique le phénomène, à condition de regarder les deux facteurs qui varient avec $n$ .

Le polynôme nodal $ω$ des nœuds équidistants décroît: son maximum vaut $6, 92 \cdot 1 0^{- 2}$ à $n = 5$ , $8, 53 \cdot 1 0^{- 3}$ à $n = 10$ , à . Divisé par , c'est même une décroissance foudroyante. Le coupable est donc nécessairement , et on peut le calculer exactement. La décomposition en éléments simples

\frac{1}{1 + 25 x ^{2}} = \frac{1}{2} (\frac{1}{1 - 5 i x} + \frac{1}{1 + 5 i x})

se dérive terme à terme et donne, au point $x = 0$ ,

f^{(k)} (0) = k! 5^{k} cos (\frac{k π}{2}), donc f^{(k)} (0) = k! 5^{k} pour k pair. (6.21)

La dérivée d'ordre $k$ contient la factorielle en facteur. Le rapport qui apparaît dans (6.17) vaut donc

\frac{M _{n + 1}}{( n + 1 )!} \approx 5^{n + 1},

et il croît géométriquement au lieu de décroître. La factorielle, qui devait tout sauver, est exactement annulée. Ce qu'il reste, $5^{n + 1} max ∣ ω ∣$ , croît: $5^{11} \cdot 8, 53 \cdot 1 0^{- 3} \approx 4, 2 \cdot 1 0^{5}$ à .

D'où vient ce facteur $5$ ? De la position des pôles complexes de $f$ . La fonction (6.20) s'annule au dénominateur en $x = \pm i /5$ : elle est parfaitement lisse sur l'axe réel, mais elle possède deux singularités à la distance $0, 2$ de l'intervalle dans le plan complexe, et c'est cette distance qui fixe la vitesse $5 = 1/0, 2$ . Une fonction réelle sans défaut visible peut donc être difficile à interpoler à cause de ce qui lui arrive hors de l'axe réel — l'un des rares endroits de ce cours où l'analyse complexe est indispensable à l'explication. La théorie correspondante (potentiel logarithmique, ellipses de Bernstein) relève d'un cours de théorie de l'approximation, et nous l'admettons ici: nous en retenons seulement le critère qualitatif, plus la singularité complexe la plus proche est près de l'intervalle, plus l'interpolation en grand degré est difficile.

Explorateur 6.1 · Le polynôme d'interpolation quand on ajoute des nœuds

Augmentez le nombre de nœuds avec les nœuds équidistants: l'erreur maximale se met à croître au lieu de décroître, et les oscillations naissent aux bords de l'intervalle. Passez ensuite aux nœuds de Tchebychev et refaites la même montée. Surveillez le dernier affichage, le rapport entre le maximum du polynôme nodal et 2 puissance moins n: il vaut exactement 1,000 pour Tchebychev quel que soit le degré — c'est la propriété minimax — et il s'envole pour les nœuds équidistants.

Degré n (soit n+1 nœuds)10

Répartition des nœudséquidistants

Raideur a de f(x) = 1/(1+ax²)25

Nombre de nœuds

Erreur max ‖f − pₙ‖∞

1,9156

max |Π(x − xᵢ)|

8,53e-3

Ce maximum divisé par 2⁻ⁿ

8,730

Question 6.7

Laquelle de ces fonctions est facile à interpoler en grand degré sur $[- 1, 1]$ avec des nœuds équidistants?

Les nœuds de Tchebychev

Puisque $ω (x) = \prod_{i} (x - x_{i})$ est le seul facteur de (6.17) que nous contrôlions, posons la question directement: quel placement des $n + 1$ nœuds rend $max_{[- 1, 1]} ∣ ω ∣$ le plus petit possible? La réponse est complète et elle est explicite.

Deux lectures de (6.23), toutes deux utiles. La première est algébrique: ce sont les racines d'un cosinus composé. La seconde est géométrique et c'est celle qu'il faut retenir: on place $n + 1$ points également espacés en angle sur le demi-cercle unité, puis on les projette verticalement sur le diamètre. La projection resserre les points près des extrémités et les écarte au centre — exactement l'inverse de ce dont on aurait spontanément envie, et exactement ce qu'il faut.

Démonstration. La seconde partie découle de la première: $ω$ est unitaire de degré $n + 1$ et s'annule aux $n + 1$ racines de $T_{n + 1}$ ; comme $T_{n + 1}$ a pour coefficient dominant , le polynôme est unitaire et possède les mêmes racines, donc lui est égal. D'où (6.25), puisque avec égalité atteinte.

La première partie, l'optimalité (6.24), est admise. La démonstration usuelle est un raisonnement par l'absurde qui exploite l'équi-oscillation: si un polynôme unitaire $q$ vérifiait $∥ q ∥_{\infty} < 2^{1 - m}$ , la différence $T_{m} - q$ , de degré au plus , changerait de signe aux points où atteint alternativement , et aurait donc racines — contradiction. Cet argument est court, mais il est le premier cas particulier du , qui gouverne toute l'approximation uniforme et dont la place est un cours de théorie de l'approximation; nous nous contentons ici du résultat et de sa lecture graphique sur la figure 6.3, où l'égalité de toutes les crêtes se voit directement.

Le passage à un intervalle quelconque

Les formules précédentes vivent sur $[- 1, 1]$ . Pour un intervalle $[a, b]$ quelconque, on utilise le changement de variable affine

x = \frac{a + b}{2} + \frac{b - a}{2} t, t \in [- 1, 1], (6.27)

qui envoie $[- 1, 1]$ sur $[a, b]$ . Les nœuds de Tchebychev de $[a, b]$ sont donc

x_{i} = \frac{a + b}{2} + \frac{b - a}{2} cos (\frac{( 2 i + 1 ) π}{2 ( n + 1 )}), i = 0, \dots, n . (6.28)

Chaque facteur $x - x_{i}$ étant multiplié par $(b - a) /2$ , la borne (6.26) devient

[a, b] max ∣ f - p_{n} ∣ \leq \frac{M _{n + 1}}{( n + 1 )!} \cdot 2 (\frac{b - a}{4})^{n + 1} . (6.29)

La colonne de droite du tableau, et la ligne qu'il ne faut pas cacher

Revenons à la fonction de Runge, maintenant interpolée aux nœuds de Tchebychev, et relisons le tableau complet.

$n$	nœuds équidistants	nœuds de Tchebychev
5	0,433	0,556
10	1,92	0,109
15	2,11	0,083

À $n = 5$ , Tchebychev est moins bon: $0, 556$ contre $0, 433$ . Ce n'est pas une erreur de calcul, c'est le comportement réel, et la raison en est parfaitement lisible. Les deux erreurs maximales sont atteintes en $x = 0$ , au sommet de la fonction, là où aucune des deux familles n'a de nœud. Or les nœuds équidistants les plus proches de $0$ sont $\pm 0, 2$ , tandis que ceux de Tchebychev sont $\pm 0, 2588$ : la répartition de Tchebychev, en resserrant les points vers les bords, dégarnit le centre, et à degré le centre est exactement ce qui compte, puisque c'est là que se trouve toute la structure de .

Cette ligne mérite d'être imprimée plutôt que cachée. L'avantage des nœuds de Tchebychev est asymptotique: il porte sur le comportement quand $n \to \infty$ , pas sur chaque degré pris isolément. À $n = 10$ le rapport est déjà de $17, 6$ en faveur de Tchebychev, à $n = 15$ de $25, 4$ , à $n = 20$ de $3900$ — mais à il est défavorable. Un cours qui ne montrerait que les deux dernières lignes enseignerait un slogan («Tchebychev, c'est mieux») au lieu d'un résultat.

Un dernier point d'honnêteté, et il est important. La borne (6.26) ne démontre pas la convergence observée. Avec (6.21), $M_{n + 1} / (n + 1)! \approx 5^{n + 1}$ , de sorte que (6.26) donne $5 \cdot (5/2)^{n}$ , qui diverge. Et pourtant l'erreur mesurée décroît: à , à , à , à . Le rapport géométrique observé entre et est de par degré. La théorie qui prédit ce nombre est celle des ellipses de Bernstein: la singularité se trouve sur l'ellipse de paramètre , et la convergence de l'interpolation de Tchebychev y est géométrique de raison . La coïncidence des deux nombres, mesuré contre prédit, est excellente — mais elle relève d'un résultat que ce cours , faute d'y disposer de l'analyse complexe nécessaire. Retenez-en que (6.18) est une , parfois très pessimiste, et qu'un majorant qui diverge ne prouve jamais qu'une suite diverge.

Question 6.8

Écrivez noeuds_tchebychev(a, b, n), qui renvoie la liste des $n + 1$ nœuds de Tchebychev de l'intervalle $[a, b]$ , d'après la formule (6.28). Le programme affiche les quatre nœuds obtenus pour $a = - 1$ , $b = 1$ et $n = 3$ , à six décimales.

tchebychev.py

Les deux autres sorties

Les nœuds de Tchebychev règlent le problème lorsqu'on choisit où mesurer. Ce n'est pas toujours le cas: une centrale d'acquisition échantillonne à intervalles réguliers, une table publiée impose ses entrées. Deux autres stratégies existent alors.

Baisser le degré et découper: l'interpolation par morceaux

La formule (6.18) contient $M_{n + 1}$ et $max ∣ ω ∣$ , tous deux pris sur l'intervalle entier. Rien n'oblige à utiliser un seul polynôme sur tout $[a, b]$ . Découpons $[a, b]$ en sous-intervalles de longueur $h$ et interpolons de degré 1 sur chacun: la borne devient $h^{2} M_{2} /8$ sur chaque morceau, donc sur l'ensemble, et elle tend vers zéro comme quand on raffine — sans aucune condition sur les dérivées d'ordre élevé de , puisque seule intervient. La divergence de Runge disparaît, parce que c'est le degré, et non le nombre de points, qui la provoquait.

Le prix est la régularité: la ligne brisée obtenue est continue mais pas dérivable aux nœuds. Les splines cubiques rétablissent la dérivabilité, et même la continuité de la dérivée seconde, en recollant des polynômes de degré 3 avec les bonnes conditions de raccord. Elles sont l'objet du chapitre 7, avec l'approximation au sens des moindres carrés — les deux réponses aux deux défauts de l'interpolation polynomiale globale, le phénomène de Runge d'une part, le bruit des données d'autre part. Sur la fonction de Runge elle-même, la spline cubique sur les mêmes nœuds équidistants est excellente là où le polynôme de degré 15 se déchaîne: c'est la comparaison qui ouvre ce chapitre 7.

Imposer aussi les dérivées: l'interpolation d'Hermite

L'autre direction consiste à demander davantage en chaque nœud, plutôt que davantage de nœuds.

En pratique on ne repart pas de zéro: il suffit de dédoubler chaque nœud dans le tableau des différences divisées, en remplaçant la différence indéfinie $f [x_{i}, x_{i}]$ par sa valeur confluente $f^{'} (x_{i})$ , conformément au théorème 6.4. Tout le reste du tableau et le schéma de Horner (6.15) fonctionnent sans modification. La formule d'erreur, établie par le même argument de fonction auxiliaire que le théorème 6.5, devient

f (x) - H (x) = \frac{f ^{(2 n + 2)} ( ξ _{x} )}{( 2 n + 2 )!} i = 0 \prod n (x - x_{i})^{2}, (6.31)

avec les carrés qui traduisent la multiplicité double des nœuds. L'ordre est donc doublé pour un même nombre de nœuds.

Ce que cela coûte: il faut connaître les dérivées $f^{'} (x_{i})$ . Sur une fonction analytique c'est souvent gratuit; sur des données mesurées, c'est le plus souvent impossible, et les estimer par différences finies revient à réintroduire l'erreur qu'on cherchait à éliminer. C'est pourquoi l'interpolation d'Hermite est surtout employée là où les pentes sont naturellement disponibles: le cas le plus important est la spline cubique d'Hermite par morceaux du chapitre 7. La même idée — un interpolant qui respecte aussi les dérivées — sert encore à reconstituer une solution continue entre les pas d'une méthode à un pas: le schéma a déjà évalué les pentes pour avancer, et il ne coûte donc rien de les réutiliser pour interpoler à l'intérieur du pas.

Question 6.9

Assemblez les briques du chapitre: écrivez erreur_max(noeuds), qui interpole la fonction de Runge sur ces nœuds et renvoie l'erreur maximale mesurée sur une grille de 2001 points de l'intervalle allant de $- 1$ à $1$ . Le programme affiche la première ligne du tableau du chapitre, pour $n = 5$ .

runge.py

Synthèse

Le problème d'interpolation impose le passage exact par $n + 1$ points de nœuds distincts, et il admet une solution polynomiale de degré $\leq n$ et une seule (théorème 6.1): la différence de deux solutions aurait $n + 1$ racines pour un degré au plus $n$ . Il se distingue de l'approximation au sens des moindres carrés, qui minimise un écart au lieu de l'annuler et qui est le sujet du chapitre 7; on interpole des données fiables, on approche des données bruitées.
La matrice de Vandermonde prouve l'existence mais ne doit jamais servir à calculer: son conditionnement $κ_{\infty}$ vaut $1, 25 \cdot 1 0^{4}$ à et à sur les nœuds équidistants de , et la résolution en double précision ne laisse plus aucun chiffre correct au-delà du degré 18.

Série d'exercices du chapitre 6Exercice 1 sur 5

Question 6.10

Combien de polynômes de degré au plus 4 passent par les cinq points $(0, 1)$ , $(1, 2)$ , $(2, 5)$ , $(3, 10)$ , $(4, 17)$ ?

Exercices

Vous pouvez afficher le corrigé directement sous chaque énoncé après avoir cherché la solution.

Exercice 6.1 · Deux formes, un seul polynôme

On considère les données $(0, 1)$ , $(1, 3)$ , $(3, 2)$ .

Écrire les trois polynômes caractéristiques de Lagrange $ℓ_{0}, ℓ_{1}, ℓ_{2}$ et en déduire sous forme développée.

Exercice 6.2 · Dimensionner une interpolation linéaire

On veut tabuler $f (x) = e^{x}$ sur $[0; 0, 5]$ à pas constant $h$ , avec interpolation linéaire entre deux entrées, et garantir une erreur inférieure ou égale à $1 0^{- 4}$ .

Exercice 6.3 · Le nœud qui manque au centre

On interpole $f (x) = 1/ (1 + 25 x^{2})$ sur $[- 1, 1]$ avec $n = 5$ , une fois sur les nœuds équidistants $- 1; - 0, 6; - 0, 2; 0, 2; 0, 6; 1$ , une fois sur les nœuds de Tchebychev , , .

Exercice 6.4 · Où s'arrête la base canonique

On reprend la matrice de Vandermonde des nœuds équidistants $x_{i} = i / n$ de $[0, 1]$ , dont le conditionnement vaut $4, 00 \cdot 1 0^{8}$ pour et pour .

Exercice 6.5 · Trois propriétés des différences divisées

Soit $f$ une fonction et $x_{0}, \dots, x_{n}$ des nœuds deux à deux distincts.

Montrer que pour tout entier $k$ avec $0 \leq k \leq n$ , pour tout . Le cas redonne (6.7).

Références

Quarteroni, A., Sacco, R. et Saleri, F., Méthodes numériques — algorithmes, analyse et applications, Springer, Milan, chap. 8 (interpolation polynomiale, différences divisées, interpolation de Tchebychev).
Rappaz, J. et Picasso, M., Introduction à l'analyse numérique, Presses polytechniques et universitaires romandes, Lausanne, chap. 4.
Burden, R. L. et Faires, J. D., Numerical Analysis, Cengage, Boston, chap. 3 (Lagrange, Newton, Hermite).
Trefethen, L. N., Approximation Theory and Approximation Practice, SIAM, Philadelphie (la référence moderne sur les nœuds de Tchebychev, les ellipses de Bernstein et la forme barycentrique).
Berrut, J.-P. et Trefethen, L. N., «Barycentric Lagrange Interpolation», SIAM Review, vol. 46, n° 3, 2004.
Schatzman, M., Analyse numérique — une approche mathématique, Dunod, Paris, chap. 5.

5	$1, 25 \cdot 1 0^{4}$	$2, 8 \cdot 1 0^{- 12}$	$7, 9 \cdot 1 0^{- 14}$
10	$4, 00 \cdot 1 0^{8}$	$8, 9 \cdot 1 0^{- 8}$	$9, 9 \cdot 1 0^{- 9}$
15	$1, 25 \cdot 1 0^{13}$	$2, 8 \cdot 1 0^{- 3}$	$4, 1 \cdot 1 0^{- 4}$
20	$4, 13 \cdot 1 0^{17}$	$9, 2 \cdot 1 0^{1}$	$2, 5$

$x_{i}$	$f [x_{i}]$	ordre 1	ordre 2	ordre 3
0,0	1,000000
0,5	0,778801	−0,442398
1,0	0,367879	−0,821843	−0,379444
1,5	0,105399	−0,524960	0,296882	0,450884

Objectifs du chapitre

Le problème d'interpolation

Existence et unicité

La forme de Lagrange

Le coût, et son défaut

La forme de Newton et les différences divisées

Le tableau triangulaire

Le coût, et l'évaluation par Horner

Les différences divisées sont des dérivées déguisées

La formule de l'erreur

Ce que chacun des trois facteurs apporte

Le phénomène de Runge

Pourquoi, au regard de la formule de l'erreur

Les nœuds de Tchebychev

Le passage à un intervalle quelconque

La colonne de droite du tableau, et la ligne qu'il ne faut pas cacher

Les deux autres sorties

Baisser le degré et découper: l'interpolation par morceaux

Imposer aussi les dérivées: l'interpolation d'Hermite

Synthèse

Exercices

Références

La forme de l'erreur pour n = 1

La borne complète

Le pas et le nombre d'entrées

Ce que gagnerait un degré de plus