TangenteX.com

Introduction

Qu'est-ce qu'une interpolation

Soit une fonction réelle continue sur un intervalle donné f(x), fonction que je suppose inconnue mais dont je connais un certain nombre de valeurs (xi,yi) dans l'intervalle.

Pour des raisons que nous verrons plus loin, je cherche à construire une fonction f*(x), telle que pour tous les points connus (xi,yi) de f(x), j'ai f*(xi) = f(xi). En d'autres termes, je veux construire une fonction f*(x) qui passe par les mêmes points que la fonction f(x), pour un ensemble fini de points donnés sur l'intervalle. Si l'ensemble était infini sur l'intervalle, nous aurions bien sur f*(x) = f(x).

Les points (xi,yi) sont nommés "pivots" ou "noeuds" dans la littérature, selon les auteurs.

Par définition, l'interpolation est donc une opération différente d'une régression. Il ne s'agit pas de construire une fonction qui passe à une distance moyenne minimum non nulle des points connus de l'intervalle, mais qui passe par tous ces points.

Dans la mesure où nous construisons la fonction f*(x), je peux la choisir de la forme qui me semble la plus simple et la plus pratique. Classiquement, on emploie la forme polynomiale, c'est à dire que je choisis \(\displaystyle  f^*(x) = \sum_{i=0}^n a_i x^{i} \), un polynôme de degré n, dont les ai sont les coefficients. On parle alors d'interpolation polynomiale. Mais il est possible de choisir d'autres formes, la plus répandue étant la forme trigonométrique.

Je vous ai parlé "d'intervalle donné". Pratiquement cet intervalle est défini par l'ensemble des points connus. L'intervalle d'interpolation est borné par le minimum des xi et par le maximum des xi.

Pourquoi interpoler une fonction

Historiquement, l'usage de l'interpolation était justifié par la faiblesse des moyens de calcul. Pour des fonction compliquées, d'usage courant, les gaussiennes ou les fonctions de Bessel par exemple, il était courant d'utiliser des tables numériques qui donnaient les valeurs de la fonction en un nombre réduit de points, qu'il fallait calculer "à la main". Et lorsqu'on voulait connaitre la valeur de la fonction entre deux points tabulés, on procédait à une interpolation.

En physique, beaucoup de grandeurs caractéristiques comme la conductitivité ou la densité des corps dépendent de plusieurs variables, température, pression par exemple. Les expérimentateurs établissaient des tables de mesures, tables qui servaient de base aux interpolations nécessaires pour connaitre la valeur d'une caractéristique en dehors des points expérimentaux.

Cet usage de l'interpolation a largement régressé, mais n'a pas totalement disparu. On trouve encore des ouvrages entiers de tables de caractéristiques et de fonctions tabulées (en stastistique en particulier). J'ai d'ailleurs utilisé cet usage de l'interpolation dans la page suivante de TangenteX.

Cependant, l'usage principal de l'interpolation en physique numérique, est un peu différent. Pour nous, il s'agit de trouver une fonction, de préférence la moins compliquée possible, qui puisse représenter le plus fidèlement possible une série de points expérimentaux. C'est indispensable lorsque nous voulons procéder à l'intégration ou à la dérivation sur une série de mesures ou encore une analyse de Fourier. Bref, lorsque nous voulons traiter mathématiquement une série de mesures, il est pratique de posséder une fonction analytique qui la représente.

Dans une autre page de TangenteX, nous avons déjà abordé ce problème, un peu différemment, en procédant à une régression linéaire. Le principe de la régression, linéaire ou autre, est d'approcher la série de mesures par une fonction qui ne passe pas forcément par tous les points de mesures, mais qui se rapproche en moyenne le plus près de ces points.

Comment établir une fonction d'interpolation

Dans la suite de cette page, nous allons découvrir plusieurs moyens de construire des fonctions d'interpolations. Ces fonctions sont des fonctions polynomiales. Nous examinerons d'abord les méthodes classiques de Lagrange et de Newton, pour lesquelles je développerai le principe et proposerai un algorithme. Puis nous aborderons les fonctions de scipy qui permettent l'interpolations des fonctions selon d'autres méthodes comme les splines.

Tous les scripts mentionnés dans cette page sont téléchargeables dans la bibliothèque de codes de TangenteX.

Les polynômes de Lagrange

Définition

Je cherche un polynôme d'interpolation unique p(x), qui passe par n+1 pivots de x0 à xn. Je peux écrire (voir la démonstration en cours d'algèbre) :

\( \displaystyle p(x) = \sum_{i=0}^n L_i(x)f(x_i) \)

où Li(x) est le polynôme élémentaire (on dit aussi nodale) de Lagrange de degré n et f(xi) l'ordonnée des points connus. Le polynôme nodale Li(x) s'annule pour tous les points x de l'intervalle, sauf pour xi.

On démontre que le polynôme d'interpolation de Lagrange s'écrit :

\( \displaystyle p(x) = \sum_{i=0}^n L_i(x)f(x_i) =  \sum_{i=0}^n f(x_i) \prod_{k=0, k \ne i} \dfrac{x - x_k}{x_i - x_k} \)

L'indice i fixe le degré du polynôme d'interpolation. Si vous possédez n points connus, le degré de votre polynôme d'interpolation sera de n-1.

Calcul numérique

Imaginons que je dispose d'un ensemble de 7 points (2,1), (3,4), (4,6), (5,7) (6,2), (7,4.5) et (8,3.5). Je veux calculer le polynôme d'interpolation sur l'intervalle [2,8].

Je vous propose un algorithme basique qui implémente "bêtement" la définition du polynôme d'interpolation de Lagrange. Il est disponible dans le script InterpolLagrangeDLE.py :

from scipy import array, linspace, zeros
from matplotlib.pylab import figure, plot, grid, title
# fonction de calcul du polynome d'interpolation de Lagrange
def InterpolLagrange(x,y,X):
    n = min(len(x), len(y))
    m = len(X)
    yh = zeros(m)
    pl = zeros(n)
    for j in range(m):
        yh[j] = 0
        for i in range(n):
            pl[i] = 1.0
            for k in range(n):
                if i != k:
                    pl[i] = pl[i]*(X[j] - x[k])/(x[i] - x[k])
            yh[j] = yh[j] + y[i]*pl[i]           
    return yh   
# définition des pivots d'interpolation
x = array([2.0,3.0,4.0,5.0,6.0,7.0,8.0])
y = array([1.0,4.0,6.0,7.0,2.0,4.5,3.5])
# calcul du polynome de Lagrange
X = linspace(min(x),max(x), 200)
Y = InterpolLagrange(x,y,X)
# tracé des résultats
fig1 = figure(figsize=(10,8))
grid()
title(u'Interpolation de Lagrange',fontsize=14)
# tracé des points
plot(x,y,'ro')
# tracé de la fonction d'interpolation
plot(X,Y,'b')

Comme vous le constatez, c'est l'implémentation basique de la définition ! Il est possible de faire mieux en termes de temsp de calcul. Ceci dit, l'algorithme d'interpolation de Lagrange n'est pas terrible en matière d'optimisation numérique. C'est d'abord un algorithme pédagogique. Dans la vraie vie, on utilise d'autres méthodes.

Avec les pivots définis ci-dessus, nous obtenons la courbe d'interpolation suivante :

Interpolation de Lagrange avec algo DLE

Nous constatons que, comme attendu, la fonction d'interpolation passe par tous les pivots.

Les limites de l'interpolation de Lagrange

L'interpolation de Lagrange possède ses limites. Tout d'abord, comme toute méthode numérique, elle comporte une marge d'erreur qu'il faut savoir estimer. Mais au delà de cette marge d'erreur, la méthode présente parfois des problèmes de convergence très ennuyeux.

Estimation de l'erreur d'interpolation

Un théorème d'analyse numérique nous dit que si je considère une fonction continue n+1 fois dérivable, avec n le degré du polynôme d'interpolation, alors si j'appelle En(x) l'erreur d'interpolation au point x appartenant à l'intervalle d'interpolation, j'ai :

\( E_n(x) = f(x) - \prod_n f(x) = \dfrac{f^{(n+1)}(\eta)}{(n + 1)!} L_{n+1} (x) \) avec \( \eta \in [a,b] \)

Un corollaire de ce théorème permet de borner cette erreur sur un intervalle [a,b] avec des pivots équirépartis par :

\( E_n(x) = max \left(|f(x) - \prod_n f(x) | \right) \le \dfrac{1}{4(n+1)} \left( \dfrac{b-a}{n}  \right)^{n+1} max \left(f^{(n+1)} (x)\right) \)

Pour illustrer l'évaluation de l'erreur d'interpolation, prenons l'exemple très classique de l'interpolation de la fonction sin(x) entre [0, 3pi]. Je vais étudier l'erreur d'interpolation commise sur cet intervalle, en observant que sur cet intervalle la valeur \( max \left(f^{(n+1)} (x)\right) \) est majorée par 1 et que donc \( E_n(x) \le \dfrac{1}{4(n+1)} \left( \dfrac{b-a}{n}  \right)^{n+1} \).

Sur le graphe ci-dessous, j'ai tracé à l'aide du script InterpoLagrangeErreur.pyl la fonction sinus entre 0 et 3pi en rouge, et les différents polynômes d'interpolation :

  • en bleu, le polynôme de degré 2,
  • en vert, le polynôme de degré 4,
  • en jaune, le polynôme de degré 6,
  • en noir, le polynôme de degré 8.

Plus le degré du polynôme augmente et plus la convergence avec la fonction sinus est bonne. D'ailleurs, vous aurez du mal à distinguer le polynôme de degré 8 (en noir) et la fonction sinus (en rouge). Cependant, méfiance, nous verrons plus loin qu'on ne peut pas augmenter impunément le degré du polynôme d'interpolation.

Interpolation fonction sinus

Pour les polynômes de degré 2, 4, 6 et 8, le script calcule la valeur max de l'écart entre le polynôme d'interpolation et la fonction sinus sur l'intervalle [0,3pi] et la majoration théorique donnée par la formule \( E_n(x) \le \dfrac{1}{4(n+1)} \left( \dfrac{b-a}{n}  \right)^{n+1} \).

Erreur Interpolation fonction sinus

La majoration théorique est tracée en rouge sur le graphe ci-dessus et l'erreur calculée en vert. Vous remarquez que l'erreur d'interpolation est très inférieure à la majoration théorique et que les deux courbes tendent vers 0 lorsque n augmente.

Pour information, voici le bout de code qui effectue les calculs d'erreurs :

# fonctions de calcul d'erreur
def ErreurTheorique(a,b,n):
    et = 1./(4*(n +1))*((b - a)/n)**(n+1)
    return et
   
def ErreurCalcul(Y,YP):
    ec = max(abs(Y - YP))
    return ec

# Evaluation de l'erreur
ErrT = [ErreurTheorique(a,b,n) for n in range(2,10,2)]
ErrC = [ErreurCalcul(Y,PInt[n](X)) for n in range(2,10,2)]
Défaut de convergence de la méthode

Vous aurez sans doute noté que dans les calculs ci-dessus, j'ai toujours utilisé des pivots équirépartis sur l'axe des abscisses. Techniquement, c'est la méthode la plus simple lorsqu'il s'agit d'interpoler une fonction, mais cette méthode peut présenter des inconvénients majeurs.

Cette équirépartition engendre une instabilité numérique lorsque le nombre de noeuds devient relativement élevé (à partir de 20 environ), car la convergence du produit \( \displaystyle \prod_{k=0, k \ne i} f \) vers f(x) n'est plus garantie. Un théorème, le théorème de Faber, permet de démontrer qu'il existe toute une variété de fonctions pour lesquelles l'erreur croît avec le nombre de pivots. Pour pallier ce défaut, on peut utiliser différentes techniques :

  • Modifier la répartition des pivots, en choisissant une répartition qui améliore la convergence sur l'intervalle d'interpolation. C'est la méthode que nous allons voir ci-après.
  • Réaliser une interpolation par morceaux sur l'intervalle, c'est à dire qu'on découpe l'intervalle en sous-intervalles, et on procède à l'interpolation de Lagrange sur ces sous-intervalles. Je vous laisse découvrir cette méthode en faisant l'exercice. L'interpolation par splines est basée sur une interpolation par morceaux.

Considérons par exemple la fonction de Cauchy \( \displaystyle f(x) = \dfrac{1}{x^2 + 1} \). C'est une fonction continue, sans problème apparent sur \( \mathbb{R} \). Essayons de l'interpoler sur l'intervalle [-5, 5 ] avec la méthode de Lagrange, en choisissant des pivots équirépartis. Le script Python RungeLagrange.py va faire ça pour nous :

# définition des pivots d'interpolation équirépartis
a = -5.0                           # borne inf de l'intervalle d'interpolation
b = 5.0                            # borne sup de l'intervalle d'interpolation

# calcul du polynome de Lagrange avec des pivots équirépartis
n = 10                             # nombre de pivots d'interpolation 
x = linspace(a,b,n)
y = f(x)
X = linspace(min(x),max(x), 200)
Y = f(X)
YL = InterpolLagrange(x,y,X)

Sur le même graphe, j'ai tracé les pivots (points rouges), la fonction de Cauchy en vert, et la fonction d'interpolation en bleu :

Interpolation pivots équirépartis

Plus on se rapproche des bords de l'intervalle et plus la fonction d'interpolation diverge de la fonction interpolée. Ce comportement est connu sous le nom de phénomène de Runge. Comme pour la fonction sinus, vous pourriez être tenté d'augmenter le nombre de pivots pour réduire la divergence : essayez pour voir ! Plus le nombre de pivots équirépartis est grand et plus la divergence est marquée.

Pour améliorer la convergence de la fonction d'interpolation, on n'utilise plus des pivots équirépartis mais des pivots répartis sur l'intervalle d'une manière particulière. Les abscisses de ces pivots sont les zéros du polynôme de Tchebyshev (Les auteurs anglo-saxons écrivent "Chebyshev") de degré n+1, dans notre exemple de degré 11. Ils sont définis par :

\( \displaystyle x_ i = \dfrac{a + b}{2} + \dfrac{b - a}{2} \cos\left(\pi \dfrac{2i + 1}{2n + 2}\right) \) pour i = 0 à n

Interpoler notre fonction de Cauchy en utilisant les pivots de Tchebyshev est simple. Il suffit de calculer les coordonnées de ces pivots à l'aide de la fourmule ci-dessus et de lancer l'interpolation de Lagrange avec ces pivots :

p = 30
xC = [((a+b)/2) + ((b-a)/2)*cos(pi*(2*i + 1)/(2*p + 2)) for i in range(0, p)]
xC = array(xC)
YC = f(xC)
YLC = InterpolLagrange(xC,YC,X)

Nous constatons sur le graphe ci-dessous que la fonction interpolée converge beaucoup mieux avec la fonction de Cauchy. Presque parfaitement, à part une divergence faiblement étendue mais importante au début de l'intervalle. L'augmentation du nombre de pivots (valeur de p) diminue cette divergence, mais au détriment du temps de calcul. C'est à vous de voir...

Interpolation pivots Chebyshev

La fonction de Cauchy est une fonction continue. Si vous tentez d'interpoler une fonction discontinue ou continue par morceaux, vous constaterez un autre phénomène, connu sous le nom de phénomène de Gibbs. Nous verrons comment nous en débarasser une autre fois.

Les polynômes de Newton

Comme je le disais plus haut, l'interpolation de Lagrange n'est pas terrible en matière de précision et de temps de traitement. Les polynômes de Newton apportent quelques améliorations en temps de calcul.

Définition

Il s'agit d'écrire autrement le polynôme de Lagrange pour faire apparaître une forme récurrente qui allégera les calculs. Avec les notations habituelles, nous cherchons à construire la relation :

\( \prod_n(x) =   \prod_{n-1}(x) + q_n(x) \)

avec qn(x), un polynôme à déterminer.

Nous avons \( q_n(x) =  \prod_n(x) - \prod_{n-1}(x) \), qui est nul pour i variant de 0 à n-1, et donc qn(x) est de la forme :

\(\displaystyle  q_n(x) = a_n  \sum_{i=0}^{n-1} (x - x_i) \)

Les coefficients an sont nommés nième différences divisées de Newton, de la forme :

\( a_n = \dfrac{f(x_n) - \prod_{n-1}(x_n)}{L_n(x_n)}\)

Si je note an = f[x0, ..., xn], an se calcule par récurrence selon la formule :

\( a_n = \dfrac{f[x_1, x_n] - f[x_0, x_{n-1}] }{x_n - x_0} \)

d'où la forme finale du polynôme d'interpolation de Newton :

\( \displaystyle \prod_{n}(x) = \sum_{i=0}^{n} a_n L_i(x) \)

L'intérêt de la forme de Newton du polynôme d'interpolation est essentiellement numérique. Lorsqu'on veut ajouter un pivot, point n'est besoin de recalculer tout le polynôme : il suffit d'ajouter un terme supplémentaire \(a_{n+1}L_{n+1} \) au polynôme \( \prod_{n}(x) \), ce qui simplifie largement l'algorithme de calcul.

Calcul numérique

Comme pour l'interpolation de Lagrange, je vous propose un script InterpolNewton.py qui implémente "bêtement" la définition donnée ci-dessus. Nous retrouvons la souplesse apportée par l'usage de la récurrence dans le calcul des coefficients an, que j'ai nommé dans le code CoeffNewton.

from scipy import copy, linspace, array
from matplotlib.pylab import figure, plot, grid, title
# fonction de calcul des coefficients du polynôme d'interpolation de Newton
def CoeffNewton(x,y):
    m = len(x)
    x = copy(x)
    yh = copy(y)
    for k in range(1,m):
        yh[k:m] = (yh[k:m] - yh[k-1])/(x[k:m] - x[k-1])
    return yh
# fonction de calcul du polynôme d'interpolation
def InterpolNewton(x, y, X):
    Y = CoeffNewton(x, y)
    n = len(x) - 1
    p = Y[n]
    for k in range(1,n+1):
        p = Y[n-k] + (X -x[n-k])*p
    return p
# définition des pivots d'interpolation
x = array([2.0,3.0,4.0,5.0,6.0,7.0,8.0])
y = array([1.0,4.0,6.0,7.0,2.0,4.5,3.5])
# calcul du polynome de Newton
X = linspace(min(x),max(x), 200)
Y = InterpolNewton(x,y,X)
# tracé des résultats
fig1 = figure(figsize=(10,8))
grid()
title(u'Interpolation de Newton',fontsize=14)
# tracé des points
plot(x,y,'ro')
# tracé de la fonction d'interpolation
plot(X,Y,'b')

Avec les mêmes pivots équirépatis utilisés plus haut, j'obtiens la courbe d'interpolation suivante :

Interpolation de Newton avec algo DLE

La courbe obtenue est identique à celle obtenue avec l'interpolation de Lagrange, mais le temps de calcul est un peu réduit, bien que cela ne soit pas très sensible.

Une remarque : le polynôme de Newton étant une autre manière d'écrire le polynôme d'interpolation de Lagrange, il souffre des mêmes problèmes de convergence avec des pivots équirépartis que ce dernier...

Quelques  méthodes Scipy

Scipy propose quelques méthodes intéressantes d'interpolation, dont une fonction d'interpolation par polynômes de Lagrange et une fonction générique qui interpole principalement par des splines.

Interpolation polynomiale de Lagrange

Scipy fournit une routine toute faite pour calculer les coefficients du polynôme d'interpolation de Lagrange pour un ensemble de points donnés. Il s'agit de la routine scipy.interpolate.lagrange, dont la documentation est disponible ici.

Le script InterpolLagrange.py me permet de calculer les coefficients du polynôme et de tracer la fonction d'interpolation :

#import des packages
from scipy import array, linspace
from numpy.polynomial.polynomial import Polynomial
from scipy.interpolate import lagrange
from matplotlib.pylab import figure, plot, grid, title
# définition des pivots d'interpolation
x = array([2.0,3.0,4.0,5.0,6.0,7.0,8.0])
y = array([1.0,4.0,6.0,7.0,2.0,4.5,3.5])
# calcul et impression des coefficients du polynome de Lagrange
P = lagrange(x,y)
Coef = Polynomial(P).coef
print Coef # dans l'ordre du coefficent le plus haut vers a0
# tracé des résultats
fig1 = figure(figsize=(10,8))
grid()
title(u'Interpolation de Lagrange',fontsize=14)
# tracé des points
plot(x,y,'ro')
# tracé de la fonction d'interpolation
X = linspace(min(x),max(x), 200)
Y = P(X)
plot(X,Y,'b')

Après exécution, le script retourne la valeur des coefficients du polynôme d'interpolation :

[ -9.23611111e-02 2.68958333e+00 -3.13715278e+01 1.87052083e+02 -6.00286111e+02 9.83508333e+02 -6.39500000e+02]

et la courbe de la fonction d'interpolation :

Interpolation de Lagrange avec scipy

Attention, comme la note de scipy le mentionne, l'interpolation de Lagrange n'est pas stable numériquement à partir d'une vingtaine de points (voir ci-dessus). Au dessus de ce nombre, il faut choisir une autre méthode parmi celles abordées ci-dessous ou procéder à une interpolation par morceaux.

La méthode interp1d

La méthode d'interpolation interp1d est documentée ici dans sa généralité et pour les détails. Il s'agit d'une méthode générale d'interpolation des fonctions réelles y = f(x), pour laquelle on peut choisir, à l'aide d'un paramètre, la méthode particulière d'interpolation :

  • linéaire (en bleue sur le graphe ci-dessous),
  • spline d'ordre 0, 1, 2 ou 3 (respectivement en rouge,vert, jaune et noir)

Le script Interpol1d.py montre l'usage de la fonction interp1d pour les différentes méthodes d'interpolation :

La méthode interp1d de scipy

L'interpolation par spline d'ordre 1 correspond à l'interpolation linéaire. L'interpolation la plus proche de Lagrange et Newton est l'interpolation par spline quadratique (ordre 2).

Revenons à notre fonction de Cauchy et essayons de l'interpoler avec la fonction interp1d. Le script CauchyInterp1d.py donne ce résultat :

Interpolation de la fonction de Cauchy avec interp1d n = 10

Cette fonction d'interpolation est obtenue avec une spline d'ordre 3 (cubique) avec n = 10. A noter que la fonction inter1d option spline n'est pas sensible à l'augmentation du nombre de pivots. Avec n = 30, j'obtiens la fonction d'interpolation :

Interpolation de la fonction de Cauchy avec interp1d n = 30

La fonction d'interpolation converge presque parfaitement avec la fonction de Cauchy.

Les précautions à prendre

Il me reste à vous mettre en garde contre les pièges de l'interpolation. Considérons par exemple le cas de la fonction \( \dfrac{\sin(x)}{1 - x} \), que je veux interpoler dans l'intervalle [-2, 3] avec 10 pivots. J'utilise la fonction interp1d avec l'option spline cubique, mais à vrai dire la méthode importe peu.

Considérons le résultat obtenu avec le script InterpolPiege.py :

Pièges de l'interpolation

Comme l'étude de la fonction le montre, on observe une asymptote en x = 1. L'équirépartition de mes pivots me donne une fonction d'interpolation complètement fausse. Si j'essaye d'interpoler avec des pivots de Tchebyshev, même résultat.

Lorsque vous suspectez que votre fonction interpolée est "pathologique", il vous appartient de définir les pivots de façon adéquate pour converger au mieux. Ce n'est pas toujours facile et parfois impossible. Aussi, prudence avec les méthodes d'interpolation....

Contenu et design par Dominique Lefebvre - www.tangenteX.com février 2020    Licence Creative Commons    Contact :