Solucion De Sistemas De Ecuaciones Lineales. Metodos Directos E Iterativos 723t36

Solución de Sistemas de Ecuaciones Lineales.

VIVIANA MARCELA BAYONA CARDENAS

UNIVERSIDAD INDUSTRIAL DE SANTANDER INGENIERIA DE PETROLEOS 2010 UNIVERSIDAD INDUSTRIAL DE SANTANDER INGENIERIA DE PETROLEOS

Pág.1

Solución de Sistemas de Ecuaciones Lineales.

1.1. Consideraciones generales 1.1.1 Introducción El objetivo de este tema es introducir al lector en la resolución de sistemas lineales de ecuaciones por métodos gaussianos. Parece conveniente, en primer lugar, establecer la notación y algunas de las bases de álgebra lineal necesarias para alcanzar el objetivo planteado. Siguiendo la notación introducida por Householder en 1964, se empleará: mayúsculas en negrita A, L, U, ∆ , Λ para las matrices, minúsculas con subíndices aij, lij, uij, δ ij, λ ij para los coeficientes de matrices, minúsculas en negrita x, y, z, b, c, d para los vectores letras griegas en minúscula α, β , γ , δ , θ , µ para los escalares. El espacio vectorial de las matrices reales m x n se denota por Rmxn; un elemento cualquiera de ese espacio, A ∈Rmxn , es una matriz rectangular de m filas y n columnas que puede escribirse como:

A=

a11 a21 . am-1,1 am,1

a12 a22 . am.1,2 am,2

...... ...... . . ......

a1,n-1 a2,n-1 . . .

a1,n a2,n . am-1,n am,n

Si m = n entonces A es cuadrada y se dice que tiene orden n. De la misma forma, C mxn es el espacio vectorial de las matrices de coeficientes complejos. En cuanto a los vectores, estos se asumen como “vectores columna”, es decir x ∈ Rn es:

x=

x1 x2 . xn-1 xn

UNIVERSIDAD INDUSTRIAL DE SANTANDER INGENIERIA DE PETROLEOS

Pág.2

donde las componentes xi son números reales. Puesto que por convención se toman los vectores como columna, objetos del tipo (x1, x2 ... xn-1, xn,) son vectores fila y se denotan por xT (T indica matriz o vector traspuesto). Además de las operaciones inherentes al espacio vectorial (suma interna y producto exterior por reales) conviene resaltar por su importancia el producto escalar de vectores. Si x e y son dos vectores de Rn, entonces xT.y = x1 y1 +x2 y2+...+xn.yn ∈ R. El producto escalar de vectores permite definir una métrica: la norma euclídea de x, que es simplemente : x =

x T .x

Del mismo modo que se ha definido xTy, también se puede definir xyT . Sin embargo, el significado de este último producto entre vectores es radicalmente distinto. Sean x e y dos vectores no nulos; la matriz xyT es de Rnxn y se escribe como :

A=

x1.y1 x2.y1 . xn-1.y1 xn.y1

x1y2 x2y2 . xn-1.y2 xn.y2

...... ...... . ...... ......

x1.yn-1 x2.yn-1 . x n-1..yn-1 x n.yn-1

x1.yn x2.yn . x n-1..yn xn.yn

y todas sus columnas (y filas) definen vectores paralelos, es decir, elementos de un espacio vectorial de dimensión uno. Por consiguiente, esta matriz es de rango uno. En realidad, toda matriz de rango uno puede expresarse como el producto de dos vectores de la forma xyT. Estas matrices son comunes en métodos numéricos y conviene saber trabajar con ellas. Por ejemplo, su almacenamiento no se hace guardando todos los coeficientes de la matriz, lo que implicaría almacenar n2 números reales; sino que se almacena únicamente las componentes de cada uno de los vectores x e y, es decir 2n números reales. Para tener una idea del ahorro computacional que esto representa basta suponer que n = 1000: mientras almacenar x e y sólo necesita de 2000 variables reales, xyT requiere un millón (es decir, 500 veces más memoria).

1.1.2 Planteamiento general A lo largo de este tema se plantea la resolución de sistemas lineales de ecuaciones del tipo: Ax = b (1.1.) donde A es una matriz de n x n coeficientes reales ai,j con i =1....n, j = 1...n. Por otra parte, UNIVERSIDAD INDUSTRIAL DE SANTANDER INGENIERIA DE PETROLEOS

Pág.3

b es el vector de términos independientes, también de coeficientes reales bi ,con i =1....n. Finalmente, xT=(x1.......xn) es el vector solución del sistema La existencia y unicidad de soluciones del sistema definido en (1.1.) es por fortuna un tema ampliamente estudiado en el álgebra lineal. Precisamente, el álgebra lineal proporciona una serie de condiciones que permiten verificar si Ax=b tiene solución y si esta es única. De estas condiciones resultan particularmente significativas: • Si las filas (columnas) son linealmente independientes, entonces el sistema A.x=b tiene solución única. • En este caso existe A-1, matriz inversa de A, tal que AA-1 = A-1A = 1 (1 matriz identidad de orden n). • La verificación de la condición de independencia lineal de filas y columnas también implica que para cualquier x ∈ Rn si Ax = 0 , entonces x = 0. •

Si el sistema tiene solución única entonces se verifica que det(A) = A≠ 0

A pesar de la indudable importancia de todas estas condiciones, en el ámbito de la resolución numérica de sistemas de ecuaciones deben ser empleadas con sumo cuidado. De acuerdo a las condiciones enunciadas, el sistema de ecuaciones (1.1.) tiene solución única si y sólo si det(A) = A≠ 0. En este caso, existirá la matriz inversa de A, A-1, que permite escribir la solución del sistema de ecuaciones como : x = A-1.b

(1.2)

Sin embargo, a pesar de tener todo el fundamento analítico necesario, calcular mediante este procedimiento la solución de (1.1) es el peor método posible desde un punto de vista numérico. Para darse cuenta de ello, basta observar sólo dos de los problemas que plantea. En primer lugar, para determinar la unicidad de la solución, debería realizarse el cálculo del determinante. Sin embargo, este trabajo puede ser bastante tedioso puesto que el determinante puede variar bruscamente con pequeños escalados de la matriz. Obsérvese, por ejemplo, que si A es de orden n, entonces det(γ A) = yn.det(A). Para ver las implicancias que esta igualdad impone bastaría tomar el caso particular de n = 100, (número de ecuaciones pequeño hoy en día), entonces: det(0.1 A) = 10-100.det(A). Es decir, dividiendo los coeficientes de A por diez, se reduce el determinante de A en un factor de 10-100. UNIVERSIDAD INDUSTRIAL DE SANTANDER INGENIERIA DE PETROLEOS

Pág.4

En síntesis, se hace muy difícil determinar numéricamente si el determinante de una matriz es realmente nulo por lo cual el uso del mismo se centra básicamente en estudios teóricos. En segundo lugar, el cálculo de la inversa de A (que presenta serios problemas asociados al almacenamiento de la matriz y a la precisión con la que obtengan los resultados), no se emplearía ni en el caso escalar (n = 1). Por ejemplo, para resolver 15x = 3 no se evaluaría primero c = 1/15 para después calcular x como x = 3c. Lo más lógico sería dividir directamente 3 por 15, x = 3/15, lo que permitiría ahorrarse una operación y un error de almacenamiento. Esta situación puede extrapolarse al caso de orden n donde la diferencia en número de operaciones es muy considerable y además los errores de redondeo pueden dar lugar a inestabilidades numéricas. 1.1.3. Resolución algebraica: MÉTODO DE CRAMER En este punto se plantea la resolución analítica de problemas muy pequeños siguiendo un posible enfoque algebraico clásico, conocido como método de Cramer. En este procedimiento, la obtención de los valores de cada una de las incógnitas que conforman el vector x, puede obtenerse a partir de la siguiente expresión:

xi =

a11

...

a1, i −1

b1

a1.i +1

...

a1, n

a 21

....

a 2, i −1

b2

a 2, i +1

...

a 2, n

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

an 1

....

an , i −1

bn A

an , i +1

...

an , n

Es interesante evaluar, para el procedimiento implícito en la ecuación anterior, el número de operaciones elementales (sumas, productos y divisiones) necesarias para obtener la solución del sistema. En primer lugar hay que evaluar n + 1 determinantes y luego realizar n divisiones. Para el cálculo de los determinantes, una de las posibles técnicas necesita de n!.n multiplicaciones y n!-1 sumas. Por consiguiente, el método de Cramer necesita de: (n + 1) (n! - 1) sumas (n + 1) n! n productos n divisiones Cada operación elemental puede tener un coste computacional distinto (por ejemplo, UNIVERSIDAD INDUSTRIAL DE SANTANDER INGENIERIA DE PETROLEOS

Pág.5

muchos ordenadores dividen empleando el método de Newton para ceros de funciones). A pesar de ello, aquí se les asignará el mismo costo computacional a todas las operaciones elementales puesto que ello permite realizar las comparaciones pertinentes. De acuerdo a lo indicado, el número de operaciones elementales con el método de Cramer es TC = (n + 1)2 n! - 1. A modo de ejemplo, la tabla siguiente muestra los valores de Tc para diferentes tamaños del sistema de ecuaciones : Operaciones elementales del método de Cramer n

TC

5

4319

10

4 x 108

100

10158

Los números presentados en la tabla anterior adquieren un mayor relieve cuando se asocian al tiempo necesario para efectuarlos. Si se dispusiera de un superordenador capaz de realizar 100 millones de operaciones en coma flotante por segundo (100 Mflops), el sistema de n = 100 necesitaría aproximadamente 3x10142 años para ser resuelto. Es evidente que el número de operaciones asociado a este método hace prohibitivo su uso, aún para sistemas pequeños. Si además se tiene en cuenta que el ENIAC (primer ordenador digital, fabricado en 1940) realizaba sólo 300 operaciones por segundo y tenía un tiempo medio entre averías de 12 horas, se comprenderá por qué la resolución de sistemas lineales de ecuaciones está en el origen del desarrollo de los métodos numéricos.

1.1.4 Resolución numérica: un enfoque global La estrategia y metodología que se aplica a la resolución numérica de sistemas lineales de ecuaciones parte de una filosofía distinta a la expuesta anteriormente. La regularidad de la matriz A no se determina por un cálculo previo de su determinante. Esta regularidad puede, en algunos problemas, ser evaluada previamente en función de su origen (por ejemplo cuando proviene de la discretización de ecuaciones diferenciales) o a partir de propiedades fácilmente computables como la dominancia diagonal. Por lo general se aplica alguno de los métodos de resolución que se verán seguidamente UNIVERSIDAD INDUSTRIAL DE SANTANDER INGENIERIA DE PETROLEOS

Pág.6

sin evaluar previamente el determinante; el cual, en muchas ocasiones (al igual que la inversa) son un subproducto de los cálculos efectuados. En realidad los algoritmos eficaces para la resolución de sistemas lineales plantean procesos con un enfoque radicalmente distinto, fundamentalmente a partir de considerar que los cálculos se realizan en ordenadores digitales. Por consiguiente, es lógico que los algoritmos se evalúen en función de su eficacia y siguiendo criterios directamente relacionados con su implementación en computadoras.. Existen tres criterios fundamentales para analizar esos algoritmos: •

Número de operaciones necesarias, íntimamente ligado al tiempo de U. Se tendrán en cuenta las operaciones elementales entre números en coma flotante (flop): +, - , /, ó *, todas a un mismo coste computacional aunque no sea exactamente cierto. El número de operaciones es obviamente un excelente indicador, aunque no debe tomarse en sentido estricto como una evaluación del tiempo de cálculo. De hecho, multiplicar el tiempo necesario para una operación por el número de operaciones siempre infravalora el tiempo necesario del algoritmo. Además del tiempo invertido en efectuar las operaciones hay una sobrecarga, debido a la gestión de la memoria, al manejo de los índices enteros, a las instrucciones lógicas en los propios bucles, etc. Aún así, el número de operaciones es un buen indicador del orden del tiempo de U porque esta sobrecarga es generalmente proporcional al número de operaciones, de forma que, aunque no se pueda predecir exactamente el tiempo final, se puede saber cómo varía (linealmente, cuadráticamente, etc.) al modificar, por ejemplo, el orden n de la matriz.

•

Necesidades de almacenamiento, que inciden clara y directamente en las limitaciones de la memoria de los diversos ordenadores; los diferentes métodos de resolución requieren almacenar las matrices de distinta forma en el ordenador y esto varía considerablemente las necesidades de memoria.

•

Rango de aplicabilidad: no todos los métodos sirven para cualquier matriz no singular; además, en función del método y de las propiedades de la matriz, la precisión de los resultados puede verse afectada dramáticamente. Como se verá más adelante, pequeños errores de redondeo pueden producir errores en la solución numérica completamente desproporcionados. No se debe olvidar que debido al enorme número de operaciones necesarias para la resolución de un sistema de ecuaciones de tamaño medio a grande, el análisis estándar de propagación de errores de redondeo no es en absoluto trivial. Conviene resaltar que cada uno de estos criterios puede ser determinante para rechazar un algoritmo. Por ejemplo, para un tipo de ordenador dado, métodos que impliquen UNIVERSIDAD INDUSTRIAL DE SANTANDER INGENIERIA DE PETROLEOS

Pág.7

exceder la memoria disponible son inutilizables por muy rápidos y precisos que resulten. Por lo tanto, el desarrollo de los algoritmos que se plantean a continuación debe tener presentes estos tres criterios en forma simultánea. Desde un punto de vista general las matrices más usuales en las ciencias aplicadas y en ingeniería pueden englobarse en dos grandes categorías: •

Matrices llenas, pero no muy grandes. Por llenas se entiende que poseen pocos elementos nulos y por no muy grandes que el número de ecuaciones es de unos pocos miles a lo sumo. Estas matrices aparecen en problemas estadísticos, matemáticos, físicos e ingenieriles.

•

Matrices vacías (esparzas) y muy grandes. En oposición al caso anterior, vacías indica que hay pocos elementos no nulos y además están situados con una cierta regularidad. En la mayoría de estos casos el número de ecuaciones supera los miles y puede llegar en ocasiones a los millones. Estas matrices son comunes en la resolución de ecuaciones diferenciales de problemas de ingeniería.

Parece lógico que los métodos para resolver sistemas lineales de ecuaciones se adecuen a las categorías de matrices anteriormente expuestas. En general los métodos directos se aplican al primer tipo de matrices, mientras que los métodos iterativos se emplean con el segundo grupo. Es importante observar que no existen reglas absolutas y que todavía en la actualidad existe cierta controversia sobre los métodos óptimos a aplicar en cada caso. En particular, la distinción establecida entre matrices llenas y vacías depende en gran medida del ordenador disponible (fundamentalmente de la memoria). De hecho, los límites han ido evolucionando a lo largo de los años a medida que también evolucionaban los ordenadores (cada vez más rápidos y con más memoria, o al menos más barata). A pesar de ello, casi nadie recomendaría métodos iterativos para matrices llenas con pocas ecuaciones; en cambio, algunos autores trabajan con métodos directos altamente sofisticados y particularizados al entorno informático disponible para resolver sistemas con varios millones de ecuaciones. Observación: En todo lo que sigue se supone que A y b son de coeficientes reales. Si los elementos de A o b son complejos, aparte de las generalizaciones de los métodos que aquí se estudian o de los algoritmos específicos para este caso, se puede replantear el problema como un sistema lineal, con matriz y término independiente reales, de 2n ecuaciones e incógnitas. Para ello se escriben la matriz y los vectores de la siguiente manera: UNIVERSIDAD INDUSTRIAL DE SANTANDER INGENIERIA DE PETROLEOS

Pág.8

A=C+iD b=e+id x=y+iz donde C y D son matrices reales nxn, y c, d, y, z son de Rn. El sistema lineal de ecuaciones original se escribe ahora como: C D

-D C

y z

=

c d

que, resuelto, da el resultado deseado.

1.2. MÉTODOS DIRECTOS 1.2.1 Introducción Los métodos directos de resolución de sistemas lineales de ecuaciones son aquellos que permiten obtener la solución después de un número finito de operaciones aritméticas. Este número de operaciones es, obviamente, función del tamaño de la matriz. Si los ordenadores pudieran almacenar y operar con todas las cifras de los números reales, es decir, si emplearan una aritmética exacta, con los métodos directos se obtendría la solución exacta del sistema en un número finito de pasos. Puesto que los ordenadores tienen una precisión finita, los errores de redondeo se propagan y la solución numérica obtenida siempre difiere de la solución exacta. La cota del error, para una matriz y término independiente dados, se asocia por lo general al número de operaciones de cada método. Se pretende, por lo tanto, obtener métodos con el mínimo número de operaciones posible. Otra particularidad de los métodos directos es que siempre conducen, después de ciertas operaciones, a la resolución de uno o varios sistemas con solución inmediata. Es decir, sistemas donde la matriz es diagonal o triangular. Los métodos para sistemas de resolución inmediata son, de hecho, métodos directos. Una clasificación habitual de estos procedimientos de resolución es aquella que considera que los mismos pertenecen a una de las dos categorías siguientes: métodos de eliminación y métodos de descomposición. En el cuadro siguiente se presenta un esquema con la clasificación de los procedimientos de cálculo más característicos.

UNIVERSIDAD INDUSTRIAL DE SANTANDER INGENIERIA DE PETROLEOS

Pág.9

Clasificación de los métodos directos •

Sistemas con solución inmediata

•

Métodos de Eliminación

 Matriz Diagonal, A = D  Matriz Triangular Superior, A = U  Matriz Triangular Inferior, A = L  Método de Gauss  Método de Gauss-Jordan

•

Métodos de Descomposición

•

Métodos de ortogonalización

    

Método de Doolittle, A = LU Método de Crout, A = LU Método de Cholesky, A = L.LT Descomposición generalizada, A=L.D.LT Método de Thomas (A Tridiagonal)

1.2.2 Sistemas con solución inmediata 1.2.2.1. Matriz Diagonal En este caso la matriz A se escribe como: d11 0 ... ... 0

A=D=

0 d22 ... ... ...

... ...  ... ...

... ... ...  ...

0 0 ... ... dnn

y la solución se obtiene directamente xi =

bi d ii

i = 1..... n

Existe solución si todos los términos de la diagonal son no nulos. Además, si se desea evaluar el determinante de A sólo es necesario calcular el producto de todos los términos de la misma. Por último, el número de operaciones necesario es de n divisiones, es decir TD = n UNIVERSIDAD INDUSTRIAL DE SANTANDER INGENIERIA DE PETROLEOS

Pág.10

operaciones elementales.

1.2.2.2. Matriz Triangular Superior

A=U=

u11 0 ... ... 0

u12 u22 ... ... ...

... ...  ... ...

... ... ... un-1,n-1 ...

u1n u2n ... un-1,n un,n

En este caso la solución de la última ecuación es trivial xn = bn/unn. Una vez conocido xn, la penúltima ecuación (la n - 1) sólo tiene una incógnita que se deduce de forma sencilla. Conocidos ahora xn y xn-1 se pasa a la ecuación anterior (la n - 2) y se resuelve para su única incógnita, Xn-2 . Retrocediendo progresivamente se obtiene el algoritmo de sustitución hacia atrás que se escribe de la siguiente forma: x n = b n / u nn n   b − u . x ∑ ij j   i j=i +1  xi =  u ii

i = n −1, n − 2,...... 1

De nuevo la solución existe si todos los términos de la diagonal de U son no nulos. En forma similar a lo indicado para matrices diagonales, el determinante se evalúa multiplicando los términos ubicados en la diagonal. El número de operaciones es ahora: 1 + 2 +.... + (n - 1) = n(n –1)/2 sumas 1 + 2 +.... + (n - 1) = n(n – 1)/2 productos n divisiones por consiguiente T∆ = n2 operaciones elementales.

UNIVERSIDAD INDUSTRIAL DE SANTANDER INGENIERIA DE PETROLEOS

Pág.11

1.2.2.3. Matriz Triangular Inferior

A=L=

l11 l21 ... ... ln1

0 l22 ... ... ...

... ...  ... ...

... ... ... ln-1,n-1 ln, n-1

0 0 ... 0 ln,n

Se aplica un algoritmo similar al anterior que se denomina de sustitución hacia adelante: x1 = b1 / u11 i −1    b − l .x  ∑ i ij j   j=1  xi =  lii

i = 2,......... n

La existencia de solución, el determinante y el número de operaciones se evalúan exactamente como en el caso anterior.

1.2.3 Métodos de eliminación

1.2.3.1. Método de Gauss En el método de eliminación de Gauss el problema original, Ax = b, se transforma mediante permutaciones adecuadas y combinaciones lineales de cada una de las ecuaciones en un sistema de la forma Ux = c donde U es una matriz triangular superior. Este nuevo sistema equivalente al original es de resolución inmediata, sólo es necesario aplicar el algoritmo de sustitución hacia atrás presentado anteriormente . Durante la transformación del sistema original al equivalente con matriz triangular, las operaciones (que sólo dependen de la matriz A) se realizan sobre la matriz y al mismo tiempo sobre el término independiente. Esto constituye una de las grandes ventajas y a la vez inconvenientes de los métodos de eliminación. Si se dispone de una serie de términos independientes, bj j = 1, - - -, m, conocidos de antemano, se efectúan sobre todos ellos, y al mismo tiempo, las operaciones necesarias para reducir el sistema y obtener una serie de vectores cj j= 1, . . . ,m. Por consiguiente, UNIVERSIDAD INDUSTRIAL DE SANTANDER INGENIERIA DE PETROLEOS

Pág.12

se calculan todos los términos independientes a la vez. Posteriormente se resuelve un sistema con matriz triangular U para cada uno de los cj. Si, por el contrario, no se conocen todos los términos independientes al iniciar los cálculos, es necesario recordar todas las transformaciones necesarias para obtener c1 partiendo de b1; seguidamente se repiten todas estas operaciones sobre los demás términos independientes hasta obtener todos los cj deseados. Hoy en día, en la mayoría de los problemas con matrices de tamaño pequeño o medio, esta propiedad de los métodos de eliminación es la determinante para su elección frente a los métodos de descomposición. Otro punto importante que conviene valorar en el método de Gauss es su importante valor pedagógico. Muchos autores denominan de manera genérica métodos gaussianos al conjunto de métodos de eliminación y descomposición, puesto que la mayoría derivan del trabajo original de Gauss escrito en 1823 (que como otros métodos fundamentales del cálculo numérico fue desarrollado mucho antes de la aparición del primer ordenador). Su implementación en un ordenador sigue siendo la más simple, y con pocas modificaciones, como ya se verá, se obtiene el método más general que existe para la resolución de sistemas lineales de ecuaciones. El algoritmo que se presenta a continuación parte de la ecuación Ax=b que se escribirá como: a11(0) a21(0) a31(0)  an1(0)

a12(0) a22(0) a32(0)  an2(0)

a13(0) a23(0) a33(0)  an3(0)

.... .... ....  ....

a1n(0) a2n(0) a3n(0)  ann(0)

x1 x2 x3  xn

=

b1(0) b2(0) b3(0)  bn(0)

donde el superíndice (0) indica coeficiente de la matriz o del término independiente originales. Si a11(0) ≠ 0, se sustrae de todas la ecuaciones, a partir de la segunda fila, la primera ecuación multiplicada por mi1=ai1(0)/a11(0) con i = 2,. .., n. Esto induce el primer sistema equivalente derivado del original donde la primera columna tiene todos los coeficientes nulos exceptuando el primer coeficiente, y el resto de los coeficientes de la matriz y del término independiente se han visto modificados a partir de la segunda fila: a11(0) 0 0 0 0

a12(0) a22(1) a32(1)  an2(1)

a13(0) a23(1) a33(1)  an3(1)

.... .... ....  ....

a1n(0) a2n(1) a3n(1)  ann(1)

x1 x2 x3  xn

con UNIVERSIDAD INDUSTRIAL DE SANTANDER INGENIERIA DE PETROLEOS

Pág.13

=

b1(0) b2(1) b3(1)  bn(1)

a i(1, j)

bi(1)

= a i(,0j)

=

− mi,1.a1(,0j)

bi(0)

− mi,1.b1(0)

= a i(,0j)

=

bi(0)

a i(,01) (0) − a1, j a1(,01)

a i(,01) (0) − b1 a1(,01)

i = 2,..... n   j = 2...... n

Ahora, si a22(1) (que ya no coincide con el coeficiente que originalmente se encontraba en su posición, a22(0)) es distinto de cero, se sustrae de todas la ecuaciones siguientes la segunda ecuación multiplicada por mi2 = ai2(1)/a22(1) con i = 3,....,n. Después de realizadas estas operaciones sobre el sistema anterior se obtiene el segundo sistema equivalente al original : a11(0)

a12(0)

a13(0)



a1n(0)

x1

b1(0)

0 0 0 0

a22(1) 0  0

a23(1) a33(2)  an3(2)

   

a2n(1) a3n(2)  ann(2)

x2 x3  xn

b2(1) b3(2)  bn(2)

=

donde la segunda columna a partir de la tercera fila sólo contiene términos nulos y los nuevos coeficientes de la matriz y término independiente se obtienen con las siguientes ecuaciones:

a i(,2j)

= a i(1, j)

− mi, 2 .a (21,)j

= a i(1, j)

a i(1, 2) (1) − a 2, j a (21,)2

a i(1, 2) (1) bi( 2) = bi(1) − mi, 2 .b (21) = bi(1) − b2 a (21,)2

i = 3,..... n   j = 3...... n

Cada paso conduce a un nuevo sistema equivalente al original con la particularidad de que la k-ésima matriz es triangular superior si sólo se miran las primeras k ecuaciones y k incógnitas.

En general, se escribe como : UNIVERSIDAD INDUSTRIAL DE SANTANDER INGENIERIA DE PETROLEOS

Pág.14

a11(0) 0  0 0  0

a12(0) a22(1)  0 0  0

      

a1,K(0) a2,K(1)  aK,K(K-1) 0  0

a1,K+1(0) a2,K+1(1)  aK,K+1(K-1) aK+1,K+1(K)  an,K+1(K)

      

a1n(0) a2n(1)  aKn(K-1) aK+1,n(K)  an,n(K)

x1 x2  xK xK+1  xn

=

b1(0) b2(1)  bK(K-1) bK+1(K)  bn(K)

que se ha obtenido a partir de las siguientes ecuaciones

a i(,kj)

bi( k )

= a i(,kj−1)

=

bi( k −1)

− mi, k .a (kk, j−1)

− mi, k .b(kk −1)

= a i(,kj−1)

= bi( k −1)

a i(,kk−1) ( k −1) − a k, j a (kk, k−1)

a i(,kk−1) ( k −1) − bk a (kk, k−1)

i = k + 1,..... n   j = k + 1...... n

Obsérvese que al pasar del (k-1)-ésimo al k-ésimo sistema es necesario realizar las siguientes operaciones: (n - k)(n - k + 1) (n - k)(n - k + 1) n-k

sumas productos divisiones

Finalmente, al deducir el (n - 1)-ésimo sistema se obtiene una matriz triangular superior: a11(0) 0  0 0

a12(0) a22(1)  0 0

    

  

a1n(0) a2n(1) 

an-1,n-1(n-2) 

an-1,n(n-2) an,n(n-1)

x1 x2  xn-1 xn

b1(0) b2(1)  bn-1(n-2) bn(n-1)

A cada uno de los términos que aparecen en la diagonal de la matriz anterior se le denomina pivote. Conviene resaltar que los pivotes no coinciden con los términos originales de la diagonal de A ; es decir, aK,K(K-1) ≠ aK,K(0) para k=1,......,n. Para desarrollar todos los pasos necesarios para obtener el sistema anterior, es necesario suponer que todos los pivotes son no nulos. Es decir, ai,i(i-1) ≠ 0, i = 1,...,n. UNIVERSIDAD INDUSTRIAL DE SANTANDER INGENIERIA DE PETROLEOS

Pág.15

A continuación se presenta el algoritmo que permite obtener la matriz y el término independiente del sistema anterior. a i(,kk−1) ( k −1) a i(,kj) = a i(,kj−1) − mi, k .a (kk, j−1) = a i(,kj−1) − a k, j a (kk, k−1)

a i(,kk−1) ( k −1) (k ) ( k −1) ( k −1) ( k −1) b i = bi − mi, k .b k = bi − bk a (kk, k−1)

k = 1,....., n − 1  i = k + 1,..... n   j = k + 1...... n

donde los términos de superíndice (0) son iguales a los originales del sistema de ecuaciones. El sistema triangular obtenido es de resolución inmediata (Véase el apartado 1.2..2.2). El número de operaciones necesarias para realizar esta primera fase de eliminación ha sido de : n −1

∑(n − k )( n − k +1) =

n ( n 2 −1) 3

sumas

∑(n − k )( n − k +1) =

n ( n 2 −1) 3

productos

k =1 n −1 k =1 n −1

∑(n − k ) =

k =1

n (n −1) 2

divisiones

Si se tienen en cuenta las operaciones correspondientes a la segunda fase de sustitución hacia atrás, el número total de operaciones elementales necesarias para el método de Gauss es: TG =

4n 3 + 9n 2 − 7 n 6

La tabla siguiente muestra el número de operaciones elementales para distintos tamaños del sistema de ecuaciones. Obviamente, se ha obtenido una importante reducción al disponer ahora de un método que crece con n3 en vez de (n! n2) correspondiente al método de Cramer.

UNIVERSIDAD INDUSTRIAL DE SANTANDER INGENIERIA DE PETROLEOS

Pág.16

Operaciones elementales del método de Gauss sin pivotamiento n

TG

5

115

10

805

100

681 550

1000

6.68 x 108

1.2.3.2. Técnicas de Pivoteo Como ya se ha comentado, se ha supuesto a lo largo de todo este análisis que los pivotes eran distintos de cero. Si durante el proceso de eliminación se obtiene un pivote nulo, por ejemplo el ak,k(k-1), se debe buscar en la parte inferior de la columna k-ésima un coeficiente no nulo, es decir de entre los ai,k(k-1) i=k+1,..., n se toma uno que sea distinto de cero. Se sustituye entonces la fila k (y su término independiente) por la fila i (y su término independiente) que se haya escogido. Si dicho coeficiente no nulo no existiera, la matriz sería singular. Esta permutación de filas no sólo tiene interés cuando el pivote es exactamente cero. Es obvio que valores pequeños del pivote pueden producir grandes errores de redondeo, ya que siempre se divide por el valor del pivote. Por consiguiente, para reducir los errores de redondeo conviene escoger el pivote máximo en valor absoluto. Para ello, hay dos técnicas posibles: 1. En el k-ésimo sistema se toma como pivote el coeficiente mayor en valor absoluto de la columna k situado por debajo de la fila k inclusive. Para ello es necesario permutar las filas k y la correspondiente al pivote escogido en la matriz y su término independiente. Esta técnica se denomina método de Gauss con pivotamiento parcial 2. En el k-ésimo sistema, se toma como pivote el coeficiente mayor en valor absoluto de la submatriz de orden n -k definida por los coeficientes que quedan por debajo de la fila k y a la derecha de la columna k. Para ello, se permuta la fila k (y el término independiente asociado) y la columna k con las correspondientes al coeficiente que cumple la condición citada. Al final del proceso deben ponerse en el orden inicial las componentes del vector solución, puesto que su posición ha sido modificada al realizar las permutaciones de columnas. Esta técnica es el método de Gauss con pivotamiento total. Estas dos últimas estrategias producen métodos numéricamente estables. El método de Gauss sin pivotamiento no es necesariamente estable. El estudio detallado de la estabilidad y propagación de errores de redondeo del método de Gauss no es trivial. UNIVERSIDAD INDUSTRIAL DE SANTANDER INGENIERIA DE PETROLEOS

Pág.17

Observación: Si la matriz A es simétrica, las matrices llenas de orden n - k sobre las que se aplica sucesivamente el algoritmo sólo permanecen simétricas si no se realiza ninguna permutación de filas o columnas. La misma observación es válida si la matriz A tiene una estructura que permite el almacenamiento en banda o en perfil. Formalización del algoritmo de Gauss para su programación en un computador. En lo que sigue se utilizará como herramienta para el desarrollo básico de los algoritmos su escritura en un metalenguaje que permita poner en evidencia las estructuras que componen el programa y la organización de sus secuencias. Tal metodología es una de las tantas que pueden utilizarse (diagramas de llaves, diagramas de flujo.. etc) y como tal debe servir para facilitar la construcción de un programa armónico, con una estructura sencilla de interpretar y revisar, en el lenguaje que se elija o se disponga. El algoritmo de Gauss es de muy sencilla implementación a partir de las expresiones detalladas en el punto 1.2.3.1. Como puede verificarse rápidamente, el siguiente desarrollo se ajusta a las expresiones indicadas allí y posee la suficiente generalidad como para su escritura final en un lenguaje de programación cualquiera. Repetir k=1,... n-1 Repetir i=k+1,... n m = aik/akk Repetir j=k+1,..., n aij = aij – m. akj Fin j bi = bi – m.bk Fin i Fin k xn = bn/ann Repetir i=n-1,.... 1 s=0 Repetir j=i+1,... n s = s + aij.xj Fin j xi = (bi-s)/aii Fin i Podrá reconocerse en su seguimiento, una primer parte a partir de la cual se obtiene la UNIVERSIDAD INDUSTRIAL DE SANTANDER INGENIERIA DE PETROLEOS

Pág.18

triangulación de la matriz coeficientes y una segunda en la cual se realiza la sustitución hacia atrás. 1.2.3.4. Matrices Simétricas. Adaptación del Método de Gauss. Ejemplo 1.2. : Resolver el siguiente sistema de ecuaciones por eliminación de Gauss. (como se advertirá, aij = aji) +2 -2 +1 +1 0

-2 +4 -4 +1 +1

+1 -4 +6 -4 +1

+1 +1 -4 +4 -2

0 +1 1 -2 +2

X1 X2 X3 X4 X5

=

+1 +2 +2 +2 +1

Siguiendo el desarrollo clásico del procedimiento de Gauss, resulta: Paso 1 m21 = -2/2 = -1

m31 = +1/2

m41 = +1/2

a21=-2 – (-1).2 =0

a31=+1 – 0.5*2 =0

a41=+1 – 0.5*2 =0

a22 = +4 – (-1).(-2)= +2

a32 = -4 – 0.5.(-2)= -3

a42 = +1 – 0.5.(-2)= +2

a23 = -4 – (-1).(+1)= -3

a33 = +6 – 0.5.(+1)= +5.5

a43 = -4 – 0.5.(+1)= -4.5

a24 = +1 – (-1).(+1)= +2

a34 = -4 – 0.5.(+1)= -4.5

a44 = +4 – 0.5.(+1)= +3.5

a25 = +1 – (-1).(0)= +1

a35 = +1 – 0.5.(0)= +1

a45 = -2 – 0.5.(0)= -2

b2 = +2 – (-1).(+1) = +3

b3 = +2 – 0.5.(+1) = +1.5

b4 = +2 – 0.5.(+1) = +1.5

•

+2

-2

+1

+1

0

+1

0 0 0 0

+2 -3 +2 +1

-3 +5.5 -4.5 +1

+2 -4.5 +3.5 -2

+1 +1 -2 +2

+3 +1.5 +1.5 +1

Notar que la submatriz punteada en rojo sigue siendo simétrica por lo cual, los coeficientes sombreados en amarillo en la tabla de procedimiento podrían no haberse calculado, y asumido directamente los términos simétricos determinados con anterioridad. Este punto constituye la diferencia fundamental que permite optimizar (y de hecho, reducir) las operaciones necesarias para llegar a triangularizar la matriz de coeficientes.

Paso 2 UNIVERSIDAD INDUSTRIAL DE SANTANDER INGENIERIA DE PETROLEOS

Pág.19

m32 = -3/2 = -1.5

m42 = +2/2 = +1

m52 = +1/2

a32=-3 – (-1.5).2 =0

a42= +2 – 1.2 =0

a52= +1 – 0.5.2 =0

a33= +5.5 – (-1.5).(-3) = +1

a43= a34 = -1.5

a53= a35 = +2.5

a34 = -4.5– (-1.5).(+2)= -1.5

a44 =+3.5– 1.(+2)= +1.5

a54 = a45 = -3

a35 = +1– (-1.5).(+1)= +2.5

a45 = -2– 1.(+1)= -3

a55 = +2– 0.5.(+1)= +1.5

b3 = +1.5 – (-1.5).(+3) = +6

b4 = +1.5 – 1.(+3) = -1.5

b5 = +1 – 0.5.(+3) = -0.5

+2 0 0 0 0

-2 +2 0 0 0

+1 -3 +1 -1.5 +2.5

+1 +2 -1.5 +1.5 -3

0 +1 +2.5 -3 +1.5

+1 +3 +6 -1.5 -0.5

Paso 3 m43 = -1.5/1 = -1.5

m53 = +2.5/1 = +2.5

a43=-1.5 – (-1.5).(1) =0

a53=+2.5 – (2.5).(1) =0

a44= +1.5 – (-1.5).(-1.5) = -0.75

a54= a45 = +0.75

a45 = -3– (-1.5).(+2.5)= +0.75

a55 = +1.5– (2.5).(+2.5)= -4.75

b4 = -1.5 – (-1.5).(+6) = +7.5

b5 = -0.5 – (2.5).(+6) = -15.5

+2 0 0 0 0

-2 +2 0 0 0

+1 -3 +1 0 0

+1 +2 -1.5 -0.75 +0.75

0 +1 +2.5 +0.75 -4.75

+1 +3 +6 +7.5 -15.5

-2 +2 0 0 0

+1 -3 +1 0 0

+1 +2 -1.5 -0.75 0

0 +1 +2.5 +0.75 -4

+1 +3 +6 +7.5 -8

Paso 4 m54 = 0.75/(-0.75) = -1 a54=+0.75 – (-1).(-0.75) =0 a55= -4.75 – (-1).(+0.75) = -4 b5 = -15.5 – (-1).(7.5) = -8

+2 0 0 0 0

Paso 5: Sustitución hacia atrás. x5 = b5/a55 = +2 x4 = (+7.5 –(0.75*2))/(-0.75) = -8 UNIVERSIDAD INDUSTRIAL DE SANTANDER INGENIERIA DE PETROLEOS

Pág.20

x3 = (+6 – ( -1.5*(-8) + 2.5*2))/(+1) = -11 x2 = (+3 – (-3*(-11) +2*(-8) +1*2))/2 = -8 x1 = (+1 – (-2*(-8) +1*(-11) +1*(-8)+0*2))/2 = +2 De lo observado en los puntos anteriores puede concluirse en que es posible prescindir de operar sobre aquellos términos ubicados por debajo de la diagonal principal porque: • o bien son coeficientes que se hacen cero (los pivotes sucesivos) • o bien el resultado a obtener ya se ha calculado previamente y se encuentra almacenado en el coeficiente simétrico por sobre la diagonal. Las conclusiones obtenidas del ejemplo anterior pueden generalizarse mediante el siguiente algoritmo para matrices simétricas: Repetir k=1,... n-1 Repetir i=k+1,... n m = aki/akk Repetir j=i,... n aij = aij – m. akj Fin j bi = bi – m.bk Fin i Fin k xn = bn/ann Repetir i=n-1,.... 1 s=0 Repetir j=i+1,... n s = s + aij.xj Fin j xi = (bi-s)/aii Fin i En donde se han destacado en color rojo las secuencias modificadas respecto del algoritmo de Gauss estándar. Se hace notar que en ningún momento se opera sobre coeficientes ubicados por debajo de la diagonal, situación que hace plantear si es realmente necesario almacenar inicialmente esos coeficientes. Cantidad de Operaciones necesarias Como se indicó anteriormente, para pasar del (k-1)esimo al K-esimo sistema es preciso realizar: UNIVERSIDAD INDUSTRIAL DE SANTANDER INGENIERIA DE PETROLEOS

Pág.21

(n - k)(n - k + 1) sumas (n - k)(n - k + 1) productos n-k divisiones sin embargo, para este caso particular, puede advertirse que, en función de las condiciones de simetría al operar sobre cada una de las n-k últimas ecuaciones, solamente en la primera (la k+1), se realizan (n-k+1) operaciones, en el resto , se continúa con la siguiente secuencia: fila

nº sumas (productos)

k+1

(n-k+1)

k+2

(n-k+1) –1

k+3

(n-k+1)-2

k+4

(n-k+1)-3





n

(n-k+1) – (n-k-1)

lo cual, implica que al pasar del (k-1) esimo sistema al K-esimo , al utilizar las propiedades de simetría de la matriz transformada se realizan: 0+1+2+.......+ (n-k-1) = (n-k)(n-k-1)/2 sumas y productos menos que aquellos que se generarían de no considerar tal propiedad, es decir que el número de operaciones resultantes para este paso sería de: (n-k).(n-k+1) - (n-k)(n-k-1)/2 sumas y productos Si ahora se considera que k= 1,....,(n-1). la cantidad de operaciones realizadas en este caso, resultaría: n −1

(n − k ) ( n − k + 3) 2 k =1

∑

n −1

(n − k ) n − k + 3) 2 k =1

∑

n −1

∑(n − k ) =

k =1

n( n −1) 2

sumas productos divisiones

En la figura siguiente se expone una comparación entre el número de operaciones a realizar al utilizar el método de Gauss estándar y aquel resultante de considerar las propiedades de simetría de la matriz. Como resulta del mismo, ya para sistemas de reducido número de ecuaciones (n>100), la cantidad de operaciones resultantes al considerar la simetría de la matriz de coeficientes es aproximadamente al mitad de UNIVERSIDAD INDUSTRIAL DE SANTANDER INGENIERIA DE PETROLEOS

Pág.22

aquella que resultaría de utilizar el método de Gauss estándar.

Comparación Gauss vs Gauss Simetrico

% Reducción de operaciones

60 50 40 30 20 10 0 0

20

40

60

80 100 120 140 160 180 200

nº de ecuaciones

1.2.3.5. Matrices simétricas y bandeadas. Adaptación del método de Gauss. Un caso de interés particular son las matrices denominadas “bandeadas”. Como se verá más adelante, tales matrices son muy comunes en problemas de ingeniería en los cuales se realiza el planteo numérico de las ecuaciones diferenciales de gobierno de un problema físico determinado. Tales matrices se caracterizan por tener un número muy importante de elementos nulos en relación al número total de coeficientes de la matriz (n 2) y porque además, esos elementos no nulos se ubican a lo ancho de una “banda” paralela a la diagonal de la matriz. Se denomina precisamente “ancho de banda” al número de columnas distribuidas en forma simétrica a la diagonal, por fuera de la cual, los coeficientes de la matriz son nulos. Ejemplo: Matriz “banda” de 8 x 8. Ancho de banda=5 a11 a21 a31 0

a12 a22 a32 a42

a13 a23 a33 a43

0 a24 a34 a44

0 0 a35 a45

0 0 0 a46

UNIVERSIDAD INDUSTRIAL DE SANTANDER INGENIERIA DE PETROLEOS

Pág.23

0 0 0 0

0 0 0 0

0 0 0 0

0 0 0 0

a53 0 0 0

a54 a64 0 0

a55 a65 a75 0

a56 a66 a76 a86

a57 a67 a77 a87

0 a68 a78 a88

donde, generalmente, aij = aji Si bien sistemas de ecuaciones con matrices de coeficientes de este tipo, denominadas habitualmente esparzas o vacías, resultan a menudo resueltos mediante técnicas iterativas, resulta de interés encarar su solución directa mediante una adaptación del algoritmo enunciado para matrices simétricas, en el punto 1.2.3.4. Esa adecuación surge de advertir que, para cada columna, los elementos ubicados por debajo de la banda ya son nulos, por lo cual, las operaciones de eliminación pueden limitarse al ancho de banda. Esto es, para obtener la matriz triangular superior: Repetir k=1,... n-1 si k+b>n fin =n sino fin =k+b Fin_si Repetir i=k+1,... fin m = aki/akk Repetir j=i,... fin aij = aij – m. akj Fin j bi = bi – m.bk Fin i Fin k en donde “b” es el semiancho de banda. En el algoritmo anterior se agregó una estructura alternativa a fin de determinar correctamente la finalización de cada uno de los lazos repetitivos. Debe notarse que, de la forma en que se ha organizado el procedimiento, la matriz triangular superior obtenida mantiene los elementos no nulos dentro de la “banda” por lo cual, la rutina de sustitución hacia atrás puede re-escribirse como: xn = bn/ann Repetir i=n-1,.... 1 si i+b>n fin =n Sino UNIVERSIDAD INDUSTRIAL DE SANTANDER INGENIERIA DE PETROLEOS

Pág.24

fin =i+b Fin_si s=0 Repetir j=i+1,... fin s = s + aij.xj Fin j xi = (bi-s)/aii Fin i Numero de operaciones El número de operaciones necesarias para la solución de un sistema determinado dependerá, además del orden de la matriz, de la magnitud del ancho de banda. En el gráfico siguiente se muestra la relación existente entre el ancho de banda, expresado como porcentaje del orden de la matriz y la relación entre el número de operaciones contempladas en el algoritmo anterior y el procedimiento de Gauss estándar. Gauss Banda Simétrico vs Gauss Estándar

En la misma se representaron las relaciones mencionadas para distintos 16 valores del orden del sistema de 14 ecuaciones, observándose que, cuando 12 n>200, la eficiencia del algoritmo se 10 hace aproximadamente constante para 8 una cierta relación entre el ancho de 6 banda y el orden de la matriz (b/n). 4 Es importante notar que en sistemas 2 provenientes de la solución numérica 0 de ecuaciones diferenciales en 0 5 10 15 20 25 30 35 40 problemas de ingeniería, es frecuente b/n [%] encontrar relaciones b/n del orden del n=1000 n=100 n=200 10 a 20 %, rango en donde la eficiencia del algoritmo GBS permite realizar entre 20 y 50 veces menos operaciones que Gauss estándar.

Nº Op. GSB/GE [%]

18

1.2.3.6. Método de Gauss-Jordan A continuación se presenta una variante del método de Gauss que conviene considerar. En este procedimiento, además de sustraer la fila k multiplicada por mik = aik(k-1)/akk(k-1) a las filas posteriores, también se sustrae a las anteriores. Es práctica común también en este caso dividir la fila k por el término diagonal a fin de que este quede unitario. UNIVERSIDAD INDUSTRIAL DE SANTANDER INGENIERIA DE PETROLEOS

Pág.25

El procedimiento general se ilustra con el siguiente ejemplo: 1 5 2

-1 -4 1

1 3 1

X1 X2 X3

=

-4 -12 11

Paso 1 m21= a21/a11 =5

m31 = a31/a11 = +2

a21 = 5- 5.1 = 0

a31 = 2- 2.(1) = 0

a22 = -4 – 5.(-1) = +1

a32 = 1 – 2.(-1) = +3

a23 = 3 – 5.1 = -2

a33 = 1 – 2 (1) = -1

b2= -12 –5..(-4) = 8

b3= 11 –2..(-4) = 19 1 0

-1 1

1 -2

X1 X2

0

3

-1

X3

=

-4 8 19

Paso 2 m12= a12/a22 =-1

m32 = a32/a22 = +3

a11 = 1- (-1).0 = 1 a12 = -1- (-1).1 = 0

a32 = 3 – 3.(1) = 0

a13 = 1 – (-1).(-2) = -1

a33 = -1 – 3 (-2) = 5

b1= -4 –(-1).(8) = 4

b3= 19 –3..(8) = -5 1 0 0

0 1 0

-1 -2 5

X1 X2 X3

0 1 0

-1 -2 1

X1 X2 X3

=

4 8 -5

=

4 8 -1

Paso 3: Dividir la 3era ecuación por a33 1 0 0 m13= a13/a33 =-1

m23 = a23/a33 = -2

a11 = 1- (-1).0 = 1 a12 = -0- (-1).0 = 0

a23 = -2 – (-2).(1) = 0

a13 = -1 – (-1).(1) = 0 UNIVERSIDAD INDUSTRIAL DE SANTANDER INGENIERIA DE PETROLEOS

Pág.26

b1= 4 –(-1).(-1) = 3

b2= 8 –(-2).(-1) = 6 1 0 0

0 1 0

0 0 1

X1 X2 X3

=

3 6 -1

Algoritmo de Gauss-Jordan Repetir k=1,... n PIV = akk Repetir j=k,... n ak,j = akj / PIV Fin j bk = bk/PIV Repetir i=1,... n Si i≠ k m = aik Repetir j=k+1,..., n aij = aij – m. akj Fin j bi = bi – m.bk Fin Si Fin i Fin k El número de operaciones de este procedimiento es T GJ = n3 + 1/2 n2 – 5/2 n +1, lo cual significa un 50 % más de operaciones que el método de Gauss estándar cuando n→∞ (n>25).

UNIVERSIDAD INDUSTRIAL DE SANTANDER INGENIERIA DE PETROLEOS

Pág.27

1.2.3.7. Técnicas de Pivoteo Supongamos el siguiente sistema de ecuaciones en donde cada número se representa en base decimal, con coma flotante y cuatro dígitos para la mantisa, el cual se resolverá en base al método de eliminación estándar de Gauss: 0.10 4.63 -3.39

3.01 -1.06 9.81

-4.88 -2.27 -4.78

X1 X2 X3

-7.70 -6.36 3.95

=

Paso 1 m21= a21/a11 =46.3 a21 = 4.63- 46.3x0.1 = 0 a22=-1.06-3.01x46.3 = -1.06-139.4=-140.5 a23 =-2.27-(-4.88x46.3)=-2.27+225.9=223.6 b2= -6.36-(-7.7x46.3)=-6.36+356.5=350.1

0.10 0 0

m31 = a31/a11 = -33.9 a31 = -3.39-(-33.9x0.1) = 0 a32 = 9.81-(-33.9x3.01) = 111.8 a33 = -4.78-(-33.9x(-4.88)) = -170.2 b3= 3.95-(-33.9x(-7.7)) = -257.1

3.01 -4.88 -140.5 223.6 111.8 -170.2

X1 X2 X3

=

-7.70 +350.1 -257.1

=

-7.70 +350.1 +21.5

Paso 2 m32= a32/a22 = -0.7957 a32 = 111.8-(m32x(-140.5))= 0 a33 = -170.2-(m32x(223.6))= -170.2+177.9=7.7 b3= -257.1-(m32x350.1)=-257.1+278.6=21.5

0.10 0 0

3.01 -140.5 0

-4.88 223.6 7.7

X1 X2 X3

X3=+2.792 (Er=2%) X2= +(350.1 –223.6x2.792)/(-140.5) =+(350.1 –624.3)/(-140.5) = +1.952 (Er=4%) X1 = (-7.70 – 3.01x1.952-(-4.88x2.792))/0.10 = (-7.7-5.876+13.62)/0.1 =0.44 (Er=9%) Los errores indicados resultan excesivos para un sistema de ecuaciones de muy reducido tamaño como el considerado. Debe notarse, que esos errores crecen en la medida que algún término diagonal resulte de un orden de magnitud notoriamente diferente al del resto de los coeficientes de la matriz. En efecto, si se rehace el problema anterior pero UNIVERSIDAD INDUSTRIAL DE SANTANDER INGENIERIA DE PETROLEOS

Pág.28

considerando a a11 = 0.01, resulta:

0.01 4.63 -3.39

3.01 -1.06 9.81

-4.88 -2.27 -4.78

X1 X2 X3

-7.70 -6.36 3.95

=

Paso 1 m21= a21/a11 =463

m31 = a31/a11 = -339

a21 = 4.63- 463x0.01 = 0

a31 = -3.39-(-339x0.01) = 0

a22=-1.06-3.01x463 = -1.06-1394=-1395

a32 = 9.81-(-339x3.01) = 1030

a23 =-2.27-(-4.88x463)=-2.27+2259=2257

a33 = -4.78-(-339x(-4.88)) = -1659

b2= -6.36-(-7.7x463)=-6.36+3565=3559

b3= 3.95-(-339x(-7.7)) = -2606

0.01 0 0

3.01 -1395 1030

-4.88 2257 -1659

X1 X2 X3

=

-7.70 3559 -2606

=

-7.70 3559 22

Paso 2 m32= a32/a22 = -0.7384 a32 = 1030-(m32x(-1395))= 0 a33 = -1659-(m32x(2257))= -1659+1667=8 b3= -2606-(m32x3559)=-2606+2628=22 0.01 0 0

3.01 -1395 0

-4.88 2257 8

X1 X2 X3

X3=+2.75 (Er=1%) X2= +(3559 –2257x2.75)/(-1395) =+(3559 –6207)/(-1395) = +1.898 (Er=2%) X1 = (-7.70 – 3.01x1.898-(-4.88x2.75))/0.01 = (-7.7-5.713+13.42)/0.01 =0.7 (Er=78%) Como resultaba de esperar, los errores relativos resultaron elevados, principalmente en el resultado obtenido para X1. Por otra parte, como ya se ha comentado, se ha supuesto a lo largo de todo este análisis UNIVERSIDAD INDUSTRIAL DE SANTANDER INGENIERIA DE PETROLEOS

Pág.29

que los pivotes eran distintos de cero. Si durante el proceso de eliminación se obtiene un pivote nulo, por ejemplo el ak,k(k-1) se debe buscar en la parte inferior de la columna késima un coeficiente no nulo, es decir de entre los al ai,k(k-1) i=k+1,..., n se toma uno que sea distinto de cero. Se sustituye entonces la fila k (y su término independiente) por la fila i (y su término independiente) que se haya escogido. Si dicho coeficiente no nulo no existiera, la matriz sería singular. Esta permutación de filas no sólo tiene interés cuando el pivote es exactamente cero. Es obvio que valores pequeños del pivote pueden producir grandes errores de redondeo, ya que siempre se divide por el valor del pivote. Por consiguiente, para reducir los errores de redondeo conviene escoger el pivote máximo en valor absoluto. Para ello, hay dos técnicas posibles: 3. En el k-ésimo sistema se toma como pivote el coeficiente mayor en valor absoluto de la columna k situado por debajo de la fila k inclusive. Para ello es necesario permutar las filas k y la correspondiente al pivote escogido en la matriz y su término independiente. Esta técnica se denomina método de Gauss con pivotamiento parcial 4. En el k-ésimo sistema, se toma como pivote el coeficiente mayor en valor absoluto de la submatriz de orden n -k definida por los coeficientes que quedan por debajo de la fila k y a la derecha de la columna k. Para ello, se permuta la fila k (y el término independiente asociado) y la columna k con las correspondientes al coeficiente que cumple la condición citada. Al final del proceso deben ponerse en el orden inicial las componentes del vector solución, puesto que su posición ha sido modificada al realizar las permutaciones de columnas. Esta técnica es el método de Gauss con pivotamiento total. Estas dos últimas estrategias producen métodos numéricamente estables. El método de Gauss sin pivotamiento no es necesariamente estable. El estudio detallado de la estabilidad y propagación de errores de redondeo del método de Gauss no es trivial. En base a los procedimientos comentados, se propone resolver el primer ejemplo utilizando una técnica de pivoteo parcial. 0.10 4.63 -3.39

3.01 -1.06 9.81

-4.88 -2.27 -4.78

X1 X2 X3

=

-7.70 -6.36 3.95

Paso 1 (k=1): Se selecciona como pivote al coeficiente de mayor valor absoluto en la columna 1, intercambiándose las filas correspondientes: UNIVERSIDAD INDUSTRIAL DE SANTANDER INGENIERIA DE PETROLEOS

Pág.30

4.63 0.10 -3.39

-1.06 3.01 9.81

-2.27 -4.88 -4.78

X1 X2 X3

-6.36 -7.70 3.95

=

Paso2 m21= a21/a11 =0.0216

m31 = a31/a11 = -0.7322

a21 = 0.10 – 0.0216x4.63 = 0

a31=-3.39 – m31x4.63 = 0

a22=3.01-(-1.06x0.0216)=3.01+0.0229=3.033

a32=9.81-(-1.06xm31)= 9.034

a23=-4.88-(-2.27x0.0216)=-4.88+0.04903=-4.831

a33=-4.78-(-2.27xm31)= -6.442

b2= -7.70-(-6.36x0.0216)=-7.70+0.1374=-7.563

b3= 3.95-(-6.36xm31)= -0.707

4.63 0 0

-1.06 3.033 9.034

-2.27 -4.831 -6.442

X1 X2 X3

=

-6.36 -7.563 -0.707

Paso 3 : Se selecciona como pivote al coeficiente de mayor valor absoluto en la columna 2 , en el rango a22..an2, intercambiándose las filas correspondientes:

4.63 0 0

-1.06 9.034 3.033

-2.27 -6.442 -4.831

X1 X2 X3

-2.27 -6.442 -2.668

X1 X2 X3

=

-6.36 -0.707 -7.563

=

-6.36 -0.707 -7.326

m32= a32/a22 = 0.3357 a32 = 3.033-(0.3357x9.034)=0 a33 = -4.831-(0.3357x(-6.442))= -2.668 b3= -7.563-(0.3357x(-0.707))= -7.326

4.63 0 0

-1.06 9.034 0

X1 = 0. 4030 , X2 = 1.88 , X3 = 2.746 En todos los casos, los errores obtenidos respecto de la solución exacta resultan inferiores al 0.1 %. 1.2.4. Métodos de Descomposición LU

UNIVERSIDAD INDUSTRIAL DE SANTANDER INGENIERIA DE PETROLEOS

Pág.31

Se ha estudiado con anterioridad aquellos procedimientos de eliminación consistentes en la aplicación directa del algoritmo de Gauss. Estos procedimientos, bastante eficientes en el manejo de problemas en donde el o los términos independientes son tratados en forma simultanea con la matriz de coeficientes, presuponen –aunque parezca obvio aclararlo- la “existencia” o conocimiento de cada uno de los vectores que conforman esos términos independientes para los cuales se realiza el cálculo. Existen algunos problemas en ingeniería, tales como aquellos vinculados a procesos de evolución en el tiempo, en los cuales los sucesivos términos independientes a plantear, dependen de los encontrados en pasos anteriores, es decir, que su conocimiento no es simultáneo ni se dispone al inicio del proceso. En razón de ello, si se utilizara la técnica de Gauss tal como se la desarrolló anteriormente, conduciría a un enorme costo computacional, toda vez que habría que proceder a realizar los mismos pasos de eliminación sobre la matriz de coeficientes para cada nuevo vector “instantáneo” de términos independientes. En este sentido, no parece muy razonable repetir una y otra vez los mismos cálculos sobre la matriz de coeficientes, siempre sobre los mismos números para obtener la misma matriz triangular superior, solo con el objeto de realizar los mismos pasos sobre el vector de términos independientes. Para ello bastaría con “recordar” las operaciones que se realizan en el proceso de eliminación para, de esta forma, repetirlas en cada paso sobre cada nuevo término independiente. De esta forma, el mayor costo estaría en el paso inicial en donde se construye una especie de “tabla de memoria” en la cual se asienten las operaciones realizadas, para después, en cada paso de solución repetir las operaciones anotadas en la “tabla” sobre cada nuevo término independiente. Este último paso resultará así mucho más económico que rehacer todas las cuentas en forma innecesaria sobre la matriz, ya que se realiza únicamente sobre los términos del nuevo vector independiente

1.2.4.1.- Descomposición LU- Método de Doolittle El problema a resolver, consiste en encontrar el vector de incógnitas X que satisface la expresión: A.x – b = 0 (1) En la siguiente presentación se supondrá que la matriz A es de 3x3 aunque el desarrollo es del todo general y puede extenderse a cualquier orden del problema. UNIVERSIDAD INDUSTRIAL DE SANTANDER INGENIERIA DE PETROLEOS

Pág.32

Supongamos que la ecuación (1) se puede expresar como un sistema triangular superior, esto es: u11 0 0

u12 u22 0

u13 u23 u33

x1 x2 x3

=

d1 d2 d3

(3)

Podrá reconocerse fácilmente que la expresión anterior coincide con el estado final de la etapa de eliminación seguida en el algoritmo de Gauss, la cual, puede resumirse en: U.x – d = 0

(4)

Ahora bien, supóngase por un momento la existencia de una matriz triangular inferior con coeficientes iguales a 1 sobre su diagonal, esto es: 1 l21 l31

L=

0 1 l32

0 0 1

(5)

la cual, supóngase también, tiene la propiedad que al premultiplicarla por la ecuación (4) su resultado es (1), es decir: L. [U.x – d] = A.x – b

(6)

en la cual, resulta evidente que: L.U = A L.d = b

(7) (8)

de resultar así, el procedimiento de cálculo puede ordenarse en base a la siguiente estrategia: a) Fase de descomposición LU: Consiste en factoriza o descomponer A en el producto de una matriz triangular inferior, L, por una superior, U. b) Fase de sustitución: en la cual se obtiene el vector solución. Para ello esta fase se compone de las siguientes dos etapas: b.1. de la ecuación (8) se obtiene d por sustitución hacia adelante. UNIVERSIDAD INDUSTRIAL DE SANTANDER INGENIERIA DE PETROLEOS

Pág.33

b.2. de la ecuación (4) se obtiene x por sustitución hacia atrás. Todo lo citado es factible de realizar a partir de la existencia de matrices U y L con las propiedades señaladas. Una forma de arribar a determinar sus coeficientes es mediante el mismo procedimiento de eliminación de Gauss ya conocido. Esta posibilidad resulta inmediata a partir de la propia forma de la matriz U, que, como se indicara con anterioridad, puede escribirse a U de la siguiente forma:

U=

a11 0 0

a12 a13 a’22 a’23 0 a”33

(9)

donde los coeficientes señalados con apóstrofes y comillas corresponden a los obtenidos en el primer y segundo paso respectivamente del proceso de eliminación de Gauss. Aunque nada de lo dicho lo sugiere, los coeficientes de la Matriz L se obtienen como “subproducto” del proceso de generación de U. En el ejemplo siguiente se pondrá este punto en evidencia aunque la metodología seguida es de carácter absolutamente general y puede ampliarse a un sistema de cualquier número de ecuaciones. Supongamos el siguiente sistema de 3 ecuaciones con 3 incógnitas:

A.x = b →

a11 a21

a12 a22

a13 a23

x1 x2

a31

a32

a33

x3

=

b1 b2 b3

Paso 1 m21= a21/a11

a11

m31 = a31/a11

a’21= a21 – m21*a11 = 0 a’31= a31 – m31*a11 = 0 a’22= a22 – m21*a12

a’32= a32 – m31*a12

a’23= a23 – m21*a13

a’33= a33 – m31*a13

→

a12

a13

0

a’22 a’23

0

a’32 a’33

Paso 2 m’32= a’32/a’22

a11

a”32= a’32 – m’32*a’22 = 0 →

0

a”33= a’33 – m’32*a’23

0

a12

a13

a’22 a’23 0

a”33

Evidentemente, para no modificar el sistema de ecuaciones original, las mismas transformaciones realizadas sobre la matriz de coeficientes debería realizarse sobre el vector de términos independientes. UNIVERSIDAD INDUSTRIAL DE SANTANDER INGENIERIA DE PETROLEOS

Pág.34

Sin embargo, supongamos que existe la posibilidad de armar una “tabla de memoria” en la cual se registren todas las operaciones a realizar sobre los coeficientes del término independiente. De esta forma, no sería necesario ir operando sobre el vector b en forma simultánea con las transformaciones sobre A, sino que las mismas podrían ser ejecutadas con posterioridad. Un posible formato para confeccionar esa “tabla de memoria” se representa a continuación:

m’21 m’31 m”32 En la misma se almacenan todos los “pivote” de cada una de las transformaciones de fila realizadas. Un punto fácilmente verificable es el hecho que, la matriz siguiente formada a partir de la “tabla de memoria”,

L=

1

0

0

m’21

1

0

m’31 m”32

1

al premultiplicar a la matriz U dada por (9), genera como resultado la matriz A original, es decir: L.U = A que es la base del procedimiento de descomposición comentado anteriormente. Una vez obtenidas las matrices indicadas, se procede a resolver los siguientes sistema de resolución directa: • L.d = b mediante sustitución hacia adelante • U.x = d por sustitución hacia atrás El ejemplo siguiente ilustra el procedimiento completo sobre un sistema de 3x3, aunque la metodología es enteramente aplicable a sistemas de cualquier orden:

A.x = b →

3 1

-1 7

-2 -3

x1 x2

3

-2

10

x3

=

+4

UNIVERSIDAD INDUSTRIAL DE SANTANDER INGENIERIA DE PETROLEOS

+7 -19

Pág.35

a) Fase de Descomposición Paso 1 m21= 1/3 = 0.333

m31 = 3/3 = 1

3

a’21= 1 – 0.333*3 = 0

a’31= 3 – 1*3 = 0

→

-1

0

a’22= 7 – 0.333*(-1) = +7.333 a’32= -2 – 1*(-1)= -1

-2

7.33 -2.33

0

-1

12

a’23= -3 – 0.333*(-2)= -2.333 a’33= 10 – 1*(-2)= 12

Paso 2 m’32= -1/7.33 = -0.136

3 →

a”32= -1 +0.136*7.33 = 0

0

a”33= 12 + 0.136*(-2.33)

L=

-1

-2

7.33 -2.33

0

0

1

0

0

0.333

1

0

1

-0.136

1

11.68 3 ;

U=

-1

0

-2

7.33 -2.33

0

0

11.68

Puede verificarse que L.U = A, excepto pequeñas diferencias producto de errores de redondeo.

b) Fase de Sustitución Como se indicó, el proceso de sustitución comprende dos etapas de resolución inmediata: a)

L.d = b , sustitución hacia adelante

d1 = b1/L11 = +7 d2 = (b2-L21*d1)/L22 = (-19-0.333*7)/1 = -21.33 d3 = (b3-L31*d1-L32*d2)/L33 = (+4-1*7 +0.136*(-21.33))/1 =-5.90 b) U.x = d , sustitución hacia atrás x3 = d3/U33 = -0.51 x2 = (d2-U23*x3)/U22 = -3.07 x1 = (d1-U12*x2-U13*x3)/U11 = +0.97

UNIVERSIDAD INDUSTRIAL DE SANTANDER INGENIERIA DE PETROLEOS

Pág.36

El procedimiento detallado puede especificarse en forma de pseudo-código a fin de facilitar su implementación en un computador: Descomposición LU Repetir k=1,... n-1 Repetir i=k+1,... n aik = aik/akk Repetir j=k+1,..., n aij = aij – aik. akj Fin j Fin i Fin k Sustitución hacia adelante d1 = b1 Repetir i=2,.... n s=0 Repetir j=1,... i-1 s = s + aij.dj Fin j di = bi-s Fin i Sustitución hacia atrás xn = dn/ann Repetir i=n-1,.... 1 s=0 Repetir j=i+1,... n s = s + aij.xj Fin j xi = (di-s)/aii Fin i

Si el algoritmo detallado se implementa para la matriz A del problema anterior, se observa que, luego de finalizada la etapa de descomposición, los coeficientes de la matriz A resultan: 3

-1

-2

0.333 7.33 -2.33 UNIVERSIDAD INDUSTRIAL DE SANTANDER INGENIERIA DE PETROLEOS

Pág.37

1

-0.136 11.68

Como se advierte, en el mismo arreglo se almacenaron los coeficientes de la matriz U (triángulo superior) y de la matriz L (triángulo inferior). En este último caso, el almacenamiento de la diagonal es un asunto trivial toda vez que todos los coeficientes ubicados sobre la misma resultan iguales a 1. Se ha optimizado así el manejo de memoria, limitando el área de trabajo a aquella necesaria solamente para el almacenamiento de la matriz de coeficientes (A). Una segunda optimización puede realizarse realmacenando los coeficientes del vector d en el proceso de sustitución hacia adelante en el mismo espacio que el vector b, cuyos coeficientes ya no resulta necesario almacenar en esa etapa del cálculo. Finalmente, un punto de interés reside en la determinación del número de operaciones realizadas por el algoritmo anteriormente desarrollado. Como puede verificarse, el número total de operaciones a realizar en el proceso de descomposición y sustitución es idéntico a aquel realizado en el algoritmo de eliminación de Gauss. Evidentemente, la gran ventaja de LU surge en el caso de tener que resolver en forma reiterada el mismo sistema pero para términos independientes diferentes. En este caso, cada nuevo sistema requiere realizar solamente un orden de n2 operaciones.

1.2.4.2. Método de Crout En el apartado anterior se desarrollaron las bases del procedimiento propuesto por Doolittle para la descomposición de sistemas de ecuaciones lineales Ax=b , a partir de la factorización de la matriz de coeficientes de la forma A = L.U. Básicamente, el procedimiento enmarcaba un reordenamiento del algoritmo de Gauss, en donde la matriz U es calculada a partir de técnicas de eliminación estándar y la matriz L se obtenía también como subproducto de ese procedimiento. El procedimiento de Crout parte de un enfoque conceptual diferente, no obstante ello, es importante adelantar, que también puede considerarse como una forma diferente de plantear el algoritmo de Gauss y por ello resultan las mismas propiedades, entre ellas, el mismo número de operaciones finales. Como todo método de descomposición posee una ventaja comparativa importante al momento de resolver sistemas de términos independientes múltiples, en donde los mismos no son conocidos en forma simultánea, como resulta característico de problemas iterativos o de evolución en el tiempo.

UNIVERSIDAD INDUSTRIAL DE SANTANDER INGENIERIA DE PETROLEOS

Pág.38

En el método de Crout la matriz A es factorizada como A= L.U en donde la matriz L es una matriz triangular inferior y U una matriz triangular superior con diagonal unitaria. Puede advertirse ya aquí una diferencia con el método de Doolittle, en el cual era la matriz L quien poseía diagonal unitaria. El método de Crout es un procedimiento del tipo recursivo, esto significa el desarrollo de un conjunto de pasos sucesivos en donde el trabajo a realizar en cada paso resulta similar o del mismo tipo pero basado en resultados obtenidos en pasos anteriores. Estas “tareas” a realizar en cada paso o “estación” consisten en la descomposición sucesiva de los menores principales de la matriz de coeficientes A. Se denomina menor principal de A de orden m, denotado como A[m], a una sub-matriz de A, de dimensión [mxm], y en donde sus coeficientes son los ai,j , i=1..m; j=1..m. Es decir: A[1] = a1,1  a1,1 A[ 2] =   a 2,1  a1,1  A[ 3] =  a 2,1  a3,1 

a1, 2  a 2,2 

a1, 2 a 2,2 a3,2

a1,3   a 2,3  a3,3 

etc....

De esta forma, puede escribirse:

en donde:

   

 A[ k ] A[ k +1] =  T  f  [ k +1]

c[ k +1]   a k +1.k +1  

 a1, k +1   a k +1,1  a  a   A2[ ,kk++11] = L[ k +1] .U [ k +1]  k +1,2  c[ k +1] =  :  ; f [ k +1] =  :      :  :     U A[ k ] c[ k +1]   L[ k ] 0 u 1]  a k ,k =+1   [ k ]  a k +[1k, + k     f[Tk +1] a k +1.k +1   l [Tk +1] l k +1.k +1   0 T 1    

a partir de lo indicado, se establece la descomposición de A[k+1] como: en donde:  u1, k +1   l k +1,1  u  l   2, k +1   k +1,2  UNIVERSIDAD INDUSTRIAL DE SANTANDER u[ k +1] =  :  ; l[ k +1] =  :  INGENIERIA DE PETROLEOS      :   :  u k , k +1  l k +1, k     

;

0  0    0 = : Pág.39   : 0 

(1)

Finalmente, de (1): L[ k ] .u [ k +1] = c[ k +1] U [Tk ] .l [ k +1] = f[ k +1]

(2) l k +1, k +1 = a k +1, k +1 − l [Tk +1] .u [ k +1]

Es importante notar que en las ecuaciones anteriores tanto u[k+1] como l[k+1] se obtienen mediante sustitución hacia adelante ya que tanto L[k] como U[k]T son matrices triangulares inferiores. Ejemplo (1): Supóngase el sistema de ecuaciones planteado como ejemplo en el apartado donde se presentó el método de Doolittle. El mismo venía expresado por:

A.x = b →

3 1

-1 7

-2 -3

x1 x2

3

-2

10

x3

+7 -19

=

+4

Se formularán a continuación cada uno de los pasos que conllevan a la factorización de la matriz de coeficientes A en el producto L.U utilizando el método de Crout. •

Paso 1: k=1

L[1] = A[1] = a11 = 3 U[1] = u11 = 1

0   L[1] U [1] L[ 2] =  T  ; U [ 2] =   l[ 2] l 2,2   0  A[1] c[ 2]  2 − 1 A[ 2] =  T =    f [ 2] a 2,2  1 7 

L[1] .u [ 2] = c[ 2]

→

U [T1] .l[ 2] = f [ 2]

→

u [ 2]  1 

3.u[ 2] = −1 1.l[ 2] = 1

∴ ∴

u[ 2] = − 1

3

l[ 2 ] = 1

l 2, 2 = a 2, 2 − l[T2] .u[ 2] = 7 −1.( − 1 ) = +7.33 3

3 L[ 2] =  1

0  1 ; U [ 2] =  7.33   0

− 0.33  1  

UNIVERSIDAD INDUSTRIAL DE SANTANDER INGENIERIA DE PETROLEOS

Pág.40

•

Paso 2: k=2

Como puede verificarse fácilmente L.U. = A, excepto pequeños errores de redondeo. 1.2.4.2.1. Formalización del Algoritmo de descomposición de Crout 0  − 2 .u[ 3] =    7.33  − 3

L[ 2] .u[ 3] = c[ 3]

→

3 1 

U [T2 ] .l[3] = f [ 3]

→

 1 − 0.33 

l3,3 = a3,3 − l[T3] .u[ 3] =10 − [3

− 0.67  u[ 3] =   − 0.32 

∴

0 3  l [ 3] =    1 − 2

∴

3  l [ 3] =   −1

− 0.67  −1].  = +11 .69 − 0.32 

Las expresiones (2) constituyen la forma básica o el núcleo del algoritmo de Crout. Asignando sucesivamente valores a k=1,2...n-1 se obtiene la factorización de la matriz de coeficientes A en el producto L.U. En el esquema siguiente se detalla, en metalenguaje, la estructura del algoritmo de L[ ] 2 L[3] =  T l  [3]

0 3 0  1 7.33 =  l3,3   l T   [3] 

0  1 − 0.33  0 0  1 ; U [3] =    l3,3 0 0    

 u[3]  A[ 2]  ; A[3] =  T   f3 1   

3 c[3]  1  = a3,3   3  

descomposición. Descomposición LU – Crout

l11=a11 , u11=1 Repetir k=1,... n-1 u1,k+1=a1,k+1 / l11 Repetir i= 2,... k s=0

3 L[3] ≡ L =  1  3

0 7.33 −1

 0   11 .69   0

;

1 U [3] ≡ U =  0  0

− 0.33 1 0

UNIVERSIDAD INDUSTRIAL DE SANTANDER INGENIERIA DE PETROLEOS

Pág.41

− 0.67  − 0.32   1  

−1 7 −2

− 2 − 3  10   

Repetir j=1,..., i-1 s = s + lij . uj,k+1 Fin j ui,k+1=(ai,k+1 – s) / lii Fin i lk+1,1 = ak+1,1 Repetir i= 2,... k s=0 Repetir j=1,..., i-1 s = s + uji . lk+1,j Fin j lk+1,i= ak+1,i – s Fin i uk+1,k+1 = 1 s=0 Repetir j=1,..., k s = s + lk+1,j . uj,k+1 Fin j lk+1,k+1= ak+1,k+1 – s Fin k El esquema anterior reproduce el procedimiento de factorización, al mismo es necesario agregarle los procedimientos de sustitución hacia atrás y hacia adelante a fin de obtener la solución del sistema de ecuaciones. No obstante, su desarrollo es idéntico al especificado para el método de Doolittle por lo cual no se repetirá innecesariamente. En lo detallado anteriormente se ha supuesto que la matriz A es tal que permite su descomposición sin pivotamiento; si esto no fuera así alguno de los coeficientes l k+1,k+1 se anularía lo cual impediría la descomposición del siguiente menor puesto que se produciría una división por cero. Las técnicas de pivotamiento también pueden emplearse en estos métodos y son muy extendidas a pesar de que complican considerablemente el algoritmo expuesto. También es importante observar que el número de operaciones necesarias para la descomposición y posterior resolución de los sistemas triangulares coincide exactamente con aquellas previstas por el método de Gauss (con las salvedades ya comentadas en caso de términos independientes múltiples disponibles en forma no simultánea). No obstante, en la forma en la cual se desarrolla el procedimiento de Crout permite que las sumatorias que se incluyen sean resueltas almacenando sus resultados parciales en variables de doble precisión (“s” en este caso), reduciendo así los errores de almacenamiento sin incrementar el volumen de memoria utilizada.

UNIVERSIDAD INDUSTRIAL DE SANTANDER INGENIERIA DE PETROLEOS

Pág.42

1.2.4.2.2. Optimización del algoritmo de Crout En la secuencia representada en la figura siguiente se esquematiza el conjunto de valores involucrados en el cálculo de un elemento ui,k+1 y lk+1,i.

L[K]

U[K]

i→

A[K]

←ui,k+1

x

←ai,k+1

=

L

↑

↑

u[k+1]

c[k+1]

U

A

i↓ L[K]

U[K] x

A[K] =

← l[k+1]

←f[k+1]

↑ lk+1,i

↑ ak+1,i

L U A Los elementos a calcular, sombreados en color negro, son obtenidos a partir de los aquellos valores ubicados en las zonas coloreadas en gris. De ello se deduce que para evaluar un elemento ui,k+1 de U, únicamente es necesario conocer el elemento correspondiente ai,k+1 de A, los elementos de U que se encuentran encima de él en su propia columna (us,k+1 s=1...i-1) y la descomposición previa del menor de A de orden k. De esto último surge que no resultan necesarios ninguno de los coeficientes de A ubicados en la columna k+1 y por encima de a i,k+1 ni del menor de A de orden k (solo es necesaria su descomposición). En resumen, el elemento ai,k+1 es necesario para calcular ui,k+1 pero no para ninguno de los siguientes elementos de U, por consiguiente, se puede almacenar el valor de ui,k+1 recientemente calculado en el mismo lugar donde se encontraba almacenado a i,k+1, es decir, los términos de U pueden almacenarse “sobre” A, es decir, sobre los “casilleros” ubicados por sobre la diagonal de A. Como los términos diagonales de U son iguales a 1 no es necesario almacenarlos. Una conclusión similar puede alcanzarse para los elementos que componen la matriz L, la UNIVERSIDAD INDUSTRIAL DE SANTANDER INGENIERIA DE PETROLEOS

Pág.43

cual puede almacenarse directamente sobre la diagonal y elementos ubicados por debajo de ella en la matriz A. En síntesis, el método de Crout está diseñado para reducir las necesidades de memoria almacenando las dos matrices triangulares en las mismas posiciones de memoria que ocupaba inicialmente la matriz de coeficientes A. Ejemplo (2): En base a la optimización comentada se rehará el proceso de descomposición sobre la matriz de coeficientes del sistema dado en el Ejemplo anterior, pero utilizando el mismo espacio de almacenamiento de A para guardar los elementos de las matrices L y U.

A.x = b →

•

Paso 1: k=1

L[1] = A[1] = a11 = 3 U[1] = u11 = 1

3 1 3

-1 7 -2

L[1] .u [ 2] = c[ 2]

→

U [T1] .l[ 2] = f [ 2]

→

-2 -3 10

x1 x2 x3

3.u[ 2] = −1 1.l[ 2] = 1

=

∴ ∴

l 2, 2 = a 2, 2 − l[T2] .u[ 2] = 7 −1.( − 1 ) = +7.33 3

UNIVERSIDAD INDUSTRIAL DE SANTANDER INGENIERIA DE PETROLEOS

Pág.44

+7 -19 +4 u[ 2] = − 1

3

l[ 2 ] = 1

Los valores obtenidos para u[2], l[2] y l2,2 pueden almacenarse directamente sobre la matriz A original, quedando la misma como:

 =

3 1 3

-1/3 -2 7.33 -3 -2 10

en donde se ha querido representar como  a las componentes del arreglo en donde se han alojado las componentes de L[2] y U[2].

•

Paso 2: k=2 L[ 2] .u [3] = c[ 3]

→

U [T2] .l[3] = f [3]

→

0  − 2 .u [3] =    7.33  − 3 0  1 3  − 0.33 1 l[3] = − 2    

3 1 

l3,3 = a3,3 −l[T3] .u[3] =10 −[3

∴

− 0.67  u [3] =   − 0.32  3  ∴ l[3] =   −1

− 0.67  −1].  = +11 .69 − 0.32 

(a) (b)

(c)

(a) Nótese que en el cálculo de u[3], todos los coeficientes no nulos necesarios para su determinación -tanto de L[2] como de c[3]- se encuentran disponibles sobre Â. Los coeficientes de u[3] calculados se almacenan sobre la matriz  en similar posición a la que ocuparían en la matriz U:

 =

3

-1/3 -0.67

1

7.33 -0.32

3

-2

10

(b) Para el cálculo de l[3] los coeficientes fuera de la diagonal no nulos de U (en este caso –1/3) se encuentra disponible sobre Â. El coeficiente que interviene en los cálculos corresponde a u1,2 el cual se encuentra almacenado en Â1,2. Los coeficientes de l[3] así determinados se almacenan sobre  en idéntica posición a la que ocuparían en la matriz L:  =

3

-1/3 -0.67

1

7.33 -0.32

3

-1

10

(c) Finalmente, el cálculo de l3,3 se realiza en base al coeficiente a3,3 (=10) y los valores de UNIVERSIDAD INDUSTRIAL DE SANTANDER INGENIERIA DE PETROLEOS

Pág.45

l[3]T y u[3] determinados en los pasos anteriores. Todos estos coeficientes se encuentran disponibles sobre la matriz Â. El valor de l3,3 así calculado es almacenado sobre  en la misma posición en la que se ubicaría en la matriz L:

 =

3

-1/3 -0.67

1

7.33 -0.32

3

-1 11.69

Como puede observarse, sobre el triángulo inferior de  se ubican los coeficientes de la matriz L, mientras que los elementos ubicados por sobre la diagonal resultan los correspondientes a la matriz U, excepto los términos diagonales, los cuales son conocidos a priori e iguales a la unidad.

1.2.4.3. Método de Cholesky Se considera aquí un procedimiento de factorización desarrollado para aquellos sistemas de ecuaciones en donde la matriz de coeficientes es simétrica y definida positiva. Este tipo de matrices resultan muy usuales en la solución de numerosos problemas de ingeniería mediante técnicas numéricas como elementos finitos, y de allí la necesidad de incluirlas de modo particular dentro del conjunto de procedimientos de resolución considerados en este capítulo. El método de Cholesky utiliza el preconocimiento de las características mencionadas para realizas una descomposición más eficaz que aquella llevada adelante por el método de Crout. En líneas generales, la secuencia paso a paso considerada por Cholesky es similar a la contemplada por Crout con la significativa diferencia que la matriz A es factorizada como A=L.LT , es decir, ahora solamente los coeficientes de L son incógnitas ya que su traspuesta hace las veces de matriz U en el proceso de descomposición.

El menor de orden 1 se obtiene nuevamente de una forma muy simple: A[1] = a11 = L[1].L[1]T = L112 ∴ L11 = a11½ La expresión de recurrencia puede obtenerse suponiendo ya desarrollada la T UNIVERSIDAD INDUSTRIAL DE A SANTANDER [ k +1] = L[ k +1] .L

[ k +1]

INGENIERIA DE PETROLEOS

 A[ k ]   fT  [ k +1]

f[ k +1]   = a k +1.k +1    

L[ k ] l [Tk +1]

  LT  [ k ] l k +1.k +1   0 T  0

Pág.46

l [ k +1]   l k +1.k +1  

(3)

descomposición de orden k, A[k]= L[k].L[k]T, de la cual puede obtenerse la descomposición para el menor de orden k+1: Nótese que, por ser A una matriz simétrica, los vectores c[k+1] y f[k+1] definidos en las ecuaciones (1) coinciden. De (3) se obtienen las ecuaciones necesarias para la descomposición del menor principal A[k+1]: L[ k ] .l [ k +1] = f[ k +1] l k +1, k +1 = ak +1, k +1 − l [Tk +1] .l [ k +1]

Debe notarse que el cálculo de l[k+1] indicado en la primera de las ecuaciones anteriores implica una sustitución hacia adelante, considerando las características de matriz triangular inferior de L. Ejemplo Factorizar la matriz de coeficientes del siguiente sistema de ecuaciones: 2

-1

0

x1

-1

2

-1

x2

0

-1

2

x3

5 =

0 0

La matriz A es simétrica y definida positiva. a) Verificación del carácter de definida positiva de A. A[1] = 2 → |A[1]| >0 A[2] =

A[3] =

→ |A[2]|=3 >0

2

-1

-1

2

2

-1

0

-1

2

-1

0

-1

2

→ |A[3]|=4 >0

c) Descomposición UNIVERSIDAD INDUSTRIAL DE SANTANDER INGENIERIA DE PETROLEOS

Pág.47

•

k=1

L[1] =

a11 =

L[1] .l[ 2] = −1 ∴ l[ 2] = l2,2 =

•

 A[1] A[ 2] =  T f  [ 2]

2 = 1.414

f[ 2]  + 2 =  a2, 2   −1 

−1  + 2 

−1 = −0.707 1.414

a2,2 − l[T2] .l[ 2] =

2 − ( −0.707 ) 2 = 1.225

 1.414 L[ 2 ] =  − 0.707

0  1.225 

k=2

L[ 2] .l[3] = f[3]

l3,3 =

 1.414 ∴  −0.707

a3,3 −l[T3] .l[3] =

 A[ 2] A[ 3] =  T f  [ 3]

 0  .l[3] =    1.225  −1 0

2 −[0

 2 −1 f[3]  −1 2 = a3,3   0 −1   

 0  → l[3] =   − 0 . 816  

 0  − 0.816 ].  =1.155 − 0 . 816  

0 −1 +2

     

 1.414 L[3] =  − 0.707   0

UNIVERSIDAD INDUSTRIAL DE SANTANDER INGENIERIA DE PETROLEOS

Pág.48

0 1.225 − 0.816

 0   1.155   0

1.2.4.3.1. Formalización del Algoritmo de descomposición de Cholesky En el esquema siguiente se detalla, en metalenguaje, la estructura del algoritmo de descomposición de Cholesky para matrices simétricas y definidas positivas. Descomposición LLT – Cholesky

l11= (a11)½ Repetir k=1,... n-1 lk+1,1=ak+1,1 / l11 Repetir i= 2,... k s=0 Repetir j=1,..., i-1 s = s + lij . lk+1,j Fin j lk+1,i=(ak+1,i – s) / lii Fin i s=0 Repetir i= 1,... k s = s + (lk+1, i)2 Fin i lk+1,k+1= (ak+1,k+1 – s) ½ Fin k Es posible demostrar que el algoritmo detallado puede utilizarse para cualquier matriz simétrica y definida positiva.

UNIVERSIDAD INDUSTRIAL DE SANTANDER INGENIERIA DE PETROLEOS

Pág.49

1.3. Técnicas de Almacenamiento de Matrices En cualquier lenguaje de programación de última generación, las matrices se almacenan en arreglos de memoria bidimensionales, en donde cada elemento es almacenado en un espacio direccionado con el nombre asignado a la matriz y el número de fila y columna en donde se ubica: A(3,2), B(1), etc. La sintaxis específica corresponde a la utilizada por cada lenguaje en particular. En el marco de estas notas se empleará la mencionada anteriormente es decir aquella en la cual cada elemento se especifica mediante el nombre del arreglo seguido, entre paréntesis, del número de fila y columna en donde se encuentra. Considerando lo natural y sencillo de lo mencionado no parece, a primera vista, necesario ni útil que deba implementarse alguna otra forma o procedimientos para almacenar las matrices en la memoria de un computador. Sin embargo imagínese por un instante la necesidad de almacenar una matriz simétrica, es evidente que solo resultaría imprescindible asignar espacios de memoria para la matriz triangular superior (o inferior, es indistinto). Supóngase, por ejemplo, que se trata de una matriz cuadrada de orden 100; en este caso al asignar un espacio de memoria para la matriz (paso imprescindible antes de su utilización en cualquier lenguaje) se reservará un arreglo de 100 filas y 100 columnas, del cual solo serán utilizados aquellos correspondientes a la diagonal y a los elementos ubicados encima –o debajo- de ella. Evidentemente, tal situación implica un malgasto de memoria que, en caso de matrices de varios miles de ecuaciones, puede limitar el tamaño del problema a resolver al no disponer el computador de una cantidad de memoria suficiente. Un caso similar, aunque más extremo aún, es el caso de las matrices bandeadas en donde los elementos no nulos se limitan a un porcentaje a menudo muy inferior que n 2 (n= orden de la matriz) razón por la cual resulta sumamente importante desarrollar metodologías de almacenamiento más eficientes. Sobre este punto en particular existe un considerable desarrollo dentro de las técnicas de cálculo numérico orientados a determinados tipos especiales de matrices, característicos de algunos problemas de la ingeniería. En los puntos siguientes se tratarán algunas de estas técnicas las cuales, claro está, son solo un pequeño segmento del numeroso campo existente sobre este particular.

UNIVERSIDAD INDUSTRIAL DE SANTANDER INGENIERIA DE PETROLEOS

Pág.50

1.3.1. Almacenamiento por columnas Consiste en almacenar en forma consecutiva en un vector de longitud mxn cada una de las columnas que conforman la matriz cuya dimensión, en general, se ha supuesto de m filas y n columnas. De esta forma, los coeficientes de la matriz serán ubicados dentro del vector en base al siguiente ordenamiento: CT = {c1, c2, ...cm, cm+1,cm+2.... cm+m, c2m+1,c2m+2...c3m, ......, c(n-1)m+1, c(n-1)m+2 ... cn.m} = = {a1,1, a2,1,...am,1, a1,2, a2,2 ,...am,2 , ......, a1,n, a2,n,... am,n} En líneas generales, en el elemento cK se almacena el término a(i,j), en donde: k = (j-i) * m + i la cual, es la ecuación de “transformación” del “sistema de posicionamiento global” (i, j), al “local”, k. 1.3.2. Almacenamiento por filas Una alternativa al almacenamiento por columnas presentado anteriormente es disponer los elementos ubicados sobre cada una de las filas de la matriz en forma sucesiva sobre un arreglo unidimensional o vector, esto es: f = {f1, f2, ..... fn, fn+1, fn+2...., f2n, f2n+1, f2n+2...., f3n..........f(m-1)*n+1, f(m-1)*n+2,..... fm*n} = = {a1,1, a1,2,... a1,n, a2,1, a2,2, ... a2,n, a3,1, a3,2,....,a3,n ........am,1, am,2,.... am,n} Una vez almacenada una matriz sobre un vector en base al ordenamiento comentado, es necesario establecer las relaciones de transformación que permitan ubicar a un elemento cualquiera a(i,j) de la matriz dentro del vector f. Es fácil comprobar que esa relación viene dada por la expresión: k = (i-1) * n + j en la cual “i” y “j” corresponden a la fila y columna respectivamente del elemento de la matriz a ubicar, y “k” la posición sobre el vector en donde los mismos se encuentran almacenados. De esta forma, en cualquiera de los algoritmos de resolución de ecuaciones UNIVERSIDAD INDUSTRIAL DE SANTANDER INGENIERIA DE PETROLEOS

Pág.51

vistos anteriormente, podrá reemplazarse cada operación sobre a(i,j) por correspondiente f(k) en forma enteramente equivalente.

el

1.3.3. Almacenamiento de matrices simétricas o triangulares. En una matriz simétrica resulta evidente que solo resulta imprescindible almacenar aquellos elementos que se ubican en la diagonal y por sobre –triangular superior- o por debajo –triangular inferior- de la misma. Una situación similar presentan las matrices triangulares (superiores o inferiores) en donde los coeficientes por debajo o por encima de la diagonal resultan nulos y por lo tanto no resulta necesario almacenarlos en la memoria. A fin de retener solamente aquellos coeficientes necesarios para el cálculo, una alternativa es disponer en un arreglo unidimensional -vector- cada una de las columnas de la matriz pero limitando su extensión a los elementos ubicados entre la primer fila y la diagonal. Se estará almacenando así una matriz triangular superior. En caso de que se requiera almacenar una matriz triangular inferior –o la mitad inferior de una matriz simétrica- las transformaciones son similares y pueden encararse sin dificultad por lo cual no serán comentadas aquí en forma explícita. A modo de ejemplo, puede suponerse el caso de una matriz triangular superior –el caso de una matriz simétrica es enteramente equivalente-: i↓ a11 0  i→ 0 A=   j→ 0  0

a12 a22       0

  



a 1i a 2i  aii     0

j↓  

 

 





a1j a2j  aij   ajj 0

 

 

a 1n a 2n  ain   ajn  ann

Los coeficientes no nulos de la misma se almacenan siguiendo un ordenamiento por columnas en un vector de acuerdo al siguiente esquema: CT = {a11, a12, a22, a13, a23, a33, ......, a1n, a2n, .... ann}

UNIVERSIDAD INDUSTRIAL DE SANTANDER INGENIERIA DE PETROLEOS

Pág.52

La relación de vinculación entre la posición de cada coeficiente sobre la matriz A y su correspondiente en el vector C es establecida mediante la siguiente relación de transformación: k =

j.( j −1) +i 2

en la cual “i” y “j” representan la fila y columna respectivamente de un elemento sobre la matriz A y “k” la posición donde éste se encuentra almacenado en el vector C. Ejemplo : utilizar la técnica de almacenamiento mencionada anteriormente dentro del algoritmo de solución de sistemas de ecuaciones con matrices simétricas de acuerdo al procedimiento de eliminación de Gauss detallado en el punto 1.2.3.4. Los coeficientes de la matriz se encuentran almacenados en un vector C con las características ya comentadas. El algoritmo especificado en aquel punto, adaptado a las condiciones de almacenamiento impuestas por la forma de la matriz, se detalla en el siguiente esquema. Gauss – Simetrico Almacenamiento por columnas de la matriz triangular superior. Definir: f[i,j] = j*(j-1)/2 + i Repetir k=1,... n-1 Repetir i=k+1,... n m = C(f[k,i]) / C(f[k,k]) Repetir j=i,... n C(f[i,j]) = C(f[i,j]) – m. C(f[k,j]) Fin j bi = bi – m.bk Fin i Fin k xn = bn/ann Repetir i=n-1,.... 1 s=0 Repetir j=i+1,... n s = s + C(f[i,j]).xj Fin j xi = (bi-s)/C(f[i,i]) Fin i UNIVERSIDAD INDUSTRIAL DE SANTANDER INGENIERIA DE PETROLEOS

Pág.53

en donde se ha hecho uso de una función predefinida (cuya sintaxis dependerá del lenguaje utilizado) la cual permite obtener la posición del coeficiente a(i,j) genérico dentro del vector C donde se encuentra almacenado. 1.3.4. Almacenamiento de matrices “bandeadas” Como se comentara en el punto 1.2.3.5., es muy frecuente que las matrices de coeficientes que forman parte de sistemas de ecuaciones lineales provenientes de la solución numérica de ecuaciones diferenciales que gobiernas problemas de ingeniería posean una forma “bandeada”. En aquel apartado se especificó un algoritmo de resolución basado en una adecuación del método de eliminación de Gauss para el caso que la matriz resulte, además de bandeada, simétrica. Ese algoritmo simplemente acotaba el “espacio” dentro del cual se realizaban las operaciones de eliminación, limitándolas a los coeficientes ubicados dentro de la banda. Es inmediato que si esto es así, no es necesario almacenar los coeficientes que se encuentran fuera de ella, punto este que será estudiado a continuación. A fin de tratar el tema con más generalidad respecto de lo indicado en el punto 1.2.3.5. se considerarán matrices “bandeadas” no necesariamente simétricas y en las cuales la “semibanda” por debajo de la diagonal, no necesariamente tiene la misma dimensión que aquella ubicada por sobre la misma. En resumen, una matriz de esta características presentaría el siguiente aspecto: ← ↑ g ↓ A=

a11 a12 a21 a22  ag+1,1 

u

→

 a1,u+1  

 

 





 

 an,n-g

an-1,n-1  an,n-1

←

g



→

en la cual: • g= semiancho inferior → a(i,j)=0 si i>j+g • u= semiancho superior → a(i,j)=0 si j>i+u UNIVERSIDAD INDUSTRIAL DE SANTANDER INGENIERIA DE PETROLEOS

Pág.54

an-u,n  an-1,n ann

↑ u ↓

Sabiendo que los coeficientes no nulos se ubican dentro de una banda de dimensión conocida, los coeficientes de la matriz A se almacenaran dentro de una matriz de dimensión ajustada al ancho de banda, es decir, una matriz rectangular de n filas y g+u+1 columnas cuya “organización” obedece al siguiente esquema:

a21  R=

a11 a22 

a12 a23



ag+1,1 ag+1,g+1        an-1,n-1 an-1,n an,n-g  an,n-1 ann ←

g

→

←

u

a1,u+1 a1,u+2  ag+1,u+g+1  an-u,n

→

Como se observa en el esquema anterior, la matriz R posee el mismo número de filas que A y un número de columnas igual al ancho de banda, esto es 1+u+g. En el caso particular (y frecuente) que u=g, los elementos diagonales de la matriz de coeficientes se almacenarán sobre la columna central de R. Evidentemente, este esquema de almacenamiento presentará ventajas sobre aquel que maneja directamente la matriz como llena, en el caso que n 2>>n(1+u+g) situación que es frecuente en los sistemas reales que responden a la solución numérica de ecuaciones diferenciales de problemas de ingeniería en donde n adquiere dimensiones de miles o superior, y el ancho de banda total (1+u+g) no supera el 10 o 20% de n. Establecida la forma en que se realizará el almacenamiento quedan por establecer las relaciones de transformación que permitan, dado un elemento genérico a(i,j) de la matriz de coeficientes del sistema, ubicarlo dentro de la tabla de almacenamiento diseñada. Es sencillo verificar que tales relaciones de transformación quedan establecidas por las siguientes expresiones: α = i β = 1+g+(j-i) las cuales permiten obtener el coeficiente R(α, β ) en donde se aloja el término A(i,j) de la matriz original. Una instancia posterior es disponer que la matriz R se almacene por filas o columnas dentro de un vector, disposición muy frecuente cuando se trabaja en lenguajes como Fortran. Para ello, no hay más que aplicar las reglas de transformación detalladas al analizar estos tipos de almacenamiento en los puntos anteriores. Por ejemplo, si se UNIVERSIDAD INDUSTRIAL DE SANTANDER INGENIERIA DE PETROLEOS

Pág.55

adopta la disposición de almacenar a la matriz R por columnas dentro de un vector V, la regla de transformación quedará finalmente como: k = (β -1).n+α = [g+(j-i)].n +i la cual permite obtener la posición “k” del vector V en donde se aloja el coeficiente A(i,j) de la matriz original. Se hace notar que al almacenar la matriz R por columnas, lo que se está realizando en realidad es un almacenamiento de A por diagonales, ya que cada columna de R se corresponde con una diagonal de la matriz original A. Esta disposición es particularmente útil cuando los coeficientes de la matriz original son generados por un programa que incluye también en su desarrollo la solución del sistema de ecuaciones. Ejemplo: Supóngase que se desea encontrar la solución del sistema de ecuaciones A.x = b donde A es una matriz de orden n, simétrica, en banda y cuyos coeficientes resultan:

↑ 2 ↓

5 -4 1

A=

← -4 6 -4 1

2→ 1 -4 6 -4 

1 -4 6 1

1 -4 -4 1

1 6 -4

 -4 5

la cual, de acuerdo a lo comentado anteriormente, puede almacenarse en una matriz rectangular con el siguiente formato:

R=

5 6 6  6 6 5

← -4 -4 -4  -4 -4

2→ 1 1 1 1

en la cual se ha prescindido de almacenar la mitad inferior de la matriz de coeficientes a partir de las características de simetría de la misma.

UNIVERSIDAD INDUSTRIAL DE SANTANDER INGENIERIA DE PETROLEOS

Pág.56

La regla de transformación que permite ubicar a un elemento genérico de la matriz A dentro de R, es según lo indicado oportunamente: α = i β = g+1+(j-i) –g = 1+j-i donde α y β representan la fila y la columna en donde el elemento a(i,j) se encuentra almacenado en R . En la expresión anterior se ha descontado el valor de “g” de la ecuación original en función de haberse almacenado únicamente la una de las “semibandas”. Por ejemplo, puede comprobarse que el coeficiente a(3,4) = –4, se encuentra almacenado en las posición R(3,2) según resulta de: α = 3 β = 1 +4 – 3 = 2 Finalmente, se reformula el algoritmo Gauss Banda Simétrico ya planteado oportunamente, considerando que solamente se almacena la mitad de la matriz en la forma compacta detallada más arriba. Definir F(i,j) = 1+j-i Repetir k=1,... n-1 fin=min(n,k+g) Repetir i=k+1,... fin m = R[k,F(k,i)]/R[k,F(k,k)] Repetir j=i,... fin R[i,F(i,j)] = R[i,F(i,j)] – m. R[k,F(k,j)] Fin j bi = bi – m.bk Fin i Fin k xn = bn/R[n,F(n,n) Repetir i=n-1,.... 1 fin =min(n,i+g) s=0 Repetir j=i+1,... fin s = s + R[i,F(i,j)].xj Fin j xi = (bi-s)/R[i,F(i,i)] Fin i

En el esquema anterior se ha denominado como “g” al ancho de la semibanda (2 para este ejemplo). UNIVERSIDAD INDUSTRIAL DE SANTANDER INGENIERIA DE PETROLEOS

Pág.57

1.4 Métodos Iterativos 1.4.1 Introducción Los métodos de solución de sistemas de ecuaciones pertenecen esencialmente a la clase de métodos directos o a la clase de métodos iterativos. Los métodos directos resuelven los sistemas de ecuaciones en un numero determinado de operaciones aritméticas (sumas, restas, multiplicaciones y divisiones), y los errores en la solución surgen enteramente de los errores de redondeo introducidos al efectuarse los cálculos. Básicamente, los métodos directos son los métodos de eliminación, de los cuales el

mas conocido es el método sistemático de eliminación gaussiano, y los

métodos de descomposición triangular. Estos últimos factorizan la matriz de coeficientes A de las ecuaciones A . x = b en A = L . U, donde L y U son matrices triangulares inferiores y superiores respectivamente; una vez que la descomposición ha sido determinada la solución es calculada resolviendo el sistema equivalente L . U . x = b, haciendo U. x = y; entonces se obtiene y resolviendo L . y = b por sustitución hacia adelante, y por último se obtiene x resolviendo U . x = y por sustitución hacia atrás. Con ambos métodos es usualmente necesario emplear pivoteo parcial para controlar el crecimiento de los errores de redondeo. Es de interés notar que cuando los cálculos relacionados con estos métodos son realizados sobre grandes computadoras modernas, los errores de redondeo introducidos durante el proceso de calculo frecuentemente tienen menos efecto sobre la solución que los errores de redondeo en los coeficientes y constantes de las ecuaciones. Si suponemos que los coeficientes y constantes de un sistema de n ecuaciones algebraicas lineales son todos menores que uno en modulo y que son conocidos con certeza en solo q lugares decimales, entonces ellos poseen errores relativos menores o iguales a ½ 10 -q. Se ha probado que si los cálculos del método de eliminación Gaussiano se llevan adelante con números cuya cantidad de lugares decimales sea mayor que (q + log 10 n) , es decir log10 n lugares decimales mas que en los coeficientes y constantes, la solución obtenida para q lugares decimales es la solución exacta de un conjunto de n ecuaciones lineales cuyos coeficientes y constantes difieren en menos de ½ 10–q de aquellos correspondientes a las ecuaciones originales (si se usa representación binaria, la mas usada en las computadoras, 10-q es reemplazado por 2-q y log10 n por log2 n). Dicho de otra forma, la solución es tan precisa como los datos lo permitan. Esto no significa que el UNIVERSIDAD INDUSTRIAL DE SANTANDER INGENIERIA DE PETROLEOS

Pág.58

numero de cifras significativas correctas en la solución sea el mismo que en el de los datos. Aquel podría ser menor, como es en el caso en que la perdida de precisión no ocurre por el método de solución utilizado sino debido a que las ecuaciones estén mal condicionadas en el sentido que pequeños cambios en los coeficientes producen grandes cambios en la solución. El numero de dígitos significativos en la solución no puede exceder el numero de cifras significativas en los primeros q lugares decimales del pívot mas pequeño. Aunque la eliminación Gaussiana y la descomposición LU son matemáticamente equivalentes, el error total del método LU puede ser hecho menor que el del método de eliminación Gaussiana si los productos internos son acumulados en variables de doble longitud, una operación que puede realizarse en algunas computadoras. Sin la acumulación de los productos internos en variables de doble longitud los dos métodos son equivalentes tanto matemática como numéricamente. Si la aproximación de la solución obtenida necesita ser mejorada el siguiente método de refinamiento es usualmente recomendado. Llamemos a la solución calculada de A . x = b como x(1). Luego, calculemos el residuo r(1) = b – A . x(1) usando variables de doble longitud, es decir de 2 . t dígitos decimales o binarios en lugar de los t dígitos estándar de las variables de simple longitud. Si hacemos x = x(1) + d, donde d es el termino de corrección exacto que se suma a la solución calculada en primera instancia, x(1), para obtener la solución exacta x, entonces A . x = b puede expresarse como A . (x(1) + d) = b. Distribuyendo el producto tendremos A . x(1) + A . d = b, y finalmente A . d = b - x(1) = r(1). La solución de esta ultima ecuación no es d exactamente sino una solución aproximada calculada que llamaremos d(1). Luego, una solución mejorada será x(2) = x(1) + d(1). Este procedimiento puede repetirse hasta lograr el grado de aproximación deseada. Con el objeto de mantener la solución de las ecuaciones en una aritmética de simple longitud, es decir t dígitos, los residuos podrían ser multiplicados por 10t, o por 2t, y los correspondientes términos de corrección d(1), d(2)......, multiplicados por 10-t, o por 2-t. Es evidente que este método de refinamiento iterativo requiere solamente una pequeña cantidad adicional de operaciones debido a que la solución de A . d = r(1) involucra solo sustituciones hacia delante y hacia atrás porque la descomposición A = L . U ha sido calculada con anterioridad. Si la matriz de coeficientes puede ser entrada totalmente en la memoria ( el almacenamiento de rápido de la computadora), entonces los métodos directos, en general, son más rápidos y aproximados que los métodos iterativos. Además, si la matriz UNIVERSIDAD INDUSTRIAL DE SANTANDER INGENIERIA DE PETROLEOS

Pág.59

de coeficientes tiene alguna propiedad o estructura en especial, en particular si son simétricas, definidas positivas o bandeadas, es posible incrementar el numero de ecuaciones a resolver mediante una programación mas eficiente que tenga en cuenta tales particularidades. Los métodos directos son ciertamente preferibles a los métodos iterativos cuando: 1. Se deben resolver varios sistemas de ecuaciones con la misma matriz de coeficientes pero diferentes vectores de términos independientes. 2. La matriz de coeficientes es aproximadamente singular. En este caso pequeños residuos (donde el residuo esta definido por r(1) = b – A . x(1)) no implican pequeños errores en la solución. Esto resulta evidente si consideramos que A-1 . r(1) = A-1 . (b – A . x(1)) = A-1 . b – x(1) = x – x(1). Entonces x – x(1) podría tener grandes componentes

cuando las componentes del vector residual r(1) son pequeñas

porque algunos de los elementos de A-1 serán grandes si A es aproximadamente singular. En la mayoría de los problemas asociados a la solución de ecuaciones diferenciales parciales las matrices son bandas y esparzas, esto es, la cantidad de elementos nulos es mucho mayor que la cantidad de elementos no nulos y están ordenados a lo largo de una estrecha banda que rodea a la diagonal principal. Para este tipo de matrices los métodos de eliminación gaussiana en su versión estándar son ineficientes ya que se efectúan una gran cantidad de operaciones innecesarias. Debe puntualizarse, como ya hemos visto, que es posible adaptar el método de Gauss a matrices bandas y esparzas . En cambio, en los métodos iterativos, combinados con técnicas eficientes de almacenamiento de matrices esparzas, hay muy pocas operaciones aritméticas asociadas con los coeficientes iguales a cero y por lo tanto se debe almacenar en la memoria de la computadora una cantidad de números considerablemente menor. Como consecuencia de esto, los métodos iterativos pueden ser usados para resolver sistemas de ecuaciones que son demasiado grandes para el uso de métodos directos. Además, la programación y el manejo de datos al aplicar métodos iterativos es mucho mas simple que cuando aplicamos métodos directos, y fundamentalmente, los métodos iterativos pueden ser aplicados a la resolución de sistemas de ecuaciones no lineales, cosa que no puede realizarse con los métodos directos. No obstante, el uso eficiente de los métodos iterativos es muy dependiente del calculo directo o estimación del valor de algunos parámetros de aceleración de la convergencia UNIVERSIDAD INDUSTRIAL DE SANTANDER INGENIERIA DE PETROLEOS

Pág.60

(ciertos métodos los poseen y otros no) y de que la matriz de coeficientes sea bien condicionada. De no ser así la convergencia será lenta y el volumen de operaciones aritméticas enorme. Con parámetros de aceleración óptimos el volumen de calculo de los métodos iterativos para grandes sistemas de ecuaciones puede ser mucho menor que en los métodos directos.

1.4.2 Fundamentos de los Métodos iterativos Un método iterativo para resolver ecuaciones es aquel en el cual una primera aproximación es usada para calcular una segunda aproximación, la cual a su vez es usada para calcular una tercera aproximación, y así sucesivamente. Se dice que el procedimiento iterativo es convergente cuando las diferencias entre la solución exacta y las sucesivas iteraciones tienden a cero al incrementarse el numero de iteraciones. En general la solución exacta nunca es obtenida en un numero finito de pasos, pero esto no es importante. Lo que importa es que las sucesivas iteraciones converjan rápidamente a valores que puedan considerarse correctos dentro de una precisión especificada. Como se mencionó previamente podríamos considerar el uso de métodos iterativos cuando, entre otros casos,

los métodos directos requieran mas espacio de

almacenamiento rápido (memoria) que el disponible y la matriz de coeficientes sea esparza pero bien condicionada. Para desarrollar los métodos de resolución en forma sencilla consideremos, por ejemplo, la resolución de un sistema general de cuatro ecuaciones: a11 ⋅ x1 + a12 ⋅ x2 + a13 ⋅ x3 + a14 ⋅ x4 = b1 a21 ⋅ x1 + a22 ⋅ x2 + a23 ⋅ x3 + a24 ⋅ x4 = b2 a31 ⋅ x1 + a32 ⋅ x2 + a33 ⋅ x3 + a34 ⋅ x4 = b3 a41 ⋅ x1 + a42 ⋅ x2 + a43 ⋅ x3 + a44 ⋅ x4 = b4 donde aii ≠ 0 para i = 1 a 4

El sistema anterior puede ser reescrito de la siguiente manera:

UNIVERSIDAD INDUSTRIAL DE SANTANDER INGENIERIA DE PETROLEOS

Pág.61

x1 =

1 ⋅ ( b1 − a12 ⋅ x2 − a13 ⋅ x3 − a14 ⋅ x4 ) a11

x2 =

1 ⋅ ( b2 − a21 ⋅ x1 − a23 ⋅ x3 − a24 ⋅ x4 ) a22

x3 =

1 ⋅ ( b3 − a31 ⋅ x1 − a32 ⋅ x2 − a34 ⋅ x4 ) a33

x4 =

1 ⋅ ( b4 − a41 ⋅ x1 − a42 ⋅ x2 − a43 ⋅ x3 ) a44

A partir de este punto plantearemos dos métodos: Jacobi y Gauss Seidel.

1.4.2.1 Método de Jacobi

Llamamos

) x(1 i

a la primera aproximación de

xi

,

) x(2 i

a la segunda aproximación, y así

sucesivamente, y asumimos que n de ellas han sido calculadas, es decir que conocemos ) x(n i

para i = 1 a 4. Entonces el método de Jacobi expresa los (n+1)-ésimos valores

iterativos exclusivamente en función de los n-ésimos valores iterativos. La iteración correspondiente al sistema de ecuaciones antes descripto es:

( ⋅ (b − a ⋅ (b − a ⋅ (b − a

x1( n + 1 ) =

1 ⋅ b1 − a12 ⋅ x(2n ) − a13 ⋅ x3( n ) − a14 ⋅ x4( n ) a11

x2( n + 1 ) =

1 a22

x3( n + 1 ) =

1 a33

x4( n + 1 ) =

1 a44

)

2

(n) (n) (n) 21 ⋅ x1 1 − a23 ⋅ x3 − a24 ⋅ x4

3

(n) (n) (n) 31 ⋅ x1 − a32 ⋅ x2 − a34 ⋅ x4

4

(n) (n) (n) 41 ⋅ x1 − a42 ⋅ x2 − a43 ⋅ x3

)

) )

En el caso general para m ecuaciones queda

xi( n +1 ) =

1 aii

i −1 m   ⋅ bi − ∑aij ⋅ x(j n ) − ∑aij ⋅ x(j n )    j =1 j =i +1  

con

i =1a m

El método de Jacobi no es incondicionalmente convergente, esto es, no siempre las sucesivas aproximaciones convergen a un valor determinado. La matriz de coeficientes del sistema a resolver debe cumplir con ciertas condiciones para asegurar la UNIVERSIDAD INDUSTRIAL DE SANTANDER INGENIERIA DE PETROLEOS

Pág.62

convergencia del método. Puede demostrarse que el método de Jacobi aplicado al sistema de ecuaciones A . x = b es convergente si la matriz A es estrictamente diagonal dominante, es decir, si en cada fila de A el módulo del elemento diagonal es mayor que la suma de los módulos de los elementos restantes de la fila.

Algoritmo del Método de Jacobi Ingresar  s Repetir i = 1,....n x1i = valores iniciales Fin repetir Repetir mientras Flag = 1 Flag = 0 Repetir i = 1,....n sum = 0 Repetir j = 1,....n Si i j sum = sum + aij . X1j Fin si Fin repetir x2i = ( bi - sum ) / aii Si x2i – x1i  > s Flag = 1 Fin si Fin repetir Repetir i = 1,....n x1i = x2i Fin repetir Fin repetir mientras

1.4.2.2 Método de Gauss Seidel

UNIVERSIDAD INDUSTRIAL DE SANTANDER INGENIERIA DE PETROLEOS

Pág.63

En este método (n+1)-ésimos valores iterativos son usados tan pronto como ellos están disponibles y la iteración correspondiente a las ecuaciones anteriores es definida por:

( ⋅ (b − a ⋅ (b − a ⋅ (b − a

x1( n +1 ) =

1 ⋅ b1 − a12 ⋅ x2( n ) − a13 ⋅ x3( n ) − a14 ⋅ x4( n ) a11

x2( n +1 ) =

1 a22

x3( n +1 ) =

1 a33

x4( n +1 ) =

1 a44

) )

2

( n +1 ) − a23 ⋅ x3( n ) − a24 ⋅ x4( n ) 21 ⋅ x1 1

3

( n +1 ) − a32 ⋅ x2( n +1 ) − a34 ⋅ x4( n ) 31 ⋅ x1

4

( n +1 ) − a42 ⋅ x(2n +1 ) − a43 ⋅ x3( n +1 ) 41 ⋅ x1

) )

En el caso general para m ecuaciones queda

xi( n +1 ) =

1 aii

i −1 m   ⋅ bi − ∑aij ⋅ x(j n +1 ) − ∑aij ⋅ x(j n )    j =1 j =i +1  

con

i =1a m

El método de Gauss – Seidel para ser convergente en la resolución de un sistema de ecuaciones A . x = b requiere que la matriz de coeficientes A cumpla con la misma condición requerida en el método de Jacobi.

Algoritmo del Método de Gauss – Seidel Ingresar 

s

Repetir i = 1,....n UNIVERSIDAD INDUSTRIAL DE SANTANDER INGENIERIA DE PETROLEOS

Pág.64

x1i = valores iniciales Fin repetir Repetir mientras Flag = 1 Flag = 0 Repetir i = 1,....n sum = 0 Repetir j = 1,....n Si i  j sum = sum + aij . X1j Fin si Fin repetir A = ( bi - sum ) / aii Si mod  A - x1i  > es Flag = 1 Fin si x1i = A A=0 Fin repetir Fin repetir mientras

Ejemplo: Sea una función u que satisface la ecuación: ∂u ∂2u = ∂t ∂x 2

para 0 < x < 1

y

t > 0.

Las condiciones iniciales son: u=1

en 0 < x < 1

para

t = 0.

x=1

para

Y las condiciones de contorno son: u=0

en x = 0

y

UNIVERSIDAD INDUSTRIAL DE SANTANDER INGENIERIA DE PETROLEOS

Pág.65

t  0.

Resolviendo esta ecuación por algún método que se vera mas adelante se llega a un sistema de ecuaciones algebraicas lineales. Las incógnitas de dicho sistema de ecuaciones son valores discretos de la función en determinados puntos del dominio. En este caso se tomaron 9 puntos interiores equidistantes en el intervalo 0 < x < 1. Entonces, en principio, las incógnitas del sistema de ecuaciones serían 9. Por cuestiones de simetría de la función solución respecto del punto central del intervalo en el sistema de ecuaciones aparecerán solo 5 incógnitas que denotaremos por ui con i = 1 a 5. Planteada la situación de esta manera el sistema resultante es: 4 ⋅ u1 − u2 = 1 − u1 + 4 ⋅ u2 − u3 = 2 − u2 + 4 ⋅ u3 − u4 = 2 − u3 + 4 ⋅ u4 − u5 = 2 − u4 + 2 ⋅ u5 = 1

Este sistema también puede expresarse como: 1 ⋅ ( u2 + 1) 4 1 = ⋅ ( u1 + u3 + 2 ) 4 1 = ⋅ ( u 2 + u4 + 2 ) 4 1 = ⋅ ( u3 + u5 + 2 ) 4 1 = ⋅ ( u4 + 1) 2

u1 = u2 u3 u4 u2

Tomaremos como valores de la primera aproximación: ui = 1

para

i = 1 a 5.

a) Resolución del sistema mediante el método de Jacobi Las ecuaciones de iteración son:

UNIVERSIDAD INDUSTRIAL DE SANTANDER INGENIERIA DE PETROLEOS

Pág.66

u1( n +1 ) = u(2n +1 ) = u3( n +1 ) = u4( n +1 ) = u5( n +1 ) =

1 ⋅ 4 1 ⋅ 4 1 ⋅ 4 1 ⋅ 4 1 ⋅ 2

(u (u (u (u (u

)

(n) 2 +1

(n) (n) 1 + u3 (n) (n) 2 + u4

(n) (n) 3 + u5

)

) + 2) + 2) +2

(n) 4 +1

Entonces, partiendo de la primera aproximación, y para encontrar la solución del sistema al cabo del primer intervalo de tiempo, iteraremos hasta que no haya variación en las primeras cuatro cifras decimales. Esto ultimo equivale a definir un cierto grado de precisión en la solución. Para la primera iteración, nos queda:

u1( 1 ) = u2( 1 ) = u3( 1 ) = u4( 1 ) = u5( 1 ) =

1 ⋅ 4 1 ⋅ 4 1 ⋅ 4 1 ⋅ 4 1 ⋅ 2

(u (u (u (u (u

)

1 ⋅ ( 1 + 1) = 0.5 4 1 (0 ) (0 ) 1 + u3 + 2 = ⋅ ( 1 + 1 + 2 ) = 1 4 1 (0 ) (0 ) 2 + u4 + 2 = ⋅ ( 1 + 1 + 2 ) = 1 4 1 (0 ) (0 ) 3 + u5 + 2 = ⋅ ( 1 + 1 + 2 ) = 1 4 1 (0 ) 4 + 1 = ⋅ ( 1 + 1) = 1 2 (0 ) 2

+1 =

) ) )

)

Para la segunda iteración: u1( 2 ) = u2( 2 ) = u3( 2 ) = u4( 2 ) = u5( 2 ) =

1 ⋅ 4 1 ⋅ 4 1 ⋅ 4 1 ⋅ 4 1 ⋅ 2

(u (u (u (u (u

)

1 ⋅ ( 1 + 1) = 0.5 4 1 (1) (1) 1 + u3 + 2 = ⋅ ( 0.5 + 1 + 2 ) = 0.875 4 1 (1) (1) 2 + u4 + 2 = ⋅ ( 1 + 1 + 2 ) = 1 4 1 (1) (1) 3 + u5 + 2 = ⋅ ( 1 + 1 + 2 ) = 1 4 1 (1) 4 + 1 = ⋅ ( 1 + 1) = 1 2 (1) 2

+1 =

)

) ) )

UNIVERSIDAD INDUSTRIAL DE SANTANDER INGENIERIA DE PETROLEOS

Pág.67

Para la tercera iteración:

u1( 3 ) = u2( 3 ) = u3( 3 ) = u4( 3 ) = u5( 3 ) =

1 ⋅ 4 1 ⋅ 4 1 ⋅ 4 1 ⋅ 4 1 ⋅ 2

(u (u (u (u (u

)

1 ⋅ ( 0.875 + 1) = 0.46875 4 1 (2) (2) 1 + u3 + 2 = ⋅ ( 0.5 + 1 + 2 ) = 0.875 4 1 (2) (2) 2 + u4 + 2 = ⋅ ( 0.875 + 1 + 2 ) = 0.96875 4 1 (2) (2) 3 + u5 + 2 = ⋅ ( 1 + 1 + 2 ) = 1 4 1 (2) 4 + 1 = ⋅ ( 1 + 1) = 1 2 (2) 2

+1 =

) ) )

)

Luego de 11 iteraciones llegamos a la solución con la aproximación requerida, lo que se resume en la siguiente tabla:

MÉTODO

DE

JACOBI

Nº de iter.

u=0

u=1

u=2

u=3

u=4

u=5

n=0

0.0000

1.0000

1.0000

1.0000

1.0000

1.0000

n=1

0.0000

0.5000

1.0000

1.0000

1.0000

1.0000

n=2

0.0000

0.5000

0.8750

1.0000

1.0000

1.0000

n=3

0.0000

0.4688

0.8750

0.9688

1.0000

1.0000

n=4

0.0000

0.4688

0.8594

0.9688

0.9922

1.0000

........

........

........

........

........

........

........

n = 10

0.0000

0.4641

0.8564

0.9614

0.9890

0.9946

n = 11

0.0000

0.4641

0.8564

0.9613

0.9890

0.9945

n = 12

0.0000

0.4641

0.8564

0.9613

0.9890

0.9945

Hasta aquí hemos encontrado la solución para el primer intervalo de tiempo. Para encontrar la solución para el próximo intervalo de tiempo utilizamos como condiciones iniciales a los resultados obtenidos en el primer intervalo y nos queda para resolver el siguiente sistema: 4 ⋅ u1 − u 2 = 0.8564 − u1 + 4 ⋅ u 2 − u 3 = 1.4254 − u 2 + 4 ⋅ u 3 − u 4 = 1.8454 − u 3 + 4 ⋅ u 4 − u 5 = 1.9558 − u 4 + 2 ⋅ u 5 = 0.9890 UNIVERSIDAD INDUSTRIAL DE SANTANDER INGENIERIA DE PETROLEOS

Pág.68

Nuevamente,

para

resolver

este

sistema

deberemos

plantear

una

primera

aproximación, que usualmente son las condiciones iniciales para este intervalo (es decir la solución del primer sistema de ecuaciones), y comenzamos con un nuevo proceso iterativo hasta encontrar la solución al segundo sistema con la precisión deseada. Y así sucesivamente.

b) Resolución del sistema mediante el método de Gauss - Seidel Las ecuaciones de iteración son:

u1( n +1 ) = u2( n +1 ) = u3( n +1 ) = u4( n +1 ) = u5( n +1 ) =

1 ⋅ 4 1 ⋅ 4 1 ⋅ 4 1 ⋅ 4 1 ⋅ 2

(u (u (u (u (u

(n) 2

)

+1

) + 2) + 2)

( n +1 ) 1

+ u3( n ) + 2

( n +1 ) 2

+ u4( n )

( n +1 ) 3

+ u5( n )

( n +1 ) 4

+1

)

De la misma manera que en el ejemplo anterior partiendo de la primera aproximación, y para encontrar la solución del sistema al cabo del primer intervalo de tiempo, iteraremos hasta que no haya variación en las primeras cuatro cifras decimales. Parara la primera iteración, nos queda:

u1( 1 ) = u2( 1 ) = u3( 1 ) = u4( 1 ) = u5( 1 ) =

1 ⋅ 4 1 ⋅ 4 1 ⋅ 4 1 ⋅ 4 1 ⋅ 2

(u (u (u (u (u

)

1 ⋅ ( 1 + 1) = 0.5 4 1 (1) (0 ) 1 + u3 + 2 = ⋅ ( 0.5 + 1 + 2 ) = 0.875 4 1 (1) (0 ) 2 + u4 + 2 = ⋅ ( 0.875 + 1 + 2 ) = 0.96875 4 1 (1) (0 ) 3 + u5 + 2 = ⋅ ( 0.96875 + 1 + 2 ) = 0.99219 4 1 (1) 4 + 1 = ⋅ ( 0.99219 + 1) = 0.99609 2 (0 ) 2

+1 =

)

) ) )

UNIVERSIDAD INDUSTRIAL DE SANTANDER INGENIERIA DE PETROLEOS

Pág.69

Para la segunda iteración:

u1( 2 ) = u(22 ) = u3( 2 ) = u4( 2 ) = u5( 2 ) =

1 ⋅ 4 1 ⋅ 4 1 ⋅ 4 1 ⋅ 4 1 ⋅ 2

(u (u (u (u (u

)

1 ⋅ ( 0.875 + 1) = 0.46875 4 1 (2) (1) 1 + u3 + 2 = ⋅ ( 0.46875 + 0.96875 + 2 ) = 0.85938 4 1 (2) (1) 2 + u4 + 2 = ⋅ ( 0.85938 + 0.99219 + 2 ) = 0.96289 4 1 (2) (1) 3 + u5 + 2 = ⋅ ( 0.96289 + 0.99609 + 2 ) = 0.98975 4 1 (2) 4 + 1 = ⋅ ( 0.98975 + 1) = 0.99488 2 (1) 2

+1 =

)

) ) )

Para la tercera iteración:

u1( 3 ) = u(23 ) = u3( 3 ) = u4( 3 ) = u5( 3 ) =

1 ⋅ 4 1 ⋅ 4 1 ⋅ 4 1 ⋅ 4 1 ⋅ 2

(u (u (u (u (u

)

1 ⋅ ( 0.85938 + 1) = 0.46485 4 1 (3) (2) 1 + u3 + 2 = ⋅ ( 0.46485 + 0.96289 + 2 ) = 0.85694 4 1 (3) (2) 2 + u4 + 2 = ⋅ ( 0.85694 + 0.98975 + 2 ) = 0.96167 4 1 (3) (2) 3 + u5 + 2 = ⋅ ( 0.96167 + 0.99488 + 2 ) = 0.98914 4 1 (3) 4 + 1 = ⋅ ( 0.98914 + 1) = 0.99457 2 (2) 2

+1 =

)

) ) )

Luego de 5 iteraciones llegamos a la solución con la aproximación requerida, lo que se resume en la siguiente tabla: MÉTODO

DE

GAUSS SEIDEL

Nº de iter.

u=0

u=1

u=2

u=3

u=4

u=5

n=0

0.0000

1.0000

1.0000

1.0000

1.0000

1.0000

n=1

0.0000

0.5000

0.8750

0.9688

0.9922

0.9961

n=2

0.0000

0.4688

0.8594

0.9629

0.9898

0.9949

n=3

0.0000

0.4649

0.8569

0.9617

0.9891

0.9946

n=4

........

........

........

........

........

........

n=5

0.0000

0.4641

0.8564

0.9613

0.9890

0.9945

n=6

0.0000

0.4641

0.8564

0.9613

0.9890

0.9945

UNIVERSIDAD INDUSTRIAL DE SANTANDER INGENIERIA DE PETROLEOS

Pág.70

Hasta aquí hemos encontrado la solución para el primer intervalo de tiempo. Para encontrar la solución para el próximo intervalo de tiempo utilizamos como condiciones iniciales a los resultados obtenidos en el primer intervalo y nos queda para resolver el mismo sistema que determinamos en el primer ejemplo: 4 ⋅ u1 − u2 = 0.8564 − u1 + 4 ⋅ u2 − u3 = 1.4254 − u2 + 4 ⋅ u3 − u4 = 1.8454 − u3 + 4 ⋅ u4 − u5 = 1.9558 − u4 + 2 ⋅ u5 = .9890

El proceso de solución continua de la misma forma en que lo describimos en el ejemplo anterior donde aplicamos el método de Jacobi.

REFERENCIAS: •

http://personal.redestb.es/ztt/tem/t6_matrices.htm

•

http://www.ditutor.com/matrices/suma_matrices.html

•

http://www.monografias.com/trabajos72/resolucion-sistemas-metodo-gaussjordan/resolucion-sistemas-metodo-gauss-jordan.shtml

UNIVERSIDAD INDUSTRIAL DE SANTANDER INGENIERIA DE PETROLEOS

Pág.71