diff --git a/tp7.py b/tp7.py index 1629b06..73f718e 100644 --- a/tp7.py +++ b/tp7.py @@ -1,4 +1,7 @@ +import random + import numpy as np +import matplotlib.pyplot as plt def trapeze_formule(f, a: float, b: float, n: int) -> float: @@ -45,6 +48,39 @@ def simpson(f, a: float, b: float, n: int) -> float: I = (dx / 3) * (f(a) + f(b) + 4 * somme_impaire + 2 * somme_paire) return I +# 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: + 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 exercice1(): f = lambda x: x ** 2 @@ -58,12 +94,57 @@ def exercice1(): simpson1 = simpson(g, a=0, b=np.pi, n=100) simpson2 = simpson(g, a=0, b=np.pi, n=200) + 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=0, b=1, c=0, d=1, n=10_000) + ig2 = monte_carlo_2d(g2, a=0, b=1, c=0, d=1, n=10_000) + print(f'I (g) = {ig}, I (g2) = {ig2}') if __name__ == '__main__': - exercice1() + #exercice1() + exercice2()