Interpolation polynomiale : algorithme de Neville

1. Introduction

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

On pose

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

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 (x_i,y_i) est :

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 (x_i,y_i),(x_i+1,y_i+1) est :

On a en effet :

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

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 :

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

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

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 y_i. Sinon, elle renvoie :

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

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