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