Méthodes d'intégration trapèze et simpson

Version imprimable

Bonjour,
je suis en train de programmer en python la méthode des trapèzes et la méthode de Simpson mais je suis confronté à un problème : Je ne retrouve pas un ordre de 2 pour la méthode des trapèzes et pas un ordre de 4 pour la méthode de Simpson. Pour trapèze, je trouve une ordre de 1,06 et pour simpson un ordre de -0,09... . De plus, pour la méthode de simpson, j'ai n'ai pas la bonne valeur de l'integrale...

Je rappelle ci-dessous les formules :

Trapèze : I = h*[f_0 + f_{N-1} + 2*(sigma_{i=1}^{N-2}(f_i))]/2 + O(h^2)

Simpson : I = h*(f_0 + f_{N-1} + 2*(sigma_{i=2}^{N-3}(f_i)) + 4*(sigma_{1}^{N-2}(f_i)))/3

Voici mon programme en python :

Code:

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
import numpy as np
import matplotlib.pyplot as plt
 
#=================================
# TD sur l'intégration numérique (transposé depuis un code matlab)
#=================================
 
# =============================
# Définition de toutes les fonctions utile
# =============================
 
def intrectangle(x,y):
    N = len(x)
    h = x[2]-x[1]
    s = 0
    for i in range(1,N):
        s+= y[i]
    return h*s
 
def inttrapeze(x,y):
    N = len(x)
    h = x[2]-x[1]
    s = 0
    for i in range(1,N-1):
        s += y[i]
    return h*(y[0]+y[N-1] + 2*s)/2.0
 
def intsimpson(x, y):
    N = len(x)
    h = x[2]-x[1]
    impaire = 0
    paire = 0
    for i in range(1, N-3, 2):
        impaire += y[i]
    for j in range(2, N-4, 2):
        paire += y[i]
    return h*(y[0]+y[N-1]+2*paire+ 4*impaire)/3
 
 
# =============================
# Initialisation de nos différents paramètres
# =============================
 
a = 0
b = 2
N = 501 # Nombre impaire pour pouvoir utiliser Simpson
x = np.arange(a, b, (b-a)/(N-1)) #Notre vecteur x par pas de h = (b-a)/(N-1)
y = np.exp(x)
I_th = np.exp(b)-np.exp(a) #Valeur théorique
 
# Test de nos différentes fonctions
 
I_rect = intrectangle(x, y)
I_trap = inttrapeze(x, y)
I_simpson = intsimpson(x, y)
 
print('Valeur exacte : ', I_th)
print('Valeur approché par Rectangle : ', I_rect)
print('Valeur approché par Trapèze : ', I_trap)
print('Valeur approché par Simpson : ', I_simpson)
 
 
# Plot de l'erreur relatif en fonction de pas h
 
H = []
err_rect = []
err_trap = []
err_simpson = []
 
for i in range(5, 3000, 20):
    X = np.arange(a, b, (b-a)/(i-1))
    Y = np.exp(X)
    H.append(X[1]-X[0])
    err_rect.append(abs(I_th-intrectangle(X, Y))/I_th)
    err_trap.append(abs(I_th-inttrapeze(X, Y))/I_th)
    err_simpson.append(abs(I_th-intsimpson(X, Y))/I_th)
 
plt.loglog(H, err_rect, H, err_trap, H, err_simpson)
plt.show()
 
# Utilisation de np.polyfit pour retrouver coefficient directeur de la droite obtenue
print(np.polyfit(np.log(H), np.log(err_rect), 1))
print(np.polyfit(np.log(H), np.log(err_trap), 1))
print(np.polyfit(np.log(H), np.log(err_simpson), 1))

:salut:

J'obtiens :

Valeur exacte :  6.38905609893065
Valeur approché par Rectangle :  6.372286505471986
Valeur approché par Trapèze :  6.359567387655495
Valeur approché par Simpson :  6.2718890256301885
[ 0.9990951  -0.43645248]
[1.06360079 0.23654477]
[0.95853913 1.23298127]

après une petite correction :

Code:

1
2
3
4
5
6
7
8
9
10
def intsimpson(x, y):
    N = len(x)
    h = x[2]-x[1]
    impaire = 0
    paire = 0
    for i in range(1, N-3, 2):
        impaire += y[i]
    for j in range(2, N-4, 2):
        paire += y[j]
    return h*(y[0]+y[N-1]+2*paire + 4*impaire)/3

22/11/2019, 23h01
cercus

Effectivement mais cela ne résous pas le problème : on trouve un ordre de 1 alors qu'on doit trouver un ordre 4...
25/11/2019, 14h28
lg_53
Outre l'erreur que daniel a relevé (un mauvais copier/coller dans lequel un indice i était resté alors que là c'était l'indice j), voici des éléments de réponse :

1) Les tests effectués ne sont pas équivalent, car vous ne maillez pas l'intervalle [1,2], mais l'intervalle [1,2[. Et oui car avec x=np.arange(a,b), b n'est pas dans x ! Il faut utiliser np.linspace(a,b,n).

2) une petite légende pour agrémenter les courbes ne serait pas de trop :
Code:

1 2 plt.legend(['rectangle', 'trapeze','simpson']) plt.show()
3) Ce code ci
Code:

1 2 3 4 5 6 7 8 9 10 def intsimpson(x,y): N=len(X) h = x[2]-x[1] ### valable que si N est impair s=0.0 for i in range(1,N-1,2): s+=y[i-1]+4*y[i]+y[i+1] return s*h/3.0
fournit bien un ordre 4 chez moi. Donc vous devez vous emmelez les pinceaux en disséquant le cas paire/impaire + les bords

Remarque: le code ne fonctionne que si x est discrétisé de manière régulière. Donc au lieu de passer x en paramètre, on devrait plutot passer directement h. Car si vous passer en entré un x qui ressemble à ca par exemple

Code:

x = [1.0 ,1.5, 2.0, 3.0, 3.2, 3.4]

votre fonction va renvoyer qqch, mais qui sera juste complètement faut

Pour debug, on peut vérifier la valeur à atteindre avec scipy.integrate.simps

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
from scipy.integrate import simps
 
def intsimpson(x, y):
    return simps(y, x)
 
 
Valeur exacte :  6.38905609893065
Valeur approché par Rectangle :  6.401842729867709
Valeur approché par Trapèze :  6.389064617669847
Valeur approché par Simpson :  6.3890560989397365
[ 1.00368864 -0.66879289]
[ 1.99985028 -2.48587882]
[ 4.00293553 -5.18001168]

Pour le plaisir :

Code:

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
from termcolor import cprint
from typing import List, Callable
import numpy as np
from scipy.integrate import simps
 
 
def rectangles(f: Callable, a: int, b: int, n: int) -> float:
    S = 0
    for i in range(0, n):
        Xi = a + (b - a) * i/float(n)
        Xj = a + (b - a) * (i + 1)/float(n)
        S += f((Xi + Xj)/2.0) * (Xj - Xi)
    return S
 
 
def trapezes(f: Callable, a: int, b: int, n: int) -> float:
    S = 0
    for i in range(0, n):
        Xi = a + (b - a) * i/float(n)
        Xj = a + (b - a) * (i + 1)/float(n)
        S += (f(Xi) + f(Xj))/2.0 * (Xj - Xi)
    return S
 
 
def simpson(f: Callable, a: int, b: int, n: int) -> float:
    S = 0
    for i in range(0, n):
        Xi = a + (b - a) * i/float(n)
        Xj = a + (b - a) * (i + 1)/float(n)
        S += (Xj - Xi) * (f(Xi) + 4.0*f((Xi + Xj)/2.0) + f(Xj))/6.0
    return S
 
 
a: int = 0
b: int = 2
n: int = 500
 
 
def fn(x: float) -> float:
    return np.exp(x)
 
 
cprint('rectangles = {}'.format(rectangles(fn, a, b, n)), 'green')
cprint('trapezes = {}'.format(trapezes(fn, a, b, n)), 'green')
cprint('simpson = {}'.format(simpson(fn, a, b, n)), 'green')
 
# debug
 
X: List[float] = []
Y: List[float] = []
h: float = (b - a)/n
 
for x in np.arange(a, b + h, h):
    X.append(x)
    Y.append(fn(x))
 
cprint('debug simpson = {}'.format(simps(Y, X)), 'red')
 
'''
rectangles = 6.38905183956191
trapezes = 6.38906461766984
simpson = 6.389056098931216
debug simpson = 6.3890560989397365
'''