Metodos Numericos-Sistemas de Ecuaciones Lineales, 1° ED. - Juan Piccini

Métodos Numéricos 2003 – IMERL - UdelaR Juan Piccini

SISTEMAS DE ECUACIONES LINEALES Notas basadas en el libro: “Numerical Analisis”, Ian Jacques, Colin Judd. METODOS ITERATIVOS

Veremos dos métodos iterativos para la solución del sistema bAx = donde

A es una matriz nxn no singular con elementos no nulos en la diagonal principal. Dichos métodos son el método de Jacobi y el método de Gauss-Seidel. La idea es partir de un punto (vector) inicial, e ir construyendo una sucesión de puntos (vectores) que converja a la solución del sistema. A efectos de facilitar el estudio , será conveniente introducir algo de notación. La

matriz A podemos escribirla como DULA ++= , donde DUL ,, son las

partes triangular inferior, triangular superior y diagonal de A respectivamente. Notemos que las matrices UL, no son las mismas de la escalerización Gaussiana. Por simplicidad, ilustraremos con un ejemplo 3x3.

Sea el sistema

3333232131

2323222121

1313212111

bxaxaxa

bxaxaxa

bxaxaxa

=++=++=++

. (1)

En este caso, tenemos que

=

=

=

33

22

11

23

1312

3231

21

00

00

00

000

00

0

0

00

000

a

a

a

Da

aa

U

aa

aL.

El primer paso es reacomodar la i-sima ecuación en (1) , de modo que

solamente el término que involucra a la variable ix quede a la izquierda del signo de igual.


Esto es,

2321313333

3231212222

3132121111

xaxabxa

xaxabxa

xaxabxa

−−=−−=

−−=

o, en notación matricial : xULbDx )( +−= (2) El siguiente paso es dividir ambos lados de la i-sima ecuación entre el elemento

diagonal iia , obteniendo

)(1

)(1

)(1

232131333

3

323121222

2

313212111

1

xaxaba

x

xaxaba

x

xaxaba

x

−−=

−−=

−−=

o, en notación matricial : ))((1 xULbDx +−= −.

Esta forma de reescribir (1) sugiere la iteración

[ ] [ ] [ ]( )[ ] [ ] [ ]( )[ ] [ ] [ ]( )mmm

mmm

mmm

xaxaba

x

xaxaba

x

xaxaba

x

232131333

13

323121222

12

313212111

11

1

1

1

−−=

−−=

−−=

+

+

+

(3)

o, en notación matricial :

[ ] [ ]( )mm xULbDx )(11 +−= −+ (4)

Este es el llamado Método de Jacobi. Para un sistema nxn, el método se define mediante


[ ] [ ]

−= ∑

=

≠=

+nj

ijj

mjiji

ii

mi xab

ax

1

1 1

Como es usual con los métodos iterativos, necesitaremos una aproximación

inicial [ ] [ ] [ ] [ ]( )00

20

10 ,...,, nxxxx = . Si no disponemos de información previa

alguna, lo usual es tomar [ ] ( )0,...,0,00 =x . Estos valores se sustituyen en los

lados derechos de (3), obteniendo [ ] [ ] [ ] [ ]( )11

21

11 ,...,, nxxxx = . Luego

sustituímos los nuevos valores en los lados derechos de (3) , obteniendo [ ] [ ] [ ] [ ]( )22

22

12 ,...,, nxxxx = y repetimos el proceso hasta que algún criterio de

parada sea satisfecho.

Por lo general, el criterio de parada suele ser del tipo [ ] [ ] ε≤−+ mm xx 1

,

donde el epsilon lo fija el usuario. Otro posible criterio puede ser [ ] ε≤− bAx m

.

Para ilustrar, veamos el siguiente ejemplo. Consideremos el sistema

742

436

625

321

321

321

=++=−+=−+

xxx

xxx

xxx

, cuya solución exacta es )1,1,1(=x .

La iteración de Jacobi viene dada por:

[ ] [ ] [ ]( )[ ] [ ] [ ]( )[ ] [ ] [ ]( )mmm

mmm

mmm

xxx

xxx

xxx

211

3

311

2

321

1

274

1

346

1

265

1

−−=

+−=

+−=

+

+

+


Como valor inicial tomamos (0,0,0), y paramos cuando

[ ] [ ]{ } 21 102

1max −+ ×≤− m

im

ii

xx

Los resultados se muestran en la siguiente tabla :

[ ] [ ] [ ]

002.1000.1001.111

999.0998.0997.010

995.0004.1000.19

010.1004.1008.18

004.1983.0989.07

962.0008.1988.06

034.1046.1053.15

084.1910.0977.04

773.0944.0860.03

983.0342.1283.12

750.1667.0200.11

0000321

mmm xxxm

Vemos que alcanzamos la solución correcta (hasta dos decimales) en once iteraciones. Consideremos ahora el mismo sistema, pero con la primera y segunda ecuaciones intercambiadas.

742

625

436

321

321

321

=++=−+=−+

xxx

xxx

xxx

La iteración de Jacobi queda ahora :

[ ] [ ] [ ]( )[ ] [ ] [ ]( )[ ] [ ] [ ]( )mmm

mmm

mmm

xxx

xxx

xxx

211

3

311

2

321

1

274

1

562

1

364

−−=

+−=

+−=

+

+

+

Partiendo de nuevo del (0,0,0), y con el mismo criterio de parada, obtenemos la siguiente tabla :


[ ] [ ] [ ]

278.83596.289625.4595

219.23547.87282.1194

656.7375.24750.373

000.1125.6750.82

750.1000.3000.41

0000321

−−−

−−−

mmm xxxm

Es claro que algo está pasando. Este ejemplo ilustra sobre la necesidad de tener alguna condición a priori para testear la convergencia (o no) del algoritmo. Resulta interesante notar que en (3), cuando usamos la segunda ecuación para

calcular [ ]12

+mx , el valor de [ ]11

+mx ya es conocido, de modo que podríamos

utilizarlo en el lugar de [ ]mx1 para calcular

[ ]12

+mx .

Como [ ]11

+mx está “más cerca” del verdadero valor 1x , es de esperar que esto

acelere la convergencia, produciendo un “mejor” valor [ ]12

+mx . Del mismo modo,

si en la tercera ecuación utilizo [ ]11

+mx y [ ]12

+mx en vez de [ ]mx1 y

[ ]mx2 para

calcular [ ]13

+mx , cabe esperar que la convergencia sea más rápida. Parecería entonces, que es mejor la iteración

[ ] [ ] [ ]( )[ ] [ ] [ ]( )[ ] [ ] [ ]( )1

2321

131333

13

3231

121122

12

313212111

11

1

1

1

+++

++

+

−−=

−−=

−−=

mmm

mmm

mmm

xaxaba

x

xaxaba

x

xaxaba

x

Esta se conoce como el método de Gauss-Seidel, y puesto que usa los valores más recientes apenas se obtienen, uno esperaría que converja más rápido que Jacobi.


Como el sistema de ecuaciones (1) puede reescribirse como

3333232131

3232222121

3132121111

bxaxaxa

xabxaxa

xaxabxa

=++−=+

−−=

en notación matricial, queda ( ) xUbxLD −=+ . La iteración de Gauss-Seidel puede escribirse entonces como,

( ) [ ] [ ]mm xUbxLD −=+ +1 o lo que es lo mismo,

[ ] ( ) [ ]( )mm UxbLDx −+= −+ 11 (5)

Para un sistema nxn, el método de Gauss-Seidel se define mediante

[ ] [ ] [ ] nixaxaba

xn

ij

mjij

i

j

mjiji

ii

mi ,...,2,1,

1

1

1

1

11 =

−−= ∑∑

+=

−

=

++

Apliquemos este método al ejemplo anterior. Veamos primero para el sistema

742

436

625

321

321

321

=++=−+=−+

xxx

xxx

xxx

Usaremos el mismo punto inicial y el mismo criterio de parada de antes. La iteración queda :

[ ] [ ] [ ]( )[ ] [ ] [ ]( )[ ] [ ] [ ]( )1

21

11

3

31

11

2

321

1

274

1

346

1

265

1

+++

++

+

−−=

+−=

+−=

mmm

mmm

mmm

xxx

xxx

xxx

Los resultados se ven en la siguiente tabla:


[ ] [ ] [ ]

001.1999.0999.07

998.0002.1003.16

004.1994.0995.05

987.0006.1024.14

019.1950.0987.03

895.0980.0220.12

033.1467.0200.11

0000321mmm xxxm

Notemos que el algoritmo converge en siete pasos (Jacobi necesitó once). Veamos ahora cuando en el sistema original intercambiamos las dos primeras ecuaciones.

742

625

436

321

321

321

=++=−+=−+

xxx

xxx

xxx

Ahora la iteración es :

[ ] [ ] [ ]( )[ ] [ ] [ ]( )[ ] [ ] [ ]( )1

21

11

3

31

11

2

321

1

274

1

562

1

364

+++

++

+

−−=

+−=

−−=

mmm

mmm

mmm

xxx

xxx

xxx

Los resultados se muestran en la siguiente tabla :

[ ] [ ] [ ]

156.95375.1894375.7603

125.750.122500.502

500.1000.7000.41

0000321

−−−

mmm xxxm


La divergencia ya se nota claramente en la tercera iteración. Antes de que pensemos que Gauss-Seidel es lo más grande que hay, mencionemos que existen ejemplos donde Jacobi converge más rápido que Gauss-Seidel, ejemplos donde uno converge y el otro no, y viceversa. Para comprender mejor qué podemos esperar, hagamos el siguiente análisis. Pongamos A M N= + , entonces la ecuación Ax b= puede escribirse como ( )M N x b+ = , o lo que es lo mismo ( )1x M b Nx−= − (6)

Esto sugiere la iteración [ ] [ ]( )1 1m mx M b Nx+ −= − (7) . Notemos que tanto Jacobi

como Gauss-Seidel son casos particulares de (7). Para Jacobi, ,M D N L U= = + y para Gauss-Seidel ,M D L N U= + = . Si a (7) le restamos (6), obtenemos

[ ] [ ]( )[ ] [ ]( )

1 1

1 1

m m

m m

x x M b Nx b Nx

x x M N x x

+ −

+ −

− = − − +

− = − −

Si llamamos [ ] [ ]1 1m me x x+ += − (el error en el paso m+1), llegamos a [ ] [ ]1 1m me M Ne+ −= − (8)

La matriz 1M N−− se llama matriz de iteración del sistema, y gobierna el cambio en el error de un paso al siguiente. Por comodidad llamemos Q a dicha matriz.

Para Jacobi y Gauss-Seidel tenemos ( ) ( ) 11 ,J G SQ D L U Q D L U−−

−= − + = − + .

Para la convergencia, necesitamos que [ ]1 0m

me +

→∞→ . Veremos una condición

necesaria y suficiente para que ello ocurra. Previamente definamos el radio espectral de una matriz, ( ) { }i

iA λρ max= es el máximo de los módulos de los

vectores propios de A. Hablando pronto y mal, el radio espectral mide el tamaño del mayor valor propio. Teorema (condición necesaria y suficiente). Supongamos que la matriz

1M N−− tiene valores propios nii ,...,1, =λ . Entonces la sucesión [ ] [ ]( )1 1m mx M b Nx+ −= − converge sii ( )1 1M Nρ −− < .

Dem. Asumiremos que existe una base de vectores propios { } ( ) nivvNMvv iiin ,...,1/,..., 1

1 =∀=− − λ

El vector de error inicial, [ ]0e puede escribirse como combinación lineal de dicha base, [ ]0

1 1 ... n ne v vα α= + + . Multiplicando ambos lados por 1Q M N−= − tenemos [ ]

nnnnn vvQvQve λαλααα ++=++= ..... 111111 Si multiplicamos de nuevo,

tendremos [ ]

nnnnnn vvQvvQe 21

211111

2 ..... λαλαλαλα ++=++=


Si repetimos esto, llegaremos a [ ]

nmnn

mm vve λαλα ++= ..111 (9).

Es claro que nisii imi ,...,110 =<→ λλ , de modo que [ ] ( )0 1me sii Qρ→ < .

Observaciones: Existen otras demostraciones para el caso en que la matriz de iteración no es diagonalizable. A partir de (9) vemos que la tasa de convergencia depende de m

iλ de modo que a menor radio espectral, más rápida será la convergencia. Como esta condición puede no ser fácil de chequear en la práctica, veremos otra condición (condición suficiente) para la convergencia, más manejable. Diremos que una matriz A es diagonal dominante cuando para cada fila de la matriz, el módulo del elemento diagonal es mayor que la suma de los módulos

del resto de la fila, o sea 1

1,...,n

ii ijjj i

a a i n=≠

> =∑

Teorema (Condición suficiente): Si la matriz de coeficientes del sistema es diagonal dominante, entonces el algoritmo converge. Dem. Si v es un vector propio de 1M N−− asociado al valor propio λ , entonces

vvNM λ=− −1 lo que puede rescribirse como ( ) 0M N vλ + = (10). Por el

teorema anterior, alcanza con probar que si A es diagonal dominante, entonces

( )1 1M Nρ −− < . Supongamos que la mayor componente del vector v es la i-sima,

esto es, j iv v j i≤ ≠ . La prueba para Jacobi y Gauss-Seidel se harán por

separado. Jacobi : Como ,M D N L U= = + , la ecuación (10) queda ( ) 0D L U vλ + + = .

La i-sima componente de dicha ecuación es 1

0n

ii i ij jjj i

a v a vλ=≠

+ =∑ , y entonces

1

1 n

ij jjii ij i

a va v

λ=≠

= − ∑ , lo cual implica 1 1

1 n ni

ij j ijj jii i ii ij i j i

va v a

a v a vλ

= =≠ ≠

= ≤∑ ∑

De modo que 1

11

n

ijjiij i

aa

λ=≠

= <∑ por ser A diagonal dominante.


Gauss-Seidel : Ahora ,M L D N U= + = , entonces la ecuación (10) queda

( ) 0D L U vλ λ+ + = . La i-sima componente de dicha ecuación es

1 1

0i n

ij j ij jj j i

a v a vλ= = +

+ =∑ ∑ , de donde 1

1

n

ij jj i

i

ij jj

a v

a vλ = +

=

=

∑

∑ (11). Para acotar esto,

cambiaremos el numerador por otro mayor, y el denominador por otro menor.

Numerador : 1 1

n n

ij j i ijj i j i

a v v a= + = +

≤∑ ∑

Denominador : 1 1

1 1 1

i i i

ij j ii i ij j ii i i ijj j j

a v a v a v a v v a− −

= = =≥ − ≥ −∑ ∑ ∑

Sustituyendo en (11), tenemos que

11

1

n

ij ij i

ii ijj

aa aλ

−= +

=

≤ −

∑∑ . Pero 1λ < es lo mismo que

1

1 1

n i

ij ii ijj i j

a a a−

= + =

< −∑ ∑ y

esto es lo mismo que ii

n

ijj

ij aa <∑≠=1

, que se cumple por hipótesis.

Metodos Numericos-Sistemas de Ecuaciones Lineales, 1° ED. - Juan Piccini

Documents

Transcript of Metodos Numericos-Sistemas de Ecuaciones Lineales, 1° ED. - Juan Piccini