1. No category

Métodos iterativos para la resolución de sistemas
Mayo de 2015.
Repasos de una situación conocida (CIN)
En CIN habéis estudiado como para resolver numéricamente una
ecuación, por ejemplo
x2 = 2,
podemos reescribir la ecuación de la forma x = f ( x ) y definir
una sucesión por la relación xn+1 = f ( xn ). Entonces si la sucesión
converge (ya la función f es continua) entonces el limite será una
solución de la ecuación. En el ejemplo anterior, en el método de
Newton, reescribimos la ecuación x2 = 2 como:
x = x − ( x2 − 2)/(2x ) = x/2 − 1/x
Luego para cada aproximación x de la solución calculamos una
nueva aproximación x 0 por medio de
x 0 = x/2 − 1/2
Elegimos algún valor inicial, por ejemplo x0 = 1. Obtenemos como
8 primeros valores
x=1.
L=[]
for i in [1..8]:
x = x/2 - 1/x
print x,
1.50000000000000 1.41666666666667 1.41421568627451 1.41421356237469
1.41421356237310 1.41421356237309 1.41421356237310 1.41421356237309
√
y observamos que la sucesión converge, efectivamente, hacía 2.
Sistemas lineales: método de Jacobi
La misma estrategia se aplica a sistema lineales compatibles
determinados (una única solución) AX = b.
Consideramos el sistema:
(
4x + 2y = 3
x + 3y = 4
En la primera ecuación despejamos la primera variable x, y en la
segunda ecuación la segunda variable y, pasando el resto al otro
lado (es el método de Jacobi):
(
4x = 3 − 2y
3y = 4 − x
métodos iterativos para la resolución de sistemas
2
Para cada aproximación ( x, y) de la solución definimos entonces
una nueva aproximación ( x 0 , y0 ) por la relación:
(
4x 0 = 3 − 2y
3y0 = 4 − x
Se elige alguna aproximación inicial, por ejemplo ( x0 , y0 ) = (0, 0).
La sucesión de vectores obtenida puede converger o no. Pero si
converge, entonces el limite es la solución del sistema.
En el ejemplo anterior:
x = 0.
y = 0.
for i in [1..10]:
xx = (3-2*y)/4
yy = (4-x)/3
(x, y) = (xx, yy)
print x, y
0.750000000000000 1.33333333333333
0.0833333333333334 1.08333333333333
0.208333333333333 1.30555555555556
0.0972222222222222 1.26388888888889
0.118055555555556 1.30092592592593
0.0995370370370371 1.29398148148148
0.103009259259259 1.30015432098765
0.0999228395061729 1.29899691358025
0.100501543209876 1.30002572016461
0.0999871399176955 1.29983281893004
observamos la convergencia de la sucesión.
Una mejora del sistema de Jacobi: el método de Gauss–Seidel
En el método de Jacobi para el sistema
(
4x + 2y = 3
x + 3y = 4
la nueva aproximación ( x 0 ; y0 ) se calcula a a partir de la aproximación anterior por medio de las formulas
(
4x 0 = 3 − 2y
3y0 = 4 − x
Pero en la formula para calcular y0 , podríamos utilizar x 0 que acaba
de ser calculada, en vez de la aproximación anterior x. Es el método
de Gauss–Seidel. Las relaciones entre cada aproximación ( x; y) y la
aproximación siguiente ( x 00 ; y00 ) queda así:
(
4x 00 = 3 − 2y
3y00 = 4 − x 00
Figura 1: Primeros pasos en el método
iterativo de Jacobi.
métodos iterativos para la resolución de sistemas
3
Calculemos los primeros términos de la sucesión de vectores obtenida de esta manera (con ( x0 ; y0 ) = (0; 0)) y observamos una
convergencia mucho más rápida:
x = 0.
y = 0.
for i in [1..10]:
x = (3-2*y)/4
y = (4-x)/3
print x, y
0.750000000000000
0.208333333333333
0.118055555555556
0.103009259259259
0.100501543209876
0.100083590534979
0.100013931755830
0.100002321959305
0.100000386993217
0.100000064498870
## comentario: aquí se utiliza el valor de x calculado justo a la línea anterior
1.08333333333333
1.26388888888889
1.29398148148148
1.29899691358025
1.29983281893004
1.29997213648834
1.29999535608139
1.29999922601357
1.29999987100226
1.29999997850038
Jacobi y Gauss–Seidel: otros ejemplos
Aquí está un ejemplo sin convergencia. Sea
"
#
" #
1 2
3
C=
,
b=
.
3 4
4
Entonces el sistema CX = b se escribe:
(
x + 2y = 3
3x + 4y = 4
Despejando x en la primera ecuación, e y en la segunda, reescribimos el sistema como:
(
x = 3 − 2y
y = (4 − 3x )/4
La sucesión que corresponde al método de Jacobi cumple:
(
xn+1 = 3 − 2yn
yn+1 = (4 − 3xn )/4
y la sucesión que corresponde al método de Gauss–Seidel cumple:
(
xn+1 = 3 − 2yn
yn+1 = (4 − 3xn+1 )/4
En este caso, ambos métodos divergen. Por ejemplo, para el método
de Gauss–Seidel, con ( x0 , y0 ) = (0, 0)
métodos iterativos para la resolución de sistemas
x=0.
y=0.
for i in [1..10]:
x = 3 -2*y
y = (4-3*x)/4
print x,y
0.650854139658888
0.217306746840088
0.108475689943148
0.0362177911400148
0.0180792816571913
0.00603629852333569
0.00301321360953181
0.00100604975388921
0.000502202268255239
0.000167674958981535
Aquí está un ejemplo con una matriz más grande. Sea


 
1
3 1 0 0
2
 1 3 1 0 

 

M=
, b =  
 0 1 3 1 
3
4
0 0 1 3
Consideramos aproximaciones ( x; y; z; t) para la solución de MX =
b. La iteración para el método de Jacobi viene dada por:

0

 x = (1 − y)/3

y0 = (2 − x − z)/3
 z0 = (3 − y − t)/3


 0
t = (4 − z)/3
y para el método de Gauss–Seidel:

x 00 = (1 − y)/3



y00 = (2 − x 00 − z)/3

z00 = (3 − y00 − t)/3


 00
t = (4 − z00 )/3
(sustituimos cada aproximación anterior por la aproximación nueva, siempre cuando ya está calculada).
Jacobi y Gauss–Seidel: ¿Cuando parar ?
Un posible criterio para parar las iteraciones es cuando una
aproximación difiere poco de la aproximación anterior. Es decir, nos
fijamos una tolerancia ε (por ejemplo 10−14 ) y paramos los cálculos
cuando obtenemos || X 0 − X || ≤ ε.
Jacobi y Gauss–Seidel: presentación con matrices
Las formulas de los métodos de Jacobi y de Gauss–Seidel para
calcular una nueva aproximación X 0 a partir de la aproximación
4
métodos iterativos para la resolución de sistemas
anterior X, son de la forma:
X 0 = MX + c
La matriz M se llama matriz del método.
Por ejemplo, para
(
4x + 2y = 3
x + 3y = 4
la iteración de Jacobi se presenta como
(
" # "
#" # "
#
x 0 = (3 − 2y)/4
x0
0
−2/4 x
3/4
,
es decir
=
+
y0 = (4 − x )/3
y0
−1/3
0
y
4/3
"
#
0
−2/4
La matriz del método de Jacobi es J =
.
−1/3
0
Para el mismo sistema, la iteración de Gauss–Seidel se presenta
como:
(
x 00 = (3 − 2y)/4
00
00
y = (4 − x )/3 = (4 − (3 − 2y)/4)/3 = 13/12 + y/6
La
del método de Gauss–Seidel para este sistema es G =
" matriz #
0 −2/4
.
0 1/6
De manera general, dado un sistema AX = b, en el método
de Jacobi cuando despejemos en cada fila i la i–esima variable,
corresponde del punto de vista matricial a separar la parte diagonal
D de A del resto de la matriz: A = D − ( A − D ) y reescribir el
sistema como
DX = −( A − D ) X + b
Por ejemplo
"
1
para A =
3
"
#
1
2
, se tiene D =
0
4
0
4
#
"
0
y A−D =
3
2
0
#
Luego se definir a partir de esta relación la iteración por
DX 0 = −( A − D ) X + b
Por tanto
X 0 = − D −1 ( A − D ) X + D −1 b
Vemos que la matriz del método de Jacobi es D −1 ( A − D ).
Para el método de Gauss–Seidel, partimos la matriz entre: su
parte triangular inferior L (incluyendo la diagonal); y el resto A − L
(es por tanto la parte triangular superior estricta) Por ejemplo
#
"
#
"
#
"
0 2
1 2
1 0
para A =
, se tiene L =
, A−L =
3 4
3 4
0 0
Entonces la iteración está definida por
LX 00 = −( A − L) X + b
y por tanto
X 00 = − L−1 ( A − L) X + L−1 b.
La matriz del método de Gauss–Seidel es por tanto − L−1 ( A − L).
5
métodos iterativos para la resolución de sistemas
Jacobi y Gauss–Seidel: garantías de convergencia
Una matriz A = ( ai,j ) tiene diagonal estrictamente dominante
por filas si para cada fila, el coeficiente de la diagonal es estrictamente superior a la suma de los valores absolutos de los otros
coeficientes de la fila:
| ai,i | > ∑ | ai,j |.
j:j6=i
Teorema 1. Si A es una matriz con diagonal estrictamente dominante
por filas entonces los métodos de Jacobi y de Gauss–Seidel convergen para
todos los sistemas AX = b.
Teorema 2. Si A es una matriz simétrica definida positiva (es decir: que
admite una descomposición de Cholesky) entonces el método de Gauss–
Seidel converge para todos los sistemas AX = b.
Estas condiciones son solamente condiciones suficientes. La convergencia puede darse incluso si no se cumplen estas condiciones.
El criterio más general es el siguiente:
Teorema 3. Sea una sucesión de vectores definida por una iteración
X 0 = Mx + c. Sea ρ el radio espectral de M, es decir el mayor valor
absoluto de los autovalores de M.
Si ρ < 1 entonces la sucesión converge, cualquier sea el valor inicialSi ρ >= 1 entonces la sucesión diverge, para casi cualquier valor
inicial.
Pero claro, en general, determinar el radio espectral es mucho
más complicado que resolver el sistema estudiado.
El método SOR (Successive OverRelaxation)
El método SOR es una mejora (en termino de velocidad de convergencia) de los métodos de Jacobi y Gauss–Seidel.
Dado un sistema AX = b, definimos a partir de cualquiera
aproximación X de la solución una nueva aproximación X 00 por el
método de Gauss–Seidel.
b definida por
En el método de SOR, la nueva aproximación es X
la relación:
b − X = ω ( X 00 − X ).
X
dónde ω es un parametro ajustable.
Ejemplo 1. Consideramos otra vez el sistema
(
4x + 2y = 3
x + 3y = 4
En el método de Gauss–Seidel, cada aproximación ( x 00 , y00 ) está
obtenida a partir de la anterior por
(
x 00 = (3 − 2y)/4
y00 = (4 − x 00 )/3
6
métodos iterativos para la resolución de sistemas
En el método SOR con parametro ω, la nueva aproximación es
( xb; yb) con
(
xb − x = ω ( x 00 − x )
yb − y = ω (y00 − y)
Por tanto:
Más explícitamente:
(
(
xb = (1 − ω ) x + ωx 00
yb = (1 − ω )y + ωy00
xb = (1 − ω ) x + ω (3 − 2 y)/4
yb = (1 − ω )y + ω (4 − x 00 )/3
Por ejemplo, para ω = 1, 5, la iteración será dada por:
(
xb = −0,5x + 1,5(3 − 2 y)/4
yb = −0,5y + 1,5(4 − x 00 )/3
Explicaremos este método con más detalle en el próximo episodio.
7