720 likes | 1.08k Views
Sistemas de lineales. Programación Numérica. Objetivos. Resolver un sistema de ecuaciones lineales acopladas Obtener la inversa de una matriz Calcular el determinante de una matriz Factorizar matrices. Solución de sistemas pequeños.
E N D
Sistemas de lineales Programación Numérica
Objetivos • Resolver un sistema de ecuaciones lineales acopladas • Obtener la inversa de una matriz • Calcular el determinante de una matriz • Factorizar matrices
Solución de sistemas pequeños La solución de un sistema lineal de 2 ecuaciones puede obtenerse fácilmente de forma gráfica a11x1 + a12x2 = b1 a21x1 + a22x2 = b2 Despejando x2 y graficando solución
Ejemplo Sea el sistema 3x1 + 2x2 = 18 –x1 + 2x2 = 2 Despejando x2 = 9 – 1.5x1 x2 = 1 + 0.5x1 (4, 3)
Tipos de soluciones No hay solución Infinidad de soluciones Sistema mal condicionado Ejemplo ½x1 + x2 = 1 ½x1 + x2 = ½ Ejemplo ½x1 + x2 = 1 x1 + 2x2 = 1 Ejemplo ½x1 + x2 = 1 0.501x1 + x2 = 0.9
Determinantes Solución de un sistema de 3x3 La matriz del sistema es: El determinante de A es:
Regla de Cramer Cada incógnita de un sistema de ecuaciones lineales algebraicas puede expresarse como la razón de dos determinantes con denominados D y con numerador obtenido a partir de D, al reemplazar la columna de coeficientes de la incógnita en cuestión por las constantes independientes b1, b2, …, bn. Por ejemplo x1 se calcula con:
Definición Un sistema lineal es un sistema de la forma: E1: a11x1 + a12x2 + a13x3 +...+ a1nxn = b1 E2: a21x1 + a22x2 + a23x3 +...+ a2nxn = b2 . En: an1x1 + an2x2 + an3x3 +...+ annxn = bn Donde x1, x2, x3.., xn son las incógnitas.
Operaciones básicas Se utilizarán tres operaciones para simplificar un sistema lineal: 1. La ecuación Ei puede multiplicarse por una constante l distinta de cero. (lEi) (Ei) 2. La ecuación Ej puede multiplicarse por una constante l distinta de cero y sumarse a otra ecuación Ei. (Ei +lEj) (Ei) 3. El orden de las ecuaciones Ei y Ei puede intercambiarse. (Ei) (Ei)
positivo Ejemplo E1: x1 + x2 + 3x4 = 4 E2: 2x1+ x2 – x3 + x4 = 1 E3: 3x1- x2 – x3 + 2x4 = -3 E4: -x1+2x2 +3x3 - x4 = 4 Las operaciones: (E2 –2E1)E2, (E3 –3E1) E3, (E4 + E1)E4 Transforman el sistema en:
Simplificar E1: x1 + x2 + 3x4 = 4 E2: - x2 – x3 - 5x4 = -7 E3: - 4x2 – x3 -7x4 =-15 E4: + 3x2 +3x3 + 2x4 = 8 Las operaciones: (E3 –4E2)E3, (E4 +3E2) E4 Transforman el sistema en:
Simplificar E1: x1 + x2 + 3x4 = 4 E2: - x2 – x3 - 5x4 = -7 E3: 3x3 +13x4 = 13 E4: -13x4 = -13 El sistema se puede resolver hacia atrás. x4 = 1 x3 = (13-13x4)/3 = 0 x2 = -(x3 + 5x4 -7) = 2 x1 = -x2 - 3x4 + 4 = -1
Comprobación E1: x1 + x2 + 3x4 = 4 E2: 2x1+ x2 – x3 + x4 = 1 E3: 3x1- x2 – x3 + 2x4 = -3 E4: -x1+2x2 +3x3 - x4 = 4 x1 = -1 x2 = 2 x3 = 0 x4 = 1 E1: (-1)+(2)+3(1)= 4 E2: 2(-1)+(2)–(0)+(1) = 1 E3: 3(-1)-(2)–(0)+2(1)= -3 E4: -(-1)+2(2)+3(0)-(1) = 4
Solución en Matlab %sistema de ecuaciones A = [1,1,0,3,4;2,1,-1,1,1;3,-1,-1,2,-3;-1,2,3,-1,4] x = [0 0 0 0]’; %operaciones %(E2 - 2E1) -> E2 A(2,:) = A(2,:) - 2*A(1,:); %(E3 - 3E1) -> E3 A(3,:) = A(3,:) - 3*A(1,:); %(E4 + E1) -> E4 A(4,:) = A(4,:) + A(1,:); %(E3 - 4E2) -> E3 A(3,:) = A(3,:) - 4*A(2,:); %(E4 + 3E2) -> E4 A(4,:) = A(4,:) + 3*A(2,:); x(4) = A(4,5)/A(4,4); x(3) = (A(3,5)-A(3,4)*x(4))/A(3,3); x(2) = (A(2,5)-A(2,3:4)*x(3:4))/A(2,2); x(1) = (A(1,5)-A(1,2:4)*x(2:4))/A(1,1);
Actividad Resolver x1 + x2 + x4 = 22x1 + x2 - x3 + x4= 1 4 x1 - x2 - 2x3 + 2 x4 = 0 3x1 - x2 - x3 + 2 x4 = -3
Matrices Una matriz es un arreglo rectangular de elementos de n renglones y m columnas.
Representación matricial El sistema de ecuaciones Se puede representar por una matriz ampliada
Ejemplo (E2 –2E1) E2, (E3 –3E1) E3, (E4 + E1) E4 (E3 –4E2) E3, (E3 +3E2) E4
Eliminación Gaussiana con sustitución hacia atrás Se forma la matriz aumentada Siempre que aii 0, realizamos la operación (Ej–(aji/aii) Ei) (Ej) Para i = 1, 2, 3, ..., n-1 y j = i + 1, i + 2, ..., n.
Sustitución hacia atrás En general
Algoritmo Gauss Entrada: número de incógnitas y de ecuaciones n, matriz aumentada (aij) donde 1<= i <= n y 1<= j <= n+1. Salida: solución x1, x2 ... xn o mensaje de que no hay solución 1. Para i = 1 hasta n-1 2. Sea p el entero más pequeño con i <= p <= n y api <> 0 3. Si no se encuentra p 4. Salida(“no existe solución”) 5. terminar
continuación 6. Si p<>i realice (Ei)<->(Ep) 7. Para j = i+1 hasta n 8. mij = aij/aii 9. Realice (Ej-mijEi) -> (Ej) 10. Si ann=0 entonces 11. Salida(“no existe solución”) 12. xn = an,n+1/ann 13. Para i = n – 1 hasta 1 14. xi = [ai,n+1-Sj=i+1 aijxj]/aii
Código en C void gauss(double A[][10],int n, double x[]){ //resuelve un sistema de ecuaciones Mx = b // A es la matriz aumentada (n x (n+1)) A = [M b] int i,j,k,p; //p es el pivote double temp,s,m; for(i = 0;i<n-1; i++){ p = i; while(p < n-1 && A[p][i] == 0) p = p + 1; if(p == n-1){ std::cout<<"No hay solución por pivoteo\n"; return; } if(p != i) for(j=0; j<n+1; j++){ temp = A[i][j]; A[i][j] = A[p][j]; A[p][j] = temp; } for(j = i+1; j<n; j++){ m = A[j][i]/A[i][i]; for(k = 0;k<n+1;k++) A[j][k] = A[j][k] - m*A[i][k]; } } if(A[n-1][n-1]==0){ std::cout<<"No solución\n"; return; } x[n-1]=A[n-1][n]/A[n-1][n-1]; for(i = n-2;i>=0; i--){ s = A[i][n]; for (j = i+1; j<n;j++) s = s - A[i][j]*x[j]; x[i] = s/A[i][i]; } }
Método de Gauss function x = ecuaciones(A) % resuelve un sistema de ecuaciones Mx = b % A es la matriz aumentada (n x (n+1)) A = [M b] [n dumy] = size(A);%n = número de renglones for i = 1 : n-1 p = i; while p < n && A(p,i) == 0; p = p + 1; end if p == n fprintf('No hay solución...'); return; end if p ~= i temp = A(i,:); A(i,:) = A(p,:); A(p,:) = temp; end for j = i+1 : n m = A(j,i)/A(i,i); for k = 1 : n+1 A(j,k) = A(j,k) - m*A(i,k); end end end
Método de Gauss if A(n,n)==0 fprintf('No hay solución...'); return; end x(n) = A(n,n+1)/A(n,n); for i = n-1 : -1: 1 s = A(i,n+1); for j = i+1 : n s = s - A(i,j)*x(j); end x(i) = s/A(i,i); end
Inversa La inversa de una matriz A de n x n es una matriz A-1 tal que A x A-1 = A-1 x A = In Una matriz que tiene inversa se le llama matriz nosingular. Una matriz que no tiene inversa es una matriz singular. Se cumple que: a. La inversa es única. b. La inversa de una matriz no singular es no singular y (A-1)-1 = A c. Si A y B son no singulares, (AB)-1 = B-1A-1
Obtención de la inversa Sea Bj la j-ésima columna de la matriz B de n x n. Si AB = C, entonces la j-esima columna de C está dada por:
Cont. Supongamos que A-1 existe y que A-1 = B = (bij). Entonces AB = I y Para encontrar B podemos resolver n ecuaciones lineales donde la j-ésima columna de la inversa es la solución del sistema lineal que en el lado derecho tiene la j-ésima columna de I.
Ejemplo Si B = A-1, entonces AB = I se tiene
Inversa en C void inversa(double B[][20],int n,double x[][20]){ double I[20][20]; double temp,m,s; int i,j,k,p; for(i = 0; i<20; i++) for(j = 0; j<20; j++) if(i==j) I[i][j] = 1; else I[i][j] = 0; double A[20][20]; for(i = 0; i<n; i++) for(j = 0; j<2*n; j++) if(j<n) A[i][j] = B[i][j]; else A[i][j] = I[i][j-n]; for(i = 0; i<20; i++) for(j = 0; j<20; j++) x[i][j] = 0; for(i = 0; i<n-1;i++){ p = i; while(p < n && A[p][i] == 0) p = p + 1; if( p == n){ std::cout << "Matriz singular...\n"; return; }
if( p != i) for( k = 0; k<2*n; k++){ temp = A[i][k]; A[i][k] = A[p][k]; A[p][k] = temp; } for(j = i+1;j<n; j++){ m = A[j][i]/A[i][i]; for( k = 0; k<2*n; k++) A[j][k] = A[j][k] - m*A[i][k]; } } if( A[n-1][n-1]==0){ std::cout << "Matriz singular...\n"; return; } for( k = 0; k<n ; k++){ x[n][k] = A[n-1][n+k]/A[n-1][n-1]; for( i = n-1; i>=0; i--){ s = A[i][n+k]; for( j = i+1 ; j<n; j++) s = s - A[i][j]*x[j][k]; x[i][k] = s/A[i][i]; } } }
Inversa en Matlab function x = inversa(B) [n dumy] = size(B); I = eye(n); A = [B,I]; x = zeros(n); for i = 1 : n-1 p = i; while p < n && A(p,i) == 0; p = p + 1; end if p == n fprintf('Matriz singular...'); return; end if p ~= i temp = A(i,:); A(i,:) = A(p,:); A(p,:) = temp; end for j = i+1 : n m = A(j,i)/A(i,i); for k = 1 : 2*n A(j,k) = A(j,k) - m*A(i,k); end end end
Inversa en Matlab if A(n,n)==0 fprintf('Matriz singular...'); return; end for k = 1: n x(n,k) = A(n,n+k)/A(n,n); for i = n-1 : -1: 1 s = A(i,n+k); for j = i+1 : n s = s - A(i,j)*x(j,k); end x(i,k) = s/A(i,i); end end
Comandos de Matlab para matrices inv(m) – invierte una matriz cuadrada m.También m^-1. / – divide, A/B es equivalente a B*inv(A) \ – divide, A\B es equivalente a inv(A)*B det(m) – calcula el determinante de m. trace(m) – suma los elementos de la diagonal de m.
Tarea 30 1. Determine cual de las siguientes matrices son no singulares y calcule su inversa. a.| 4 2 6| b.|1 2 0| c.|4 0 0| | 3 0 7| |2 1 -1| |0 0 0| |-2 -1 -3| |3 1 1| |0 0 3|d.|1 1 -1 1| e.|4 0 0 0| f.|2 0 1 2| |1 2 -4 -2| |6 7 0 0| |1 1 0 2| |2 1 1 5| |9 11 1 0| |2 -1 3 1| |-1 0 -2 -4| |5 4 1 1| |3 -1 4 3|
Álgebra lineal Definición: dos matrices son iguales si tienen el mismo tamaño n x m y si aij = bij para todo i = 1, 2, .., n y j = 1, 2, ..., m. Definición: la suma de dos matrices A y B del mismo tamaño n x m es igual a una matriz C donde cij = aij + bij para todo i = 1, 2, .., n y j = 1, 2, ..., m. Definición: la multiplicación de una matriz A de tamaño n x m por una escalar l es una matriz de tamaño n x m cuyos elementos son laij para todo i = 1, 2, .., n y j = 1, 2, ..., m.
Propiedades de las matrices Si A, B y C son matrices de tamaño n x m, y l y m son números reales, se cumplen las siguientes propiedades: a. A + B = B + A e. l(A + B) = lB + lA b. (A + B) + C = A + (B + C) f. (l + m)A = lA + mA c. A + 0 = 0 + A = A g. l (mA) = (lm)A d. A + (-A) = -A + A = 0 h. 1A = A
Multiplicación de matrices Sea A una matriz de n x m y B una matriz de m x p. El producto de la matriz A por la matriz B se define como la matriz C de n x p donde los elementos de C se calculan por:
Definiciones Una matriz cuadrada es la que tiene igual número de renglones que de columnas. Una matriz diagonal es una matriz cuadrada con D =(dij) con dij = 0 simpre que i<>j. La matriz identidad de orden n, In = (dij) es una matriz diagonal con dij = 1 si i = j y dij =0 si i<>j. Una matriz triangular superior de n x n U = (uij) tiene, para toda j = 1, 2, ..., n, los elementos uij = 0, para cada i = j +1, j + 2, ..., n Una matriz triangular inferior de n x n L = (lij) tiene, para toda j = 1, 2, ..., n, los elementos lij = 0, para cada i = 1, 2, ..., j – 1
Traspuesta La traspuesta de una matriz A de n x m es una matriz At tal que la i_ésima columna de At es la misma que el i-ésimo renglón de A. Se cumple que a. (At)t = A b. (A + B)t = At + Bt c. (AB)t = Bt At d. Si A-1 existe, (A-1)t = (At)-1
Determinante Si A es una matriz de 1 x 1, entonces det(A) = x. Si a es de orden mayor que 1, calcule el determinante de A como sigue: Escoja cualquier fila o columna. para cada elemento A[i, j] en esa fila o columna, forme el producto (-1)(i+ j)*A[i, j]*det(menor(A[i, j])) Donde det(menor(A[i, j])) es el determinante del menor de A[i, j]. det(A) = suma de todos los productos para la columna o fila seleccionada.
Propiedades a. Si un renglón o columna tiene solo ceros, el determinante es cero. b. Si se intercambian 2 renglones o columnas, el signo del determinante cambia c. Si dos columnas o renglones son iguales, el determinante es cero. d. Si se multiplica un renglón o columna por un numero real el determinante se multiplica por ese número real. e. Si se suma un múltiplo de un renglón o columna a otro renglón o columna, el determinante no se altera. f. El determinante de un producto de matrices es igual al producto de los determinantes de cada una. g. El determinante de la inversa es el inverso del determinante de la matriz original.
Cálculo eficiente del determinante Convertir la matriz en una matriz triangular superior y calcular el determinante de ésta. El determinante de una matriz triangular superior es igual al producto de los elementos de la diagonal principal. La matriz triangular se obtiene mediante la eliminación gaussiana.
Determinante en Matlab for j = i+1 : n m = A(j,i)/A(i,i); for k = 1 : n A(j,k) = A(j,k) - m*A(i,k); end end end if A(n,n)==0 x = 0; return; end p = 1; for i = 1:n p = p*A(i,i); end x = p; function x = determinante(A) [n dumy] = size(A); for i = 1 : n-1 p = i; while p < n && A(p,i) == 0; p = p + 1; end if p == n x = 0; return; end if p ~= i temp = A(i,:); A(i,:) = A(p,:); A(p,:) = temp; end
Descomposición LU La descomposición LU transforma una matriz a en el producto de dos matrices triangulares A = LU (1) Donde L es triangular inferior y U es triangular superior. Un sistema Ax = b se puede escribir como LUx = b (2) Si hacemos Ux = d (3) Obtenemos Ld = b (4) Con (4) se obtiene z y después con (3) se obtiene x. En los casos en que se tengan varios sistemas con la misma matriz de coeficientes es más eficiente aplicar este método que resolver cada uno mediante Gauss.
Supongamos un sistema de 3x3 En el primer paso de eliminación gaussiana se multiplica el renglón 1 por f 21 = a21/a11 y se le resta al segundo para eliminar a21. Luego se multiplica el renglón 1 por f 31 = a31/a11 y se le resta al tercer para eliminar a31. El paso final es multiplicar el segundo renglón por f 32 = a’32/a’22 y restar el resultado al tercer renglón. Los valores de f 21, f 31 y f 32 se pueden retener en a21, a31 y a32 . Después de la eliminación
Esta matriz representa un almacenamiento eficiente de la descomposición LU de A. A = L U donde
Ejemplo Los valores de f 21, f 31 son f 21 = 0.1/3 = 0.033333 y f 31 = 0.3/3 = 0.1 El valore de f 32 es f 32 = 0.19/7.003333 = 0.027130