Métodos Iteractivos para Sistemas de Ecuaciones Lineales2
Métodos Iteractivos para Sistemas de Ecuaciones Lineales2
Métodos Iteractivos para Sistemas de Ecuaciones Lineales2
FACULTAD DE INGENIERÍA
ANÁLISIS NUMÉRICO
TECNICAS ITERATIVAS PARA RESOLVER SISTEMAS LINEALES 2
Docente: Mg. Alvaro Espinosa Pérez
INTRODUCCIÓN
Una técnica iterativa para resolver un sistema lineal 𝐴𝑥 = 𝑏 de 𝑛 𝑥 𝑛 empieza con una aproximación inicial
𝑥 (0) a la solución 𝑥, y genera una sucesión de vectores {𝑥 𝑘 }∞𝑘=0 que converge a 𝑥. La mayoría de estas
técnicas iterativas involucran un proceso que convierte el sistema 𝐴𝑥 = 𝑏 en un sistema equivalente de la
forma 𝑥 = 𝑇 𝑥 + 𝑐 para alguna matriz 𝑇 de 𝑛 𝑥 𝑛 y un vector 𝑐. Ya seleccionado el vector inicial 𝑥 (0) la
sucesión de vectores de solución aproximada se genera calculando:
Las técnicas iterativas se emplean raras veces para resolver sistemas lineales de dimensió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.
MÉTODO DE JACOBI
El método de Jacobi es un método iterativo con el cual se resuelve el sistema lineal 𝐴𝑥 = 𝑏. Se despeja la
primera variable de la primera ecuación, la segunda variable de la segunda ecuación y así sucesivamente. Si
tenemos la k-ésima aproximación 𝑥 (𝑘) entonces substituimos en los términos de la derecha y obtenemos el
valor de 𝑥 𝑘+1 .
Como valor inicial se suele tomar 𝑥 (0) = (0, 0, 0) aunque puede elegirse cualquier otro punto. Miremos un
ejemplo que nos permita comprender el método de Jacobi, el cual es usado para encontrar la solución de un
sistema cuadrado de ecuaciones lineales.
Este consiste en resolver la i-ésima ecuación de 𝐴 𝑥 = 𝑏 para 𝑥𝑖 para obtener, siempre y cuando 𝑎𝑖𝑖 ≠ 0, que:
𝑛
𝑎𝑖𝑗 𝑥𝑗 𝑏𝑖
𝑥𝑖 = ∑ (− )+
𝑎𝑖𝑖 𝑎𝑖𝑖
𝑗=1
𝑗≠𝑖
(𝑘) (𝑘+1)
Y generar cada 𝑥𝑖 de las componentes de 𝑥𝑖 para 𝑘 > 1 con:
𝑛
(𝑘) 1 (𝑘−1)
𝑥𝑖 = ∑(−𝑎𝑖𝑘 𝑥𝑗 ) + 𝑏𝑖 𝑝𝑎𝑟𝑎 𝑖 = 1,2, … , 𝑛
𝑎𝑖𝑖
𝑗=1
[ 𝑗≠𝑖 ]
El método también puede escribirse en la forma 𝑥 (𝑘) = 𝑇 𝑥 (𝑘−1) + 𝑐 dividiendo a 𝐴 en su parte diagonal
y no-diagonal. Para ver esto, sean 𝐷 la matriz diagonal cuya diagonal es la misma que la diagonal de 𝐴, −𝐿 la
parte triangular estrictamente inferior de 𝐴, y −𝑈 la parte triangular estrictamente superior de A.
Con esta notación, se separa en:
𝑥 = 𝐷 −1 (𝐿 + 𝑈)𝑥 + 𝐷 −1 𝑏
Donde:
1
𝐷 −1 = [ ]
𝑎𝑖𝑖
𝑥 (𝑘) = 𝑇𝑗 𝑥 (𝑘−1) + 𝑐𝑗
10 + 𝑥2 − 𝑥3
𝑥1 =
5
11 − 2 𝑥1 + 𝑥3
𝑥2 =
8
3 + 𝑥1 − 𝑥2
𝑥3 =
4
(0) (0)
(1) 10 + 𝑥2 − 𝑥3 10 + 0 − 0
𝑥1 = = =2
5 5
(0) (0)
(1) 11 − 2𝑥1 + 𝑥3 11 − 2(0) + 0
𝑥2 = = = 1,375
8 8
(0) (0)
(1) 3 + 𝑥1 − 𝑥2 3 +0−0
𝑥3 = = = 0,75
4 4
(2) (1) (1) (1)
Para hallar la segunda iteración 𝑥𝑖 con 𝑥 (1) = (𝑥1 , 𝑥2 , 𝑥3 ) = (2, 1,375, 0,75):
(1) (1)
(2) 10 + 𝑥2 − 𝑥3 10 + 1,375 − 0,75
𝑥1 = = = 2,125
5 5
(1) (1)
(2) 11 − 2𝑥1 + 𝑥3 11 − 2(2) + 0,75
𝑥2 = = = 0,96875
8 8
(1) (1)
(2) 3 + 𝑥1 − 𝑥2 3 + 2 − 1,375
𝑥3 = = = 1,0625
4 4
Y continuamos el proceso iterativo.
5 −1 1 𝑥1 10
(2 8 −1 ) (𝑥 2 ) = (11)
−1 1 4 𝑥3 3
Tenemos:
5 −1 1 10
𝐴=( 2 8 −1) , 𝑏 = (11)
−1 1 4 3
5 0 0 0 0 0 0 1 −1
𝐷 = (0 8 0) , 𝐿 = (−2 0 0) 𝑦 𝑈 = (0 0 1)
0 0 4 1 −1 0 0 0 0
Como:
𝑥 (𝑘) = 𝐷 −1 (𝐿 + 𝑈)𝑥 (𝑘−1) + 𝐷 −1 𝑏 = 𝑇𝑗 𝑥 (𝑘−1) + 𝑐𝑗 , 𝑘 = 1,2, …
0 0,2 −0,2 0 2 2
𝑥 (𝑘) = 𝑇𝑗 𝑥 (𝑘−1) + 𝑐𝑗 ⟹ 𝑥 (1) = 𝑇𝑗 𝑥 (0) + 𝑐𝑗 = (−0,25 0 0,125) (0) + (1,375) = (1,375)
0,25 −0,25 0 0 0,75 0,75
k (𝒌)
𝒙𝟏
(𝒌)
𝒙𝟐
(𝒌)
𝒙𝟑
0 0 0 0
1 2 1,375 0,75
2 2,125 0,96875 1,0625
3 1,98125 0,9765625 1,015625
4 1,9921875 1,00664063 0,99140625
5 2,00304688 1,00087891 1,00019531
6 2,00013672 0,9992627 1,00071289
7 1,99970996 1,00005493 0,99985596
8 2,00003979 1,0000545 0,9999635
9 2,0000182 0,99998549 1,00001907
10 1,99999328 0,99999783 0,99999978
11 1,99999961 1,00000165 0,99999838
12 2,00000066 0,99999989 1,00000031
13 1,99999992 0,99999987 1,00000009
14 1,99999996 1,00000003 0,99999996
15 2,00000001 1,00000001 1
16 2 1 1
Como en todo método iterativo, debe establecerse un criterio de paro para las iteraciones, el cual, en este
contexto, puede ser:
(𝑘) (𝑘−1)
𝑥𝑖 − 𝑥𝑖
| (𝑘)
| <∈
𝑥𝑖
Uno de los principales problemas de los métodos iterativos es la garantía de que el método va a converger, es
decir, va a producir una sucesión de aproximaciones cada vez efectivamente más próximas a la solución. En el
caso del método de Jacobi no existe una condición exacta para la convergencia. Lo mejor es una condición que
garantiza la convergencia, pero en caso de no cumplirse puede o no haberla es la siguiente: Si la matriz de
coeficientes original del sistema de ecuaciones es diagonalmente dominante, el método de Jacobi seguro
converge.
𝑥 (𝑘) = 𝑇𝑔 𝑥 (𝑘−1) + 𝑐𝑔
Para que la matriz triangular inferior (𝐷 − 𝐿) sea no singular, es necesario y suficiente que 𝑎𝑖𝑖 ≠ 0 para
cada 𝑖 = 1,2, … , 𝑛.
10 + 𝑥2 − 𝑥3
𝑥1 =
5
11 − 2 𝑥1 + 𝑥3
𝑥2 =
8
3 + 𝑥1 − 𝑥2
𝑥3 =
4
Tenemos que la forma de realizar las iteraciones es:
(𝑘−1) (𝑘−1)
(𝑘) 10 + 𝑥2 − 𝑥3
𝑥1 =
5
(𝑘) (𝑘−1)
(𝑘) 11 − 2𝑥1 + 𝑥3
𝑥2 =
8
(𝑘) (𝑘)
(𝑘) 3 + 𝑥1 − 𝑥2
𝑥3 =
4
(0) (0) (0)
A partir del punto inicial 𝑥 (0) = (𝑥1 , 𝑥2 , 𝑥3 ) = (0, 0, 0), tenemos:
(0) (0)
(1) 10 + 𝑥2 − 𝑥3 10 + 0 − 0
𝑥1 = = =2
5 5
(1) (0)
(1) 11 − 2𝑥1 + 𝑥3 11 − 2(2) + 0
𝑥2 = = = 0,875
8 8
(1) (1)
(1) 3 + 𝑥1 − 𝑥2 3 + 2 − 0,875
𝑥3 = = = 1,03125
4 4
(1) (1) (1) (1)
Para 𝑥𝑖 = (𝑥1 , 𝑥2 , 𝑥3 ) = (2, 0,875, 1,03125), podemos calcular 𝑥𝑖(2) así:
(1) (1)
(2) 10 + 𝑥2 − 𝑥3 10 + 2 − 1,03125
𝑥1 = = = 2,193
5 5
(2) (1)
(2) 11 − 2𝑥1 + 𝑥3 11 − 2(2,193) + 1,03125
𝑥2 = = = 0,995
8 8
(2) (2)
(2) 3 + 𝑥1 − 𝑥2 3 + 2,193 − 0,995
𝑥3 = = = 1,0495
4 4
5 −1 1 10
𝐴=( 2 8 −1) , 𝑏 = (11)
−1 1 4 3
Con:
5 0 0 0 0 0 0 1 −1
𝐷 = (0 8 0) , 𝐿 = (−2 0 0) 𝑦 𝑈 = (0 0 1)
0 0 4 1 −1 0 0 0 0
5 0 0 −1 0 1 −1 0,2 0 0 0 1 −1
𝑇𝐺 = (𝐷 − 𝐿)−1 𝑈 = ( 2 8 0) (0 0 1 ) = ( −0,05 0,125 0 ) (0 0 1)
−1 1 4 0 0 0 0,0625 −0,03125 0,25 0 0 0
0 0,2 −0,2
𝑇𝐺 = (0 −0,05 0,175 )
0 0,0625 −0,09375
0,2 0 0 10 2
𝑐𝐺 = (𝐷 − 𝐿)−1 𝑏 = ( −0,05 0,125 0 ) (11) = ( 0,875 )
0,0625 −0,03125 0,25 3 1,03125
Entonces, la primera iteración es:
0 0,2 −0,2 0 2 2
𝑥 (1) = (0 −0,05 0,175 ) (0) + ( 0,875 ) = ( 0,875 )
0 0,0625 −0,09375 0 1,03125 1,03125
La segunda iteración:
0 0,2 −0,2 2 2
𝑥 (2) = 𝑇𝐺 𝑥 (1) + 𝑐𝐺 = (0 −0,05 0,175 ) ( 0,875 ) + ( 0,875 )
0 0,0625 −0,09375 1,03125 1,03125
1,96875
𝑥 (2) = (1,0117188)
0,9892578
k (𝒌)
𝒙𝟏
(𝒌)
𝒙𝟐
(𝒌)
𝒙𝟑
0 0 0 0
1 2 0,875 1,03125
2 1,96875 1,01171875 0,98925781
3 2,00449219 0,99753418 1,0017395
4 1,99915894 1,0004277 0,99968281
5 2,00014898 0,99992311 1,00005647
6 1,99997333 1,00001373 0,9999899
7 2,00000477 0,99999755 1,0000018
8 1,99999915 1,00000044 0,99999968
9 2,00000015 0,99999992 1,00000006
10 1,99999997 1,00000001 0,99999999
11 2 1 1
Como podemos notar, el método de Gauss-Seidel nos da la respuesta en 11 iteraciones, en tanto que con
Jacobi se obtuvo en 16, lo cual habla de la rapidez de convergencia del uno comparado con el otro.
Este método es consistente para cualquier 𝑤 ≠ 0. Cuando 𝑤 = 1 coincide con el método de Gauss-Seidel.
En particular cuando 0 < 𝑤 < 1 el método se llama sub-relajación mientras que si 𝑤 > 1 se llama sobre
relajación.
Para matrices arbitrarias 𝐴 ∈ 𝑅 𝑛 𝑥 𝑛 en general, no es posible determinar exactamente un valor óptimo para
el párametro de relajación 𝑤.
Tal como en los métodos anteriores, para encontrar las componentes del vector de aproximación, debemos
utilizar la siguiente igualdad:
𝑖−1 𝑛
(𝑘) 𝑤 (𝑘) (𝑘−1) (𝑘−1)
𝑥𝑖 = [− ∑(𝑎𝑖𝑗 𝑥𝑗 ) − ∑ (𝑎𝑖𝑗 𝑥𝑗 ) + 𝑏𝑖 ] + (1 − 𝑤)𝑥𝑖
𝑎𝑖𝑖
𝑗=1 𝑗=𝑖+1
Si 𝑤 ∈ 𝑅, una condición necesaria para que el método converja es que 0 < 𝑤 < 2. Esta condición es también
suficiente si A es simétrica definida positiva.
A partir del punto inicial (𝑥1(0) , 𝑥2(0), 𝑥3(0) ) = (0, 0, 0), tenemos:
(0) (0)
(1) 10 + 𝑥2 − 𝑥3 10 + 0 − 0
𝑥1 = 1,25 [ ] + (1 − 1,25)𝑥1(0) = 1,25 [ ] + (−0,25)(0) = 2,5
5 5
(1) (0)
(1) 11 − 2𝑥1 + 𝑥3 11 − 2(2,5) + 0
𝑥2 = 1,25 [ ] + (1 − 1,25)𝑥2(0) = 1,25 [ ] + (−0,25)(0) = 0,9375
8 8
(1) (1)
(1) 3 + 𝑥1 − 𝑥2 3 + 2,5 − 0,9375
𝑥3 = 1,25 [ ] + (1 − 1,25)𝑥3(0) = 1,25 [ ] + (−0,25)(0) = 1,4258
4 4
Los resultados del proceso iterativo se muestran en la siguiente tabla:
0 𝟎 𝟎 𝟎
1 2 0,875 1,03125
2 1,96875 1,01171875 0,98925781
3 2,00449219 0,99753418 1,0017395
4 1,99915894 1,0004277 0,99968281
5 2,00014898 0,99992311 1,00005647
6 1,99997333 1,00001373 0,9999899
7 2,00000477 0,99999755 1,0000018
8 1,99999915 1,00000044 0,99999968
9 2,00000015 0,99999992 1,00000006
10 1,99999997 1,00000001 0,99999999
11 2 1 1
Un problema que se presenta al usar el método SOR, es cómo escoger el valor apropiado de 𝑤. Aun cuando
no se conoce una respuesta completa a esta pregunta para un sistema lineal general 𝑛 × 𝑛, los siguientes
resultados pueden usarse en ciertas situaciones:
TEOREMA (Kahan): Si 𝑎𝑖𝑖 ≠ 0 para cada 𝑖 = 1, 2, . . . , 𝑛, entonces 𝜌(𝑇𝑤 ) ≥ |𝑤 − 1|. Esto implica que
𝜌(𝑇𝑤 ) < 1 sólo si 0 < 𝑤 < 2, donde 𝑇𝑤 = (𝐷 − 𝑤 𝐿)−1 [(1 − 𝑤) 𝐷 + 𝑤𝑈] es la matriz de iteración del
método SOR.
TEOREMA (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 𝑥 (0) del vector solución.
TEOREMA: Si A es una matriz positiva definida y tridiagonal, entonces 𝜌(𝑇𝐺 ) = [𝜌(𝑇𝑗 )]2 < 1, la elección
óptima de 𝑤 para el método SOR es:
1
𝑤=
1 + √1 − [𝜌(𝑇𝑗 )]2
y con este valor de 𝑤, 𝜌(𝑇𝑤) = 𝑤 − 1.
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: Si ‖𝑇‖ < 1, para cualquier norma matricial natural, entonces la sucesión definida en la ecuación
𝑥 𝑘 = 𝑇 𝑥 𝑘−1 + 𝑐 , {𝑥 𝑘 }∞𝑘=0 , converge para cualquier 𝑥
(0) 𝜖𝑅 𝑛 , a un vector 𝑥 ∈ 𝑅 𝑛 , y se satisfacen las
Con:
‖𝑥 (1) − 𝑥 (0) ‖ = ‖𝑇𝑐 ‖ ≤ ‖𝑇‖‖𝑐 ‖
Y nos queda:
‖𝑇‖𝑘+1
‖𝑥 − 𝑥 𝑘 ‖ ≤ ‖𝑐 ‖
1 − ‖𝑇‖
Ejemplo: Considere el siguiente sistema:
5𝑥1 − 𝑥2 + 𝑥3 = 10
2𝑥1 + 8 𝑥2 − 𝑥3 = 11
−𝑥1 + 𝑥2 + 4 𝑥3 = 3
Vamos a demostrar que el método de Jacobi es convergente para el sistema; reduciendo el sistema a la forma
especial para la iteración de Jacobi, tenemos:
10 + 𝑥2 − 𝑥3
𝑥1 = = 0 + 0,2 𝑥2 − 0,2 𝑥3 + 2
5
11 − 2 𝑥1 + 𝑥3
𝑥2 = = −0,25 𝑥1 + 0 + 0,125 𝑥3 + 1,375
8
3 + 𝑥1 − 𝑥2
𝑥3 = = 0,25 𝑥1 − 0,25 𝑥2 + 0 + 0,75
4
En este ejemplo,
0 0,2 −0,2 2
𝑇𝑗 = (−0,25 0 0,125) 𝑦 𝑐𝑗 = (1,375)
0,25 −0,25 0 0,75
Utilizando, la norma 𝑙1 , tenemos:
En consecuencia, el proceso de iteración para el sistema dado es convergente. Si ahora consideramos como
aproximación inicial de la raíz 𝑥 el vector:
3
(0)
𝑥 = 𝑐 = (2, 1,375, 0.75) ⇒ ‖𝑐 ‖1 = ∑|𝑥𝑖 | = 2 + 1,375 + 0.75 = 4,125
𝑖=1
Sea ahora 𝑘 el número de iteraciones requeridas (Ej. 𝐸 = 10−4 ) para conseguir la exactitud especificada,
tenemos:
𝑘
‖𝑇‖1 𝑘+1 0,5𝑘+1 𝑥 4,125
‖𝑥 − 𝑥 ‖1 ≤ ‖𝑐 ‖1 = < 10−4
1 − ‖𝑇‖1 0.5
0,5𝑘+1 < 1,2121 𝑥 10−5 ⇒ log10 0,5𝑘+1 < log10 (1,2121 𝑥 10−5 )
(𝑘 + 1)log10 (0,5) < −4,9165
−(𝑘 + 1)0,3010 < −4,9165
𝑘 > 16,33 − 1 ≈ 15,33
Se puede tomar 𝑘 = 16. Nótese que la estimación teórica del número de iteraciones necesarias para asegurar
la exactitud especificada es excesivamente alta. A menudo se obtiene la exactitud deseada en un número
menor de iteraciones.
TEOREMA: Sea 𝐴 una matriz real de orden 𝑛, y sea 𝑏 una matriz columna real del mismo orden. Si la matriz
𝐴 es diagonalmente dominante (por filas o por columnas) en sentido estricto, entonces la iteración de Jacobi
converge a la solución del sistema 𝐴𝑥 = 𝑏 para cualquier punto de partida 𝑥0 ∈ 𝑅 𝑛 .
TEOREMA: Si 𝐴 es una matriz simétrica real definida positiva de orden 𝑛, y 𝑏 es una matriz columna real del
mismo orden, entonces la iteración de Gauss-Seidel converge a la única solución del sistema 𝐴𝑥 = 𝑏 para
cualquier punto de partida 𝑥0 ∈ 𝑅 𝑛 .