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

Le polynôme a₀ + a₁⋅x + a₂⋅x² 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 r₁, r₂ et r₃, c'est-a-dire le polynôme (x – r₁)⋅(x – r₂)⋅(x – r₃), 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 a₀ + a₁⋅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

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 a₀ = 0 :

p = nppol.polyfit(x, y, [1, 2]) # le degré 0 est absent de la liste des degrés

Classe Polynomial

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))

Notes et références

---

Manipulation de matrices < ↑ > Statistiques