I. Introduction
Dans un précédent billet nous avons réalisé une classe en Python qui nous a permis d'effectuer des opérations de base entre polynômes.
On se propose maintenant d'utiliser cette classe pour développer et réduire un polynôme d'interpolation de Newton.
Après avoir défini l'expression générale du polynôme d'interpolation, on écrira des fonctions en Python permettant d'obtenir sa forme développée et réduite (forme usuelle des polynômes).
II. Polynôme d'interpolation de Newton
II-A. Définition
Le polynôme d'interpolation de Newton N(x) associé à k+1 points (x0, y0), …,(xk, yk) est défini par :
N(x) = [y0] + [y0, y1](x - x0) + [y0, y1, y2](x - x0)(x - x1) + ... + [y0, y1, ..., yk](x - x0)(x - x1)...(x - xk-1)
Où [y0], [y0, y1], [y0, y1, y2], ..., [y0, y1, ..., yk] représentent les différences divisées à l'ordre 0, 1, 2, ..., k.
Par exemple, les différences divisées d'ordre 0, 1 et 2 sont définies comme suit :
631183
On obtient ainsi le schéma récursif suivant :
631184
III. Développement d'un polynôme d'interpolation de Newton
III-A. Classe Polynome
Elle permet de créer des objets Polynome et de redéfinir les opérateurs d'addition et de multiplication pour ces nouveaux objets.
631204
Module de la Classe Polynome :
class Polynome:
def __init__(self, liste_coefs=): # méthode constructeur de la classe
self.coefs = liste_coefs # on définit la liste des coefficients du polynôme [a0, a1, ..., an]
self.reduire() # suppression si nécessaire des zéros en queue de liste de coefficients. Exemple : [2, 3, 1, 0, 0] -> [2, 3, 1]
def __str__(self): # permet d'afficher le polynôme sous la forme 1 + 2x + 3x^2
s="" # initialisation de la chaîne de caractères
# on vérifie d’abord si le degré du polynôme est nul
if (len(self.coefs)-1==0):
return str(self.coefs)
else: # sinon
if self.coefs!=0:
s=str(self.coefs) + " + "
for i in range(1, len(self.coefs)): # parcours des indices des coefficients du polynôme : [a1, a2, ..., an]
if self.coefs!=0: # si le coefficient de degré i n'est pas nul
if self.coefs!=1: # si le coefficient de degré i est différent de 1
s+="{}.X^{} + ".format(self.coefs,i)
else: s+="X^{} + ".format(i)
return s[:-3] # on retourne l'expression du polynôme
def eval_degre(): # retourne le degré du polynôme
return (len(self.coefs)-1)
def __add__(self, other): # méthode permettant de redéfinir l'opérateur « + » pour 2 polynômes : (1 + 2x + x^2) + (1 + x) = 2 + 3x + x^2
# p1 = self, p2 = other
if len(other.coefs) >len(self.coefs): # si degré de p2 > degré de p1
liste_coefs = other.coefs[:]; n = len(self.coefs) # on copie les coefs du polynôme de degré le plus élevé et la longueur de la liste de coefs la plus petite.
else: liste_coefs = self.coefs[:]; n = len(other.coefs) # sinon, ...
for i in range(n): # parcours des indices de liste_coefs
liste_coefs = self.coefs + other.coefs # addition des coefficients de degré i
return Polynome(liste_coefs) # renvoie le polynôme résultat de l'addition
def reduire(self):
# tant que le dernier élément de la liste est nul
while self.coefs[-1] == 0 and len(self.coefs)>1:
self.coefs.pop() # supprimer le dernier élément
def __mul__(self, other): # méthode permettant de redéfinir l'opérateur « * » pour 2 polynômes : (1 + x) * (1 + 2x) = 1 + 3x + 2x^2
# initialisation de la liste des coefficients qu'avec des zéros
liste_coefs=*(len(self.coefs)+len(other.coefs)-1) # exemple : [0, 0, 0]
for i1 in range(len(self.coefs)): # parcours des indices des coefs du polynôme n°1
for i2 in range(len(other.coefs)): # parcours des indices des coefs du polynôme n°2
# multiplication des coefficients d'indices i1 et i2
liste_coefs[i1+i2] = liste_coefs[i1+i2] + self.coefs*other.coefs
poly=Polynome(liste_coefs) # création de l'objet Polynome basé sur la liste
return poly # renvoie le polynôme résultat de la multiplication
def __pow__(self, n): # méthode permettant de redéfinir l'opérateur de puissance : self ** n
# on crée l'objet poly à partir d'une liste contenant un seul élément 1, élément neutre pour la multiplication des polynômes
poly = Polynome(,self.domaine)
for i in range(n): # on multiplie n fois poly par self à l'aide de l'opérateur *
poly = poly*self # équivalent à : poly = poly.__mul__(self)
return poly # renvoie le polynôme résultat de l'opération (self ** n)
def __eq__(poly1, other): # méthode permettant de redéfinir l'opérateur « == » pour 2 polynômes
return (poly1.coefs==other.coefs) # renvoie True si les 2 liste de coefficients sont égales
def eval(self,x): # méthode permettant d'évaluer le polynôme en fonction de x
# intialisation des variables
valeur_polynome = self.coefs # valeur_polynome = a0
power=1
for coef in self.coefs[1:]: # parcours des coefficients du polynôme à partir de a1 : a1, a2, ..., an
power = power*x # calcul de la puissance de x pour ai : power = x^i
valeur_polynome += coef*power # valeur_polynome = valeur_polynome + ai*x^i
return valeur_polynome # renvoie la valeur du polynôme
def eval_horner(self,x): # méthode permettant d'évaluer le polynôme en fonction de x
# intialisation de la variable
valeur_polynome = self.coefs[-1] # valeur_polynome = an
for coef in reversed(self.coefs[:-1]): # parcours des coefficients du polynôme à partir de an-1 : an-1, ..., a1, a0
valeur_polynome = valeur_polynome*x + coef # valeur_polynome = valeur_polynome*x + ai
return valeur_polynome # renvoie la valeur du polynôme
Note importante : on va utiliser par la suite les opérateurs « + » et « * » entre polynômes pour développer et réduire progressivement le polynôme de Newton.
III-B. Evaluation des différences divisées
Pour éviter de recalculer à chaque fois les différences divisées successives [y0], [y0, y1], [y0, y1, y2], ..., on doit d'abord les mémoriser à l'aide d'une fonction Python ayant comme arguments les 2 listes de valeurs x et y :
def differences_divisees(x,y):
n=len(x)-1 # nombre de valeurs de la série - 1, il correspond au degré du polynôme, exemple pour n=2 : a2x^2 + a1x + a0, où a2, a1 et a0 représentent les coefficients
dy = y.copy() # on copie la liste y dans dy : différences divisées à l'ordre 0
diff_divisees = [dy] # on initialise la liste permettant d'enregistrer les différences divisées :
for k in range(1,n+1): # on parcourt les ordres suivants des différences divisées : ordre 1, 2, ..., n
for i in range(n-k+1): # on parcourt les indices de x et dy jusqu'à n-k
h = (x[i+k] - x) # on détermine la différence entre xi+k et xi pour le calcul de la différence divisée à l'ordre k
dyi = (dy[i+1] - dy)/h # on évalue la différence divisée à l'ordre k : [yi, ..., yi+k])
dy = dyi # on remplace l'élément d'indice i de la liste dy par la valeur de la différence divisée dyi
diff_divisees = diff_divisees + [dy] # on ajoute à la liste la différence divisée d'ordre k.
return diff_divisees # on renvoie la liste des différences divisées
Déterminons maintenant le polynôme développé et réduit à partir des 2 listes de valeurs.
III-C. Développement du polynôme suivant un schéma classique
Nous avons vu précédemment que le polynôme d'interpolation de Newton N(x) associé à k+1 points (x0, y0), …,(xk, yk) est défini par :
N(x) = [y0] + [y0, y1](x - x0) + [y0, y1, y2](x - x0)(x - x1) + ... + [y0, y1, ..., yk](x - x0)(x - x1)...(x - xk-1)
On se propose d'écrire une fonction qui permette de développer et de réduire cette expression de façon à obtenir un polynôme exprimé sous la forme usuelle :
N(x) = a0 + a1x + a2x2 + ... + akxk
Où a0, a1, ..., ak désignent les coefficients réels du polynôme.
La fonction developper_polynome_newton prend en arguments 2 listes x et y, représentant la série de points (xi,yi), et renvoie l'objet Polynome recherché avec ses coefficients :
def developper_polynome_newton(x, y):
diff_divisees= differences_divisees(x,y) # on génère la liste des différences divisées : [, [y0, y1], [y0, y1, y2], ...]
k=len(x)-1 # indice maxi de la liste des x
y0 = diff_divisees # évaluation de la différence divisée
polynome_developpe = Polynome() # on initialise le polynôme avec -> P(x) = y0
produit_facteurs = Polynome() # on initialise le produit de facteurs avec -> p = 1
for i in range(k): # parcours des indices : 0, 1, ..., k-1
facteur_x = Polynome([-x,1]) # on crée le polynôme correspondant au facteur (x - xi)
diff_div = Polynome([diff_divisees[i+1]]) # polynôme associé à la différence divisée : p = diff_divisees[i+1]
produit_facteurs = produit_facteurs*facteur_x # on génère le polynôme résultat du développement du produit de facteurs : (x -x0)(x-x1)...
polynome_developpe = polynome_developpe + diff_div*produit_facteurs # on développe et réduit progressivement le polynôme de Newton
return polynome_developpe # renvoie le polynôme développé et réduit
Note importante : dans ce cas nous avons besoin au mieux de réaliser 2k-1 multiplications pour développer notre polynôme.
III-D. Développement du polynôme suivant un schéma de Horner
Le polynôme d'interpolation de Newton peut également s'écrire suivant un schéma de Horner :
N(x) = ((...([y0, y1, ..., yk](x - xk-1) + [y0, y1, ..., yk-1])(x - xk-2) + ... + [y0, y1, y2])(x - x1) + [y0, y1])(x - x0) + [y0]
On peut identifier la relation de récurrence valable pour i compris entre 0 et k-1 :
Pi+1(x) = Pi(x)(x - xk-i-1) + [y0, y1, ..., yk-i-1]
Avec :
P0(x) = [y0, y1, ..., yk]
Et :
Pk(x) = N(x)
Nous pouvons ainsi développer notre polynôme avec seulement k multiplications.
La fonction developper_polynome_newton_horner basée sur ce schéma prend en arguments 2 listes x et y, représentant la série de points (xi,yi), et renvoie le polynôme recherché avec sa liste de coefficients :
def developper_polynome_newton_horner(x, y):
diff_divisees = differences_divisees(x,y) # on génère la liste des différences divisées : [, [y0, y1], [y0, y1, y2], ...]
k=len(x)-1 # indice maxi de la liste des x
polynome_developpe = Polynome([diff_divisees]) # on initialise le polynôme avec [y0, y1, ..., yk] -> P(x) = [y0, y1, ..., yk]
for i in range(k-1,-1,-1): # parcours des indices : k-1, k-2, ..., 1, 0
facteur_x = Polynome([-x,1]) # on crée le polynôme correspondant au facteur (x - xi)
diff_div = Polynome([diff_divisees]) # polynôme associé à la différence divisée : p = diff_divisees
polynome_developpe = polynome_developpe*facteur_x + diff_div # on développe et réduit progressivement le polynôme de Newton
return polynome_developpe # renvoie le polynôme développé et réduit
III-E. Test des fonctions
Module complet de test contenant la classe Polynome avec :
x=[1, 2, 4]
y=[6, 17, 57]
class Polynome:
def __init__(self, liste_coefs=): # méthode constructeur de la classe
self.coefs = liste_coefs # on définit la liste des coefficients du polynôme [a0, a1, ..., an]
self.reduire() # suppression si nécessaire des zéros en queue de liste de coefficients. Exemple : [2, 3, 1, 0, 0] -> [2, 3, 1]
def __str__(self): # permet d'afficher le polynôme sous la forme 1 + 2x + 3x^2
s="" # initialisation de la chaîne de caractères
# on vérifie d’abord si le degré du polynôme est nul
if (len(self.coefs)-1==0):
return str(self.coefs)
else: # sinon
if self.coefs!=0:
s=str(self.coefs) + " + "
for i in range(1, len(self.coefs)): # parcours des indices des coefficients du polynôme : [a1, a2, ..., an]
if self.coefs!=0: # si le coefficient de degré i n'est pas nul
if self.coefs!=1: # si le coefficient de degré i est différent de 1
s+="{}.X^{} + ".format(self.coefs,i)
else: s+="X^{} + ".format(i)
return s[:-3] # on retourne l'expression du polynôme
def eval_degre(self): # retourne le degré du polynôme
return (len(self.coefs)-1)
def __add__(self, other): # méthode permettant de redéfinir l'opérateur « + » pour 2 polynômes : (1 + 2x + x^2) + (1 + x) = 2 + 3x + x^2
# p1 = self, p2 = other
if len(other.coefs) >len(self.coefs): # si degré de p2 > degré de p1
liste_coefs = other.coefs[:]; n = len(self.coefs) # on copie les coefs du polynôme de degré le plus élevé et la longueur de la liste de coefs la plus petite.
else: liste_coefs = self.coefs[:]; n = len(other.coefs) # sinon, ...
for i in range(n): # parcours des indices de liste_coefs
liste_coefs = self.coefs + other.coefs # addition des coefficients de degré i
return Polynome(liste_coefs) # renvoie le polynôme résultat de l'addition
def reduire(self):
# tant que le dernier élément de la liste est nul
while self.coefs[-1] == 0 and len(self.coefs)>1:
self.coefs.pop() # supprimer le dernier élément
def __mul__(self, other): # méthode permettant de redéfinir l'opérateur « * » pour 2 polynômes : (1 + x) * (1 + 2x) = 1 + 3x + 2x^2
# initialisation de la liste des coefficients qu'avec des zéros
liste_coefs=*(len(self.coefs)+len(other.coefs)-1) # exemple : [0, 0, 0]
for i1 in range(len(self.coefs)): # parcours des indices des coefs du polynôme n°1
for i2 in range(len(other.coefs)): # parcours des indices des coefs du polynôme n°2
# multiplication des coefficients d'indices i1 et i2
liste_coefs[i1+i2] = liste_coefs[i1+i2] + self.coefs*other.coefs
poly=Polynome(liste_coefs) # création de l'objet Polynome basé sur la liste
return poly # renvoie le polynôme résultat de la multiplication
def __pow__(self, n): # méthode permettant de redéfinir l'opérateur de puissance : self ** n
# on crée l'objet poly à partir d'une liste contenant un seul élément 1, élément neutre pour la multiplication des polynômes
poly = Polynome(,self.domaine)
for i in range(n): # on multiplie n fois poly par self à l'aide de l'opérateur *
poly = poly*self # équivalent à : poly = poly.__mul__(self)
return poly # renvoie le polynôme résultat de l'opération (self ** n)
def __eq__(poly1, other): # méthode permettant de redéfinir l'opérateur « == » pour 2 polynômes
return (poly1.coefs==other.coefs) # renvoie True si les 2 liste de coefficients sont égales
def eval(self,x): # méthode permettant d'évaluer le polynôme en fonction de x
# intialisation des variables
valeur_polynome = self.coefs # valeur_polynome = a0
power=1
for coef in self.coefs[1:]: # parcours des coefficients du polynôme à partir de a1 : a1, a2, ..., an
power = power*x # calcul de la puissance de x pour ai : power = x^i
valeur_polynome += coef*power # valeur_polynome = valeur_polynome + ai*x^i
return valeur_polynome # renvoie la valeur du polynôme
def eval_horner(self,x): # méthode permettant d'évaluer le polynôme en fonction de x
# intialisation de la variable
valeur_polynome = self.coefs[-1] # valeur_polynome = an
for coef in reversed(self.coefs[:-1]): # parcours des coefficients du polynôme à partir de an-1 : an-1, ..., a1, a0
valeur_polynome = valeur_polynome*x + coef # valeur_polynome = valeur_polynome*x + ai
return valeur_polynome # renvoie la valeur du polynôme
def differences_divisees(x,y):
n=len(x)-1 # nombre de valeurs de la série - 1, il correspond au degré du polynôme, exemple pour n=2 : a2x^2 + a1x + a0, où a2, a1 et a0 représentent les coefficients
dy = y.copy() # on copie la liste y dans dy : différences divisées à l'ordre 0
diff_divisees = [dy] # on initialise la liste permettant d'enregistrer les différences divisées avec
for k in range(1,n+1): # on parcourt les ordres suivants des différences divisées : ordre 1, 2, ..., n
for i in range(n-k+1): # on parcourt les indices de x et dy jusqu'à n-k
h = (x[i+k] - x) # on détermine la différence entre xi+k et xi pour le calcul de la différence divisée à l'ordre k
dyi = (dy[i+1] - dy)/h # on évalue la différence divisée à l'ordre k : [yi, ..., yi+k])
dy = dyi # on remplace l'élément d'indice i de la liste dy par la valeur de la différence divisée dyi
diff_divisees = diff_divisees + [dy] # on ajoute à la liste la différence divisée d'ordre k.
return diff_divisees # on renvoie la liste des différences divisées
def expression_polynome_newton(x, y):
diff_divisees= differences_divisees(x,y) # on génère la liste des différences divisées : [, [y0, y1], [y0, y1, y2], ...]
k=len(x)-1 # indice maxi de la liste des x
y0 = diff_divisees # évaluation de la différence divisée
polynome_newton = str(y0) # on initialise le polynôme avec -> p(x) = y0
produit_facteurs = '' # on initialise le produit de facteur avec une chaîne vide
for i in range(k): # parcours des indices : 0, 1, ..., k-1
if x<0:
facteur_x = "(x + " + str(-x) + ")" # facteur de x : (x - xi)
else : facteur_x = "(x - " + str(x) + ")"
diff_div = str(diff_divisees[i+1]) # différence divisée à l'ordre i+1 : diff_divisees[i+1]
produit_facteurs = produit_facteurs + facteur_x # on génère le polynôme associé au produit de facteurs : (x -x0)(x-x1)...
polynome_newton = polynome_newton + " + " + diff_div + str(".") + produit_facteurs # on génère progressivement l'expression du polynôme de Newton
return polynome_newton # renvoie l'expression du polynôme de Newton
def eval_polynome_newton(x, y, valeur_x):
diff_divisees = differences_divisees(x,y) # on génère la liste des différences divisées : [, [y0, y1], [y0, y1, y2], ...]
k=len(x)-1 # indice maxi de la liste des x
valeur_polynome_newton = diff_divisees # on initialise la variable : valeur_polynome_newton = [y0, y1, ..., yk]
for i in range(k-1,-1,-1): # parcours des indices : k-1, k-2, ..., 1, 0
facteur_x = (valeur_x -x) # on évalue le facteur (x - xi)
valeur_polynome_newton = valeur_polynome_newton*facteur_x + diff_divisees # on évalue progressivement le polynôme de Newton
return valeur_polynome_newton # renvoie la valeur du polynôme
def developper_polynome_newton(x, y):
diff_divisees= differences_divisees(x,y) # on génère la liste des différences divisées : [, [y0, y1], [y0, y1, y2], ...]
k=len(x)-1 # indice maxi de la liste des x
y0 = diff_divisees # évaluation de la différence divisée
polynome_developpe = Polynome() # on initialise le polynôme avec -> P(x) = y0
produit_facteurs = Polynome() # on initialise le produit de facteurs avec -> p = 1
for i in range(k): # parcours des indices : 0, 1, ..., k-1
facteur_x = Polynome([-x,1]) # on crée le polynôme correspondant au facteur (x - xi)
diff_div = Polynome([diff_divisees[i+1]]) # polynôme associé à la différence divisée : p = diff_divisees[i+1]
produit_facteurs = produit_facteurs*facteur_x # on génère le polynôme résultat du développement du produit de facteurs : (x -x0)(x-x1)...
polynome_developpe = polynome_developpe + diff_div*produit_facteurs # on développe et réduit progressivement le polynôme de Newton
return polynome_developpe # renvoie le polynôme développé et réduit
def developper_polynome_newton_horner(x, y):
diff_divisees = differences_divisees(x,y) # on génère la liste des différences divisées : [, [y0, y1], [y0, y1, y2], ...]
k=len(x)-1 # indice maxi de la liste des x
polynome_developpe = Polynome([diff_divisees]) # on initialise le polynôme avec [y0, y1, ..., yk] -> P(x) = [y0, y1, ..., yk]
for i in range(k-1,-1,-1): # parcours des indices : k-1, k-2, ..., 1, 0
facteur_x = Polynome([-x,1]) # on crée le polynôme correspondant au facteur (x - xi)
diff_div = Polynome([diff_divisees]) # polynôme associé à la différence divisée : p = diff_divisees
polynome_developpe = polynome_developpe*facteur_x + diff_div # on développe et réduit progressivement le polynôme de Newton
return polynome_developpe # renvoie le polynôme développé et réduit
# Initialisation des 2 listes de valeurs x et y
x=[1, 2, 4]
y=[6, 17, 57]
print("Forme classique du polynôme de Newton :")
# Affiche l'expression classique du polynôme de Newton
print('N(x) = ' + expression_polynome_newton(x, y))
# Affiche la valeur de N(3)
print('N(3) = ' + str(eval_polynome_newton(x, y, 3)))
print("")
# Détermination du polynôme développé et réduit
polynome_developpe = developper_polynome_newton(x, y)
print("Forme développée et réduite du polynôme :")
# Affiche l'expression du polynôme sous sa forme développée et réduite
print('N(x) = ' + polynome_developpe.__str__())
# Affiche la valeur de N(3)
print('N(3) = ' + str(polynome_developpe.eval(3)))
print("")
# Détermination du polynôme développé et réduit (méthode de Horner)
polynome_developpe = developper_polynome_newton_horner(x, y)
print("Forme développée et réduite du polynôme - méthode de Horner :")
# Affiche l'expression du polynôme sous sa forme développée et réduite
print('N(x) = ' + polynome_developpe.__str__())
# Affiche la valeur de N(3)
print('N(3) = ' + str(polynome_developpe.eval_horner(3)))
Le code affiche :
Forme classique du polynôme de Newton :
N(x) = 6 + 11.0.(x - 1) + 3.0.(x - 1)(x - 2)
N(3) = 34.0
Forme développée et réduite du polynôme :
N(x) = 1.0 + 2.0.X^1 + 3.0.X^2
N(3) = 34.0
Forme développée et réduite du polynôme - méthode de Horner :
N(x) = 1.0 + 2.0.X^1 + 3.0.X^2
N(3) = 34.0
IV. Conclusion
Après avoir donné l'expression générale du polynôme d'interpolation de Newton, nous avons pu utiliser notre classe Polynome pour obtenir sa forme développée et réduite.
Chacun pourra ensuite librement développer d'autres types de polynômes à l'aide de cette classe.
Sources :
https://fr.wikipedia.org/wiki/Interpolation_newtonienne
https://en.wikipedia.org/wiki/Newton_polynomial
https://fr.wikipedia.org/wiki/Diff%C3%A9rences_divis%C3%A9es
Soutenez le club developpez.com en souscrivant un abonnement pour que nous puissions continuer à vous proposer des publications.