Referencia: Chapra 10.2 p287, pdf312
Se plantea resolver el sistema de ecuaciones usando matrices triangulares
A = \begin{pmatrix} 3 & -0.1 & -0.2 \\ 0.1 & 7 & -0.3 \\0.3 & -0.2 & 10 \end{pmatrix} B = [7.85,-19.3,71.4]Para encontrar la solución al sistema de ecuaciones, se plantea resolver:
– sustitución hacia adelante de LY=B
– sustitución hacia atras para UX=Y
Algoritmo en Python
Al algoritmo de la sección anterior se añade los procedimientos correspondientes con los que se obtiene la solución:
Matriz aumentada [[ 3. -0.1 -0.2 7.85] [ 0.1 7. -0.3 -19.3 ] [ 0.3 -0.2 10. 71.4 ]] Pivoteo parcial: Pivoteo por filas NO requerido matriz L: [[ 1. 0. 0. ] [ 0.03333333 1. 0. ] [ 0.1 -0.02712994 1. ]] Matriz U: [[ 3. -0.1 -0.2 ] [ 0. 7.00333333 -0.29333333] [ 0. 0. 10.01204188]] B1 : [[ 7.85] [-19.3 ] [ 71.4 ]] Y Sustitución hacia adelante [[ 7.85 ] [-19.56166667] [ 70.08429319]] X Sustitución hacia atras [ 3. -2.5 7. ] >>>
Instrucciones en Python
# Solución usando Matrices L y U # continuación de ejercicio anterior import numpy as np def pivoteafila(A,B,vertabla=False): ''' Pivotea parcial por filas Si hay ceros en diagonal es matriz singular, Tarea: Revisar si diagonal tiene ceros ''' A = np.array(A,dtype=float) B = np.array(B,dtype=float) # Matriz aumentada nB = len(np.shape(B)) if nB == 1: B = np.transpose([B]) AB = np.concatenate((A,B),axis=1) if vertabla==True: print('Matriz aumentada') print(AB) print('Pivoteo parcial:') # Pivoteo por filas AB tamano = np.shape(AB) n = tamano[0] m = tamano[1] # Para cada fila en AB pivoteado = 0 for i in range(0,n-1,1): # columna desde diagonal i en adelante columna = np.abs(AB[i:,i]) dondemax = np.argmax(columna) # dondemax no es en diagonal if (dondemax != 0): # intercambia filas temporal = np.copy(AB[i,:]) AB[i,:] = AB[dondemax+i,:] AB[dondemax+i,:] = temporal pivoteado = pivoteado + 1 if vertabla==True: print(' ',pivoteado, 'intercambiar filas: ',i,'y', dondemax+i) if vertabla==True: if pivoteado==0: print(' Pivoteo por filas NO requerido') else: print(AB) return(AB) # PROGRAMA ------------ # INGRESO A = [[ 3. , -0.1, -0.2], [ 0.1, 7. , -0.3], [ 0.3, -0.2, 10. ]] B = [7.85,-19.3,71.4] # PROCEDIMIENTO AB = pivoteafila(A,B,vertabla=True) tamano = np.shape(AB) n = tamano[0] m = tamano[1] AB0 = np.copy(AB) B1 = np.copy(AB[:,m-1]) L = np.identity(n,dtype=float) # Inicializa L # eliminacion hacia adelante for i in range(0,n-1,1): pivote = AB[i,i] adelante = i+1 for k in range(adelante,n,1): factor = AB[k,i]/pivote AB[k,:] = AB[k,:] - AB[i,:]*factor L[k,i] = factor # llena L U = np.copy(AB[:,:n]) # Resolver LY = B B1 = np.transpose([B1]) AB =np.concatenate((L,B1),axis=1) # sustitución hacia adelante Y = np.zeros(n,dtype=float) Y[0] = AB[0,n] for i in range(1,n,1): suma = 0 for j in range(0,i,1): suma = suma + AB[i,j]*Y[j] b = AB[i,n] Y[i] = (b-suma)/AB[i,i] Y = np.transpose([Y]) # Resolver UX = Y AB =np.concatenate((U,Y),axis=1) # sustitución hacia atrás ultfila = n-1 ultcolumna = m-1 X = np.zeros(n,dtype=float) for i in range(ultfila,0-1,-1): suma = 0 for j in range(i+1,ultcolumna,1): suma = suma + AB[i,j]*X[j] b = AB[i,ultcolumna] X[i] = (b-suma)/AB[i,i] # SALIDA print('matriz L: ') print(L) print('Matriz U: ') print(U) print('B1 :') print(B1) print('Y Sustitución hacia adelante') print(Y) print('X Sustitución hacia atras') print(X)