User

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.

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 x_i = a + ih, pour i = 0, … , n et h = (b – a)/n. La formule d'intégration de degré n est définie ainsi :

où les w_i 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 (x_i, f(x_i)) et l_i(x) le polynôme de Lagrange tel que :

Alors :

Ainsi :

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 l_i(x) pour obtenir les coefficients de quadrature w_i.

Enfin, on obtient la valeur de l'intégrale de f sur l'intervalle [a, b] en utilisant la formule :

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) = 2x² + x + 1 = 1 + x + 2x²

à l'aide de sa liste de coefficients [1, 1, 2] :

Code Python :

Sélectionner tout - Visualiser dans une fenêtre à part

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 :

Code Python :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
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 :

p_i = p_i*(x - x_j)/(x_i - x_j)

Equivaut en utilisant des objets Polynome à :

Code Python :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
dj = (x[i] - x[j])
pi = pi * Polynome([-x[j]/dj, 1/dj])

L'objet Polynome([-x[j], 1]) représente en fait le facteur (x_j - x) ou (x - x_j), et donc l'objet Polynome([-x[j]/dj, 1/dj]) représente logiquement le rapport (x - x_j)/(x_i- x_j).

On peut ainsi développer les polynômes de Lagrange l_i(x) comme des produits de plusieurs facteurs :

Code Python :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
...
    # 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([1]) 
 
        # 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[i] - x[j])
                pi = pi * Polynome([-x[j]/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 :

Code Python :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
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=[cst]*(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]/(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.x² + x + 1

La primitive est donnée par :

P(x) = (2/3).x³ + (1/2).x² + 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 :

Code Python :

Sélectionner tout - Visualiser dans une fenêtre à part

Pi = pi.primitive()

La valeur de l'intégrale entre a et b égale au coefficient de quadrature w_i est ensuite donnée par :

Code Python :

Sélectionner tout - Visualiser dans une fenêtre à part

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 :

Comme on sait déjà évaluer les coefficients de quadrature w_i, 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[i] = a + i*h, pour i = 0, …, n ;
• y : liste des valeurs de y, avec y[i] = f(x[i]).

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 :

Code Python :

Sélectionner tout - Visualiser dans une fenêtre à part

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
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[0]; 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([1]) 
 
        # 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[i] - x[j])
                pi = pi * Polynome([-x[j]/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[i] à la variable valeur_integrale
        valeur_integrale  += wi*y[i] # 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 :

Elle permet d'obtenir la formule exacte de Newton-Cotes qui est un peu plus compliquée à traduire en Python :

Code Python :

Sélectionner tout - Visualiser dans une fenêtre à part

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
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[0]; 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([1]) 
 
        # 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[j]-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[i] à la variable valeur_integrale
        valeur_integrale  += wi*y[i] # 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 :

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 [x_k, x_k+m-1] ;
• k : indice de la 1^re 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 :

Code Python :

Sélectionner tout - Visualiser dans une fenêtre à part

valeur_integrale = valeur_integrale + integration_newton_cotes(x[k:k+m], y[k:k+m])

On obtient finalement la fonction générale :

Code Python :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
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[k], 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 :

Code Python :

Sélectionner tout - Visualiser dans une fenêtre à part

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
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
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=[0], 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[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 __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[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]) 
 
        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 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=[cst]*(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]/(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[0]; 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([1]) 
 
        # 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[i] - x[j])
                pi = pi * Polynome([-x[j]/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[i] à la variable valeur_integrale
        valeur_integrale  += wi*y[i] # 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[0]; 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([1]) 
 
        # 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[j]-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[i] à la variable valeur_integrale
        valeur_integrale  += wi*y[i] # 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[k], 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 :

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 l_i(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...int%C3%A9grale
https://fr.wikipedia.org/wiki/Interp...n_lagrangienne
https://fr.wikipedia.org/wiki/Ph%C3%...%A8ne_de_Runge
https://fr.wikipedia.org/wiki/Calcul...les_composites
https://www.developpez.net/forums/bl...lynome-python/
https://docs.python.org/fr/3/library/decimal.html
https://docs.python.org/fr/3/tutoria...tingpoint.html