Feat: Ajoute newton_cotes et monte_carlo_2d

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