93

5 Métodos iterativos para resolver sistemas de ecuaciones lineales

La resolución de sistemas de ecuaciones lineales también puede hacerse con fórmulas iterativas
que permiten acercarse a la respuesta mediante aproximaciones sucesivas, sin embargo desde el
punto de vista práctico es preferible usar métodos directos que con el soporte computacional actual
resuelven grandes sistemas en forma eficiente y con mucha precisión, a diferencia de los sistemas
de ecuaciones no-lineales cuya solución no se puede obtener mediante métodos directos.

Las fórmulas iterativas no siempre convergen, su análisis puede ser complicado, la convergencia
es de primer orden y se requiere elegir algún vector inicial para comenzar el proceso iterativo. En la
época actual el estudio de estos métodos iterativos se puede considerar principalmente como de
interés teórico matemático, excepto para resolver grandes sistemas de ecuaciones lineales con
matrices esparcidas y cuya convergencia se puede determinar fácilmente.
Para definir un método iterativo, se expresa el sistema AX = B en la forma X = G(X) con el mismo
fundamento descrito en el método del Punto Fijo para resolver ecuaciones no lineales.

5.1 Método de Jacobi

5.1.1 Formulación matemática
Dado un sistema de ecuaciones lineales
a1,1x1 + a1,2x2 + ... + a1,nxn = b1
a2,1x1 + a2,2x2 + ... + a2,nxn = b2
. . . .
an,1x1 + an,2x2 + ... + an,nxn = bn

Expresado abreviadamente en notación matricial: AX = B

Una forma simple para obtener la forma X = G(X) consiste en re-escribir el sistema despejando de
la ecuación i la variable xi a condición de que ai,i sea diferente de 0:

x1 = 1/a1,1 (b1 – a1,2x2 – a1,3x3 - ... – a1,nxn)

x2 = 1/a2,2 (b2 – a2,1x1 – a2,3x3 - ... – a2,nxn)
. . . .
xn = 1/an,n (bn – an,1x1 – an,2x2 - ... – an,n-1xn-1)

El cual puede escribirse con la notación de sumatoria:

i −1 n n
1 1
xi = (bi - ∑ ai,j x j − ∑ ai,j x j ) = (bi - ∑ ai,j x j ); i = 1, 2, ..., n;
ai,i j= 1 j= i + 1 ai,i =j 1,j≠i
El sistema está en la forma X = G(X) la cual sugiere su uso iterativo.

Utilizamos un índice para indicar iteración:

X(k+1) = G(X(k)), k=0, 1, 2, .... (iteraciones)

Fórmula iterativa de Jacobi:

n
1
xi(k + 1) = (bi - ∑ ai,j x (k)
j ); i = 1, 2, ..., n; k = 0, 1, 2, ...
ai,i =j 1,j≠i

X(0) es el vector inicial. A partir de este vector se obtienen sucesivamente los vectores X(1), X(2), ...
Si el método converge, entonces X(k) tiende a la solución X a medida que k crece:
X (k) → X
k →∞
 94

Ejemplo. Dadas las ecuaciones:

5x1 – 3x2 + x3 = 5
2x1 + 4x2 – x3 = 6
2x1 – 3x2 + 8x3 = 4

Formule un sistema iterativo con el método de Jacobi:

x1(k+1) = 1/5 (5 + 3x2(k) – x3(k))
x2(k+1) = 1/4 (6 – 2x1(k) + x3(k))
x3(k+1) = 1/8 (4 – 2x1(k) + 3x2(k)) k = 0, 1, 2, ...

Realice dos iteraciones, comenzando con los valores iniciales: x1(0) = x2(0) = x3(0) = 1

k=0: x1(1) = 1/5 (5 + 3x2(0) – x3(0)) = 1/5 (5 + 3(1) – (1)) = 1/5 (7) = 1.4
x2(1) = 1/4 (6 – 2x1(0) + x3(0)) = 1/4 (6 – 2(1) + (1)) = 1/4 (5) = 1.25
x3(1) = 1/8 (4 – 2x1(0) + 3x2(0)) = 1/8 (4 – 2(1) + 3(1)) = 1/8 (5) = 0.625
k=1:
x1(2) = 1/5 (5 + 3x2(1) – x3(1)) = 1/5 (5 + 3(1.25) – (0.625)) = 1.6250
x2(2) = 1/4 (6 – 2x1(1) + x3(1)) = 1/4 (6 – 2(1.4) + (0.625)) = 0.9563
x3(2) = 1/8 (4 – 2x1(1) + 3x2(1)) = 1/8 (4 – 2(1.4) + 3(1.25)) = 0.6188

5.1.2 Manejo computacional de la fórmula de Jacobi

La siguiente función en MATLAB recibe un vector X y entrega un nuevo vector X calculado en
cada iteración

function x = jacobi(a,b,x)
n=length(x);
t=x; % t es asignado con el vector X ingresado
for i=1:n
s=a(i,1:i-1)*t(1:i-1)+a(i,i+1:n)*t(i+1:n);
x(i) = (b(i) - s)/a(i,i);
end

Uso de la función Jacobi para el ejemplo anterior:

>> a= [5, -3, 1; 2, 4, -1;2, -3, 8];

>> b= [5; 6; 4];
>> x= [1; 1; 1]; Vector inicial
>> x=jacobi(a,b,x) Repita este comando y observe la convergencia
x=
1.4000
1.2500
0.6250
>> x=jacobi(a,b,x)
x=
1.6250
0.9563
0.6188
>> x=jacobi(a,b,x)
x=
1.4500
0.8422
0.4523 etc.
 95

5.1.3 Algoritmo de Jacobi

a: matriz de coeficientes, b: vector de constantes del sistema de n ecuaciones lineales

e: estimación del error para la solución, m: máximo de iteraciones permitidas
x: vector inicial para la solución, k es el conteo de iteraciones
t←x
Para k = 1, 2, . . ., m
Calcular el vector x con la fórmula de Jacobi
Si la norma del vector x – t es menor que e
x es el vector solución con precisión e
Terminar
sino
t←x
fin
fin
Al exceder el máximo de iteraciones entregar un vector x nulo

5.1.4 Instrumentación computacional del método de Jacobi

La siguiente función en MATLAB recibe la matriz de coeficientes a y el vector de constantes b de
un sistema lineal. Adicionalmente recibe un vector inicial x, la estimación del error e y el máximo
de iteraciones permitidas m. La función entrega el vector x calculado y el número de iteraciones
realizadas k. Si el método no converge, x contendrá un vector nulo y el número de iteraciones k
será igual al máximo m.

function [x,k] = jacobim(a,b,x,e,m)

n=length(x);
for k=1:m
t=x;
for i=1:n
s=a(i,1:i-1)*t(1:i-1)+a(i,i+1:n)*t(i+1:n);
x(i) = (b(i) - s)/a(i,i);
end
if norm((x-t),inf)<e
return
end
end
x=[ ];
k=m;

Ejemplo
Use la función jacobim para encontrar el vector solución del ejemplo anterior con una precisión
de 0.0001 Determine cuántas iteraciones se realizaron. Comenzar con el vector inicial x=[1; 1; 1]

>> a=[5, -3, 1; 2, 4, -1; 2, -3, 8];

>> b=[5; 6; 4];
>> x=[1; 1; 1];
>> [x,k]=jacobim(a,b,x,0.0001,20)
x=
1.4432
0.8973
0.4757
k=
12
 96

5.1.5 Forma matricial del método de Jacobi

Dado el sistema de ecuaciones lineales
AX = B
La matriz A se re-escribe como la suma de tres matrices:
A=L+D+U
D es una matriz diagonal con elementos iguales a los elementos de la diagonal principal de A
L es una matriz triangular inferior con ceros en la diagonal principal y los otros elementos iguales a
los elementos respectivos de la matriz A
U es una matriz triangular superior con ceros en la diagonal principal y los otros elementos iguales
a los elementos respectivos de la matriz A
En forma desarrollada:

 a1,1 a1,2 ... a1,n   0 0 ... 0   a1,1 0 0 0   0 a1,2 ... a1,n 

a   0   0 a 2,2 0 0   0 0
2,1 a 2,2 ... a2,n  a2,1 0 ... ... a2,n 
A= =L + D + U =  + +
 ... ... ... ...   ... ... ... ...  ... ... ... ...  ... ... ... ... 
       
 an,1 an,2 ... an,n   an,1 an,2 ... 0   0 0 0 an,n   0 0 0 0 
Sustituyendo en la ecuación:
(L + D + U)X = B
LX + DX + UX = B
X = D-1B - D-1LX - D-1UX, Siempre que D-1 exista

Ecuación recurrente del método de Jacobi según el método del Punto Fijo X = G(X)
X = D-1B - D-1 (L + U)X
En donde
 1
D −1 =  
 ai,i  nxn
Ecuación recurrente desarrollada
 x1   b1 / a1,1   0 a1,2 / a1,1 a1,3 / a1,1 ... a1,n / a1,1   x1 
  b / a   a / a 0 a2,3 / a2,2 ... a2,n / a2,2   x 2 
 x2 
=  2 2,2  −  2,1 2,2  
 :   :   : : : ... :  : 
      
 x2   n n,n   n,1 an,n
b / a a / an,2 / an,n an,3 / an,n ... 0   xn 

Fórmula recurrente iterativa:

X (k + 1) =D−1B − D−1 (L + U)X (k) , k =0,1,2,... (iteraciones)
Fórmula iterativa con las matrices desarrolladas
 x1(k + 1)   b1 / a1,1   0 a1,2 / a1,1 a1,3 / a1,1 ... a1,n / a1,1   x1 
(k)
 (k + 1)  b / a   a / a  
 x2   2 2,2  −  2,1 2,2 0 a2,3 / a2,2 ... a2,n / a2,2   x (k)
2 
=
   :   : : : ... :  : 
 :      
 x (k + 1)   n n,n   n,1 an,n
b / a a / an,2 / an,n an,3 / an,n ... 0   xn(k) 
 n 

X(0) es el vector inicial. A partir de este vector se obtienen sucesivamente los vectores X(1), X(2), ...
 97

Si el método converge, entonces la sucesión { X(k) }k=0,1,2,.. tiende al vector solución X

5.1.6 Práctica computacional con la notación matricial

Resuelva el ejemplo anterior usando la ecuación de recurrencia matricial directamente desde la
ventana de commandos de MATLAB

>> a=[5 -3 1; 2 4 -1; 2 -3 8]

a=
5 -3 1
2 4 -1
2 -3 8
>> b=[5;6;4]
b=
5
6
4
>> d=[5 0 0; 0 4 0; 0 0 8] Matriz diagonal
d=
5 0 0
0 4 0
0 0 8
>> l=[0 0 0; 2 0 0; 2 -3 0] Matriz triangular inferior con ceros en la diagonal
l=
0 0 0
2 0 0
2 -3 0
>> u=[0 -3 1; 0 0 -1; 0 0 0] Matriz triangular superior con ceros en la diagonal
u=
0 -3 1
0 0 -1
0 0 0
>> x=[1;1;1] Vector inicial
x=
1
1
1
>> x=inv(d)*b-inv(d)*(l+u)*x Formula iterativa de Jacobi en forma matricial
x=
1.4000
1.2500
0.6250
>> x=inv(d)*b-inv(d)*(l+u)*x
x=
1.6250
0.9562
0.6187
>> x=inv(d)*b-inv(d)*(l+u)*x
x=
1.4500
0.8422
0.4523

Etc
 98

5.2 Método de Gauss-Seidel

Se diferencia del método anterior al usar los valores más recientes del vector X, es decir aquellos
que ya están calculados, en lugar de los valores de la iteración anterior como en el método de
Jacobi. Por este motivo, podemos suponer que el método de Gauss-Seidel en general converge o
diverge más rápidamente que el método de Jacobi.

5.2.1 Formulación matemática

La fórmula de Gauss-Seidel se la obtiene directamente de la fórmula de Jacobi separando la
sumatoria en dos partes: los componentes que aún no han sido calculados se los toma de la
iteración anterior k, mientras que los que ya están calculados, se los toma de la iteración k+1:

i −1 n
1
xi(k + 1) =
ai,i
(bi - ∑ ai,jx(kj +1) - ∑ ai,jx(k)
j ); i = 1, 2, ..., n; k = 0, 1, 2, ...
j= 1 j= i + 1

En general, el método de Gauss-Seidel requiere menos iteraciones que el método de Jacobi en

caso de que converja

Ejemplo. Dadas las ecuaciones:

5x1 – 3x2 + x3 = 5
2x1 + 4x2 – x3 = 6
2x1 – 3x2 + 8x3 = 4

Formule un sistema iterativo con el método de Gauss-Seidel:

x1= 1/5 (5 + 3x2 – x3 )
x2= 1/4 (6 – 2x1+ x3 )
x3= 1/8 (4 – 2x1+ 3x2)

Fórmula iterativa:
x1(k+1) = 1/5 (5 + 3x2(k) – x3(k))
x2(k+1) = 1/4 (6 – 2x1(k+1) + x3(k))
x3(k+1) = 1/8 (4 – 2x1(k+1) + 3x2(k+1))

Comenzando con el vector inicial: x1(0) = x2(0) = x3(0) = 1, realice dos iteraciones:

k=0: x1(1) = 1/5 (5 + 3x2(0) – x3(0)) = 1/5 (5 + 3(1) – (1)) = 1.4

x2(1) = 1/4 (6 – 2x1(1) + x3(0)) = 1/4 (6 – 2(1.4) + (1)) = 1.05
x3(1) = 1/8 (4 – 2x1(1) + 3x2(1)) = 1/8 (4 – 2(1.4) + 3(1.05)) = 0.5438
k=1:
x1(2) = 1/5 (5 + 3x2(1) – x3(1)) = 1/5 (5 + 3(1.05) – (0.5438))=1.5212
x2(2) = 1/4 (6 – 2x1(2) + x3(1)) = 1/4 (6 – 2(1.5212) + (0.5438)) = 0.8753
x3(2) = 1/8 (4 – 2x1(2) + 3x2(2)) =1/8 (4 – 2(1.5212) + 3(0.8753))=0.4479
 99

5.2.2 Manejo computacional de la fórmula de Gauss-Seidel

La siguiente función en MATLAB recibe un vector X y entrega un nuevo vector X calculado en
cada iteración

function x = gaussseidel(a,b,x)
n=length(x);
for i=1:n
s=a(i,1:i-1)*x(1:i-1)+a(i,i+1:n)*x(i+1:n); %Usa el vector x actualizado
x(i) = (b(i) - s)/a(i,i);
end

Resuelva el ejemplo anterior usando la función gaussseidel:

>> a = [5, -3, 1; 2, 4, -1;2, -3, 8];

>> b = [5; 6; 4];
>> x = [1; 1; 1]; Vector inicial
>> x=gaussseidel(a,b,x) Repetir este comando y observar la convergencia
x=
1.4000
1.0500
0.5438
>> x=gaussseidel(a,b,x)
x=
1.5213
0.8753
0.4479
>> x=gaussseidel(a,b,x)
x=
1.4356
0.8942
0.4764
>> x=gaussseidel(a,b,x)
x=
1.4412
0.8985
0.4766 etc.

En general, la convergencia es más rápida que con el método de Jacobi

5.2.3 Instrumentación computacional del método de Gauss-Seidel

Similar al método de Jacobi , conviene instrumentar el método incluyendo el control de iteraciones
dentro de la función, suministrando los parámetros apropiados.

Los datos para la función son la matriz de coeficientes a y el vector de constantes b de un

sistema lineal. Adicionalmente recibe un vector inicial x, el criterio de error e y el máximo de
iteraciones permitidas m. La función entrega el vector x calculado y el número de iteraciones
realizadas k. Si el método no converge, x contendrá un vector nulo y el número de iteraciones k
será igual al máximo m.
 100

function [x,k] = gaussseidelm(a,b,x,e,m)

n=length(x);
for k=1:m
t=x;
for i=1:n
s=a(i,1:i-1)*x(1:i-1)+a(i,i+1:n)*x(i+1:n);
x(i) = (b(i) - s)/a(i,i);
end
if norm((x-t),inf)<e
return
end
end
x=[ ];
k=m;

Ejemplo. Use la función gaussseidelm para encontrar el vector solución del ejemplo anterior con
una precisión de 0.0001 y determine cuántas iteraciones se realizaron si se comienza con el
vector inicial x=[1; 1; 1]

>> a=[5, -3, 1; 2, 4, -1; 2, -3, 8];

>> b=[5; 6; 4];
>> x=[1; 1; 1];
>> [x,k]=gaussseidelm(a,b,x,0.0001,20)
x=
1.4432
0.8973
0.4757
k=
7

5.2.4 Forma matricial del método de Gauss-Seidel

Dado el sistema de ecuaciones lineales
AX = B
La matriz A se re-escribe como la suma de tres matrices:
A=L+D+U
D es una matriz diagonal con elementos iguales a los elementos de la diagonal principal de A
L es una matriz triangular inferior con ceros en la diagonal principal y los otros elementos iguales a
los elementos respectivos de la matriz A
U es una matriz triangular superior con ceros en la diagonal principal y los otros elementos iguales
a los elementos respectivos de la matriz A
Sustituyendo en la ecuación:
(L + D + U)X = B
LX + DX + UX = B
X = D-1B - D-1LX - D-1UX, Siempre que D-1 exista
 101

Ecuación recurrente con las matrices desarrolladas

 x1   b1 / a1,1   0 0 0 ... 0   x1  0 a1,2 / a1,1 a1,3 / a1,1 ... a1,n / a1,1   x1 
 x  b / a   a / a 0 0 ... 0   x 2  0 0 a2,3 / a2,2 ... a2,n / a2,2   x 2 
 2 =  2 2,2  −  2,1 2,2  −
 :   :   : : : ... :  :  : : : ... :  : 
         
 x 2   n n,n   n,1 n,n
b / a a / a a n,2 / a n,n a n,3 / an,n ... 0   xn  0 0 0 ... 0   xn 

El sistema está en la forma recurrente del punto fijo X = G(X) que sugiere su uso iterativo. En el
método de Gauss-Seidel se utilizan los valores recientemente calculados del vector X.

Fórmula iterativa

X (k + 1) =
D−1B − D−1LX (k + 1) − D−1UX (k) , k =
0,1,2,...

X(0) es el vector inicial. A partir de este vector se obtienen sucesivamente los vectores X(1), X(2), ...

Si el método converge, entonces la sucesión { X(k) }k=0,1,2,.. tiende al vector solución X

Fórmula matricial iterativa del método de Gauss-Seidel:

 x1(k +1)   b1 / a1,1   0 0 0

(k + 1)
... 0   x1  0 a1,2 / a1,1 a1,3 / a1,1 ... a1,n / a1,1   x1 
(k)
 (k + 1)    a / a   (k + 1)     (k) 
 x 2   2 2,2   2,1 2,2
b / a 0 0 ... 0   x 2  0 0 a2,3 / a2,2 ... a2,n / a2,2   x 2 
 = −   −
 :   : : : ... :  :  : : : ... :  : 
 :         (k) 
 x (k + 1)  bn / an,n   an,1 / an,n an,2 / an,n an,3 / an,n ... 0   x (k + 1)  0 0 0 ... 0   xn 
 n   n 

5.3 Método de relajación

Es un dispositivo para acelerar la convergencia (o divergencia) de los métodos iterativos

5.3.1 Formulación matemática

Se la obtiene de la fórmula de Gauss-Seidel incluyendo un factor de convergencia
i −1
ω n
xi(k + 1) = xi(k) + (bi - ∑ ai,j x (k
j
+ 1)
- ∑ ai,jx(k)
j ); i = 1, 2, ..., n; k = 0, 1, 2, ...
ai,i j= 1 j= i

0 < ω < 2 es el factor de relajación

Si ω =1 la fórmula se reduce a la fórmula iterativa de Gauss-Seidel

Ejemplo. Dadas las ecuaciones:

5x1 – 3x2 + x3 = 5
2x1 + 4x2 – x3 = 6
2x1 – 3x2 + 8x3 = 4

Formule un sistema iterativo con el método de Relajación:

x1(k+1) = x1(k) + ω /5 (5 – 5x1(k) + 3x2 (k) – x3 (k))
x2(k+1) = x2(k) + ω /4 (6 – 2x1(k+1) – 4x2(k) – x3(k))
x3 (k+1)= x3(k) + ω /8 (4 – 2x1(k+1) + 3x2(k+1) – 8x3(k))