[Link] Técnicas iterativas para resolver sistemas lineales — Cap.

XVII

CAPITULO XVII. TECNICAS ITERATIVAS PARA RESOLVER

SISTEMAS LINEALES

1. INTRODUCCION Y METODO
Una técnica iterativa para resolver un sistema lineal A x = b de n × n empieza con
una aproximación inicial x(0) a la solución x, y genera una sucesión de vectores {x(k) }∞
k=0
que converge a x. La mayorı́a de estas técnicas iterativas involucran un proceso que
convierte el sistema A x = b en un sistema equivalente de la forma x = T x + c para
alguna matriz T de n × n y un vector c. Ya seleccionado el vector inicial x(0) la sucesión
de vectores de solución aproximada se genera calculando
x(k) = T x(k−1) + c (XV II.1)
para cada k = 1, 2, 3, . . .. Este tipo de procedimiento nos recuerda a la iteración del punto
fijo estudiada en la tercera parte.
Las técnicas iterativas se emplean raras veces para resolver sistemas lineales de di-
mensión pequeña ya que el tiempo requerido para lograr una precisión suficiente excede
al de las técnicas directas como el método de eliminación Gaussiana. Sin embargo, para
sistemas grandes con un gran porcentaje de ceros, estas técnicas son eficientes en términos
de almacenamiento en la computadora y del tiempo requerido. Los sistemas de este tipo
surgen frecuentemente en la solución numérica de problemas de valores en la frontera y
de ecuaciones diferenciales parciales.
Ejemplo 1.
El sistema lineal A x = b dado por
E1 : 10 x1 − x2 + 2 x3 = 6,
E2 : − x1 + 11 x2 − x3 + 3 x4 = 25 ,
E3 : 2 x1 − x2 + 10 x3 − x4 = −11 ,
E4 : 3 x2 − x3 + 8 x4 = 15 ,
tiene por solución a x = (1, 2, −1, 1)t . Para convertir A x = b a la forma x = T x + c,
resolvemos la ecuación Ei para cada i = 1, 2, 3, 4, obteniendo:
1
x1 = 10 x2 − 15 x3 + 3
5 ,

1 1 3 25
x2 = 11 x1 + 11 x3 − 11 x4 + 11 ,

x3 = − 15 x1 + 1
10 x2 + 1
10 x4 − 11
10 ,

3 1 15
x4 = − 8 x2 + 8 x3 + 8 .
En este ejemplo,
 1   
0 10 − 15 0 3
5
   
 1 1 3   25 
 11 0 11 − 11   11 
   
T =  y c=  .
 1 1 1   11 
 −5 10 0 10   − 10 
   
0 − 38 1
8 0 15
8

229
[Link] Técnicas iterativas para resolver sistemas lineales — Cap. XVII

Como una aproximación inicial tomemos a x(0) = (0, 0, 0, 0)t y generemos x(1) mediante:
(1) 1 (0) 1 (0) 3
x1 = 10 x2 − 5 x3 + 5 = 0.6000 ,
(1) 1 (0) 1 (0) 3 (0) 25
x2 = 11 x1 + 11 x3 − 11 x4 + 11 = 2.2727 ,
(1) (0) 1 (0) 1 (0)
x3 = − 15 x1 + 10 x2 + 10 x4 − 11
10 = −1.1000 ,
(1) 3 (0) 1 (0) 15
x4 = − 8 x2 + 8 x3 + 8 = 1.8750 .

(k) (k) (k) (k)

Las iteraciones adicionales x(k) = (x1 , x2 , x3 , x4 )t , se generan de manera similar y
se presentan en la tabla siguiente.
Tabla 1
(k) (k) (k) (k)
k x1 x2 x3 x4
0 0.0000 0.0000 0.0000 0.0000
1 0.6000 2.2727 −1.1000 1.8750
2 1.0473 1.7159 −0.80523 0.88524
3 0.93264 2.0533 −1.0493 1.1309
4 1.0152 1.9537 −0.96811 0.97385
5 0.98899 2.0114 −1.0103 1.0213
6 1.0032 1.9923 −0.99453 0.99444
7 0.99814 2.0023 −1.0020 1.0036
8 1.0006 1.9987 −0.99904 0.99889
9 0.99968 2.0004 −1.0004 1.0006
10 1.0001 1.9998 −0.99984 0.99980

La decisión de parar después de diez iteraciones está basada en el hecho de que

||x(10) − x(9) ||∞ 8.0 × 10−4

(10)
= < 10−3 .
||x ||∞ 1.9998

En realidad, ||x(10) − x||∞ = 0.0002.

El método del ejemplo 1 se llama método iterativo de Jacobi. Este consiste en
resolver la i−ésima ecuación de A x = b para xi para obtener, siempre y cuando aii 6= 0,
que
n ³
X aij xj ´ bi
xi = − + para i = 1, 2, . . . , n (XV II.2)
j=1
a ii aii
j6=i

(k)
y generar cada xi de las componentes de x(k−1) para k ≥ 1 con
n
(k) 1 X (k−1)
xi = [ (−aij xj ) + bi ] para i = 1, 2, . . . , n . (XV II.3)
aii j=1
j6=i

El método puede escribirse en la forma x(k) = T x(k−1) + c dividiendo a A en su parte

diagonal y no-diagonal. Para ver esto, sean D la matriz diagonal cuya diagonal es la

230
[Link] Técnicas iterativas para resolver sistemas lineales — Cap. XVII

misma que la diagonal de A, −L la parte triangular estrictamente inferior de A, y −U la

parte triangular estrictamente superior de A. Con esta notación, se separa en
   
a11 a12 . . . a1n a11 0 . . . 0
a a22 . . . a2n   0 a22 . . . 0 
A =  21 = +
... ... ... ... ... ... ... ...
an1 an2 . . . ann 0 0 . . . ann
   
0 0 ... 0 0 −a12 . . . −a1n
 −a21 0 ... 0   0 0 ... ... 
− − =
... ... ... ... . . . . . . . . . −an−1,n
−an1 . . . −an,n−1 0 0 0 ... 0
= D − L − U .

La ecuación A x = b ó (D − L − U ) x = b se transforma entonces en D x = (L + U ) x + b,

y finalmente
x = D−1 (L + U ) x + D−1 b . (XV II.4)

Esto da lugar a la forma matricial de la técnica iterativa de Jacobi:

x(k) = D−1 (L + U ) x(k−1) + D−1 b , k = 1, 2, . . . . (XV II.5)

En la práctica, la ecuación (XV II.3) es la que se usa para los cálculos, reservando a la
ecuación (XV II.5) para propósitos teóricos.

2. LOS ALGORITMOS DE JACOBI Y DE GAUSS-SEIDEL

Para resumir el método iterativo de Jacobi, presentamos el siguiente algoritmo:

Algoritmo iterativo de Jacobi.
==================================================
Para resolver el sistema lineal A x = b con una aproximación inicial dada x(0) .
Entrada: número de incógnitas y de ecuaciones n; las componentes de la matriz A = (aij )
donde 1 ≤ i, j ≤ n; las componentes bi , con 1 ≤ i ≤ n, del término no homogéneo b; las
componentes XOi , con 1 ≤ i ≤ n, de la aproximación inicial XO = x(0) ; la tolerancia
TOL; el número máximo de iteraciones N0 .
Salida: solución aproximada x1 , x2 , . . . , xn ó mensaje de que el número de iteraciones
fue excedido.
Paso 1: Tomar k = 1.
Paso 2: Mientras que k ≤ N0 seguir los pasos 3–6.
Paso 3: Para i = 1, 2, . . . , n tomar

Xn
1
xi = [− (aij XOj ) + bi ] .
aii j=1
j6=i

Paso 4: Si ||x − XO|| < T OL entonces SALIDA (x1 , x2 , . . . , xn );

(procedimiento completado satisfactoriamente) PARAR.

231
[Link] Técnicas iterativas para resolver sistemas lineales — Cap. XVII

Paso 5: Tomar k = k + 1.
Paso 6: Para i = 1, 2, . . . , n tomar XOi = xi .
Paso 7: SALIDA (número máximo de iteraciones excedido);
(procedimiento completado sin éxito) PARAR.
==================================================
El paso 3 del algoritmo requiere que aii 6= 0 para cada i = 1, 2, . . . , n. Si éste no es
el caso, se puede realizar un reordenamiento de las ecuaciones para que ningún aii = 0,
a menos que el sistema sea singular. Se sugiere que las ecuaciones sean arregladas de tal
manera que aii sea lo más grande posible para acelerar la convergencia.
En el paso 4, el criterio de paro ha sido ||x − XO|| < T OL; otro criterio de paro es
iterar hasta que
||x(k) − x(k−1) ||
||x(k) ||
sea menor que alguna tolerancia predeterminada ε > 0. Para este propósito, se puede
usar cualquier norma conveniente; la que más se usa es la norma l∞ .
Un análisis de la ecuación (XV II.3) sugiere una posible mejora en el algoritmo
(k)
iterativo de Jacobi. Para calcular xi , se usan las componentes de x(k−1) . Como para i >
(k) (k) (k)
1, x1 , x2 , . . ., xi−1 ya han sido calculadas y supuestamente son mejores aproximaciones
(k) (k) (k) (k)
a la solución real x1 , x2 , . . ., xi−1 que x1 , x2 , . . ., xi−1 , parece razonable calcular xi
usando los valores calculados más recientemente; es decir,
i−1
X Xn
(k) 1 (k) (k−1)
xi = [− (aij xj ) − (aij xj ) + bi ] , (XV II.6)
aii j=1 j=i+1

para cada i = 1, 2, . . . , n en vez de la ecuación (XV II.3).

Ejemplo 2.
El sistema lineal A x = b dado por

E1 : 10 x1 − x2 + 2 x3 = 6,
E2 : − x1 + 11 x2 − x3 + 3 x4 = 25 ,
E3 : 2 x1 − x2 + 10 x3 − x4 = −11 ,
E4 : 3 x2 − x3 + 8 x4 = 15 ,

fue resuelto en el ejemplo 1 con el método iterativo de Jacobi. Incorporando la ecuación

(XV II.6) en el algoritmo iterativo de Jacobi, se obtienen las ecuaciones que se usarán
para cada k = 1, 2, . . .:
(k) 1 (k−1) 1 (k−1) 3
x1 = 10 x2 − 5 x3 + 5 ,
(k) 1 (k) 1 (k−1) 3 (k−1) 25
x2 = 11 x1 + 11 x3 − 11 x4 + 11 ,
(k) (k) 1 (k) 1 (k−1)
x3 = − 15 x1 + 10 x2 + 10 x4 − 11
10 ,
(k) 3 (k) 1 (k) 15
x4 = − 8 x2 + 8 x3 + 8 .

232
 [Link] Técnicas iterativas para resolver sistemas lineales — Cap. XVII

Tomando x(0) = (0, 0, 0, 0)t , generamos los vectores iterados de la tabla 2.

Tabla 2
(k) (k) (k) (k)
k x1 x2 x3 x4
0 0.0000 0.0000 0.0000 0.0000
1 0.6000 2.3273 −0.98727 0.87885
2 1.0302 2.0369 −1.0145 0.98435
3 1.0066 2.0035 −1.0025 0.99838
4 1.0009 2.0003 −1.0003 0.99985
5 1.0001 2.0000 −1.0000 1.0000

Ya que
||x(5) − x(4) ||∞ 0.0008
(4)
= = 4 × 10−4 ,
||x ||∞ 2.000
se acepta x(5) como una aproximación razonable a la solución. Es interesante notar que el
método de Jacobi en el ejemplo 1 requiere el doble de iteraciones para la misma precisión.
La técnica presentada en el ejemplo 2 se llama método iterativo de Gauss-Seidel.
Para escribir este método en la forma matricial (XV II.1) se multiplican ambos lados de
la ecuación (XV II.6) por aii y se recolectan todos los k−ésimos términos iterados para
dar
(k) (k) (k) (k−1)
ai1 x1 + ai2 x3 + . . . + aii xi = −ai,i+1 xi+1 − . . . − ain x(k−1)
n + bi ,

para cada i = 1, 2, . . . , n. Escribiendo las n ecuaciones tenemos:

(k) (k−1) (k−1)
a11 x1 = −a12 x2 − a13 x3 − . . . − a1n x(k−1)
n + b1 ,
(k) (k) (k−1)
a21 x1 + a22 x2 = −a23 x3 . . . − a2n x(k−1)
n + b2 ,

... ... ... ... ... ...

(k) (k)
an1 x1 + an2 x2 + . . . + ann x(k)
n = bn ,

y se sigue que, en forma matricial, el método de Gauss-Seidel puede ser representado

como (D − L) x(k) = U x(k−1) + b, ó

x(k) = (D − L)−1 U x(k−1) + (D − L)−1 b . (XV II.7)

Para que la matriz triangular inferior (D − L) sea no singular, es necesario y suficiente

que aii 6= 0 para cada i = 1, 2, . . . , n.
Para resumir el método iterativo de Gauss-Seidel, presentamos el siguiente algoritmo:
Algoritmo iterativo de Gauss-Seidel.
=====================================================
Para resolver el sistema lineal A x = b con una aproximación inicial dada x(0) .
Entrada: número de incógnitas y de ecuaciones n; las componentes de la matriz A = (aij )
donde 1 ≤ i, j ≤ n; las componentes bi , con 1 ≤ i ≤ n, del término no homogéneo b; las

233
[Link] Técnicas iterativas para resolver sistemas lineales — Cap. XVII

componentes XOi , con 1 ≤ i ≤ n, de la aproximación inicial XO = x(0) ; la tolerancia

TOL; el número máximo de iteraciones N0 .
Salida: solución aproximada x1 , x2 , . . . , xn ó mensaje de que el número de iteraciones
fue excedido.
Paso 1: Tomar k = 1.
Paso 2: Mientras que k ≤ N0 seguir los pasos 3–6.
Paso 3: Para i = 1, 2, . . . , n tomar

i−1
X Xn
1
xi = [− (aij xj ) − (aij XOj ) + bi ] .
aii j=1 j=i+1

Paso 4: Si ||x − XO|| < T OL entonces SALIDA (x1 , x2 , . . . , xn );

(procedimiento completado satisfactoriamente) PARAR.
Paso 5: Tomar k = k + 1.
Paso 6: Para i = 1, 2, . . . , n tomar XOi = xi .
Paso 7: SALIDA (número máximo de iteraciones excedido);
(procedimiento completado sin éxito) PARAR.
==================================================
Los resultados de los ejemplos 1 y 2 parecen implicar que el método de Gauss-Seidel
es superior al método de Jacobi. Este es generalmente cierto, pero no siempre. En
realidad, hay sistemas lineales para los cuales el método de Jacobi converge y el método
de Gauss-Seidel no, y viceversa.

3. CONVERGENCIA DE LOS PROCESOS ITERATIVOS

Para estudiar la convergencia de las técnicas generales de iteración, consideramos la

fórmula (XV II.1)
x(k) = T x(k−1) + c

para cada k = 1, 2, . . ., donde x(0) es arbitrario. Este estudio requerirá del siguiente lema:
Lema XVII.1
Si el radio espectral ρ(T ) satisface que ρ(T ) < 1, ó si la norma de la matriz T satisface
que ||T || < 1, entonces (I − T )−1 existe y

(I − T )−1 = I + T + T 2 + . . . .

Teorema XVII.2
Para cualquier x(0) ∈ Rn , la sucesión {x(k) }∞
k=0 definida por (XV II.1)

x(k) = T x(k−1) + c

para cada k ≥ 1 y c 6= 0, converge a la solución única de x = T x + c si y sólo si ρ(T ) < 1.

234
[Link] Técnicas iterativas para resolver sistemas lineales — Cap. XVII

Demostración: de la ecuación (XV II.1), se tiene que

x(k) = T x(k−1) + c =
= T (T x(k−2) + c) + c =
= T 2 x(k−2) + (T + I) c =
...
= T k x(0) + (T k−1 + . . . + T + I) c .

Suponiendo que ρ(T ) < 1, podemos usar el Teorema XIII.15 y el Lema XVII.1 para
obtener
(k) k (0)
¡ k−1
X ¢
lim x = lim T x + lim Tj c
k→∞ k→∞ k→∞
j=0
(0) −1
= 0·x + (I − T ) c = (I − T )−1 c .

De (XV II.1) x = lim x(k) = (I − T )−1 c será la solución única de x = T x + c.

k→∞
Para probar el recı́proco, sea {x(k) }∞
k=0 convergente a x para cualquier x
(0)
. De la
ecuación (XV II.1) sigue que x = T x + c, ası́ que para cada k,

x − x(k) = T (x − x(k−1) ) = . . . = T k (x − x(0) ) .

Por lo tanto, para cualquier vector x(0) ,

lim T k (x − x(0) ) = lim x − x(k) = 0 .

k→∞ k→∞

Consecuentemente, si z es un vector arbitrario y x(0) = x − z, entonces

lim T k z = lim T k [x − (x − z)] = 0 ,

k→∞ k→∞

lo cual, por el Teorema XIII.15, implica que ρ(T ) < 1. c.q.d.

Un Teorema parecido nos dará condiciones de suficiencia para la convergencia de los
procesos de iteración usando las normas en lugar del radio espectral.
Teorema XVII.3
Si ||T || < 1, para cualquier norma matricial natural, entonces la sucesión definida en
la ecuación (XV II.1), {x(k) }∞
k=0 , converge para cualquier x
(0)
∈ Rn , a un vector x ∈ Rn ,
y se satisfacen las siguientes cotas de error:

||x − x(k) || ≤ ||T ||k ||x(0) − x|| , (XV II.8)

y
||T ||k
||x − x(k) || ≤ ||x(1) − x(0) || . (XV II.9)
1 − ||T ||

235
[Link] Técnicas iterativas para resolver sistemas lineales — Cap. XVII

Demostración: comenzando con un vector arbitrario x(0) , formaremos una secuencia de

aproximaciones
x(1) = T x(0) + c ,

x(2) = T x(1) + c ,

... ... ... ...

x(k) = T x(k−1) + c ,

de donde
x(k) = T k x(0) + (T k−1 + . . . + T + I) c .

Como para ||T || < 1 tenemos ||T k || → 0 cuando k → ∞, se deduce que

∞
X
k 2 k−1
lim T = 0 y lim (I + T + T + . . . + T )= T k = (I − T )−1 .
k→∞ k→∞
k=0

Y por tanto, pasando al lı́mite cuando k → ∞, tenemos

x = lim x(k) = (I − T )−1 c .

k→∞

Esto prueba la convergencia del proceso iterativo. Además, tenemos (I − T ) x = c ó

x = T x + c, lo cual quiere decir que el vector x en el lı́mite es una solución del sistema.
Como la matriz (I − T ) no es singular, la solución x es única. Hemos ası́ demostrado la
primera parte del Teorema.
Demostramos ahora la cota de error (XV II.8). Supongamos que x(k+p) y x(k) son
dos aproximaciones de la solución del sistema lineal x = T x+c; de la ecuación (XV II.1),
tenemos:

||x(k+p) − x(k) || = ||T x(k+p−1) − T x(k−1) || = ||T (x(k+p−1) − x(k−1) )|| = . . .

= ||T k (x(p) − x(0) )|| ≤ ||T ||k ||x(p) − x(0) || .

Ahora pasando al lı́mite cuando p → ∞, obtenemos

lim ||x(k+p) − x(k) || ≤ lim ||T ||k ||x(p) − x(0) || = ||T ||k lim ||x(p) − x(0) ||
p→∞ p→∞ p→∞

y entonces
||x − x(k) || ≤ ||T ||k ||x − x(0) || ,

que es la cota de error (XV II.8)

Finalmente demostramos la cota de error (XV II.9). Como antes, supongamos que
(k+p)
x y x(k) son dos aproximaciones de la solución del sistema lineal x = T x + c.
Tenemos

||x(k+p) − x(k) || ≤ ||x(k+1) − x(k) || + ||x(k+2) − x(k+1) || + . . . + ||x(k+p) − x(k+p−1) || .

236
[Link] Técnicas iterativas para resolver sistemas lineales — Cap. XVII

Por lo visto antes:

||x(m+1) − x(m) || ≤ ||T || ||x(m) − x(m−1) || ≤ ||T ||m−k ||x(k+1) − x(k) || ,

para m > k ≥ 1. Entonces tenemos:

||x(p+k) − x(k) || ≤ ||x(k+1) − x(k) || + ||T || ||x(k+1) − x(k) || + . . . +

1
+ ||T ||p−1 ||x(k+1) − x(k) || ≤ ||x(k+1) − x(k) || ≤
1 − ||T ||
||T || ||T ||k
≤ ||x(k) − x(k−1) || ≤ . . . ≤ ||x(1) − x(0) || ,
1 − ||T || 1 − ||T ||

de donde se deduce la cota de error (XV II.9). c.q.d.

Notése que si en particular elegimos x(0) = c, entonces x(1) = T c + c y

||x(1) − x(0) || = ||T c|| ≤ ||T || ||c|| ,

y la cota (XV II.9) nos da:

||T ||k+1
||x − x(k) || ≤ ||c|| . (XV II.90 )
1 − ||T ||

Ejemplo 3.
Demostrar que el proceso de iteración de Jacobi es convergente para el sistema lineal
siguiente:
E1 : 10 x1 − x2 + 2 x3 − 3 x4 = 0,
E2 : x1 + 10 x2 − x3 + 2 x4 = 5,
E3 : 2 x1 + 3 x2 + 20 x3 − x4 = −10 ,
E4 : 3 x1 + 2 x2 + x3 + 20 x4 = 15 .
¿Cuántas iteraciones han de efectuarse para hallar las raı́ces del sistema con un error
menor de 10−4 ?
Reduciendo el sistema a la forma especial para la iteración de Jacobi, tenemos

x1 = 0.1 x2 − 0.2 x3 + 0.3 x4 ,

x2 = −0.1 x1 + 0.1 x3 − 0.2 x4 + 0.5 ,

x3 = −0.1 x1 − 0.15 x2 + 0.05 x4 − 0.5 ,

x4 = −0.15 x1 − 0.1 x2 − 0.05 x3 + 0.75 .

Entonces la matriz del sistema es:

 
0 0.1 −0.2 0.3
 −0.1 0 0.1 −0.2 
T =  .
−0.1 −0.15 0 0.05
−0.15 −0.1 −0.05 0

237
[Link] Técnicas iterativas para resolver sistemas lineales — Cap. XVII

Utilizando, por ejemplo, la norma l1 , tenemos:

||T ||1 = max{0.35, 0.35, 0.35, 0.55} = 0.55 < 1 .

En consecuencia el proceso de iteración para el sistema dado es convergente. Si conside-

ramos como aproximación inicial de la raı́z x el vector

x(0) = c = (0.0, 0.5, −0.5, 0.75)t ,

entonces
||c||1 = 0.0 + 0.5 + 0.5 + 0.75 = 1.75 .

Sea ahora k el número de iteraciones requeridas para conseguir la exactitud especificada.

Utilizando la fórmula (XV II.90 ), tenemos:

(k) ||T ||k+1 0.55k+1 × 1.75

||x − x ||1 ≤ 1
||c||1 = < 10−4 .
1 − ||T ||1 0.45

De aquı́,
45
0.55k+1 < 10−4
175
o sea
(k + 1) log10 0.55 < log10 45 − log10 175 − 4

−(k + 1) 0.25964 < 1.65321 − 2.24304 − 4 = −4.58983

y consecuentemente

4.58983
k+1> ≈ 17.7 =⇒ k > 16.7 .
0.25964
Podemos tomar k = 17. Notése que la estimación teórica del número de iteraciones
necesarias para asegurar la exactitud especificada es excesivamente alto. A menudo se
obtiene la exactitud deseada en un número menor de iteraciones.
Para aplicar los resultados de arriba a las técnicas iterativas de Jacobi o Gauss-
Seidel, necesitamos escribir las matrices de iteración del método de Jacobi, TJ , dadas en
(XV II.5) y del método de Gauss-Seidel, TGS , dadas en (XV II.7), como

TJ = D−1 (L + U ) y TGS = (D − L)−1 U .

De ser ρ(TJ ) ó ρ(TGS ) menores que uno, es claro que la sucesión {x(k) }∞
k=0 converge a
la solución x de A x = b. Por ejemplo, el esquema de Jacobi (ver ecuación (XV II.5))
tiene:
x(k) = D−1 (L + U ) x(k−1) + D−1 b ,

y si {x(k) }∞
k=0 converge a x, entonces

x = D−1 (L + U ) x + D−1 b .

238
[Link] Técnicas iterativas para resolver sistemas lineales — Cap. XVII

Esto implica que

D x = (L + U ) x + b y (D − L − U ) x = b .

Ya que D − L − U = A, luego x satisface A x = b. De manera parecida se procede con

el esquema de Gauss-Seidel dado por la ecuación (XV II.7).
Podemos dar ahora condiciones de suficiencia fáciles de verificar para la convergencia
de los métodos de Jacobi y de Gauss-Seidel.
Teorema XVII.4
Si A es una matriz estrictamente dominante diagonalmente, entonces, para cualquier
elección de x(0) ∈ Rn ambos métodos, el de Jacobi o el de Gauss-Seidel, dan lugar a
sucesiones {x(k) }∞
k=0 que convergen a la solución de A x = b.

La relación entre la rapidez de convergencia y el radio espectral de la matriz de

iteración T se puede ver de la desigualdad (XV II.8). Como (XV II.8) se satisface para
cualquier norma matricial natural se sigue, de la afirmación que siguió al Teorema XIII.14,
que
||x(k) − x|| ≈ ρ(T )k ||x(0) − x|| . (XV II.10)

Supongamos que ρ(T ) < 1 y que se va a usar x(0) = 0 en una técnica iterativa para
aproximar x con un error relativo máximo de 10−t . Por la estimación (XV II.10), el
error relativo después de k iteraciones es aproximadamente ρ(T )k , ası́ que se espera una
precisión de 10−t si
ρ(T )k ≤ 10−t ,

esto es, si
t
k≥ .
− log10 ρ(T )
Por lo tanto, es deseable escoger la técnica iterativa con el menor ρ(T ) < 1 para el sistema
particular A x = b.
En general no se conoce cuál de las dos técnicas, la de Jacobi o la de Gauss-Seidel,
debe usarse. Sin embargo, en un caso especial, sı́ se conoce la respuesta.
Teorema XVII.5 (Stein-Rosenberg)
Si aij ≤ 0 para cada i 6= j y aii > 0 para cada i = 1, 2, . . . , n, entonces se satisface
una y solamente una de las siguientes afirmaciones:
a) 0 < ρ(TGS ) < ρ(TJ ) < 1;
b) 1 < ρ(TJ ) < ρ(TGS );
c) ρ(TGS ) = ρ(TJ ) = 0;
d) ρ(TJ ) = ρ(TGS ) = 1;
Para el caso especial descrito en el Teorema XVII.5, vemos que cuando un método
converge, entonces ambos convergen, siendo el método de Gauss-Seidel más rápido que el
método de Jacobi.

239
[Link] Técnicas iterativas para resolver sistemas lineales — Cap. XVII

4. LOS METODOS DE RELAJACION

Como la razón de convergencia de un procedimiento depende del radio espectral de
la matriz asociada con el método, una manera de seleccionar un procedimiento que nos
lleve a una convergencia acelerada consiste en escoger un método cuya matriz asociada
tenga un radio espectral mı́nimo. Estos procedimientos nos llevan a los métodos de
relajación. Pero antes de formular la teorı́a de los métodos de relajación, veamos las ideas
fundamentales de la forma más simple. Supongamos que se dispone de un sistema de
ecuaciones lineales

E1 : a11 x1 + a12 x2 + . . . + a1n xn = b1 ,

E2 : a21 x1 + a22 x2 + . . . + a2n xn = b2 ,
(XV II.11)
... ... ... ... ... ...
En : an1 x1 + an2 x2 + . . . + ann xn = bn .

Transformaremos este sistema de la manera siguiente: pondremos los términos constantes

a la izquierda y dividiremos la primera ecuación por −a11 , la segunda por −a22 , etc.
Obtendremos entonces un sistema que está listo para la relajación:

E1 : − x1 + b12 x2 + . . . + b1n xn + c1 = 0 ,
E2 : b21 x1 − x2 + . . . + b2n xn + c2 = 0 ,
(XV II.12)
... ... ... ... ... ... ...
En : bn1 x1 + bn2 x2 + . . . − xn + cn = 0 ,

donde
aij bi
bij = − (i 6= j) y ci = . (XV II.13)
aii aii
(0) (0)
Supongamos que x(0) = (x1 , . . . , xn ) es la aproximación inicial a la solución del sistema
dado. Sustituyendo estos valores en el sistema tendremos los restos
n
X
(0) (0) (0) (1) (0)
R1 = c1 − x1 + b1j xj = x1 − x1 ,
j=2

... ...
... ...
X n
(0) (0) (0) (1) (0)
Rk = ck − xk + bkj xj = xk − xk ,
(XV II.14)
j=1
j6=k

... ... ... ...

n−1
X (0)
Rn(0) = cn − x(0)
n + bnj xj = x(1) (0)
n − xn .
j=1

(0) (0) (0)

Si damos un incremento δxs a una de las incógnitas xs , el resto correspondiente Rs
(0) (0)
quederá disminuido en δxs y todos los otros restos Ri (i 6= s) quedarán aumentados en

240
[Link] Técnicas iterativas para resolver sistemas lineales — Cap. XVII
(0) (1)
bis δxs . De este modo, para hacer que desaparezca el resto siguiente Ri es suficiente
(1) (1) (0)
dar a xs un incremento δxs = Rs y tendremos
(1) (0)
Rs(1) = 0 y Ri = Ri + bis δx(0)
s para i 6= s . (XV II.15)

Ası́ el método de relajación, en su forma más simple, consiste en reducir el resto

numéricamente más elevado a cero, en cada etapa, cambiando el valor del componente
apropiado de la aproximación. El proceso acaba cuando todos los restos del último sistema
transformado son iguales a cero con la exactitud requerida.
Ejemplo 4.
Vamos a resolver el sistema lineal A x = b dado por
E1 : 10 x1 − 2 x2 − 2 x3 = 6,
E2 : − x1 + 10 x2 − 2 x3 = 7,
E3 : − x1 − x2 + 10 x3 = 8,
con el método de relajación y aritmética de dos dı́gitos.
Reduzcamos el sistema a una forma conveniente para la relajación:
E1 : − x1 + 0.2 x2 + 0.2 x3 + 0.6 = 0 ,
E2 : 0.1 x1 − x2 + 0.2 x3 + 0.7 = 0 ,
E3 : 0.1 x1 + 0.1 x2 − x3 + 0.8 = 0 .
(0) (0) (0)
Eligiendo los valores x1 = x2 = x3 = 0 como aproximaciones iniciales de las raı́ces,
obtenemos los restos:
(0) (0) (0)
R1 = 0.60, R2 = 0.70, R3 = 0.80.
(0)
Siendo R3 = 0.80 el resto más grande, consideramos
(0)
δx3 = 0.80

de donde obtendremos los restos

(1) (0)
R1 = R1 + 0.2 × 0.8 = 0.60 + 0.16 = 0.76 ,
(1) (0)
R2 = R2 + 0.2 × 0.8 = 0.70 + 0.16 = 0.86 ,
(1) (0)
R3 = R3 − 0.8 = 0.80 − 0.80 = 0.0 .

Establezcamos ahora
(1)
δx2 = 0.86
y ası́ sucesivamente. Los resultados de los cálculos se dan en la tabla 3.
(k)
Sumando todos los incrementos δxi (i = 1, 2, 3; k = 0, 1, . . .), tendremos los valores
de las raı́ces:
x1 = 0.0 + 0.93 + 0.062 + 0.007 = 0.999 = 1.0 ,
x2 = 0.0 + 0.86 + 0.13 + 0.016 = 1.006 = 1.0 ,
x3 = 0.0 + 0.80 + 0.18 + 0.019 = 0.999 = 1.0 .

241
[Link] Técnicas iterativas para resolver sistemas lineales — Cap. XVII

En este caso el sistema se ha resuelto exactamente.

Tabla 3

k x1 R1 x2 R2 x3 R3
0 0.0 0.60 0.0 0.70 0.0 0.80
0.16 0.16 0.80 −0.80
0.76 0.86 0.0

1 0.17 0.86 −0.86 0.086

0.93 0.0 0.086

2 0.93 −0.93 0.093 0.093

0.0 0.093 0.18

3 0.036 0.036 0.18 −0.18

0.036 0.13 0.0

4 0.026 0.13 −0.13 0.013

0.062 0.0 0.013

5 0.062 −0.062 0.012 0.0062

0.0 0.012 0.019

6 0.0038 0.0038 0.019 −0.019

0.0038 0.016 0.0

7 0.0032 0.016 −0.016 0.0016

0.007 0.0 0.0016

8 0.007 −0.007 0.0007 0.0007

0.0 0.0007 0.0023
P
1.0 1.0 1.0

Vamos ahora a describir los métodos de relajación. Antes de describir un proce-

dimiento para seleccionar tales métodos, necesitamos introducir una manera nueva de
medir la cantidad por la cual una aproximación a la solución de un sistema lineal difiere
de la solución real del sistema. El método hace uso del denominado vector residual.

Definición. Si x̃ ∈ Rn es una aproximación a la solución del sistema lineal definido por

A x = b, el vector residual de x̃ con respecto a este sistema se define como r = b − A x̃.

En procedimientos como los métodos de Jacobi o de Gauss-Seidel se asocia un vector

residual con cada cálculo de una componente aproximada del vector solución. El objetivo
del método consiste en generar una sucesión de aproximaciones que hagan que los vectores
residuales asociados converjan a cero. Supongamos que tomamos

(k) (k) (k) (k)

ri = (r1i , r2i , . . . , rni )t

242
[Link] Técnicas iterativas para resolver sistemas lineales — Cap. XVII

para denotar al vector residual para el método de Gauss-Seidel correspondiente al vector

solución aproximado
(k) (k) (k) (k−1)
(x1 , x2 , . . . , xi−1 , xi , . . . , x(k−1)
n )t .

(k)
La m−ésima componente de ri es

i−1
X n
X
(k) (k) (k−1)
rmi = bm − amj xj − amj xj (XV II.16)
j=1 j=i

ó
i−1
X n
X
(k) (k) (k−1) (k−1)
rmi = bm − amj xj − amj xj − ami xi
j=1 j=i+1

(k)
para cada m = 1, 2, . . . , n. En particular, la i−ésima componente de ri es

i−1
X n
X
(k) (k) (k−1) (k−1)
rii = bi − aij xj − aij xj − aii xi ;
j=1 j=i+1

ası́ que
i−1
X n
X
(k−1) (k) (k) (k−1)
aii xi + rii = bi − aij xj − aij xj . (XV II.17)
j=1 j=i+1

(k)
Recuérdese, sin embargo, que en el método de Gauss-Seidel xi se escoge como

P
i−1
(k) P
n
(k−1)
− (aij xj ) − (aij xj ) + bi
(k) j=1 j=i+1
xi = , (XV II.6)
aii
(k−1) (k) (k)
ası́ que la ecuación (XV II.17) puede escribirse como aii xi + rii = aii xi ó
(k)
(k) (k−1) r
xi = xi + ii . (XV II.18)
aii

Podemos derivar otra conexión entre los vectores residuales y la técnica de Gauss-
(k)
Seidel. De (XV II.16), la i−ésima componente de ri+1 es

i
X n
X
(k) (k) (k−1)
ri,i+1 = bi − aij xj − aij xj
j=1 j=i+1
(XV II.19)
i−1
X n
X
(k) (k−1) (k)
= bi − aij xj − aij xj − aii xi .
j=1 j=i+1

(k)
La ecuación (XV II.6) implica que ri,i+1 = 0. Entonces, en cierto sentido, la técnica de
(k)
Gauss-Seidel está ideada para requerir que la i−ésima componente de ri+1 sea cero.

243
[Link] Técnicas iterativas para resolver sistemas lineales — Cap. XVII

Reducir una coordenada del vector residual a cero, sin embargo, no es necesariamente
(k)
la manera más eficiente de reducir la norma del vector ri+1 . En realidad, modificando el
procedimiento de Gauss-Seidel en la forma de la ecuación (XV II.18) a:
(k)
(k) (k−1) rii
xi = xi +ω (XV II.20)
aii

para ciertas elecciones de ω positivo nos llevará a una convergencia significativamente

más rápida.
Los métodos que emplean la ecuación (XV II.20) se conocen como métodos de
relajación. Para 0 < ω < 1, los procedimientos se llaman métodos de sub-relajación
y se pueden emplear para obtener la convergencia de algunos sistemas que no son conver-
gentes por el método de Gauss-Seidel. Para ω > 1, los procedimientos se llaman métodos
de sobre-relajación y se pueden usar para acelerar la convergencia de sistemas que son
convergentes por el método de Gauss-Seidel. Estos métodos se abrevian frecuentemente
como SOR (de Successive Over-Relaxation) y son particularmente útiles para re-
solver los sistemas lineales que aparecen en la solución numérica de ciertas ecuaciones
diferenciales parciales.
Antes de ilustrar las ventajas del método SOR notamos que usando la ecuación
(XV II.17), la ecuación (XV II.20) se puede reformular para propósitos de cómputo como

ωh i
i−1
X Xn
(k) (k−1) (k) (k−1)
xi = (1 − ω) xi + bi − aij xj − aij xj . (XV II.21)
aii j=1 j=i+1

Para determinar la forma matricial del método SOR reescribimos (XV II.21) como

i−1
X n
X
(k) (k) (k−1) (k−1)
aii xi +ω aij xj = (1 − ω) aii xi −ω aij xj + ω bi
j=1 j=i+1

ası́ que
(D − ω L) x(k) = [(1 − ω) D + ω U ] x(k−1) + ω b

ó
x(k) = (D − ω L)−1 [(1 − ω) D + ω U ] x(k−1) + ω (D − ω L)−1 b .

Algoritmo iterativo Successive Over-Relaxation (SOR).

==================================================
Para resolver el sistema lineal A x = b dados el parámetro ω y una aproximación inicial
x(0) .
Entrada: número de incógnitas y de ecuaciones n; las componentes de la matriz A = (aij )
donde 1 ≤ i, j ≤ n; las componentes bi , con 1 ≤ i ≤ n, del término no homogéneo b; las
componentes XOi , con 1 ≤ i ≤ n, de la aproximación inicial XO = x(0) ; el parámetro ω;
la tolerancia TOL; el número máximo de iteraciones N0 .
Salida: solución aproximada x1 , x2 , . . . , xn ó mensaje de que el número de iteraciones
fue excedido.

244
[Link] Técnicas iterativas para resolver sistemas lineales — Cap. XVII

Paso 1: Tomar k = 1.
Paso 2: Mientras que k ≤ N0 seguir los pasos 3–6.
Paso 3: Para i = 1, 2, . . . , n tomar
i−1 n
ω£ X X ¤
xi = (1 − ω) XOi + − (aij xj ) − (aij XOj ) + bi .
aii j=1 j=i+1

Paso 4: Si ||x − XO|| < T OL entonces SALIDA (x1 , x2 , . . . , xn );

(procedimiento completado satisfactoriamente) PARAR.
Paso 5: Tomar k = k + 1.
Paso 6: Para i = 1, 2, . . . , n tomar XOi = xi .
Paso 7: SALIDA (número máximo de iteraciones excedido);
(procedimiento completado sin éxito) PARAR.
==================================================
Ejemplo 5.
El sistema lineal A x = b dado por
E1 : 4 x1 + 3 x2 = 24 ,
E2 : 3 x1 + 4 x2 − x3 = 30 ,
E3 : − x2 + 4 x3 = −24 ,
t
tiene por solución x = (3, 4, −5) . Se usarán los métodos de Gauss-Seidel y el SOR con
ω = 1.25 para resolver este sistema usando x(0) = (1, 1, 1)t para ambos métodos. Las
ecuaciones para el método de Gauss-Seidel son
(k) (k−1)
x1 = − 0.75 x2 +6 ,
(k) (k) (k−1)
x2 = − 0.75 x1 + 0.25 x3 + 7.5 ,
(k) (k)
x3 = 0.25 x2 − 6 ,
para cada k = 1, 2, . . ., y las ecuaciones para el método SOR con ω = 1.25 son
(k) (k−1) (k−1)
x1 = − 0.25 x1 − 0.9375 x2 + 7.5 ,
(k) (k) (k−1) (k−1)
x2 = − 0.9375 x1 − 0.25 x2 + 0.3125 x3 + 9.375 ,
(k) (k) (k−1)
x3 = 0.3125 x2 − 0.25 x3 − 7.5 .
Las primeras siete iteraciones de cada método se muestran en las tablas 4 y 5.
Para obtener una precisión de siete lugares decimales el método de Gauss-Seidel
requiere de 34 iteraciones en contra de las 14 que se necesitan en el método de sobre-
relajación con ω = 1.25.
Tabla 4
(k) (k) (k)
k x1 x2 x3
0 1.000000 1.000000 1.000000
1 5.250000 3.812500 −5.046875
2 3.1406250 3.8828125 −5.0292969
3 3.0878906 3.9267578 −5.0183105
4 3.0549317 3.9542236 −5.0114441
5 3.0343323 3.9713898 −5.0071526
6 3.0214577 3.9821186 −5.0044703
7 3.0134111 3.9888241 −5.0027940

245
[Link] Técnicas iterativas para resolver sistemas lineales — Cap. XVII

Tabla 5
(k) (k) (k)
k x1 x2 x3
0 1.000000 1.000000 1.000000
1 6.312500 3.5195313 −6.6501465
2 2.6223144 3.9585266 −4.6004238
3 3.1333027 4.0102646 −5.0966864
4 2.9570513 4.0074838 −4.9734897
5 3.0037211 4.0029250 −5.0057135
6 2.9963275 4.0009263 −4.9982822
7 3.0000498 4.0002586 −5.0003486

Un problema que se presenta al usar el método SOR, es cómo escoger el valor apro-
piado de ω. Aún cuando no se conoce una respuesta completa a esta pregunta para un
sistema lineal general n × n, los siguientes resultados pueden usarse en ciertas situaciones.
Teorema XVII.6 (Kahan)
Si aii 6= 0 para cada i = 1, 2, . . . , n, entonces ρ(Tω ) ≥ |ω − 1|. Esto implica que
ρ(Tω ) < 1 sólo si 0 < ω < 2, donde Tω = (D − ω L)−1 [(1 − ω) D + ω U ] es la matriz de
iteración del método SOR.
Teorema XVII.7 (Ostrowski-Reich)
Si A es una matriz positiva definida y 0 < ω < 2, entonces el método SOR converge
para cualquier elección de la aproximación inicial x(0) del vector solución.
Teorema XVII.8
Si A es una matriz positiva definida y tridiagonal, entonces ρ(TGS ) = [ρ(TJ )]2 < 1,
la elección óptima de ω para el método SOR es
2
ω= p , (XV II.22)
1 + 1 − [ρ(TJ )]2

y con este valor de ω, ρ(Tω ) = ω − 1.

Ejemplo 6.
En el ejemplo 5 la matriz A estaba dada por
 
4 3 0
A = 3 4 −1  .
0 −1 4

Esta matriz es positiva definida y tridiagonal, ası́ que se aplica el Teorema XVII.8. Como

TJ = D−1 (L + U ) =
   
1/4 0 0 0 −3 0
=  0 1/4 0   −3 0 1 =
0 0 1/4 0 1 0
 
0 −0.75 0
=  −0.75 0 0.25  ,
0 0.25 0

246
[Link] Técnicas iterativas para resolver sistemas lineales — Cap. XVII

tenemos que  
−λ −0.75 0

TJ − λI = −0.75 −λ 0.25  ,
0 0.25 −λ
con lo que
det(TJ − λI) = −λ(λ2 − 0.625) .

Por lo tanto,
√
ρ(TJ ) = 0.625 ≈ 0.79057

y la ecuación (XV II.22) nos da

2 2 2
ω= p = p = √ ≈ 1.24 .
1+ 1 − ρ(TGS ) 1+ 1 − [ρ(TJ )]2 1+ 1 − 0.625

Esto explica la rápida convergencia obtenida usando ω = 1.25 en el ejemplo 5.

5. ELECCION DEL METODO PARA RESOLVER SISTEMAS LINEALES

Cuando el sistema lineal es lo suficientemente pequeño para que sea fácilmente aco-
modado en la memoria principal de un ordenador, es en general más eficaz usar una
técnica directa que minimice el efecto del error de redondeo. Especı́ficamente, es ade-
cuado el algoritmo de eliminación Gaussiana con pivoteo escalado de columna.
Los sistemas lineales grandes cuyos coeficientes son entradas básicamente de ceros y
que aparecen en patrones regulares se pueden resolver generalmente de una manera efi-
ciente usando un procedimiento iterativo como el discutido en este capı́tulo. Los sistemas
de este tipo aparecen naturalmente, por ejemplo, cuando se usan técnicas de diferencias
finitas para resolver problemas de valor en la frontera, una aplicación común en la solución
numérica de ecuaciones diferenciales parciales.

247
[Link] Técnicas iterativas para resolver sistemas lineales — Cap. XVII

EJERCICIOS.

1. Encontrar las dos primeras iteraciones del método de Jacobi, del método de
Gauss-Seidel y del método SOR para los siguientes sistemas lineales, usando
x(0) = 0.

a) 2 x1 − x2 + x3 = −1 ,
3 x1 + 3 x2 + 9 x3 = 0,
3 x1 + 3 x2 + 5 x3 = 4.

b) 2 x2 + 4 x3 = 0,
x1 − x2 − x3 = 0.375 ,
x1 − x2 + 2 x3 = 0.

c) 10 x1 − x2 = 9,
− x1 + 10 x2 − 2 x3 = 7,
− 2 x2 + 10 x3 = 6.

d) 2 x1 = 3,
x1 + 1.5 x2 = 4.5 ,
− 3 x2 + 0.5 x3 = −6.6 ,
2 x1 − 2 x2 + x3 + x4 = 0.8 .

e) 2 x1 − x2 + 10 x3 = −11 ,
3 x2 − x3 + 8 x4 = −11 ,
10 x1 − x2 + 2 x3 = 6,
− x1 + 11 x2 − x3 + 3 x4 = 25 .

f) 10 x1 − x2 + 2 x3 = 6,
− x1 + 11 x2 − x3 + 3 x4 = 25 ,
2 x1 − x2 + 10 x3 = −11 ,
3 x2 − x3 + 8 x4 = −11 .

g) 4 x1 − 2 x2 = 0,
−2 x1 + 5 x2 − x3 = 2,
− x2 + 4 x3 + 2 x4 = 3,
2 x3 + 3 x4 = −2 .

2. Aplicar, si es posible, los algoritmos iterativos de Jacobi, de Gauss-Seidel y el

algoritmo SOR, con ω = 1.2, para resolver los sistemas lineales del ejercicio
1. Usar T OL = 10−2 y el número máximo de iteraciones N = 25.

248