Factorisation LU Python

Discutez d'informatique ici !
acoustica
Membre Irrationnel
Messages: 1043
Enregistré le: 08 Juil 2008, 11:00

Factorisation LU Python

par acoustica » 12 Aoû 2012, 14:46

Bonjour à tous,

EDIT : PROBLEME RESOLU. Le programme avec espaces en début de ligne s'obtient en faisant "éditer le message"

Source du problème : l'inversion de matrices par la méthode des cofacteurs devient ingérable à partir d'une matrice trop grande, les erreurs d'arrondis deviennent énormes après un calcul de l'ordre de n!.


J'ai un petit soucis avec un programme Python, il s'agit d'une décomposition LU (avec tout ce qui doit être programmé autour, à savoir le produit de matrices, l'inversion, etc...). Bref, le programme semble marcher au début, mais ça a commencé à foirer quand j'ai voulu aussi avoir la matrice "low", puisqu'il faut à ce moment-là inverser (ce qui réclame plein de procédures parallèles).

Voici un exemple de factorisation LU :

Image

Et voici l'algorithme :

Image


Le problème que je rencontre est que ça indique "non inversible" à chaque fois. Avant la décomposition "up" marchait, mais je ne sais pas, j'ai dû toucher à un truc qu'il ne fallait pas toucher...
Quand on retire tout ce qui a trait à l'inversibilité (pas besoin pour calculer "up", on n'a besoin d'inverser que pour "low", là ça marche, en supprimant tout les "def Det", "def com"...) là ça marche et on obtient bien une matrice triangulaire supérieure. C'est pour l'inférieure que là il y a besoin d'inverser et que systématiquement ça indique non inversible. Pourtant je ne vois vraiment pas d'où vient le problème....

Ci-dessous le programme, mais le forum a la fâcheuse particularité de ne pas tenir compte des espace en début de ligne. Donc pour un éventuel copié-collé, les espaces sont dans le programme en lien :


n=int(input("dimension de la matrice carrée ? "))
import math

A=[]
T=[]


def produitmatrices(C,D):
P=[]
p=0

for i in range(n):
P.append([])

for i in range(n):
for j in range(n):
P[j].append(0)

for i in range(n):
for j in range(n):
for m in range(n):
p=p+float(C[i][m])*float(D[m][j])
P[i][j]=p
p=0
return P

def Mssij(M,i,j):
l=len(M)
Mwthij=[]
stMwthij=[]
for y in range(l):
if y!=i:
for z in range(l):
if z!=j:
Mwthij.append(M[y][z])
for u in range(0,len(Mwthij),l-1):
stMwthij.append(Mwthij[u:u+l-1])
return stMwthij

def Det(T):
if len(T)==1:
return T[0][0]
if len(T)==2:
d=T[0][0]*T[1][1]-T[0][1]*T[1][0]
return d
else:
s=0
for j in range(len(T)):
N=Mssij(T,j,0)
if j//2==j/2:
s=s+T[j][0]*Det(N)
else:
s=s-T[j][0]*Det(N)
return s

def com(W):
w=len(W)
X=[]

for i in range(w):
X.append([])

for g in range(w):
for h in range(w):
v=Det(Mssij(W,g,h)*((-1)**(g+h)))
X[g].append(v)
return X


def KfoisMat(E,f):
for o in range(len(E)):
for r in range(len(E)):
E[o][r]=f*E[o][r]
return E

def transpose(O):
X=[]
for i in range(len(O)):
X.append([])
for o in range(len(O)):
for r in range(len(O)):
X[r].append(O[o][r])
return X

def inversematrice(M):
if Det(M)==0:
print("PAS INVERSIBLE")
else :
return(transpose(KfoisMat(com(M),1/Det(M))))


for i in range(n):
A.append([])

for i in range(n):
for j in range(n):
A[j].append(0)

for i in range(n-1):
T.append([])

for i in range(n-1):
for j in range(n):
T[i].append([])

for k in range(n-1):
for i in range(n):
for j in range(n):
if i==j:
T[k][i].append(1)
else:
T[k][i].append(0)

for i in range(1,n+1):
for j in range(1,n+1):
print("a",i,j,"? :")
A[i-1][j-1]=float(input(" "))

S=A

for k in range(n-1):
if A[k][k]!=0:
pivot=float(A[k][k])
else:
for b in range(n):
if A[b][k]!=0:
pivot=float(A[b][k])
break
for i in range(k+1,n):
T[k][i][k]=-float(A[i][k])/float(pivot)
A=produitmatrices(T[k],A)

for dep in range(n):
print(S[dep])

print("=")
for k in range(n-1):
T[k]=inversematrice(T[k])
for k in range(1,n-1):
T[k]=produitmatrices(T[k-1],T[k])

for q in range(n):
print(T[k][q])

print("")
print("&*")
print("")
for q in range(n):
print(A[q])



 

Retourner vers ϟ Informatique

Qui est en ligne

Utilisateurs parcourant ce forum : Aucun utilisateur enregistré et 2 invités

Tu pars déja ?



Fais toi aider gratuitement sur Maths-forum !

Créé un compte en 1 minute et pose ta question dans le forum ;-)
Inscription gratuite

Identification

Pas encore inscrit ?

Ou identifiez-vous :

Inscription gratuite