Table des matières PDFPython

Interpolation polynomiale : algorithme de Neville

1. Introduction

On considère une fonction f(x) dont on connaît les valeurs aux N points

x0,x1,...xN-1

On pose

yi=f(xi)(1)

Pour calculer f(x) pour x quelconque, on peut effectuer une interpolation polynomiale, qui consiste à calculer f(x) en utilisant l'unique polynôme de degré N-1 qui passe par les points

(x0,y0),(x1,y1),...(xN-1,yN-1)(2)

On présente ici l'algorithme de Neville, qui permet de faire efficacement une interpolation polynomiale lorsque le nombre de valeurs de x est faible.

2. Algorithme de Neville

L'algorithme de Neville ([1]) construit les polynômes par récurrence. On commence par les polynômes de degré 0. Le polynôme de degré 0 qui passe par le point (xi,yi) est :

Pi(x)=yi(3)

Partant de ces N polynômes de degré 0, on construit les polynômes de degré 1. Le polynôme de degré 1 passant par les points (xi,yi),(xi+1,yi+1) est :

Pi,i+1(x)=(x-xi+1)Pi(x)+(xi-x)Pi+1(x)xi-xi+1(4)

On a en effet :

Pi,i+1(xi)=Pi(xi)=yi(5)Pi,i+1(xi+1)=Pi+1(xi+1)=yi+1(6)

La relation de récurrence permettant de passer de deux polynômes de degré m-1 à un polynôme de degré m est :

Pi,i+1,...i+m(x)=(x-xi+m)Pi,i+1,...,i+m-1(x)+(xi-x)Pi+1,i+2,...i+m(x)xi-xi+m(7)

Pour le démontrer, on raisonne par récurrence, en supposant que les deux polynômes de degré m-1 apparaissant dans la relation ci-dessus passent par les m points figurant en indice. On vérifie alors que le polynôme de degré m passe par les m+1 points d'indices i,i+1,...i+m. Pour le premier point, on a en effet :

Pi,i+1,...i+m(xi)=(xi-xi+m)Pi,i+1,...,i+m-1(xi)xi-xi+m=yi(8)

La vérification est similaire pour le point xi+m. Pour un point xi+k (commun aux deux polynômes de degré m-1), on a :

Pi,i+1,...i+m(xi+k)=(xi+k-xi+m)yi+k+(xi-xi+k)yi+kxi-xi+m=yi+k(9)

Le tableau suivant montre l'enchaînement des calculs jusqu'au polynôme de degré N-1, dans le cas N=4 :

P0 P01 P1 P012 P12 P0123P2 P123 P23 P3 (10)

3. Implémentation récursive

L'algorithme de Neville est aisément programmé de manière récursive. Pour cela, on définit une fonction r(i,m,x) qui calcule la valeur du polynôme de degré m en fonction des valeurs des deux polynômes de degré m-1

Si m=0, la fonction renvoie yi. Sinon, elle renvoie :

(x-xi+m)r(i,m-1,x)+(xi-x)r(i+1,m-1,x)xi-xi+m(11)

On définit une classe dont le constructeur prend en argument les valeurs de x et y sous forme de listes.

neville.py
class Neville:
    def __init__(self,x,y):
        self.x = x
        self.y = y
        self.N = len(x)
             

La fonction suivante définit la récurrence, en s'appelant récursivement deux fois :

    def recurrence(self,i,m,x):
        if m==0:
            return self.y[i]
        else:
            return ((x-self.x[i+m])*self.recurrence(i,m-1,x)+\
            (self.x[i]-x)*self.recurrence(i+1,m-1,x))/(self.x[i]-self.x[i+m])
             

Voici la fonction qui effectue l'interpolation :

    def interpoler(self,x):
        return self.recurrence(0,self.N-1,x)   
             

Pour tester l'algorithme, la fonction suivante effectue l'évaluation d'un polynôme dont les coefficients sont donnés dans une liste c (rangés par degré croissant) :

    def eval_polynome(self,c,x):
        p = len(c)
        i = p-1
        y = c[i]
        while i>0:
            i -= 1
            y = y*x+c[i]
        return y
              

La fonction suivante calcule les valeurs des listes x et y à partir d'un polynôme donné, sur un intervalle donné :

    def polynome(self,c,xmin,xmax):
        p = len(c)
        dx = (xmax-xmin)/p
        for k in range(p):
            self.x[k] = xmin+k*dx
            self.y[k] = self.eval_polynome(c,self.x[k])
              

Faisons un test avec un polynôme de degré 4 :

import neville
x=[0,0,0,0]
y=[0,0,0,0]
nev = neville.Neville(x,y)
c = [1,2,-3,2]
nev.polynome(c,0.0,10.0)
x0 = 3.564
y0 = nev.interpoler(x0)
              
print(y0)
--> 60.562252287999996
print(y0-nev.eval_polynome(c,x0))
--> -1.4210854715202004e-14
Références
[1]  W.H. Press, S.A. Teukolsky, W.T. Vetterling, B.P. Flannery,  Numerical recipes, the art of scientific computing,  (Cambridge University Press, 2007)
Creative Commons LicenseTextes et figures sont mis à disposition sous contrat Creative Commons.