I. Introduction
En analyse numérique, la méthode de Newton-Cotes permet d'évaluer une intégrale sur un intervalle réel [a, b], ceci à l’aide d’une interpolation polynomiale de la fonction en des points répartis uniformément.
638755
On propose dans ce billet une implémentation en Python de cette méthode à l'aide de la classe Polynome.
Note : les définitions données dans ce billet sont toutes issues des pages wikipedia Formule de Newton-Cotes et Calcul numérique d'une intégrale.
II. Description de la méthode
La fonction f est évaluée en des points équidistants xi = a + ih, pour i = 0, … , n et h = (b – a)/n. La formule d'intégration de degré n est définie ainsi :
638693
où les wi sont appelés les coefficients de quadrature. Ils se déduisent d'une base de polynômes de Lagrange et sont indépendants de la fonction f.
Plus précisément, si L(x) est l'interpolation lagrangienne aux points (xi, f(xi)) et li(x) le polynôme de Lagrange tel que :
638694
Alors :
638695
Ainsi :
638696
III. Implémentation de la méthode de Newton-Cotes
On doit d'abord développer, puis intégrer sur un intervalle [a, b] chaque polynôme de Lagrange li(x) pour obtenir les coefficients de quadrature wi.
Enfin, on obtient la valeur de l'intégrale de f sur l'intervalle [a, b] en utilisant la formule :
638693
III-A. Développement et réduction des polynômes de Lagrange à l'aide de la classe Polynome
Comme on l'a déjà vu la classe Polynome permet de représenter un polynôme de la forme :
p(x) = 2x2 + x + 1 = 1 + x + 2x2
à l'aide de sa liste de coefficients [1, 1, 2] :
p = Polynome([1, 1, 2])
On peut également réaliser des opérations de base entre ces objets Polynome, comme la multiplication entre deux polynômes :
p1 = Polynome([1, 1, 2])
p2 = Polynome([2, 1, 1])
p = p1*p2
On obtient comme résultat un objet p représentant un nouveau polynôme exprimé sous sa forme développée et réduite.
Ainsi, la multiplication :
pi = pi*(x - xj)/(xi - xj)
Equivaut en utilisant des objets Polynome à :
dj = (x - x)
pi = pi * Polynome([-x/dj, 1/dj])
L'objet Polynome([-x, 1]) représente en fait le facteur (xj - x) ou (x - xj), et donc l'objet Polynome([-x/dj, 1/dj]) représente logiquement le rapport (x - xj)/(xi- xj).
On peut ainsi développer les polynômes de Lagrange li(x) comme des produits de plusieurs facteurs :
...
# parcours des indices i des valeurs de x : 0 -> len(x)-1
for i in range(len(x)):
# création du polynôme initial : p(x) = 1
pi = Polynome()
# parcours des indices j des valeurs de x : 0 -> len(x)-1
for j in range(len(x)):
if j!=i: # si j différent de i
# p = p*(x - xj)/(xi - xj)
dj = (x - x)
pi = pi * Polynome([-x/dj,1/dj])
...
III-B. Intégration des polynômes développés et réduits
On va ensuite ajouter une méthode à notre classe Polynome pour nous permettre de déterminer la primitive d'une fonction polynomiale :
638782
class Polynome:
...
def primitive(self, cst=0): # renvoie la primitive de la fonction polynomiale f(x)
# initialisation de la liste des coefficients avec des constantes
liste_coefs=*(len(self.coefs)+1) # exemple : [cst, cst, cst, cst, cst]
# parcours des indices des coefficients du polynôme
for i in range(len(self.coefs)):
# calcul du coefficient du terme de degré i+1
liste_coefs[i+1] = self.coefs/(i+1)
# renvoie l'objet Polynome représentant la primitive
return Polynome(liste_coefs)
La primitive d'une fonction polynomiale est obtenue en ajoutant 1 au degré de chaque terme et en divisant chaque terme par son nouveau degré.
Par exemple, pour un polynôme de degré 2 :
p(x) = 2.x2 + x + 1
La primitive est donnée par :
P(x) = (2/3).x3 + (1/2).x2 + x + Cst
Cst désignant une constante arbitraire
En Python, pour obtenir la primitive Pi d'un objet pi représentant par exemple un polynôme de Lagrange développé et réduit, on fait simplement :
Pi = pi.primitive()
La valeur de l'intégrale entre a et b égale au coefficient de quadrature wi est ensuite donnée par :
wi = Pi.eval(b) - Pi.eval(a)
III-C. Intégration complète
La formule d'intégration de degré n est définie ainsi :
638693
Comme on sait déjà évaluer les coefficients de quadrature wi, on peut en déduire facilement le code complet de la fonction Python permettant de calculer l'intégrale numérique de f entre a et b.
La fonction integration_newton_cotes() a comme arguments :
- x : liste des valeurs de x, avec x = a + i*h, pour i = 0, …, n ;
- y : liste des valeurs de y, avec y = f(x).
Elle renvoie une valeur numérique donnant une estimation de la surface située en dessous de la courbe représentative de la fonction f entre a et b :
def integration_newton_cotes(x, y):
# bornes inférieure et supérieure de l'intervalle d'intégration [a, b] : première et dernière valeurs de x
a= x; b = x[-1]
# initialisation de la variable à retourner : valeur de l'intégrale
valeur_integrale = 0
# parcours des indices i des valeurs de x : 0 -> len(x)-1
for i in range(len(x)):
# création du polynôme initial : p(x) = 1
pi = Polynome()
# parcours des indices j des valeurs de x : 0 -> len(x)-1
for j in range(len(x)):
if j!=i: # si j différent de i
# p = p*(x - xj)/(xi - xj)
dj = (x - x)
pi = pi * Polynome([-x/dj,1/dj])
# primitive du polynôme p obtenu
Pi = pi.primitive()
# intégrale entre a et b
wi = (Pi.eval(b) - Pi.eval(a))
# ajout du produit wi*y à la variable valeur_integrale
valeur_integrale += wi*y # valeur_integrale = valeur_integrale + wi*yi
# on retourne la valeur résultat de l'intégration par la formule de Newton-Cotes
return valeur_integrale
En fait, avec le changement de variable u =(x - a)/h, l'expression des coefficients de quadrature peut également s'écrire :
638960
Elle permet d'obtenir la formule exacte de Newton-Cotes qui est un peu plus compliquée à traduire en Python :
def integration_newton_cotes2(x, y):
# bornes inférieure et supérieure de l'intervalle d'intégration [a, b] : première et dernière valeurs de x
a= x; b = x[-1]
# nombre de sous-intervalles [a+ih, a+(i+1)h] avec i ϵ {0, 1, … , n} et h=(b-a)/n
n = len(x)-1
# pas d'intégration sur [a, b]
h = (b-a)/n
# initialisation de la variable à retourner : valeur de l'intégrale
valeur_integrale = 0
# parcours des indices i des valeurs de x : 0 -> n
for i in range(n+1):
# création du polynôme initial : p(x) = 1
pi = Polynome()
# parcours des indices j des valeurs de x : 0 -> n
for j in range(n+1):
if j!=i: # si j différent de i
# changement de variable : uj = (xj -a)/h
uj = (x-a)/h
# pi = pi*(u - uj)
pi = pi * Polynome([-uj,1])
# primitive du polynôme p obtenu
Pi = pi.primitive()
# intégrale entre 0 et n : formule de Newton-Cotes
wi = h*pow(-1,(n-i))/(math.factorial(i)*math.factorial(n-i))*(Pi.eval(float(n)) - Pi.eval(0.0))
# ajout du produit wi*y à la variable valeur_integrale
valeur_integrale += wi*y # valeur_integrale = valeur_integrale + wi*yi
# on retourne la valeur résultat de l'intégration par la formule de Newton-Cotes
return valeur_integrale
III.D Amélioration de la méthode
Bien qu'une formule de Newton-Cotes puisse être établie pour n'importe quel degré, l'utilisation de degrés supérieurs peut causer des erreurs d'arrondi, et la convergence n’est pas assurée lorsque le degré augmente à cause du phénomène de Runge.
Pour cette raison, il est généralement préférable de se restreindre aux premiers degrés, et d'utiliser des formules composites pour améliorer la précision de la formule de quadrature.
L’erreur peut ainsi être réduite en découpant simplement l’intervalle [a, b] en sous-intervalles de m points, ceci dans le but de déterminer une valeur approchée de l’intégrale sur chacun d’eux.
L’intégrale sur [a, b] est estimée par la somme des valeurs ainsi calculées :
638756
La fonction f peut donc être interpolée sur chaque sous-intervalle à l’aide de son évaluation en m points équidistants par un polynôme de degré m – 1 issu d’une base de polynômes de Lagrange et dont l’intégrale est obtenue de la même façon que précédemment.
En notant :
- m : nombre de points équidistants des sous-intervalles [xk, xk+m-1] ;
- k : indice de la 1re valeur de x pour chaque groupe de m points.
k prenant successivement les valeurs 0, (m-1), 2(m-1), .., (n-m+1).
On peut ainsi utiliser la fonction Python précédente sur des sous-listes de longueur m :
valeur_integrale = valeur_integrale + integration_newton_cotes(x[k:k+m], y[k:k+m])
On obtient finalement la fonction générale :
def integration_numerique(x, y, m):
# nombre de sous-intervalles [a+ih, a+(i+1)h] avec i ϵ {0, 1, … , n} et h=(b-a)/n
n = len(x) - 1 # nombre de valeurs de x - 1
if (n % (m-1) != 0): # si n n'est pas un multiple de (m-1)
print(str(n) + " n'est pas un multiple de " + str(m-1) + " !")
return None # on sort de la fonction
# initialisation de la variable à retourner : valeur de l'intégrale
valeur_integrale = 0
# parcours des indices des valeurs de x avec un pas de m-1 : 0, (m-1), 2(m-1), ..., (n-m+1).
for k in range(0, n, m-1):
# calcul de l'intégrale sur un sous-intervalle [x, x[k+m-1]] puis ajout du résultat à la variable valeur_integrale
valeur_integrale += integration_newton_cotes(x[k:k+m], y[k:k+m])
# on retourne la valeur résultat de l'intégration numérique
return valeur_integrale
III.E Module de test
Enfin, on donne le code complet contenant la classe Polynome et les différentes fonctions pour effectuer les tests :
import math # module utilisé pour disposer des fonctions mathématiques : cos(x), sin(x), exp(x), etc.
from decimal import Decimal, getcontext # module utilisé pour convertir les float en decimal : xi = Decimal(a + h*i); yi = Decimal(f(xi))
class Polynome:
def __init__(self, liste_coefs=, domaine=None): # 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 __sub__(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()
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 primitive(self, cst=0): # renvoie la primitive de la fonction polynomiale f(x)
# initialisation de la liste des coefficients qu'avec des constantes
liste_coefs=*(len(self.coefs)+1) # exemple : [cst, cst, cst, cst, cst]
# parcours des indices des coefficients du polynôme
for i in range(len(self.coefs)):
# calcul du coefficient du terme de degré i+1
liste_coefs[i+1] = self.coefs/(i+1)
# renvoie l'objet Polynome représentant la primitive
return Polynome(liste_coefs)
def integration_newton_cotes(x, y):
# bornes inférieure et supérieure de l'intervalle d'intégration [a, b] : première et dernière valeurs de x
a= x; b = x[-1]
# initialisation de la variable à retourner : valeur de l'intégrale
valeur_integrale = 0
# parcours des indices i des valeurs de x : 0 -> len(x)-1
for i in range(len(x)):
# création du polynôme initial : p(x) = 1
pi = Polynome()
# parcours des indices j des valeurs de x : 0 -> len(x)-1
for j in range(len(x)):
if j!=i: # si j différent de i
# p = p*(x - xj)/(xi - xj)
dj = (x - x)
pi = pi * Polynome([-x/dj,1/dj])
# primitive du polynôme p obtenu
Pi = pi.primitive()
# intégrale entre a et b
wi = (Pi.eval(b) - Pi.eval(a))
# ajout du produit wi*y à la variable valeur_integrale
valeur_integrale += wi*y # valeur_integrale = valeur_integrale + wi*yi
# on retourne la valeur résultat de l'intégration par la formule de Newton-Cotes
return valeur_integrale
def integration_newton_cotes2(x, y):
# bornes inférieure et supérieure de l'intervalle d'intégration [a, b] : première et dernière valeurs de x
a= x; b = x[-1]
# nombre de sous-intervalles [a+ih, a+(i+1)h] avec i ϵ {0, 1, … , n} et h=(b-a)/n
n = len(x)-1
# pas d'intégration sur [a, b]
h = (b-a)/n
# initialisation de la variable à retourner : valeur de l'intégrale
valeur_integrale = 0
# parcours des indices i des valeurs de x : 0 -> n
for i in range(n+1):
# création du polynôme initial : p(x) = 1
pi = Polynome()
# parcours des indices j des valeurs de x : 0 -> n
for j in range(n+1):
if j!=i: # si j différent de i
# changement de variable : uj = (xj -a)/h
uj = (x-a)/h
# pi = pi*(u - uj)
pi = pi * Polynome([-uj,1])
# primitive du polynôme p obtenu
Pi = pi.primitive()
# intégrale entre 0 et n : formule de Newton-Cotes
wi = h*pow(-1,(n-i))/(math.factorial(i)*math.factorial(n-i))*(Pi.eval(float(n)) - Pi.eval(0.0))
# ajout du produit wi*y à la variable valeur_integrale
valeur_integrale += wi*y # valeur_integrale = valeur_integrale + wi*yi
# on retourne la valeur résultat de l'intégration par la formule de Newton-Cotes
return valeur_integrale
def integration_numerique(x, y, m):
# nombre de sous-intervalles [a+ih, a+(i+1)h] avec i ϵ {0, 1, … , n} et h=(b-a)/n
n = len(x) - 1 # nombre de valeurs de x - 1
if (n % (m-1) != 0): # si n n'est pas un multiple de (m-1)
print(str(n) + " n'est pas un multiple de " + str(m-1) + " !")
return None # on sort de la fonction
# initialisation de la variable à retourner : valeur de l'intégrale
valeur_integrale = 0
# parcours des indices des valeurs de x avec un pas de m-1 : 0, (m-1), 2(m-1), ..., (n-m+1).
for k in range(0, n, m-1):
# calcul de l'intégrale sur un sous-intervalle [x, x[k+m-1]] puis ajout du résultat à la variable valeur_integrale
valeur_integrale += integration_newton_cotes2(x[k:k+m], y[k:k+m])
# on retourne la valeur résultat de l'intégration numérique
return valeur_integrale
# donne la précision des nombres décimaux utilisés dans les calculs
#getcontext().prec = 32
# fonction à intégrer sur l'intervalle [a, b]
f = lambda x: math.exp(-x/2)
# limite inférieure et supérieure de l'intervalle d'intégration : intégration de la fonction f sur [a, b]
a= 0.0; b = 10.0
# nombre de sous-intervalles [a+ih, a+(i+1)h] avec i ϵ {0, 1, … , n} et h=(b-a)/n
n = 20
# longueur de chaque sous-intervalle
h = (b-a)/n
# intialisation des listes de valeurs x et y
x=[]; y=[]
# parcours des indices des n+1 valeurs de x et y
for i in range(n+1):
xi = a + h*i # valeur de xi, pour convertir un float en decimal : xi = Decimal(a + h*i)
yi = f(xi) # valeur de yi, pour convertir un float en decimal : yi = Decimal(f(xi))
x.append(xi) # ajout de xi à la liste x
y.append(yi) # ajout de yi à la liste y
print("I. Intégration numérique de la fonction f(x)={0} par la formule de Newton-Cotes".format("exp(-x/2)"))
print()
print("Intégration numérique sur l'intervalle [{0}, {1}] avec {2} points :".format(a,b,n+1))
print()
# valeur approchée résultat de l'intégration numérique de la fonction f sur l'intervalle [a, b]
valeur_approchee = integration_newton_cotes(x, y)
print("Valeur approchée de l'intégrale : " + str(valeur_approchee))
# valeur réelle de l'intégrale obtenue à l'aide de la primitive de la fonction f(x)
valeur_reelle = (-2*math.exp(-b/2) + 2*math.exp(-a/2))
print("Valeur réelle de l'intégrale : " + str(valeur_reelle))
print();print()
print("II. Intégration numérique de la fonction f(x)={0}".format("exp(-x/2) à l'aide de formules composites"))
print()
# découpage du segment [a, b] en sous-intervalles de longueur m*h : on prendra (n multiple de (m-1)
m = 3 # nombre de points équidistants : méthode de Simpson pour m=3 points
print("Intégration numérique sur l'intervalle [{0}, {1}] avec {2} points et des sous-intervalles de {3} points équidistants :".format(a,b,n+1,m))
print()
# valeur approchée résultat de l'intégration numérique de la fonction f sur l'intervalle [a, b] avec des sous-intervalles de m points équidistants
valeur_approchee = integration_numerique(x, y, m)
print("Valeur approchée de l'intégrale : " + str(valeur_approchee))
print("Valeur réelle de l'intégrale : " + str(valeur_reelle))
Note : il utilise en plus la librairie decimal pour éviter les erreurs dans les opérations de nombres à virgule flottante.
Le code affiche :
638949
IV. Conclusion
Après avoir présenté la méthode d'intégration de Newton-Cotes, nous avons pu l'implémenter en Python à l'aide de la classe Polynome. Enfin, nous avons proposé une amélioration de cette intégration numérique en utilisant des formules composites.
La principale difficulté est donc d'obtenir les polynômes de Lagrange li(x) sous leur forme développée et réduite, ensuite leur intégration sur chaque sous-intervalle est assez simple à réaliser.
Comme on a pu le noter, on peut également utiliser la formule de Newton-Cotes obtenue après un changement de variable. Elle donne des résultats plus précis, mais elle est aussi plus compliquée à mettre en œuvre.
Sources :
https://fr.wikipedia.org/wiki/Formule_de_Newton-Cotes
https://fr.wikipedia.org/wiki/Calcul_num%C3%A9rique_d%27une_int%C3%A9grale
https://fr.wikipedia.org/wiki/Interpolation_lagrangienne
https://fr.wikipedia.org/wiki/Ph%C3%A9nom%C3%A8ne_de_Runge
https://fr.wikipedia.org/wiki/Calcul_num%C3%A9rique_d%27une_int%C3%A9grale#Formules_composites
https://www.developpez.net/forums/blogs/44027-user/b10406/creer-classe-polynome-python/
https://docs.python.org/fr/3/library/decimal.html
https://docs.python.org/fr/3/tutorial/floatingpoint.html
Soutenez le club developpez.com en souscrivant un abonnement pour que nous puissions continuer à vous proposer des publications.