Método de eliminación de Gauss con pivotamiento

1. Método de eliminación de Gauss con pivotamiento Métodos Matemáticos de Especialidad (Mecánica-Máquinas) Madrid, 10 de octubre de 2005 Javier García de Jalón ETSII - Departamento de Matemática Aplicadaa la Ingeniería Industrial

2. Introducción • La resolución de sistemas de ecuaciones lineales es uno de los problemas matemáticos más importantes en Ingeniería • Hasta la llegada de los computadores digitales (segunda mitad del s. XX) la capacidad de resolver sistemas de ecuaciones estaba muy limitada, no por la complejidad del problema, sino por el número de operaciones aritméticas • Ahora se puede resolver con un PC un sistema 1000×1000 en menos de 1 seg. • Con programas especiales que aprovechan la estructura de la matriz se pueden resolver de forma rutinaria con PCs, sistemas de decenas ó cientos de miles de ecuaciones lineales • Muchos métodos matemáticos (cálculo de valores y vectores propios, integración de ecuaciones diferenciales, optimización, ...) se reducen a la resolución repetida de sistemas de ecuaciones lineales • La resolución de sistemas de ecuaciones lineales tiene además un importante valor didáctico • para los métodos numéricos en general • para la programación de ordenadores

3. b A(:,1)x1 A(:,3)x3 A(:,1) A(:,3) A(:,2) A(:,2)x2 Interpretaciones del sistema Ax=b • Intersección de "planos" • La solución es el punto que satisface las ecuaciones de todos los hiperplanos, es decir, su intersección. • La intersección puede ser unpunto, un subespacio de dimensiónn−r, o el conjunto vacío. • Combinación lineal de vectores • El vector término independiente b es una combinación lineal de las columnas de A, cuyoscoeficientes son los valores de x:

4. m=3, n=2, r=2 m=3, n=2, r=2 m=2, n=2, r=2 bIm(A) bIm(A) m=1, n=2, r=1 m=2, n=2, r=1 m=2, n=2, r=1 bIm(A) bIm(A) Algunos casos posibles del sistema Ax=b • Las siguientes figuras muestran algunas de las posibles situaciones en un sistema de dos ecuaciones con dos incógnitas

5. El método de eliminación de Gauss (1/4) • Se hace un cero en la posición (2,1) multiplicando por la matriz P21 donde • Se obtienen ceros en el resto de la 1ª columna pre-multiplicando por las matrices P31 y P41 : • Obsérvese la forma en que se multiplican las matrices Pi1 • Si la matriz A es simétrica, la submatriz de P1A en la que hay que seguir haciendo ceros también lo es.

6. El método de eliminación de Gauss (2/4) • De la misma forma se hacen ceros, debajo de los elementos de la diagonal, en las columnas 2 y 3: • Para llegar a este resultado se ha supuesto que los elementos que aparecen sobre la diagonal −los pivots− son distintos de cero, pues aparecen en los denominadores de los factoes mji. • Las matrices elementales Pi tienen propiedades especiales.

7. Las matrices Pj obtenidas operan de una forma particular: Por una parte, el producto P3P2 no ofrece ninguna forma especial: La inversa de las matrices Pj se obtiene cambiando el signo a los factores mij Finalmente, el producto de las inversas de P3P2−en orden opuesto− se realiza por simple “superposición”: El método de eliminación de Gauss (3/4)

8. El método de eliminación de Gauss (4/4) • Finalmente se obtiene la matriz Atriangularizada (o en forma de escalera, como se verá al hablar del pivotamiento) • Pre-multiplicando por las inversas de las matrices Pi y teniendo en cuenta la forma que adopta dicho producto de inversas: • Conclusión: • El método de eliminación de Gauss equivale a descomponer la matriz A en el producto de una matriz triangular inferior L con “unos” en la diagonal, por una matriz triangular superior U: • La factorización LU tiene una gran importancia en álgebra lineal numérica.

9. k j k k = i i 0 k i k k * * * = 0 Programa de Gauss básico n×n (1/2) • Algoritmo básico con Matlab • La triangularización se basa en hacer un cero en (i,k) combinando las filas i y k: • Para calcular las incógnitas mediante la vuelta atrás, el proceso es el siguiente: for k=1:n-1 for i=k+1:n m=A(i,k)/A(k,k); for j=k+1:n A(i,j)=A(i,j)-m*A(k,j); end b(i)=b(i)-m*b(k); end end x(n)=b(n)/A(n,n); for k=n-1:-1:1 s=0; for i=k+1:n s=s+A(k,i)*x(i); end x(k)=(b(k)-s)/A(k,k); end

10. k 0 k k 1 * * 1 = k i k k * * * = 0 Programa de Gauss básico n×n (2/2) • Programa para resolver un sistema en la forma LUx=Ly=b % factorización A=LU for k=1:n-1 for i=k+1:n m=A(i,k)/A(k,k); for j=k+1:n A(i,j)=A(i,j)-m*A(k,j); end % se almacena m en la posición del 0 A(i,k)=m; end end % Resolución del sistema Ly=b y(1)=b(1); for k=2:n s=0; for j=1:k-1 s=s+A(k,j)*y(j); end y(k)=b(k)-s; end % Resolución del sistema Ux=y x(n)=y(n)/A(n,n); for k=n-1:-1:1 s=0; for j=k+1:n s=s+A(k,j)*x(j); end x(k)=(y(k)-s)/A(k,k); end

11. Método de Gauss con pivotamiento (1/7) • En el proceso de eliminación de Gauss es necesario dividir por los elementos que aparecen en la diagonal, que se denominan “pivots” • En ocasiones puede aparecer un cero (o un elemento de valor muy pequeño) en la posición de un pívot: en ese caso el proceso no puede continuar. • Si debajo del pívot hay elementos 0 se pueden permutar las filas correspondientes y proseguir la eliminación (la solución del sistema no varía): • La permutación de filas puede hacerse pre-multiplicando por una matriz P: • La matriz P puede construirse a partir de un vector p con el orden de los pivots

12. Método de Gauss con pivotamiento (2/7) • Puede suceder que algunas columnas carezcan de pívot válido • Por ejemplo, cuando al hacer ceros en la columna 2 los elementos correspondientes de la columna 3 se han hecho también cero: el proceso normal falla y hay dos posibilidades: • La eliminación puede proseguir en las columnas posteriores y llegar hasta un último pívot en la última fila (x3 es una variable libre o independiente, que puede tomar un valor arbitrario). • Le eliminación se detiene porque debajo del último pívot encontrado todas las filas restantes se han hecho cero (unafila dependiente de las demás, y dos variables libres: x3 y x5).

13. Método de Gauss con pivotamiento (3/7) • Método de Gauss con pivotamiento por columnas: • Se busca el mayor elemento en valor absoluto de la primera columna. Este elemento servirá como pívot para hacer ceros en la primera columna. • El segundo pívot es el máximo elemento en valor absoluto de la 2ª columna, sin contar el elemento de la fila del primer pívot. Con este pívot se hacen ceros en los elementos de la 2ª columna que no pertenecen a una fila de un pívot anterior. • De modo análogo se calculan los pivots en las restantes columnas y se hacen cero los elementos de cada columna en filas en las que no hay pivots. • En el sistema final es fácil hallar las incógnitas xn, ..., x1, a partir de las filas en que han aparecido el último pívot, hasta la del primero. • En la práctica, no es necesario intercambiar filas: basta conocer las filas dónde han aparecido los sucesivos pivots:

14. Método de Gauss con pivotamiento (4/7) • Posibles dificultades del pivotamiento por columnas: • Si una ecuación del sistema (o una fila de la matriz) se multiplica por un número muy grande, la solución del sistema no varía, pero puede alterar el orden en el que se eligen los pivots, y con dicho orden la precisión de la solución calculada. • Por ello es conveniente que todas las filas de la matriz sean vectores que tengan aproximadamente la misma norma, o bien que, al elegir los pivots se comparen los elementos de una columna dividido cada uno de ellos por la norma de su fila. • La norma más fácil de utilizar es la norma sub-infinito: ||A||. • A este proceso se le denomina escalamiento o scaling por filas. • Ver programas de Matlab: • pivotcol1.m: Selección directa de pivots en cada columna. • pivotcol2.m: Se comparan los elementos divididos por la norma de su fila.

15. % fichero pivotcol1.m % Gauss y pivotamiento por columnas n=5; A=round(rand(n,n)*20-10); b=round(rand(n,1)*20-10); x=zeros(n,1); % se guardan los datos de partida AA=A; bb=b; % vector con el orden de los pivots p=[1:n]; for k=1:n-1 % buscar el pívot en columna k [pk,ik]=max(abs(A(p(k:n),k))); % intercambiar elementos en p temp=p(k); p(k)=p(ik+k-1); p(ik+k-1)=temp; % hacer ceros en columna pívot for i=p(k+1:n) m=A(i,k)/A(p(k),k); A(i,k:n)=A(i,k:n)- ... m*A(p(k),k:n); b(i)=b(i)-m*b(p(k)); %A(i,k)=m; end end % vuelta atrás x(n)=b(p(n))/A(p(n),n); for i=n-1:-1:1 s=A(p(i),i+1:n)*x(i+1:n); x(i)=(b(p(i))-s)/A(p(i),i); end xe=AA\bb; disp('Solución de Matlab'); disp(xe); disp('Solución calculada'); disp(x); %I=eye(n); %P=I(p,:); %B=P*A; disp('Ya he terminado') Método de Gauss con pivotamiento (5/7) • Implementación del método de Gauss con pivotamiento por columnas

16. Método de Gauss con pivotamiento (6/7) • Método de Gauss con pivotamiento total: • En primer lugar se busca el mayor elemento en valor absoluto de toda la matriz. Este elemento será el primer pívot (fila 4, columna 2). • Se hacen ceros todos los elementos de la columna 2 (menos el pívot). • Se busca el máximo elemento en valor absoluto excluyendo la fila 4 y la columna 2. Este elemento será el 2º pívot (fila 3, columna 4). • Se hacen ceros en la columna 4 excepto en las filas de pivots anteriores. • Se prosigue de la misma forma hasta que todas las filas y/o columnas tienen un pívot. • En la vuelta atrás, las incógnitas se calculan en orden decreciente de columnas en las que han aparecido los pivots, y se calculan en el orden decreciente de filas en las que han aparecido pivots. • Ver programas de Matlab pivototal1.m, pivototal3.m, pivototalMain.m y pivototal.m.

17. Método de Gauss con pivotamiento (7/7) • Método de Gauss con pivotamiento total (cont.): • El pivotamiento total mejora algo la estabilidad del proceso de eliminación de Gauss, pero en general se considera que esta ganancia no merece la pena por el esfuerzo adicional que lleva consigo. • Se utiliza en ocasiones cuando el objetivo es determinar el rango de la matriz o una partición óptima entre variables básicas y variables libres. • Con pivotamiento total puede ser conveniente utilizar escalamiento por filas y por columnas, según la expresión (D1 y D2 perm. diag.): • El escalamiento por columnas implica un cambio de unidades en las incógnitas: en algunos casos puede ser conveniente no mezclar metros con micras y en otros casos puede ser necesario hacerlo. • No hay criterios absolutos para hacer un buen escalamiento por columnas. • Un buen criterio podría ser hacer que el máximo elemento de cada fila y columna esté entre ½ y 1. • Si al hacer el scaling se multiplica y divide siempre por potencias de 2, no se introducen errores de redondeo adicionales.

18. Sistema de m ecuaciones con n incógnitas (1/2) • Resumen del método de eliminación de Gauss • Mediante el método de Gauss con permutación de filas un sistema m×nse puede reducir a la forma de “escalera” siguiente: • Aparecen r pivots (distintos de cero), que no pueden tener elementos distintos de cero debajo y a la izquierda. El rango de A es r. • Conclusiones de la forma de “escalera” PA=LU • La matriz P indica posibles permutaciones de filas para encontrar pivots. • Sólo r de las m ecuaciones eran independientes y sólo r de los m vectores fila (en Rn) eran independientes. Hay m−r filas que se han hecho cero. • Las columnas (y las incógnitas) aparecen divididas en dos grupos: las que tienen pivots (variables básicas o dependientes) y las que no tienen pivots (variables libres o independientes).

19. Sistema de m ecuaciones con n incógnitas (2/2) • Caso r=m, m<n: filas indep.; menos ecuaciones que incógnitas • El vector b siempre pertenece a Im(A): hay solución. • Hay n−mvariables libres: sistema indeterminado con infinitas soluciones. • Caso r<m, mn: filas no independientes; no más filas que columnas • No hay m columnas independientes: puede no haber solución si bIm(A). • Como r<n, hay variables libres: si hay solución, hay infinitas soluciones. • Caso rn, n<m: más ecuaciones que incógnitas • Las columnas pueden ser independientes (r=n) o no (r<n). • Sólo habrá solución si bIm(A). En este caso habrá infinitas soluciones si hay variables libres (r<n). • Caso r=m=n: tantas ecuaciones independientes como incógnitas • En este caso existe solución y ésta es única. • Forma de escalera: • La forma de escalera permite identificar fácilmente estos casos: • Existe solución si las filas que se hacen cero en el primer miembro −en la forma de escalera− se hagan también cero en el segundo miembro.

20. Algunas funciones de Matlab • Operador barra invertida (\) • Por lo general, la operación A\b produce el mismo resultado que inv(A)*b, con algunas diferencias entre las que está la forma de realizar los cálculos. • Si A es una matriz cuadrada no singular, A\b es la solución de Ax=b calculada por eliminación de Gauss. • Si A y b son tales que el sistema tiene infinitas soluciones o ninguna, A\b da la solución de variables libres nulas o la de mínimos cuadrados, respectivamente. • El operador \ es "inteligente", en el sentido de que es capaz de descubrir y aprovechar si A es simétrica y definida-positiva, o si tiene forma triangular. • Función [L,U,P]=lu(A) • Calcula la factorización LU de la matriz PA, donde P es una matriz de permutación que representa el orden de filas en que han aparecido los pivots. • Función [R,c]=rref(A,tol) • Reduce A a forma de escalera con ceros también en la parte superior. El vector c es un vector de enteros cuyo número de elementos indica el rango de la matriz y cuyos elementos indican las columnas de los pivots (variables básicas o dependientes).