Untitled

895 days ago by pub

def mP(M): """Contruimos una matriz P de pivoteo para el método de Doolittle""" # Creamos matriz identidad m =M.nrows() P=identity_matrix(m) # Contruimos la matriz de transfomación de tal forma que coloque # en la diagonal principal de M el mayor elemento de la primera # columna de las sucesivas submatrices principales for j in range(m): fila = max(range(j, m), key=lambda i: M[i,j]) if j != fila: # intercambiamos filas P[j], P[fila] = P[fila], P[j] return P def mLU(A): """Efectuamos una Descomposición LU de A (A matiz cuadrada) en PA = LU. La función devuelve P, L y U.""" n = A.nrows() # Creamos matrices LU con ceros L = matrix(QQ,n,n) U = matrix(QQ,n,n) # Creamos la matriz P de pivoteo obtenemos la matriz PA P = mP(A) PA = P*A # Efectuamos la decomposición LU for s in range(n): # Todos los elementos de la diagonal de L son unitarios L[s,s] = 1 #construimos la fila s de U for j in range(s,n): U[s,j] = PA[s,j] - sum(U[k,j] * L[s,k] for k in range(s)) #construimos la columna j de L for i in range(s,n): L[i,s] = (PA[i,s] - sum(U[k,s] * L[i,k] for k in range(s))) / U[s,s] return (P, L, U) # Metodo de sustitucion regresiva o subida para # sistemas de ecuaciones lineales triangulares superiores def subida(A, b): # comprobamos validez del sistema m,n = A.nrows(),A.ncols() tam = b.nrows() if ((m <> n) or (tam <> m)): print 'Las dimensiones del sistema no son válidas' return True if A.det()==0: print 'El sistema no tiene solución única o es incompatible' return True # Vector de soluciones x = matrix(QQ,n,1) x[n-1] = b[n-1]/A[n-1,n-1] for k in range(n-2,-1,-1): suma = sum(x[i]*A[k,i] for i in range(n-1,k,-1)) x[k] = (b[k]-suma)/A[k,k] return x # Metodo de sustitucion progresiva o bajada para # sistemas de ecuaciones lineales triangulares inferiores def bajada(A, b): # comprobamos validez del sistema m,n = A.nrows(),A.ncols() tam = b.nrows() if ((m <> n) or (tam <> m)): print 'Las dimensiones del sistema no son válidas' return True if A.det()==0: print 'El sistema no tiene solución única o es incompatible' return True # Vector de soluciones x = matrix(QQ,n,1) x[0] = b[0]/A[0,0] for k in range(1,n): suma = sum(x[i]*A[k,i] for i in range(0,k)) x[k] = (b[k]-suma)/A[k,k] return x 
       
b=matrix(QQ,4,1,[1,1,1,1]) A=matrix(QQ,4,4,[5,2,1,1,1,4,2,1,1,1,3,1,1,2,3,7]) P,L,U= mLU(A) show(html("$"+latex(P)+latex(A)+"="+latex(L)+latex(U)+"$")) #caculamos Ly=Pb y=bajada(L,P*b) show(html("$\\text{y = }" +latex(y.n()) + "$")) #calulamos Ux=y x=subida(U,y) show(html("$\\text{x = }" +latex(x.n()) + "$")) #solución SAGE sol=A\b show(sol.n()) show("error relativo: " + ((sol-x).norm(oo)/sol.norm(oo)).str())