Transcript Sistema de Ecuaciones Lineales Metodos Iterativos

Metodos de Solución Iterativos




Empezar con una aproximación inicial para el
vector solución (x0)
Actualizar en cada iteración el vector x usando
el sistema Ax=b
Cada iteración involucra el producto matrizvector.
Si A es esparcida este producto es realizado
eficientemente.
1
Procedimiento de solución Iterativa




Escribir el sistema Ax=b en una forma
equivalente x=Tx+c
Empezando con x0, genere una secuencia de
aproximaciones {xk} iterativamente por
xk+1=Txk+c
Representación de T y c dependen del tipo de
método usado.
Pero para cada método T y c son obtenidas a
partir de A y b.
2
Convergencia




Cuando k, la secuencia {xk} converge a un
vector solución bajo algunas condiciones en la
Matriz T.
Esto impone condiciones diferentes en la matriz
A para diferentes métodos.
Para la misma matriz A, un método puede
converger mientras que otro puede divergir.
Por lo tanto para cada método la relación entre
A y T deben ser encontradas para decidir la
convergencia.
3
Diferentes metodos Iterativos



Iteración de Jacobi
Iteración de Gauss-Seidel
Successive Over Relaxation (S.O.R)


SOR es un método usado para acelerar la
convergencia.
La iteración de Gauss-Seidel es un caso especial del
método SOR.
4
Iteración de Jacobi
a11 x1  a12 x2    a1n xn  b1
a21 x1  a22 x2    a2 n xn  b2

an1 x1  an 2 x2    ann xn  bn
 x10 
 0
x2 
0

x 

 0
 xn 
1
1 
(b1  a12 x20    a1n xn0 )
k 1
xi  bi 
a11
aii 
1
x12 
(b2  a21 x10  a23 x30    a2 n xn0 )
a22
1
x1n 
(bn  an1 x10  an 2 x20    ann 1 xn01 )
ann
x11 

aij x   aij x 

j 1
j i 1

i 1
n
k
j
k
j
5
Método de Jacobi. Forma Matricial
 Descomponiendo A = D - L - U.
-U=triu(A)-D
-U
D
=
-L
-L=tril(A)-D
D=diag(diag(A))
6
xk+1=Txk+c - iteración por el método de
Jacobi
Se puede escribir como A=D-L-U (No es una
factorización)
 a11 a12
a
 21 a22
a31 a32
a13  a11 0
a23    0 a22
a33   0
0
0  0
0
0 0  a12
0    a21
0
0  0
0
a33   a31  a32 0 0
0
Ax=b  (D-L-U)x=b
1
k 1
xi 
aii
Dxk+1
i 1
n


k
k
bi   aij x j   aij x j 
j 1
j i 1


 
Lx k
Uxk
 a13 
 a23 
0 
Dxk+1 = (L+U)xk+b
xk+1=D-1(L+U)xk+D-1b
T=D-1(L+U)
c=D-1b
7
iteración Gauss-Seidel (GS)
Use lo último
al actualizar
a11 x1  a12 x2    a1n xn  b1
a21 x1  a22 x2    a2 n xn  b2

an1 x1  an 2 x2    ann xn  bn
1
1 
(b1  a12 x20    a1n xn0 )
k 1
xi  bi 
a11
aii 
1
x12 
(b2  a21 x11  a23 x30    a2 n xn0 )
a22
1
x1n 
(bn  an1 x11  an 2 x12    ann 1 x1n1 )
ann
x11 
 x10 
 0
x2 
0

x 

 0
 xn 
i 1
a x
j 1
ij
k 1
j

  aij x 
j i 1

n
k
j
8
x(k+1)=Tx(k)+x iteración de Gauss-Seidel
Ax=b  (D-L-U)x=b
k 1
i
x
1

aii
Dxk+1
i 1
n


k 1
k
bi   aij x j   aij x j 
j 1
j i 1


 
Lx k 1
(D-L)xk+1 =Uxk+b
Uxk
xk+1=(D-L)-1Uxk+(D-L)-1b
Tgs=(D-L)-1U
cgs=(D-L)-1b
9
Comparación



İteración de Gauss-Seidel converge más
rápidamente que la iteración de Jacobi desde
que este usa la última actualización.
Pero existen algunos casos que la iteración de
Jacobi converge pero Gauss-Seidel no.
El método de sobre relajación sucesiva es
usada para acelerar la convergencia del
método de Gauss-Seidel.
10
Metodo Sobre Relajación Sucesiva
(SOR)

Puede ser escrita como sigue
k 1
i
x
1
x 
aii
k
i
i 1
n


k 1
k
bi   aij x j   aij x j 
j 1
j i


xik 1  xik   ik
término Corrector
i2
xi3
xi2
xi1
0
i
x

 i1
 i0
i1
2
i
Multiplicando por
 1
Converge más
rápido
i0
11
SOR
xik 1  xik   ik
k 1
i
x
xik 1
xik 1
i 1
n


k 1
k
bi   aij x j   aij x j 
j 1
j i


i 1
n


1
k
k 1
k
 (1   ) xi  
bi   aij x j   aij x j 
aii 
j 1
j  i 1

 (1   ) xik   ~
xik 1
1
 x 
aii
k
i
Donde el ultimo termino es la estimación de Gauss-Seidel
1<<2 Sobre-relajación (convergencia rápida)
0<<1 Sub-relajación (convergencia más lenta)
Existe un valor óptimo para 
Encontrarlo por prueba y error
12
x(k+1)=Tx(k)+c iteración para SOR
k 1
i
x
1
 (1   ) x  
aii
k
i
i 1
n


k 1
k
bi   aij x j   aij x j 
j 1
j i 1


Dxk+1=(1-)Dxk+b+Lxk+1+Uxk
(D- L)xk+1=[(1-)D+U]xk+b
T=(D- L)-1[(1-)D+U]
c= (D- L)-1b
13
Convergencia de los métodos iterativos
Define el vector solución como
Define el vector error como
xˆ
ek
x  e  xˆ
k
k
Substituye esto en
x k 1  Tx k  c
ek 1  xˆ  T (ek  xˆ)  c  Txˆ  c  Tek
e
k 1
 Te  TTe
k
k 1
 TTTe
k 2
T
( k 1) 0
e
14
Convergencia de los Métodos Iterativos
iteración
e
k 1
T
( k 1) 0
e T
( k 1)
e
0
potencia
El método iterativo convergería para cualquier vector
inicial arbitrario si la siguiente condición es satisfecha
Condición de Convergencia
Lim e k 1  0 si Lim T ( k 1)  0
k 
k 
15
Norma de un vector
La norma de un vector debe satisfacer estas
condiciones:
x  0 Para cualquier vect orno nulo x
x  0 si y solo si x es un vector nulo
αx   x
Para un escalar α
x y  x  y
Las normas Vectoriales pueden ser definidas de
diferentes formas en tanto que la definición de
norma sea satisfecha.
16
Normas de vectores
Comunmente usadas
norma Suma o norma ℓ1
x 1  x1  x2   xn
norma Euclideana ó norma ℓ2
x 2  x  x  x
2
1
2
2
2
n
norma Máxima o norma ℓ
x   maxi xi
17
Norma de una matriz
La norma de una matriz debe
satisfacer estas condiciones:
A 0
A  0 si y solo si A es una matriznula
αA   A
para α escalar
A B  A  B
Importante identidad
Ax  A x
x es un vector
18
Normas de matrices mas usadas
Norma Máxima suma_columna o norma ℓ1
m
A 1  max  aij
1 j  n
i 1
Norma Espectral o norma ℓ2
A 2  maximo valor propio de A A
T
Norma Maxima suma_fila o norma ℓ
n
A   max aij
1i  m
j 1
19
Ejemplo

Calcule las normas ℓ1 y ℓ de la matriz
3 9 5
7 2 4 


6 8 1 
16
19
17  A 
13
15
10
 A1
20
Condición de Convergencia
lim e
k 
k 1
 0 si lim T
k 
( k 1)
0
Expresar T en terminos de matriz modal P y 
: Matriz Diagonal con valores propios de T en la diagonal
T  PP
1
T ( k 1)  PP 1 PP 1  PP 1
T ( k 1)  P( k 1) P 1
1k 1



k 1

2

k 1  




k 1 
n 

lim T ( k 1)  0  lim P( k 1) P 1  0  lim ( k 1)  0
k 
k 
k 
 lim ki 1  0   i  1 for i  1,2,...,n
k 
21
Condición Suficiente para convergencia
Si la magnitud de todos los valores propios de la
Matriz de iteración T es menor que 1 entonces la
iteración es convergente
Los valores propios son mas fácil de calcular que la
norma de una matriz
Tx  x
Tx   x 

   x  T x    T   (T )  T
Tx  T x 

 (T )  1
condición suficiente para convergencia
22
Convergencia de la iteración de Jacobi
T=D-1(L+U)

 0

 a21
 a22
T  



 an1
 a
 nn
a12

a11
0



a23

a22







ann 1

ann
a1n 

a11 

a2 n 

a22 
 
an 1n 


an 1n 1 

0 

23
Convergencia de la iteración de Jacobi
Evaluar la norma infinita (suma máxima fila) de T
T

n
aij
j 1
i j
aii
1 
n
 1 P ara i  1,2,...,n
 aii   aij
j 1
i j
Matriz Diagonal
estrictamente
Dominante
Si A es una matriz con diagonal estrictamente
dominante, entonces la iteración de Jacobi
converge para cualquier valor inicial
24
Criterios de Parada

Ax=b

En cualquier iteración k, el término residual es
rk=b-Axk

Verificar la norma del término residual
||b-Axk||

Si esto es menor que la cota del valor de
parada
25
Ejemplo 1 (Iteración de Jacobi)
 4  1 1  x1   7 
 4  8 1  x    21

 2  

 2 1 5  x3   15 
0 
x 0  0
0
b  Ax 0
2
 26.7395
Matriz Diagonal estrictamente dominante
7  x20  x30
7
x 
  1.75
4
4
0
0
21

4
x

x
21
1
3
x12 
  2.625
8
8
0
0
15

2
x

x
15
1
2
x31 
  3.0
5
5
1
1
b  Ax 1  10.0452
2
26
Ejemplo 1 continuación...
7  x12  x31
7  2.625 3
x 

 1.65625
4
4
21 4 x11  x31 21 4 1.75  3
2
x2 

 3.875
8
8
15  2 x11  x12 15  2 1.75  2.625
2
x3 

 4.225
5
5
2
1
7  3.875 4.225
 1.6625
4
21 4 1.65625 4.225
x23 
 3.98125
8
15  2 1.65625 3.875
x33 
 2.8875
5
b  Ax 2
2
 6.7413
x13 
b  Ax 3
2
 1.9534
Matriz es diagonal estrictamente dominante, las iteraciones de
Jacobi son convergentes.
27
Ejemplo 2
 2 1 5  x1   15 
 4  8 1  x    21

 2  

 4  1 1  x3   7 
0 
x 0  0
0
b  Ax 0
2
 26.7395
La matriz no es diagonal estrictamente dominante
 15  x20  5 x30  15
x 

 7.5
2
2
21 4 x10  x30
21
1
x2 
  2.625
8
8
x31  7  4 x10  x20
 7.0
1
1
b  Ax 1  54.8546
2
28
Ejemplo 2 continuación...
 15  2.625 5  7
 11.3125
2
21 4  7.5  7
x12 
 0.25
8
x31  7  4  7.5  2.625  39.625
x11 
b  Ax 2
2
 208 .3761
El término del residual aumenta en cada iteración,
de tal forma que las iteraciones divergen.
Note que la matriz no es diagonalmente
estrictamente dominante
Cuando la matriz no tiene diagonal estrictamente
dominante, puede converger como no.
29
Convergencia de la iteración de
Gauss-Seidel

Iteración GS converge para cualquier vector
inicial si A es una matriz diagonal
estrictamente dominante

Iteración GS converge para cualquier vector
inicial si A es una matriz simétrica y definida
positiva – La matriz A es definida positiva si
xTAx>0 para cualquier vector x no nulo.
30
Ejemplo1 (Iteración de Gauss-Seidel)
 4  1 1  x1   7 
 4  8 1  x    21

 2  

 2 1 5  x3   15 
0 
x 0  0
0
b  Ax 0
2
 26.7395
Matriz Diagonal estrictamente dominante
7
7  x20  x30

 1.75
x 
4
4
21 4 x11  x30 21 4 1.75
1

 3.5
x2 
8
8
15  2 x11  x12 15  2 1.75  3.5
1

 3.0
x3 
5
5
1
1
b  Ax 1
2
 3.0414
b  Ax 1  10.0452
2
İteración de
Jacobi
31
Ejemplo 1 continuación...
7  x12  x31
7  3.5  3
x 

 1.875
4
4
21 4 x12  x31 21 4 1.875 3
2
x2 

 3.9375
8
8
15  2 x12  x22 15  2 1.875 3.9375
2
x3 

 2.9625
5
5
2
1
b  Ax 2
b  Ax 2
2
2
 0.4765
 6.7413
Iteración de Jacobi
Cuando ambos métodos de Jacobi y GaussSeidel convergen, Gauss-Seidel converge más
rápido.
32
Convergencia del método SOR



Si 0<<2, método SOR converge para
cualquier valor inicial si A es una matriz
simétrica y definida positiva.
Si >2, método SOR diverge
Si 0<<1, SOR método converge pera la
velocidad de convergencia es mas lenta que el
método de Gauss-Seidel.
33
Conteo de operaciones




El # de operaciones para la Eliminación gaussiana o la
descomposición LU es de 0 (n3), orden de n3
Para los métodos iterativos, el número de
multiplicaciones escalares es 0 (n2) en cada iteración.
Si el número total de las iteraciones requeridas para la
convergencia es mucho menos que n, entonces los
métodos iterativos son más eficiente que métodos
directos.
Los Métodos iterativos también se satisfacen bien para
las matrices esparcidas.
34
Formas Matriciales. Resumen
La solución del sistema A x = b se obtiene
mediante la siguiente expresión recursiva:
x ( k ) = Tx ( k-1 ) + c
Método
Jacobi
Gauss-Seidel
SOR
A= D - L - U
T
c
D-1 (L+U)
D-1 b
( D -L)-1 U
( D -L)-1 b
(D-w L)-1 [(1-w) D + w U ]
w(D-w L)-1 b
35
Problema 1
Resolver el siguiente sistema por el método SOR,
considere ω=1.25.
4 x1  x2  2
 x1  4 x2  x3  6
 x2  4 x3  2
x  0 x2  0 x3  0
0
1
0
0
Aplicamos el metodo de SOR:
xik 1  (1   ) xik   ~
xik 1
36
Problema 1
2  x2
20
1
~
x1 

 0.5
4
4
1
0
x1   ~
x11  1    x1  1.25 x0.5  1  1.25x0  0.625
0
6  x1  x3
6  0.625  0
1
~
x2 

 1.65625
4
4
1
0
x2   ~
x21  1    x2  1.25 x1.65625  1  1.25x0 
1
0
 2.0703125
2  x2
2  2.0703125
~
x31 

 1.017578125
4
4
1
0
x3   ~
x31  1    x3 
1
 1.25 x1.017578125  1  1.25x0  1.24197265625
37
Problema 1
k
x1
x2
x3
0
0
0
0
1
0.625
2.0703125
1.2719727
2
1.1157227
2.1035767
0.9643745
3
1.003437
1.9640469
0.997671
4
0.9879054
2.0044809
1.0019825
5
1.0044239
2.0008818
0.9997799
38
Problema 2
Sea el sistema A x = b :
 Para k=-1, es la matriz A definida positiva?
 Para que valores de k el sistema converge, al usar el
método de Gauss-Seidel?
 Hacer 03 iteraciones de Gauss-Seidel para k=-3
 2 k   x1  6
  1 3   x   9 

 2   
39
Problema 2
A es definida positiva si:
xT Ax  0, para t odo vect orcolumna x no nulo
x1
 2  1  x1 
2
2
2
x2 

2
x

(
x

x
)

2
x
para t odo x no nulo.
1
1
2
2 0



  1 3   x2 
Observese que también satisface el criterio de
Silvester
40
Problema 2
TG  ( D  L) 1U
 2 0
( D  L)  

  1 3
1 / 2 0 
( D  L) 1  

1
/
6
1
/
3


0  0  k  0  k / 2 
 1/ 2
TG  
 0 0   0  k / 6 

1
/
6
1
/
3


 

k
Det TG  I  (   )( )
6
1  0 2   k / 6  (TG )  max   k / 6
Exist econvergencia cuando :  (TG )  1
Est o se cumplesiempreque : - 6  k  6
41
Problema 2
P ara Gauss - Seidel (k  -3)
x1( n 1)  3  1.5 x2( n )
x2( n 1)  3  x1( n 1) / 3
n
x1
x2
0
0
0
1
3
4
2
9
6
3
12
7
4
13.5
7.5
5
14.25
7.75
6
14.625
7.875
7
14.8125
7.9375
42