Python pour le calcul scientifique/Interpolation, extrapolation et lissage

Un livre de Wikilivres.
Sauter à la navigation Sauter à la recherche

Rappelons que dorénavant les programmes commencent tous par :

import numpy as np
import matplotlib.pyplot as plt

Interpolation à une dimension[modifier | modifier le wikicode]

Trois procédés d'interpolation.
Interpolation linéaire avec trois manières d'extrapoler :
  • avec une valeur constante ;
  • extrapolation périodique ;
  • extrapolation linéaire.

Nous disposons de n points M de coordonnées (xpi, ypi)0 ≤ in – 1. Nous supposons que ces points décrivent une fonction continue ƒ et nous voulons connaître une approximation de ƒ(x) pour un point x quelconque. Les méthodes d'interpolation utilisent des fonctions qui passent exactement par les points (xpi, ypi). Entre ces points, c'est-à-dire pour une valeur de x comprise entre deux valeurs xpi et xpi + 1, on utilise typiquement une des solutions suivantes :

  • on détermine le xpj le plus proche de x et l'on prend y = ypj ce qui donne une fonction en escalier ;
  • on trace un segment de droite entre (xpi, ypi) et (xpi + 1, ypi + 1), le point (x, y) est situé sur ce segment de droite ; c'est la méthode d'interpolation linéaire, la fonction est affine par parties, son graphe est un polygone ouvert ;
  • on utilise le polynôme de degré n – 1 passant pour tous les points, le point (x, y) est situé sur ce polynôme ;
  • on utilise des splines (cerces), des polynômes par partie, en général de degré 2 ou 3, dont les dérivées sont raccordées.

En dehors de l'intervalle [xp 0 ; xp n – 1], on peut extrapoler les valeurs de la manière suivante :

  • pas d'extrapolation, la fonction retourne des valeurs NaN ;
  • extrapolation par une valeur constante, soit la valeur la plus proche (xp 0 à gauche, xp n – 1 à droite), soit une valeur fixée arbitrairement (éventuellement 0) ;
  • soit en prolongeant la fonction avec la loi sur le segment [xp 0 ; xp 1] à gauche, [xp n – 2 ; xp n – 1] à droite ;
  • soit en considérant que la fonction est périodique de période xp n – 1xp 0.

La fonction numpy.interp()[modifier | modifier le wikicode]

La fonction np.interp() permet de faire des interpolations linéaires :

y = np.interp(x, xp, yp)

  • xp et yp sont des séquences (vecteurs, listes, n-uplets) de réels, la liste des points décrivant la fonction ;
  • x est la séquence de valeurs auxquelles on s'intéresse.

Par défaut, l'extrapolation se fait par des valeurs constantes, yp 0 et yp n – 1. On peut indiquer que la fonction a une période avec l'argument period. Par exemple :

import numpy as np
import matplotlib.pyplot as plt

# **********************
# * numpy.interp()     *
# * avec option period *
# **********************

xp = np.linspace(0, 2*np.pi, 5) # Cinq points
yp = np.sin(xp)

x = np.linspace(-2*np.pi, 4*np.pi, 25)

y = np.interp(x, xp, yp, period=2*np.pi) # interpolation linéaire

plt.plot(x, y, "y-")
plt.plot(xp, yp, "r+")

Le module scipy.interpolate[modifier | modifier le wikicode]

Quatre procédures d'interpolation polynomiale.

Le module SciPy fournit d'autres méthodes. Pour cela, il faut charger le module :

from scipy import interpolate

La fonction interpolate.interp1d() (interpolation à une dimension, l'avant-dernière lettre est le chiffre 1) crée une fonction d'interpolation. Pour avoir une interpolation linéaire comme précédemment, on écrit :

f = interpolate.interp1d(xp, yp)
y = f(x)

On peut autoriser l'extrapolation linéaire avec l'option fill_value="extrapolate". En absence de ce paramètre, une valeur de x en dehors de l'intervalle des xp génère une erreur.

Cette fonction propose d'autre méthodes d'interpolation : avec l'option kind= :

  • "linear" : méthode par défaut ;
  • "nearest" : renvoie la valeur de yp correspondant au xp le plus proche ;
  • "zero" : fait une interpolation par une spline d'ordre 0 ;
  • "slinear" : " par une spline d'ordre 1 ;
  • "quadratic" : " par une spline d'ordre 2 ;
  • "cubic" : " par une spline d'ordre 3.

Par exemple :

interpolate.interp1d(xp, yp, fill_value="extrapolate", kind="nearest")

La classe interpolate.BarycentricInterpolator correspond à une fonction polynôme passant par un ensemble de points (xp, yp) :

f = interpolate.BarycentricInterpolator(xp, yp)
y = f.__call__(x)

On a le même résultat avec la fonction interpolate.barycentric_interpolate() :

y = interpolate.barycentric_interpolate(xp, yp, x)

La classe interpolate.KroghInterpolator et la fonction interpolate.krogh_interpolate() permettent de fixer les dérivées successives du polynôme : pour définir la dérivée au point i, il suffit de mettre deux fois la valeur xpi dans le vecteur xp, et de mettre ypi, ypi dans le vecteur yp. Si une valeur xpi figure n fois, on définit les dérivées jusqu'à l'ordre n – 1.

f = interpolate.KroghInterpolator(xp, yp)
y = f.__call__(x)

y = interpolate.Krogh_interpolate(xp, yp, x)

On peut avoir la dérivée d'ordre n avec la syntaxe :

f = interpolate.KroghInterpolator(xp, yp)
y_n = f.derivative(x, n) # dérivée n-ième
y_N = f.derivatives(x, N) # dérivées première à n-ième ;
                          # notez le « s » au nom de la fonction

On peut avoir une interpolation par un polynôme cubique par parties PCHIP (piecewise cubic hermite interpolating polynomial) avec la classe interpolate.PchipInterpolator ou la fonction interpolate.pchip_interpolate() :

f = interpolate.PchipInterpolator(xp, yp) # renvoie NaN en cas d'extrapolation
f = interpolate.PchipInterpolator(xp, yp, extrapolate=True) # autorise l'extrapolation
y = f.__call__(x) # on peut ajouter extrapolate=True lors de l'évaluation
y_n = f.__call__(x, n) # valeurs de la dérivée n-ième
f_n = f.derivative(n) # construit la fonction dérivée n-ième
F_n = f.antiderivative(n) # construit la n-ième primitive
racines = f.roots() # racines du polynôme

y = interpolate.pchip_interpolate(xi, yi, x)
y_n = interpolate.pchip_interpolate(xi, yi, x, n) # dérivée n-ième

La classe interpolate.Akima1DInterpolator définit une fonction d'interpolation selon la méthode Akima : elle construit des sous-splines à partir de polynômes cubiques par partie :

f = interpolate.Akima1DInterpolator(xp, yp)
y = f.__call__(x)
f_n = f.derivative(n) # construit la fonction dérivée n-ième
F_n = f.antiderivative(n) # construit la n-ième primitive
racines = f.roots() # racines du polynôme

Voici à titre d'exemple le code source de l'image ci-contre.

Notes et références[modifier | modifier le wikicode]


Statistiques < > Régression et optimisation