Mathématiques et Python - Développer et réduire un polynôme d'interpolation de Newton à l'aide de la classe Polynome

Un billet blog de Denis Hulo

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
class Polynome: 
  
    def __init__(self, liste_coefs=[0]): # 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[0]) 
        else: # sinon 
            if self.coefs[0]!=0: 
                s=str(self.coefs[0]) + " + " 
            for i in range(1, len(self.coefs)): # parcours des indices des coefficients du polynôme : [a1, a2, ..., an] 
                if self.coefs[i]!=0: # si le coefficient de degré i n'est pas nul 
                    if self.coefs[i]!=1: # si le coefficient de degré i est différent de 1 
                        s+="{}.X^{} + ".format(self.coefs[i],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[i] = self.coefs[i] + other.coefs[i] # 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=[0]*(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[i1]*other.coefs[i2] 
  
        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([1],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[0] # 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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
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[0]] # on initialise la liste permettant d'enregistrer les différences divisées : [y0] 
  
    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[i]) # 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[i])/h # on évalue la différence divisée à l'ordre k : [yi, ..., yi+k]) 
  
            dy[i] = 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[0]] # on ajoute à la liste la différence divisée d'ordre k. 
  
    return diff_divisees # on renvoie la liste des différences divisées

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
def developper_polynome_newton(x, y): 
  
    diff_divisees= differences_divisees(x,y) # on génère la liste des différences divisées : [[y0], [y0, y1], [y0, y1, y2], ...] 
    k=len(x)-1 # indice maxi de la liste des x 
    y0 = diff_divisees[0] # évaluation de la différence divisée [y0] 
  
    polynome_developpe = Polynome([y0]) # on initialise le polynôme avec [y0] -> P(x) = y0 
    produit_facteurs = Polynome([1]) # on initialise le produit de facteurs avec [1] -> p = 1  
  
    for i in range(k): # parcours des indices : 0, 1, ..., k-1 
        facteur_x = Polynome([-x[i],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

1
2
3
4
5
6
7
8
9
10
11
12
13
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], [y0, y1], [y0, y1, y2], ...] 
    k=len(x)-1 # indice maxi de la liste des x 
  
    polynome_developpe = Polynome([diff_divisees[k]]) # 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[i],1]) # on crée le polynôme correspondant au facteur (x - xi) 
        diff_div = Polynome([diff_divisees[i]]) # polynôme associé à la différence divisée : p = diff_divisees[i] 
        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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
class Polynome: 
  
    def __init__(self, liste_coefs=[0]): # 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[0]) 
        else: # sinon 
            if self.coefs[0]!=0: 
                s=str(self.coefs[0]) + " + " 
            for i in range(1, len(self.coefs)): # parcours des indices des coefficients du polynôme : [a1, a2, ..., an] 
                if self.coefs[i]!=0: # si le coefficient de degré i n'est pas nul 
                    if self.coefs[i]!=1: # si le coefficient de degré i est différent de 1 
                        s+="{}.X^{} + ".format(self.coefs[i],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[i] = self.coefs[i] + other.coefs[i] # 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=[0]*(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[i1]*other.coefs[i2] 
  
        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([1],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[0] # 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[0]] # on initialise la liste permettant d'enregistrer les différences divisées avec [y0] 
  
    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[i]) # 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[i])/h # on évalue la différence divisée à l'ordre k : [yi, ..., yi+k]) 
  
            dy[i] = 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[0]] # 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], [y0, y1], [y0, y1, y2], ...] 
    k=len(x)-1 # indice maxi de la liste des x 
    y0 = diff_divisees[0] # évaluation de la différence divisée [y0] 
    polynome_newton = str(y0) # on initialise le polynôme avec [y0] -> 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[i]<0: 
            facteur_x = "(x + " + str(-x[i]) + ")" # facteur de x : (x - xi) 
        else : facteur_x = "(x - " + str(x[i]) + ")" 
        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], [y0, y1], [y0, y1, y2], ...] 
    k=len(x)-1 # indice maxi de la liste des x 
  
    valeur_polynome_newton = diff_divisees[k] # 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[i]) # on évalue le facteur (x - xi) 
        valeur_polynome_newton = valeur_polynome_newton*facteur_x + diff_divisees[i] # 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], [y0, y1], [y0, y1, y2], ...] 
    k=len(x)-1 # indice maxi de la liste des x 
    y0 = diff_divisees[0] # évaluation de la différence divisée [y0] 
  
    polynome_developpe = Polynome([y0]) # on initialise le polynôme avec [y0] -> P(x) = y0 
    produit_facteurs = Polynome([1]) # on initialise le produit de facteurs avec [1] -> p = 1  
  
    for i in range(k): # parcours des indices : 0, 1, ..., k-1 
        facteur_x = Polynome([-x[i],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], [y0, y1], [y0, y1, y2], ...] 
    k=len(x)-1 # indice maxi de la liste des x 
  
    polynome_developpe = Polynome([diff_divisees[k]]) # 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[i],1]) # on crée le polynôme correspondant au facteur (x - xi) 
        diff_div = Polynome([diff_divisees[i]]) # polynôme associé à la différence divisée : p = diff_divisees[i] 
        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)))