Python pour le calcul scientifique/Polynômes
Les fonctions spécifiques aux polynômes sont dans le sous-module numpy.polynomial.polynomial. Comme ce préfixe est fastidieux à taper, nous l'abrégeons en nppol. L'en-tête du programme est donc :
#!/usr/bin/python3
import numpy as np
import numpy.polynomial.polynomial as nppol
import matplotlib.pyplot as plt
Par la suite, on peut manipuler les polynômes de deux manières : soit sous la forme d'un vecteur de coefficients, soit sous la forme d'une classe.
Vecteur de coefficients modifier
Le polynôme a0 + a1⋅x + a2⋅x2 est décrit par le vecteur [a0, a1, a2] (ou par le n-uplet (a0, a1, a2)). Pour évaluer le polynôme avec une valeur x donnée, on utilise :
nppol.polyval(x, [a0, a1, a2])
où x peut être une matrice, le résultat étant alors la matrice des résultats. Cela utilise la méthode de Horner.
Pour connaître les racine d'un polynôme :
nppol.polyroots([a0, a1, a2])
À l'inverse, si l'on veut avoir les coefficients du polynôme ayant pour racines r1, r2 et r3, c'est-a-dire le polynôme (x – r1)⋅(x – r2)⋅(x – r3), on utilise :
nppol.polyfromroots([r1, r2, r3])
On peut aussi évaluer un polynôme défini directement par ses racines :
nppol.polyalfromroots(x, [r1, r2, r3])
Le vecteur du affine polynôme a0 + a1⋅x, donc le vecteur [a0, a1], s'obtient avec
nppol.polyline(a0, a1)
NumPy propose des polynômes élémentaires :
nppol.polyzero # p(x) = 0
nppol.polyone # p(x) = 1
nppol.polyx # p(x) = x
On peut effectuer les opérations polynomiales suivantes :
nppol.polyadd(p1, p2) # somme de deux polynômes, p1 + p2
nppol.polysub(p1, p2) # p1 - p2
nppol.polymul(p1, p2) # produit des polynômes p1 × p2
nppol.polymulx(p) # produit de p × ''x''
nppol.polydiv(p1, p2) # division euclidienne p1/p2
nppol.polypol(p, n) # élévation à la puissance p^n
nppol.polyder(p) # dérivation p'
nppol.polyder(p, n) # dérivée n-ième
nppol.polyint(p) # primitive
nppol.polyint(p, 1, k) # primitive avec la constante d'intégration k
nppol.polyint(p, n) # n-ième primitive
nppol.polyint(p, 2, [k1, k2]) #
Régression polynomiale modifier
Si l'on a un nuage de points M(x, y), les variables x et y étant des vecteurs de même taille, alors on peut déterminer le « meilleur polynôme » de degré d décrivant le nuage de points M c'est-à-dire la fonction polynôme ayant le plus petit écart quadratique. Cette régression polynomiale se fait de la manière suivante :
nppol.polyfit(x, y, d)
Par exemple :
x = np.arange(10)
y = x*x + np.random.randn(10) # valeurs perturbées
p = nppol.polyfit(x, y, 2) # régression polynomiale de degré 2
print(p) # devrait être proche de [0, 0, 1]
plt.plot(x, y, "+")
plt.plot(x, nppol.polyval(x, p))
Si l'on veut forcer certains coefficients à zéro, on n'indique non pas le degré du polynôme mais la liste des degrés qu'il contient. Par exemple, pour imposer a0 = 0 :
p = nppol.polyfit(x, y, [1, 2]) # le degré 0 est absent de la liste des degrés
Classe Polynomial
modifier
Le module propose la classe Polynomial. Le polynôme a donc une représentation propre et dispose de méthodes optimisées.
Les exemples ci-dessus deviennent :
p = nppol.Polynomial([a0, a1, a2]) # ou nppol.Polynomial((a0, a1, a2)) ; p(x) = a0 + a1⋅x + a2⋅x^2
p.__call(x)__ # évaluation de p(x), x étant un nombre ou un vecteur
p.roots() # racines de p
p.degree() # degré de p
p = nppol.Polynomial.fromroots([r0, r1, r2]) # ou .fromroots((r0, r1, r2)) ; p(x) = (x - r1)⋅(x - r2)⋅(x - r3)
p = nppol.Polynomial.fit(x, y, 2) # régression polynomiale de degré 2 sur (x, y)
p = nppol.Polynomial.fit(x, y, [1, 2]) # idem mais a0 = 0 (puisque le degré 0 est absent du vecteur des degrés)
Les opérations classiques s'appliquent aux polynômes : +, -, *, //, %, divmod(), **. La dérivation et l'intégration s'obtiennent avec :
p.deriv() # dérivée
p.deriv(n) # n-ième dérivée
p.integ() # primitive
p.integ(1, k) # primitive avec la constante d'intégration k
p.integ(n) # n-ième primitive
p.integ(2, [k1, k2]) # 2e primitive avec les constantes d'intégration k1 et k2
Pour obtenir le vecteur des coefficients, on peut utiliser :
p.convert().coef
Pour reprendre l'exemple de régression polynomiale ci-dessus :
x = np.arange(10)
y = x*x + np.random.randn(10) # valeurs perturbées
p = nppol.Polynomial.fit(x, y, 2) # régression polynomiale de degré 2
print(p.convert().coef) # devrait être proche de [0, 0, 1]
plt.plot(x, y, "+")
plt.plot(x, p.__call__(x))