[TP0] Méthodes numériques

2020-09-01 17:53:50 +02:00 · 2020-09-01 17:53:50 +02:00 · 49d542c19d
commit 49d542c19d
parent ed8a5bed4f
1 changed files with 107 additions and 0 deletions
--- a/NUMERIQUES/main.py
+++ b/NUMERIQUES/main.py
@ -0,0 +1,107 @@
+import numpy as np
+
+class Question:
+    def __init__(self,i=1, l=0, n=''):
+        self.name=n
+        self.level=l
+        self.number=i
+    def __enter__(self):
+        print('\n' + (self.level*2)*' ' + f"-> {self.number}. : {self.name} -- Début")
+        return self
+    def __exit__(self, exc_type, exc_value, exc_traceback):
+        print((self.level*2)*' ' + f"<- {self.number}. : {self.name} -- Fin\n")
+
+
+with Question(1):
+    with Question(1,1):
+        def dichotomie(f, n, a=0, b=2):
+            iterations = 0
+            while abs(b-a) >= 2*10**(-n):
+                iterations += 1
+                c = (a+b) / 2
+                if f(c) == 0:
+                    return c
+                elif f(a)*f(c) > 0:
+                    a = c
+                else:
+                    b = c
+            return c, iterations
+        f = lambda x: x**3 - 3*x**2 + 1
+        print(dichotomie(f,5,0,1.5))
+        print(dichotomie(np.sin, 3, 3, 4))
+        print(dichotomie(np.sin, 10, 3, 4))
+
+
+    with Question(2,1):
+        def newton(f, fp, x0, nbiter, epsilon=1e-10):
+            for i in range(nbiter):
+                x1 = x0 - f(x0) / fp(x0)
+                if abs(x1 - x0) <= epsilon:
+                    break
+                x0 = x1
+            return x1
+        print(newton(lambda x: x**3 - 3*x**2 + 1, lambda x: 3*x**2 - 6*x, 1.5, 200, epsilon=1e-100))
+        print(newton(np.sin, np.cos, 3, 100))
+        print(abs(newton(np.sin, np.cos, 3, 100) - dichotomie(np.sin, 10, 3, 4)[0]))
+
+
+with Question(2):
+    from math import cos, exp, sin, log
+    from scipy.integrate import quad
+    from time import time
+    with Question(1,2):
+        def rectangles(f, a, b, n):
+            res = 0
+            for k in range(n):
+                res += f(a + k*(b - a)/n)
+            return res*(b - a)/n
+        def f1(x):
+            return x
+        def f2(x):
+            return x**10
+        f3 = cos
+        f4 = exp
+        print(rectangles(f1, 0, 1, 1000), quad(f1, 0, 1)[0])
+        print(rectangles(f2, 0, 1, 1000), quad(f2, 0, 1)[0])
+        print(rectangles(f3, 0, np.pi/2, 1000), quad(f3, 0, np.pi/2)[0])
+        print(rectangles(f4, -3, 3, 1000), quad(f4, -3, 3)[0])
+
+
+    with Question(2,2):
+        def trapezes(f, a, b, n):
+            res = 0
+            for k in range(n):
+                res += f(a + k*(b - a)/n) + f(a + (k+1)*(b - a)/n)
+            return res*(b - a)/(2*n)
+        print(trapezes(f1, 0, 1, 1000), quad(f1, 0, 1)[0])
+        print(trapezes(f2, 0, 1, 1000), quad(f2, 0, 1)[0])
+        print(trapezes(f3, 0, np.pi/2, 1000), quad(f3, 0, np.pi/2)[0])
+        print(trapezes(f4, -3, 3, 1000), quad(f4, -3, 3)[0])
+
+
+    with Question(3, 2):
+        def simpson(f, a, b, n):
+            res = 0
+            for k in range(n):
+                res += f(a + k*(b - a)/n) + 4*f(((a+k*(b - a)/n) + (a + (k+1)*(b - a)/n))/2) + f(a + (k+1)*(b - a)/n)
+            return res*(b - a)/(6*n)
+        print(simpson(f1, 0, 1, 1000), quad(f1, 0, 1)[0])
+        print(simpson(f2, 0, 1, 1000), quad(f2, 0, 1)[0])
+        print(simpson(f3, 0, np.pi/2, 1000), quad(f3, 0, np.pi/2)[0])
+        print(simpson(f4, -3, 3, 1000), quad(f4, -3, 3)[0])
+
+
+with Question(3):
+    import matplotlib.pyplot as plt
+    with Question(1, 2):
+        def euler(F, t0, y0, t1, n):
+            res = [y0]*(n+1)
+            for i in range(1,n+1):
+                res[i] = res[i-1] + (t1-t0)/n * F(t0 + i*(t1-t0)/n, res[i-1])
+            return res
+        def expo_euler(t, yt):
+            return yt
+        n = 100
+        plt.plot(np.linspace(0, 1, num=n+1), euler(expo_euler, 0, 1, 1, n))
+        plt.plot(np.linspace(0, 1, num=n+1), np.exp(np.linspace(0, 1, num=n+1))) 
+        plt.show()