Sistema de Ecuaciones Lineales Metodos Iterativos
Download
Report
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 xn01 )
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 x1n1 )
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
1i 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 PP
1
T ( k 1) PP 1 PP 1 PP 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
20
1
~
x1
0.5
4
4
1
0
x1 ~
x11 1 x1 1.25 x0.5 1 1.25x0 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.25x0
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.25x0 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