commit bdeaa7e415e1cfbfcafcef3ea663a103178beb0d Author: Namu Date: Wed Oct 15 15:32:14 2025 +0200 Feat: Add tp5 diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..fb192ea --- /dev/null +++ b/.gitignore @@ -0,0 +1,5 @@ +# jetbrains +.idea/ + +# python venv +.venv/ diff --git a/ex1.py b/ex1.py new file mode 100644 index 0000000..d22a1ac --- /dev/null +++ b/ex1.py @@ -0,0 +1,42 @@ +import numpy as np + + +def resolve_up(A: np.array, b: np.array) -> np.array: + x = np.zeros(len(b)) + for i in range(len(A) - 1, -1, -1): + right_of_diagonal = 0.0 + for k in range(i+1, len(A)): + right_of_diagonal += A[i, k] * x[k] + x[i] = (b[i] - right_of_diagonal) / A[i,i] + return x + + +def resolve_down(A: np.array, b: np.array) -> np.array: + x = np.zeros(len(b)) + for i in range(len(A)): + sum = 0.0 + for k in range(i): + sum += A[i,k] * x[k] + x[i] = (b[i] - sum) / A[i, i] + return x + + +if __name__ == '__main__': + A = np.array([ + [2, 4, -6], + [0, -1, 1], + [0, 0, -2], + ], dtype=float) + + b = np.array([2, 3, -7], dtype=float) + + A_down = np.array([ + [2, 0, 0], + [4, -1, 0], + [-6, 1, -2], + ]) + + b_down = np.array([2, 3, -7]) + + print(resolve_up(A, b)) + print(resolve_down(A_down, b_down)) diff --git a/ex2.py b/ex2.py new file mode 100644 index 0000000..d3b0dc1 --- /dev/null +++ b/ex2.py @@ -0,0 +1,54 @@ +""" +Ici, on implémente gauss +""" +import numpy as np +from ex1 import resolve_up + + +def gauss(A: np.array, b: np.array) -> (np.array, np.array): + n = len(A) + for k in range(n): # k = pivot + for i in range(k+1, n): # i = ligne en dessous du pivot + g = A[i, k] / A[k, k] # multiplication pour éliminer A[i, k] + A[i, k:] -= g * A[k, k:] # row operation on A + b[i] -= g * b[k] # same row operation on b + + return A, b # matrice triangularisé et le b modifié + + +def cholesky(A: np.array) -> np.array: + L = np.zeros(A) + L = np.sqrt(A[0, 0]) + m = len(A) + # first column + for i in range(1, m): + L[i, 0] = A[i, 0] / L[0, 0] + + # généralisation k = 2, ..., m + for k in range(1, m): + # diagonale + L[k, k] = np.sqrt(A[k, k] - np.sum(L[k, :k]**2)) + for i in range(k+1, m): + for j in range(k-1): + L[i,k] = (A[i,k] - np.sum(L[i, j] * L[k, j])) / L[k, k] + return L + + +if __name__ == '__main__': + A = np.array([ + [1, 1, 1, 1], + [1, 5, 5, 5], + [1, 5, 14, 14], + [1, 5, 14, 15] + ], dtype=np.float64) + + b = np.array([1, 0, 0, 0], dtype=np.float64) + + A_triangle, b_gauss = gauss(A, b) + print(f"A triangulsarisée:\n {A_triangle}") + + x = resolve_up(A, b_gauss) + print(f"x={x}") + + L = cholesky(A) + print(f"L\n {L}") diff --git a/ex3.py b/ex3.py new file mode 100644 index 0000000..e1f6623 --- /dev/null +++ b/ex3.py @@ -0,0 +1,136 @@ +import numpy as np + + +# L'inverse d'une diagonale +# D_inv = np.diag(1 / np.diag(A)) + + +def to_D(A: np.array) -> np.array: + D = np.zeros_like(A) + for i in range(len(A)): + D[i, i] = A[i, i] + return D + + +def to_L(A: np.array) -> np.array: + L = np.zeros_like(A) + for i in range(len(A)): + for j in range(len(A)): + if i < j: + L[i, j] = A[i, j] + return L + + +def to_U(A: np.array) -> np.array: + U = np.zeros_like(A) + for i in range(len(A)): + for j in range(len(A)): + if i > j: + U[i, j] = A[i, j] + return U + + +def diag_strict_dominante(A) -> bool: + diag_sum = 0 + for i in range(len(A)): + diag_sum += A[i, i] + other_sum = 0 + for i in range(len(A)): + for j in range(len(A)): + if i != j: + other_sum += A[i, j] + return diag_sum > other_sum + + +def jacobi(A, b): + if not diag_strict_dominante(A): + raise Exception('A doit être à diagnonale strictement dominante') + + L = to_L(A) + U = to_U(A) + x0 = np.array([0,0,0]) + epsilon = 1e-6 + max_iter = 100_000 + + x = x0 + for k in range(max_iter): + x_new = np.diag(1 / np.diag(A)) @ ((L + U) @ x) + np.diag(1 / np.diag(A)) @ b + if np.linalg.norm(x_new - x, ord=2) < epsilon or np.linalg.norm(b - A @ x_new, ord=2) < epsilon: + break + x = x_new + + return x + + +def gauss_seidel(A, b): + x0 = np.array([0, 0, 0]) + + D = to_D(A) + L = to_L(A) + U = to_U(A) + epsilon = 1e-6 + done = False + + x = x0 + + while not done: + x_new = np.linalg.inv(D - L) @ U @ x + np.linalg.inv(D - L) @ b + + done: bool = np.linalg.norm(x_new - x, ord=2) < epsilon + x = x_new + + return x + + +def relaxation(A, b, omega=1.0, epsilon=1e-6, max_iter=100_000): + D = np.diag(np.diag(A)) + L = np.tril(A, k=-1) + U = np.triu(A, k=1) + x = np.zeros_like(b, dtype=float) + + # Pré-calculer (D - ωL)^(-1) une seule fois + inv_D_omega_L = np.linalg.inv(D - omega * L) + + if omega == 1: + return gauss_seidel(A, b) + + for _ in range(max_iter): + x_new = inv_D_omega_L @ ((1 - omega) * D @ x + omega * (U @ x + b)) + if np.linalg.norm(x_new - x, ord=2) < epsilon: + return x_new + x = x_new + + raise RuntimeError("La méthode de relaxation n'a pas convergé.") + + + +if __name__ == '__main__': + A = np.array([ + [8,4,1], + [1,6,-5], + [1,-2,-6] + ]) + + b = np.array([1,0,0]) + + D = to_D(A) + L = to_L(A) + U = to_U(A) + + res_jacobi = jacobi(A, b) + print(res_jacobi) + + res_gauss = gauss_seidel(A, b) + print(res_gauss) + + # je la commente car omega = 1 utilise la ^m fonction que gauss_seidel + #res_relaxation_1 = relaxation(A, b, 1) + #print(res_relaxation_1) + + res_relaxation_less_1 = relaxation(A, b, 0) + print(res_relaxation_less_1) + + res_relaxation_2 = relaxation(A, b, 2) + print(res_relaxation_2) + + print('fini') diff --git a/ex4.py b/ex4.py new file mode 100644 index 0000000..a576cb7 --- /dev/null +++ b/ex4.py @@ -0,0 +1,37 @@ +import numpy as np +import matplotlib.pyplot as plt + +from ex2 import gauss +from ex1 import resolve_up + +if __name__ == '__main__': + A = np.array([ + [-1, 5.6, 20.4], + [-1, 28, 26], + [-1, 8, 30] + ], dtype=np.float32) + + b = np.array([84.84, 361, 228.75], dtype=np.float32) + + A_res, b_res = gauss(A, b) + print(f'A: {A_res}') + print(f'b: {b_res}') + + x = resolve_up(A_res, b_res) + print(f'x: {x}') + + fig, ax = plt.subplots() + cercles = [(2.8, 10.2, 5.2), (14, 13, 7), (4, 15, 3.5)] + + for x, y, r in cercles: + cercle = plt.Circle((x, y), r, color='blue', fill=False) + ax.add_patch(cercle) + + # Ajuster les limites du graphique + ax.set_xlim(-4, 22) + ax.set_ylim(4, 21) + ax.set_aspect('equal') + + plt.title("Cercles tracés avec matplotlib") + plt.grid(True) + plt.show() diff --git a/ex5.py b/ex5.py new file mode 100644 index 0000000..eb6496e --- /dev/null +++ b/ex5.py @@ -0,0 +1,31 @@ + +class Point: + def __init__(self, x: float, y: float): + self.x = x + self.y = y + + +def interpolation_lagrange(points: list[Point], x: float) -> float: + interpolation: float = 0. + n = len(points) + for i in range(n): + term: float = 1. # 0.0 will always give 0, it's a product ! So lets put 1 + # Lagrange L(x) + for k in range(n): + if k != i: + term *= ((x - points[k].x) / (points[i].x - points[k].x)) + # interpolation with a piece of Lagrange (y_i) + interpolation += term * points[i].y + + return interpolation + + +def interpolation_newton(points: list[Point], x: float) -> float: + interpolation = 0. # alpha sum + n = len(points) + for i in range(n): + a_i = 1. + for k in range(n-1): + a_i *= (points[k+1].y - points[k].y) / (points[k+1].x - points[k].x) + + return interpolation diff --git a/requirements.txt b/requirements.txt new file mode 100644 index 0000000..71b334a Binary files /dev/null and b/requirements.txt differ