258 lines
7.3 KiB
Python
258 lines
7.3 KiB
Python
import random
|
|
|
|
import numpy as np
|
|
import matplotlib.pyplot as plt
|
|
|
|
|
|
def trapeze_formule(f, a: float, b: float, n: int) -> float:
|
|
"""
|
|
|
|
:param f: la fonction f (bien passer une fonction ou une lambda)
|
|
:param a: la petite borne
|
|
:param b: la grande borne
|
|
:param n: nombre de "tranche" dans l'intervalle
|
|
:return: l'approximation de la valeur de l'intervalle
|
|
"""
|
|
|
|
# largeur de chaque sous-intervalle
|
|
dx = (b - a) / n
|
|
|
|
# la somme des valeurs de f aux points entre a et b, f(a) f(b) exclus
|
|
somme_inferieure = 0.
|
|
|
|
for i in range(1, n):
|
|
xi = a + i * dx
|
|
somme_inferieure += f(xi)
|
|
|
|
# I = (delta X / 2) * [f(a) + f(b) + 2 * somme_inferieure ]
|
|
I = (dx / 2) * (f(a) + f(b) + 2 * somme_inferieure)
|
|
return I
|
|
|
|
def trapeze_version_numpy(f, a: float, b: float, n: int) -> float:
|
|
"""
|
|
|
|
:param f: la fonction f (bien passer une fonction ou une lambda)
|
|
:param a: la petite borne
|
|
:param b: la grande borne
|
|
:param n: nombre de "tranche" dans l'intervalle
|
|
:return: l'approximation de la valeur de l'intervalle
|
|
"""
|
|
|
|
# largeur de chaque sous-intervalle
|
|
dx = (b - a) / n
|
|
x = np.linspace(a, b, n+1)
|
|
y = f(x)
|
|
return dx * (np.sum(y) - 0.5*(y[0] + y[-1]))
|
|
|
|
|
|
def simpson(f, a: float, b: float, n: int) -> float:
|
|
"""
|
|
|
|
:param f: la fonction f(x)
|
|
:param a: la petite borne
|
|
:param b: la grande borne
|
|
:param n: nombre de "tranche" dans l'intervalle
|
|
:return: l'approximation de la valeur de l'intervalle
|
|
"""
|
|
|
|
if n % 2 != 0:
|
|
raise ValueError(f'n must be even, got {n}')
|
|
|
|
# largeur de chaque sous-intervalle
|
|
dx = (b - a) / n
|
|
|
|
somme_paire = 0.
|
|
somme_impaire = 0.
|
|
|
|
for i in range(1, n):
|
|
xi = a + i * dx
|
|
if i % 2 == 0:
|
|
somme_paire += f(xi)
|
|
else:
|
|
somme_impaire += f(xi)
|
|
|
|
I = (dx / 3) * (f(a) + f(b) + 4 * somme_impaire + 2 * somme_paire)
|
|
return I
|
|
|
|
|
|
def simpson_numpy(f, a: float, b: float, n: int) -> float:
|
|
if n % 2 != 0:
|
|
raise ValueError(f'n must be even, got {n}')
|
|
|
|
# largeur de chaque sous-intervalle
|
|
dx = (b - a) / n
|
|
x = np.linspace(a, b, n + 1)
|
|
y = f(x)
|
|
|
|
coefficients = np.ones(n + 1) # Crée un tableau de 1.0 de taille n+1
|
|
coefficients[1:n:2] = 4 # Remplace les indices impairs (1, 3, 5, ...) par 4
|
|
coefficients[2:n:2] = 2 # Remplace les indices pairs (2, 4, 6, ...) par 2
|
|
|
|
return (dx / 3) * np.sum(coefficients * y)
|
|
|
|
|
|
# Les b possibles dans Newton Cotes
|
|
POIDS = {
|
|
2: [1/2,1/2],
|
|
3: [1/6, 4/6, 1/6],
|
|
4: [1/8, 3/8, 3/8, 1/8],
|
|
5: [
|
|
7 / 90,
|
|
32 / 90,
|
|
12 / 90,
|
|
32 / 90,
|
|
7 / 90,
|
|
],
|
|
6: [19/288, 75/288, 50/288, 50/288, 75/255, 19/255],
|
|
7: [41/840, 216/840, 27/840, 272/840, 27/840, 216/840, 41/840]
|
|
}
|
|
|
|
def newton_cotes(f, a: float, b: float, s: int, n: int) -> float:
|
|
"""
|
|
|
|
:param f: la fonction f(x)
|
|
:param a: la petite borne
|
|
:param b: la grande borne
|
|
:param s: l'ordre (faire +1 pour savoir quel poid utiliser) Cela détermine comme on approxime la courbe de la fonction
|
|
Une courbe (méthode des trapèzes), un polynôme de deg. 2 (Simpson), etc.
|
|
:param n: Nombre de "tranche" dans l'intervalle
|
|
:return: l'approximation de la valeur de l'intervalle
|
|
"""
|
|
|
|
weights = POIDS[s+1] # Poids pour l'ordre "s"
|
|
h = (b - a) / n
|
|
Aj = 0.
|
|
|
|
for j in range(n):
|
|
somme = 0.
|
|
x_j = a + j * h
|
|
for i in range(s):
|
|
x = x_j + (i / (s - 1)) * h # points à l'intérieur du sous-intervalle
|
|
somme += weights[i] * f(x)
|
|
somme += h / (s-1) # facteur h/(s-1) selon Newton-Cotes
|
|
Aj += somme * h
|
|
|
|
return Aj
|
|
|
|
|
|
def newton_cotes_numpy(f, a: float, b: float, s: int, n: int) -> float:
|
|
"""
|
|
Approximation de l'intégrale de f sur [a, b] par la formule de Newton-Cotes d'ordre s.
|
|
|
|
:param f: Fonction à intégrer.
|
|
:param a: Borne inférieure.
|
|
:param b: Borne supérieure.
|
|
:param s: Ordre de la formule (1: trapèzes, 2: Simpson, etc.).
|
|
:param n: Nombre de sous-intervalles (doit être divisible par s pour les formules fermées).
|
|
:return: Approximation de l'intégrale.
|
|
"""
|
|
|
|
if n % s != 0:
|
|
raise ValueError(f'n must be a multiple of s, got s = {s} and n = {n}')
|
|
|
|
# b_i dans le tableau
|
|
weights_by_orders = {
|
|
2: [1 / 2, 1 / 2], # trapèzes
|
|
3: [1 / 6, 4 / 6, 1 / 6], # simpson
|
|
4: [1 / 8, 3 / 8, 3 / 8, 1 / 8], # simpson 3/8
|
|
5: [
|
|
7 / 90,
|
|
32 / 90,
|
|
12 / 90,
|
|
32 / 90,
|
|
7 / 90,
|
|
], # bool
|
|
6: [19 / 288, 75 / 288, 50 / 288, 50 / 288, 75 / 255, 19 / 255],
|
|
7: [41 / 840, 216 / 840, 27 / 840, 272 / 840, 27 / 840, 216 / 840, 41 / 840]
|
|
}
|
|
|
|
weights = weights_by_orders[s+1]
|
|
h = (b - a) / n
|
|
x = np.linspace(a, b, n * s + 1) # Tous les points d'évaluation
|
|
y = f(x)
|
|
|
|
# Applique les poids de manière glissante sur chaque sous-intervalle
|
|
return np.sum(
|
|
np.convolve(y, weights[::-1], mode='valid') * (h / (s + 1)) # Convolution pour appliquer les poids
|
|
)
|
|
|
|
|
|
def exercice1():
|
|
f = lambda x: x ** 2
|
|
approx = trapeze_formule(f, a=0, b=1, n=10)
|
|
print(approx)
|
|
|
|
g = lambda x : np.sin(x)
|
|
trapeze1 = trapeze_formule(g, a=0, b=np.pi, n=100)
|
|
trapeze2 = trapeze_formule(g, a=0, b=np.pi, n=200)
|
|
|
|
giga_trapeze1 = trapeze_version_numpy(g, 0, np.pi, 100)
|
|
giga_trapeze2 = trapeze_version_numpy(g, 0, np.pi, 200)
|
|
|
|
print(f'VERSION MOI {trapeze1}, VERSION CHAT {giga_trapeze1}')
|
|
print(f'VERSION MOI {trapeze2}, VERSION CHAT {giga_trapeze2}')
|
|
|
|
simpson1 = simpson(g, a=0, b=np.pi, n=100)
|
|
simpson2 = simpson(g, a=0, b=np.pi, n=200)
|
|
|
|
giga_simpson1 = simpson_numpy(g, 0, np.pi, 100)
|
|
giga_simpson2 = simpson_numpy(g, 0, np.pi, 200)
|
|
|
|
print(f'MOI: {simpson1}, CHAT : {giga_simpson1}')
|
|
print(f'MOI: {simpson2}, CHAT : {giga_simpson2}')
|
|
|
|
newton1 = newton_cotes(g, a=0, b=np.pi, n=100, s=4)
|
|
newton2 = newton_cotes(g, a=0, b=np.pi, n=200, s=4)
|
|
|
|
# le vrai résultat est 2.
|
|
print(f'{trapeze1}, {trapeze2} erreur: {2-trapeze1} {2-trapeze2}')
|
|
print(f'{simpson1}, {simpson2} erreur: {2 - simpson1} {2 - simpson2}')
|
|
print(f'{newton1}, {newton2} erreur: {2 - newton1} {2 - newton2}')
|
|
|
|
|
|
def monte_carlo_2d(f, a: float, b: float, c: float, d: float, n: int) -> float:
|
|
"""
|
|
|
|
:param f: la fonction
|
|
:param a: borne a (pour x)
|
|
:param b: borne b (pour x)
|
|
:param c: borne c (pour y)
|
|
:param d: borne d (pour y)
|
|
:param n: nombre de tirages aléatoire (N grand → epsilon petit)
|
|
:return: I
|
|
"""
|
|
sum = 0.
|
|
points_x = []
|
|
points_y = []
|
|
for i in range(1, n):
|
|
x = a + (b - a) * random.random()
|
|
y = c + (d - c) * random.random()
|
|
if a < x < b and c < y < d:
|
|
sum += f(x, y)
|
|
points_x.append(x)
|
|
points_y.append(y)
|
|
|
|
I = (sum / n) * (b - a) * (d - c)
|
|
epsilon = np.abs(.58 - I)
|
|
print(f'Err absolue : {epsilon}')
|
|
|
|
plt.scatter(points_x, points_y, color='blue')
|
|
plt.title('Surface fonction')
|
|
plt.show()
|
|
|
|
return I
|
|
|
|
|
|
def exercice2():
|
|
g = lambda x, y: 1
|
|
g2 = lambda x, y: np.exp(-x ** 2 - y ** 2)
|
|
|
|
ig = monte_carlo_2d(g, a=-1, b=1, c=-1, d=1, n=10_000)
|
|
ig2 = monte_carlo_2d(g2, a=-1, b=1, c=-1, d=1, n=10_000)
|
|
print(f'I (g) = {ig}, I (g2) = {ig2}')
|
|
|
|
|
|
if __name__ == '__main__':
|
|
exercice1()
|
|
#exercice2()
|