From 128e07d420f3eacdd450d35da7bfba792762b378 Mon Sep 17 00:00:00 2001 From: Namu Date: Sun, 19 Oct 2025 19:03:43 +0200 Subject: [PATCH] Fix: Correct jacobi, gauss, cholesky --- ei_blanc.py | 0 ex1.py | 22 ++++++++++--- ex2.py | 72 +++++++++++++++++++++++++++------------- ex3.py | 95 +++++++++++++++++++++-------------------------------- ex5.py | 40 +++++++++++++++------- 5 files changed, 133 insertions(+), 96 deletions(-) create mode 100644 ei_blanc.py diff --git a/ei_blanc.py b/ei_blanc.py new file mode 100644 index 0000000..e69de29 diff --git a/ex1.py b/ex1.py index d22a1ac..68a2d13 100644 --- a/ex1.py +++ b/ex1.py @@ -1,22 +1,34 @@ import numpy as np -def resolve_up(A: np.array, b: np.array) -> np.array: +def resolve_up(A, b) -> list[float]: + """ + + :param A: Une matrice + :param b: Un vecteur + :return: + """ 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)): + 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] + x[i] = (b[i] - right_of_diagonal) / A[i, i] return x -def resolve_down(A: np.array, b: np.array) -> np.array: +def resolve_down(A, b) -> list[float]: + """ + + :param A: Une matrice + :param b: Un vecteur + :return: + """ 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] + sum += A[i, k] * x[k] x[i] = (b[i] - sum) / A[i, i] return x diff --git a/ex2.py b/ex2.py index d3b0dc1..27b7fe2 100644 --- a/ex2.py +++ b/ex2.py @@ -1,11 +1,29 @@ -""" -Ici, on implémente gauss -""" import numpy as np -from ex1 import resolve_up +from ex1 import resolve_up, resolve_down -def gauss(A: np.array, b: np.array) -> (np.array, np.array): +def est_symetrique_definie_positive(A) -> bool: + """ + Vérifie que A = A^t et que toutes les valeurs propres soit positive + :param A: + :return: + """ + + # Symétrie (avec tolérance pour les flottants) + sym = np.allclose(A, A.T) + # Définie positive : toutes les valeurs propres > 0 + valeurs_propres = np.linalg.eigvals(A) + definie_positive = np.all(valeurs_propres > 0) + return sym and definie_positive + + +def gauss(A, b): + """ + Pas de condition, mais pas optimal non plus. + :param A: Une matrice + :param b: Un vecteur + :return: A triangularisée et b modifié en conséquence. + """ n = len(A) for k in range(n): # k = pivot for i in range(k+1, n): # i = ligne en dessous du pivot @@ -13,24 +31,25 @@ def gauss(A: np.array, b: np.array) -> (np.array, np.array): 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é + return A, b # matrice triangularisée et le vecteur 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] +def cholesky(A): + if not est_symetrique_definie_positive(A): + raise ValueError('Matrice non symétrique définie positive') + + m = A.shape[0] # Taille de la matrice (m x m) + L = np.zeros_like(A, dtype=np.float32) # Matrice L initialisée à zéro + + for j in range(m): # Parcourt chaque colonne j de 0 à m-1 + for i in range(j, m): # Parcourt chaque ligne i de j à m-1 + s = sum(L[i, k] * L[j, k] for k in range(j)) # Somme des produits L[i,k] * L[j,k] pour k < j + + if i == j: # Élément diagonal + L[i, j] = np.sqrt(A[i, i] - s) # Calcul de la racine carrée pour la diagonale + else: # Élément hors diagonale + L[i, j] = (A[i, j] - s) / L[j, j] # Calcul des éléments hors diagonale return L @@ -47,8 +66,17 @@ if __name__ == '__main__': A_triangle, b_gauss = gauss(A, b) print(f"A triangulsarisée:\n {A_triangle}") - x = resolve_up(A, b_gauss) + x = resolve_up(A_triangle, b_gauss) print(f"x={x}") + A = np.array([ + [4, 1, 1, 1], + [1, 5, 2, 1], + [1, 2, 6, 2], + [1, 1, 2, 7] + ]) + L = cholesky(A) - print(f"L\n {L}") + print(f"L: {L}") + x = resolve_down(L, b) + print(f'x: {x}') diff --git a/ex3.py b/ex3.py index e1f6623..e3fe251 100644 --- a/ex3.py +++ b/ex3.py @@ -5,81 +5,62 @@ import numpy as np # 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 + """ + Pour chaque ligne, il faut que la somme des coefs non-diagonaux soit inférieure au coef diagonal + TOUT EST EN VALEUR ABSOLUE. + :param A: + :return: + """ + A = np.array(A) + n = len(A) + for i in range(n): + diag = abs(A[i,i]) + row_sum = np.sum(np.abs(A[i,:])) - diag # On fait la somme de toute la ligne puis on soustrait le coef diagonal + if diag <= row_sum: + return False + return True -def jacobi(A, b): +def jacobi(A, b, epsilon=1e-6, max_iter=100_000): if not diag_strict_dominante(A): - raise Exception('A doit être à diagnonale strictement dominante') + raise ValueError("A doit être à diagonale strictement dominante") - L = to_L(A) - U = to_U(A) - x0 = np.array([0,0,0]) - epsilon = 1e-6 - max_iter = 100_000 + x = np.zeros_like(b, dtype=float) + L = np.tril(A, k=-1) # Partie triangulaire inférieure (sans diagonale) + U = np.triu(A, k=1) # Partie triangulaire supérieure (sans diagonale) - 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: + for _ in range(max_iter): + x_new = np.diag(1 / np.diag(A)) @ (b - (L + U) @ x) + if np.linalg.norm(x_new - x, ord=2) < epsilon: break x = x_new - return x + return x_new def gauss_seidel(A, b): x0 = np.array([0, 0, 0]) - D = to_D(A) - L = to_L(A) - U = to_U(A) + D = np.diag(np.diag(A)) + L = np.tril(A, k=-1) + U = np.triu(A, k=1) epsilon = 1e-6 - done = False + max_iter = 100_000 + + # Pré-calcul de l'inverse de (D - L) pour éviter de le recalculer à chaque itération + inv_D_minus_L = np.linalg.inv(D - L) x = x0 - while not done: - x_new = np.linalg.inv(D - L) @ U @ x + np.linalg.inv(D - L) @ b + for _ in range(max_iter): + x_new = inv_D_minus_L @ (-U @ x) + inv_D_minus_L @ b - done: bool = np.linalg.norm(x_new - x, ord=2) < epsilon + if np.linalg.norm(x_new - x, ord=2) < epsilon: + return x_new x = x_new - return x + raise RuntimeError('Gauss-Seidel') def relaxation(A, b, omega=1.0, epsilon=1e-6, max_iter=100_000): @@ -113,9 +94,9 @@ if __name__ == '__main__': b = np.array([1,0,0]) - D = to_D(A) - L = to_L(A) - U = to_U(A) + D = np.diag(A) + L = np.tril(A) + U = np.triu(A) res_jacobi = jacobi(A, b) print(res_jacobi) diff --git a/ex5.py b/ex5.py index 4bb0dfc..8e83bdc 100644 --- a/ex5.py +++ b/ex5.py @@ -4,6 +4,7 @@ import numpy as np dich_steps = [] new_raph_steps = [] + def f(x) -> float: return np.exp(x) - 2 * np.cos(x) @@ -26,18 +27,33 @@ def dichotomie(x_0, a, b, epsilon, n=0) -> float: def f_prime(x) -> float: - return 1 + np.sin(x) + return np.exp(x) + 2 * np.sin(x) -def newton_raphson(x_0, epsilon) -> float: - x_1 = x_0 - (f(x_0) / f_prime(x_0)) - new_raph_steps.append(x_1) +def newton_raphson(x_0, epsilon, max_iter=1_000) -> float | None: + """ + Condition : connaître la dérivée de f(x) + :param x_0: + :param epsilon: + :param max_iter: + :return: + """ + x = x_0 - if np.abs(x_1 - x_0) < epsilon: - return x_1 - elif np.abs(f(x_0)) < epsilon: - return x_1 - return newton_raphson(x_1, epsilon) + for _ in range(max_iter): + x_prev = x + fx = f(x) + fpx = f_prime(x) + if np.abs(fx) < epsilon: + return x + elif fpx == 0: + return None + x = x - (fx / fpx) + new_raph_steps.append(x) + if np.abs(x - x_prev) < epsilon: + return x + + return None def ex1() -> None: @@ -71,10 +87,10 @@ def ex1() -> None: plt.show() """ - La dichotomie à moins d'étape que NR, mais NR à l'air plus rapide, - comme Flash McQueen + La NR à moins d'étape que la dichotomie, NR à l'air rapide, + comme Flash McQueen (Kachow) """ if __name__ == '__main__': - ex1() \ No newline at end of file + ex1()