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