83 lines
2.4 KiB
Python
83 lines
2.4 KiB
Python
import numpy as np
|
|
from ex1 import resolve_up, resolve_down
|
|
|
|
|
|
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
|
|
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ée et le vecteur b modifié
|
|
|
|
|
|
|
|
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
|
|
|
|
|
|
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_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: {L}")
|
|
x = resolve_down(L, b)
|
|
print(f'x: {x}')
|