1 Tópicos Especiales en Computación NuméricaLa Universidad del Zulia Facultad de Ingeniería División de Estudios para Graduados Asignatura: Tópicos Especiales en Computación Numérica Ecuaciones algebraicas lineales Eliminación de Gauss Descomposición LU e inversión de matrices Matrices especiales y método de Gauss-Siedel Prof. Luis Zerpa, M.Sc.

2 Motivación En el tema anterior (Raíces de ecuaciones) se estudiaron métodos para determinar el valor de x que satisface a una sola ecuación, f(x) = 0 Ahora, nos ocuparemos del caso para determinar los valores de x1, x2,… xn que en forma simultánea satisfacen a un conjunto de ecuaciones algebraicas lineales

3 Motivación Estudiaremos métodos para resolver conjuntos de ecuaciones algebraicas lineales que son de la forma general, a: coeficientes constantes b: constantes n: No. de ecuaciones

4 Métodos empleados antes de la era de las computadorasPara pocas ecuaciones (n ≤ 3), las ecuaciones lineales se pueden resolver con rapidez mediante técnicas simples Método gráfico Para n = 2  la solución corresponde a la intersección de líneas rectas Para n = 3  cada ecuación representa un plano. La solución corresponde al punto donde se intersectan los 3 planos Para n > 3  los métodos gráficos no funcionan Caso que pueden ocasionar problemas Líneas paralelas que jamás se cruzan Dos líneas coinciden Número infinito de soluciones Sistemas singulares Sistemas cercanos a situación singular (mal condicionados) Sensibles a errores de redondeo

5 Métodos empleados antes de la era de las computadorasDeterminantes y Regla de Cramer Determinante de una matriz de coeficientes 22 Determinante de una matriz de coeficientes 33

6 Métodos empleados antes de la era de las computadoraslos sistemas singulares tienen determinante igual a cero Los sistemas mal condicionados (casi singulares) tienen determinante cercano a cero

7 Métodos empleados antes de la era de las computadorasRegla de Cramer Cada incógnita en un sistema de ecuaciones algebraicas lineales puede ser expresada como una fracción de dos determinantes El denominador es el determinante de la matriz de coeficientes El numerador es el determinante de una matriz en la cual se reemplaza la columna de coeficientes de la incógnita por el vector de constantes {B} Para 3 ecuaciones resulta práctica Para n > 3 no es muy eficiente porque la evaluación de determinantes consume tiempo

8 Ecuaciones algebraicas lineales y práctica de la ingenieríaMuchas de las ecuaciones fundamentales de ingeniería están basadas en leyes de conservación En términos matemáticos, esos principios conducen a ecuaciones de balance que relacionan el comportamiento del sistema con las propiedades o características y los estímulos externos que actúan sobre el sistema En el capítulo anterior se trabajó con sistemas de un solo componente que resultaba en una sola ecuación que se resolvía con métodos de cálculo de raíces Los sistemas multicomponentes resultan en un conjunto agrupado de ecuaciones matemáticas que deben ser resueltas simultáneamente

9 Ecuaciones algebraicas lineales y práctica de la ingenieríaProblemas multicomponentes surgen tanto de modelos matemáticos de variables agrupadas como de variables distribuidas Los problemas de variables agrupadas involucran componentes finitos acoplados Armaduras Reactores Circuitos eléctricos Los problemas de variables distribuidas intentan describir detalles espaciales de los sistemas sobre una base continua Las ecuaciones diferenciales derivadas a partir de leyes de conservación especifican la distribución de la variable dependiente para tales sistemas Esas ecuaciones se pueden resolver numéricamente al convertirlas en un sistema equivalente de ecuaciones algebraicas simultáneas

10 Ecuaciones algebraicas lineales y práctica de la ingenieríaAdemás de los problemas físicos, las ecuaciones algebraicas lineales simultáneas surgen también en diferentes contextos de problemas matemáticos Algunas técnicas numéricas de uso general que emplean ecuaciones simultáneas son el análisis de regresión y la interpolación segmentaria

11 Antecedentes matemáticos necesariosNotación matricial En la solución de ecuaciones algebraicas lineales No. de ecuaciones = No. de filas No. de variables = No. de columnas Filas Columnas Debe ser igual para que una solución única sea posible

12 Antecedentes matemáticos necesariosOperaciones de matrices Suma Multiplicación de una matriz por un escalar: se multiplica cada elemento de la matriz por el escalar Multiplicación de dos matrices Es imposible multiplicar dos matrices si el número de columnas de la primera no es igual al número de filas de la segunda

13 Antecedentes matemáticos necesariosOperaciones de matrices División: la división de una matriz no es una operación definida sin embargo, si una matriz A es cuadrada y no singular, existe otra matriz A-1, llamada inversa de A para la cual AA-1 = A-1A=I Transpuesta de una matriz: consiste en transformar sus filas en columnas y viceversa Traza de una matriz: es la suma de los elementos de su diagonal principal

14 Representación de ecuaciones algebraicas lineales en forma matricialUna manera formal para obtener la solución usando algebra matricial es multiplicando cada lado de la ecuación por la inversa de A No resulta muy eficiente por la obtención de la inversa de la matriz Es necesario métodos numéricos

15 Métodos numéricos para la solución de ecuaciones algebraicas linealesEliminación de Gauss Descomposición LU  valiosa para casos donde se necesita evaluar muchos vectores del lado derecho. Permite hacer eficiente el cálculo de la matriz inversa Técnicas eficientes para la solución de sistemas tridiagonales (matrices en banda) Método de Gauss-Seidel  método iterativo

16 Eliminación de Gauss Este método involucra una combinación de ecuaciones para eliminar las incógnitas Es uno de los métodos más antiguos y sigue siendo uno de los algoritmos de mayor importancia

17 Eliminación de Gauss Eliminación de incógnitasLa estrategia básica es multiplicar las ecuaciones por constantes, de tal forma que se elimine una de las incógnitas cuando se combinen las ecuaciones El resultado es una sola ecuación que se puede resolver para la incógnita restante Este valor se sustituye en las ecuaciones originales para calcular la otra variable Este método representa la base para la eliminación de Gauss Se puede extender a grandes sistemas de ecuaciones desarrollando un esquema sistemático para eliminar incógnitas y sustituir hacia atrás

18 Eliminación de Gauss Eliminación de Gauss simple  el método consiste en dos fases Eliminación de incógnitas Reduce el conjunto de ecuaciones a un sistema triangular superior Primero se elimina la primera incógnita, x1, desde la segunda hasta la n-enésima fila, multiplicando por a21/a11 a la primera ecuación, luego restando ésta a la segunda El procedimiento es repetido para las ecuaciones restantes Para este paso la ecuación 1 es la ecuación pivote y a11 es el coeficiente pivote Solución por sustitución hacia atrás Al finalizar la eliminación, la ecuación n puede resolverse para xn Para evaluar las x restantes

19 Eliminación de Gauss Seudo código DO k = 1, n-1 DO i = k+1, nFactor = ai,k/ak,k DO j = k+1, n ai,j = ai,j - factor * ak,j END DO bi = bi -factor * bk xn =bn / an,n DO i = n-1, 1, -1 sum = 0 DO j = i+1, n sum = sum + ai,j * xj xi = (bi - sum) / ai,i

20 Eliminación de Gauss Ejemplo Sustitución Eliminación X ={3.0000 ]; B = [7.85 -19.3 71.4]; A =[ ] B ={7.8500 } A =[ ] B ={7.8500 } A =[ ] B ={7.8500 } X ={3.0000 7.0000} Sustitución big) big = dummy. p = ii. END IF. END DO. IF (p ≠ k) DO jj = k, n. dummy = ap,jj. ap,jj = ak,jj. ak,jj = dummy. dummy = bp. bp = bk. bk = dummy.", "width": "800" } 32 Técnicas para mejorar las soluciones del método de eliminación de GaussEscalamiento El escalamiento revela si el pivoteo es necesario Si es necesario, se pivotea pero se retienen los coeficientes originales de la ecuación escalada El escalamiento se usa para calcular los valores escalados de los coeficientes que servirán como un criterio de pivoteo Ejemplo Escalamiento como criterio de pivoteo Manteniendo coeficientes originales Problema original Escalando y pivoteando { "@context": "http://schema.org", "@type": "ImageObject", "contentUrl": "http://slideplayer.es/4267461/14/images/32/T%C3%A9cnicas+para+mejorar+las+soluciones+del+m%C3%A9todo+de+eliminaci%C3%B3n+de+Gauss.jpg", "name": "Técnicas para mejorar las soluciones del método de eliminación de Gauss", "description": "Escalamiento. El escalamiento revela si el pivoteo es necesario. Si es necesario, se pivotea pero se retienen los coeficientes originales de la ecuación escalada. El escalamiento se usa para calcular los valores escalados de los coeficientes que servirán como un criterio de pivoteo. Ejemplo. Escalamiento como criterio de pivoteo. Manteniendo coeficientes originales. Problema original. Escalando y pivoteando.", "width": "800" } 33 Descomposición LU e Inversión de matricesEl principal atractivo de este método es que el paso de eliminación, que consume tiempo, se puede reformular de tal manera que involucre sólo operaciones sobre los elementos de la matriz de coeficientes, A De esta forma, es muy adecuado para aquellas situaciones donde se debe evaluar muchos vectores {B} El método de eliminación de Gauss puede implementarse como una descomposición LU La descomposición LU proporciona un medio eficaz para calcular la matriz inversa, la cual a su vez permite evaluar la condición de un sistema { "@context": "http://schema.org", "@type": "ImageObject", "contentUrl": "http://slideplayer.es/4267461/14/images/33/Descomposici%C3%B3n+LU+e+Inversi%C3%B3n+de+matrices.jpg", "name": "Descomposición LU e Inversión de matrices", "description": "El principal atractivo de este método es que el paso de eliminación, que consume tiempo, se puede reformular de tal manera que involucre sólo operaciones sobre los elementos de la matriz de coeficientes, A. De esta forma, es muy adecuado para aquellas situaciones donde se debe evaluar muchos vectores {B} El método de eliminación de Gauss puede implementarse como una descomposición LU. La descomposición LU proporciona un medio eficaz para calcular la matriz inversa, la cual a su vez permite evaluar la condición de un sistema.", "width": "800" } 34 Descomposición LU Partiendo de un sistema de ecuaciones lineales de la forma, Este se puede ordenar como, El primer paso de la eliminación de Gauss resulta en un sistema con una matriz triángular superior Que puede ser expresada como, Ahora, suponga que existe una matriz triangular inferior con números 1 sobre la diagonal que tiene la siguiente propiedad si esta propiedad se cumple, de las reglas de multiplicación de matrices se obtiene, { "@context": "http://schema.org", "@type": "ImageObject", "contentUrl": "http://slideplayer.es/4267461/14/images/34/Descomposici%C3%B3n+LU+Partiendo+de+un+sistema+de+ecuaciones+lineales+de+la+forma%2C+Este+se+puede+ordenar+como%2C.jpg", "name": "Descomposición LU Partiendo de un sistema de ecuaciones lineales de la forma, Este se puede ordenar como,", "description": "El primer paso de la eliminación de Gauss resulta en un sistema con una matriz triángular superior. Que puede ser expresada como, Ahora, suponga que existe una matriz triangular inferior con números 1 sobre la diagonal. que tiene la siguiente propiedad. si esta propiedad se cumple, de las reglas de multiplicación de matrices se obtiene,", "width": "800" } 35 Descomposición LU Estrategia para resolver el sistemaPaso de descomposición LU: la matriz [A], se factoriza o descompone en matrices triangulares inferior [L] y superior [U] Paso de sustitución: [L] y [U] se usan par determinar una solución {X} para un vector {B}. Este paso consta de dos subpasos: Se determina el vector intermedio {D} resolviendo [L]{D}={B} por sustitución hacia delante, debido a que [L] es una matriz triangular inferior Se determina {X} resolviendo [U]{X}={D} por sustitución hacia atrás { "@context": "http://schema.org", "@type": "ImageObject", "contentUrl": "http://slideplayer.es/4267461/14/images/35/Descomposici%C3%B3n+LU+Estrategia+para+resolver+el+sistema.jpg", "name": "Descomposición LU Estrategia para resolver el sistema", "description": "Paso de descomposición LU: la matriz [A], se factoriza o descompone en matrices triangulares inferior [L] y superior [U] Paso de sustitución: [L] y [U] se usan par determinar una solución {X} para un vector {B}. Este paso consta de dos subpasos: Se determina el vector intermedio {D} resolviendo [L]{D}={B} por sustitución hacia delante, debido a que [L] es una matriz triangular inferior. Se determina {X} resolviendo [U]{X}={D} por sustitución hacia atrás.", "width": "800" } 36 Descomposición LU Descomposición LU con base en la eliminación de Gauss Partiendo de una matriz de coeficientes, se llega a una matriz triangular superior Para llegar a esta matriz [U] La matriz triangular inferior que tiene la propiedad requerida para la descomposición LU es se multiplicó la fila 1 por el factor f21 = a21/a11 y restando el resultado a la fila 2 se eliminó a21 se multiplicó la fila 1 por el factor f31 = a31/a11 y restando el resultado a la fila 3 se eliminó a31 se multiplicó la fila 2 por el factor f32 = a32’/a22’ y restando el resultado a la fila 3 se eliminó a32’ { "@context": "http://schema.org", "@type": "ImageObject", "contentUrl": "http://slideplayer.es/4267461/14/images/36/Descomposici%C3%B3n+LU+Descomposici%C3%B3n+LU+con+base+en+la+eliminaci%C3%B3n+de+Gauss..jpg", "name": "Descomposición LU Descomposición LU con base en la eliminación de Gauss.", "description": "Partiendo de una matriz de coeficientes, se llega a una matriz triangular superior. Para llegar a esta matriz [U] La matriz triangular inferior que tiene la propiedad requerida para la descomposición LU es. se multiplicó la fila 1 por el factor f21 = a21/a11 y restando el resultado a la fila 2 se eliminó a21. se multiplicó la fila 1 por el factor f31 = a31/a11 y restando el resultado a la fila 3 se eliminó a31. se multiplicó la fila 2 por el factor f32 = a32’/a22’ y restando el resultado a la fila 3 se eliminó a32’", "width": "800" } 37 Descomposición LU Descomposición LU con base en la eliminación de Gauss Paso de descomposición LU: la matriz [A], se factoriza o descompone en matrices triangulares inferior [L] y superior [U] Paso de sustitución: [L] y [U] se usan par determinar una solución {X} para un vector {B}. Este paso consta de dos subpasos: Se determina el vector intermedio {D} resolviendo [L]{D}={B} por sustitución hacia delante Se determina {X} resolviendo [U]{X}={D} por sustitución hacia atrás { "@context": "http://schema.org", "@type": "ImageObject", "contentUrl": "http://slideplayer.es/4267461/14/images/37/Descomposici%C3%B3n+LU+Descomposici%C3%B3n+LU+con+base+en+la+eliminaci%C3%B3n+de+Gauss..jpg", "name": "Descomposición LU Descomposición LU con base en la eliminación de Gauss.", "description": "Paso de descomposición LU: la matriz [A], se factoriza o descompone en matrices triangulares inferior [L] y superior [U] Paso de sustitución: [L] y [U] se usan par determinar una solución {X} para un vector {B}. Este paso consta de dos subpasos: Se determina el vector intermedio {D} resolviendo [L]{D}={B} por sustitución hacia delante. Se determina {X} resolviendo [U]{X}={D} por sustitución hacia atrás.", "width": "800" } 38 Seudo código para descomposición LUSUB Ludecomp(a, b, n, tol, x, er) DIM on, sn er = 0 CALL Decompose(a, n, tol, o, s, er) IF er <> -1 THEN CALL Substitute(a, o, n, b, x) END IF END Ludecomp { "@context": "http://schema.org", "@type": "ImageObject", "contentUrl": "http://slideplayer.es/4267461/14/images/38/Seudo+c%C3%B3digo+para+descomposici%C3%B3n+LU.jpg", "name": "Seudo código para descomposición LU", "description": "SUB Ludecomp(a, b, n, tol, x, er) DIM on, sn. er = 0. CALL Decompose(a, n, tol, o, s, er) IF er -1 THEN. CALL Substitute(a, o, n, b, x) END IF. END Ludecomp.", "width": "800" } 39 Seudo código para descomposición LUSUB Decompose(a, n, tol, o, s, er) DO i = 1,n oi = i si = ABS(ai,1) DO j = 2,n IF ABS(ai,j)> si THEN si = ABS(ai,j) END DO DO k = 1,n-1 CALL Pivot(a, o, s, n, k) IF ABS(ao(k),k/so(k)) < tol THEN er = -1 PRINT ao(k),k/so(k) EXIT DO END IF DO i = k+1,n factor = ao(i),k/ao(k),k ao(i),k = factor DO j = k+1,n ao(i),j = ao(i),j - factor*a IF ABS(ao(k),k/so(k))< tol THEN END Decompose { "@context": "http://schema.org", "@type": "ImageObject", "contentUrl": "http://slideplayer.es/4267461/14/images/39/Seudo+c%C3%B3digo+para+descomposici%C3%B3n+LU.jpg", "name": "Seudo código para descomposición LU", "description": "SUB Decompose(a, n, tol, o, s, er) DO i = 1,n. oi = i. si = ABS(ai,1) DO j = 2,n. IF ABS(ai,j)> si THEN si = ABS(ai,j) END DO. DO k = 1,n-1. CALL Pivot(a, o, s, n, k) IF ABS(ao(k),k/so(k)) < tol THEN. er = -1. PRINT ao(k),k/so(k) EXIT DO. END IF. DO i = k+1,n. factor = ao(i),k/ao(k),k. ao(i),k = factor. DO j = k+1,n. ao(i),j = ao(i),j - factor*a. IF ABS(ao(k),k/so(k))< tol THEN. END Decompose.", "width": "800" } 40 Seudo código para descomposición LUSUB Pivot(a, o, s, n, k) p = k big = ABS(ao(k),k/so(k)) DO ii = k+1,n dummy = ABS(ao(ii),k/so(ii)) IF dummy > big THEN big = dummy p = ii END IF END DO dummy = op op = ok ok = dummy END Pivot SUB Substitute(a, o, n, b, x) DO i = 2,n sum = bo(i) DO 1,i-1 sum = sum - ao(i),j*bo(j) END DO bo(i) = sum xn = bo(n)/ao(n),n DO i = n-1,1,-1 sum = 0 DO j = i+1,n sum = sum + ao(i),j*xj xi =(bo(i) - sum)/ao(i),i END Substitute { "@context": "http://schema.org", "@type": "ImageObject", "contentUrl": "http://slideplayer.es/4267461/14/images/40/Seudo+c%C3%B3digo+para+descomposici%C3%B3n+LU.jpg", "name": "Seudo código para descomposición LU", "description": "SUB Pivot(a, o, s, n, k) p = k. big = ABS(ao(k),k/so(k)) DO ii = k+1,n. dummy = ABS(ao(ii),k/so(ii)) IF dummy > big THEN. big = dummy. p = ii. END IF. END DO. dummy = op. op = ok. ok = dummy. END Pivot. SUB Substitute(a, o, n, b, x) DO i = 2,n. sum = bo(i) DO 1,i-1. sum = sum - ao(i),j*bo(j) END DO. bo(i) = sum. xn = bo(n)/ao(n),n. DO i = n-1,1,-1. sum = 0. DO j = i+1,n. sum = sum + ao(i),j*xj. xi =(bo(i) - sum)/ao(i),i. END Substitute.", "width": "800" } 41 Matriz inversa Para una matriz cuadrada [A], hay otra matriz [A]-1 conocida como la inversa de [A], para la cual se cumple, [A] [A]-1 = [A]-1[A] = I La matriz inversa se puede calcular en una forma de columna por columna a partir de vectores unitarios como vector de constantes del sistema de ecuaciones lineales algebraicas Por ejemplo, para determinar la primera columna de la matriz inversa se resuelve el sistema con el vector de constantes B=[1 0 0]T para determinar la segunda columna se usa B=[0 1 0]T y así sucesivamente { "@context": "http://schema.org", "@type": "ImageObject", "contentUrl": "http://slideplayer.es/4267461/14/images/41/Matriz+inversa+Para+una+matriz+cuadrada+%5BA%5D%2C+hay+otra+matriz+%5BA%5D-1+conocida+como+la+inversa+de+%5BA%5D%2C+para+la+cual+se+cumple%2C.jpg", "name": "Matriz inversa Para una matriz cuadrada [A], hay otra matriz [A]-1 conocida como la inversa de [A], para la cual se cumple,", "description": "[A] [A]-1 = [A]-1[A] = I. La matriz inversa se puede calcular en una forma de columna por columna a partir de vectores unitarios como vector de constantes del sistema de ecuaciones lineales algebraicas. Por ejemplo, para determinar la primera columna de la matriz inversa se resuelve el sistema con el vector de constantes B=[1 0 0]T. para determinar la segunda columna se usa B=[0 1 0]T. y así sucesivamente.", "width": "800" } 42 Matriz inversa La descomposición LU representa la mejor forma para implementar el cálculo de la matriz inversa, ya que una vez obtenida la descomposición LU de la matriz A se puede calcular su inversa resolviendo cada columna con los vectores unitarios como constantes Ejemplo: Determinar la inversa de A = [ ] 1. Descomposición LU de la matriz A LU(A) = [ ] B1 = [1] [0] [0] B2 = [0] [1] [0] B3 = [0] [0] [1] 2. Sustitución, InvA(1) = [0.3325] [ ] [ ] InvA(2) = [0.0049] [0.1429] [0.0027] InvA(3) = [0.0068] [0.0042] [0.0999] { "@context": "http://schema.org", "@type": "ImageObject", "contentUrl": "http://slideplayer.es/4267461/14/images/42/Matriz+inversa.jpg", "name": "Matriz inversa", "description": "La descomposición LU representa la mejor forma para implementar el cálculo de la matriz inversa, ya que una vez obtenida la descomposición LU de la matriz A se puede calcular su inversa resolviendo cada columna con los vectores unitarios como constantes. Ejemplo: Determinar la inversa de. A = [ 3 -0.1 -0.2 0.1 7 -0.3 0.3 -0.2 10] 1. Descomposición LU de la matriz A. LU(A) = [ 3 -0.1 -0.2 0.0333 7.0033 -0.2933 0.1 -0.0271 10.0120] B1 = [1] [0] [0] B2 = [0] [1] [0] B3 = [0] [0] [1] 2. Sustitución, InvA(1) = [0.3325] [-0.0052] [-0.0101] InvA(2) = [0.0049] [0.1429] [0.0027] InvA(3) = [0.0068] [0.0042] [0.0999]", "width": "800" } 43 Análisis de error y condición del sistemaLa matriz inversa permite determinar si un sistema está mal condicionado, para esto existen 3 métodos: Escalar la matriz de coeficientes [A], de tal manera que el elemento más grande en cada fila sea 1. Si al invertir la matriz escalada existen elementos de la inversa [A]-1 que sean varios ordenes de magnitud mayores que la unidad, es probable que el sistema esté mal condicionado Multiplicar la inversa por la matriz de coeficientes original y verificar que [A][A]-1  I. Si no es así, indica que el sistema está mal condicionado Invertir la matriz inversa y verificar que el resultado está lo suficientemente cercano a la matriz original. Si no es así, indica que el sistema está mal condicionado { "@context": "http://schema.org", "@type": "ImageObject", "contentUrl": "http://slideplayer.es/4267461/14/images/43/An%C3%A1lisis+de+error+y+condici%C3%B3n+del+sistema.jpg", "name": "Análisis de error y condición del sistema", "description": "La matriz inversa permite determinar si un sistema está mal condicionado, para esto existen 3 métodos: Escalar la matriz de coeficientes [A], de tal manera que el elemento más grande en cada fila sea 1. Si al invertir la matriz escalada existen elementos de la inversa [A]-1 que sean varios ordenes de magnitud mayores que la unidad, es probable que el sistema esté mal condicionado. Multiplicar la inversa por la matriz de coeficientes original y verificar que [A][A]-1  I. Si no es así, indica que el sistema está mal condicionado. Invertir la matriz inversa y verificar que el resultado está lo suficientemente cercano a la matriz original. Si no es así, indica que el sistema está mal condicionado.", "width": "800" } 44 Número de condición de una matrizEste número mide la sensibilidad de la solución de un sistema de ecuaciones lineales a errores en los datos Valores cercanos a 1 indican que el sistema está bien condicionado Valores grandes indican que la matriz es casi singular { "@context": "http://schema.org", "@type": "ImageObject", "contentUrl": "http://slideplayer.es/4267461/14/images/44/N%C3%BAmero+de+condici%C3%B3n+de+una+matriz.jpg", "name": "Número de condición de una matriz", "description": "Este número mide la sensibilidad de la solución de un sistema de ecuaciones lineales a errores en los datos. Valores cercanos a 1 indican que el sistema está bien condicionado. Valores grandes indican que la matriz es casi singular.", "width": "800" } 45 Número de condición de una matrizNorma Euclidiana - Para vector - Para matriz Norma columna-suma Norma 2, o normal espectral max es el eigenvalor más grande de [A]T[A]. Esta es la la norma mínima, por lo tanto proporciona la medida de tamaño más ajustada Norma infinita / matriz uniforme / fila-suma Norma de Frobenius { "@context": "http://schema.org", "@type": "ImageObject", "contentUrl": "http://slideplayer.es/4267461/14/images/45/N%C3%BAmero+de+condici%C3%B3n+de+una+matriz.jpg", "name": "Número de condición de una matriz", "description": "Norma Euclidiana. - Para vector. - Para matriz. Norma columna-suma. Norma 2, o normal espectral. max es el eigenvalor más grande de [A]T[A]. Esta es la la norma mínima, por lo tanto proporciona la medida de tamaño más ajustada. Norma infinita / matriz uniforme / fila-suma. Norma de Frobenius.", "width": "800" } 46 Matrices especiales Matrices banda Matrices simétricasUna matriz banda es una matriz cuadrada en la que todos sus elementos son cero, con excepción de una banda centrada sobre la diagonal principal Las dimensiones de un sistema de banda se pueden cuantificar con dos parámetros: Ancho de banda, BW Ancho de media banda, HBW Estos se relacionan por BW = 2HBW + 1 En general un sistema de banda es aquel para el cual aij = 0 si |i – j| > HBW HBW+1 Diagonal HBW BW { "@context": "http://schema.org", "@type": "ImageObject", "contentUrl": "http://slideplayer.es/4267461/14/images/46/Matrices+especiales+Matrices+banda+Matrices+sim%C3%A9tricas.jpg", "name": "Matrices especiales Matrices banda Matrices simétricas", "description": "Una matriz banda es una matriz cuadrada en la que todos sus elementos son cero, con excepción de una banda centrada sobre la diagonal principal. Las dimensiones de un sistema de banda se pueden cuantificar con dos parámetros: Ancho de banda, BW. Ancho de media banda, HBW. Estos se relacionan por BW = 2HBW + 1. En general un sistema de banda es aquel para el cual aij = 0 si |i – j| > HBW. HBW+1. Diagonal. HBW. BW.", "width": "800" } 47 Matrices especiales La eliminación de Gauss o la descomposición LU pueden emplearse para resolver sistemas de banda, pero si el pivoteo no es necesario resultan ineficientes, porque se utilizaría tiempo y espacio innecesario en el almacenamiento y manejo de ceros Si se sabe de antemano que el pivoteo es innecesario, se pueden desarrollar algoritmos muy eficientes que no involucren los elementos cero fuera de la banda { "@context": "http://schema.org", "@type": "ImageObject", "contentUrl": "http://slideplayer.es/4267461/14/images/47/Matrices+especiales.jpg", "name": "Matrices especiales", "description": "La eliminación de Gauss o la descomposición LU pueden emplearse para resolver sistemas de banda, pero si el pivoteo no es necesario resultan ineficientes, porque se utilizaría tiempo y espacio innecesario en el almacenamiento y manejo de ceros. Si se sabe de antemano que el pivoteo es innecesario, se pueden desarrollar algoritmos muy eficientes que no involucren los elementos cero fuera de la banda.", "width": "800" } 48 Sistemas tridiagonalesUn sistema tridiagonal (ancho de banda = 3) se puede expresar como Se cambia la notación a  e, f, g b  r para evitar guardar ceros, ahorrando espacio { "@context": "http://schema.org", "@type": "ImageObject", "contentUrl": "http://slideplayer.es/4267461/14/images/48/Sistemas+tridiagonales.jpg", "name": "Sistemas tridiagonales", "description": "Un sistema tridiagonal (ancho de banda = 3) se puede expresar como. Se cambia la notación. a  e, f, g. b  r. para evitar guardar ceros, ahorrando espacio.", "width": "800" } 49 Sistemas tridiagonalesEjemplo B = [ ] A = [ ] f = [ ] g = [ ] e = [ ] { "@context": "http://schema.org", "@type": "ImageObject", "contentUrl": "http://slideplayer.es/4267461/14/images/49/Sistemas+tridiagonales.jpg", "name": "Sistemas tridiagonales", "description": "Ejemplo. B = [14.5000 -4 -6 -8 -10 -12 -14 -16 -30] A = [-5 1.5 0 0 0 0 0 0 0 2.5 -5 1.5 0 0 0 0 0 0 0 2.5 -5 1.5 0 0 0 0 0 0 0 2.5 -5 1.5 0 0 0 0 0 0 0 2.5 -5 1.5 0 0 0 0 0 0 0 2.5 -5 1.5 0 0 0 0 0 0 0 2.5 -5 1.5 0 0 0 0 0 0 0 2.5 -5 1.5 0 0 0 0 0 0 0 2.5 -5] f = [ -5 -5 -5 -5 -5 -5 -5 -5 -5] g = [ 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 0] e = [ 0 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5]", "width": "800" } 50 Algoritmo de Thomas (TDMA, tridiagonal matrix algorithm)El algoritmo consiste en tres pasos: Descomposición Sustitución hacia delante Sustitución hacia atrás Manteniendo todas las ventajas de la descomposición LU { "@context": "http://schema.org", "@type": "ImageObject", "contentUrl": "http://slideplayer.es/4267461/14/images/50/Algoritmo+de+Thomas+%28TDMA%2C+tridiagonal+matrix+algorithm%29.jpg", "name": "Algoritmo de Thomas (TDMA, tridiagonal matrix algorithm)", "description": "El algoritmo consiste en tres pasos: Descomposición. Sustitución hacia delante. Sustitución hacia atrás. Manteniendo todas las ventajas de la descomposición LU.", "width": "800" } 51 Seudo código para algoritmo de Thomasa) Descomposición DO k = 2, n ek = ek/fk-1 fk = fk - ek*gk-1 END DO b) Sustitución hacia adelante DO k = 2,n rk = rk - ek*rk-1 c) Sustitución hacia atrás xn = rn/fn DO k = n-1,1,-1 xk = (rk - gk*xk+1)/fk { "@context": "http://schema.org", "@type": "ImageObject", "contentUrl": "http://slideplayer.es/4267461/14/images/51/Seudo+c%C3%B3digo+para+algoritmo+de+Thomas.jpg", "name": "Seudo código para algoritmo de Thomas", "description": "a) Descomposición. DO k = 2, n. ek = ek/fk-1. fk = fk - ek*gk-1. END DO. b) Sustitución hacia adelante. DO k = 2,n. rk = rk - ek*rk-1. c) Sustitución hacia atrás. xn = rn/fn. DO k = n-1,1,-1. xk = (rk - gk*xk+1)/fk.", "width": "800" } 52 Descomposición de CholeskyLa descomposición o factorización de Cholesky expresa una matriz simétrica como el producto de una matriz triangular y su transpuesta A = L∙LT  L: matriz triangular inferior No todas las matrices simétricas se pueden factorizar de esta forma Las matrices que tienen este tipo de factorización son las matrices simétricas definidas positivas. Esto implica que todos los elementos de la diagonal sean positivos y que los elementos fuera de la diagonal no sean muy grandes { "@context": "http://schema.org", "@type": "ImageObject", "contentUrl": "http://slideplayer.es/4267461/14/images/52/Descomposici%C3%B3n+de+Cholesky.jpg", "name": "Descomposición de Cholesky", "description": "La descomposición o factorización de Cholesky expresa una matriz simétrica como el producto de una matriz triangular y su transpuesta. A = L∙LT  L: matriz triangular inferior. No todas las matrices simétricas se pueden factorizar de esta forma. Las matrices que tienen este tipo de factorización son las matrices simétricas definidas positivas. Esto implica que todos los elementos de la diagonal sean positivos y que los elementos fuera de la diagonal no sean muy grandes.", "width": "800" } 53 Descomposición de CholeskyLos términos de la descomposición se pueden multiplicar entre si. El resultado se puede expresar en forma simple por relaciones recurrentes Para la fila k { "@context": "http://schema.org", "@type": "ImageObject", "contentUrl": "http://slideplayer.es/4267461/14/images/53/Descomposici%C3%B3n+de+Cholesky.jpg", "name": "Descomposición de Cholesky", "description": "Los términos de la descomposición se pueden multiplicar entre si. El resultado se puede expresar en forma simple por relaciones recurrentes. Para la fila k.", "width": "800" } 54 Descomposición de CholeskyEjemplo, Para k = 1 Para k = 2 Para k = 3 Sustitución hacia adelante Sustitución hacia atrás { "@context": "http://schema.org", "@type": "ImageObject", "contentUrl": "http://slideplayer.es/4267461/14/images/54/Descomposici%C3%B3n+de+Cholesky.jpg", "name": "Descomposición de Cholesky", "description": "Ejemplo, Para k = 1. Para k = 2. Para k = 3. Sustitución hacia adelante. Sustitución hacia atrás.", "width": "800" } 55 Seudo código para la descomposición de Choleskyfor k = 1:n for i = 1:k-1 sum = 0; for j = 1:i-1 sum = sum + A(i,j)*A(k,j); end A(k,i) = (A(k,i) - sum)/A(i,i); for j = 1:k-1 sum = sum + A(k,j)^2; A(k,k) = sqrt(A(k,k) - sum); { "@context": "http://schema.org", "@type": "ImageObject", "contentUrl": "http://slideplayer.es/4267461/14/images/55/Seudo+c%C3%B3digo+para+la+descomposici%C3%B3n+de+Cholesky.jpg", "name": "Seudo código para la descomposición de Cholesky", "description": "for k = 1:n. for i = 1:k-1. sum = 0; for j = 1:i-1. sum = sum + A(i,j)*A(k,j); end. A(k,i) = (A(k,i) - sum)/A(i,i); for j = 1:k-1. sum = sum + A(k,j)^2; A(k,k) = sqrt(A(k,k) - sum);", "width": "800" } 56 Método de Gauss-SeidelEste es un método iterativo Dado un conjunto de ecuaciones, AX = B Si los elementos de la diagonal son diferentes de cero, se puede resolver la ecuación i para la variable i, donde i = 1…n Se puede empezar el proceso de solución al escoger los valores iniciales de las variables x (xi = 0) { "@context": "http://schema.org", "@type": "ImageObject", "contentUrl": "http://slideplayer.es/4267461/14/images/56/M%C3%A9todo+de+Gauss-Seidel.jpg", "name": "Método de Gauss-Seidel", "description": "Este es un método iterativo. Dado un conjunto de ecuaciones, AX = B. Si los elementos de la diagonal son diferentes de cero, se puede resolver la ecuación i para la variable i, donde i = 1…n. Se puede empezar el proceso de solución al escoger los valores iniciales de las variables x (xi = 0)", "width": "800" } 57 Método de Gauss-SeidelLos valores iniciales se sustituyen en la primera ecuación para calcular un nuevo valor para x1 Este nuevo valor de x1 junto con los demás valores iniciales se sustituyen en la segunda ecuación para calcular un nuevo valor para x2 Este proceso se repite hasta calcular los nuevos valores de las n variables Después se regresa a la primera ecuación y se repite todo el procedimiento hasta que la solución converja a la solución real La convergencia se puede verificar usando el criterio, { "@context": "http://schema.org", "@type": "ImageObject", "contentUrl": "http://slideplayer.es/4267461/14/images/57/M%C3%A9todo+de+Gauss-Seidel.jpg", "name": "Método de Gauss-Seidel", "description": "Los valores iniciales se sustituyen en la primera ecuación para calcular un nuevo valor para x1. Este nuevo valor de x1 junto con los demás valores iniciales se sustituyen en la segunda ecuación para calcular un nuevo valor para x2. Este proceso se repite hasta calcular los nuevos valores de las n variables. Después se regresa a la primera ecuación y se repite todo el procedimiento hasta que la solución converja a la solución real. La convergencia se puede verificar usando el criterio,", "width": "800" } 58 Criterio de convergencia del método de Gauss-SeidelEste método es similar en esencia al método de iteración de punto de fijo que se usa para el cálculo de raíces de una ecuación Presenta las mismas desventajas: En algunos casos no converge En algunos casos la convergencia es lenta Las condiciones suficientes para la convergencia de dos ecuaciones no lineales también aplican para ecuaciones lineales cuando se usa Gauss-Seidel { "@context": "http://schema.org", "@type": "ImageObject", "contentUrl": "http://slideplayer.es/4267461/14/images/58/Criterio+de+convergencia+del+m%C3%A9todo+de+Gauss-Seidel.jpg", "name": "Criterio de convergencia del método de Gauss-Seidel", "description": "Este método es similar en esencia al método de iteración de punto de fijo que se usa para el cálculo de raíces de una ecuación. Presenta las mismas desventajas: En algunos casos no converge. En algunos casos la convergencia es lenta. Las condiciones suficientes para la convergencia de dos ecuaciones no lineales también aplican para ecuaciones lineales cuando se usa Gauss-Seidel.", "width": "800" } 59 Criterio de convergencia del método de Gauss-SeidelEn el caso de dos ecuaciones el método de Gauss-Seidel se expresa como Las derivadas parciales de estas ecuaciones con respecto a las variables son Para que se cumplan las condiciones suficientes de convergencia { "@context": "http://schema.org", "@type": "ImageObject", "contentUrl": "http://slideplayer.es/4267461/14/images/59/Criterio+de+convergencia+del+m%C3%A9todo+de+Gauss-Seidel.jpg", "name": "Criterio de convergencia del método de Gauss-Seidel", "description": "En el caso de dos ecuaciones el método de Gauss-Seidel se expresa como. Las derivadas parciales de estas ecuaciones con respecto a las variables son. Para que se cumplan las condiciones suficientes de convergencia.", "width": "800" } 60 Criterio de convergencia del método de Gauss-SeidelEl valor absoluto de la pendiente de las ecuaciones rectas debe ser menor que la unidad para asegurar convergencia Rerformulando, El elemento diagonal debe ser mayor que el elemento fuera de la diagonal para cada fila (sistemas diagonal dominantes) Generalizando para n ecuaciones El criterio es suficiente pero no necesario para convergencia { "@context": "http://schema.org", "@type": "ImageObject", "contentUrl": "http://slideplayer.es/4267461/14/images/60/Criterio+de+convergencia+del+m%C3%A9todo+de+Gauss-Seidel.jpg", "name": "Criterio de convergencia del método de Gauss-Seidel", "description": "El valor absoluto de la pendiente de las ecuaciones rectas debe ser menor que la unidad para asegurar convergencia. Rerformulando, El elemento diagonal debe ser mayor que el elemento fuera de la diagonal para cada fila (sistemas diagonal dominantes) Generalizando para n ecuaciones. El criterio es suficiente pero no necesario para convergencia.", "width": "800" } 61 Mejoras a la convergencia por medio de relajaciónLa relajación representa una ligera modificación al método de Gauss-Seidel y está diseñada para mejorar la convergencia Después de calcular cada nuevo valor de x, ese valor se modifica por un promedio ponderado de los resultados de las iteraciones anterior y actual: : es el coeficiente de relajación que tiene un valor entre 0 y 2 Si  = 1  el resultado no se modifica Si 0 <  < 1  el resultado es un promedio ponderado de xinuevo y xianterior (subrelajación), se usa para hacer que un sistema no convergente, converja o converja más rápido al amortiguar sus oscilaciones Si 1 <  < 2  se le da una ponderación extra al valor actual (sobrerelajación), acelera la convergencia de un sistema que ya es convergente. También es conocida como sobrerelajación simultánea o sucesiva, SOR { "@context": "http://schema.org", "@type": "ImageObject", "contentUrl": "http://slideplayer.es/4267461/14/images/61/Mejoras+a+la+convergencia+por+medio+de+relajaci%C3%B3n.jpg", "name": "Mejoras a la convergencia por medio de relajación", "description": "La relajación representa una ligera modificación al método de Gauss-Seidel y está diseñada para mejorar la convergencia. Después de calcular cada nuevo valor de x, ese valor se modifica por un promedio ponderado de los resultados de las iteraciones anterior y actual: : es el coeficiente de relajación que tiene un valor entre 0 y 2. Si  = 1  el resultado no se modifica. Si 0 <  < 1  el resultado es un promedio ponderado de xinuevo y xianterior (subrelajación), se usa para hacer que un sistema no convergente, converja o converja más rápido al amortiguar sus oscilaciones. Si 1 <  < 2  se le da una ponderación extra al valor actual (sobrerelajación), acelera la convergencia de un sistema que ya es convergente. También es conocida como sobrerelajación simultánea o sucesiva, SOR.", "width": "800" } 62 Seudo código para el método de Gauss-Seidel con relajaciónSUBROUTINE Gseid(a,b,n,x,imax,es,lambda) for i = 1:n dummy = a(i,i); for j = 1:n a(i,j) = a(i,j)/dummy; end b(i) = b(i)/dummy; sum = b(i); if i ~= j sum = sum - a(i,j)*x(j); x(i) = sum; iter = 1; while iter < maxIteraciones sentinel = 1; for i = 1:n old = x(i); sum = b(i); for j = 1:n if i ~= j sum = sum - a(i,j)*x(j); end x(i) = lambda * sum +(1.0-lambda)*old; if sentinel == 1 AND x(1) ~= 0.0 ea = abs((x(i)-old)/x(i))*100; if ea > es sentinel = 0; iter = iter + 1; if sentinel == 1 break

21 Eliminación de Gauss Número de operaciones de punto flotante para Gauss simple Para un sistema que se hace cada vez más grande, el tiempo de cálculo se incrementa considerablemente La mayor parte del esfuerzo ocurre en el paso de la eliminación. Por lo que se hace necesario hacer más eficiente el procedimiento

22 Desventajas del método de eliminación de GaussDivisión entre cero Durante las fases de eliminación y sustitución es posible que ocurra división entre cero También se pueden presentar problemas cuando el coeficiente es muy cercano a cero Para evitar estos problemas se utiliza una técnica de pivoteo La normalización de la primera fila involucra una división entre cero, a11 = 0

23 Desventajas del método de eliminación de GaussErrores de redondeo Debido a que las computadoras manejan sólo un número limitado de cifras significativas, pueden ocurrir errores de redondeo y se deben considerar al evaluar los resultados Estos errores pueden ser importantes para sistemas con un gran número de ecuaciones Debido a que cada resultado depende del anterior, el error de los primeros pasos tiende a propagarse Una regla general es la de suponer que los errores de redondeo son importantes cuando n ≥ 100 Siempre se debe sustituir los resultados en las ecuaciones originales y verificar si ha ocurrido un error sustancial

24 Desventajas del método de eliminación de GaussSistemas mal condicionados Sistemas bien condicionados son aquellos en los que un pequeño cambio en uno o más coeficientes provoca un pequeño cambio en la solución Sistemas mal condicionados son aquellos en donde pequeños cambios en los coeficientes generan grandes cambios en la solución Es decir un amplio rango de soluciones puede satisfacer las ecuaciones en forma aproximada Los errores de redondeo pueden inducir pequeños cambios en los coeficientes, si el sistema está mal condicionado estos cambios artificiales pueden generar grandes errores en la solución

25 Desventajas del método de eliminación de GaussSistemas mal condicionados Ejemplo Cambiando a21 de 1.1 a 1.05 sustituyendo en las ecuaciones originales Solución Solución La solución es aproximada

26 Desventajas del método de eliminación de GaussSistemas mal condicionados Esta situación se puede caracterizar de forma matemática, escribiendo las ecuaciones en su forma general Arreglando las ecuaciones en un formato de líneas rectas Si las pendientes son casi iguales Un sistema mal condicionado es aquel en el cual su determinante es cercano a cero Si el determinante es exactamente cero, las pendientes son idénticas, y el sistema no tiene solución o hay un número infinito de soluciones

27 Desventajas del método de eliminación de GaussSistemas singulares Son aquellos donde dos o más ecuaciones son iguales En el caso donde dos ecuaciones son iguales se pierde un grado de libertad siendo imposible resolver el problema de n-1 ecuaciones con n incógnitas Tales casos podrían no ser obvios cuando se trabaja con grandes conjuntos de ecuaciones Se hace necesario tener una forma que de manera automática detecte la singularidad del sistema Esto se logra debido al hecho de que el determinante de un sistema singular es cero Durante el proceso de eliminación se chequea si un elemento de la diagonal es cero, al descubrir uno se puede terminar inmediatamente y generar una excepción o mensaje de error

28 Técnicas para mejorar las soluciones del método de eliminación de GaussUso de más cifras significativas La solución más simple para el mal condicionamiento es usar más cifras significativas en los cálculos El uso de la precisión expandida tiene un precio que se eleva en forma de tiempo de cálculo y cantidad de memoria Pivoteo Ocurren problemas de división por cero cuando el coeficiente pivote es cero Cuando el coeficiente pivote es cercano a cero se pueden introducir errores de redondeo, porque su magnitud puede ser muy pequeña al compararla con la de los demás coeficientes Para evitar esto se utiliza el pivoteo parcial

29 Técnicas para mejorar las soluciones del método de eliminación de GaussPivoteo parcial Antes de normalizar cada fila, se determina el mayor coeficiente pivote disponible en la columna que está por debajo del elemento pivote Las filas se intercambian de manera tal que el coeficiente más grande sea el pivote Ventajas del pivoteo parcial Evita la división entre cero Minimiza el error de redondeo

30 Técnicas para mejorar las soluciones del método de eliminación de GaussPivoteo parcial Ejemplo Problema original Con pivoteo a11 es cercano a cero Restando a la 2da ecuación Restando a la 2da ecuación Se restan dos números casi iguales

31 Seudo código para implementar el pivoteo parcialSe puede usar como una subrutina que podría ser llamada directamente después del inicio del primer ciclo de eliminación Aquí se intercambian de forma física las filas. Para grandes matrices esto puede consumir mucho tiempo Lo que se hace es no intercambiar las filas sino que se guarda el orden de los pivotes en un vector, y según este orden se llevan a cabo las operaciones de eliminación y sustitución SUB pivot(a, b, n, k) p = k big = abs(ak,k) DO ii = k+1, n dummy = abs(aii,k) IF (dummy > big) big = dummy p = ii END IF END DO IF (p ≠ k) DO jj = k, n dummy = ap,jj ap,jj = ak,jj ak,jj = dummy dummy = bp bp = bk bk = dummy

32 Técnicas para mejorar las soluciones del método de eliminación de GaussEscalamiento El escalamiento revela si el pivoteo es necesario Si es necesario, se pivotea pero se retienen los coeficientes originales de la ecuación escalada El escalamiento se usa para calcular los valores escalados de los coeficientes que servirán como un criterio de pivoteo Ejemplo Escalamiento como criterio de pivoteo Manteniendo coeficientes originales Problema original Escalando y pivoteando

33 Descomposición LU e Inversión de matricesEl principal atractivo de este método es que el paso de eliminación, que consume tiempo, se puede reformular de tal manera que involucre sólo operaciones sobre los elementos de la matriz de coeficientes, A De esta forma, es muy adecuado para aquellas situaciones donde se debe evaluar muchos vectores {B} El método de eliminación de Gauss puede implementarse como una descomposición LU La descomposición LU proporciona un medio eficaz para calcular la matriz inversa, la cual a su vez permite evaluar la condición de un sistema

34 Descomposición LU Partiendo de un sistema de ecuaciones lineales de la forma, Este se puede ordenar como, El primer paso de la eliminación de Gauss resulta en un sistema con una matriz triángular superior Que puede ser expresada como, Ahora, suponga que existe una matriz triangular inferior con números 1 sobre la diagonal que tiene la siguiente propiedad si esta propiedad se cumple, de las reglas de multiplicación de matrices se obtiene,

35 Descomposición LU Estrategia para resolver el sistemaPaso de descomposición LU: la matriz [A], se factoriza o descompone en matrices triangulares inferior [L] y superior [U] Paso de sustitución: [L] y [U] se usan par determinar una solución {X} para un vector {B}. Este paso consta de dos subpasos: Se determina el vector intermedio {D} resolviendo [L]{D}={B} por sustitución hacia delante, debido a que [L] es una matriz triangular inferior Se determina {X} resolviendo [U]{X}={D} por sustitución hacia atrás

36 Descomposición LU Descomposición LU con base en la eliminación de Gauss Partiendo de una matriz de coeficientes, se llega a una matriz triangular superior Para llegar a esta matriz [U] La matriz triangular inferior que tiene la propiedad requerida para la descomposición LU es se multiplicó la fila 1 por el factor f21 = a21/a11 y restando el resultado a la fila 2 se eliminó a21 se multiplicó la fila 1 por el factor f31 = a31/a11 y restando el resultado a la fila 3 se eliminó a31 se multiplicó la fila 2 por el factor f32 = a32’/a22’ y restando el resultado a la fila 3 se eliminó a32’

37 Descomposición LU Descomposición LU con base en la eliminación de Gauss Paso de descomposición LU: la matriz [A], se factoriza o descompone en matrices triangulares inferior [L] y superior [U] Paso de sustitución: [L] y [U] se usan par determinar una solución {X} para un vector {B}. Este paso consta de dos subpasos: Se determina el vector intermedio {D} resolviendo [L]{D}={B} por sustitución hacia delante Se determina {X} resolviendo [U]{X}={D} por sustitución hacia atrás

38 Seudo código para descomposición LUSUB Ludecomp(a, b, n, tol, x, er) DIM on, sn er = 0 CALL Decompose(a, n, tol, o, s, er) IF er <> -1 THEN CALL Substitute(a, o, n, b, x) END IF END Ludecomp

39 Seudo código para descomposición LUSUB Decompose(a, n, tol, o, s, er) DO i = 1,n oi = i si = ABS(ai,1) DO j = 2,n IF ABS(ai,j)> si THEN si = ABS(ai,j) END DO DO k = 1,n-1 CALL Pivot(a, o, s, n, k) IF ABS(ao(k),k/so(k)) < tol THEN er = -1 PRINT ao(k),k/so(k) EXIT DO END IF DO i = k+1,n factor = ao(i),k/ao(k),k ao(i),k = factor DO j = k+1,n ao(i),j = ao(i),j - factor*a IF ABS(ao(k),k/so(k))< tol THEN END Decompose

40 Seudo código para descomposición LUSUB Pivot(a, o, s, n, k) p = k big = ABS(ao(k),k/so(k)) DO ii = k+1,n dummy = ABS(ao(ii),k/so(ii)) IF dummy > big THEN big = dummy p = ii END IF END DO dummy = op op = ok ok = dummy END Pivot SUB Substitute(a, o, n, b, x) DO i = 2,n sum = bo(i) DO 1,i-1 sum = sum - ao(i),j*bo(j) END DO bo(i) = sum xn = bo(n)/ao(n),n DO i = n-1,1,-1 sum = 0 DO j = i+1,n sum = sum + ao(i),j*xj xi =(bo(i) - sum)/ao(i),i END Substitute

41 Matriz inversa Para una matriz cuadrada [A], hay otra matriz [A]-1 conocida como la inversa de [A], para la cual se cumple, [A] [A]-1 = [A]-1[A] = I La matriz inversa se puede calcular en una forma de columna por columna a partir de vectores unitarios como vector de constantes del sistema de ecuaciones lineales algebraicas Por ejemplo, para determinar la primera columna de la matriz inversa se resuelve el sistema con el vector de constantes B=[1 0 0]T para determinar la segunda columna se usa B=[0 1 0]T y así sucesivamente

42 Matriz inversa La descomposición LU representa la mejor forma para implementar el cálculo de la matriz inversa, ya que una vez obtenida la descomposición LU de la matriz A se puede calcular su inversa resolviendo cada columna con los vectores unitarios como constantes Ejemplo: Determinar la inversa de A = [ ] 1. Descomposición LU de la matriz A LU(A) = [ ] B1 = [1] [0] [0] B2 = [0] [1] [0] B3 = [0] [0] [1] 2. Sustitución, InvA(1) = [0.3325] [ ] [ ] InvA(2) = [0.0049] [0.1429] [0.0027] InvA(3) = [0.0068] [0.0042] [0.0999]

43 Análisis de error y condición del sistemaLa matriz inversa permite determinar si un sistema está mal condicionado, para esto existen 3 métodos: Escalar la matriz de coeficientes [A], de tal manera que el elemento más grande en cada fila sea 1. Si al invertir la matriz escalada existen elementos de la inversa [A]-1 que sean varios ordenes de magnitud mayores que la unidad, es probable que el sistema esté mal condicionado Multiplicar la inversa por la matriz de coeficientes original y verificar que [A][A]-1  I. Si no es así, indica que el sistema está mal condicionado Invertir la matriz inversa y verificar que el resultado está lo suficientemente cercano a la matriz original. Si no es así, indica que el sistema está mal condicionado

44 Número de condición de una matrizEste número mide la sensibilidad de la solución de un sistema de ecuaciones lineales a errores en los datos Valores cercanos a 1 indican que el sistema está bien condicionado Valores grandes indican que la matriz es casi singular

45 Número de condición de una matrizNorma Euclidiana - Para vector - Para matriz Norma columna-suma Norma 2, o normal espectral max es el eigenvalor más grande de [A]T[A]. Esta es la la norma mínima, por lo tanto proporciona la medida de tamaño más ajustada Norma infinita / matriz uniforme / fila-suma Norma de Frobenius

46 Matrices especiales Matrices banda Matrices simétricasUna matriz banda es una matriz cuadrada en la que todos sus elementos son cero, con excepción de una banda centrada sobre la diagonal principal Las dimensiones de un sistema de banda se pueden cuantificar con dos parámetros: Ancho de banda, BW Ancho de media banda, HBW Estos se relacionan por BW = 2HBW + 1 En general un sistema de banda es aquel para el cual aij = 0 si |i – j| > HBW HBW+1 Diagonal HBW BW

47 Matrices especiales La eliminación de Gauss o la descomposición LU pueden emplearse para resolver sistemas de banda, pero si el pivoteo no es necesario resultan ineficientes, porque se utilizaría tiempo y espacio innecesario en el almacenamiento y manejo de ceros Si se sabe de antemano que el pivoteo es innecesario, se pueden desarrollar algoritmos muy eficientes que no involucren los elementos cero fuera de la banda

48 Sistemas tridiagonalesUn sistema tridiagonal (ancho de banda = 3) se puede expresar como Se cambia la notación a  e, f, g b  r para evitar guardar ceros, ahorrando espacio

49 Sistemas tridiagonalesEjemplo B = [ ] A = [ ] f = [ ] g = [ ] e = [ ]

50 Algoritmo de Thomas (TDMA, tridiagonal matrix algorithm)El algoritmo consiste en tres pasos: Descomposición Sustitución hacia delante Sustitución hacia atrás Manteniendo todas las ventajas de la descomposición LU

51 Seudo código para algoritmo de Thomasa) Descomposición DO k = 2, n ek = ek/fk-1 fk = fk - ek*gk-1 END DO b) Sustitución hacia adelante DO k = 2,n rk = rk - ek*rk-1 c) Sustitución hacia atrás xn = rn/fn DO k = n-1,1,-1 xk = (rk - gk*xk+1)/fk

52 Descomposición de CholeskyLa descomposición o factorización de Cholesky expresa una matriz simétrica como el producto de una matriz triangular y su transpuesta A = L∙LT  L: matriz triangular inferior No todas las matrices simétricas se pueden factorizar de esta forma Las matrices que tienen este tipo de factorización son las matrices simétricas definidas positivas. Esto implica que todos los elementos de la diagonal sean positivos y que los elementos fuera de la diagonal no sean muy grandes

53 Descomposición de CholeskyLos términos de la descomposición se pueden multiplicar entre si. El resultado se puede expresar en forma simple por relaciones recurrentes Para la fila k

54 Descomposición de CholeskyEjemplo, Para k = 1 Para k = 2 Para k = 3 Sustitución hacia adelante Sustitución hacia atrás

55 Seudo código para la descomposición de Choleskyfor k = 1:n for i = 1:k-1 sum = 0; for j = 1:i-1 sum = sum + A(i,j)*A(k,j); end A(k,i) = (A(k,i) - sum)/A(i,i); for j = 1:k-1 sum = sum + A(k,j)^2; A(k,k) = sqrt(A(k,k) - sum);

56 Método de Gauss-SeidelEste es un método iterativo Dado un conjunto de ecuaciones, AX = B Si los elementos de la diagonal son diferentes de cero, se puede resolver la ecuación i para la variable i, donde i = 1…n Se puede empezar el proceso de solución al escoger los valores iniciales de las variables x (xi = 0)

57 Método de Gauss-SeidelLos valores iniciales se sustituyen en la primera ecuación para calcular un nuevo valor para x1 Este nuevo valor de x1 junto con los demás valores iniciales se sustituyen en la segunda ecuación para calcular un nuevo valor para x2 Este proceso se repite hasta calcular los nuevos valores de las n variables Después se regresa a la primera ecuación y se repite todo el procedimiento hasta que la solución converja a la solución real La convergencia se puede verificar usando el criterio,

58 Criterio de convergencia del método de Gauss-SeidelEste método es similar en esencia al método de iteración de punto de fijo que se usa para el cálculo de raíces de una ecuación Presenta las mismas desventajas: En algunos casos no converge En algunos casos la convergencia es lenta Las condiciones suficientes para la convergencia de dos ecuaciones no lineales también aplican para ecuaciones lineales cuando se usa Gauss-Seidel

59 Criterio de convergencia del método de Gauss-SeidelEn el caso de dos ecuaciones el método de Gauss-Seidel se expresa como Las derivadas parciales de estas ecuaciones con respecto a las variables son Para que se cumplan las condiciones suficientes de convergencia

60 Criterio de convergencia del método de Gauss-SeidelEl valor absoluto de la pendiente de las ecuaciones rectas debe ser menor que la unidad para asegurar convergencia Rerformulando, El elemento diagonal debe ser mayor que el elemento fuera de la diagonal para cada fila (sistemas diagonal dominantes) Generalizando para n ecuaciones El criterio es suficiente pero no necesario para convergencia

61 Mejoras a la convergencia por medio de relajaciónLa relajación representa una ligera modificación al método de Gauss-Seidel y está diseñada para mejorar la convergencia Después de calcular cada nuevo valor de x, ese valor se modifica por un promedio ponderado de los resultados de las iteraciones anterior y actual: : es el coeficiente de relajación que tiene un valor entre 0 y 2 Si  = 1  el resultado no se modifica Si 0 <  < 1  el resultado es un promedio ponderado de xinuevo y xianterior (subrelajación), se usa para hacer que un sistema no convergente, converja o converja más rápido al amortiguar sus oscilaciones Si 1 <  < 2  se le da una ponderación extra al valor actual (sobrerelajación), acelera la convergencia de un sistema que ya es convergente. También es conocida como sobrerelajación simultánea o sucesiva, SOR

62 Seudo código para el método de Gauss-Seidel con relajaciónSUBROUTINE Gseid(a,b,n,x,imax,es,lambda) for i = 1:n dummy = a(i,i); for j = 1:n a(i,j) = a(i,j)/dummy; end b(i) = b(i)/dummy; sum = b(i); if i ~= j sum = sum - a(i,j)*x(j); x(i) = sum; iter = 1; while iter < maxIteraciones sentinel = 1; for i = 1:n old = x(i); sum = b(i); for j = 1:n if i ~= j sum = sum - a(i,j)*x(j); end x(i) = lambda * sum +(1.0-lambda)*old; if sentinel == 1 AND x(1) ~= 0.0 ea = abs((x(i)-old)/x(i))*100; if ea > es sentinel = 0; iter = iter + 1; if sentinel == 1 break