Capitulo 3 SOLUCION NUMERICA DE SISTEMAS DE ECUACIONES 105
II E1 II = Max I - 10 I 110 I = 10 II R1 II = Max I 0 I I - 5 I = 5
asf que un vector error residual pequeno (relativo al vector de terminos independientes
b = (1) II b 11 1 = 23 II b 112 21095 II b II = 21 ) corresponde a un vector error
relativamente grande
Para el sistema (34) X2 = (10) es la soluci6n exacta y si consideramos como una soluci6n00
aproximada a X- 2 = (34) que es la soluci6n exacta del sistema perturbado (34) entonces97
el vector error de X2 con respecto a X2 es
E2 =(34) _(10) =(-66)97 0 97
y el vector error residual correspondiente a la soluci6n aproximada X2 es
R2 =(411) _( 41) =(01J 97 97 0
Como
entonces nuevamente un vector error residual pequeno no corresponde a un vector error pequeno bull
EI ejemplo anterior pone de manifiesto que II R II pequeno no necesariamente implica que
II E II tambien sea pequeno Sin embargo a partir del siguiente teorema podremos probar
que satisfecha cierta condici6n 1111 ii pequeno implica ~ tambien pequenoshy
Teorema 33 Sea A E Rn xn una matriz no-singular y X la soluci6n exacta del sistema
AX = b b - O Si X es una soluci6n aproximada del sistema AX = b entonces para cualquier norma matricial inducida se tiene que
106 METODOS NUMERICOS
(37)
Demostraci6n Como R = AX - b = AX - AX = A(X - x) = AE YA es invertible entonces
E = A -1R b = AX X = A -1b Y aplicando la desigualdad (3 6) se obtiene
II R II $ 11~ 1111 E II es decir II II $ 11 E II y II E II $ 11 A -1 1111 R II
de donde
(38)1111 $11E II $11A -111 11R II
Aplicando la misma desigualdad (3 6) se tiene que
II b II $11 A 1111 X II es decir Ii 1111 $11 X II y II X II $11 A -1 1111 b II
de donde
1111 Iii $ 11 X II $ 11 A -11111 b IIbull o equivalentemente
-----1---- lt _1_ lt _II A_ II (39) II A -1 1111 b II - II X II - II b II
Combinando (3 8) y (3 9) obtenemos las siguientes cotas para el error relativo II ~ II en
terminos del error residual relativo ~ II b II
II R II 1 lt ~~JL II A 1111 A -111~ (310) II b 1111 A 11 11 A - 111 - 11 X 11 - II b II
que era 10 que querra demostrarse V
Dp acuerdo con este teorema 33 si se satisface la condici6n II A 11 11A-111 1 entonces ~ II b II
y I ~ II son mas 0 menos del mismo tamano As que si Iii 1111 es pequeno tambiEm 10 sera
II ~ II y si 11 ill es grande tambien 10 sera II ~ II por 10 tanto si II A 11 11A -111 1 podremos11
distinguir una soluci6n aproximada X buena de una
II R IIrelatlvo fblf
EI numero Cond (A) = II A 1111 A -1 II se lIamarcl NUMERO
de la matriz no-singular A relativo a la norma matricial
depende de la norma matricial usada sin embargo
matricial inducida pues
De acuerdo con la relaci6n (37) dada en el
entonces el error relativo 1 ~ II y el error resid
mismo tamano y podremos distinguir una observando el error residual relativo pero enfrI
De 10 anterior se espera que A tenga un buen residual relativo pequero implique cornesponcll
de AX = b si Cond (A) 1 caso en el ~
sistema AX = b esta bien condiciondo) 51 comportamiento en el sentido que un error una soluci6n aproximada mala ydiremos esta mal condicionado)
es poder determinar cuando una SOIUCuJ4I
tratar de distinguir si el sistema AX middotb
Para la matriz
del ejemplo 31 tenemos
luego
(37)
(38)
(39)
(310)
Capitulo 3 SOLUCION NUMERICA DE SISTEMAS DE ECUACIONES 107
distinguir una soluci6n aproximada X buena de una mala observando el error residual
relativo ~ b
EI numero Cond (A) = II A 1111 A -1 II se IIamara NUMERO DE CONDICION 0 CONDICIONAL
de la matriz no-singular A relativo a la norma matricial usada Aunque el valor de Cond(A)
depende de la norma matricial usada sin embargo Cond (A) ~ 1 cualquiera sea la norma
matricial inducida pues
In = AA -1 In II ~ II A 1111 A - 1 II y In II =~~o II In II ~ II =~~d ~ II =1
De acuerdo con la relaci6n (3 7) dada en el teorema 33 vemos que si Cond (A) 1
entonces el error relativo II ~ II y el error residual relativo II ii son mas 0 menos del
mismo tamano y podremos distinguir una soluci6n aproximada buena de una mala
observando el error residual relativo pero entre mas grande sea Cond(A) menor es la
informaci6n que se puede obtener del error relativo a partir del error residual relativo
De 10 anterior se espera que A tenga un buen comportamiento en el sentido de que un error residual relativo pequeno implique correspondientemente una buena soluci6n aproximada
de AX =b si Cond (A) 1 caso en el cual diremos que A esta bien condicionada (el
sistema AX =b esta bien condicionado) Si Cond (A) gtgt 1 es posible que A tenga un mal
comportamiento en el sentido que un error residual relativo pequeno puede corresponder a una soluci6n aproximada mala y diremos que A esta mal condicionada (el sistema AX = b esta mal condicionado)
A pesar de las definiciones anteriores no debemos olvidar que 10 que realmente nos interesa
es poder determinar cuando una soluci6n aproximada Xde un sistema AX = b es buena y tratar de distinguir si el sistema AX = b esta bien 0 mal condicionado
Para la matriz
del ejemplo 31 tenemos
A- 1 _ 1 ( 10 -1J - 05 - 1005 1
I A =Max 111 + 111 11005 1 + 1 10 1 = Max 2 2005 = 2005
II A-1 II = _1 II ( 10 - 1J II =_1 1105 = 221 05 - 1005 1 05
co
luego
108 METODOS NUMERICOS
Cond - (A) == II A 1100 II A - 110) = (2005)(221) == 443105 raquo 1
pequeno puedeEste numero de cond icion nos dice que un error residual relativo
muy grande asi que A puede considerarse mal
cond icionada
Veamos que puede decirse en este caso de la calidad de la soluci6n aproximada
- ( 10)
corresponder a un error relativo
X == - 8 del sistema
X +y == 2 1005x+10y==21
Para este ejemplo tenemos
II RII 5 y Cond (A) == 443105
II b II en 21
asi que la desigualdad (37) dada en el teorema 33 se convierte en
esto es
_ Ilx - x1 11 537 x10 6 ~ II x 110) 0) ~ 1055
5 10 que indica que aunque el error residual relativo es pequeno 21 el numero de condicion
es tan grande (443105) que hace que la solucion calculada pueda tener un error relativo de
hasta 1055 asi que nada puede decirse de la cercania entre X y X bull
Instrucciones en DERIVE
NORMAJNF(A) Simplifica en la norma del maximo de la matriz A II A II
COND_INF(A) Simplifica en el numero de condici6n relativo a la norma del maximo de la
matriz A es decir simplifica en el numero Cond O) (A) == II f-II00 11 A - 1 11 00 0
Existe otro numero asociado con una matriz al cual se Ie denomina tambien numero de condicion A continuaci6n nos referiremos a tal numero
Del teorema 32 se sabe que p (A)sl AI
I
pero como los valores propios de A tiene que
con cr(A) == A ECAes 1
EI numero Cond(A
Para la matriz A ( 1~
numero de
Capitulo 3 SOLUCION NUMERICA DE SISTEMAS DE ECUACIONES 109
Del teorema 32 se sabe que p (A) II A II para toda norma matricial inducida asl que
Cond (A) = II A IIII A -1 ~p(A) p(A-1)
pero como los valores propios de A - 1 son los reciprocos de los valores propios de A se tiene que
Max I A I ) Ecr(A) ( )
Cond (A) ~ I I == Condo AMm A
AEcr( A)
con a(A) =AE C A esvalorpropiodeA espectro de A (Recuerde que
p(A -1) = Max I A I = 1 I ) ecr(A-1) Mm A
AEcr( A)
Max I A I EI numero Cond (A)= AEcr(A) se denomina numero de condici6n espectral de A
bull Min I A I AEcr( A)
Segun se acaba de probar Cond(A) ~ Cond(A)
(1 1J Para la matnz A = se bene que 1005 10
det(A_AI) = 1- A 1 I =A2- 11A- 05 1005 10 - A
as que los valores propios de A son A1 1100454358 A 2 -45435778 x 10 -3 Y por tanto
Condo (A) 1100454358 2421999592 raquo 1 bull 45435778 x 10-3
Dado un sistema AX = b si 8A Y 8b denotan perturbaciones en A y b respectivamente el siguiente teorema cuya demostracion puede ser consultada en Ortega 1990 paginas 32 y
33 establece una cota para el error relativo II X- X II en terminos de las perturbaciones
IIXII
11 MII 118b II relatlvas lfAl lfbf y Cond(A) donde X es la soluclon exacta de AX = b Y X es la
solucion exacta del sistema perturbado (A + 8 A)X = b + 8 b
110 METODOS NUMERICOS
Teorema 34 Supongase que A es no-singular y que 11 8 A II lt II A~ II (esta hipotes is asegura
que A + 8 A es invertible y que 1- Cond(A) li SA II gt 0) Si X es la solucion exacta del IIA II
sistema perturbado (A + 0 A)X = b + 8 b entonces X aproxima a la solucion exacta X del
sistema AX = b b cI 0 con la siguiente estimacion de error
II X - X II ~ Cond (A) ( IL~~JL~ J (311)
IIXII 1_cond(A)[ 13 All l llb ll IIAII
II AII
La desigualdad (3 11) dice que si la matriz A esta bien condicionada es decir si
Cond (A) 1 entonces cambios pequerios en A y b producen correspondientemente
cambios pequerios en la soluci6n del sistema (el sistema AX = b esta bien condicionado) Por otro lado si A esta mal condicionada entonces cam bios pequenos en A y b pueden producir grandes cam bios en la solucion del sistema (el sistema AX = b esta mal condicionado)
Ejercicio 31 Estime la cota de error dada en el teorema 34 para los sistemas (34) y (34) del ejemplo 31 bull
Ejercicio 32 a) Calcule Cond (A) usando 11middot11 2 1111 y II to para las siguientes matrices
456 218) ( 279 138
b) Que puede decir del condicionamiento de los siguientes sistemas de ecuaciones lineales
f39x + 16x2 = 55 456X l + 218x2 = 674 I) ii) bulll 68x + 29x2 = 97 279x + 138x 2 = 417
ESTABILIDAD NUMERICA EN LA ELiMINACION GAUSSIANA
Volvamos al metodo de eliminacion Gaussiana (simple) y consideremos el siguiente ejemplo
Ejemplo 33 Resuelva el siguiente sistema de ecuaciones lineales usando eliminacion Gaussiana con sustitucion regresiva y aritmetica (decimal) con redondeo a tres dfgitos
El 03x + 589x2 = 592 E2 531x - 610x2 = 470
Usando elimnaci6n Gaussiana obtenemos
03 589 592) ~~- (A b) = ( 531 - 610 470
y por sustituci6n regresiva
_ 592 - 589(101 592 Xl = 03
luego la solucion calculada es X (~2 1 =(shy
Instrucci6n en DERIVE
PIVOT(A i j) Usa oper~cione~e tiene matriz obtenida de la matnz A q
Que puede decir de la calidad de la
Para intentar responder esta
dadas par el teorema 3 3
03 589) Como A = ( 531 -610
los calculos se obtienen III
v
(311 )
Capitulo 3 SOlUCION NUMERICA DE SISTEMAS DE ECUACIONES 111
(A b) = ( 03 589 592) __( 5__) _m 2 _=1_77 ) ~) ( 03 589E2_ ~~ _( _ _ 592 ) 531 -610 470 0 - 10400 -10500
y por sustituci6n regresiva
x = -10500 =101 2 - 10400
)(1 = 592 - 589(101) = 592 - 595 = - 3 = - 100 D3 D3 D3
- (X1 J (-1 00)luego la soluci6n calculada es X = _ = 101x2
Instrucci6n en DERIVE
PIVOT(A i j) Usa operaciones elemetales de fila para Simplificar (0 aproXimar) en una matriz obtenida de la matriz A que tiene ceros en la columna j y por debajo de la fila i 0
Que puede decir de la calidad de la soluci6n aproximada X
IIX -XIIPara intentar responder esta pregunta encontremos las cotas para el error relativo
IIXII dadas por el teorema 33
03 589)Como A = entonces usando aritmetica con redondeo a tres digitos para todos( 531 - 610
los calculos se obtienen las siguientes aproximaciones
- 1 1 (-610 - 589)A =-shy-313 - 531 03
II A II = Max 589 114 = 589
A - 1 II = _ 1_ Max 650 534 = _ 1_ 650 = 208II 313 313
entonces Cond (A) = II A II II A -111 = (589)(208) = 123 que no es muy grande
comparado con uno asi que la matriz A puede considerarse bien condicionada
(Por ciertas consideraciones te6ricas sobre el numero de condici6n de una matriz A las cuales pueden ser consultadas en Burden 1985 paginas 481 y 482 cuando se trabaja en
aritmetica finita (decimal) con redondeo a t-digitos y Cond(A) ~ 10t se espera un mal
112 METOOOS NUMERICOS
comportamiento de A con respecto a la solucion de AX = b Y A se considera mal
condicionada En este ejemplo Cond(A) = 123 lt 103 )
Ahora
03 589 ) (-100) ( 592)AX-b e ( 531 - 610 101 470
Estas operaciones se realizari en doble preci sion (6 digitos)
59189) ( 592) ( - 011 )(= -59261 - 470 = -106261
(Para evitar la perdida de cifras sign ificativas se debe calcular el vector erro r res idual
R = AX - b en doble precision)
Conv irtiendo este ultimo resultado a tres digltos usando redondeo se obtiene
R = (-011) - 106
Entonces II R II =106 y II b II = 592 Y por tanto
II R la 1 106 1 11gtlt - X II 106 II R II en
-II-11- ( )= - - =145 ~ II II ~ 2202= 123- = Cond (A)-II - lI shyb co Cond A 592 123 X 592 b DO
II R pero como Cond (A) ~ 1 entonces se espera que fbiC Y X- X sean mas 0 menosII xL
II X- X L sea tamble ndel mismo tamano y ya que es grande se espera que II X II
grande +
En situaciones como a observada en este ejemplo se sugiere hacer un refinamiento
iterativo sobre la soluci6n caculada X 0 usar esta soluci6n calcuada como aproximaci6n inicial en un metodo iterativ~ con et proposito de tratar de mejorar la soluci6n aproximada y ograr que el error residua l relativo sea mas pequeno EI metodo de reflna81iento iterativo puede ser consutado en Kincaid 1972 pag inas 174-176
La soucion exacta del sistema en consideracion es X = (J= ( 0) A que se debe a
diferencia entre la soucion exacta y a solucion calculada
Observe que e error en e calculo de )(2 con respecto a x2 fue de solo 01 (un error relativo
de 1) Y este error fue multiplicado par un factor de aproximadamente - 2000 al obtener
)(1 debido a orden en que se realiz6 la el iminaci6n Gaussiana
Instruccion en DERIVE
RESUELVA_1(Ab) Simplifica en a SOUCI6n entra como un vector fila 0
En el ejemplo anterior la elimina~i6n sistema de ecuaciones lineales ble~ nII1II
del algoritmo de eliminaci6n Gausslana pequeno) Hay sin embargo SIlUC -
Gaussiana es numencamente estable consultarse en Burden 1985 paginas366 y
Teorema 35 Si A =(ajiLn (EDD) par filas es decir si
entonces A es invertible (no-singular intercambio de filas en cualquier calculos son estables con respeclO
Notese que como consecuencia
por filas entonces A t~ene unos en su diagonal pnnClpalyU
~~a~ (vease Burden 19a5
satisface
Capitulo 3 SOLUCION NUMERIC A DE SISTEMAS DE ECUACIONES 113
Instrucci6n en DERIVE
RESUELVA_1(Ab) Simplifica en la soluci6n exacta X del sistema AX=b EI vector b se entra como un vector fila 0
En el ejemplo anterior la eliminaci6n Gaussiana condujo a una respuesta defectuosa de un sistema de ecuaciones lineales bien condicionado Esto muestra la inestabilidad numerica del algoritmo de eliminaci6n Gaussiana (consecuencia de la divisi6n por un nOmero (pivote) pequero) Hay sin embargo situaciones en las cuales el algoritmo de eliminaci6n Gaussiana es numericamente estable EI siguiente teorema cuya demostraci6n puede consultarse en Burden 1985 paginas 366 y 367 se refiere a una de tales situaciones
Teorema 35 Si A = (ai j )n xn es una matriz estrictamente dominante diagonalmente
(EDD) por filas es decir si
n Iaii Igt II ai j I para cada i == 12 n
)= 1 ) 1
entonces A es invertible (no-singular) Ademas se puede rea lizar eliminaci6n Gaussiana sin intercambio de filas en cualquier sistema AX = b para obtener su Onica soluci6n y los calculos son estables con respecto al crecimiento de los errores de redondeo V
N6tese que como consecuencia del teorema anterior se tiene que Si A = (a l ) )n xn es EDD
por filas entonces A tiene factorizaci6n LU es decir A == LU con L triangular inferior con unos en su diagonal principal y U triangular superior (escalonada)
sea tambien
Observe que la matriz de coeficientes del ejemplo 33 anterior no es EDD por filas
sectJeor9~~am9i~ -valioo=paramatrj simetricas
(vease Bllrden 19 5 ~ina 368) Una matriz A E Rnxn ~ i~etrica se dice definida
positiva si satisface una cualquiera de las sig uientes condiciones (Ias- cualesson equival~tes) a _-----
- - -- shy
i) ~AX gt 0 paralodo 2CERn X 0 J
iiXTodos los valores propios de A son pos1tlvOS7
iii) Todos los pivotes obtenidos en la eliminaci6n Gaussiana sobre A sin intercambio de filas - ~on positiv~s
iv) Todas las submatrices principales de A tienen determinante positiv~
(Las submatrices principales de la matriz A == (ai J) son las matrices - n ~ n
1k a ] k == 12bull f) )
akk v
114 METODOS NUMERICOS
N6tese nuevamente que Si A E R nxn es simetrica y definida positiva entonces A tiene
factorizaci6n A = LU con L triangular inferior con sus componentes sabre la diagonal principal iguales a uno y U triangular superior (ver mas adelante factorizaci6n de Choleski)
Observe que la matriz de coeficientes de l ejemplo 33 anterior no es simetrica
34 ESTRATEGIAS DE PIVOTEO
EI eJemplo 33 anterior muestra una de las dificultades que pueden surgir al aplicar el metodo
de eliminaci6n Gaussiana cuando el pivote ajj-1) es pequeno comparado con algunos
(j-1) e ementos I a t para J ~ I t ~ n
Para tratar de evitar tales dificultades se introduce en el metodo de eliminaci6n Gaussiana una estrategia IIamada de pivoteo la cual consiste en seleccionar el pivote de acuerdo con un cierto criter io Nosotros usaremos dos estrategias la estrategia de pivoteo maximo por columna 0 pivoteo parcial y la estrategia de pivoteo escalado de fila 0 escalamiento
341 Pivoteo maximo por columna 0 pivoteo parcial Esta estrategia difiere de
eliminaci6n Gaussiana simple unicamente en la escogencia del pivote ajj-1) la cual se hace
ahora as
Para j = 12 n - 1 se determina el menor entero k J ~ k ~ n tal que
y
es decir seleccionamos el primer elemento diferente de cero sobre la columna j-esima a
partir de la j-esima fila y que tenga mayor valor absoluto (para j = 1 a~j j1 ) = al~) ak 1 )
Si tal k no existe el sistema no tiene soluci6n unica y el proceso se puede terminar
Si tal k existe y k ot j entonces hacemos intercambio de las ecuaciones j-esima y k-esima
y continuamos con la eliminaci6n Gaussiana l
Ilustremos esta estrategia para resolver el sistema
E1 03X1 + 589x2 = 592 E2 531x1 - 61 OX 2 = 470
entonces k = 2 ot 1= j asl que
Por sustituci6n regresiva
Observe que en este caSO
SWAP(A i j)
que es el mismo del ejemplo 33 usando aritmetica con redondeo a tres digitos
Como para j = 1 se tiene que
Capitulo 3 SOLUCION NUMERICA DE SISTEMAS DE ECUACIONES 115
Max Iall I Ia21 I = Max 03 531 = 531 = Ia21 I 0
entonces k = 2 1= j as i que intercambiamos El y E2 Y continuamos con la el imi nac i6n
531 - 610 (A b)= ( 03 589 592) 470) (531 - 610 470 03 589 592
E2I (~) ( m21=5 65 ~ 1O - l) 531 - 610 5 31 470) ------------------~) ( o 589 589
Por sustituci6n regresiva obtenemos
x = 589 = 100 Xl = 470 +610(100) = 531 = 100 2 589 531 531
Observe que en este caso X= (J es la soluci6n exacta del sistema dado bull
Instrucci6n en DERIVE
SWAP(A i j) Intercambia las filas (0 elementos) i y j de la matriz A (de un vector) 0
Nota En el procedimiento de pivoteo maximo por columna (p ivoteo parcial) cada multiplicador mi j es tal que
y aunque esta estrategia permite resolver satisfactoriamente muchos sistemas de ecuaciones lineales hay casos donde fracasa como se ilustra en el siguiente ejemplo
Ejemplo 34 Consideremos el sistema
El 300Xl + 58900x2 59200 E2 531x1 -610X2 - 470
el cual es un sistema equivalente al del ejemplo 33 (los coeficientes de la primera ecuaci6n
en el sistema del ejemplo 33 han sido multiplicados por 10 3) EI pivoteo maximo por
columna con aritmetica de redondeo a tres dig itos nos lIeva a los siguientes resultados
59200) E2 1 (-~pound) (m 2 1 =177) ) ( 300 58900 59200 ) (A b) = ( 300 58900 531 - 610 470 o - 10400 -10500
y por sustituci6n regresiva x2 = 101 Y Xl = - 100 que es la misma solucion que se obtiene si ~ usamos eliminaci6n Gaussiana simple
116 METODOS NUMERICOS
En casas como el de este ejemplo donde un pivote es mucho mas pequeno que alguno de los coeficientes de la ecuacion que el encabeza se recomienda la tecnica conocida como pivoteo escalado de fila 0 escalamiento la cual es nuestra segunda estrategia
342 Pivoteo escalado de fila Esta tecnica solo difiere de la eliminaci6n Gaussiana simple al igual que el pivoteo parcial en la escogencia del pivote
lEsta vez el pivote alr ) se escoge como se indica a continuacion
Para j = 12 n - 1 hacemos 10 siguiente
a) Para i = j j + 1 n calculamos
S = Max Ia(l- l) I Factor de escala I j s i s n II
Si Sj = 0 entonces el sistema no tiene solucion unica y el proceso se puede terminar
b) Para i = jj+1 n calculamos
c) Encontramos el menor entero k con j o k 0 n tal que
Ia~jjl ) I IaH- l) I -------- = Max --------
Sk j S j $ n SI
Si tal k no existe entonces el sistema no tiene solucion unica y el proceso se puede terminar
Si tal k existe y kF j entonces hacemos intercambio de las ecuaciones j-esima y kshyesima
y continuamos con la eliminacion Gaussiana V
Apliquemos esta estrategia para resolver el sistema del ejemplo 34 usando aritmetica con redondeo a tres dfgitos
Para j = 1
a) Sl = Max Iall I I a12 I = Max 300 58900 = 58900F 0 Y
Max 53~ 610 IshyS2 =Max a21 a n -shy
b) Ahora ~= 51
la0
l 300
S2
531 al l I a21 -Max - 0 -
c) Max -S -S- - 5890 NIl 1 2
las ecuaciolllltanto intercamblamos
Gaussiana
300 58900 (A b) = ( 531 -610
Por sustitucion regresiva
Xl tOO
Observe que en este GaSO X
35 FACTORIZACION
Capitulo 3 SOLUCION NUMERICA DE SISTEMAS DE ECUACIONES 117
8 2 =Max I a21 I I a22 1 =Max 531 610 = 610 0
b) Ahora
Ia I 300 - 8- = 58900 Y
1
Ia 2 I 531 - - = - shy
6108 2
c) Max ~ ~ I=Max ~ 531 = 531 =~ 0 as que k = 2 1= j y por8 8 2 58900 60 60 8 2
tanto intercambiamos las ecuaciones E y E2 Y continuamos con la eliminacion
Gaussiana
(A b) = (300 58900 59200) p ( 531 -610 470 ) 2 )
531 - 610 470 300 58900 59200
E2 (300) (m2 =565 ) ( 531 -6 1 0531 47 0 )
--------------~) o 58900 58900
Por sustituci6n regresiva
_ 470 +610(100) 531X2 = 100 Xl = =-- = 100
531 531
Observe que en este caso X= ( ~ J = ( 100) es la solucion exacta del sistema bull 100x2
35 FACTORIZACION TRIANGULAR
Consideremos un sistema AX = b con A no-singular y b 0 Con respecto a la matriz A se sabe que existen matrices P de permutaci6n L triangular inferior con sus componentes sobre la diagonal principal iguales a uno y U triangular superior (escalon ada) tales que PA = LU Una forma eficiente computacionalmente de encontrar P L Y U usando eliminaci6n Gaussiana se muestra en el siguiente ejemplo Ejemplo 35 Resuelva el siguiente sistema de ecuaciones lineales usando eliminacion Gaussiana con pivoteo parcial y obtenga una factorizaci6n PA =LU para la matriz A de coeficientes asociada con este metodo
- x + 2X2 - 3X3 = -2 3x1 - 3x2 - X 3 = -4
x + x2 = 3IEmpezamos introduciendo un vector p =(PP2P3) T el cual inicializamos_con P = i i = 123 Y
donde se almacenaran los intercambios necesarios en el proceso de eliminacion Gaussiana
118 METODOS NUMERICOS
con pivoteo parcial (el numero de componentes del vector p coincide con el orden n del sistema a resolver)
- 1 2 -3 -2] -4] - 2 P = (123) T (Ab) = ~ - 3 -1 - 4
[ o 3 3
3 - 3 - 1 - 4(m=-H (m3l = ~ )
E( - ~) E3(~) 10 10(-i) 3 3
2(i) 1 13 3 3
(Observe que cada multiplicador mij es almacenado en la posici6n correspondiente (i j) en
la matriz de trabajo)
3 - 3 -1 - 4
Max 121 =2 p=( 23 1) T bull (i) 2 1
-
3 13 3
(-i) 10 3
10 3
(Observe que la permutaci6n se hace para las filas 2 y 3 completas es decir incluyendo los multiplicadores)
3 - 3 - 1 --4(mJ2 = ~ ) E3( ~)
2(i) 13 3 3
(- i) (i) 7 11 2 2
La eficiencia en el metoda indicado se debe a que en la misma matriz de trabajo se almacenan los multiplicadores que van a conformar la matriz L (en el ejemplo son los numeros que se encuentran dentro de parentesis) 10 que significa un ahorro de memoria y como los intercambios necesarios afectan simultaneamente a las matrices L y U se evita tener que volver a la matriz original a realizar los intercambios ya observados y repetir la eliminaci6n Gaussiana con pivoteo parcial De esta manera al terminar el proceso de eliminaci6n podemos leer en la matriz final la parte estrictamente triangular inferior de L (son los numeros entre parentesis) y la matriz triangular superior U (que es la parte triangular superior de la matriz final) y en el vector p final quedan almacenados los intercambios realizados que se usan para producir la matriz de permutaci6n P
Para el ejemplo 35
3 -31 0 0 1 0 2U=0L= 3 1 1 0 0-3 2
(Verifique que PA =LU)
Para obtener la soluci6n del sistema
reducido
yobtenemos
C6mo se resuelve el
Capitulo 3 SOLUCION NUMERICA DE SISTEMAS DE ECUACIONES 119
Para el ejemplo 35
1 0 0 3 -3 - 1 1 1 1
L= - 0 U= 0 2 - Y como P~ (231) entonces [~ 0 lp3 3 1 1 7 0
-- - 0 0 3 2 2
(Verifique que PA = LU )
Para obtener la soluci6n del sistema original usamos sustituci6n regresiva en el sistema reducido
-4 13
UX = 3
11
2 ~
c y obtenemos
asi que la soluci6n (exacta) del sistema dado es X = (23 40 ~)T21 21 7
C6mo se resuelve el sistema AX = b a partir de la factorizaci6n PA = LU obtenida bull
351 Algunas aplicaciones de la factorizaci6n PA = LU La factorizaci6n PA = LU es utilizada eficientemente en aquellos casos donde se trabaja repetfdamente con la misma matriz A Dos de esos casos se presentan a continuaci6n
1) Resolver varios sistemas AX =b con la misma matriz de coeficientes A ya que en p L Y U esta almacenado todo el proceso de eliminaci6n Gaussiana EI algoritmo se bas a en la siguiente equivalencia
AX = b lt=gt PAX = Pb lt=gt LUX = Pb
lt=gt UX = L- 1Pb
lt=gtUXoo c y c oo L- 1Pb
lt=gt UX == C Y Lc = Pb lt=gt Lc = Pb Y UX = c
Los pasos a seguir son
Paso 1 Calcular Pb Paso 2 Resolver para c Lc = Pb por sustituci6n progresiva Paso 3 Resolver para X UX = c por sustituci6n regresiva
Como ejercicio resuelva el sistema del ejemplo anterior usando este algoritmo
120 METODOS NUMERICOS
2) Encontrar la matriz inversa A -1 de una matriz invertible A resolviendo los n-sistemas
AX =eU) j = 12 n
donde e(1) =(O O1 o 0f ERn
t posicion j
La solucion X del sistema AX = e(n j = 12 n produce la correspondiente columna j shy
esima de la matriz A -1
Como ejercicio compare el numero de operaciones para encontrar A - usando el metodo de Gauss-Jordan con el numero de operaciones resolviendo los n-sistemas indicados antes V
Ejercicio 33 Calcule la inversa de la matriz A de coeficientes del sistema del ejemplo 35 usando el metodo de Gauss-Jordan y tambien usando la factorizacion PA = LU bull
36 SISTEMAS TRIDIAGONALES
Un caso muy importante de sistemas de ecuaciones lineales que requiere un tratamiento especial es el de los sistemas tridiagonales Tales sistemas aparecen en diversas aplicaciones como por ejemplo al utilizar metodos de diferencias flnitas en la solucion de problemas con valores en la frontera para ecuaciones diferenciales ordinarias y como veremos mas adelante en el problema de la interpolacion segmentaria cubica
Un sistema tridiagonal es de la forma
dx1 + C1X2
a 2 x + d2 x2 + C 2 X 3
a3 x2 + d3 x3 + C 3 X 4
an-1Xn- 2 + dn_ xn_ + cn_xn = bn_1
anxn_1 + dnxn = bn
La matriz de coeficientes del sistema es
C 0 0 od1
0a 2 d2 c2
o a3 d3 c 3
A =
o
la cual se dice una matriz TRIDIAGONAL
En general
ij=12 n
suponiendo que la matriz A de podemos usar eliminaclon Gausslana eliminacion Gaussiana simple es adecuadl EDD por filas) Otra forma de resolver
EDD por filas es a partir de la componentes de las matrices L Y U as triangular inferior con todos sus superior Para encontrar tales malnce5
o 0 0
Y2 o o o o
L =
o Yn t
o o
Como Cl Ct
Y2 Cl Y2C +az
0 h a2
LU =
0
0
d
=
106 METODOS NUMERICOS
(37)
Demostraci6n Como R = AX - b = AX - AX = A(X - x) = AE YA es invertible entonces
E = A -1R b = AX X = A -1b Y aplicando la desigualdad (3 6) se obtiene
II R II $ 11~ 1111 E II es decir II II $ 11 E II y II E II $ 11 A -1 1111 R II
de donde
(38)1111 $11E II $11A -111 11R II
Aplicando la misma desigualdad (3 6) se tiene que
II b II $11 A 1111 X II es decir Ii 1111 $11 X II y II X II $11 A -1 1111 b II
de donde
1111 Iii $ 11 X II $ 11 A -11111 b IIbull o equivalentemente
-----1---- lt _1_ lt _II A_ II (39) II A -1 1111 b II - II X II - II b II
Combinando (3 8) y (3 9) obtenemos las siguientes cotas para el error relativo II ~ II en
terminos del error residual relativo ~ II b II
II R II 1 lt ~~JL II A 1111 A -111~ (310) II b 1111 A 11 11 A - 111 - 11 X 11 - II b II
que era 10 que querra demostrarse V
Dp acuerdo con este teorema 33 si se satisface la condici6n II A 11 11A-111 1 entonces ~ II b II
y I ~ II son mas 0 menos del mismo tamano As que si Iii 1111 es pequeno tambiEm 10 sera
II ~ II y si 11 ill es grande tambien 10 sera II ~ II por 10 tanto si II A 11 11A -111 1 podremos11
distinguir una soluci6n aproximada X buena de una
II R IIrelatlvo fblf
EI numero Cond (A) = II A 1111 A -1 II se lIamarcl NUMERO
de la matriz no-singular A relativo a la norma matricial
depende de la norma matricial usada sin embargo
matricial inducida pues
De acuerdo con la relaci6n (37) dada en el
entonces el error relativo 1 ~ II y el error resid
mismo tamano y podremos distinguir una observando el error residual relativo pero enfrI
De 10 anterior se espera que A tenga un buen residual relativo pequero implique cornesponcll
de AX = b si Cond (A) 1 caso en el ~
sistema AX = b esta bien condiciondo) 51 comportamiento en el sentido que un error una soluci6n aproximada mala ydiremos esta mal condicionado)
es poder determinar cuando una SOIUCuJ4I
tratar de distinguir si el sistema AX middotb
Para la matriz
del ejemplo 31 tenemos
luego
(37)
(38)
(39)
(310)
Capitulo 3 SOLUCION NUMERICA DE SISTEMAS DE ECUACIONES 107
distinguir una soluci6n aproximada X buena de una mala observando el error residual
relativo ~ b
EI numero Cond (A) = II A 1111 A -1 II se IIamara NUMERO DE CONDICION 0 CONDICIONAL
de la matriz no-singular A relativo a la norma matricial usada Aunque el valor de Cond(A)
depende de la norma matricial usada sin embargo Cond (A) ~ 1 cualquiera sea la norma
matricial inducida pues
In = AA -1 In II ~ II A 1111 A - 1 II y In II =~~o II In II ~ II =~~d ~ II =1
De acuerdo con la relaci6n (3 7) dada en el teorema 33 vemos que si Cond (A) 1
entonces el error relativo II ~ II y el error residual relativo II ii son mas 0 menos del
mismo tamano y podremos distinguir una soluci6n aproximada buena de una mala
observando el error residual relativo pero entre mas grande sea Cond(A) menor es la
informaci6n que se puede obtener del error relativo a partir del error residual relativo
De 10 anterior se espera que A tenga un buen comportamiento en el sentido de que un error residual relativo pequeno implique correspondientemente una buena soluci6n aproximada
de AX =b si Cond (A) 1 caso en el cual diremos que A esta bien condicionada (el
sistema AX =b esta bien condicionado) Si Cond (A) gtgt 1 es posible que A tenga un mal
comportamiento en el sentido que un error residual relativo pequeno puede corresponder a una soluci6n aproximada mala y diremos que A esta mal condicionada (el sistema AX = b esta mal condicionado)
A pesar de las definiciones anteriores no debemos olvidar que 10 que realmente nos interesa
es poder determinar cuando una soluci6n aproximada Xde un sistema AX = b es buena y tratar de distinguir si el sistema AX = b esta bien 0 mal condicionado
Para la matriz
del ejemplo 31 tenemos
A- 1 _ 1 ( 10 -1J - 05 - 1005 1
I A =Max 111 + 111 11005 1 + 1 10 1 = Max 2 2005 = 2005
II A-1 II = _1 II ( 10 - 1J II =_1 1105 = 221 05 - 1005 1 05
co
luego
108 METODOS NUMERICOS
Cond - (A) == II A 1100 II A - 110) = (2005)(221) == 443105 raquo 1
pequeno puedeEste numero de cond icion nos dice que un error residual relativo
muy grande asi que A puede considerarse mal
cond icionada
Veamos que puede decirse en este caso de la calidad de la soluci6n aproximada
- ( 10)
corresponder a un error relativo
X == - 8 del sistema
X +y == 2 1005x+10y==21
Para este ejemplo tenemos
II RII 5 y Cond (A) == 443105
II b II en 21
asi que la desigualdad (37) dada en el teorema 33 se convierte en
esto es
_ Ilx - x1 11 537 x10 6 ~ II x 110) 0) ~ 1055
5 10 que indica que aunque el error residual relativo es pequeno 21 el numero de condicion
es tan grande (443105) que hace que la solucion calculada pueda tener un error relativo de
hasta 1055 asi que nada puede decirse de la cercania entre X y X bull
Instrucciones en DERIVE
NORMAJNF(A) Simplifica en la norma del maximo de la matriz A II A II
COND_INF(A) Simplifica en el numero de condici6n relativo a la norma del maximo de la
matriz A es decir simplifica en el numero Cond O) (A) == II f-II00 11 A - 1 11 00 0
Existe otro numero asociado con una matriz al cual se Ie denomina tambien numero de condicion A continuaci6n nos referiremos a tal numero
Del teorema 32 se sabe que p (A)sl AI
I
pero como los valores propios de A tiene que
con cr(A) == A ECAes 1
EI numero Cond(A
Para la matriz A ( 1~
numero de
Capitulo 3 SOLUCION NUMERICA DE SISTEMAS DE ECUACIONES 109
Del teorema 32 se sabe que p (A) II A II para toda norma matricial inducida asl que
Cond (A) = II A IIII A -1 ~p(A) p(A-1)
pero como los valores propios de A - 1 son los reciprocos de los valores propios de A se tiene que
Max I A I ) Ecr(A) ( )
Cond (A) ~ I I == Condo AMm A
AEcr( A)
con a(A) =AE C A esvalorpropiodeA espectro de A (Recuerde que
p(A -1) = Max I A I = 1 I ) ecr(A-1) Mm A
AEcr( A)
Max I A I EI numero Cond (A)= AEcr(A) se denomina numero de condici6n espectral de A
bull Min I A I AEcr( A)
Segun se acaba de probar Cond(A) ~ Cond(A)
(1 1J Para la matnz A = se bene que 1005 10
det(A_AI) = 1- A 1 I =A2- 11A- 05 1005 10 - A
as que los valores propios de A son A1 1100454358 A 2 -45435778 x 10 -3 Y por tanto
Condo (A) 1100454358 2421999592 raquo 1 bull 45435778 x 10-3
Dado un sistema AX = b si 8A Y 8b denotan perturbaciones en A y b respectivamente el siguiente teorema cuya demostracion puede ser consultada en Ortega 1990 paginas 32 y
33 establece una cota para el error relativo II X- X II en terminos de las perturbaciones
IIXII
11 MII 118b II relatlvas lfAl lfbf y Cond(A) donde X es la soluclon exacta de AX = b Y X es la
solucion exacta del sistema perturbado (A + 8 A)X = b + 8 b
110 METODOS NUMERICOS
Teorema 34 Supongase que A es no-singular y que 11 8 A II lt II A~ II (esta hipotes is asegura
que A + 8 A es invertible y que 1- Cond(A) li SA II gt 0) Si X es la solucion exacta del IIA II
sistema perturbado (A + 0 A)X = b + 8 b entonces X aproxima a la solucion exacta X del
sistema AX = b b cI 0 con la siguiente estimacion de error
II X - X II ~ Cond (A) ( IL~~JL~ J (311)
IIXII 1_cond(A)[ 13 All l llb ll IIAII
II AII
La desigualdad (3 11) dice que si la matriz A esta bien condicionada es decir si
Cond (A) 1 entonces cambios pequerios en A y b producen correspondientemente
cambios pequerios en la soluci6n del sistema (el sistema AX = b esta bien condicionado) Por otro lado si A esta mal condicionada entonces cam bios pequenos en A y b pueden producir grandes cam bios en la solucion del sistema (el sistema AX = b esta mal condicionado)
Ejercicio 31 Estime la cota de error dada en el teorema 34 para los sistemas (34) y (34) del ejemplo 31 bull
Ejercicio 32 a) Calcule Cond (A) usando 11middot11 2 1111 y II to para las siguientes matrices
456 218) ( 279 138
b) Que puede decir del condicionamiento de los siguientes sistemas de ecuaciones lineales
f39x + 16x2 = 55 456X l + 218x2 = 674 I) ii) bulll 68x + 29x2 = 97 279x + 138x 2 = 417
ESTABILIDAD NUMERICA EN LA ELiMINACION GAUSSIANA
Volvamos al metodo de eliminacion Gaussiana (simple) y consideremos el siguiente ejemplo
Ejemplo 33 Resuelva el siguiente sistema de ecuaciones lineales usando eliminacion Gaussiana con sustitucion regresiva y aritmetica (decimal) con redondeo a tres dfgitos
El 03x + 589x2 = 592 E2 531x - 610x2 = 470
Usando elimnaci6n Gaussiana obtenemos
03 589 592) ~~- (A b) = ( 531 - 610 470
y por sustituci6n regresiva
_ 592 - 589(101 592 Xl = 03
luego la solucion calculada es X (~2 1 =(shy
Instrucci6n en DERIVE
PIVOT(A i j) Usa oper~cione~e tiene matriz obtenida de la matnz A q
Que puede decir de la calidad de la
Para intentar responder esta
dadas par el teorema 3 3
03 589) Como A = ( 531 -610
los calculos se obtienen III
v
(311 )
Capitulo 3 SOlUCION NUMERICA DE SISTEMAS DE ECUACIONES 111
(A b) = ( 03 589 592) __( 5__) _m 2 _=1_77 ) ~) ( 03 589E2_ ~~ _( _ _ 592 ) 531 -610 470 0 - 10400 -10500
y por sustituci6n regresiva
x = -10500 =101 2 - 10400
)(1 = 592 - 589(101) = 592 - 595 = - 3 = - 100 D3 D3 D3
- (X1 J (-1 00)luego la soluci6n calculada es X = _ = 101x2
Instrucci6n en DERIVE
PIVOT(A i j) Usa operaciones elemetales de fila para Simplificar (0 aproXimar) en una matriz obtenida de la matriz A que tiene ceros en la columna j y por debajo de la fila i 0
Que puede decir de la calidad de la soluci6n aproximada X
IIX -XIIPara intentar responder esta pregunta encontremos las cotas para el error relativo
IIXII dadas por el teorema 33
03 589)Como A = entonces usando aritmetica con redondeo a tres digitos para todos( 531 - 610
los calculos se obtienen las siguientes aproximaciones
- 1 1 (-610 - 589)A =-shy-313 - 531 03
II A II = Max 589 114 = 589
A - 1 II = _ 1_ Max 650 534 = _ 1_ 650 = 208II 313 313
entonces Cond (A) = II A II II A -111 = (589)(208) = 123 que no es muy grande
comparado con uno asi que la matriz A puede considerarse bien condicionada
(Por ciertas consideraciones te6ricas sobre el numero de condici6n de una matriz A las cuales pueden ser consultadas en Burden 1985 paginas 481 y 482 cuando se trabaja en
aritmetica finita (decimal) con redondeo a t-digitos y Cond(A) ~ 10t se espera un mal
112 METOOOS NUMERICOS
comportamiento de A con respecto a la solucion de AX = b Y A se considera mal
condicionada En este ejemplo Cond(A) = 123 lt 103 )
Ahora
03 589 ) (-100) ( 592)AX-b e ( 531 - 610 101 470
Estas operaciones se realizari en doble preci sion (6 digitos)
59189) ( 592) ( - 011 )(= -59261 - 470 = -106261
(Para evitar la perdida de cifras sign ificativas se debe calcular el vector erro r res idual
R = AX - b en doble precision)
Conv irtiendo este ultimo resultado a tres digltos usando redondeo se obtiene
R = (-011) - 106
Entonces II R II =106 y II b II = 592 Y por tanto
II R la 1 106 1 11gtlt - X II 106 II R II en
-II-11- ( )= - - =145 ~ II II ~ 2202= 123- = Cond (A)-II - lI shyb co Cond A 592 123 X 592 b DO
II R pero como Cond (A) ~ 1 entonces se espera que fbiC Y X- X sean mas 0 menosII xL
II X- X L sea tamble ndel mismo tamano y ya que es grande se espera que II X II
grande +
En situaciones como a observada en este ejemplo se sugiere hacer un refinamiento
iterativo sobre la soluci6n caculada X 0 usar esta soluci6n calcuada como aproximaci6n inicial en un metodo iterativ~ con et proposito de tratar de mejorar la soluci6n aproximada y ograr que el error residua l relativo sea mas pequeno EI metodo de reflna81iento iterativo puede ser consutado en Kincaid 1972 pag inas 174-176
La soucion exacta del sistema en consideracion es X = (J= ( 0) A que se debe a
diferencia entre la soucion exacta y a solucion calculada
Observe que e error en e calculo de )(2 con respecto a x2 fue de solo 01 (un error relativo
de 1) Y este error fue multiplicado par un factor de aproximadamente - 2000 al obtener
)(1 debido a orden en que se realiz6 la el iminaci6n Gaussiana
Instruccion en DERIVE
RESUELVA_1(Ab) Simplifica en a SOUCI6n entra como un vector fila 0
En el ejemplo anterior la elimina~i6n sistema de ecuaciones lineales ble~ nII1II
del algoritmo de eliminaci6n Gausslana pequeno) Hay sin embargo SIlUC -
Gaussiana es numencamente estable consultarse en Burden 1985 paginas366 y
Teorema 35 Si A =(ajiLn (EDD) par filas es decir si
entonces A es invertible (no-singular intercambio de filas en cualquier calculos son estables con respeclO
Notese que como consecuencia
por filas entonces A t~ene unos en su diagonal pnnClpalyU
~~a~ (vease Burden 19a5
satisface
Capitulo 3 SOLUCION NUMERIC A DE SISTEMAS DE ECUACIONES 113
Instrucci6n en DERIVE
RESUELVA_1(Ab) Simplifica en la soluci6n exacta X del sistema AX=b EI vector b se entra como un vector fila 0
En el ejemplo anterior la eliminaci6n Gaussiana condujo a una respuesta defectuosa de un sistema de ecuaciones lineales bien condicionado Esto muestra la inestabilidad numerica del algoritmo de eliminaci6n Gaussiana (consecuencia de la divisi6n por un nOmero (pivote) pequero) Hay sin embargo situaciones en las cuales el algoritmo de eliminaci6n Gaussiana es numericamente estable EI siguiente teorema cuya demostraci6n puede consultarse en Burden 1985 paginas 366 y 367 se refiere a una de tales situaciones
Teorema 35 Si A = (ai j )n xn es una matriz estrictamente dominante diagonalmente
(EDD) por filas es decir si
n Iaii Igt II ai j I para cada i == 12 n
)= 1 ) 1
entonces A es invertible (no-singular) Ademas se puede rea lizar eliminaci6n Gaussiana sin intercambio de filas en cualquier sistema AX = b para obtener su Onica soluci6n y los calculos son estables con respecto al crecimiento de los errores de redondeo V
N6tese que como consecuencia del teorema anterior se tiene que Si A = (a l ) )n xn es EDD
por filas entonces A tiene factorizaci6n LU es decir A == LU con L triangular inferior con unos en su diagonal principal y U triangular superior (escalonada)
sea tambien
Observe que la matriz de coeficientes del ejemplo 33 anterior no es EDD por filas
sectJeor9~~am9i~ -valioo=paramatrj simetricas
(vease Bllrden 19 5 ~ina 368) Una matriz A E Rnxn ~ i~etrica se dice definida
positiva si satisface una cualquiera de las sig uientes condiciones (Ias- cualesson equival~tes) a _-----
- - -- shy
i) ~AX gt 0 paralodo 2CERn X 0 J
iiXTodos los valores propios de A son pos1tlvOS7
iii) Todos los pivotes obtenidos en la eliminaci6n Gaussiana sobre A sin intercambio de filas - ~on positiv~s
iv) Todas las submatrices principales de A tienen determinante positiv~
(Las submatrices principales de la matriz A == (ai J) son las matrices - n ~ n
1k a ] k == 12bull f) )
akk v
114 METODOS NUMERICOS
N6tese nuevamente que Si A E R nxn es simetrica y definida positiva entonces A tiene
factorizaci6n A = LU con L triangular inferior con sus componentes sabre la diagonal principal iguales a uno y U triangular superior (ver mas adelante factorizaci6n de Choleski)
Observe que la matriz de coeficientes de l ejemplo 33 anterior no es simetrica
34 ESTRATEGIAS DE PIVOTEO
EI eJemplo 33 anterior muestra una de las dificultades que pueden surgir al aplicar el metodo
de eliminaci6n Gaussiana cuando el pivote ajj-1) es pequeno comparado con algunos
(j-1) e ementos I a t para J ~ I t ~ n
Para tratar de evitar tales dificultades se introduce en el metodo de eliminaci6n Gaussiana una estrategia IIamada de pivoteo la cual consiste en seleccionar el pivote de acuerdo con un cierto criter io Nosotros usaremos dos estrategias la estrategia de pivoteo maximo por columna 0 pivoteo parcial y la estrategia de pivoteo escalado de fila 0 escalamiento
341 Pivoteo maximo por columna 0 pivoteo parcial Esta estrategia difiere de
eliminaci6n Gaussiana simple unicamente en la escogencia del pivote ajj-1) la cual se hace
ahora as
Para j = 12 n - 1 se determina el menor entero k J ~ k ~ n tal que
y
es decir seleccionamos el primer elemento diferente de cero sobre la columna j-esima a
partir de la j-esima fila y que tenga mayor valor absoluto (para j = 1 a~j j1 ) = al~) ak 1 )
Si tal k no existe el sistema no tiene soluci6n unica y el proceso se puede terminar
Si tal k existe y k ot j entonces hacemos intercambio de las ecuaciones j-esima y k-esima
y continuamos con la eliminaci6n Gaussiana l
Ilustremos esta estrategia para resolver el sistema
E1 03X1 + 589x2 = 592 E2 531x1 - 61 OX 2 = 470
entonces k = 2 ot 1= j asl que
Por sustituci6n regresiva
Observe que en este caSO
SWAP(A i j)
que es el mismo del ejemplo 33 usando aritmetica con redondeo a tres digitos
Como para j = 1 se tiene que
Capitulo 3 SOLUCION NUMERICA DE SISTEMAS DE ECUACIONES 115
Max Iall I Ia21 I = Max 03 531 = 531 = Ia21 I 0
entonces k = 2 1= j as i que intercambiamos El y E2 Y continuamos con la el imi nac i6n
531 - 610 (A b)= ( 03 589 592) 470) (531 - 610 470 03 589 592
E2I (~) ( m21=5 65 ~ 1O - l) 531 - 610 5 31 470) ------------------~) ( o 589 589
Por sustituci6n regresiva obtenemos
x = 589 = 100 Xl = 470 +610(100) = 531 = 100 2 589 531 531
Observe que en este caso X= (J es la soluci6n exacta del sistema dado bull
Instrucci6n en DERIVE
SWAP(A i j) Intercambia las filas (0 elementos) i y j de la matriz A (de un vector) 0
Nota En el procedimiento de pivoteo maximo por columna (p ivoteo parcial) cada multiplicador mi j es tal que
y aunque esta estrategia permite resolver satisfactoriamente muchos sistemas de ecuaciones lineales hay casos donde fracasa como se ilustra en el siguiente ejemplo
Ejemplo 34 Consideremos el sistema
El 300Xl + 58900x2 59200 E2 531x1 -610X2 - 470
el cual es un sistema equivalente al del ejemplo 33 (los coeficientes de la primera ecuaci6n
en el sistema del ejemplo 33 han sido multiplicados por 10 3) EI pivoteo maximo por
columna con aritmetica de redondeo a tres dig itos nos lIeva a los siguientes resultados
59200) E2 1 (-~pound) (m 2 1 =177) ) ( 300 58900 59200 ) (A b) = ( 300 58900 531 - 610 470 o - 10400 -10500
y por sustituci6n regresiva x2 = 101 Y Xl = - 100 que es la misma solucion que se obtiene si ~ usamos eliminaci6n Gaussiana simple
116 METODOS NUMERICOS
En casas como el de este ejemplo donde un pivote es mucho mas pequeno que alguno de los coeficientes de la ecuacion que el encabeza se recomienda la tecnica conocida como pivoteo escalado de fila 0 escalamiento la cual es nuestra segunda estrategia
342 Pivoteo escalado de fila Esta tecnica solo difiere de la eliminaci6n Gaussiana simple al igual que el pivoteo parcial en la escogencia del pivote
lEsta vez el pivote alr ) se escoge como se indica a continuacion
Para j = 12 n - 1 hacemos 10 siguiente
a) Para i = j j + 1 n calculamos
S = Max Ia(l- l) I Factor de escala I j s i s n II
Si Sj = 0 entonces el sistema no tiene solucion unica y el proceso se puede terminar
b) Para i = jj+1 n calculamos
c) Encontramos el menor entero k con j o k 0 n tal que
Ia~jjl ) I IaH- l) I -------- = Max --------
Sk j S j $ n SI
Si tal k no existe entonces el sistema no tiene solucion unica y el proceso se puede terminar
Si tal k existe y kF j entonces hacemos intercambio de las ecuaciones j-esima y kshyesima
y continuamos con la eliminacion Gaussiana V
Apliquemos esta estrategia para resolver el sistema del ejemplo 34 usando aritmetica con redondeo a tres dfgitos
Para j = 1
a) Sl = Max Iall I I a12 I = Max 300 58900 = 58900F 0 Y
Max 53~ 610 IshyS2 =Max a21 a n -shy
b) Ahora ~= 51
la0
l 300
S2
531 al l I a21 -Max - 0 -
c) Max -S -S- - 5890 NIl 1 2
las ecuaciolllltanto intercamblamos
Gaussiana
300 58900 (A b) = ( 531 -610
Por sustitucion regresiva
Xl tOO
Observe que en este GaSO X
35 FACTORIZACION
Capitulo 3 SOLUCION NUMERICA DE SISTEMAS DE ECUACIONES 117
8 2 =Max I a21 I I a22 1 =Max 531 610 = 610 0
b) Ahora
Ia I 300 - 8- = 58900 Y
1
Ia 2 I 531 - - = - shy
6108 2
c) Max ~ ~ I=Max ~ 531 = 531 =~ 0 as que k = 2 1= j y por8 8 2 58900 60 60 8 2
tanto intercambiamos las ecuaciones E y E2 Y continuamos con la eliminacion
Gaussiana
(A b) = (300 58900 59200) p ( 531 -610 470 ) 2 )
531 - 610 470 300 58900 59200
E2 (300) (m2 =565 ) ( 531 -6 1 0531 47 0 )
--------------~) o 58900 58900
Por sustituci6n regresiva
_ 470 +610(100) 531X2 = 100 Xl = =-- = 100
531 531
Observe que en este caso X= ( ~ J = ( 100) es la solucion exacta del sistema bull 100x2
35 FACTORIZACION TRIANGULAR
Consideremos un sistema AX = b con A no-singular y b 0 Con respecto a la matriz A se sabe que existen matrices P de permutaci6n L triangular inferior con sus componentes sobre la diagonal principal iguales a uno y U triangular superior (escalon ada) tales que PA = LU Una forma eficiente computacionalmente de encontrar P L Y U usando eliminaci6n Gaussiana se muestra en el siguiente ejemplo Ejemplo 35 Resuelva el siguiente sistema de ecuaciones lineales usando eliminacion Gaussiana con pivoteo parcial y obtenga una factorizaci6n PA =LU para la matriz A de coeficientes asociada con este metodo
- x + 2X2 - 3X3 = -2 3x1 - 3x2 - X 3 = -4
x + x2 = 3IEmpezamos introduciendo un vector p =(PP2P3) T el cual inicializamos_con P = i i = 123 Y
donde se almacenaran los intercambios necesarios en el proceso de eliminacion Gaussiana
118 METODOS NUMERICOS
con pivoteo parcial (el numero de componentes del vector p coincide con el orden n del sistema a resolver)
- 1 2 -3 -2] -4] - 2 P = (123) T (Ab) = ~ - 3 -1 - 4
[ o 3 3
3 - 3 - 1 - 4(m=-H (m3l = ~ )
E( - ~) E3(~) 10 10(-i) 3 3
2(i) 1 13 3 3
(Observe que cada multiplicador mij es almacenado en la posici6n correspondiente (i j) en
la matriz de trabajo)
3 - 3 -1 - 4
Max 121 =2 p=( 23 1) T bull (i) 2 1
-
3 13 3
(-i) 10 3
10 3
(Observe que la permutaci6n se hace para las filas 2 y 3 completas es decir incluyendo los multiplicadores)
3 - 3 - 1 --4(mJ2 = ~ ) E3( ~)
2(i) 13 3 3
(- i) (i) 7 11 2 2
La eficiencia en el metoda indicado se debe a que en la misma matriz de trabajo se almacenan los multiplicadores que van a conformar la matriz L (en el ejemplo son los numeros que se encuentran dentro de parentesis) 10 que significa un ahorro de memoria y como los intercambios necesarios afectan simultaneamente a las matrices L y U se evita tener que volver a la matriz original a realizar los intercambios ya observados y repetir la eliminaci6n Gaussiana con pivoteo parcial De esta manera al terminar el proceso de eliminaci6n podemos leer en la matriz final la parte estrictamente triangular inferior de L (son los numeros entre parentesis) y la matriz triangular superior U (que es la parte triangular superior de la matriz final) y en el vector p final quedan almacenados los intercambios realizados que se usan para producir la matriz de permutaci6n P
Para el ejemplo 35
3 -31 0 0 1 0 2U=0L= 3 1 1 0 0-3 2
(Verifique que PA =LU)
Para obtener la soluci6n del sistema
reducido
yobtenemos
C6mo se resuelve el
Capitulo 3 SOLUCION NUMERICA DE SISTEMAS DE ECUACIONES 119
Para el ejemplo 35
1 0 0 3 -3 - 1 1 1 1
L= - 0 U= 0 2 - Y como P~ (231) entonces [~ 0 lp3 3 1 1 7 0
-- - 0 0 3 2 2
(Verifique que PA = LU )
Para obtener la soluci6n del sistema original usamos sustituci6n regresiva en el sistema reducido
-4 13
UX = 3
11
2 ~
c y obtenemos
asi que la soluci6n (exacta) del sistema dado es X = (23 40 ~)T21 21 7
C6mo se resuelve el sistema AX = b a partir de la factorizaci6n PA = LU obtenida bull
351 Algunas aplicaciones de la factorizaci6n PA = LU La factorizaci6n PA = LU es utilizada eficientemente en aquellos casos donde se trabaja repetfdamente con la misma matriz A Dos de esos casos se presentan a continuaci6n
1) Resolver varios sistemas AX =b con la misma matriz de coeficientes A ya que en p L Y U esta almacenado todo el proceso de eliminaci6n Gaussiana EI algoritmo se bas a en la siguiente equivalencia
AX = b lt=gt PAX = Pb lt=gt LUX = Pb
lt=gt UX = L- 1Pb
lt=gtUXoo c y c oo L- 1Pb
lt=gt UX == C Y Lc = Pb lt=gt Lc = Pb Y UX = c
Los pasos a seguir son
Paso 1 Calcular Pb Paso 2 Resolver para c Lc = Pb por sustituci6n progresiva Paso 3 Resolver para X UX = c por sustituci6n regresiva
Como ejercicio resuelva el sistema del ejemplo anterior usando este algoritmo
120 METODOS NUMERICOS
2) Encontrar la matriz inversa A -1 de una matriz invertible A resolviendo los n-sistemas
AX =eU) j = 12 n
donde e(1) =(O O1 o 0f ERn
t posicion j
La solucion X del sistema AX = e(n j = 12 n produce la correspondiente columna j shy
esima de la matriz A -1
Como ejercicio compare el numero de operaciones para encontrar A - usando el metodo de Gauss-Jordan con el numero de operaciones resolviendo los n-sistemas indicados antes V
Ejercicio 33 Calcule la inversa de la matriz A de coeficientes del sistema del ejemplo 35 usando el metodo de Gauss-Jordan y tambien usando la factorizacion PA = LU bull
36 SISTEMAS TRIDIAGONALES
Un caso muy importante de sistemas de ecuaciones lineales que requiere un tratamiento especial es el de los sistemas tridiagonales Tales sistemas aparecen en diversas aplicaciones como por ejemplo al utilizar metodos de diferencias flnitas en la solucion de problemas con valores en la frontera para ecuaciones diferenciales ordinarias y como veremos mas adelante en el problema de la interpolacion segmentaria cubica
Un sistema tridiagonal es de la forma
dx1 + C1X2
a 2 x + d2 x2 + C 2 X 3
a3 x2 + d3 x3 + C 3 X 4
an-1Xn- 2 + dn_ xn_ + cn_xn = bn_1
anxn_1 + dnxn = bn
La matriz de coeficientes del sistema es
C 0 0 od1
0a 2 d2 c2
o a3 d3 c 3
A =
o
la cual se dice una matriz TRIDIAGONAL
En general
ij=12 n
suponiendo que la matriz A de podemos usar eliminaclon Gausslana eliminacion Gaussiana simple es adecuadl EDD por filas) Otra forma de resolver
EDD por filas es a partir de la componentes de las matrices L Y U as triangular inferior con todos sus superior Para encontrar tales malnce5
o 0 0
Y2 o o o o
L =
o Yn t
o o
Como Cl Ct
Y2 Cl Y2C +az
0 h a2
LU =
0
0
d
=
(37)
(38)
(39)
(310)
Capitulo 3 SOLUCION NUMERICA DE SISTEMAS DE ECUACIONES 107
distinguir una soluci6n aproximada X buena de una mala observando el error residual
relativo ~ b
EI numero Cond (A) = II A 1111 A -1 II se IIamara NUMERO DE CONDICION 0 CONDICIONAL
de la matriz no-singular A relativo a la norma matricial usada Aunque el valor de Cond(A)
depende de la norma matricial usada sin embargo Cond (A) ~ 1 cualquiera sea la norma
matricial inducida pues
In = AA -1 In II ~ II A 1111 A - 1 II y In II =~~o II In II ~ II =~~d ~ II =1
De acuerdo con la relaci6n (3 7) dada en el teorema 33 vemos que si Cond (A) 1
entonces el error relativo II ~ II y el error residual relativo II ii son mas 0 menos del
mismo tamano y podremos distinguir una soluci6n aproximada buena de una mala
observando el error residual relativo pero entre mas grande sea Cond(A) menor es la
informaci6n que se puede obtener del error relativo a partir del error residual relativo
De 10 anterior se espera que A tenga un buen comportamiento en el sentido de que un error residual relativo pequeno implique correspondientemente una buena soluci6n aproximada
de AX =b si Cond (A) 1 caso en el cual diremos que A esta bien condicionada (el
sistema AX =b esta bien condicionado) Si Cond (A) gtgt 1 es posible que A tenga un mal
comportamiento en el sentido que un error residual relativo pequeno puede corresponder a una soluci6n aproximada mala y diremos que A esta mal condicionada (el sistema AX = b esta mal condicionado)
A pesar de las definiciones anteriores no debemos olvidar que 10 que realmente nos interesa
es poder determinar cuando una soluci6n aproximada Xde un sistema AX = b es buena y tratar de distinguir si el sistema AX = b esta bien 0 mal condicionado
Para la matriz
del ejemplo 31 tenemos
A- 1 _ 1 ( 10 -1J - 05 - 1005 1
I A =Max 111 + 111 11005 1 + 1 10 1 = Max 2 2005 = 2005
II A-1 II = _1 II ( 10 - 1J II =_1 1105 = 221 05 - 1005 1 05
co
luego
108 METODOS NUMERICOS
Cond - (A) == II A 1100 II A - 110) = (2005)(221) == 443105 raquo 1
pequeno puedeEste numero de cond icion nos dice que un error residual relativo
muy grande asi que A puede considerarse mal
cond icionada
Veamos que puede decirse en este caso de la calidad de la soluci6n aproximada
- ( 10)
corresponder a un error relativo
X == - 8 del sistema
X +y == 2 1005x+10y==21
Para este ejemplo tenemos
II RII 5 y Cond (A) == 443105
II b II en 21
asi que la desigualdad (37) dada en el teorema 33 se convierte en
esto es
_ Ilx - x1 11 537 x10 6 ~ II x 110) 0) ~ 1055
5 10 que indica que aunque el error residual relativo es pequeno 21 el numero de condicion
es tan grande (443105) que hace que la solucion calculada pueda tener un error relativo de
hasta 1055 asi que nada puede decirse de la cercania entre X y X bull
Instrucciones en DERIVE
NORMAJNF(A) Simplifica en la norma del maximo de la matriz A II A II
COND_INF(A) Simplifica en el numero de condici6n relativo a la norma del maximo de la
matriz A es decir simplifica en el numero Cond O) (A) == II f-II00 11 A - 1 11 00 0
Existe otro numero asociado con una matriz al cual se Ie denomina tambien numero de condicion A continuaci6n nos referiremos a tal numero
Del teorema 32 se sabe que p (A)sl AI
I
pero como los valores propios de A tiene que
con cr(A) == A ECAes 1
EI numero Cond(A
Para la matriz A ( 1~
numero de
Capitulo 3 SOLUCION NUMERICA DE SISTEMAS DE ECUACIONES 109
Del teorema 32 se sabe que p (A) II A II para toda norma matricial inducida asl que
Cond (A) = II A IIII A -1 ~p(A) p(A-1)
pero como los valores propios de A - 1 son los reciprocos de los valores propios de A se tiene que
Max I A I ) Ecr(A) ( )
Cond (A) ~ I I == Condo AMm A
AEcr( A)
con a(A) =AE C A esvalorpropiodeA espectro de A (Recuerde que
p(A -1) = Max I A I = 1 I ) ecr(A-1) Mm A
AEcr( A)
Max I A I EI numero Cond (A)= AEcr(A) se denomina numero de condici6n espectral de A
bull Min I A I AEcr( A)
Segun se acaba de probar Cond(A) ~ Cond(A)
(1 1J Para la matnz A = se bene que 1005 10
det(A_AI) = 1- A 1 I =A2- 11A- 05 1005 10 - A
as que los valores propios de A son A1 1100454358 A 2 -45435778 x 10 -3 Y por tanto
Condo (A) 1100454358 2421999592 raquo 1 bull 45435778 x 10-3
Dado un sistema AX = b si 8A Y 8b denotan perturbaciones en A y b respectivamente el siguiente teorema cuya demostracion puede ser consultada en Ortega 1990 paginas 32 y
33 establece una cota para el error relativo II X- X II en terminos de las perturbaciones
IIXII
11 MII 118b II relatlvas lfAl lfbf y Cond(A) donde X es la soluclon exacta de AX = b Y X es la
solucion exacta del sistema perturbado (A + 8 A)X = b + 8 b
110 METODOS NUMERICOS
Teorema 34 Supongase que A es no-singular y que 11 8 A II lt II A~ II (esta hipotes is asegura
que A + 8 A es invertible y que 1- Cond(A) li SA II gt 0) Si X es la solucion exacta del IIA II
sistema perturbado (A + 0 A)X = b + 8 b entonces X aproxima a la solucion exacta X del
sistema AX = b b cI 0 con la siguiente estimacion de error
II X - X II ~ Cond (A) ( IL~~JL~ J (311)
IIXII 1_cond(A)[ 13 All l llb ll IIAII
II AII
La desigualdad (3 11) dice que si la matriz A esta bien condicionada es decir si
Cond (A) 1 entonces cambios pequerios en A y b producen correspondientemente
cambios pequerios en la soluci6n del sistema (el sistema AX = b esta bien condicionado) Por otro lado si A esta mal condicionada entonces cam bios pequenos en A y b pueden producir grandes cam bios en la solucion del sistema (el sistema AX = b esta mal condicionado)
Ejercicio 31 Estime la cota de error dada en el teorema 34 para los sistemas (34) y (34) del ejemplo 31 bull
Ejercicio 32 a) Calcule Cond (A) usando 11middot11 2 1111 y II to para las siguientes matrices
456 218) ( 279 138
b) Que puede decir del condicionamiento de los siguientes sistemas de ecuaciones lineales
f39x + 16x2 = 55 456X l + 218x2 = 674 I) ii) bulll 68x + 29x2 = 97 279x + 138x 2 = 417
ESTABILIDAD NUMERICA EN LA ELiMINACION GAUSSIANA
Volvamos al metodo de eliminacion Gaussiana (simple) y consideremos el siguiente ejemplo
Ejemplo 33 Resuelva el siguiente sistema de ecuaciones lineales usando eliminacion Gaussiana con sustitucion regresiva y aritmetica (decimal) con redondeo a tres dfgitos
El 03x + 589x2 = 592 E2 531x - 610x2 = 470
Usando elimnaci6n Gaussiana obtenemos
03 589 592) ~~- (A b) = ( 531 - 610 470
y por sustituci6n regresiva
_ 592 - 589(101 592 Xl = 03
luego la solucion calculada es X (~2 1 =(shy
Instrucci6n en DERIVE
PIVOT(A i j) Usa oper~cione~e tiene matriz obtenida de la matnz A q
Que puede decir de la calidad de la
Para intentar responder esta
dadas par el teorema 3 3
03 589) Como A = ( 531 -610
los calculos se obtienen III
v
(311 )
Capitulo 3 SOlUCION NUMERICA DE SISTEMAS DE ECUACIONES 111
(A b) = ( 03 589 592) __( 5__) _m 2 _=1_77 ) ~) ( 03 589E2_ ~~ _( _ _ 592 ) 531 -610 470 0 - 10400 -10500
y por sustituci6n regresiva
x = -10500 =101 2 - 10400
)(1 = 592 - 589(101) = 592 - 595 = - 3 = - 100 D3 D3 D3
- (X1 J (-1 00)luego la soluci6n calculada es X = _ = 101x2
Instrucci6n en DERIVE
PIVOT(A i j) Usa operaciones elemetales de fila para Simplificar (0 aproXimar) en una matriz obtenida de la matriz A que tiene ceros en la columna j y por debajo de la fila i 0
Que puede decir de la calidad de la soluci6n aproximada X
IIX -XIIPara intentar responder esta pregunta encontremos las cotas para el error relativo
IIXII dadas por el teorema 33
03 589)Como A = entonces usando aritmetica con redondeo a tres digitos para todos( 531 - 610
los calculos se obtienen las siguientes aproximaciones
- 1 1 (-610 - 589)A =-shy-313 - 531 03
II A II = Max 589 114 = 589
A - 1 II = _ 1_ Max 650 534 = _ 1_ 650 = 208II 313 313
entonces Cond (A) = II A II II A -111 = (589)(208) = 123 que no es muy grande
comparado con uno asi que la matriz A puede considerarse bien condicionada
(Por ciertas consideraciones te6ricas sobre el numero de condici6n de una matriz A las cuales pueden ser consultadas en Burden 1985 paginas 481 y 482 cuando se trabaja en
aritmetica finita (decimal) con redondeo a t-digitos y Cond(A) ~ 10t se espera un mal
112 METOOOS NUMERICOS
comportamiento de A con respecto a la solucion de AX = b Y A se considera mal
condicionada En este ejemplo Cond(A) = 123 lt 103 )
Ahora
03 589 ) (-100) ( 592)AX-b e ( 531 - 610 101 470
Estas operaciones se realizari en doble preci sion (6 digitos)
59189) ( 592) ( - 011 )(= -59261 - 470 = -106261
(Para evitar la perdida de cifras sign ificativas se debe calcular el vector erro r res idual
R = AX - b en doble precision)
Conv irtiendo este ultimo resultado a tres digltos usando redondeo se obtiene
R = (-011) - 106
Entonces II R II =106 y II b II = 592 Y por tanto
II R la 1 106 1 11gtlt - X II 106 II R II en
-II-11- ( )= - - =145 ~ II II ~ 2202= 123- = Cond (A)-II - lI shyb co Cond A 592 123 X 592 b DO
II R pero como Cond (A) ~ 1 entonces se espera que fbiC Y X- X sean mas 0 menosII xL
II X- X L sea tamble ndel mismo tamano y ya que es grande se espera que II X II
grande +
En situaciones como a observada en este ejemplo se sugiere hacer un refinamiento
iterativo sobre la soluci6n caculada X 0 usar esta soluci6n calcuada como aproximaci6n inicial en un metodo iterativ~ con et proposito de tratar de mejorar la soluci6n aproximada y ograr que el error residua l relativo sea mas pequeno EI metodo de reflna81iento iterativo puede ser consutado en Kincaid 1972 pag inas 174-176
La soucion exacta del sistema en consideracion es X = (J= ( 0) A que se debe a
diferencia entre la soucion exacta y a solucion calculada
Observe que e error en e calculo de )(2 con respecto a x2 fue de solo 01 (un error relativo
de 1) Y este error fue multiplicado par un factor de aproximadamente - 2000 al obtener
)(1 debido a orden en que se realiz6 la el iminaci6n Gaussiana
Instruccion en DERIVE
RESUELVA_1(Ab) Simplifica en a SOUCI6n entra como un vector fila 0
En el ejemplo anterior la elimina~i6n sistema de ecuaciones lineales ble~ nII1II
del algoritmo de eliminaci6n Gausslana pequeno) Hay sin embargo SIlUC -
Gaussiana es numencamente estable consultarse en Burden 1985 paginas366 y
Teorema 35 Si A =(ajiLn (EDD) par filas es decir si
entonces A es invertible (no-singular intercambio de filas en cualquier calculos son estables con respeclO
Notese que como consecuencia
por filas entonces A t~ene unos en su diagonal pnnClpalyU
~~a~ (vease Burden 19a5
satisface
Capitulo 3 SOLUCION NUMERIC A DE SISTEMAS DE ECUACIONES 113
Instrucci6n en DERIVE
RESUELVA_1(Ab) Simplifica en la soluci6n exacta X del sistema AX=b EI vector b se entra como un vector fila 0
En el ejemplo anterior la eliminaci6n Gaussiana condujo a una respuesta defectuosa de un sistema de ecuaciones lineales bien condicionado Esto muestra la inestabilidad numerica del algoritmo de eliminaci6n Gaussiana (consecuencia de la divisi6n por un nOmero (pivote) pequero) Hay sin embargo situaciones en las cuales el algoritmo de eliminaci6n Gaussiana es numericamente estable EI siguiente teorema cuya demostraci6n puede consultarse en Burden 1985 paginas 366 y 367 se refiere a una de tales situaciones
Teorema 35 Si A = (ai j )n xn es una matriz estrictamente dominante diagonalmente
(EDD) por filas es decir si
n Iaii Igt II ai j I para cada i == 12 n
)= 1 ) 1
entonces A es invertible (no-singular) Ademas se puede rea lizar eliminaci6n Gaussiana sin intercambio de filas en cualquier sistema AX = b para obtener su Onica soluci6n y los calculos son estables con respecto al crecimiento de los errores de redondeo V
N6tese que como consecuencia del teorema anterior se tiene que Si A = (a l ) )n xn es EDD
por filas entonces A tiene factorizaci6n LU es decir A == LU con L triangular inferior con unos en su diagonal principal y U triangular superior (escalonada)
sea tambien
Observe que la matriz de coeficientes del ejemplo 33 anterior no es EDD por filas
sectJeor9~~am9i~ -valioo=paramatrj simetricas
(vease Bllrden 19 5 ~ina 368) Una matriz A E Rnxn ~ i~etrica se dice definida
positiva si satisface una cualquiera de las sig uientes condiciones (Ias- cualesson equival~tes) a _-----
- - -- shy
i) ~AX gt 0 paralodo 2CERn X 0 J
iiXTodos los valores propios de A son pos1tlvOS7
iii) Todos los pivotes obtenidos en la eliminaci6n Gaussiana sobre A sin intercambio de filas - ~on positiv~s
iv) Todas las submatrices principales de A tienen determinante positiv~
(Las submatrices principales de la matriz A == (ai J) son las matrices - n ~ n
1k a ] k == 12bull f) )
akk v
114 METODOS NUMERICOS
N6tese nuevamente que Si A E R nxn es simetrica y definida positiva entonces A tiene
factorizaci6n A = LU con L triangular inferior con sus componentes sabre la diagonal principal iguales a uno y U triangular superior (ver mas adelante factorizaci6n de Choleski)
Observe que la matriz de coeficientes de l ejemplo 33 anterior no es simetrica
34 ESTRATEGIAS DE PIVOTEO
EI eJemplo 33 anterior muestra una de las dificultades que pueden surgir al aplicar el metodo
de eliminaci6n Gaussiana cuando el pivote ajj-1) es pequeno comparado con algunos
(j-1) e ementos I a t para J ~ I t ~ n
Para tratar de evitar tales dificultades se introduce en el metodo de eliminaci6n Gaussiana una estrategia IIamada de pivoteo la cual consiste en seleccionar el pivote de acuerdo con un cierto criter io Nosotros usaremos dos estrategias la estrategia de pivoteo maximo por columna 0 pivoteo parcial y la estrategia de pivoteo escalado de fila 0 escalamiento
341 Pivoteo maximo por columna 0 pivoteo parcial Esta estrategia difiere de
eliminaci6n Gaussiana simple unicamente en la escogencia del pivote ajj-1) la cual se hace
ahora as
Para j = 12 n - 1 se determina el menor entero k J ~ k ~ n tal que
y
es decir seleccionamos el primer elemento diferente de cero sobre la columna j-esima a
partir de la j-esima fila y que tenga mayor valor absoluto (para j = 1 a~j j1 ) = al~) ak 1 )
Si tal k no existe el sistema no tiene soluci6n unica y el proceso se puede terminar
Si tal k existe y k ot j entonces hacemos intercambio de las ecuaciones j-esima y k-esima
y continuamos con la eliminaci6n Gaussiana l
Ilustremos esta estrategia para resolver el sistema
E1 03X1 + 589x2 = 592 E2 531x1 - 61 OX 2 = 470
entonces k = 2 ot 1= j asl que
Por sustituci6n regresiva
Observe que en este caSO
SWAP(A i j)
que es el mismo del ejemplo 33 usando aritmetica con redondeo a tres digitos
Como para j = 1 se tiene que
Capitulo 3 SOLUCION NUMERICA DE SISTEMAS DE ECUACIONES 115
Max Iall I Ia21 I = Max 03 531 = 531 = Ia21 I 0
entonces k = 2 1= j as i que intercambiamos El y E2 Y continuamos con la el imi nac i6n
531 - 610 (A b)= ( 03 589 592) 470) (531 - 610 470 03 589 592
E2I (~) ( m21=5 65 ~ 1O - l) 531 - 610 5 31 470) ------------------~) ( o 589 589
Por sustituci6n regresiva obtenemos
x = 589 = 100 Xl = 470 +610(100) = 531 = 100 2 589 531 531
Observe que en este caso X= (J es la soluci6n exacta del sistema dado bull
Instrucci6n en DERIVE
SWAP(A i j) Intercambia las filas (0 elementos) i y j de la matriz A (de un vector) 0
Nota En el procedimiento de pivoteo maximo por columna (p ivoteo parcial) cada multiplicador mi j es tal que
y aunque esta estrategia permite resolver satisfactoriamente muchos sistemas de ecuaciones lineales hay casos donde fracasa como se ilustra en el siguiente ejemplo
Ejemplo 34 Consideremos el sistema
El 300Xl + 58900x2 59200 E2 531x1 -610X2 - 470
el cual es un sistema equivalente al del ejemplo 33 (los coeficientes de la primera ecuaci6n
en el sistema del ejemplo 33 han sido multiplicados por 10 3) EI pivoteo maximo por
columna con aritmetica de redondeo a tres dig itos nos lIeva a los siguientes resultados
59200) E2 1 (-~pound) (m 2 1 =177) ) ( 300 58900 59200 ) (A b) = ( 300 58900 531 - 610 470 o - 10400 -10500
y por sustituci6n regresiva x2 = 101 Y Xl = - 100 que es la misma solucion que se obtiene si ~ usamos eliminaci6n Gaussiana simple
116 METODOS NUMERICOS
En casas como el de este ejemplo donde un pivote es mucho mas pequeno que alguno de los coeficientes de la ecuacion que el encabeza se recomienda la tecnica conocida como pivoteo escalado de fila 0 escalamiento la cual es nuestra segunda estrategia
342 Pivoteo escalado de fila Esta tecnica solo difiere de la eliminaci6n Gaussiana simple al igual que el pivoteo parcial en la escogencia del pivote
lEsta vez el pivote alr ) se escoge como se indica a continuacion
Para j = 12 n - 1 hacemos 10 siguiente
a) Para i = j j + 1 n calculamos
S = Max Ia(l- l) I Factor de escala I j s i s n II
Si Sj = 0 entonces el sistema no tiene solucion unica y el proceso se puede terminar
b) Para i = jj+1 n calculamos
c) Encontramos el menor entero k con j o k 0 n tal que
Ia~jjl ) I IaH- l) I -------- = Max --------
Sk j S j $ n SI
Si tal k no existe entonces el sistema no tiene solucion unica y el proceso se puede terminar
Si tal k existe y kF j entonces hacemos intercambio de las ecuaciones j-esima y kshyesima
y continuamos con la eliminacion Gaussiana V
Apliquemos esta estrategia para resolver el sistema del ejemplo 34 usando aritmetica con redondeo a tres dfgitos
Para j = 1
a) Sl = Max Iall I I a12 I = Max 300 58900 = 58900F 0 Y
Max 53~ 610 IshyS2 =Max a21 a n -shy
b) Ahora ~= 51
la0
l 300
S2
531 al l I a21 -Max - 0 -
c) Max -S -S- - 5890 NIl 1 2
las ecuaciolllltanto intercamblamos
Gaussiana
300 58900 (A b) = ( 531 -610
Por sustitucion regresiva
Xl tOO
Observe que en este GaSO X
35 FACTORIZACION
Capitulo 3 SOLUCION NUMERICA DE SISTEMAS DE ECUACIONES 117
8 2 =Max I a21 I I a22 1 =Max 531 610 = 610 0
b) Ahora
Ia I 300 - 8- = 58900 Y
1
Ia 2 I 531 - - = - shy
6108 2
c) Max ~ ~ I=Max ~ 531 = 531 =~ 0 as que k = 2 1= j y por8 8 2 58900 60 60 8 2
tanto intercambiamos las ecuaciones E y E2 Y continuamos con la eliminacion
Gaussiana
(A b) = (300 58900 59200) p ( 531 -610 470 ) 2 )
531 - 610 470 300 58900 59200
E2 (300) (m2 =565 ) ( 531 -6 1 0531 47 0 )
--------------~) o 58900 58900
Por sustituci6n regresiva
_ 470 +610(100) 531X2 = 100 Xl = =-- = 100
531 531
Observe que en este caso X= ( ~ J = ( 100) es la solucion exacta del sistema bull 100x2
35 FACTORIZACION TRIANGULAR
Consideremos un sistema AX = b con A no-singular y b 0 Con respecto a la matriz A se sabe que existen matrices P de permutaci6n L triangular inferior con sus componentes sobre la diagonal principal iguales a uno y U triangular superior (escalon ada) tales que PA = LU Una forma eficiente computacionalmente de encontrar P L Y U usando eliminaci6n Gaussiana se muestra en el siguiente ejemplo Ejemplo 35 Resuelva el siguiente sistema de ecuaciones lineales usando eliminacion Gaussiana con pivoteo parcial y obtenga una factorizaci6n PA =LU para la matriz A de coeficientes asociada con este metodo
- x + 2X2 - 3X3 = -2 3x1 - 3x2 - X 3 = -4
x + x2 = 3IEmpezamos introduciendo un vector p =(PP2P3) T el cual inicializamos_con P = i i = 123 Y
donde se almacenaran los intercambios necesarios en el proceso de eliminacion Gaussiana
118 METODOS NUMERICOS
con pivoteo parcial (el numero de componentes del vector p coincide con el orden n del sistema a resolver)
- 1 2 -3 -2] -4] - 2 P = (123) T (Ab) = ~ - 3 -1 - 4
[ o 3 3
3 - 3 - 1 - 4(m=-H (m3l = ~ )
E( - ~) E3(~) 10 10(-i) 3 3
2(i) 1 13 3 3
(Observe que cada multiplicador mij es almacenado en la posici6n correspondiente (i j) en
la matriz de trabajo)
3 - 3 -1 - 4
Max 121 =2 p=( 23 1) T bull (i) 2 1
-
3 13 3
(-i) 10 3
10 3
(Observe que la permutaci6n se hace para las filas 2 y 3 completas es decir incluyendo los multiplicadores)
3 - 3 - 1 --4(mJ2 = ~ ) E3( ~)
2(i) 13 3 3
(- i) (i) 7 11 2 2
La eficiencia en el metoda indicado se debe a que en la misma matriz de trabajo se almacenan los multiplicadores que van a conformar la matriz L (en el ejemplo son los numeros que se encuentran dentro de parentesis) 10 que significa un ahorro de memoria y como los intercambios necesarios afectan simultaneamente a las matrices L y U se evita tener que volver a la matriz original a realizar los intercambios ya observados y repetir la eliminaci6n Gaussiana con pivoteo parcial De esta manera al terminar el proceso de eliminaci6n podemos leer en la matriz final la parte estrictamente triangular inferior de L (son los numeros entre parentesis) y la matriz triangular superior U (que es la parte triangular superior de la matriz final) y en el vector p final quedan almacenados los intercambios realizados que se usan para producir la matriz de permutaci6n P
Para el ejemplo 35
3 -31 0 0 1 0 2U=0L= 3 1 1 0 0-3 2
(Verifique que PA =LU)
Para obtener la soluci6n del sistema
reducido
yobtenemos
C6mo se resuelve el
Capitulo 3 SOLUCION NUMERICA DE SISTEMAS DE ECUACIONES 119
Para el ejemplo 35
1 0 0 3 -3 - 1 1 1 1
L= - 0 U= 0 2 - Y como P~ (231) entonces [~ 0 lp3 3 1 1 7 0
-- - 0 0 3 2 2
(Verifique que PA = LU )
Para obtener la soluci6n del sistema original usamos sustituci6n regresiva en el sistema reducido
-4 13
UX = 3
11
2 ~
c y obtenemos
asi que la soluci6n (exacta) del sistema dado es X = (23 40 ~)T21 21 7
C6mo se resuelve el sistema AX = b a partir de la factorizaci6n PA = LU obtenida bull
351 Algunas aplicaciones de la factorizaci6n PA = LU La factorizaci6n PA = LU es utilizada eficientemente en aquellos casos donde se trabaja repetfdamente con la misma matriz A Dos de esos casos se presentan a continuaci6n
1) Resolver varios sistemas AX =b con la misma matriz de coeficientes A ya que en p L Y U esta almacenado todo el proceso de eliminaci6n Gaussiana EI algoritmo se bas a en la siguiente equivalencia
AX = b lt=gt PAX = Pb lt=gt LUX = Pb
lt=gt UX = L- 1Pb
lt=gtUXoo c y c oo L- 1Pb
lt=gt UX == C Y Lc = Pb lt=gt Lc = Pb Y UX = c
Los pasos a seguir son
Paso 1 Calcular Pb Paso 2 Resolver para c Lc = Pb por sustituci6n progresiva Paso 3 Resolver para X UX = c por sustituci6n regresiva
Como ejercicio resuelva el sistema del ejemplo anterior usando este algoritmo
120 METODOS NUMERICOS
2) Encontrar la matriz inversa A -1 de una matriz invertible A resolviendo los n-sistemas
AX =eU) j = 12 n
donde e(1) =(O O1 o 0f ERn
t posicion j
La solucion X del sistema AX = e(n j = 12 n produce la correspondiente columna j shy
esima de la matriz A -1
Como ejercicio compare el numero de operaciones para encontrar A - usando el metodo de Gauss-Jordan con el numero de operaciones resolviendo los n-sistemas indicados antes V
Ejercicio 33 Calcule la inversa de la matriz A de coeficientes del sistema del ejemplo 35 usando el metodo de Gauss-Jordan y tambien usando la factorizacion PA = LU bull
36 SISTEMAS TRIDIAGONALES
Un caso muy importante de sistemas de ecuaciones lineales que requiere un tratamiento especial es el de los sistemas tridiagonales Tales sistemas aparecen en diversas aplicaciones como por ejemplo al utilizar metodos de diferencias flnitas en la solucion de problemas con valores en la frontera para ecuaciones diferenciales ordinarias y como veremos mas adelante en el problema de la interpolacion segmentaria cubica
Un sistema tridiagonal es de la forma
dx1 + C1X2
a 2 x + d2 x2 + C 2 X 3
a3 x2 + d3 x3 + C 3 X 4
an-1Xn- 2 + dn_ xn_ + cn_xn = bn_1
anxn_1 + dnxn = bn
La matriz de coeficientes del sistema es
C 0 0 od1
0a 2 d2 c2
o a3 d3 c 3
A =
o
la cual se dice una matriz TRIDIAGONAL
En general
ij=12 n
suponiendo que la matriz A de podemos usar eliminaclon Gausslana eliminacion Gaussiana simple es adecuadl EDD por filas) Otra forma de resolver
EDD por filas es a partir de la componentes de las matrices L Y U as triangular inferior con todos sus superior Para encontrar tales malnce5
o 0 0
Y2 o o o o
L =
o Yn t
o o
Como Cl Ct
Y2 Cl Y2C +az
0 h a2
LU =
0
0
d
=
108 METODOS NUMERICOS
Cond - (A) == II A 1100 II A - 110) = (2005)(221) == 443105 raquo 1
pequeno puedeEste numero de cond icion nos dice que un error residual relativo
muy grande asi que A puede considerarse mal
cond icionada
Veamos que puede decirse en este caso de la calidad de la soluci6n aproximada
- ( 10)
corresponder a un error relativo
X == - 8 del sistema
X +y == 2 1005x+10y==21
Para este ejemplo tenemos
II RII 5 y Cond (A) == 443105
II b II en 21
asi que la desigualdad (37) dada en el teorema 33 se convierte en
esto es
_ Ilx - x1 11 537 x10 6 ~ II x 110) 0) ~ 1055
5 10 que indica que aunque el error residual relativo es pequeno 21 el numero de condicion
es tan grande (443105) que hace que la solucion calculada pueda tener un error relativo de
hasta 1055 asi que nada puede decirse de la cercania entre X y X bull
Instrucciones en DERIVE
NORMAJNF(A) Simplifica en la norma del maximo de la matriz A II A II
COND_INF(A) Simplifica en el numero de condici6n relativo a la norma del maximo de la
matriz A es decir simplifica en el numero Cond O) (A) == II f-II00 11 A - 1 11 00 0
Existe otro numero asociado con una matriz al cual se Ie denomina tambien numero de condicion A continuaci6n nos referiremos a tal numero
Del teorema 32 se sabe que p (A)sl AI
I
pero como los valores propios de A tiene que
con cr(A) == A ECAes 1
EI numero Cond(A
Para la matriz A ( 1~
numero de
Capitulo 3 SOLUCION NUMERICA DE SISTEMAS DE ECUACIONES 109
Del teorema 32 se sabe que p (A) II A II para toda norma matricial inducida asl que
Cond (A) = II A IIII A -1 ~p(A) p(A-1)
pero como los valores propios de A - 1 son los reciprocos de los valores propios de A se tiene que
Max I A I ) Ecr(A) ( )
Cond (A) ~ I I == Condo AMm A
AEcr( A)
con a(A) =AE C A esvalorpropiodeA espectro de A (Recuerde que
p(A -1) = Max I A I = 1 I ) ecr(A-1) Mm A
AEcr( A)
Max I A I EI numero Cond (A)= AEcr(A) se denomina numero de condici6n espectral de A
bull Min I A I AEcr( A)
Segun se acaba de probar Cond(A) ~ Cond(A)
(1 1J Para la matnz A = se bene que 1005 10
det(A_AI) = 1- A 1 I =A2- 11A- 05 1005 10 - A
as que los valores propios de A son A1 1100454358 A 2 -45435778 x 10 -3 Y por tanto
Condo (A) 1100454358 2421999592 raquo 1 bull 45435778 x 10-3
Dado un sistema AX = b si 8A Y 8b denotan perturbaciones en A y b respectivamente el siguiente teorema cuya demostracion puede ser consultada en Ortega 1990 paginas 32 y
33 establece una cota para el error relativo II X- X II en terminos de las perturbaciones
IIXII
11 MII 118b II relatlvas lfAl lfbf y Cond(A) donde X es la soluclon exacta de AX = b Y X es la
solucion exacta del sistema perturbado (A + 8 A)X = b + 8 b
110 METODOS NUMERICOS
Teorema 34 Supongase que A es no-singular y que 11 8 A II lt II A~ II (esta hipotes is asegura
que A + 8 A es invertible y que 1- Cond(A) li SA II gt 0) Si X es la solucion exacta del IIA II
sistema perturbado (A + 0 A)X = b + 8 b entonces X aproxima a la solucion exacta X del
sistema AX = b b cI 0 con la siguiente estimacion de error
II X - X II ~ Cond (A) ( IL~~JL~ J (311)
IIXII 1_cond(A)[ 13 All l llb ll IIAII
II AII
La desigualdad (3 11) dice que si la matriz A esta bien condicionada es decir si
Cond (A) 1 entonces cambios pequerios en A y b producen correspondientemente
cambios pequerios en la soluci6n del sistema (el sistema AX = b esta bien condicionado) Por otro lado si A esta mal condicionada entonces cam bios pequenos en A y b pueden producir grandes cam bios en la solucion del sistema (el sistema AX = b esta mal condicionado)
Ejercicio 31 Estime la cota de error dada en el teorema 34 para los sistemas (34) y (34) del ejemplo 31 bull
Ejercicio 32 a) Calcule Cond (A) usando 11middot11 2 1111 y II to para las siguientes matrices
456 218) ( 279 138
b) Que puede decir del condicionamiento de los siguientes sistemas de ecuaciones lineales
f39x + 16x2 = 55 456X l + 218x2 = 674 I) ii) bulll 68x + 29x2 = 97 279x + 138x 2 = 417
ESTABILIDAD NUMERICA EN LA ELiMINACION GAUSSIANA
Volvamos al metodo de eliminacion Gaussiana (simple) y consideremos el siguiente ejemplo
Ejemplo 33 Resuelva el siguiente sistema de ecuaciones lineales usando eliminacion Gaussiana con sustitucion regresiva y aritmetica (decimal) con redondeo a tres dfgitos
El 03x + 589x2 = 592 E2 531x - 610x2 = 470
Usando elimnaci6n Gaussiana obtenemos
03 589 592) ~~- (A b) = ( 531 - 610 470
y por sustituci6n regresiva
_ 592 - 589(101 592 Xl = 03
luego la solucion calculada es X (~2 1 =(shy
Instrucci6n en DERIVE
PIVOT(A i j) Usa oper~cione~e tiene matriz obtenida de la matnz A q
Que puede decir de la calidad de la
Para intentar responder esta
dadas par el teorema 3 3
03 589) Como A = ( 531 -610
los calculos se obtienen III
v
(311 )
Capitulo 3 SOlUCION NUMERICA DE SISTEMAS DE ECUACIONES 111
(A b) = ( 03 589 592) __( 5__) _m 2 _=1_77 ) ~) ( 03 589E2_ ~~ _( _ _ 592 ) 531 -610 470 0 - 10400 -10500
y por sustituci6n regresiva
x = -10500 =101 2 - 10400
)(1 = 592 - 589(101) = 592 - 595 = - 3 = - 100 D3 D3 D3
- (X1 J (-1 00)luego la soluci6n calculada es X = _ = 101x2
Instrucci6n en DERIVE
PIVOT(A i j) Usa operaciones elemetales de fila para Simplificar (0 aproXimar) en una matriz obtenida de la matriz A que tiene ceros en la columna j y por debajo de la fila i 0
Que puede decir de la calidad de la soluci6n aproximada X
IIX -XIIPara intentar responder esta pregunta encontremos las cotas para el error relativo
IIXII dadas por el teorema 33
03 589)Como A = entonces usando aritmetica con redondeo a tres digitos para todos( 531 - 610
los calculos se obtienen las siguientes aproximaciones
- 1 1 (-610 - 589)A =-shy-313 - 531 03
II A II = Max 589 114 = 589
A - 1 II = _ 1_ Max 650 534 = _ 1_ 650 = 208II 313 313
entonces Cond (A) = II A II II A -111 = (589)(208) = 123 que no es muy grande
comparado con uno asi que la matriz A puede considerarse bien condicionada
(Por ciertas consideraciones te6ricas sobre el numero de condici6n de una matriz A las cuales pueden ser consultadas en Burden 1985 paginas 481 y 482 cuando se trabaja en
aritmetica finita (decimal) con redondeo a t-digitos y Cond(A) ~ 10t se espera un mal
112 METOOOS NUMERICOS
comportamiento de A con respecto a la solucion de AX = b Y A se considera mal
condicionada En este ejemplo Cond(A) = 123 lt 103 )
Ahora
03 589 ) (-100) ( 592)AX-b e ( 531 - 610 101 470
Estas operaciones se realizari en doble preci sion (6 digitos)
59189) ( 592) ( - 011 )(= -59261 - 470 = -106261
(Para evitar la perdida de cifras sign ificativas se debe calcular el vector erro r res idual
R = AX - b en doble precision)
Conv irtiendo este ultimo resultado a tres digltos usando redondeo se obtiene
R = (-011) - 106
Entonces II R II =106 y II b II = 592 Y por tanto
II R la 1 106 1 11gtlt - X II 106 II R II en
-II-11- ( )= - - =145 ~ II II ~ 2202= 123- = Cond (A)-II - lI shyb co Cond A 592 123 X 592 b DO
II R pero como Cond (A) ~ 1 entonces se espera que fbiC Y X- X sean mas 0 menosII xL
II X- X L sea tamble ndel mismo tamano y ya que es grande se espera que II X II
grande +
En situaciones como a observada en este ejemplo se sugiere hacer un refinamiento
iterativo sobre la soluci6n caculada X 0 usar esta soluci6n calcuada como aproximaci6n inicial en un metodo iterativ~ con et proposito de tratar de mejorar la soluci6n aproximada y ograr que el error residua l relativo sea mas pequeno EI metodo de reflna81iento iterativo puede ser consutado en Kincaid 1972 pag inas 174-176
La soucion exacta del sistema en consideracion es X = (J= ( 0) A que se debe a
diferencia entre la soucion exacta y a solucion calculada
Observe que e error en e calculo de )(2 con respecto a x2 fue de solo 01 (un error relativo
de 1) Y este error fue multiplicado par un factor de aproximadamente - 2000 al obtener
)(1 debido a orden en que se realiz6 la el iminaci6n Gaussiana
Instruccion en DERIVE
RESUELVA_1(Ab) Simplifica en a SOUCI6n entra como un vector fila 0
En el ejemplo anterior la elimina~i6n sistema de ecuaciones lineales ble~ nII1II
del algoritmo de eliminaci6n Gausslana pequeno) Hay sin embargo SIlUC -
Gaussiana es numencamente estable consultarse en Burden 1985 paginas366 y
Teorema 35 Si A =(ajiLn (EDD) par filas es decir si
entonces A es invertible (no-singular intercambio de filas en cualquier calculos son estables con respeclO
Notese que como consecuencia
por filas entonces A t~ene unos en su diagonal pnnClpalyU
~~a~ (vease Burden 19a5
satisface
Capitulo 3 SOLUCION NUMERIC A DE SISTEMAS DE ECUACIONES 113
Instrucci6n en DERIVE
RESUELVA_1(Ab) Simplifica en la soluci6n exacta X del sistema AX=b EI vector b se entra como un vector fila 0
En el ejemplo anterior la eliminaci6n Gaussiana condujo a una respuesta defectuosa de un sistema de ecuaciones lineales bien condicionado Esto muestra la inestabilidad numerica del algoritmo de eliminaci6n Gaussiana (consecuencia de la divisi6n por un nOmero (pivote) pequero) Hay sin embargo situaciones en las cuales el algoritmo de eliminaci6n Gaussiana es numericamente estable EI siguiente teorema cuya demostraci6n puede consultarse en Burden 1985 paginas 366 y 367 se refiere a una de tales situaciones
Teorema 35 Si A = (ai j )n xn es una matriz estrictamente dominante diagonalmente
(EDD) por filas es decir si
n Iaii Igt II ai j I para cada i == 12 n
)= 1 ) 1
entonces A es invertible (no-singular) Ademas se puede rea lizar eliminaci6n Gaussiana sin intercambio de filas en cualquier sistema AX = b para obtener su Onica soluci6n y los calculos son estables con respecto al crecimiento de los errores de redondeo V
N6tese que como consecuencia del teorema anterior se tiene que Si A = (a l ) )n xn es EDD
por filas entonces A tiene factorizaci6n LU es decir A == LU con L triangular inferior con unos en su diagonal principal y U triangular superior (escalonada)
sea tambien
Observe que la matriz de coeficientes del ejemplo 33 anterior no es EDD por filas
sectJeor9~~am9i~ -valioo=paramatrj simetricas
(vease Bllrden 19 5 ~ina 368) Una matriz A E Rnxn ~ i~etrica se dice definida
positiva si satisface una cualquiera de las sig uientes condiciones (Ias- cualesson equival~tes) a _-----
- - -- shy
i) ~AX gt 0 paralodo 2CERn X 0 J
iiXTodos los valores propios de A son pos1tlvOS7
iii) Todos los pivotes obtenidos en la eliminaci6n Gaussiana sobre A sin intercambio de filas - ~on positiv~s
iv) Todas las submatrices principales de A tienen determinante positiv~
(Las submatrices principales de la matriz A == (ai J) son las matrices - n ~ n
1k a ] k == 12bull f) )
akk v
114 METODOS NUMERICOS
N6tese nuevamente que Si A E R nxn es simetrica y definida positiva entonces A tiene
factorizaci6n A = LU con L triangular inferior con sus componentes sabre la diagonal principal iguales a uno y U triangular superior (ver mas adelante factorizaci6n de Choleski)
Observe que la matriz de coeficientes de l ejemplo 33 anterior no es simetrica
34 ESTRATEGIAS DE PIVOTEO
EI eJemplo 33 anterior muestra una de las dificultades que pueden surgir al aplicar el metodo
de eliminaci6n Gaussiana cuando el pivote ajj-1) es pequeno comparado con algunos
(j-1) e ementos I a t para J ~ I t ~ n
Para tratar de evitar tales dificultades se introduce en el metodo de eliminaci6n Gaussiana una estrategia IIamada de pivoteo la cual consiste en seleccionar el pivote de acuerdo con un cierto criter io Nosotros usaremos dos estrategias la estrategia de pivoteo maximo por columna 0 pivoteo parcial y la estrategia de pivoteo escalado de fila 0 escalamiento
341 Pivoteo maximo por columna 0 pivoteo parcial Esta estrategia difiere de
eliminaci6n Gaussiana simple unicamente en la escogencia del pivote ajj-1) la cual se hace
ahora as
Para j = 12 n - 1 se determina el menor entero k J ~ k ~ n tal que
y
es decir seleccionamos el primer elemento diferente de cero sobre la columna j-esima a
partir de la j-esima fila y que tenga mayor valor absoluto (para j = 1 a~j j1 ) = al~) ak 1 )
Si tal k no existe el sistema no tiene soluci6n unica y el proceso se puede terminar
Si tal k existe y k ot j entonces hacemos intercambio de las ecuaciones j-esima y k-esima
y continuamos con la eliminaci6n Gaussiana l
Ilustremos esta estrategia para resolver el sistema
E1 03X1 + 589x2 = 592 E2 531x1 - 61 OX 2 = 470
entonces k = 2 ot 1= j asl que
Por sustituci6n regresiva
Observe que en este caSO
SWAP(A i j)
que es el mismo del ejemplo 33 usando aritmetica con redondeo a tres digitos
Como para j = 1 se tiene que
Capitulo 3 SOLUCION NUMERICA DE SISTEMAS DE ECUACIONES 115
Max Iall I Ia21 I = Max 03 531 = 531 = Ia21 I 0
entonces k = 2 1= j as i que intercambiamos El y E2 Y continuamos con la el imi nac i6n
531 - 610 (A b)= ( 03 589 592) 470) (531 - 610 470 03 589 592
E2I (~) ( m21=5 65 ~ 1O - l) 531 - 610 5 31 470) ------------------~) ( o 589 589
Por sustituci6n regresiva obtenemos
x = 589 = 100 Xl = 470 +610(100) = 531 = 100 2 589 531 531
Observe que en este caso X= (J es la soluci6n exacta del sistema dado bull
Instrucci6n en DERIVE
SWAP(A i j) Intercambia las filas (0 elementos) i y j de la matriz A (de un vector) 0
Nota En el procedimiento de pivoteo maximo por columna (p ivoteo parcial) cada multiplicador mi j es tal que
y aunque esta estrategia permite resolver satisfactoriamente muchos sistemas de ecuaciones lineales hay casos donde fracasa como se ilustra en el siguiente ejemplo
Ejemplo 34 Consideremos el sistema
El 300Xl + 58900x2 59200 E2 531x1 -610X2 - 470
el cual es un sistema equivalente al del ejemplo 33 (los coeficientes de la primera ecuaci6n
en el sistema del ejemplo 33 han sido multiplicados por 10 3) EI pivoteo maximo por
columna con aritmetica de redondeo a tres dig itos nos lIeva a los siguientes resultados
59200) E2 1 (-~pound) (m 2 1 =177) ) ( 300 58900 59200 ) (A b) = ( 300 58900 531 - 610 470 o - 10400 -10500
y por sustituci6n regresiva x2 = 101 Y Xl = - 100 que es la misma solucion que se obtiene si ~ usamos eliminaci6n Gaussiana simple
116 METODOS NUMERICOS
En casas como el de este ejemplo donde un pivote es mucho mas pequeno que alguno de los coeficientes de la ecuacion que el encabeza se recomienda la tecnica conocida como pivoteo escalado de fila 0 escalamiento la cual es nuestra segunda estrategia
342 Pivoteo escalado de fila Esta tecnica solo difiere de la eliminaci6n Gaussiana simple al igual que el pivoteo parcial en la escogencia del pivote
lEsta vez el pivote alr ) se escoge como se indica a continuacion
Para j = 12 n - 1 hacemos 10 siguiente
a) Para i = j j + 1 n calculamos
S = Max Ia(l- l) I Factor de escala I j s i s n II
Si Sj = 0 entonces el sistema no tiene solucion unica y el proceso se puede terminar
b) Para i = jj+1 n calculamos
c) Encontramos el menor entero k con j o k 0 n tal que
Ia~jjl ) I IaH- l) I -------- = Max --------
Sk j S j $ n SI
Si tal k no existe entonces el sistema no tiene solucion unica y el proceso se puede terminar
Si tal k existe y kF j entonces hacemos intercambio de las ecuaciones j-esima y kshyesima
y continuamos con la eliminacion Gaussiana V
Apliquemos esta estrategia para resolver el sistema del ejemplo 34 usando aritmetica con redondeo a tres dfgitos
Para j = 1
a) Sl = Max Iall I I a12 I = Max 300 58900 = 58900F 0 Y
Max 53~ 610 IshyS2 =Max a21 a n -shy
b) Ahora ~= 51
la0
l 300
S2
531 al l I a21 -Max - 0 -
c) Max -S -S- - 5890 NIl 1 2
las ecuaciolllltanto intercamblamos
Gaussiana
300 58900 (A b) = ( 531 -610
Por sustitucion regresiva
Xl tOO
Observe que en este GaSO X
35 FACTORIZACION
Capitulo 3 SOLUCION NUMERICA DE SISTEMAS DE ECUACIONES 117
8 2 =Max I a21 I I a22 1 =Max 531 610 = 610 0
b) Ahora
Ia I 300 - 8- = 58900 Y
1
Ia 2 I 531 - - = - shy
6108 2
c) Max ~ ~ I=Max ~ 531 = 531 =~ 0 as que k = 2 1= j y por8 8 2 58900 60 60 8 2
tanto intercambiamos las ecuaciones E y E2 Y continuamos con la eliminacion
Gaussiana
(A b) = (300 58900 59200) p ( 531 -610 470 ) 2 )
531 - 610 470 300 58900 59200
E2 (300) (m2 =565 ) ( 531 -6 1 0531 47 0 )
--------------~) o 58900 58900
Por sustituci6n regresiva
_ 470 +610(100) 531X2 = 100 Xl = =-- = 100
531 531
Observe que en este caso X= ( ~ J = ( 100) es la solucion exacta del sistema bull 100x2
35 FACTORIZACION TRIANGULAR
Consideremos un sistema AX = b con A no-singular y b 0 Con respecto a la matriz A se sabe que existen matrices P de permutaci6n L triangular inferior con sus componentes sobre la diagonal principal iguales a uno y U triangular superior (escalon ada) tales que PA = LU Una forma eficiente computacionalmente de encontrar P L Y U usando eliminaci6n Gaussiana se muestra en el siguiente ejemplo Ejemplo 35 Resuelva el siguiente sistema de ecuaciones lineales usando eliminacion Gaussiana con pivoteo parcial y obtenga una factorizaci6n PA =LU para la matriz A de coeficientes asociada con este metodo
- x + 2X2 - 3X3 = -2 3x1 - 3x2 - X 3 = -4
x + x2 = 3IEmpezamos introduciendo un vector p =(PP2P3) T el cual inicializamos_con P = i i = 123 Y
donde se almacenaran los intercambios necesarios en el proceso de eliminacion Gaussiana
118 METODOS NUMERICOS
con pivoteo parcial (el numero de componentes del vector p coincide con el orden n del sistema a resolver)
- 1 2 -3 -2] -4] - 2 P = (123) T (Ab) = ~ - 3 -1 - 4
[ o 3 3
3 - 3 - 1 - 4(m=-H (m3l = ~ )
E( - ~) E3(~) 10 10(-i) 3 3
2(i) 1 13 3 3
(Observe que cada multiplicador mij es almacenado en la posici6n correspondiente (i j) en
la matriz de trabajo)
3 - 3 -1 - 4
Max 121 =2 p=( 23 1) T bull (i) 2 1
-
3 13 3
(-i) 10 3
10 3
(Observe que la permutaci6n se hace para las filas 2 y 3 completas es decir incluyendo los multiplicadores)
3 - 3 - 1 --4(mJ2 = ~ ) E3( ~)
2(i) 13 3 3
(- i) (i) 7 11 2 2
La eficiencia en el metoda indicado se debe a que en la misma matriz de trabajo se almacenan los multiplicadores que van a conformar la matriz L (en el ejemplo son los numeros que se encuentran dentro de parentesis) 10 que significa un ahorro de memoria y como los intercambios necesarios afectan simultaneamente a las matrices L y U se evita tener que volver a la matriz original a realizar los intercambios ya observados y repetir la eliminaci6n Gaussiana con pivoteo parcial De esta manera al terminar el proceso de eliminaci6n podemos leer en la matriz final la parte estrictamente triangular inferior de L (son los numeros entre parentesis) y la matriz triangular superior U (que es la parte triangular superior de la matriz final) y en el vector p final quedan almacenados los intercambios realizados que se usan para producir la matriz de permutaci6n P
Para el ejemplo 35
3 -31 0 0 1 0 2U=0L= 3 1 1 0 0-3 2
(Verifique que PA =LU)
Para obtener la soluci6n del sistema
reducido
yobtenemos
C6mo se resuelve el
Capitulo 3 SOLUCION NUMERICA DE SISTEMAS DE ECUACIONES 119
Para el ejemplo 35
1 0 0 3 -3 - 1 1 1 1
L= - 0 U= 0 2 - Y como P~ (231) entonces [~ 0 lp3 3 1 1 7 0
-- - 0 0 3 2 2
(Verifique que PA = LU )
Para obtener la soluci6n del sistema original usamos sustituci6n regresiva en el sistema reducido
-4 13
UX = 3
11
2 ~
c y obtenemos
asi que la soluci6n (exacta) del sistema dado es X = (23 40 ~)T21 21 7
C6mo se resuelve el sistema AX = b a partir de la factorizaci6n PA = LU obtenida bull
351 Algunas aplicaciones de la factorizaci6n PA = LU La factorizaci6n PA = LU es utilizada eficientemente en aquellos casos donde se trabaja repetfdamente con la misma matriz A Dos de esos casos se presentan a continuaci6n
1) Resolver varios sistemas AX =b con la misma matriz de coeficientes A ya que en p L Y U esta almacenado todo el proceso de eliminaci6n Gaussiana EI algoritmo se bas a en la siguiente equivalencia
AX = b lt=gt PAX = Pb lt=gt LUX = Pb
lt=gt UX = L- 1Pb
lt=gtUXoo c y c oo L- 1Pb
lt=gt UX == C Y Lc = Pb lt=gt Lc = Pb Y UX = c
Los pasos a seguir son
Paso 1 Calcular Pb Paso 2 Resolver para c Lc = Pb por sustituci6n progresiva Paso 3 Resolver para X UX = c por sustituci6n regresiva
Como ejercicio resuelva el sistema del ejemplo anterior usando este algoritmo
120 METODOS NUMERICOS
2) Encontrar la matriz inversa A -1 de una matriz invertible A resolviendo los n-sistemas
AX =eU) j = 12 n
donde e(1) =(O O1 o 0f ERn
t posicion j
La solucion X del sistema AX = e(n j = 12 n produce la correspondiente columna j shy
esima de la matriz A -1
Como ejercicio compare el numero de operaciones para encontrar A - usando el metodo de Gauss-Jordan con el numero de operaciones resolviendo los n-sistemas indicados antes V
Ejercicio 33 Calcule la inversa de la matriz A de coeficientes del sistema del ejemplo 35 usando el metodo de Gauss-Jordan y tambien usando la factorizacion PA = LU bull
36 SISTEMAS TRIDIAGONALES
Un caso muy importante de sistemas de ecuaciones lineales que requiere un tratamiento especial es el de los sistemas tridiagonales Tales sistemas aparecen en diversas aplicaciones como por ejemplo al utilizar metodos de diferencias flnitas en la solucion de problemas con valores en la frontera para ecuaciones diferenciales ordinarias y como veremos mas adelante en el problema de la interpolacion segmentaria cubica
Un sistema tridiagonal es de la forma
dx1 + C1X2
a 2 x + d2 x2 + C 2 X 3
a3 x2 + d3 x3 + C 3 X 4
an-1Xn- 2 + dn_ xn_ + cn_xn = bn_1
anxn_1 + dnxn = bn
La matriz de coeficientes del sistema es
C 0 0 od1
0a 2 d2 c2
o a3 d3 c 3
A =
o
la cual se dice una matriz TRIDIAGONAL
En general
ij=12 n
suponiendo que la matriz A de podemos usar eliminaclon Gausslana eliminacion Gaussiana simple es adecuadl EDD por filas) Otra forma de resolver
EDD por filas es a partir de la componentes de las matrices L Y U as triangular inferior con todos sus superior Para encontrar tales malnce5
o 0 0
Y2 o o o o
L =
o Yn t
o o
Como Cl Ct
Y2 Cl Y2C +az
0 h a2
LU =
0
0
d
=
numero de
Capitulo 3 SOLUCION NUMERICA DE SISTEMAS DE ECUACIONES 109
Del teorema 32 se sabe que p (A) II A II para toda norma matricial inducida asl que
Cond (A) = II A IIII A -1 ~p(A) p(A-1)
pero como los valores propios de A - 1 son los reciprocos de los valores propios de A se tiene que
Max I A I ) Ecr(A) ( )
Cond (A) ~ I I == Condo AMm A
AEcr( A)
con a(A) =AE C A esvalorpropiodeA espectro de A (Recuerde que
p(A -1) = Max I A I = 1 I ) ecr(A-1) Mm A
AEcr( A)
Max I A I EI numero Cond (A)= AEcr(A) se denomina numero de condici6n espectral de A
bull Min I A I AEcr( A)
Segun se acaba de probar Cond(A) ~ Cond(A)
(1 1J Para la matnz A = se bene que 1005 10
det(A_AI) = 1- A 1 I =A2- 11A- 05 1005 10 - A
as que los valores propios de A son A1 1100454358 A 2 -45435778 x 10 -3 Y por tanto
Condo (A) 1100454358 2421999592 raquo 1 bull 45435778 x 10-3
Dado un sistema AX = b si 8A Y 8b denotan perturbaciones en A y b respectivamente el siguiente teorema cuya demostracion puede ser consultada en Ortega 1990 paginas 32 y
33 establece una cota para el error relativo II X- X II en terminos de las perturbaciones
IIXII
11 MII 118b II relatlvas lfAl lfbf y Cond(A) donde X es la soluclon exacta de AX = b Y X es la
solucion exacta del sistema perturbado (A + 8 A)X = b + 8 b
110 METODOS NUMERICOS
Teorema 34 Supongase que A es no-singular y que 11 8 A II lt II A~ II (esta hipotes is asegura
que A + 8 A es invertible y que 1- Cond(A) li SA II gt 0) Si X es la solucion exacta del IIA II
sistema perturbado (A + 0 A)X = b + 8 b entonces X aproxima a la solucion exacta X del
sistema AX = b b cI 0 con la siguiente estimacion de error
II X - X II ~ Cond (A) ( IL~~JL~ J (311)
IIXII 1_cond(A)[ 13 All l llb ll IIAII
II AII
La desigualdad (3 11) dice que si la matriz A esta bien condicionada es decir si
Cond (A) 1 entonces cambios pequerios en A y b producen correspondientemente
cambios pequerios en la soluci6n del sistema (el sistema AX = b esta bien condicionado) Por otro lado si A esta mal condicionada entonces cam bios pequenos en A y b pueden producir grandes cam bios en la solucion del sistema (el sistema AX = b esta mal condicionado)
Ejercicio 31 Estime la cota de error dada en el teorema 34 para los sistemas (34) y (34) del ejemplo 31 bull
Ejercicio 32 a) Calcule Cond (A) usando 11middot11 2 1111 y II to para las siguientes matrices
456 218) ( 279 138
b) Que puede decir del condicionamiento de los siguientes sistemas de ecuaciones lineales
f39x + 16x2 = 55 456X l + 218x2 = 674 I) ii) bulll 68x + 29x2 = 97 279x + 138x 2 = 417
ESTABILIDAD NUMERICA EN LA ELiMINACION GAUSSIANA
Volvamos al metodo de eliminacion Gaussiana (simple) y consideremos el siguiente ejemplo
Ejemplo 33 Resuelva el siguiente sistema de ecuaciones lineales usando eliminacion Gaussiana con sustitucion regresiva y aritmetica (decimal) con redondeo a tres dfgitos
El 03x + 589x2 = 592 E2 531x - 610x2 = 470
Usando elimnaci6n Gaussiana obtenemos
03 589 592) ~~- (A b) = ( 531 - 610 470
y por sustituci6n regresiva
_ 592 - 589(101 592 Xl = 03
luego la solucion calculada es X (~2 1 =(shy
Instrucci6n en DERIVE
PIVOT(A i j) Usa oper~cione~e tiene matriz obtenida de la matnz A q
Que puede decir de la calidad de la
Para intentar responder esta
dadas par el teorema 3 3
03 589) Como A = ( 531 -610
los calculos se obtienen III
v
(311 )
Capitulo 3 SOlUCION NUMERICA DE SISTEMAS DE ECUACIONES 111
(A b) = ( 03 589 592) __( 5__) _m 2 _=1_77 ) ~) ( 03 589E2_ ~~ _( _ _ 592 ) 531 -610 470 0 - 10400 -10500
y por sustituci6n regresiva
x = -10500 =101 2 - 10400
)(1 = 592 - 589(101) = 592 - 595 = - 3 = - 100 D3 D3 D3
- (X1 J (-1 00)luego la soluci6n calculada es X = _ = 101x2
Instrucci6n en DERIVE
PIVOT(A i j) Usa operaciones elemetales de fila para Simplificar (0 aproXimar) en una matriz obtenida de la matriz A que tiene ceros en la columna j y por debajo de la fila i 0
Que puede decir de la calidad de la soluci6n aproximada X
IIX -XIIPara intentar responder esta pregunta encontremos las cotas para el error relativo
IIXII dadas por el teorema 33
03 589)Como A = entonces usando aritmetica con redondeo a tres digitos para todos( 531 - 610
los calculos se obtienen las siguientes aproximaciones
- 1 1 (-610 - 589)A =-shy-313 - 531 03
II A II = Max 589 114 = 589
A - 1 II = _ 1_ Max 650 534 = _ 1_ 650 = 208II 313 313
entonces Cond (A) = II A II II A -111 = (589)(208) = 123 que no es muy grande
comparado con uno asi que la matriz A puede considerarse bien condicionada
(Por ciertas consideraciones te6ricas sobre el numero de condici6n de una matriz A las cuales pueden ser consultadas en Burden 1985 paginas 481 y 482 cuando se trabaja en
aritmetica finita (decimal) con redondeo a t-digitos y Cond(A) ~ 10t se espera un mal
112 METOOOS NUMERICOS
comportamiento de A con respecto a la solucion de AX = b Y A se considera mal
condicionada En este ejemplo Cond(A) = 123 lt 103 )
Ahora
03 589 ) (-100) ( 592)AX-b e ( 531 - 610 101 470
Estas operaciones se realizari en doble preci sion (6 digitos)
59189) ( 592) ( - 011 )(= -59261 - 470 = -106261
(Para evitar la perdida de cifras sign ificativas se debe calcular el vector erro r res idual
R = AX - b en doble precision)
Conv irtiendo este ultimo resultado a tres digltos usando redondeo se obtiene
R = (-011) - 106
Entonces II R II =106 y II b II = 592 Y por tanto
II R la 1 106 1 11gtlt - X II 106 II R II en
-II-11- ( )= - - =145 ~ II II ~ 2202= 123- = Cond (A)-II - lI shyb co Cond A 592 123 X 592 b DO
II R pero como Cond (A) ~ 1 entonces se espera que fbiC Y X- X sean mas 0 menosII xL
II X- X L sea tamble ndel mismo tamano y ya que es grande se espera que II X II
grande +
En situaciones como a observada en este ejemplo se sugiere hacer un refinamiento
iterativo sobre la soluci6n caculada X 0 usar esta soluci6n calcuada como aproximaci6n inicial en un metodo iterativ~ con et proposito de tratar de mejorar la soluci6n aproximada y ograr que el error residua l relativo sea mas pequeno EI metodo de reflna81iento iterativo puede ser consutado en Kincaid 1972 pag inas 174-176
La soucion exacta del sistema en consideracion es X = (J= ( 0) A que se debe a
diferencia entre la soucion exacta y a solucion calculada
Observe que e error en e calculo de )(2 con respecto a x2 fue de solo 01 (un error relativo
de 1) Y este error fue multiplicado par un factor de aproximadamente - 2000 al obtener
)(1 debido a orden en que se realiz6 la el iminaci6n Gaussiana
Instruccion en DERIVE
RESUELVA_1(Ab) Simplifica en a SOUCI6n entra como un vector fila 0
En el ejemplo anterior la elimina~i6n sistema de ecuaciones lineales ble~ nII1II
del algoritmo de eliminaci6n Gausslana pequeno) Hay sin embargo SIlUC -
Gaussiana es numencamente estable consultarse en Burden 1985 paginas366 y
Teorema 35 Si A =(ajiLn (EDD) par filas es decir si
entonces A es invertible (no-singular intercambio de filas en cualquier calculos son estables con respeclO
Notese que como consecuencia
por filas entonces A t~ene unos en su diagonal pnnClpalyU
~~a~ (vease Burden 19a5
satisface
Capitulo 3 SOLUCION NUMERIC A DE SISTEMAS DE ECUACIONES 113
Instrucci6n en DERIVE
RESUELVA_1(Ab) Simplifica en la soluci6n exacta X del sistema AX=b EI vector b se entra como un vector fila 0
En el ejemplo anterior la eliminaci6n Gaussiana condujo a una respuesta defectuosa de un sistema de ecuaciones lineales bien condicionado Esto muestra la inestabilidad numerica del algoritmo de eliminaci6n Gaussiana (consecuencia de la divisi6n por un nOmero (pivote) pequero) Hay sin embargo situaciones en las cuales el algoritmo de eliminaci6n Gaussiana es numericamente estable EI siguiente teorema cuya demostraci6n puede consultarse en Burden 1985 paginas 366 y 367 se refiere a una de tales situaciones
Teorema 35 Si A = (ai j )n xn es una matriz estrictamente dominante diagonalmente
(EDD) por filas es decir si
n Iaii Igt II ai j I para cada i == 12 n
)= 1 ) 1
entonces A es invertible (no-singular) Ademas se puede rea lizar eliminaci6n Gaussiana sin intercambio de filas en cualquier sistema AX = b para obtener su Onica soluci6n y los calculos son estables con respecto al crecimiento de los errores de redondeo V
N6tese que como consecuencia del teorema anterior se tiene que Si A = (a l ) )n xn es EDD
por filas entonces A tiene factorizaci6n LU es decir A == LU con L triangular inferior con unos en su diagonal principal y U triangular superior (escalonada)
sea tambien
Observe que la matriz de coeficientes del ejemplo 33 anterior no es EDD por filas
sectJeor9~~am9i~ -valioo=paramatrj simetricas
(vease Bllrden 19 5 ~ina 368) Una matriz A E Rnxn ~ i~etrica se dice definida
positiva si satisface una cualquiera de las sig uientes condiciones (Ias- cualesson equival~tes) a _-----
- - -- shy
i) ~AX gt 0 paralodo 2CERn X 0 J
iiXTodos los valores propios de A son pos1tlvOS7
iii) Todos los pivotes obtenidos en la eliminaci6n Gaussiana sobre A sin intercambio de filas - ~on positiv~s
iv) Todas las submatrices principales de A tienen determinante positiv~
(Las submatrices principales de la matriz A == (ai J) son las matrices - n ~ n
1k a ] k == 12bull f) )
akk v
114 METODOS NUMERICOS
N6tese nuevamente que Si A E R nxn es simetrica y definida positiva entonces A tiene
factorizaci6n A = LU con L triangular inferior con sus componentes sabre la diagonal principal iguales a uno y U triangular superior (ver mas adelante factorizaci6n de Choleski)
Observe que la matriz de coeficientes de l ejemplo 33 anterior no es simetrica
34 ESTRATEGIAS DE PIVOTEO
EI eJemplo 33 anterior muestra una de las dificultades que pueden surgir al aplicar el metodo
de eliminaci6n Gaussiana cuando el pivote ajj-1) es pequeno comparado con algunos
(j-1) e ementos I a t para J ~ I t ~ n
Para tratar de evitar tales dificultades se introduce en el metodo de eliminaci6n Gaussiana una estrategia IIamada de pivoteo la cual consiste en seleccionar el pivote de acuerdo con un cierto criter io Nosotros usaremos dos estrategias la estrategia de pivoteo maximo por columna 0 pivoteo parcial y la estrategia de pivoteo escalado de fila 0 escalamiento
341 Pivoteo maximo por columna 0 pivoteo parcial Esta estrategia difiere de
eliminaci6n Gaussiana simple unicamente en la escogencia del pivote ajj-1) la cual se hace
ahora as
Para j = 12 n - 1 se determina el menor entero k J ~ k ~ n tal que
y
es decir seleccionamos el primer elemento diferente de cero sobre la columna j-esima a
partir de la j-esima fila y que tenga mayor valor absoluto (para j = 1 a~j j1 ) = al~) ak 1 )
Si tal k no existe el sistema no tiene soluci6n unica y el proceso se puede terminar
Si tal k existe y k ot j entonces hacemos intercambio de las ecuaciones j-esima y k-esima
y continuamos con la eliminaci6n Gaussiana l
Ilustremos esta estrategia para resolver el sistema
E1 03X1 + 589x2 = 592 E2 531x1 - 61 OX 2 = 470
entonces k = 2 ot 1= j asl que
Por sustituci6n regresiva
Observe que en este caSO
SWAP(A i j)
que es el mismo del ejemplo 33 usando aritmetica con redondeo a tres digitos
Como para j = 1 se tiene que
Capitulo 3 SOLUCION NUMERICA DE SISTEMAS DE ECUACIONES 115
Max Iall I Ia21 I = Max 03 531 = 531 = Ia21 I 0
entonces k = 2 1= j as i que intercambiamos El y E2 Y continuamos con la el imi nac i6n
531 - 610 (A b)= ( 03 589 592) 470) (531 - 610 470 03 589 592
E2I (~) ( m21=5 65 ~ 1O - l) 531 - 610 5 31 470) ------------------~) ( o 589 589
Por sustituci6n regresiva obtenemos
x = 589 = 100 Xl = 470 +610(100) = 531 = 100 2 589 531 531
Observe que en este caso X= (J es la soluci6n exacta del sistema dado bull
Instrucci6n en DERIVE
SWAP(A i j) Intercambia las filas (0 elementos) i y j de la matriz A (de un vector) 0
Nota En el procedimiento de pivoteo maximo por columna (p ivoteo parcial) cada multiplicador mi j es tal que
y aunque esta estrategia permite resolver satisfactoriamente muchos sistemas de ecuaciones lineales hay casos donde fracasa como se ilustra en el siguiente ejemplo
Ejemplo 34 Consideremos el sistema
El 300Xl + 58900x2 59200 E2 531x1 -610X2 - 470
el cual es un sistema equivalente al del ejemplo 33 (los coeficientes de la primera ecuaci6n
en el sistema del ejemplo 33 han sido multiplicados por 10 3) EI pivoteo maximo por
columna con aritmetica de redondeo a tres dig itos nos lIeva a los siguientes resultados
59200) E2 1 (-~pound) (m 2 1 =177) ) ( 300 58900 59200 ) (A b) = ( 300 58900 531 - 610 470 o - 10400 -10500
y por sustituci6n regresiva x2 = 101 Y Xl = - 100 que es la misma solucion que se obtiene si ~ usamos eliminaci6n Gaussiana simple
116 METODOS NUMERICOS
En casas como el de este ejemplo donde un pivote es mucho mas pequeno que alguno de los coeficientes de la ecuacion que el encabeza se recomienda la tecnica conocida como pivoteo escalado de fila 0 escalamiento la cual es nuestra segunda estrategia
342 Pivoteo escalado de fila Esta tecnica solo difiere de la eliminaci6n Gaussiana simple al igual que el pivoteo parcial en la escogencia del pivote
lEsta vez el pivote alr ) se escoge como se indica a continuacion
Para j = 12 n - 1 hacemos 10 siguiente
a) Para i = j j + 1 n calculamos
S = Max Ia(l- l) I Factor de escala I j s i s n II
Si Sj = 0 entonces el sistema no tiene solucion unica y el proceso se puede terminar
b) Para i = jj+1 n calculamos
c) Encontramos el menor entero k con j o k 0 n tal que
Ia~jjl ) I IaH- l) I -------- = Max --------
Sk j S j $ n SI
Si tal k no existe entonces el sistema no tiene solucion unica y el proceso se puede terminar
Si tal k existe y kF j entonces hacemos intercambio de las ecuaciones j-esima y kshyesima
y continuamos con la eliminacion Gaussiana V
Apliquemos esta estrategia para resolver el sistema del ejemplo 34 usando aritmetica con redondeo a tres dfgitos
Para j = 1
a) Sl = Max Iall I I a12 I = Max 300 58900 = 58900F 0 Y
Max 53~ 610 IshyS2 =Max a21 a n -shy
b) Ahora ~= 51
la0
l 300
S2
531 al l I a21 -Max - 0 -
c) Max -S -S- - 5890 NIl 1 2
las ecuaciolllltanto intercamblamos
Gaussiana
300 58900 (A b) = ( 531 -610
Por sustitucion regresiva
Xl tOO
Observe que en este GaSO X
35 FACTORIZACION
Capitulo 3 SOLUCION NUMERICA DE SISTEMAS DE ECUACIONES 117
8 2 =Max I a21 I I a22 1 =Max 531 610 = 610 0
b) Ahora
Ia I 300 - 8- = 58900 Y
1
Ia 2 I 531 - - = - shy
6108 2
c) Max ~ ~ I=Max ~ 531 = 531 =~ 0 as que k = 2 1= j y por8 8 2 58900 60 60 8 2
tanto intercambiamos las ecuaciones E y E2 Y continuamos con la eliminacion
Gaussiana
(A b) = (300 58900 59200) p ( 531 -610 470 ) 2 )
531 - 610 470 300 58900 59200
E2 (300) (m2 =565 ) ( 531 -6 1 0531 47 0 )
--------------~) o 58900 58900
Por sustituci6n regresiva
_ 470 +610(100) 531X2 = 100 Xl = =-- = 100
531 531
Observe que en este caso X= ( ~ J = ( 100) es la solucion exacta del sistema bull 100x2
35 FACTORIZACION TRIANGULAR
Consideremos un sistema AX = b con A no-singular y b 0 Con respecto a la matriz A se sabe que existen matrices P de permutaci6n L triangular inferior con sus componentes sobre la diagonal principal iguales a uno y U triangular superior (escalon ada) tales que PA = LU Una forma eficiente computacionalmente de encontrar P L Y U usando eliminaci6n Gaussiana se muestra en el siguiente ejemplo Ejemplo 35 Resuelva el siguiente sistema de ecuaciones lineales usando eliminacion Gaussiana con pivoteo parcial y obtenga una factorizaci6n PA =LU para la matriz A de coeficientes asociada con este metodo
- x + 2X2 - 3X3 = -2 3x1 - 3x2 - X 3 = -4
x + x2 = 3IEmpezamos introduciendo un vector p =(PP2P3) T el cual inicializamos_con P = i i = 123 Y
donde se almacenaran los intercambios necesarios en el proceso de eliminacion Gaussiana
118 METODOS NUMERICOS
con pivoteo parcial (el numero de componentes del vector p coincide con el orden n del sistema a resolver)
- 1 2 -3 -2] -4] - 2 P = (123) T (Ab) = ~ - 3 -1 - 4
[ o 3 3
3 - 3 - 1 - 4(m=-H (m3l = ~ )
E( - ~) E3(~) 10 10(-i) 3 3
2(i) 1 13 3 3
(Observe que cada multiplicador mij es almacenado en la posici6n correspondiente (i j) en
la matriz de trabajo)
3 - 3 -1 - 4
Max 121 =2 p=( 23 1) T bull (i) 2 1
-
3 13 3
(-i) 10 3
10 3
(Observe que la permutaci6n se hace para las filas 2 y 3 completas es decir incluyendo los multiplicadores)
3 - 3 - 1 --4(mJ2 = ~ ) E3( ~)
2(i) 13 3 3
(- i) (i) 7 11 2 2
La eficiencia en el metoda indicado se debe a que en la misma matriz de trabajo se almacenan los multiplicadores que van a conformar la matriz L (en el ejemplo son los numeros que se encuentran dentro de parentesis) 10 que significa un ahorro de memoria y como los intercambios necesarios afectan simultaneamente a las matrices L y U se evita tener que volver a la matriz original a realizar los intercambios ya observados y repetir la eliminaci6n Gaussiana con pivoteo parcial De esta manera al terminar el proceso de eliminaci6n podemos leer en la matriz final la parte estrictamente triangular inferior de L (son los numeros entre parentesis) y la matriz triangular superior U (que es la parte triangular superior de la matriz final) y en el vector p final quedan almacenados los intercambios realizados que se usan para producir la matriz de permutaci6n P
Para el ejemplo 35
3 -31 0 0 1 0 2U=0L= 3 1 1 0 0-3 2
(Verifique que PA =LU)
Para obtener la soluci6n del sistema
reducido
yobtenemos
C6mo se resuelve el
Capitulo 3 SOLUCION NUMERICA DE SISTEMAS DE ECUACIONES 119
Para el ejemplo 35
1 0 0 3 -3 - 1 1 1 1
L= - 0 U= 0 2 - Y como P~ (231) entonces [~ 0 lp3 3 1 1 7 0
-- - 0 0 3 2 2
(Verifique que PA = LU )
Para obtener la soluci6n del sistema original usamos sustituci6n regresiva en el sistema reducido
-4 13
UX = 3
11
2 ~
c y obtenemos
asi que la soluci6n (exacta) del sistema dado es X = (23 40 ~)T21 21 7
C6mo se resuelve el sistema AX = b a partir de la factorizaci6n PA = LU obtenida bull
351 Algunas aplicaciones de la factorizaci6n PA = LU La factorizaci6n PA = LU es utilizada eficientemente en aquellos casos donde se trabaja repetfdamente con la misma matriz A Dos de esos casos se presentan a continuaci6n
1) Resolver varios sistemas AX =b con la misma matriz de coeficientes A ya que en p L Y U esta almacenado todo el proceso de eliminaci6n Gaussiana EI algoritmo se bas a en la siguiente equivalencia
AX = b lt=gt PAX = Pb lt=gt LUX = Pb
lt=gt UX = L- 1Pb
lt=gtUXoo c y c oo L- 1Pb
lt=gt UX == C Y Lc = Pb lt=gt Lc = Pb Y UX = c
Los pasos a seguir son
Paso 1 Calcular Pb Paso 2 Resolver para c Lc = Pb por sustituci6n progresiva Paso 3 Resolver para X UX = c por sustituci6n regresiva
Como ejercicio resuelva el sistema del ejemplo anterior usando este algoritmo
120 METODOS NUMERICOS
2) Encontrar la matriz inversa A -1 de una matriz invertible A resolviendo los n-sistemas
AX =eU) j = 12 n
donde e(1) =(O O1 o 0f ERn
t posicion j
La solucion X del sistema AX = e(n j = 12 n produce la correspondiente columna j shy
esima de la matriz A -1
Como ejercicio compare el numero de operaciones para encontrar A - usando el metodo de Gauss-Jordan con el numero de operaciones resolviendo los n-sistemas indicados antes V
Ejercicio 33 Calcule la inversa de la matriz A de coeficientes del sistema del ejemplo 35 usando el metodo de Gauss-Jordan y tambien usando la factorizacion PA = LU bull
36 SISTEMAS TRIDIAGONALES
Un caso muy importante de sistemas de ecuaciones lineales que requiere un tratamiento especial es el de los sistemas tridiagonales Tales sistemas aparecen en diversas aplicaciones como por ejemplo al utilizar metodos de diferencias flnitas en la solucion de problemas con valores en la frontera para ecuaciones diferenciales ordinarias y como veremos mas adelante en el problema de la interpolacion segmentaria cubica
Un sistema tridiagonal es de la forma
dx1 + C1X2
a 2 x + d2 x2 + C 2 X 3
a3 x2 + d3 x3 + C 3 X 4
an-1Xn- 2 + dn_ xn_ + cn_xn = bn_1
anxn_1 + dnxn = bn
La matriz de coeficientes del sistema es
C 0 0 od1
0a 2 d2 c2
o a3 d3 c 3
A =
o
la cual se dice una matriz TRIDIAGONAL
En general
ij=12 n
suponiendo que la matriz A de podemos usar eliminaclon Gausslana eliminacion Gaussiana simple es adecuadl EDD por filas) Otra forma de resolver
EDD por filas es a partir de la componentes de las matrices L Y U as triangular inferior con todos sus superior Para encontrar tales malnce5
o 0 0
Y2 o o o o
L =
o Yn t
o o
Como Cl Ct
Y2 Cl Y2C +az
0 h a2
LU =
0
0
d
=
110 METODOS NUMERICOS
Teorema 34 Supongase que A es no-singular y que 11 8 A II lt II A~ II (esta hipotes is asegura
que A + 8 A es invertible y que 1- Cond(A) li SA II gt 0) Si X es la solucion exacta del IIA II
sistema perturbado (A + 0 A)X = b + 8 b entonces X aproxima a la solucion exacta X del
sistema AX = b b cI 0 con la siguiente estimacion de error
II X - X II ~ Cond (A) ( IL~~JL~ J (311)
IIXII 1_cond(A)[ 13 All l llb ll IIAII
II AII
La desigualdad (3 11) dice que si la matriz A esta bien condicionada es decir si
Cond (A) 1 entonces cambios pequerios en A y b producen correspondientemente
cambios pequerios en la soluci6n del sistema (el sistema AX = b esta bien condicionado) Por otro lado si A esta mal condicionada entonces cam bios pequenos en A y b pueden producir grandes cam bios en la solucion del sistema (el sistema AX = b esta mal condicionado)
Ejercicio 31 Estime la cota de error dada en el teorema 34 para los sistemas (34) y (34) del ejemplo 31 bull
Ejercicio 32 a) Calcule Cond (A) usando 11middot11 2 1111 y II to para las siguientes matrices
456 218) ( 279 138
b) Que puede decir del condicionamiento de los siguientes sistemas de ecuaciones lineales
f39x + 16x2 = 55 456X l + 218x2 = 674 I) ii) bulll 68x + 29x2 = 97 279x + 138x 2 = 417
ESTABILIDAD NUMERICA EN LA ELiMINACION GAUSSIANA
Volvamos al metodo de eliminacion Gaussiana (simple) y consideremos el siguiente ejemplo
Ejemplo 33 Resuelva el siguiente sistema de ecuaciones lineales usando eliminacion Gaussiana con sustitucion regresiva y aritmetica (decimal) con redondeo a tres dfgitos
El 03x + 589x2 = 592 E2 531x - 610x2 = 470
Usando elimnaci6n Gaussiana obtenemos
03 589 592) ~~- (A b) = ( 531 - 610 470
y por sustituci6n regresiva
_ 592 - 589(101 592 Xl = 03
luego la solucion calculada es X (~2 1 =(shy
Instrucci6n en DERIVE
PIVOT(A i j) Usa oper~cione~e tiene matriz obtenida de la matnz A q
Que puede decir de la calidad de la
Para intentar responder esta
dadas par el teorema 3 3
03 589) Como A = ( 531 -610
los calculos se obtienen III
v
(311 )
Capitulo 3 SOlUCION NUMERICA DE SISTEMAS DE ECUACIONES 111
(A b) = ( 03 589 592) __( 5__) _m 2 _=1_77 ) ~) ( 03 589E2_ ~~ _( _ _ 592 ) 531 -610 470 0 - 10400 -10500
y por sustituci6n regresiva
x = -10500 =101 2 - 10400
)(1 = 592 - 589(101) = 592 - 595 = - 3 = - 100 D3 D3 D3
- (X1 J (-1 00)luego la soluci6n calculada es X = _ = 101x2
Instrucci6n en DERIVE
PIVOT(A i j) Usa operaciones elemetales de fila para Simplificar (0 aproXimar) en una matriz obtenida de la matriz A que tiene ceros en la columna j y por debajo de la fila i 0
Que puede decir de la calidad de la soluci6n aproximada X
IIX -XIIPara intentar responder esta pregunta encontremos las cotas para el error relativo
IIXII dadas por el teorema 33
03 589)Como A = entonces usando aritmetica con redondeo a tres digitos para todos( 531 - 610
los calculos se obtienen las siguientes aproximaciones
- 1 1 (-610 - 589)A =-shy-313 - 531 03
II A II = Max 589 114 = 589
A - 1 II = _ 1_ Max 650 534 = _ 1_ 650 = 208II 313 313
entonces Cond (A) = II A II II A -111 = (589)(208) = 123 que no es muy grande
comparado con uno asi que la matriz A puede considerarse bien condicionada
(Por ciertas consideraciones te6ricas sobre el numero de condici6n de una matriz A las cuales pueden ser consultadas en Burden 1985 paginas 481 y 482 cuando se trabaja en
aritmetica finita (decimal) con redondeo a t-digitos y Cond(A) ~ 10t se espera un mal
112 METOOOS NUMERICOS
comportamiento de A con respecto a la solucion de AX = b Y A se considera mal
condicionada En este ejemplo Cond(A) = 123 lt 103 )
Ahora
03 589 ) (-100) ( 592)AX-b e ( 531 - 610 101 470
Estas operaciones se realizari en doble preci sion (6 digitos)
59189) ( 592) ( - 011 )(= -59261 - 470 = -106261
(Para evitar la perdida de cifras sign ificativas se debe calcular el vector erro r res idual
R = AX - b en doble precision)
Conv irtiendo este ultimo resultado a tres digltos usando redondeo se obtiene
R = (-011) - 106
Entonces II R II =106 y II b II = 592 Y por tanto
II R la 1 106 1 11gtlt - X II 106 II R II en
-II-11- ( )= - - =145 ~ II II ~ 2202= 123- = Cond (A)-II - lI shyb co Cond A 592 123 X 592 b DO
II R pero como Cond (A) ~ 1 entonces se espera que fbiC Y X- X sean mas 0 menosII xL
II X- X L sea tamble ndel mismo tamano y ya que es grande se espera que II X II
grande +
En situaciones como a observada en este ejemplo se sugiere hacer un refinamiento
iterativo sobre la soluci6n caculada X 0 usar esta soluci6n calcuada como aproximaci6n inicial en un metodo iterativ~ con et proposito de tratar de mejorar la soluci6n aproximada y ograr que el error residua l relativo sea mas pequeno EI metodo de reflna81iento iterativo puede ser consutado en Kincaid 1972 pag inas 174-176
La soucion exacta del sistema en consideracion es X = (J= ( 0) A que se debe a
diferencia entre la soucion exacta y a solucion calculada
Observe que e error en e calculo de )(2 con respecto a x2 fue de solo 01 (un error relativo
de 1) Y este error fue multiplicado par un factor de aproximadamente - 2000 al obtener
)(1 debido a orden en que se realiz6 la el iminaci6n Gaussiana
Instruccion en DERIVE
RESUELVA_1(Ab) Simplifica en a SOUCI6n entra como un vector fila 0
En el ejemplo anterior la elimina~i6n sistema de ecuaciones lineales ble~ nII1II
del algoritmo de eliminaci6n Gausslana pequeno) Hay sin embargo SIlUC -
Gaussiana es numencamente estable consultarse en Burden 1985 paginas366 y
Teorema 35 Si A =(ajiLn (EDD) par filas es decir si
entonces A es invertible (no-singular intercambio de filas en cualquier calculos son estables con respeclO
Notese que como consecuencia
por filas entonces A t~ene unos en su diagonal pnnClpalyU
~~a~ (vease Burden 19a5
satisface
Capitulo 3 SOLUCION NUMERIC A DE SISTEMAS DE ECUACIONES 113
Instrucci6n en DERIVE
RESUELVA_1(Ab) Simplifica en la soluci6n exacta X del sistema AX=b EI vector b se entra como un vector fila 0
En el ejemplo anterior la eliminaci6n Gaussiana condujo a una respuesta defectuosa de un sistema de ecuaciones lineales bien condicionado Esto muestra la inestabilidad numerica del algoritmo de eliminaci6n Gaussiana (consecuencia de la divisi6n por un nOmero (pivote) pequero) Hay sin embargo situaciones en las cuales el algoritmo de eliminaci6n Gaussiana es numericamente estable EI siguiente teorema cuya demostraci6n puede consultarse en Burden 1985 paginas 366 y 367 se refiere a una de tales situaciones
Teorema 35 Si A = (ai j )n xn es una matriz estrictamente dominante diagonalmente
(EDD) por filas es decir si
n Iaii Igt II ai j I para cada i == 12 n
)= 1 ) 1
entonces A es invertible (no-singular) Ademas se puede rea lizar eliminaci6n Gaussiana sin intercambio de filas en cualquier sistema AX = b para obtener su Onica soluci6n y los calculos son estables con respecto al crecimiento de los errores de redondeo V
N6tese que como consecuencia del teorema anterior se tiene que Si A = (a l ) )n xn es EDD
por filas entonces A tiene factorizaci6n LU es decir A == LU con L triangular inferior con unos en su diagonal principal y U triangular superior (escalonada)
sea tambien
Observe que la matriz de coeficientes del ejemplo 33 anterior no es EDD por filas
sectJeor9~~am9i~ -valioo=paramatrj simetricas
(vease Bllrden 19 5 ~ina 368) Una matriz A E Rnxn ~ i~etrica se dice definida
positiva si satisface una cualquiera de las sig uientes condiciones (Ias- cualesson equival~tes) a _-----
- - -- shy
i) ~AX gt 0 paralodo 2CERn X 0 J
iiXTodos los valores propios de A son pos1tlvOS7
iii) Todos los pivotes obtenidos en la eliminaci6n Gaussiana sobre A sin intercambio de filas - ~on positiv~s
iv) Todas las submatrices principales de A tienen determinante positiv~
(Las submatrices principales de la matriz A == (ai J) son las matrices - n ~ n
1k a ] k == 12bull f) )
akk v
114 METODOS NUMERICOS
N6tese nuevamente que Si A E R nxn es simetrica y definida positiva entonces A tiene
factorizaci6n A = LU con L triangular inferior con sus componentes sabre la diagonal principal iguales a uno y U triangular superior (ver mas adelante factorizaci6n de Choleski)
Observe que la matriz de coeficientes de l ejemplo 33 anterior no es simetrica
34 ESTRATEGIAS DE PIVOTEO
EI eJemplo 33 anterior muestra una de las dificultades que pueden surgir al aplicar el metodo
de eliminaci6n Gaussiana cuando el pivote ajj-1) es pequeno comparado con algunos
(j-1) e ementos I a t para J ~ I t ~ n
Para tratar de evitar tales dificultades se introduce en el metodo de eliminaci6n Gaussiana una estrategia IIamada de pivoteo la cual consiste en seleccionar el pivote de acuerdo con un cierto criter io Nosotros usaremos dos estrategias la estrategia de pivoteo maximo por columna 0 pivoteo parcial y la estrategia de pivoteo escalado de fila 0 escalamiento
341 Pivoteo maximo por columna 0 pivoteo parcial Esta estrategia difiere de
eliminaci6n Gaussiana simple unicamente en la escogencia del pivote ajj-1) la cual se hace
ahora as
Para j = 12 n - 1 se determina el menor entero k J ~ k ~ n tal que
y
es decir seleccionamos el primer elemento diferente de cero sobre la columna j-esima a
partir de la j-esima fila y que tenga mayor valor absoluto (para j = 1 a~j j1 ) = al~) ak 1 )
Si tal k no existe el sistema no tiene soluci6n unica y el proceso se puede terminar
Si tal k existe y k ot j entonces hacemos intercambio de las ecuaciones j-esima y k-esima
y continuamos con la eliminaci6n Gaussiana l
Ilustremos esta estrategia para resolver el sistema
E1 03X1 + 589x2 = 592 E2 531x1 - 61 OX 2 = 470
entonces k = 2 ot 1= j asl que
Por sustituci6n regresiva
Observe que en este caSO
SWAP(A i j)
que es el mismo del ejemplo 33 usando aritmetica con redondeo a tres digitos
Como para j = 1 se tiene que
Capitulo 3 SOLUCION NUMERICA DE SISTEMAS DE ECUACIONES 115
Max Iall I Ia21 I = Max 03 531 = 531 = Ia21 I 0
entonces k = 2 1= j as i que intercambiamos El y E2 Y continuamos con la el imi nac i6n
531 - 610 (A b)= ( 03 589 592) 470) (531 - 610 470 03 589 592
E2I (~) ( m21=5 65 ~ 1O - l) 531 - 610 5 31 470) ------------------~) ( o 589 589
Por sustituci6n regresiva obtenemos
x = 589 = 100 Xl = 470 +610(100) = 531 = 100 2 589 531 531
Observe que en este caso X= (J es la soluci6n exacta del sistema dado bull
Instrucci6n en DERIVE
SWAP(A i j) Intercambia las filas (0 elementos) i y j de la matriz A (de un vector) 0
Nota En el procedimiento de pivoteo maximo por columna (p ivoteo parcial) cada multiplicador mi j es tal que
y aunque esta estrategia permite resolver satisfactoriamente muchos sistemas de ecuaciones lineales hay casos donde fracasa como se ilustra en el siguiente ejemplo
Ejemplo 34 Consideremos el sistema
El 300Xl + 58900x2 59200 E2 531x1 -610X2 - 470
el cual es un sistema equivalente al del ejemplo 33 (los coeficientes de la primera ecuaci6n
en el sistema del ejemplo 33 han sido multiplicados por 10 3) EI pivoteo maximo por
columna con aritmetica de redondeo a tres dig itos nos lIeva a los siguientes resultados
59200) E2 1 (-~pound) (m 2 1 =177) ) ( 300 58900 59200 ) (A b) = ( 300 58900 531 - 610 470 o - 10400 -10500
y por sustituci6n regresiva x2 = 101 Y Xl = - 100 que es la misma solucion que se obtiene si ~ usamos eliminaci6n Gaussiana simple
116 METODOS NUMERICOS
En casas como el de este ejemplo donde un pivote es mucho mas pequeno que alguno de los coeficientes de la ecuacion que el encabeza se recomienda la tecnica conocida como pivoteo escalado de fila 0 escalamiento la cual es nuestra segunda estrategia
342 Pivoteo escalado de fila Esta tecnica solo difiere de la eliminaci6n Gaussiana simple al igual que el pivoteo parcial en la escogencia del pivote
lEsta vez el pivote alr ) se escoge como se indica a continuacion
Para j = 12 n - 1 hacemos 10 siguiente
a) Para i = j j + 1 n calculamos
S = Max Ia(l- l) I Factor de escala I j s i s n II
Si Sj = 0 entonces el sistema no tiene solucion unica y el proceso se puede terminar
b) Para i = jj+1 n calculamos
c) Encontramos el menor entero k con j o k 0 n tal que
Ia~jjl ) I IaH- l) I -------- = Max --------
Sk j S j $ n SI
Si tal k no existe entonces el sistema no tiene solucion unica y el proceso se puede terminar
Si tal k existe y kF j entonces hacemos intercambio de las ecuaciones j-esima y kshyesima
y continuamos con la eliminacion Gaussiana V
Apliquemos esta estrategia para resolver el sistema del ejemplo 34 usando aritmetica con redondeo a tres dfgitos
Para j = 1
a) Sl = Max Iall I I a12 I = Max 300 58900 = 58900F 0 Y
Max 53~ 610 IshyS2 =Max a21 a n -shy
b) Ahora ~= 51
la0
l 300
S2
531 al l I a21 -Max - 0 -
c) Max -S -S- - 5890 NIl 1 2
las ecuaciolllltanto intercamblamos
Gaussiana
300 58900 (A b) = ( 531 -610
Por sustitucion regresiva
Xl tOO
Observe que en este GaSO X
35 FACTORIZACION
Capitulo 3 SOLUCION NUMERICA DE SISTEMAS DE ECUACIONES 117
8 2 =Max I a21 I I a22 1 =Max 531 610 = 610 0
b) Ahora
Ia I 300 - 8- = 58900 Y
1
Ia 2 I 531 - - = - shy
6108 2
c) Max ~ ~ I=Max ~ 531 = 531 =~ 0 as que k = 2 1= j y por8 8 2 58900 60 60 8 2
tanto intercambiamos las ecuaciones E y E2 Y continuamos con la eliminacion
Gaussiana
(A b) = (300 58900 59200) p ( 531 -610 470 ) 2 )
531 - 610 470 300 58900 59200
E2 (300) (m2 =565 ) ( 531 -6 1 0531 47 0 )
--------------~) o 58900 58900
Por sustituci6n regresiva
_ 470 +610(100) 531X2 = 100 Xl = =-- = 100
531 531
Observe que en este caso X= ( ~ J = ( 100) es la solucion exacta del sistema bull 100x2
35 FACTORIZACION TRIANGULAR
Consideremos un sistema AX = b con A no-singular y b 0 Con respecto a la matriz A se sabe que existen matrices P de permutaci6n L triangular inferior con sus componentes sobre la diagonal principal iguales a uno y U triangular superior (escalon ada) tales que PA = LU Una forma eficiente computacionalmente de encontrar P L Y U usando eliminaci6n Gaussiana se muestra en el siguiente ejemplo Ejemplo 35 Resuelva el siguiente sistema de ecuaciones lineales usando eliminacion Gaussiana con pivoteo parcial y obtenga una factorizaci6n PA =LU para la matriz A de coeficientes asociada con este metodo
- x + 2X2 - 3X3 = -2 3x1 - 3x2 - X 3 = -4
x + x2 = 3IEmpezamos introduciendo un vector p =(PP2P3) T el cual inicializamos_con P = i i = 123 Y
donde se almacenaran los intercambios necesarios en el proceso de eliminacion Gaussiana
118 METODOS NUMERICOS
con pivoteo parcial (el numero de componentes del vector p coincide con el orden n del sistema a resolver)
- 1 2 -3 -2] -4] - 2 P = (123) T (Ab) = ~ - 3 -1 - 4
[ o 3 3
3 - 3 - 1 - 4(m=-H (m3l = ~ )
E( - ~) E3(~) 10 10(-i) 3 3
2(i) 1 13 3 3
(Observe que cada multiplicador mij es almacenado en la posici6n correspondiente (i j) en
la matriz de trabajo)
3 - 3 -1 - 4
Max 121 =2 p=( 23 1) T bull (i) 2 1
-
3 13 3
(-i) 10 3
10 3
(Observe que la permutaci6n se hace para las filas 2 y 3 completas es decir incluyendo los multiplicadores)
3 - 3 - 1 --4(mJ2 = ~ ) E3( ~)
2(i) 13 3 3
(- i) (i) 7 11 2 2
La eficiencia en el metoda indicado se debe a que en la misma matriz de trabajo se almacenan los multiplicadores que van a conformar la matriz L (en el ejemplo son los numeros que se encuentran dentro de parentesis) 10 que significa un ahorro de memoria y como los intercambios necesarios afectan simultaneamente a las matrices L y U se evita tener que volver a la matriz original a realizar los intercambios ya observados y repetir la eliminaci6n Gaussiana con pivoteo parcial De esta manera al terminar el proceso de eliminaci6n podemos leer en la matriz final la parte estrictamente triangular inferior de L (son los numeros entre parentesis) y la matriz triangular superior U (que es la parte triangular superior de la matriz final) y en el vector p final quedan almacenados los intercambios realizados que se usan para producir la matriz de permutaci6n P
Para el ejemplo 35
3 -31 0 0 1 0 2U=0L= 3 1 1 0 0-3 2
(Verifique que PA =LU)
Para obtener la soluci6n del sistema
reducido
yobtenemos
C6mo se resuelve el
Capitulo 3 SOLUCION NUMERICA DE SISTEMAS DE ECUACIONES 119
Para el ejemplo 35
1 0 0 3 -3 - 1 1 1 1
L= - 0 U= 0 2 - Y como P~ (231) entonces [~ 0 lp3 3 1 1 7 0
-- - 0 0 3 2 2
(Verifique que PA = LU )
Para obtener la soluci6n del sistema original usamos sustituci6n regresiva en el sistema reducido
-4 13
UX = 3
11
2 ~
c y obtenemos
asi que la soluci6n (exacta) del sistema dado es X = (23 40 ~)T21 21 7
C6mo se resuelve el sistema AX = b a partir de la factorizaci6n PA = LU obtenida bull
351 Algunas aplicaciones de la factorizaci6n PA = LU La factorizaci6n PA = LU es utilizada eficientemente en aquellos casos donde se trabaja repetfdamente con la misma matriz A Dos de esos casos se presentan a continuaci6n
1) Resolver varios sistemas AX =b con la misma matriz de coeficientes A ya que en p L Y U esta almacenado todo el proceso de eliminaci6n Gaussiana EI algoritmo se bas a en la siguiente equivalencia
AX = b lt=gt PAX = Pb lt=gt LUX = Pb
lt=gt UX = L- 1Pb
lt=gtUXoo c y c oo L- 1Pb
lt=gt UX == C Y Lc = Pb lt=gt Lc = Pb Y UX = c
Los pasos a seguir son
Paso 1 Calcular Pb Paso 2 Resolver para c Lc = Pb por sustituci6n progresiva Paso 3 Resolver para X UX = c por sustituci6n regresiva
Como ejercicio resuelva el sistema del ejemplo anterior usando este algoritmo
120 METODOS NUMERICOS
2) Encontrar la matriz inversa A -1 de una matriz invertible A resolviendo los n-sistemas
AX =eU) j = 12 n
donde e(1) =(O O1 o 0f ERn
t posicion j
La solucion X del sistema AX = e(n j = 12 n produce la correspondiente columna j shy
esima de la matriz A -1
Como ejercicio compare el numero de operaciones para encontrar A - usando el metodo de Gauss-Jordan con el numero de operaciones resolviendo los n-sistemas indicados antes V
Ejercicio 33 Calcule la inversa de la matriz A de coeficientes del sistema del ejemplo 35 usando el metodo de Gauss-Jordan y tambien usando la factorizacion PA = LU bull
36 SISTEMAS TRIDIAGONALES
Un caso muy importante de sistemas de ecuaciones lineales que requiere un tratamiento especial es el de los sistemas tridiagonales Tales sistemas aparecen en diversas aplicaciones como por ejemplo al utilizar metodos de diferencias flnitas en la solucion de problemas con valores en la frontera para ecuaciones diferenciales ordinarias y como veremos mas adelante en el problema de la interpolacion segmentaria cubica
Un sistema tridiagonal es de la forma
dx1 + C1X2
a 2 x + d2 x2 + C 2 X 3
a3 x2 + d3 x3 + C 3 X 4
an-1Xn- 2 + dn_ xn_ + cn_xn = bn_1
anxn_1 + dnxn = bn
La matriz de coeficientes del sistema es
C 0 0 od1
0a 2 d2 c2
o a3 d3 c 3
A =
o
la cual se dice una matriz TRIDIAGONAL
En general
ij=12 n
suponiendo que la matriz A de podemos usar eliminaclon Gausslana eliminacion Gaussiana simple es adecuadl EDD por filas) Otra forma de resolver
EDD por filas es a partir de la componentes de las matrices L Y U as triangular inferior con todos sus superior Para encontrar tales malnce5
o 0 0
Y2 o o o o
L =
o Yn t
o o
Como Cl Ct
Y2 Cl Y2C +az
0 h a2
LU =
0
0
d
=
v
(311 )
Capitulo 3 SOlUCION NUMERICA DE SISTEMAS DE ECUACIONES 111
(A b) = ( 03 589 592) __( 5__) _m 2 _=1_77 ) ~) ( 03 589E2_ ~~ _( _ _ 592 ) 531 -610 470 0 - 10400 -10500
y por sustituci6n regresiva
x = -10500 =101 2 - 10400
)(1 = 592 - 589(101) = 592 - 595 = - 3 = - 100 D3 D3 D3
- (X1 J (-1 00)luego la soluci6n calculada es X = _ = 101x2
Instrucci6n en DERIVE
PIVOT(A i j) Usa operaciones elemetales de fila para Simplificar (0 aproXimar) en una matriz obtenida de la matriz A que tiene ceros en la columna j y por debajo de la fila i 0
Que puede decir de la calidad de la soluci6n aproximada X
IIX -XIIPara intentar responder esta pregunta encontremos las cotas para el error relativo
IIXII dadas por el teorema 33
03 589)Como A = entonces usando aritmetica con redondeo a tres digitos para todos( 531 - 610
los calculos se obtienen las siguientes aproximaciones
- 1 1 (-610 - 589)A =-shy-313 - 531 03
II A II = Max 589 114 = 589
A - 1 II = _ 1_ Max 650 534 = _ 1_ 650 = 208II 313 313
entonces Cond (A) = II A II II A -111 = (589)(208) = 123 que no es muy grande
comparado con uno asi que la matriz A puede considerarse bien condicionada
(Por ciertas consideraciones te6ricas sobre el numero de condici6n de una matriz A las cuales pueden ser consultadas en Burden 1985 paginas 481 y 482 cuando se trabaja en
aritmetica finita (decimal) con redondeo a t-digitos y Cond(A) ~ 10t se espera un mal
112 METOOOS NUMERICOS
comportamiento de A con respecto a la solucion de AX = b Y A se considera mal
condicionada En este ejemplo Cond(A) = 123 lt 103 )
Ahora
03 589 ) (-100) ( 592)AX-b e ( 531 - 610 101 470
Estas operaciones se realizari en doble preci sion (6 digitos)
59189) ( 592) ( - 011 )(= -59261 - 470 = -106261
(Para evitar la perdida de cifras sign ificativas se debe calcular el vector erro r res idual
R = AX - b en doble precision)
Conv irtiendo este ultimo resultado a tres digltos usando redondeo se obtiene
R = (-011) - 106
Entonces II R II =106 y II b II = 592 Y por tanto
II R la 1 106 1 11gtlt - X II 106 II R II en
-II-11- ( )= - - =145 ~ II II ~ 2202= 123- = Cond (A)-II - lI shyb co Cond A 592 123 X 592 b DO
II R pero como Cond (A) ~ 1 entonces se espera que fbiC Y X- X sean mas 0 menosII xL
II X- X L sea tamble ndel mismo tamano y ya que es grande se espera que II X II
grande +
En situaciones como a observada en este ejemplo se sugiere hacer un refinamiento
iterativo sobre la soluci6n caculada X 0 usar esta soluci6n calcuada como aproximaci6n inicial en un metodo iterativ~ con et proposito de tratar de mejorar la soluci6n aproximada y ograr que el error residua l relativo sea mas pequeno EI metodo de reflna81iento iterativo puede ser consutado en Kincaid 1972 pag inas 174-176
La soucion exacta del sistema en consideracion es X = (J= ( 0) A que se debe a
diferencia entre la soucion exacta y a solucion calculada
Observe que e error en e calculo de )(2 con respecto a x2 fue de solo 01 (un error relativo
de 1) Y este error fue multiplicado par un factor de aproximadamente - 2000 al obtener
)(1 debido a orden en que se realiz6 la el iminaci6n Gaussiana
Instruccion en DERIVE
RESUELVA_1(Ab) Simplifica en a SOUCI6n entra como un vector fila 0
En el ejemplo anterior la elimina~i6n sistema de ecuaciones lineales ble~ nII1II
del algoritmo de eliminaci6n Gausslana pequeno) Hay sin embargo SIlUC -
Gaussiana es numencamente estable consultarse en Burden 1985 paginas366 y
Teorema 35 Si A =(ajiLn (EDD) par filas es decir si
entonces A es invertible (no-singular intercambio de filas en cualquier calculos son estables con respeclO
Notese que como consecuencia
por filas entonces A t~ene unos en su diagonal pnnClpalyU
~~a~ (vease Burden 19a5
satisface
Capitulo 3 SOLUCION NUMERIC A DE SISTEMAS DE ECUACIONES 113
Instrucci6n en DERIVE
RESUELVA_1(Ab) Simplifica en la soluci6n exacta X del sistema AX=b EI vector b se entra como un vector fila 0
En el ejemplo anterior la eliminaci6n Gaussiana condujo a una respuesta defectuosa de un sistema de ecuaciones lineales bien condicionado Esto muestra la inestabilidad numerica del algoritmo de eliminaci6n Gaussiana (consecuencia de la divisi6n por un nOmero (pivote) pequero) Hay sin embargo situaciones en las cuales el algoritmo de eliminaci6n Gaussiana es numericamente estable EI siguiente teorema cuya demostraci6n puede consultarse en Burden 1985 paginas 366 y 367 se refiere a una de tales situaciones
Teorema 35 Si A = (ai j )n xn es una matriz estrictamente dominante diagonalmente
(EDD) por filas es decir si
n Iaii Igt II ai j I para cada i == 12 n
)= 1 ) 1
entonces A es invertible (no-singular) Ademas se puede rea lizar eliminaci6n Gaussiana sin intercambio de filas en cualquier sistema AX = b para obtener su Onica soluci6n y los calculos son estables con respecto al crecimiento de los errores de redondeo V
N6tese que como consecuencia del teorema anterior se tiene que Si A = (a l ) )n xn es EDD
por filas entonces A tiene factorizaci6n LU es decir A == LU con L triangular inferior con unos en su diagonal principal y U triangular superior (escalonada)
sea tambien
Observe que la matriz de coeficientes del ejemplo 33 anterior no es EDD por filas
sectJeor9~~am9i~ -valioo=paramatrj simetricas
(vease Bllrden 19 5 ~ina 368) Una matriz A E Rnxn ~ i~etrica se dice definida
positiva si satisface una cualquiera de las sig uientes condiciones (Ias- cualesson equival~tes) a _-----
- - -- shy
i) ~AX gt 0 paralodo 2CERn X 0 J
iiXTodos los valores propios de A son pos1tlvOS7
iii) Todos los pivotes obtenidos en la eliminaci6n Gaussiana sobre A sin intercambio de filas - ~on positiv~s
iv) Todas las submatrices principales de A tienen determinante positiv~
(Las submatrices principales de la matriz A == (ai J) son las matrices - n ~ n
1k a ] k == 12bull f) )
akk v
114 METODOS NUMERICOS
N6tese nuevamente que Si A E R nxn es simetrica y definida positiva entonces A tiene
factorizaci6n A = LU con L triangular inferior con sus componentes sabre la diagonal principal iguales a uno y U triangular superior (ver mas adelante factorizaci6n de Choleski)
Observe que la matriz de coeficientes de l ejemplo 33 anterior no es simetrica
34 ESTRATEGIAS DE PIVOTEO
EI eJemplo 33 anterior muestra una de las dificultades que pueden surgir al aplicar el metodo
de eliminaci6n Gaussiana cuando el pivote ajj-1) es pequeno comparado con algunos
(j-1) e ementos I a t para J ~ I t ~ n
Para tratar de evitar tales dificultades se introduce en el metodo de eliminaci6n Gaussiana una estrategia IIamada de pivoteo la cual consiste en seleccionar el pivote de acuerdo con un cierto criter io Nosotros usaremos dos estrategias la estrategia de pivoteo maximo por columna 0 pivoteo parcial y la estrategia de pivoteo escalado de fila 0 escalamiento
341 Pivoteo maximo por columna 0 pivoteo parcial Esta estrategia difiere de
eliminaci6n Gaussiana simple unicamente en la escogencia del pivote ajj-1) la cual se hace
ahora as
Para j = 12 n - 1 se determina el menor entero k J ~ k ~ n tal que
y
es decir seleccionamos el primer elemento diferente de cero sobre la columna j-esima a
partir de la j-esima fila y que tenga mayor valor absoluto (para j = 1 a~j j1 ) = al~) ak 1 )
Si tal k no existe el sistema no tiene soluci6n unica y el proceso se puede terminar
Si tal k existe y k ot j entonces hacemos intercambio de las ecuaciones j-esima y k-esima
y continuamos con la eliminaci6n Gaussiana l
Ilustremos esta estrategia para resolver el sistema
E1 03X1 + 589x2 = 592 E2 531x1 - 61 OX 2 = 470
entonces k = 2 ot 1= j asl que
Por sustituci6n regresiva
Observe que en este caSO
SWAP(A i j)
que es el mismo del ejemplo 33 usando aritmetica con redondeo a tres digitos
Como para j = 1 se tiene que
Capitulo 3 SOLUCION NUMERICA DE SISTEMAS DE ECUACIONES 115
Max Iall I Ia21 I = Max 03 531 = 531 = Ia21 I 0
entonces k = 2 1= j as i que intercambiamos El y E2 Y continuamos con la el imi nac i6n
531 - 610 (A b)= ( 03 589 592) 470) (531 - 610 470 03 589 592
E2I (~) ( m21=5 65 ~ 1O - l) 531 - 610 5 31 470) ------------------~) ( o 589 589
Por sustituci6n regresiva obtenemos
x = 589 = 100 Xl = 470 +610(100) = 531 = 100 2 589 531 531
Observe que en este caso X= (J es la soluci6n exacta del sistema dado bull
Instrucci6n en DERIVE
SWAP(A i j) Intercambia las filas (0 elementos) i y j de la matriz A (de un vector) 0
Nota En el procedimiento de pivoteo maximo por columna (p ivoteo parcial) cada multiplicador mi j es tal que
y aunque esta estrategia permite resolver satisfactoriamente muchos sistemas de ecuaciones lineales hay casos donde fracasa como se ilustra en el siguiente ejemplo
Ejemplo 34 Consideremos el sistema
El 300Xl + 58900x2 59200 E2 531x1 -610X2 - 470
el cual es un sistema equivalente al del ejemplo 33 (los coeficientes de la primera ecuaci6n
en el sistema del ejemplo 33 han sido multiplicados por 10 3) EI pivoteo maximo por
columna con aritmetica de redondeo a tres dig itos nos lIeva a los siguientes resultados
59200) E2 1 (-~pound) (m 2 1 =177) ) ( 300 58900 59200 ) (A b) = ( 300 58900 531 - 610 470 o - 10400 -10500
y por sustituci6n regresiva x2 = 101 Y Xl = - 100 que es la misma solucion que se obtiene si ~ usamos eliminaci6n Gaussiana simple
116 METODOS NUMERICOS
En casas como el de este ejemplo donde un pivote es mucho mas pequeno que alguno de los coeficientes de la ecuacion que el encabeza se recomienda la tecnica conocida como pivoteo escalado de fila 0 escalamiento la cual es nuestra segunda estrategia
342 Pivoteo escalado de fila Esta tecnica solo difiere de la eliminaci6n Gaussiana simple al igual que el pivoteo parcial en la escogencia del pivote
lEsta vez el pivote alr ) se escoge como se indica a continuacion
Para j = 12 n - 1 hacemos 10 siguiente
a) Para i = j j + 1 n calculamos
S = Max Ia(l- l) I Factor de escala I j s i s n II
Si Sj = 0 entonces el sistema no tiene solucion unica y el proceso se puede terminar
b) Para i = jj+1 n calculamos
c) Encontramos el menor entero k con j o k 0 n tal que
Ia~jjl ) I IaH- l) I -------- = Max --------
Sk j S j $ n SI
Si tal k no existe entonces el sistema no tiene solucion unica y el proceso se puede terminar
Si tal k existe y kF j entonces hacemos intercambio de las ecuaciones j-esima y kshyesima
y continuamos con la eliminacion Gaussiana V
Apliquemos esta estrategia para resolver el sistema del ejemplo 34 usando aritmetica con redondeo a tres dfgitos
Para j = 1
a) Sl = Max Iall I I a12 I = Max 300 58900 = 58900F 0 Y
Max 53~ 610 IshyS2 =Max a21 a n -shy
b) Ahora ~= 51
la0
l 300
S2
531 al l I a21 -Max - 0 -
c) Max -S -S- - 5890 NIl 1 2
las ecuaciolllltanto intercamblamos
Gaussiana
300 58900 (A b) = ( 531 -610
Por sustitucion regresiva
Xl tOO
Observe que en este GaSO X
35 FACTORIZACION
Capitulo 3 SOLUCION NUMERICA DE SISTEMAS DE ECUACIONES 117
8 2 =Max I a21 I I a22 1 =Max 531 610 = 610 0
b) Ahora
Ia I 300 - 8- = 58900 Y
1
Ia 2 I 531 - - = - shy
6108 2
c) Max ~ ~ I=Max ~ 531 = 531 =~ 0 as que k = 2 1= j y por8 8 2 58900 60 60 8 2
tanto intercambiamos las ecuaciones E y E2 Y continuamos con la eliminacion
Gaussiana
(A b) = (300 58900 59200) p ( 531 -610 470 ) 2 )
531 - 610 470 300 58900 59200
E2 (300) (m2 =565 ) ( 531 -6 1 0531 47 0 )
--------------~) o 58900 58900
Por sustituci6n regresiva
_ 470 +610(100) 531X2 = 100 Xl = =-- = 100
531 531
Observe que en este caso X= ( ~ J = ( 100) es la solucion exacta del sistema bull 100x2
35 FACTORIZACION TRIANGULAR
Consideremos un sistema AX = b con A no-singular y b 0 Con respecto a la matriz A se sabe que existen matrices P de permutaci6n L triangular inferior con sus componentes sobre la diagonal principal iguales a uno y U triangular superior (escalon ada) tales que PA = LU Una forma eficiente computacionalmente de encontrar P L Y U usando eliminaci6n Gaussiana se muestra en el siguiente ejemplo Ejemplo 35 Resuelva el siguiente sistema de ecuaciones lineales usando eliminacion Gaussiana con pivoteo parcial y obtenga una factorizaci6n PA =LU para la matriz A de coeficientes asociada con este metodo
- x + 2X2 - 3X3 = -2 3x1 - 3x2 - X 3 = -4
x + x2 = 3IEmpezamos introduciendo un vector p =(PP2P3) T el cual inicializamos_con P = i i = 123 Y
donde se almacenaran los intercambios necesarios en el proceso de eliminacion Gaussiana
118 METODOS NUMERICOS
con pivoteo parcial (el numero de componentes del vector p coincide con el orden n del sistema a resolver)
- 1 2 -3 -2] -4] - 2 P = (123) T (Ab) = ~ - 3 -1 - 4
[ o 3 3
3 - 3 - 1 - 4(m=-H (m3l = ~ )
E( - ~) E3(~) 10 10(-i) 3 3
2(i) 1 13 3 3
(Observe que cada multiplicador mij es almacenado en la posici6n correspondiente (i j) en
la matriz de trabajo)
3 - 3 -1 - 4
Max 121 =2 p=( 23 1) T bull (i) 2 1
-
3 13 3
(-i) 10 3
10 3
(Observe que la permutaci6n se hace para las filas 2 y 3 completas es decir incluyendo los multiplicadores)
3 - 3 - 1 --4(mJ2 = ~ ) E3( ~)
2(i) 13 3 3
(- i) (i) 7 11 2 2
La eficiencia en el metoda indicado se debe a que en la misma matriz de trabajo se almacenan los multiplicadores que van a conformar la matriz L (en el ejemplo son los numeros que se encuentran dentro de parentesis) 10 que significa un ahorro de memoria y como los intercambios necesarios afectan simultaneamente a las matrices L y U se evita tener que volver a la matriz original a realizar los intercambios ya observados y repetir la eliminaci6n Gaussiana con pivoteo parcial De esta manera al terminar el proceso de eliminaci6n podemos leer en la matriz final la parte estrictamente triangular inferior de L (son los numeros entre parentesis) y la matriz triangular superior U (que es la parte triangular superior de la matriz final) y en el vector p final quedan almacenados los intercambios realizados que se usan para producir la matriz de permutaci6n P
Para el ejemplo 35
3 -31 0 0 1 0 2U=0L= 3 1 1 0 0-3 2
(Verifique que PA =LU)
Para obtener la soluci6n del sistema
reducido
yobtenemos
C6mo se resuelve el
Capitulo 3 SOLUCION NUMERICA DE SISTEMAS DE ECUACIONES 119
Para el ejemplo 35
1 0 0 3 -3 - 1 1 1 1
L= - 0 U= 0 2 - Y como P~ (231) entonces [~ 0 lp3 3 1 1 7 0
-- - 0 0 3 2 2
(Verifique que PA = LU )
Para obtener la soluci6n del sistema original usamos sustituci6n regresiva en el sistema reducido
-4 13
UX = 3
11
2 ~
c y obtenemos
asi que la soluci6n (exacta) del sistema dado es X = (23 40 ~)T21 21 7
C6mo se resuelve el sistema AX = b a partir de la factorizaci6n PA = LU obtenida bull
351 Algunas aplicaciones de la factorizaci6n PA = LU La factorizaci6n PA = LU es utilizada eficientemente en aquellos casos donde se trabaja repetfdamente con la misma matriz A Dos de esos casos se presentan a continuaci6n
1) Resolver varios sistemas AX =b con la misma matriz de coeficientes A ya que en p L Y U esta almacenado todo el proceso de eliminaci6n Gaussiana EI algoritmo se bas a en la siguiente equivalencia
AX = b lt=gt PAX = Pb lt=gt LUX = Pb
lt=gt UX = L- 1Pb
lt=gtUXoo c y c oo L- 1Pb
lt=gt UX == C Y Lc = Pb lt=gt Lc = Pb Y UX = c
Los pasos a seguir son
Paso 1 Calcular Pb Paso 2 Resolver para c Lc = Pb por sustituci6n progresiva Paso 3 Resolver para X UX = c por sustituci6n regresiva
Como ejercicio resuelva el sistema del ejemplo anterior usando este algoritmo
120 METODOS NUMERICOS
2) Encontrar la matriz inversa A -1 de una matriz invertible A resolviendo los n-sistemas
AX =eU) j = 12 n
donde e(1) =(O O1 o 0f ERn
t posicion j
La solucion X del sistema AX = e(n j = 12 n produce la correspondiente columna j shy
esima de la matriz A -1
Como ejercicio compare el numero de operaciones para encontrar A - usando el metodo de Gauss-Jordan con el numero de operaciones resolviendo los n-sistemas indicados antes V
Ejercicio 33 Calcule la inversa de la matriz A de coeficientes del sistema del ejemplo 35 usando el metodo de Gauss-Jordan y tambien usando la factorizacion PA = LU bull
36 SISTEMAS TRIDIAGONALES
Un caso muy importante de sistemas de ecuaciones lineales que requiere un tratamiento especial es el de los sistemas tridiagonales Tales sistemas aparecen en diversas aplicaciones como por ejemplo al utilizar metodos de diferencias flnitas en la solucion de problemas con valores en la frontera para ecuaciones diferenciales ordinarias y como veremos mas adelante en el problema de la interpolacion segmentaria cubica
Un sistema tridiagonal es de la forma
dx1 + C1X2
a 2 x + d2 x2 + C 2 X 3
a3 x2 + d3 x3 + C 3 X 4
an-1Xn- 2 + dn_ xn_ + cn_xn = bn_1
anxn_1 + dnxn = bn
La matriz de coeficientes del sistema es
C 0 0 od1
0a 2 d2 c2
o a3 d3 c 3
A =
o
la cual se dice una matriz TRIDIAGONAL
En general
ij=12 n
suponiendo que la matriz A de podemos usar eliminaclon Gausslana eliminacion Gaussiana simple es adecuadl EDD por filas) Otra forma de resolver
EDD por filas es a partir de la componentes de las matrices L Y U as triangular inferior con todos sus superior Para encontrar tales malnce5
o 0 0
Y2 o o o o
L =
o Yn t
o o
Como Cl Ct
Y2 Cl Y2C +az
0 h a2
LU =
0
0
d
=
112 METOOOS NUMERICOS
comportamiento de A con respecto a la solucion de AX = b Y A se considera mal
condicionada En este ejemplo Cond(A) = 123 lt 103 )
Ahora
03 589 ) (-100) ( 592)AX-b e ( 531 - 610 101 470
Estas operaciones se realizari en doble preci sion (6 digitos)
59189) ( 592) ( - 011 )(= -59261 - 470 = -106261
(Para evitar la perdida de cifras sign ificativas se debe calcular el vector erro r res idual
R = AX - b en doble precision)
Conv irtiendo este ultimo resultado a tres digltos usando redondeo se obtiene
R = (-011) - 106
Entonces II R II =106 y II b II = 592 Y por tanto
II R la 1 106 1 11gtlt - X II 106 II R II en
-II-11- ( )= - - =145 ~ II II ~ 2202= 123- = Cond (A)-II - lI shyb co Cond A 592 123 X 592 b DO
II R pero como Cond (A) ~ 1 entonces se espera que fbiC Y X- X sean mas 0 menosII xL
II X- X L sea tamble ndel mismo tamano y ya que es grande se espera que II X II
grande +
En situaciones como a observada en este ejemplo se sugiere hacer un refinamiento
iterativo sobre la soluci6n caculada X 0 usar esta soluci6n calcuada como aproximaci6n inicial en un metodo iterativ~ con et proposito de tratar de mejorar la soluci6n aproximada y ograr que el error residua l relativo sea mas pequeno EI metodo de reflna81iento iterativo puede ser consutado en Kincaid 1972 pag inas 174-176
La soucion exacta del sistema en consideracion es X = (J= ( 0) A que se debe a
diferencia entre la soucion exacta y a solucion calculada
Observe que e error en e calculo de )(2 con respecto a x2 fue de solo 01 (un error relativo
de 1) Y este error fue multiplicado par un factor de aproximadamente - 2000 al obtener
)(1 debido a orden en que se realiz6 la el iminaci6n Gaussiana
Instruccion en DERIVE
RESUELVA_1(Ab) Simplifica en a SOUCI6n entra como un vector fila 0
En el ejemplo anterior la elimina~i6n sistema de ecuaciones lineales ble~ nII1II
del algoritmo de eliminaci6n Gausslana pequeno) Hay sin embargo SIlUC -
Gaussiana es numencamente estable consultarse en Burden 1985 paginas366 y
Teorema 35 Si A =(ajiLn (EDD) par filas es decir si
entonces A es invertible (no-singular intercambio de filas en cualquier calculos son estables con respeclO
Notese que como consecuencia
por filas entonces A t~ene unos en su diagonal pnnClpalyU
~~a~ (vease Burden 19a5
satisface
Capitulo 3 SOLUCION NUMERIC A DE SISTEMAS DE ECUACIONES 113
Instrucci6n en DERIVE
RESUELVA_1(Ab) Simplifica en la soluci6n exacta X del sistema AX=b EI vector b se entra como un vector fila 0
En el ejemplo anterior la eliminaci6n Gaussiana condujo a una respuesta defectuosa de un sistema de ecuaciones lineales bien condicionado Esto muestra la inestabilidad numerica del algoritmo de eliminaci6n Gaussiana (consecuencia de la divisi6n por un nOmero (pivote) pequero) Hay sin embargo situaciones en las cuales el algoritmo de eliminaci6n Gaussiana es numericamente estable EI siguiente teorema cuya demostraci6n puede consultarse en Burden 1985 paginas 366 y 367 se refiere a una de tales situaciones
Teorema 35 Si A = (ai j )n xn es una matriz estrictamente dominante diagonalmente
(EDD) por filas es decir si
n Iaii Igt II ai j I para cada i == 12 n
)= 1 ) 1
entonces A es invertible (no-singular) Ademas se puede rea lizar eliminaci6n Gaussiana sin intercambio de filas en cualquier sistema AX = b para obtener su Onica soluci6n y los calculos son estables con respecto al crecimiento de los errores de redondeo V
N6tese que como consecuencia del teorema anterior se tiene que Si A = (a l ) )n xn es EDD
por filas entonces A tiene factorizaci6n LU es decir A == LU con L triangular inferior con unos en su diagonal principal y U triangular superior (escalonada)
sea tambien
Observe que la matriz de coeficientes del ejemplo 33 anterior no es EDD por filas
sectJeor9~~am9i~ -valioo=paramatrj simetricas
(vease Bllrden 19 5 ~ina 368) Una matriz A E Rnxn ~ i~etrica se dice definida
positiva si satisface una cualquiera de las sig uientes condiciones (Ias- cualesson equival~tes) a _-----
- - -- shy
i) ~AX gt 0 paralodo 2CERn X 0 J
iiXTodos los valores propios de A son pos1tlvOS7
iii) Todos los pivotes obtenidos en la eliminaci6n Gaussiana sobre A sin intercambio de filas - ~on positiv~s
iv) Todas las submatrices principales de A tienen determinante positiv~
(Las submatrices principales de la matriz A == (ai J) son las matrices - n ~ n
1k a ] k == 12bull f) )
akk v
114 METODOS NUMERICOS
N6tese nuevamente que Si A E R nxn es simetrica y definida positiva entonces A tiene
factorizaci6n A = LU con L triangular inferior con sus componentes sabre la diagonal principal iguales a uno y U triangular superior (ver mas adelante factorizaci6n de Choleski)
Observe que la matriz de coeficientes de l ejemplo 33 anterior no es simetrica
34 ESTRATEGIAS DE PIVOTEO
EI eJemplo 33 anterior muestra una de las dificultades que pueden surgir al aplicar el metodo
de eliminaci6n Gaussiana cuando el pivote ajj-1) es pequeno comparado con algunos
(j-1) e ementos I a t para J ~ I t ~ n
Para tratar de evitar tales dificultades se introduce en el metodo de eliminaci6n Gaussiana una estrategia IIamada de pivoteo la cual consiste en seleccionar el pivote de acuerdo con un cierto criter io Nosotros usaremos dos estrategias la estrategia de pivoteo maximo por columna 0 pivoteo parcial y la estrategia de pivoteo escalado de fila 0 escalamiento
341 Pivoteo maximo por columna 0 pivoteo parcial Esta estrategia difiere de
eliminaci6n Gaussiana simple unicamente en la escogencia del pivote ajj-1) la cual se hace
ahora as
Para j = 12 n - 1 se determina el menor entero k J ~ k ~ n tal que
y
es decir seleccionamos el primer elemento diferente de cero sobre la columna j-esima a
partir de la j-esima fila y que tenga mayor valor absoluto (para j = 1 a~j j1 ) = al~) ak 1 )
Si tal k no existe el sistema no tiene soluci6n unica y el proceso se puede terminar
Si tal k existe y k ot j entonces hacemos intercambio de las ecuaciones j-esima y k-esima
y continuamos con la eliminaci6n Gaussiana l
Ilustremos esta estrategia para resolver el sistema
E1 03X1 + 589x2 = 592 E2 531x1 - 61 OX 2 = 470
entonces k = 2 ot 1= j asl que
Por sustituci6n regresiva
Observe que en este caSO
SWAP(A i j)
que es el mismo del ejemplo 33 usando aritmetica con redondeo a tres digitos
Como para j = 1 se tiene que
Capitulo 3 SOLUCION NUMERICA DE SISTEMAS DE ECUACIONES 115
Max Iall I Ia21 I = Max 03 531 = 531 = Ia21 I 0
entonces k = 2 1= j as i que intercambiamos El y E2 Y continuamos con la el imi nac i6n
531 - 610 (A b)= ( 03 589 592) 470) (531 - 610 470 03 589 592
E2I (~) ( m21=5 65 ~ 1O - l) 531 - 610 5 31 470) ------------------~) ( o 589 589
Por sustituci6n regresiva obtenemos
x = 589 = 100 Xl = 470 +610(100) = 531 = 100 2 589 531 531
Observe que en este caso X= (J es la soluci6n exacta del sistema dado bull
Instrucci6n en DERIVE
SWAP(A i j) Intercambia las filas (0 elementos) i y j de la matriz A (de un vector) 0
Nota En el procedimiento de pivoteo maximo por columna (p ivoteo parcial) cada multiplicador mi j es tal que
y aunque esta estrategia permite resolver satisfactoriamente muchos sistemas de ecuaciones lineales hay casos donde fracasa como se ilustra en el siguiente ejemplo
Ejemplo 34 Consideremos el sistema
El 300Xl + 58900x2 59200 E2 531x1 -610X2 - 470
el cual es un sistema equivalente al del ejemplo 33 (los coeficientes de la primera ecuaci6n
en el sistema del ejemplo 33 han sido multiplicados por 10 3) EI pivoteo maximo por
columna con aritmetica de redondeo a tres dig itos nos lIeva a los siguientes resultados
59200) E2 1 (-~pound) (m 2 1 =177) ) ( 300 58900 59200 ) (A b) = ( 300 58900 531 - 610 470 o - 10400 -10500
y por sustituci6n regresiva x2 = 101 Y Xl = - 100 que es la misma solucion que se obtiene si ~ usamos eliminaci6n Gaussiana simple
116 METODOS NUMERICOS
En casas como el de este ejemplo donde un pivote es mucho mas pequeno que alguno de los coeficientes de la ecuacion que el encabeza se recomienda la tecnica conocida como pivoteo escalado de fila 0 escalamiento la cual es nuestra segunda estrategia
342 Pivoteo escalado de fila Esta tecnica solo difiere de la eliminaci6n Gaussiana simple al igual que el pivoteo parcial en la escogencia del pivote
lEsta vez el pivote alr ) se escoge como se indica a continuacion
Para j = 12 n - 1 hacemos 10 siguiente
a) Para i = j j + 1 n calculamos
S = Max Ia(l- l) I Factor de escala I j s i s n II
Si Sj = 0 entonces el sistema no tiene solucion unica y el proceso se puede terminar
b) Para i = jj+1 n calculamos
c) Encontramos el menor entero k con j o k 0 n tal que
Ia~jjl ) I IaH- l) I -------- = Max --------
Sk j S j $ n SI
Si tal k no existe entonces el sistema no tiene solucion unica y el proceso se puede terminar
Si tal k existe y kF j entonces hacemos intercambio de las ecuaciones j-esima y kshyesima
y continuamos con la eliminacion Gaussiana V
Apliquemos esta estrategia para resolver el sistema del ejemplo 34 usando aritmetica con redondeo a tres dfgitos
Para j = 1
a) Sl = Max Iall I I a12 I = Max 300 58900 = 58900F 0 Y
Max 53~ 610 IshyS2 =Max a21 a n -shy
b) Ahora ~= 51
la0
l 300
S2
531 al l I a21 -Max - 0 -
c) Max -S -S- - 5890 NIl 1 2
las ecuaciolllltanto intercamblamos
Gaussiana
300 58900 (A b) = ( 531 -610
Por sustitucion regresiva
Xl tOO
Observe que en este GaSO X
35 FACTORIZACION
Capitulo 3 SOLUCION NUMERICA DE SISTEMAS DE ECUACIONES 117
8 2 =Max I a21 I I a22 1 =Max 531 610 = 610 0
b) Ahora
Ia I 300 - 8- = 58900 Y
1
Ia 2 I 531 - - = - shy
6108 2
c) Max ~ ~ I=Max ~ 531 = 531 =~ 0 as que k = 2 1= j y por8 8 2 58900 60 60 8 2
tanto intercambiamos las ecuaciones E y E2 Y continuamos con la eliminacion
Gaussiana
(A b) = (300 58900 59200) p ( 531 -610 470 ) 2 )
531 - 610 470 300 58900 59200
E2 (300) (m2 =565 ) ( 531 -6 1 0531 47 0 )
--------------~) o 58900 58900
Por sustituci6n regresiva
_ 470 +610(100) 531X2 = 100 Xl = =-- = 100
531 531
Observe que en este caso X= ( ~ J = ( 100) es la solucion exacta del sistema bull 100x2
35 FACTORIZACION TRIANGULAR
Consideremos un sistema AX = b con A no-singular y b 0 Con respecto a la matriz A se sabe que existen matrices P de permutaci6n L triangular inferior con sus componentes sobre la diagonal principal iguales a uno y U triangular superior (escalon ada) tales que PA = LU Una forma eficiente computacionalmente de encontrar P L Y U usando eliminaci6n Gaussiana se muestra en el siguiente ejemplo Ejemplo 35 Resuelva el siguiente sistema de ecuaciones lineales usando eliminacion Gaussiana con pivoteo parcial y obtenga una factorizaci6n PA =LU para la matriz A de coeficientes asociada con este metodo
- x + 2X2 - 3X3 = -2 3x1 - 3x2 - X 3 = -4
x + x2 = 3IEmpezamos introduciendo un vector p =(PP2P3) T el cual inicializamos_con P = i i = 123 Y
donde se almacenaran los intercambios necesarios en el proceso de eliminacion Gaussiana
118 METODOS NUMERICOS
con pivoteo parcial (el numero de componentes del vector p coincide con el orden n del sistema a resolver)
- 1 2 -3 -2] -4] - 2 P = (123) T (Ab) = ~ - 3 -1 - 4
[ o 3 3
3 - 3 - 1 - 4(m=-H (m3l = ~ )
E( - ~) E3(~) 10 10(-i) 3 3
2(i) 1 13 3 3
(Observe que cada multiplicador mij es almacenado en la posici6n correspondiente (i j) en
la matriz de trabajo)
3 - 3 -1 - 4
Max 121 =2 p=( 23 1) T bull (i) 2 1
-
3 13 3
(-i) 10 3
10 3
(Observe que la permutaci6n se hace para las filas 2 y 3 completas es decir incluyendo los multiplicadores)
3 - 3 - 1 --4(mJ2 = ~ ) E3( ~)
2(i) 13 3 3
(- i) (i) 7 11 2 2
La eficiencia en el metoda indicado se debe a que en la misma matriz de trabajo se almacenan los multiplicadores que van a conformar la matriz L (en el ejemplo son los numeros que se encuentran dentro de parentesis) 10 que significa un ahorro de memoria y como los intercambios necesarios afectan simultaneamente a las matrices L y U se evita tener que volver a la matriz original a realizar los intercambios ya observados y repetir la eliminaci6n Gaussiana con pivoteo parcial De esta manera al terminar el proceso de eliminaci6n podemos leer en la matriz final la parte estrictamente triangular inferior de L (son los numeros entre parentesis) y la matriz triangular superior U (que es la parte triangular superior de la matriz final) y en el vector p final quedan almacenados los intercambios realizados que se usan para producir la matriz de permutaci6n P
Para el ejemplo 35
3 -31 0 0 1 0 2U=0L= 3 1 1 0 0-3 2
(Verifique que PA =LU)
Para obtener la soluci6n del sistema
reducido
yobtenemos
C6mo se resuelve el
Capitulo 3 SOLUCION NUMERICA DE SISTEMAS DE ECUACIONES 119
Para el ejemplo 35
1 0 0 3 -3 - 1 1 1 1
L= - 0 U= 0 2 - Y como P~ (231) entonces [~ 0 lp3 3 1 1 7 0
-- - 0 0 3 2 2
(Verifique que PA = LU )
Para obtener la soluci6n del sistema original usamos sustituci6n regresiva en el sistema reducido
-4 13
UX = 3
11
2 ~
c y obtenemos
asi que la soluci6n (exacta) del sistema dado es X = (23 40 ~)T21 21 7
C6mo se resuelve el sistema AX = b a partir de la factorizaci6n PA = LU obtenida bull
351 Algunas aplicaciones de la factorizaci6n PA = LU La factorizaci6n PA = LU es utilizada eficientemente en aquellos casos donde se trabaja repetfdamente con la misma matriz A Dos de esos casos se presentan a continuaci6n
1) Resolver varios sistemas AX =b con la misma matriz de coeficientes A ya que en p L Y U esta almacenado todo el proceso de eliminaci6n Gaussiana EI algoritmo se bas a en la siguiente equivalencia
AX = b lt=gt PAX = Pb lt=gt LUX = Pb
lt=gt UX = L- 1Pb
lt=gtUXoo c y c oo L- 1Pb
lt=gt UX == C Y Lc = Pb lt=gt Lc = Pb Y UX = c
Los pasos a seguir son
Paso 1 Calcular Pb Paso 2 Resolver para c Lc = Pb por sustituci6n progresiva Paso 3 Resolver para X UX = c por sustituci6n regresiva
Como ejercicio resuelva el sistema del ejemplo anterior usando este algoritmo
120 METODOS NUMERICOS
2) Encontrar la matriz inversa A -1 de una matriz invertible A resolviendo los n-sistemas
AX =eU) j = 12 n
donde e(1) =(O O1 o 0f ERn
t posicion j
La solucion X del sistema AX = e(n j = 12 n produce la correspondiente columna j shy
esima de la matriz A -1
Como ejercicio compare el numero de operaciones para encontrar A - usando el metodo de Gauss-Jordan con el numero de operaciones resolviendo los n-sistemas indicados antes V
Ejercicio 33 Calcule la inversa de la matriz A de coeficientes del sistema del ejemplo 35 usando el metodo de Gauss-Jordan y tambien usando la factorizacion PA = LU bull
36 SISTEMAS TRIDIAGONALES
Un caso muy importante de sistemas de ecuaciones lineales que requiere un tratamiento especial es el de los sistemas tridiagonales Tales sistemas aparecen en diversas aplicaciones como por ejemplo al utilizar metodos de diferencias flnitas en la solucion de problemas con valores en la frontera para ecuaciones diferenciales ordinarias y como veremos mas adelante en el problema de la interpolacion segmentaria cubica
Un sistema tridiagonal es de la forma
dx1 + C1X2
a 2 x + d2 x2 + C 2 X 3
a3 x2 + d3 x3 + C 3 X 4
an-1Xn- 2 + dn_ xn_ + cn_xn = bn_1
anxn_1 + dnxn = bn
La matriz de coeficientes del sistema es
C 0 0 od1
0a 2 d2 c2
o a3 d3 c 3
A =
o
la cual se dice una matriz TRIDIAGONAL
En general
ij=12 n
suponiendo que la matriz A de podemos usar eliminaclon Gausslana eliminacion Gaussiana simple es adecuadl EDD por filas) Otra forma de resolver
EDD por filas es a partir de la componentes de las matrices L Y U as triangular inferior con todos sus superior Para encontrar tales malnce5
o 0 0
Y2 o o o o
L =
o Yn t
o o
Como Cl Ct
Y2 Cl Y2C +az
0 h a2
LU =
0
0
d
=
Capitulo 3 SOLUCION NUMERIC A DE SISTEMAS DE ECUACIONES 113
Instrucci6n en DERIVE
RESUELVA_1(Ab) Simplifica en la soluci6n exacta X del sistema AX=b EI vector b se entra como un vector fila 0
En el ejemplo anterior la eliminaci6n Gaussiana condujo a una respuesta defectuosa de un sistema de ecuaciones lineales bien condicionado Esto muestra la inestabilidad numerica del algoritmo de eliminaci6n Gaussiana (consecuencia de la divisi6n por un nOmero (pivote) pequero) Hay sin embargo situaciones en las cuales el algoritmo de eliminaci6n Gaussiana es numericamente estable EI siguiente teorema cuya demostraci6n puede consultarse en Burden 1985 paginas 366 y 367 se refiere a una de tales situaciones
Teorema 35 Si A = (ai j )n xn es una matriz estrictamente dominante diagonalmente
(EDD) por filas es decir si
n Iaii Igt II ai j I para cada i == 12 n
)= 1 ) 1
entonces A es invertible (no-singular) Ademas se puede rea lizar eliminaci6n Gaussiana sin intercambio de filas en cualquier sistema AX = b para obtener su Onica soluci6n y los calculos son estables con respecto al crecimiento de los errores de redondeo V
N6tese que como consecuencia del teorema anterior se tiene que Si A = (a l ) )n xn es EDD
por filas entonces A tiene factorizaci6n LU es decir A == LU con L triangular inferior con unos en su diagonal principal y U triangular superior (escalonada)
sea tambien
Observe que la matriz de coeficientes del ejemplo 33 anterior no es EDD por filas
sectJeor9~~am9i~ -valioo=paramatrj simetricas
(vease Bllrden 19 5 ~ina 368) Una matriz A E Rnxn ~ i~etrica se dice definida
positiva si satisface una cualquiera de las sig uientes condiciones (Ias- cualesson equival~tes) a _-----
- - -- shy
i) ~AX gt 0 paralodo 2CERn X 0 J
iiXTodos los valores propios de A son pos1tlvOS7
iii) Todos los pivotes obtenidos en la eliminaci6n Gaussiana sobre A sin intercambio de filas - ~on positiv~s
iv) Todas las submatrices principales de A tienen determinante positiv~
(Las submatrices principales de la matriz A == (ai J) son las matrices - n ~ n
1k a ] k == 12bull f) )
akk v
114 METODOS NUMERICOS
N6tese nuevamente que Si A E R nxn es simetrica y definida positiva entonces A tiene
factorizaci6n A = LU con L triangular inferior con sus componentes sabre la diagonal principal iguales a uno y U triangular superior (ver mas adelante factorizaci6n de Choleski)
Observe que la matriz de coeficientes de l ejemplo 33 anterior no es simetrica
34 ESTRATEGIAS DE PIVOTEO
EI eJemplo 33 anterior muestra una de las dificultades que pueden surgir al aplicar el metodo
de eliminaci6n Gaussiana cuando el pivote ajj-1) es pequeno comparado con algunos
(j-1) e ementos I a t para J ~ I t ~ n
Para tratar de evitar tales dificultades se introduce en el metodo de eliminaci6n Gaussiana una estrategia IIamada de pivoteo la cual consiste en seleccionar el pivote de acuerdo con un cierto criter io Nosotros usaremos dos estrategias la estrategia de pivoteo maximo por columna 0 pivoteo parcial y la estrategia de pivoteo escalado de fila 0 escalamiento
341 Pivoteo maximo por columna 0 pivoteo parcial Esta estrategia difiere de
eliminaci6n Gaussiana simple unicamente en la escogencia del pivote ajj-1) la cual se hace
ahora as
Para j = 12 n - 1 se determina el menor entero k J ~ k ~ n tal que
y
es decir seleccionamos el primer elemento diferente de cero sobre la columna j-esima a
partir de la j-esima fila y que tenga mayor valor absoluto (para j = 1 a~j j1 ) = al~) ak 1 )
Si tal k no existe el sistema no tiene soluci6n unica y el proceso se puede terminar
Si tal k existe y k ot j entonces hacemos intercambio de las ecuaciones j-esima y k-esima
y continuamos con la eliminaci6n Gaussiana l
Ilustremos esta estrategia para resolver el sistema
E1 03X1 + 589x2 = 592 E2 531x1 - 61 OX 2 = 470
entonces k = 2 ot 1= j asl que
Por sustituci6n regresiva
Observe que en este caSO
SWAP(A i j)
que es el mismo del ejemplo 33 usando aritmetica con redondeo a tres digitos
Como para j = 1 se tiene que
Capitulo 3 SOLUCION NUMERICA DE SISTEMAS DE ECUACIONES 115
Max Iall I Ia21 I = Max 03 531 = 531 = Ia21 I 0
entonces k = 2 1= j as i que intercambiamos El y E2 Y continuamos con la el imi nac i6n
531 - 610 (A b)= ( 03 589 592) 470) (531 - 610 470 03 589 592
E2I (~) ( m21=5 65 ~ 1O - l) 531 - 610 5 31 470) ------------------~) ( o 589 589
Por sustituci6n regresiva obtenemos
x = 589 = 100 Xl = 470 +610(100) = 531 = 100 2 589 531 531
Observe que en este caso X= (J es la soluci6n exacta del sistema dado bull
Instrucci6n en DERIVE
SWAP(A i j) Intercambia las filas (0 elementos) i y j de la matriz A (de un vector) 0
Nota En el procedimiento de pivoteo maximo por columna (p ivoteo parcial) cada multiplicador mi j es tal que
y aunque esta estrategia permite resolver satisfactoriamente muchos sistemas de ecuaciones lineales hay casos donde fracasa como se ilustra en el siguiente ejemplo
Ejemplo 34 Consideremos el sistema
El 300Xl + 58900x2 59200 E2 531x1 -610X2 - 470
el cual es un sistema equivalente al del ejemplo 33 (los coeficientes de la primera ecuaci6n
en el sistema del ejemplo 33 han sido multiplicados por 10 3) EI pivoteo maximo por
columna con aritmetica de redondeo a tres dig itos nos lIeva a los siguientes resultados
59200) E2 1 (-~pound) (m 2 1 =177) ) ( 300 58900 59200 ) (A b) = ( 300 58900 531 - 610 470 o - 10400 -10500
y por sustituci6n regresiva x2 = 101 Y Xl = - 100 que es la misma solucion que se obtiene si ~ usamos eliminaci6n Gaussiana simple
116 METODOS NUMERICOS
En casas como el de este ejemplo donde un pivote es mucho mas pequeno que alguno de los coeficientes de la ecuacion que el encabeza se recomienda la tecnica conocida como pivoteo escalado de fila 0 escalamiento la cual es nuestra segunda estrategia
342 Pivoteo escalado de fila Esta tecnica solo difiere de la eliminaci6n Gaussiana simple al igual que el pivoteo parcial en la escogencia del pivote
lEsta vez el pivote alr ) se escoge como se indica a continuacion
Para j = 12 n - 1 hacemos 10 siguiente
a) Para i = j j + 1 n calculamos
S = Max Ia(l- l) I Factor de escala I j s i s n II
Si Sj = 0 entonces el sistema no tiene solucion unica y el proceso se puede terminar
b) Para i = jj+1 n calculamos
c) Encontramos el menor entero k con j o k 0 n tal que
Ia~jjl ) I IaH- l) I -------- = Max --------
Sk j S j $ n SI
Si tal k no existe entonces el sistema no tiene solucion unica y el proceso se puede terminar
Si tal k existe y kF j entonces hacemos intercambio de las ecuaciones j-esima y kshyesima
y continuamos con la eliminacion Gaussiana V
Apliquemos esta estrategia para resolver el sistema del ejemplo 34 usando aritmetica con redondeo a tres dfgitos
Para j = 1
a) Sl = Max Iall I I a12 I = Max 300 58900 = 58900F 0 Y
Max 53~ 610 IshyS2 =Max a21 a n -shy
b) Ahora ~= 51
la0
l 300
S2
531 al l I a21 -Max - 0 -
c) Max -S -S- - 5890 NIl 1 2
las ecuaciolllltanto intercamblamos
Gaussiana
300 58900 (A b) = ( 531 -610
Por sustitucion regresiva
Xl tOO
Observe que en este GaSO X
35 FACTORIZACION
Capitulo 3 SOLUCION NUMERICA DE SISTEMAS DE ECUACIONES 117
8 2 =Max I a21 I I a22 1 =Max 531 610 = 610 0
b) Ahora
Ia I 300 - 8- = 58900 Y
1
Ia 2 I 531 - - = - shy
6108 2
c) Max ~ ~ I=Max ~ 531 = 531 =~ 0 as que k = 2 1= j y por8 8 2 58900 60 60 8 2
tanto intercambiamos las ecuaciones E y E2 Y continuamos con la eliminacion
Gaussiana
(A b) = (300 58900 59200) p ( 531 -610 470 ) 2 )
531 - 610 470 300 58900 59200
E2 (300) (m2 =565 ) ( 531 -6 1 0531 47 0 )
--------------~) o 58900 58900
Por sustituci6n regresiva
_ 470 +610(100) 531X2 = 100 Xl = =-- = 100
531 531
Observe que en este caso X= ( ~ J = ( 100) es la solucion exacta del sistema bull 100x2
35 FACTORIZACION TRIANGULAR
Consideremos un sistema AX = b con A no-singular y b 0 Con respecto a la matriz A se sabe que existen matrices P de permutaci6n L triangular inferior con sus componentes sobre la diagonal principal iguales a uno y U triangular superior (escalon ada) tales que PA = LU Una forma eficiente computacionalmente de encontrar P L Y U usando eliminaci6n Gaussiana se muestra en el siguiente ejemplo Ejemplo 35 Resuelva el siguiente sistema de ecuaciones lineales usando eliminacion Gaussiana con pivoteo parcial y obtenga una factorizaci6n PA =LU para la matriz A de coeficientes asociada con este metodo
- x + 2X2 - 3X3 = -2 3x1 - 3x2 - X 3 = -4
x + x2 = 3IEmpezamos introduciendo un vector p =(PP2P3) T el cual inicializamos_con P = i i = 123 Y
donde se almacenaran los intercambios necesarios en el proceso de eliminacion Gaussiana
118 METODOS NUMERICOS
con pivoteo parcial (el numero de componentes del vector p coincide con el orden n del sistema a resolver)
- 1 2 -3 -2] -4] - 2 P = (123) T (Ab) = ~ - 3 -1 - 4
[ o 3 3
3 - 3 - 1 - 4(m=-H (m3l = ~ )
E( - ~) E3(~) 10 10(-i) 3 3
2(i) 1 13 3 3
(Observe que cada multiplicador mij es almacenado en la posici6n correspondiente (i j) en
la matriz de trabajo)
3 - 3 -1 - 4
Max 121 =2 p=( 23 1) T bull (i) 2 1
-
3 13 3
(-i) 10 3
10 3
(Observe que la permutaci6n se hace para las filas 2 y 3 completas es decir incluyendo los multiplicadores)
3 - 3 - 1 --4(mJ2 = ~ ) E3( ~)
2(i) 13 3 3
(- i) (i) 7 11 2 2
La eficiencia en el metoda indicado se debe a que en la misma matriz de trabajo se almacenan los multiplicadores que van a conformar la matriz L (en el ejemplo son los numeros que se encuentran dentro de parentesis) 10 que significa un ahorro de memoria y como los intercambios necesarios afectan simultaneamente a las matrices L y U se evita tener que volver a la matriz original a realizar los intercambios ya observados y repetir la eliminaci6n Gaussiana con pivoteo parcial De esta manera al terminar el proceso de eliminaci6n podemos leer en la matriz final la parte estrictamente triangular inferior de L (son los numeros entre parentesis) y la matriz triangular superior U (que es la parte triangular superior de la matriz final) y en el vector p final quedan almacenados los intercambios realizados que se usan para producir la matriz de permutaci6n P
Para el ejemplo 35
3 -31 0 0 1 0 2U=0L= 3 1 1 0 0-3 2
(Verifique que PA =LU)
Para obtener la soluci6n del sistema
reducido
yobtenemos
C6mo se resuelve el
Capitulo 3 SOLUCION NUMERICA DE SISTEMAS DE ECUACIONES 119
Para el ejemplo 35
1 0 0 3 -3 - 1 1 1 1
L= - 0 U= 0 2 - Y como P~ (231) entonces [~ 0 lp3 3 1 1 7 0
-- - 0 0 3 2 2
(Verifique que PA = LU )
Para obtener la soluci6n del sistema original usamos sustituci6n regresiva en el sistema reducido
-4 13
UX = 3
11
2 ~
c y obtenemos
asi que la soluci6n (exacta) del sistema dado es X = (23 40 ~)T21 21 7
C6mo se resuelve el sistema AX = b a partir de la factorizaci6n PA = LU obtenida bull
351 Algunas aplicaciones de la factorizaci6n PA = LU La factorizaci6n PA = LU es utilizada eficientemente en aquellos casos donde se trabaja repetfdamente con la misma matriz A Dos de esos casos se presentan a continuaci6n
1) Resolver varios sistemas AX =b con la misma matriz de coeficientes A ya que en p L Y U esta almacenado todo el proceso de eliminaci6n Gaussiana EI algoritmo se bas a en la siguiente equivalencia
AX = b lt=gt PAX = Pb lt=gt LUX = Pb
lt=gt UX = L- 1Pb
lt=gtUXoo c y c oo L- 1Pb
lt=gt UX == C Y Lc = Pb lt=gt Lc = Pb Y UX = c
Los pasos a seguir son
Paso 1 Calcular Pb Paso 2 Resolver para c Lc = Pb por sustituci6n progresiva Paso 3 Resolver para X UX = c por sustituci6n regresiva
Como ejercicio resuelva el sistema del ejemplo anterior usando este algoritmo
120 METODOS NUMERICOS
2) Encontrar la matriz inversa A -1 de una matriz invertible A resolviendo los n-sistemas
AX =eU) j = 12 n
donde e(1) =(O O1 o 0f ERn
t posicion j
La solucion X del sistema AX = e(n j = 12 n produce la correspondiente columna j shy
esima de la matriz A -1
Como ejercicio compare el numero de operaciones para encontrar A - usando el metodo de Gauss-Jordan con el numero de operaciones resolviendo los n-sistemas indicados antes V
Ejercicio 33 Calcule la inversa de la matriz A de coeficientes del sistema del ejemplo 35 usando el metodo de Gauss-Jordan y tambien usando la factorizacion PA = LU bull
36 SISTEMAS TRIDIAGONALES
Un caso muy importante de sistemas de ecuaciones lineales que requiere un tratamiento especial es el de los sistemas tridiagonales Tales sistemas aparecen en diversas aplicaciones como por ejemplo al utilizar metodos de diferencias flnitas en la solucion de problemas con valores en la frontera para ecuaciones diferenciales ordinarias y como veremos mas adelante en el problema de la interpolacion segmentaria cubica
Un sistema tridiagonal es de la forma
dx1 + C1X2
a 2 x + d2 x2 + C 2 X 3
a3 x2 + d3 x3 + C 3 X 4
an-1Xn- 2 + dn_ xn_ + cn_xn = bn_1
anxn_1 + dnxn = bn
La matriz de coeficientes del sistema es
C 0 0 od1
0a 2 d2 c2
o a3 d3 c 3
A =
o
la cual se dice una matriz TRIDIAGONAL
En general
ij=12 n
suponiendo que la matriz A de podemos usar eliminaclon Gausslana eliminacion Gaussiana simple es adecuadl EDD por filas) Otra forma de resolver
EDD por filas es a partir de la componentes de las matrices L Y U as triangular inferior con todos sus superior Para encontrar tales malnce5
o 0 0
Y2 o o o o
L =
o Yn t
o o
Como Cl Ct
Y2 Cl Y2C +az
0 h a2
LU =
0
0
d
=
114 METODOS NUMERICOS
N6tese nuevamente que Si A E R nxn es simetrica y definida positiva entonces A tiene
factorizaci6n A = LU con L triangular inferior con sus componentes sabre la diagonal principal iguales a uno y U triangular superior (ver mas adelante factorizaci6n de Choleski)
Observe que la matriz de coeficientes de l ejemplo 33 anterior no es simetrica
34 ESTRATEGIAS DE PIVOTEO
EI eJemplo 33 anterior muestra una de las dificultades que pueden surgir al aplicar el metodo
de eliminaci6n Gaussiana cuando el pivote ajj-1) es pequeno comparado con algunos
(j-1) e ementos I a t para J ~ I t ~ n
Para tratar de evitar tales dificultades se introduce en el metodo de eliminaci6n Gaussiana una estrategia IIamada de pivoteo la cual consiste en seleccionar el pivote de acuerdo con un cierto criter io Nosotros usaremos dos estrategias la estrategia de pivoteo maximo por columna 0 pivoteo parcial y la estrategia de pivoteo escalado de fila 0 escalamiento
341 Pivoteo maximo por columna 0 pivoteo parcial Esta estrategia difiere de
eliminaci6n Gaussiana simple unicamente en la escogencia del pivote ajj-1) la cual se hace
ahora as
Para j = 12 n - 1 se determina el menor entero k J ~ k ~ n tal que
y
es decir seleccionamos el primer elemento diferente de cero sobre la columna j-esima a
partir de la j-esima fila y que tenga mayor valor absoluto (para j = 1 a~j j1 ) = al~) ak 1 )
Si tal k no existe el sistema no tiene soluci6n unica y el proceso se puede terminar
Si tal k existe y k ot j entonces hacemos intercambio de las ecuaciones j-esima y k-esima
y continuamos con la eliminaci6n Gaussiana l
Ilustremos esta estrategia para resolver el sistema
E1 03X1 + 589x2 = 592 E2 531x1 - 61 OX 2 = 470
entonces k = 2 ot 1= j asl que
Por sustituci6n regresiva
Observe que en este caSO
SWAP(A i j)
que es el mismo del ejemplo 33 usando aritmetica con redondeo a tres digitos
Como para j = 1 se tiene que
Capitulo 3 SOLUCION NUMERICA DE SISTEMAS DE ECUACIONES 115
Max Iall I Ia21 I = Max 03 531 = 531 = Ia21 I 0
entonces k = 2 1= j as i que intercambiamos El y E2 Y continuamos con la el imi nac i6n
531 - 610 (A b)= ( 03 589 592) 470) (531 - 610 470 03 589 592
E2I (~) ( m21=5 65 ~ 1O - l) 531 - 610 5 31 470) ------------------~) ( o 589 589
Por sustituci6n regresiva obtenemos
x = 589 = 100 Xl = 470 +610(100) = 531 = 100 2 589 531 531
Observe que en este caso X= (J es la soluci6n exacta del sistema dado bull
Instrucci6n en DERIVE
SWAP(A i j) Intercambia las filas (0 elementos) i y j de la matriz A (de un vector) 0
Nota En el procedimiento de pivoteo maximo por columna (p ivoteo parcial) cada multiplicador mi j es tal que
y aunque esta estrategia permite resolver satisfactoriamente muchos sistemas de ecuaciones lineales hay casos donde fracasa como se ilustra en el siguiente ejemplo
Ejemplo 34 Consideremos el sistema
El 300Xl + 58900x2 59200 E2 531x1 -610X2 - 470
el cual es un sistema equivalente al del ejemplo 33 (los coeficientes de la primera ecuaci6n
en el sistema del ejemplo 33 han sido multiplicados por 10 3) EI pivoteo maximo por
columna con aritmetica de redondeo a tres dig itos nos lIeva a los siguientes resultados
59200) E2 1 (-~pound) (m 2 1 =177) ) ( 300 58900 59200 ) (A b) = ( 300 58900 531 - 610 470 o - 10400 -10500
y por sustituci6n regresiva x2 = 101 Y Xl = - 100 que es la misma solucion que se obtiene si ~ usamos eliminaci6n Gaussiana simple
116 METODOS NUMERICOS
En casas como el de este ejemplo donde un pivote es mucho mas pequeno que alguno de los coeficientes de la ecuacion que el encabeza se recomienda la tecnica conocida como pivoteo escalado de fila 0 escalamiento la cual es nuestra segunda estrategia
342 Pivoteo escalado de fila Esta tecnica solo difiere de la eliminaci6n Gaussiana simple al igual que el pivoteo parcial en la escogencia del pivote
lEsta vez el pivote alr ) se escoge como se indica a continuacion
Para j = 12 n - 1 hacemos 10 siguiente
a) Para i = j j + 1 n calculamos
S = Max Ia(l- l) I Factor de escala I j s i s n II
Si Sj = 0 entonces el sistema no tiene solucion unica y el proceso se puede terminar
b) Para i = jj+1 n calculamos
c) Encontramos el menor entero k con j o k 0 n tal que
Ia~jjl ) I IaH- l) I -------- = Max --------
Sk j S j $ n SI
Si tal k no existe entonces el sistema no tiene solucion unica y el proceso se puede terminar
Si tal k existe y kF j entonces hacemos intercambio de las ecuaciones j-esima y kshyesima
y continuamos con la eliminacion Gaussiana V
Apliquemos esta estrategia para resolver el sistema del ejemplo 34 usando aritmetica con redondeo a tres dfgitos
Para j = 1
a) Sl = Max Iall I I a12 I = Max 300 58900 = 58900F 0 Y
Max 53~ 610 IshyS2 =Max a21 a n -shy
b) Ahora ~= 51
la0
l 300
S2
531 al l I a21 -Max - 0 -
c) Max -S -S- - 5890 NIl 1 2
las ecuaciolllltanto intercamblamos
Gaussiana
300 58900 (A b) = ( 531 -610
Por sustitucion regresiva
Xl tOO
Observe que en este GaSO X
35 FACTORIZACION
Capitulo 3 SOLUCION NUMERICA DE SISTEMAS DE ECUACIONES 117
8 2 =Max I a21 I I a22 1 =Max 531 610 = 610 0
b) Ahora
Ia I 300 - 8- = 58900 Y
1
Ia 2 I 531 - - = - shy
6108 2
c) Max ~ ~ I=Max ~ 531 = 531 =~ 0 as que k = 2 1= j y por8 8 2 58900 60 60 8 2
tanto intercambiamos las ecuaciones E y E2 Y continuamos con la eliminacion
Gaussiana
(A b) = (300 58900 59200) p ( 531 -610 470 ) 2 )
531 - 610 470 300 58900 59200
E2 (300) (m2 =565 ) ( 531 -6 1 0531 47 0 )
--------------~) o 58900 58900
Por sustituci6n regresiva
_ 470 +610(100) 531X2 = 100 Xl = =-- = 100
531 531
Observe que en este caso X= ( ~ J = ( 100) es la solucion exacta del sistema bull 100x2
35 FACTORIZACION TRIANGULAR
Consideremos un sistema AX = b con A no-singular y b 0 Con respecto a la matriz A se sabe que existen matrices P de permutaci6n L triangular inferior con sus componentes sobre la diagonal principal iguales a uno y U triangular superior (escalon ada) tales que PA = LU Una forma eficiente computacionalmente de encontrar P L Y U usando eliminaci6n Gaussiana se muestra en el siguiente ejemplo Ejemplo 35 Resuelva el siguiente sistema de ecuaciones lineales usando eliminacion Gaussiana con pivoteo parcial y obtenga una factorizaci6n PA =LU para la matriz A de coeficientes asociada con este metodo
- x + 2X2 - 3X3 = -2 3x1 - 3x2 - X 3 = -4
x + x2 = 3IEmpezamos introduciendo un vector p =(PP2P3) T el cual inicializamos_con P = i i = 123 Y
donde se almacenaran los intercambios necesarios en el proceso de eliminacion Gaussiana
118 METODOS NUMERICOS
con pivoteo parcial (el numero de componentes del vector p coincide con el orden n del sistema a resolver)
- 1 2 -3 -2] -4] - 2 P = (123) T (Ab) = ~ - 3 -1 - 4
[ o 3 3
3 - 3 - 1 - 4(m=-H (m3l = ~ )
E( - ~) E3(~) 10 10(-i) 3 3
2(i) 1 13 3 3
(Observe que cada multiplicador mij es almacenado en la posici6n correspondiente (i j) en
la matriz de trabajo)
3 - 3 -1 - 4
Max 121 =2 p=( 23 1) T bull (i) 2 1
-
3 13 3
(-i) 10 3
10 3
(Observe que la permutaci6n se hace para las filas 2 y 3 completas es decir incluyendo los multiplicadores)
3 - 3 - 1 --4(mJ2 = ~ ) E3( ~)
2(i) 13 3 3
(- i) (i) 7 11 2 2
La eficiencia en el metoda indicado se debe a que en la misma matriz de trabajo se almacenan los multiplicadores que van a conformar la matriz L (en el ejemplo son los numeros que se encuentran dentro de parentesis) 10 que significa un ahorro de memoria y como los intercambios necesarios afectan simultaneamente a las matrices L y U se evita tener que volver a la matriz original a realizar los intercambios ya observados y repetir la eliminaci6n Gaussiana con pivoteo parcial De esta manera al terminar el proceso de eliminaci6n podemos leer en la matriz final la parte estrictamente triangular inferior de L (son los numeros entre parentesis) y la matriz triangular superior U (que es la parte triangular superior de la matriz final) y en el vector p final quedan almacenados los intercambios realizados que se usan para producir la matriz de permutaci6n P
Para el ejemplo 35
3 -31 0 0 1 0 2U=0L= 3 1 1 0 0-3 2
(Verifique que PA =LU)
Para obtener la soluci6n del sistema
reducido
yobtenemos
C6mo se resuelve el
Capitulo 3 SOLUCION NUMERICA DE SISTEMAS DE ECUACIONES 119
Para el ejemplo 35
1 0 0 3 -3 - 1 1 1 1
L= - 0 U= 0 2 - Y como P~ (231) entonces [~ 0 lp3 3 1 1 7 0
-- - 0 0 3 2 2
(Verifique que PA = LU )
Para obtener la soluci6n del sistema original usamos sustituci6n regresiva en el sistema reducido
-4 13
UX = 3
11
2 ~
c y obtenemos
asi que la soluci6n (exacta) del sistema dado es X = (23 40 ~)T21 21 7
C6mo se resuelve el sistema AX = b a partir de la factorizaci6n PA = LU obtenida bull
351 Algunas aplicaciones de la factorizaci6n PA = LU La factorizaci6n PA = LU es utilizada eficientemente en aquellos casos donde se trabaja repetfdamente con la misma matriz A Dos de esos casos se presentan a continuaci6n
1) Resolver varios sistemas AX =b con la misma matriz de coeficientes A ya que en p L Y U esta almacenado todo el proceso de eliminaci6n Gaussiana EI algoritmo se bas a en la siguiente equivalencia
AX = b lt=gt PAX = Pb lt=gt LUX = Pb
lt=gt UX = L- 1Pb
lt=gtUXoo c y c oo L- 1Pb
lt=gt UX == C Y Lc = Pb lt=gt Lc = Pb Y UX = c
Los pasos a seguir son
Paso 1 Calcular Pb Paso 2 Resolver para c Lc = Pb por sustituci6n progresiva Paso 3 Resolver para X UX = c por sustituci6n regresiva
Como ejercicio resuelva el sistema del ejemplo anterior usando este algoritmo
120 METODOS NUMERICOS
2) Encontrar la matriz inversa A -1 de una matriz invertible A resolviendo los n-sistemas
AX =eU) j = 12 n
donde e(1) =(O O1 o 0f ERn
t posicion j
La solucion X del sistema AX = e(n j = 12 n produce la correspondiente columna j shy
esima de la matriz A -1
Como ejercicio compare el numero de operaciones para encontrar A - usando el metodo de Gauss-Jordan con el numero de operaciones resolviendo los n-sistemas indicados antes V
Ejercicio 33 Calcule la inversa de la matriz A de coeficientes del sistema del ejemplo 35 usando el metodo de Gauss-Jordan y tambien usando la factorizacion PA = LU bull
36 SISTEMAS TRIDIAGONALES
Un caso muy importante de sistemas de ecuaciones lineales que requiere un tratamiento especial es el de los sistemas tridiagonales Tales sistemas aparecen en diversas aplicaciones como por ejemplo al utilizar metodos de diferencias flnitas en la solucion de problemas con valores en la frontera para ecuaciones diferenciales ordinarias y como veremos mas adelante en el problema de la interpolacion segmentaria cubica
Un sistema tridiagonal es de la forma
dx1 + C1X2
a 2 x + d2 x2 + C 2 X 3
a3 x2 + d3 x3 + C 3 X 4
an-1Xn- 2 + dn_ xn_ + cn_xn = bn_1
anxn_1 + dnxn = bn
La matriz de coeficientes del sistema es
C 0 0 od1
0a 2 d2 c2
o a3 d3 c 3
A =
o
la cual se dice una matriz TRIDIAGONAL
En general
ij=12 n
suponiendo que la matriz A de podemos usar eliminaclon Gausslana eliminacion Gaussiana simple es adecuadl EDD por filas) Otra forma de resolver
EDD por filas es a partir de la componentes de las matrices L Y U as triangular inferior con todos sus superior Para encontrar tales malnce5
o 0 0
Y2 o o o o
L =
o Yn t
o o
Como Cl Ct
Y2 Cl Y2C +az
0 h a2
LU =
0
0
d
=
Capitulo 3 SOLUCION NUMERICA DE SISTEMAS DE ECUACIONES 115
Max Iall I Ia21 I = Max 03 531 = 531 = Ia21 I 0
entonces k = 2 1= j as i que intercambiamos El y E2 Y continuamos con la el imi nac i6n
531 - 610 (A b)= ( 03 589 592) 470) (531 - 610 470 03 589 592
E2I (~) ( m21=5 65 ~ 1O - l) 531 - 610 5 31 470) ------------------~) ( o 589 589
Por sustituci6n regresiva obtenemos
x = 589 = 100 Xl = 470 +610(100) = 531 = 100 2 589 531 531
Observe que en este caso X= (J es la soluci6n exacta del sistema dado bull
Instrucci6n en DERIVE
SWAP(A i j) Intercambia las filas (0 elementos) i y j de la matriz A (de un vector) 0
Nota En el procedimiento de pivoteo maximo por columna (p ivoteo parcial) cada multiplicador mi j es tal que
y aunque esta estrategia permite resolver satisfactoriamente muchos sistemas de ecuaciones lineales hay casos donde fracasa como se ilustra en el siguiente ejemplo
Ejemplo 34 Consideremos el sistema
El 300Xl + 58900x2 59200 E2 531x1 -610X2 - 470
el cual es un sistema equivalente al del ejemplo 33 (los coeficientes de la primera ecuaci6n
en el sistema del ejemplo 33 han sido multiplicados por 10 3) EI pivoteo maximo por
columna con aritmetica de redondeo a tres dig itos nos lIeva a los siguientes resultados
59200) E2 1 (-~pound) (m 2 1 =177) ) ( 300 58900 59200 ) (A b) = ( 300 58900 531 - 610 470 o - 10400 -10500
y por sustituci6n regresiva x2 = 101 Y Xl = - 100 que es la misma solucion que se obtiene si ~ usamos eliminaci6n Gaussiana simple
116 METODOS NUMERICOS
En casas como el de este ejemplo donde un pivote es mucho mas pequeno que alguno de los coeficientes de la ecuacion que el encabeza se recomienda la tecnica conocida como pivoteo escalado de fila 0 escalamiento la cual es nuestra segunda estrategia
342 Pivoteo escalado de fila Esta tecnica solo difiere de la eliminaci6n Gaussiana simple al igual que el pivoteo parcial en la escogencia del pivote
lEsta vez el pivote alr ) se escoge como se indica a continuacion
Para j = 12 n - 1 hacemos 10 siguiente
a) Para i = j j + 1 n calculamos
S = Max Ia(l- l) I Factor de escala I j s i s n II
Si Sj = 0 entonces el sistema no tiene solucion unica y el proceso se puede terminar
b) Para i = jj+1 n calculamos
c) Encontramos el menor entero k con j o k 0 n tal que
Ia~jjl ) I IaH- l) I -------- = Max --------
Sk j S j $ n SI
Si tal k no existe entonces el sistema no tiene solucion unica y el proceso se puede terminar
Si tal k existe y kF j entonces hacemos intercambio de las ecuaciones j-esima y kshyesima
y continuamos con la eliminacion Gaussiana V
Apliquemos esta estrategia para resolver el sistema del ejemplo 34 usando aritmetica con redondeo a tres dfgitos
Para j = 1
a) Sl = Max Iall I I a12 I = Max 300 58900 = 58900F 0 Y
Max 53~ 610 IshyS2 =Max a21 a n -shy
b) Ahora ~= 51
la0
l 300
S2
531 al l I a21 -Max - 0 -
c) Max -S -S- - 5890 NIl 1 2
las ecuaciolllltanto intercamblamos
Gaussiana
300 58900 (A b) = ( 531 -610
Por sustitucion regresiva
Xl tOO
Observe que en este GaSO X
35 FACTORIZACION
Capitulo 3 SOLUCION NUMERICA DE SISTEMAS DE ECUACIONES 117
8 2 =Max I a21 I I a22 1 =Max 531 610 = 610 0
b) Ahora
Ia I 300 - 8- = 58900 Y
1
Ia 2 I 531 - - = - shy
6108 2
c) Max ~ ~ I=Max ~ 531 = 531 =~ 0 as que k = 2 1= j y por8 8 2 58900 60 60 8 2
tanto intercambiamos las ecuaciones E y E2 Y continuamos con la eliminacion
Gaussiana
(A b) = (300 58900 59200) p ( 531 -610 470 ) 2 )
531 - 610 470 300 58900 59200
E2 (300) (m2 =565 ) ( 531 -6 1 0531 47 0 )
--------------~) o 58900 58900
Por sustituci6n regresiva
_ 470 +610(100) 531X2 = 100 Xl = =-- = 100
531 531
Observe que en este caso X= ( ~ J = ( 100) es la solucion exacta del sistema bull 100x2
35 FACTORIZACION TRIANGULAR
Consideremos un sistema AX = b con A no-singular y b 0 Con respecto a la matriz A se sabe que existen matrices P de permutaci6n L triangular inferior con sus componentes sobre la diagonal principal iguales a uno y U triangular superior (escalon ada) tales que PA = LU Una forma eficiente computacionalmente de encontrar P L Y U usando eliminaci6n Gaussiana se muestra en el siguiente ejemplo Ejemplo 35 Resuelva el siguiente sistema de ecuaciones lineales usando eliminacion Gaussiana con pivoteo parcial y obtenga una factorizaci6n PA =LU para la matriz A de coeficientes asociada con este metodo
- x + 2X2 - 3X3 = -2 3x1 - 3x2 - X 3 = -4
x + x2 = 3IEmpezamos introduciendo un vector p =(PP2P3) T el cual inicializamos_con P = i i = 123 Y
donde se almacenaran los intercambios necesarios en el proceso de eliminacion Gaussiana
118 METODOS NUMERICOS
con pivoteo parcial (el numero de componentes del vector p coincide con el orden n del sistema a resolver)
- 1 2 -3 -2] -4] - 2 P = (123) T (Ab) = ~ - 3 -1 - 4
[ o 3 3
3 - 3 - 1 - 4(m=-H (m3l = ~ )
E( - ~) E3(~) 10 10(-i) 3 3
2(i) 1 13 3 3
(Observe que cada multiplicador mij es almacenado en la posici6n correspondiente (i j) en
la matriz de trabajo)
3 - 3 -1 - 4
Max 121 =2 p=( 23 1) T bull (i) 2 1
-
3 13 3
(-i) 10 3
10 3
(Observe que la permutaci6n se hace para las filas 2 y 3 completas es decir incluyendo los multiplicadores)
3 - 3 - 1 --4(mJ2 = ~ ) E3( ~)
2(i) 13 3 3
(- i) (i) 7 11 2 2
La eficiencia en el metoda indicado se debe a que en la misma matriz de trabajo se almacenan los multiplicadores que van a conformar la matriz L (en el ejemplo son los numeros que se encuentran dentro de parentesis) 10 que significa un ahorro de memoria y como los intercambios necesarios afectan simultaneamente a las matrices L y U se evita tener que volver a la matriz original a realizar los intercambios ya observados y repetir la eliminaci6n Gaussiana con pivoteo parcial De esta manera al terminar el proceso de eliminaci6n podemos leer en la matriz final la parte estrictamente triangular inferior de L (son los numeros entre parentesis) y la matriz triangular superior U (que es la parte triangular superior de la matriz final) y en el vector p final quedan almacenados los intercambios realizados que se usan para producir la matriz de permutaci6n P
Para el ejemplo 35
3 -31 0 0 1 0 2U=0L= 3 1 1 0 0-3 2
(Verifique que PA =LU)
Para obtener la soluci6n del sistema
reducido
yobtenemos
C6mo se resuelve el
Capitulo 3 SOLUCION NUMERICA DE SISTEMAS DE ECUACIONES 119
Para el ejemplo 35
1 0 0 3 -3 - 1 1 1 1
L= - 0 U= 0 2 - Y como P~ (231) entonces [~ 0 lp3 3 1 1 7 0
-- - 0 0 3 2 2
(Verifique que PA = LU )
Para obtener la soluci6n del sistema original usamos sustituci6n regresiva en el sistema reducido
-4 13
UX = 3
11
2 ~
c y obtenemos
asi que la soluci6n (exacta) del sistema dado es X = (23 40 ~)T21 21 7
C6mo se resuelve el sistema AX = b a partir de la factorizaci6n PA = LU obtenida bull
351 Algunas aplicaciones de la factorizaci6n PA = LU La factorizaci6n PA = LU es utilizada eficientemente en aquellos casos donde se trabaja repetfdamente con la misma matriz A Dos de esos casos se presentan a continuaci6n
1) Resolver varios sistemas AX =b con la misma matriz de coeficientes A ya que en p L Y U esta almacenado todo el proceso de eliminaci6n Gaussiana EI algoritmo se bas a en la siguiente equivalencia
AX = b lt=gt PAX = Pb lt=gt LUX = Pb
lt=gt UX = L- 1Pb
lt=gtUXoo c y c oo L- 1Pb
lt=gt UX == C Y Lc = Pb lt=gt Lc = Pb Y UX = c
Los pasos a seguir son
Paso 1 Calcular Pb Paso 2 Resolver para c Lc = Pb por sustituci6n progresiva Paso 3 Resolver para X UX = c por sustituci6n regresiva
Como ejercicio resuelva el sistema del ejemplo anterior usando este algoritmo
120 METODOS NUMERICOS
2) Encontrar la matriz inversa A -1 de una matriz invertible A resolviendo los n-sistemas
AX =eU) j = 12 n
donde e(1) =(O O1 o 0f ERn
t posicion j
La solucion X del sistema AX = e(n j = 12 n produce la correspondiente columna j shy
esima de la matriz A -1
Como ejercicio compare el numero de operaciones para encontrar A - usando el metodo de Gauss-Jordan con el numero de operaciones resolviendo los n-sistemas indicados antes V
Ejercicio 33 Calcule la inversa de la matriz A de coeficientes del sistema del ejemplo 35 usando el metodo de Gauss-Jordan y tambien usando la factorizacion PA = LU bull
36 SISTEMAS TRIDIAGONALES
Un caso muy importante de sistemas de ecuaciones lineales que requiere un tratamiento especial es el de los sistemas tridiagonales Tales sistemas aparecen en diversas aplicaciones como por ejemplo al utilizar metodos de diferencias flnitas en la solucion de problemas con valores en la frontera para ecuaciones diferenciales ordinarias y como veremos mas adelante en el problema de la interpolacion segmentaria cubica
Un sistema tridiagonal es de la forma
dx1 + C1X2
a 2 x + d2 x2 + C 2 X 3
a3 x2 + d3 x3 + C 3 X 4
an-1Xn- 2 + dn_ xn_ + cn_xn = bn_1
anxn_1 + dnxn = bn
La matriz de coeficientes del sistema es
C 0 0 od1
0a 2 d2 c2
o a3 d3 c 3
A =
o
la cual se dice una matriz TRIDIAGONAL
En general
ij=12 n
suponiendo que la matriz A de podemos usar eliminaclon Gausslana eliminacion Gaussiana simple es adecuadl EDD por filas) Otra forma de resolver
EDD por filas es a partir de la componentes de las matrices L Y U as triangular inferior con todos sus superior Para encontrar tales malnce5
o 0 0
Y2 o o o o
L =
o Yn t
o o
Como Cl Ct
Y2 Cl Y2C +az
0 h a2
LU =
0
0
d
=
116 METODOS NUMERICOS
En casas como el de este ejemplo donde un pivote es mucho mas pequeno que alguno de los coeficientes de la ecuacion que el encabeza se recomienda la tecnica conocida como pivoteo escalado de fila 0 escalamiento la cual es nuestra segunda estrategia
342 Pivoteo escalado de fila Esta tecnica solo difiere de la eliminaci6n Gaussiana simple al igual que el pivoteo parcial en la escogencia del pivote
lEsta vez el pivote alr ) se escoge como se indica a continuacion
Para j = 12 n - 1 hacemos 10 siguiente
a) Para i = j j + 1 n calculamos
S = Max Ia(l- l) I Factor de escala I j s i s n II
Si Sj = 0 entonces el sistema no tiene solucion unica y el proceso se puede terminar
b) Para i = jj+1 n calculamos
c) Encontramos el menor entero k con j o k 0 n tal que
Ia~jjl ) I IaH- l) I -------- = Max --------
Sk j S j $ n SI
Si tal k no existe entonces el sistema no tiene solucion unica y el proceso se puede terminar
Si tal k existe y kF j entonces hacemos intercambio de las ecuaciones j-esima y kshyesima
y continuamos con la eliminacion Gaussiana V
Apliquemos esta estrategia para resolver el sistema del ejemplo 34 usando aritmetica con redondeo a tres dfgitos
Para j = 1
a) Sl = Max Iall I I a12 I = Max 300 58900 = 58900F 0 Y
Max 53~ 610 IshyS2 =Max a21 a n -shy
b) Ahora ~= 51
la0
l 300
S2
531 al l I a21 -Max - 0 -
c) Max -S -S- - 5890 NIl 1 2
las ecuaciolllltanto intercamblamos
Gaussiana
300 58900 (A b) = ( 531 -610
Por sustitucion regresiva
Xl tOO
Observe que en este GaSO X
35 FACTORIZACION
Capitulo 3 SOLUCION NUMERICA DE SISTEMAS DE ECUACIONES 117
8 2 =Max I a21 I I a22 1 =Max 531 610 = 610 0
b) Ahora
Ia I 300 - 8- = 58900 Y
1
Ia 2 I 531 - - = - shy
6108 2
c) Max ~ ~ I=Max ~ 531 = 531 =~ 0 as que k = 2 1= j y por8 8 2 58900 60 60 8 2
tanto intercambiamos las ecuaciones E y E2 Y continuamos con la eliminacion
Gaussiana
(A b) = (300 58900 59200) p ( 531 -610 470 ) 2 )
531 - 610 470 300 58900 59200
E2 (300) (m2 =565 ) ( 531 -6 1 0531 47 0 )
--------------~) o 58900 58900
Por sustituci6n regresiva
_ 470 +610(100) 531X2 = 100 Xl = =-- = 100
531 531
Observe que en este caso X= ( ~ J = ( 100) es la solucion exacta del sistema bull 100x2
35 FACTORIZACION TRIANGULAR
Consideremos un sistema AX = b con A no-singular y b 0 Con respecto a la matriz A se sabe que existen matrices P de permutaci6n L triangular inferior con sus componentes sobre la diagonal principal iguales a uno y U triangular superior (escalon ada) tales que PA = LU Una forma eficiente computacionalmente de encontrar P L Y U usando eliminaci6n Gaussiana se muestra en el siguiente ejemplo Ejemplo 35 Resuelva el siguiente sistema de ecuaciones lineales usando eliminacion Gaussiana con pivoteo parcial y obtenga una factorizaci6n PA =LU para la matriz A de coeficientes asociada con este metodo
- x + 2X2 - 3X3 = -2 3x1 - 3x2 - X 3 = -4
x + x2 = 3IEmpezamos introduciendo un vector p =(PP2P3) T el cual inicializamos_con P = i i = 123 Y
donde se almacenaran los intercambios necesarios en el proceso de eliminacion Gaussiana
118 METODOS NUMERICOS
con pivoteo parcial (el numero de componentes del vector p coincide con el orden n del sistema a resolver)
- 1 2 -3 -2] -4] - 2 P = (123) T (Ab) = ~ - 3 -1 - 4
[ o 3 3
3 - 3 - 1 - 4(m=-H (m3l = ~ )
E( - ~) E3(~) 10 10(-i) 3 3
2(i) 1 13 3 3
(Observe que cada multiplicador mij es almacenado en la posici6n correspondiente (i j) en
la matriz de trabajo)
3 - 3 -1 - 4
Max 121 =2 p=( 23 1) T bull (i) 2 1
-
3 13 3
(-i) 10 3
10 3
(Observe que la permutaci6n se hace para las filas 2 y 3 completas es decir incluyendo los multiplicadores)
3 - 3 - 1 --4(mJ2 = ~ ) E3( ~)
2(i) 13 3 3
(- i) (i) 7 11 2 2
La eficiencia en el metoda indicado se debe a que en la misma matriz de trabajo se almacenan los multiplicadores que van a conformar la matriz L (en el ejemplo son los numeros que se encuentran dentro de parentesis) 10 que significa un ahorro de memoria y como los intercambios necesarios afectan simultaneamente a las matrices L y U se evita tener que volver a la matriz original a realizar los intercambios ya observados y repetir la eliminaci6n Gaussiana con pivoteo parcial De esta manera al terminar el proceso de eliminaci6n podemos leer en la matriz final la parte estrictamente triangular inferior de L (son los numeros entre parentesis) y la matriz triangular superior U (que es la parte triangular superior de la matriz final) y en el vector p final quedan almacenados los intercambios realizados que se usan para producir la matriz de permutaci6n P
Para el ejemplo 35
3 -31 0 0 1 0 2U=0L= 3 1 1 0 0-3 2
(Verifique que PA =LU)
Para obtener la soluci6n del sistema
reducido
yobtenemos
C6mo se resuelve el
Capitulo 3 SOLUCION NUMERICA DE SISTEMAS DE ECUACIONES 119
Para el ejemplo 35
1 0 0 3 -3 - 1 1 1 1
L= - 0 U= 0 2 - Y como P~ (231) entonces [~ 0 lp3 3 1 1 7 0
-- - 0 0 3 2 2
(Verifique que PA = LU )
Para obtener la soluci6n del sistema original usamos sustituci6n regresiva en el sistema reducido
-4 13
UX = 3
11
2 ~
c y obtenemos
asi que la soluci6n (exacta) del sistema dado es X = (23 40 ~)T21 21 7
C6mo se resuelve el sistema AX = b a partir de la factorizaci6n PA = LU obtenida bull
351 Algunas aplicaciones de la factorizaci6n PA = LU La factorizaci6n PA = LU es utilizada eficientemente en aquellos casos donde se trabaja repetfdamente con la misma matriz A Dos de esos casos se presentan a continuaci6n
1) Resolver varios sistemas AX =b con la misma matriz de coeficientes A ya que en p L Y U esta almacenado todo el proceso de eliminaci6n Gaussiana EI algoritmo se bas a en la siguiente equivalencia
AX = b lt=gt PAX = Pb lt=gt LUX = Pb
lt=gt UX = L- 1Pb
lt=gtUXoo c y c oo L- 1Pb
lt=gt UX == C Y Lc = Pb lt=gt Lc = Pb Y UX = c
Los pasos a seguir son
Paso 1 Calcular Pb Paso 2 Resolver para c Lc = Pb por sustituci6n progresiva Paso 3 Resolver para X UX = c por sustituci6n regresiva
Como ejercicio resuelva el sistema del ejemplo anterior usando este algoritmo
120 METODOS NUMERICOS
2) Encontrar la matriz inversa A -1 de una matriz invertible A resolviendo los n-sistemas
AX =eU) j = 12 n
donde e(1) =(O O1 o 0f ERn
t posicion j
La solucion X del sistema AX = e(n j = 12 n produce la correspondiente columna j shy
esima de la matriz A -1
Como ejercicio compare el numero de operaciones para encontrar A - usando el metodo de Gauss-Jordan con el numero de operaciones resolviendo los n-sistemas indicados antes V
Ejercicio 33 Calcule la inversa de la matriz A de coeficientes del sistema del ejemplo 35 usando el metodo de Gauss-Jordan y tambien usando la factorizacion PA = LU bull
36 SISTEMAS TRIDIAGONALES
Un caso muy importante de sistemas de ecuaciones lineales que requiere un tratamiento especial es el de los sistemas tridiagonales Tales sistemas aparecen en diversas aplicaciones como por ejemplo al utilizar metodos de diferencias flnitas en la solucion de problemas con valores en la frontera para ecuaciones diferenciales ordinarias y como veremos mas adelante en el problema de la interpolacion segmentaria cubica
Un sistema tridiagonal es de la forma
dx1 + C1X2
a 2 x + d2 x2 + C 2 X 3
a3 x2 + d3 x3 + C 3 X 4
an-1Xn- 2 + dn_ xn_ + cn_xn = bn_1
anxn_1 + dnxn = bn
La matriz de coeficientes del sistema es
C 0 0 od1
0a 2 d2 c2
o a3 d3 c 3
A =
o
la cual se dice una matriz TRIDIAGONAL
En general
ij=12 n
suponiendo que la matriz A de podemos usar eliminaclon Gausslana eliminacion Gaussiana simple es adecuadl EDD por filas) Otra forma de resolver
EDD por filas es a partir de la componentes de las matrices L Y U as triangular inferior con todos sus superior Para encontrar tales malnce5
o 0 0
Y2 o o o o
L =
o Yn t
o o
Como Cl Ct
Y2 Cl Y2C +az
0 h a2
LU =
0
0
d
=
Capitulo 3 SOLUCION NUMERICA DE SISTEMAS DE ECUACIONES 117
8 2 =Max I a21 I I a22 1 =Max 531 610 = 610 0
b) Ahora
Ia I 300 - 8- = 58900 Y
1
Ia 2 I 531 - - = - shy
6108 2
c) Max ~ ~ I=Max ~ 531 = 531 =~ 0 as que k = 2 1= j y por8 8 2 58900 60 60 8 2
tanto intercambiamos las ecuaciones E y E2 Y continuamos con la eliminacion
Gaussiana
(A b) = (300 58900 59200) p ( 531 -610 470 ) 2 )
531 - 610 470 300 58900 59200
E2 (300) (m2 =565 ) ( 531 -6 1 0531 47 0 )
--------------~) o 58900 58900
Por sustituci6n regresiva
_ 470 +610(100) 531X2 = 100 Xl = =-- = 100
531 531
Observe que en este caso X= ( ~ J = ( 100) es la solucion exacta del sistema bull 100x2
35 FACTORIZACION TRIANGULAR
Consideremos un sistema AX = b con A no-singular y b 0 Con respecto a la matriz A se sabe que existen matrices P de permutaci6n L triangular inferior con sus componentes sobre la diagonal principal iguales a uno y U triangular superior (escalon ada) tales que PA = LU Una forma eficiente computacionalmente de encontrar P L Y U usando eliminaci6n Gaussiana se muestra en el siguiente ejemplo Ejemplo 35 Resuelva el siguiente sistema de ecuaciones lineales usando eliminacion Gaussiana con pivoteo parcial y obtenga una factorizaci6n PA =LU para la matriz A de coeficientes asociada con este metodo
- x + 2X2 - 3X3 = -2 3x1 - 3x2 - X 3 = -4
x + x2 = 3IEmpezamos introduciendo un vector p =(PP2P3) T el cual inicializamos_con P = i i = 123 Y
donde se almacenaran los intercambios necesarios en el proceso de eliminacion Gaussiana
118 METODOS NUMERICOS
con pivoteo parcial (el numero de componentes del vector p coincide con el orden n del sistema a resolver)
- 1 2 -3 -2] -4] - 2 P = (123) T (Ab) = ~ - 3 -1 - 4
[ o 3 3
3 - 3 - 1 - 4(m=-H (m3l = ~ )
E( - ~) E3(~) 10 10(-i) 3 3
2(i) 1 13 3 3
(Observe que cada multiplicador mij es almacenado en la posici6n correspondiente (i j) en
la matriz de trabajo)
3 - 3 -1 - 4
Max 121 =2 p=( 23 1) T bull (i) 2 1
-
3 13 3
(-i) 10 3
10 3
(Observe que la permutaci6n se hace para las filas 2 y 3 completas es decir incluyendo los multiplicadores)
3 - 3 - 1 --4(mJ2 = ~ ) E3( ~)
2(i) 13 3 3
(- i) (i) 7 11 2 2
La eficiencia en el metoda indicado se debe a que en la misma matriz de trabajo se almacenan los multiplicadores que van a conformar la matriz L (en el ejemplo son los numeros que se encuentran dentro de parentesis) 10 que significa un ahorro de memoria y como los intercambios necesarios afectan simultaneamente a las matrices L y U se evita tener que volver a la matriz original a realizar los intercambios ya observados y repetir la eliminaci6n Gaussiana con pivoteo parcial De esta manera al terminar el proceso de eliminaci6n podemos leer en la matriz final la parte estrictamente triangular inferior de L (son los numeros entre parentesis) y la matriz triangular superior U (que es la parte triangular superior de la matriz final) y en el vector p final quedan almacenados los intercambios realizados que se usan para producir la matriz de permutaci6n P
Para el ejemplo 35
3 -31 0 0 1 0 2U=0L= 3 1 1 0 0-3 2
(Verifique que PA =LU)
Para obtener la soluci6n del sistema
reducido
yobtenemos
C6mo se resuelve el
Capitulo 3 SOLUCION NUMERICA DE SISTEMAS DE ECUACIONES 119
Para el ejemplo 35
1 0 0 3 -3 - 1 1 1 1
L= - 0 U= 0 2 - Y como P~ (231) entonces [~ 0 lp3 3 1 1 7 0
-- - 0 0 3 2 2
(Verifique que PA = LU )
Para obtener la soluci6n del sistema original usamos sustituci6n regresiva en el sistema reducido
-4 13
UX = 3
11
2 ~
c y obtenemos
asi que la soluci6n (exacta) del sistema dado es X = (23 40 ~)T21 21 7
C6mo se resuelve el sistema AX = b a partir de la factorizaci6n PA = LU obtenida bull
351 Algunas aplicaciones de la factorizaci6n PA = LU La factorizaci6n PA = LU es utilizada eficientemente en aquellos casos donde se trabaja repetfdamente con la misma matriz A Dos de esos casos se presentan a continuaci6n
1) Resolver varios sistemas AX =b con la misma matriz de coeficientes A ya que en p L Y U esta almacenado todo el proceso de eliminaci6n Gaussiana EI algoritmo se bas a en la siguiente equivalencia
AX = b lt=gt PAX = Pb lt=gt LUX = Pb
lt=gt UX = L- 1Pb
lt=gtUXoo c y c oo L- 1Pb
lt=gt UX == C Y Lc = Pb lt=gt Lc = Pb Y UX = c
Los pasos a seguir son
Paso 1 Calcular Pb Paso 2 Resolver para c Lc = Pb por sustituci6n progresiva Paso 3 Resolver para X UX = c por sustituci6n regresiva
Como ejercicio resuelva el sistema del ejemplo anterior usando este algoritmo
120 METODOS NUMERICOS
2) Encontrar la matriz inversa A -1 de una matriz invertible A resolviendo los n-sistemas
AX =eU) j = 12 n
donde e(1) =(O O1 o 0f ERn
t posicion j
La solucion X del sistema AX = e(n j = 12 n produce la correspondiente columna j shy
esima de la matriz A -1
Como ejercicio compare el numero de operaciones para encontrar A - usando el metodo de Gauss-Jordan con el numero de operaciones resolviendo los n-sistemas indicados antes V
Ejercicio 33 Calcule la inversa de la matriz A de coeficientes del sistema del ejemplo 35 usando el metodo de Gauss-Jordan y tambien usando la factorizacion PA = LU bull
36 SISTEMAS TRIDIAGONALES
Un caso muy importante de sistemas de ecuaciones lineales que requiere un tratamiento especial es el de los sistemas tridiagonales Tales sistemas aparecen en diversas aplicaciones como por ejemplo al utilizar metodos de diferencias flnitas en la solucion de problemas con valores en la frontera para ecuaciones diferenciales ordinarias y como veremos mas adelante en el problema de la interpolacion segmentaria cubica
Un sistema tridiagonal es de la forma
dx1 + C1X2
a 2 x + d2 x2 + C 2 X 3
a3 x2 + d3 x3 + C 3 X 4
an-1Xn- 2 + dn_ xn_ + cn_xn = bn_1
anxn_1 + dnxn = bn
La matriz de coeficientes del sistema es
C 0 0 od1
0a 2 d2 c2
o a3 d3 c 3
A =
o
la cual se dice una matriz TRIDIAGONAL
En general
ij=12 n
suponiendo que la matriz A de podemos usar eliminaclon Gausslana eliminacion Gaussiana simple es adecuadl EDD por filas) Otra forma de resolver
EDD por filas es a partir de la componentes de las matrices L Y U as triangular inferior con todos sus superior Para encontrar tales malnce5
o 0 0
Y2 o o o o
L =
o Yn t
o o
Como Cl Ct
Y2 Cl Y2C +az
0 h a2
LU =
0
0
d
=
118 METODOS NUMERICOS
con pivoteo parcial (el numero de componentes del vector p coincide con el orden n del sistema a resolver)
- 1 2 -3 -2] -4] - 2 P = (123) T (Ab) = ~ - 3 -1 - 4
[ o 3 3
3 - 3 - 1 - 4(m=-H (m3l = ~ )
E( - ~) E3(~) 10 10(-i) 3 3
2(i) 1 13 3 3
(Observe que cada multiplicador mij es almacenado en la posici6n correspondiente (i j) en
la matriz de trabajo)
3 - 3 -1 - 4
Max 121 =2 p=( 23 1) T bull (i) 2 1
-
3 13 3
(-i) 10 3
10 3
(Observe que la permutaci6n se hace para las filas 2 y 3 completas es decir incluyendo los multiplicadores)
3 - 3 - 1 --4(mJ2 = ~ ) E3( ~)
2(i) 13 3 3
(- i) (i) 7 11 2 2
La eficiencia en el metoda indicado se debe a que en la misma matriz de trabajo se almacenan los multiplicadores que van a conformar la matriz L (en el ejemplo son los numeros que se encuentran dentro de parentesis) 10 que significa un ahorro de memoria y como los intercambios necesarios afectan simultaneamente a las matrices L y U se evita tener que volver a la matriz original a realizar los intercambios ya observados y repetir la eliminaci6n Gaussiana con pivoteo parcial De esta manera al terminar el proceso de eliminaci6n podemos leer en la matriz final la parte estrictamente triangular inferior de L (son los numeros entre parentesis) y la matriz triangular superior U (que es la parte triangular superior de la matriz final) y en el vector p final quedan almacenados los intercambios realizados que se usan para producir la matriz de permutaci6n P
Para el ejemplo 35
3 -31 0 0 1 0 2U=0L= 3 1 1 0 0-3 2
(Verifique que PA =LU)
Para obtener la soluci6n del sistema
reducido
yobtenemos
C6mo se resuelve el
Capitulo 3 SOLUCION NUMERICA DE SISTEMAS DE ECUACIONES 119
Para el ejemplo 35
1 0 0 3 -3 - 1 1 1 1
L= - 0 U= 0 2 - Y como P~ (231) entonces [~ 0 lp3 3 1 1 7 0
-- - 0 0 3 2 2
(Verifique que PA = LU )
Para obtener la soluci6n del sistema original usamos sustituci6n regresiva en el sistema reducido
-4 13
UX = 3
11
2 ~
c y obtenemos
asi que la soluci6n (exacta) del sistema dado es X = (23 40 ~)T21 21 7
C6mo se resuelve el sistema AX = b a partir de la factorizaci6n PA = LU obtenida bull
351 Algunas aplicaciones de la factorizaci6n PA = LU La factorizaci6n PA = LU es utilizada eficientemente en aquellos casos donde se trabaja repetfdamente con la misma matriz A Dos de esos casos se presentan a continuaci6n
1) Resolver varios sistemas AX =b con la misma matriz de coeficientes A ya que en p L Y U esta almacenado todo el proceso de eliminaci6n Gaussiana EI algoritmo se bas a en la siguiente equivalencia
AX = b lt=gt PAX = Pb lt=gt LUX = Pb
lt=gt UX = L- 1Pb
lt=gtUXoo c y c oo L- 1Pb
lt=gt UX == C Y Lc = Pb lt=gt Lc = Pb Y UX = c
Los pasos a seguir son
Paso 1 Calcular Pb Paso 2 Resolver para c Lc = Pb por sustituci6n progresiva Paso 3 Resolver para X UX = c por sustituci6n regresiva
Como ejercicio resuelva el sistema del ejemplo anterior usando este algoritmo
120 METODOS NUMERICOS
2) Encontrar la matriz inversa A -1 de una matriz invertible A resolviendo los n-sistemas
AX =eU) j = 12 n
donde e(1) =(O O1 o 0f ERn
t posicion j
La solucion X del sistema AX = e(n j = 12 n produce la correspondiente columna j shy
esima de la matriz A -1
Como ejercicio compare el numero de operaciones para encontrar A - usando el metodo de Gauss-Jordan con el numero de operaciones resolviendo los n-sistemas indicados antes V
Ejercicio 33 Calcule la inversa de la matriz A de coeficientes del sistema del ejemplo 35 usando el metodo de Gauss-Jordan y tambien usando la factorizacion PA = LU bull
36 SISTEMAS TRIDIAGONALES
Un caso muy importante de sistemas de ecuaciones lineales que requiere un tratamiento especial es el de los sistemas tridiagonales Tales sistemas aparecen en diversas aplicaciones como por ejemplo al utilizar metodos de diferencias flnitas en la solucion de problemas con valores en la frontera para ecuaciones diferenciales ordinarias y como veremos mas adelante en el problema de la interpolacion segmentaria cubica
Un sistema tridiagonal es de la forma
dx1 + C1X2
a 2 x + d2 x2 + C 2 X 3
a3 x2 + d3 x3 + C 3 X 4
an-1Xn- 2 + dn_ xn_ + cn_xn = bn_1
anxn_1 + dnxn = bn
La matriz de coeficientes del sistema es
C 0 0 od1
0a 2 d2 c2
o a3 d3 c 3
A =
o
la cual se dice una matriz TRIDIAGONAL
En general
ij=12 n
suponiendo que la matriz A de podemos usar eliminaclon Gausslana eliminacion Gaussiana simple es adecuadl EDD por filas) Otra forma de resolver
EDD por filas es a partir de la componentes de las matrices L Y U as triangular inferior con todos sus superior Para encontrar tales malnce5
o 0 0
Y2 o o o o
L =
o Yn t
o o
Como Cl Ct
Y2 Cl Y2C +az
0 h a2
LU =
0
0
d
=
Capitulo 3 SOLUCION NUMERICA DE SISTEMAS DE ECUACIONES 119
Para el ejemplo 35
1 0 0 3 -3 - 1 1 1 1
L= - 0 U= 0 2 - Y como P~ (231) entonces [~ 0 lp3 3 1 1 7 0
-- - 0 0 3 2 2
(Verifique que PA = LU )
Para obtener la soluci6n del sistema original usamos sustituci6n regresiva en el sistema reducido
-4 13
UX = 3
11
2 ~
c y obtenemos
asi que la soluci6n (exacta) del sistema dado es X = (23 40 ~)T21 21 7
C6mo se resuelve el sistema AX = b a partir de la factorizaci6n PA = LU obtenida bull
351 Algunas aplicaciones de la factorizaci6n PA = LU La factorizaci6n PA = LU es utilizada eficientemente en aquellos casos donde se trabaja repetfdamente con la misma matriz A Dos de esos casos se presentan a continuaci6n
1) Resolver varios sistemas AX =b con la misma matriz de coeficientes A ya que en p L Y U esta almacenado todo el proceso de eliminaci6n Gaussiana EI algoritmo se bas a en la siguiente equivalencia
AX = b lt=gt PAX = Pb lt=gt LUX = Pb
lt=gt UX = L- 1Pb
lt=gtUXoo c y c oo L- 1Pb
lt=gt UX == C Y Lc = Pb lt=gt Lc = Pb Y UX = c
Los pasos a seguir son
Paso 1 Calcular Pb Paso 2 Resolver para c Lc = Pb por sustituci6n progresiva Paso 3 Resolver para X UX = c por sustituci6n regresiva
Como ejercicio resuelva el sistema del ejemplo anterior usando este algoritmo
120 METODOS NUMERICOS
2) Encontrar la matriz inversa A -1 de una matriz invertible A resolviendo los n-sistemas
AX =eU) j = 12 n
donde e(1) =(O O1 o 0f ERn
t posicion j
La solucion X del sistema AX = e(n j = 12 n produce la correspondiente columna j shy
esima de la matriz A -1
Como ejercicio compare el numero de operaciones para encontrar A - usando el metodo de Gauss-Jordan con el numero de operaciones resolviendo los n-sistemas indicados antes V
Ejercicio 33 Calcule la inversa de la matriz A de coeficientes del sistema del ejemplo 35 usando el metodo de Gauss-Jordan y tambien usando la factorizacion PA = LU bull
36 SISTEMAS TRIDIAGONALES
Un caso muy importante de sistemas de ecuaciones lineales que requiere un tratamiento especial es el de los sistemas tridiagonales Tales sistemas aparecen en diversas aplicaciones como por ejemplo al utilizar metodos de diferencias flnitas en la solucion de problemas con valores en la frontera para ecuaciones diferenciales ordinarias y como veremos mas adelante en el problema de la interpolacion segmentaria cubica
Un sistema tridiagonal es de la forma
dx1 + C1X2
a 2 x + d2 x2 + C 2 X 3
a3 x2 + d3 x3 + C 3 X 4
an-1Xn- 2 + dn_ xn_ + cn_xn = bn_1
anxn_1 + dnxn = bn
La matriz de coeficientes del sistema es
C 0 0 od1
0a 2 d2 c2
o a3 d3 c 3
A =
o
la cual se dice una matriz TRIDIAGONAL
En general
ij=12 n
suponiendo que la matriz A de podemos usar eliminaclon Gausslana eliminacion Gaussiana simple es adecuadl EDD por filas) Otra forma de resolver
EDD por filas es a partir de la componentes de las matrices L Y U as triangular inferior con todos sus superior Para encontrar tales malnce5
o 0 0
Y2 o o o o
L =
o Yn t
o o
Como Cl Ct
Y2 Cl Y2C +az
0 h a2
LU =
0
0
d
=
120 METODOS NUMERICOS
2) Encontrar la matriz inversa A -1 de una matriz invertible A resolviendo los n-sistemas
AX =eU) j = 12 n
donde e(1) =(O O1 o 0f ERn
t posicion j
La solucion X del sistema AX = e(n j = 12 n produce la correspondiente columna j shy
esima de la matriz A -1
Como ejercicio compare el numero de operaciones para encontrar A - usando el metodo de Gauss-Jordan con el numero de operaciones resolviendo los n-sistemas indicados antes V
Ejercicio 33 Calcule la inversa de la matriz A de coeficientes del sistema del ejemplo 35 usando el metodo de Gauss-Jordan y tambien usando la factorizacion PA = LU bull
36 SISTEMAS TRIDIAGONALES
Un caso muy importante de sistemas de ecuaciones lineales que requiere un tratamiento especial es el de los sistemas tridiagonales Tales sistemas aparecen en diversas aplicaciones como por ejemplo al utilizar metodos de diferencias flnitas en la solucion de problemas con valores en la frontera para ecuaciones diferenciales ordinarias y como veremos mas adelante en el problema de la interpolacion segmentaria cubica
Un sistema tridiagonal es de la forma
dx1 + C1X2
a 2 x + d2 x2 + C 2 X 3
a3 x2 + d3 x3 + C 3 X 4
an-1Xn- 2 + dn_ xn_ + cn_xn = bn_1
anxn_1 + dnxn = bn
La matriz de coeficientes del sistema es
C 0 0 od1
0a 2 d2 c2
o a3 d3 c 3
A =
o
la cual se dice una matriz TRIDIAGONAL
En general
ij=12 n
suponiendo que la matriz A de podemos usar eliminaclon Gausslana eliminacion Gaussiana simple es adecuadl EDD por filas) Otra forma de resolver
EDD por filas es a partir de la componentes de las matrices L Y U as triangular inferior con todos sus superior Para encontrar tales malnce5
o 0 0
Y2 o o o o
L =
o Yn t
o o
Como Cl Ct
Y2 Cl Y2C +az
0 h a2
LU =
0
0
d
=