Tema 2: Sistemas Lineales

Autor: Enrique Mérida Casermeiro

Un sistema lineal de ecuaciones puede expresarse:
a1,1 x1 + a1, 2 x2 + ... + a1,n xn = b1
a2,1 x1 + a2, 2 x2 + ... + a2,n xn = b2
. . .
am ,1 x1 + am , 2 x2 + ... + am ,n xn = bm
o bien matricialmente mediante: Ax=b con A matriz
mxn, b vector mx1, x vector nx1.
Los objetivos del tema será resolver numéricamente
estos sistemas, cálculo del determinante de una matriz, de
su inversa, autovalores y autovectores.
2
Conocimientos previos: Muchos conceptos relaciona-
dos con el tema se suponen conocidos de Álgebra, entre
ellos:
- Según el número de soluciones los sistemas lineales
pueden ser:
SCD: Sistema compatible determinado (solución
única).
SCI: Sistema compatible indeterminado (infinitas
soluciones).
SI: Sistema incompatible (no tienen solución).
- Rango de una matriz.
- Regla de Cramer.
- Teorema de Rouché-Frobernius.
- Cálculo de la matriz inversa mediante: A−1 = Adj ( A' )
| A3 |
- Matriz inversible = tiene inversa = determinante no cero.
- Matriz singular = no tiene inversa = determinante cero.
- Matriz diagonal.
- Matriz triangular superior.
- Matriz triangular inferior.
- Matriz traspuesta.
- Matriz simétrica.
- Matriz antisimétrica: A=-A’
- Matriz simétrica definida positiva: xAx’>0
- Matriz simétrica semidefinida positiva xAx’ mayor o igual
que cero.
Se caracterizan mediante los sucesivos determinantes que
deben ser mayores que 0 (def. positiva) o mayores o iguales
que cero (semidef. positiva).
4
 - Autovalores de una matriz. Calculo.
- Autovectores de una matriz. Cálculo.
- Matrices equivalentes: tienen los mismos autovalores.

Ineficiencia del método de Cramer:

Si A es una matriz cuadrada y |A| es distinto de cero, la
solución del sistema Ax=b puede hacerse mediante:
∆i
xi =
∆
luego la solución del sistema necesita del cálculo de n+1
determinantes y n divisiones.
Cada determinante requiere: n!-1 sumas y n!(n-1)
productos. El sistema necesita n(n+1)!-1 operaciones. 5
Norma vectorial: Dado un espacio vectorial V se define
una norma vectorial como una aplicación ||.||:V Æ R que
verifica:
1) Para todo vector x se verifica ||x|| ≥ 0
2) ||x|| = 0 si y solo si x = 0
3) Para todo real λ y para todo vector x, se verifica : λ x = λ x
4) Para cualesquiera vectores x e y, se verifica :||x + y|| ≤ ||x|| + ||y||

La norma de un vector es una medida de su tamaño.

Puede interpretarse como la distancia del punto x al origen.
Análogamente la distancia entre dos vectores puede consi-
derarse como la norma del vector diferencia entre ellos:
d(x,y)=||x-y||.

6
Ejemplos de normas vectoriales:
Norma k:
Dado un vector x se define su norma k mediante:
n k

∑x
k k k
x k =k i = k
x1 + x2 + ... + xn
i =1

Norma 1: n
x 1 = ∑ xi = x1 + x2 + ... + xn
i =1

Norma 2 o euclídea:
n

∑x
2 2 2 2
x2= i = x1 + x2 + ... + xn
i =1
7
 Norma infinita:
x ∞
= max i xi

Ejemplo:
Dado el vector x=(2, -3, 4, -5) resulta:
x 1 = 2 + 3 + 4 + 5 = 14
x 2 = 4 + 9 + 16 + 25 = 54
x ∞
=5

y para la norma 4:
x 4 = 4 16 + 81 + 256 + 625 = 4 978

8
 Norma Matricial:
Es una aplicación de las matrices en R, que verifica las 4
propiedades de norma y además:
5) Para cualesquiera matrices A y B se verifica : AB ≤ A B

Normas inducidas:
Para cada norma vectorial existe una norma matricial
inducida definida por:
A N
= max x =1 Ax

Las normas inducidas verifican ||I||=1.

Radio espectral de una matriz A es el módulo del
autovalor con mayor módulo. ρ ( A) = max i λ i 9
 Ejemplos de normas matriciales:
Norma 1: Es la inducida por la norma 1 vectorial y se
calcula mediante:
A 1 = max j ∑ | ai , j |
i

Norma infinito: Es la inducida por la norma infinito

vectorial y se calcula mediante:
A ∞ = max i ∑ | ai , j |
j

Norma 2 o espectral: Es la inducida por la norma 2 vecto-

rial y se calcula mediante:
A 2 = ( ρ ( A' A) )
1/ 2

10
Norma de Frobernius:
Es un ejemplo de norma no inducida por una norma
vectorial ya que:
|| I || F = n (n orden de la matriz)

Se define por:

1/ 2
 
=  ∑∑ ai , j
2
AF 

 i j 

11
Ejemplo: Dadas las matrices A y B hallar sus distintas
normas matriciales.  2 0 - 1
   - 2 3 0
A= 0 1 3 B =  
- 2 3 5   0 − 1 1
 

A 1 = max{4,4,9} = 9, A ∞ = max{3,4,10} = 10,

A F = 4 + 1 + 1 + A9 =+max{4,4,9}
4 + 9=+9, 25A == max{3,4,10
1 53 } = 10,
∞

A F = 4 + 1 + 1 + 9 + 4 + 9 + 25 = 53

8 6 8  8 - λ 6 8
 
A' A =  6 10 18  6 10 - λ 18 = 0 ⇒
 8 18 35  8 18 35 - λ
 
p(λ ) = −λ 3 + 53λ 2 − 286λ + 36 = 0 que tiene las raices :
0.1289 5.9501 y 46.921, por tanto ρ ( A'A) = 46.921 y
A = 46.921 ≈ 6.8499 12
  − 2 3 0
B =   ⇒ B 1 = max{2,4,−1}, B = max{5,2} = 5
 0 −1 1 
∞

B F
= 4 + 9 + 1 + 1 = 15
 4 −B6=B =− 2 − 023 30  0⇒ ⇒B4B=−=maxmax
λ −6 0
  0 0 − 1 − 1  1 
{2{
, 4
2
1 1
, −
4,1}
− =
1},4, B B= =
max {
max
5,
∞
{
2}
5∞
,=2}5= 5

B ' B =  − 6 10B B = =−4 +419++91++⇒ − 6 10 − λ − 1 = 0 ⇒

 
11+=1 =1515

 0 − 1  414- 6−06 0 4-λ 04-−6 λ 0− 6 −0 1 1 − λ

F F

 6 10
B ' BB=' B -=6 −10 - 1 − 1-6⇒ 10
− -λ
6 10
 0  0- 1 −11 1 0
-1 − λ= 0 ⇒ −1 = 0 ⇒
-01 1-λλ −1 1− λ
2  3 3 2 2 
p(λ ) = -λ 3 + 15λ - 17λ = 0 que tiene las raices :
p ( λp(
) =λ )−=λ -λ+ 15+ 15λ
λ − 17 λ ==00 que
- 17λ quetiene
tienelas
lasraices
raices::
0, 1.235
0, 1.235y 13.765,
y 13.765,
por tanto ρ(B'B )ρ(B'
por tanto = 13B) = 13.765
.765 y y
0, 1.235 y 13.765, por
B B= =13.13
765 ≈ 3≈
.765
2
3tanto
.7101
2
.7101 ρ(B' B) = 13.765 y
B 2 = 13.765 ≈ 3.7101

13
Propiedades de las normas matriciales:

1) Toda norma matricial verifica: ρ ( A) ≤ A

2) El radio espectral de A es la inferior de sus normas:

∀ε > 0, existe una norma de A tal que : A < ρ ( A) + ε

3) La sucesión de potencias de una matriz A converge

a la matriz 0 si y solo si el radio espectral de A es
menor que 1.
{ A } → 0 si y solo si ρ (A) < 1
n

14
Definiciones:
Limite de una sucesión matricial:
{ An } → A ⇔ lim n →∞ An − A = 0

Matriz de diagonal dominante por filas:

∀i, ai,i > ∑ ai , j

j ≠i

Matriz de diagonal dominante por columnas:

∀i, ai,i > ∑ a j ,i

j ≠i

15
Ejemplos: 3 0 0 −1
 
5 3 1   −1 3 1 0 
  B=
A = 1 − 4 − 2 0 2 −4 1 
4 2   
7 2 0 − 2 − 5 
  

La matriz A es estrictamente dominante por filas:

|5|>|3|+|1| |-4|>|1|+|-2| |7|>|4|+|2|
pero no por columnas al fallar las columnas 1 y 2:
|5|=|1|+|4| |-4|<|3|+|2|

La matriz B es también estrictamente dominante por

filas pero no por columnas:
filas: 3>0+0+1, 3>1+1+0, 4>0+2+0 y 5>2+0+2
columnas: 3=1+0+2 (falla la primera columna)
16
 Métodos directos de resolución de sistemas:
Se denominan métodos directos a los que proporcionan la
solución exacta en un número finito de pasos, como alter-
nativa a los métodos iterativos que solo pretenden en cada
paso aproximarse a la solución y en general necesitan infi-
nitos pasos para obtener la solución exacta.
El problema con los métodos directos es que al trabajar
con números inexactos (punto flotante) propagan y aumen-
tan los errores llegando a invalidar la solución obtenida en
matrices de tamaño mediano.
Existen matrices muy mal condicionadas que impiden la
solución por estos métodos para n=2 ó 3, pero en general y
con técnicas especiales pueden resolver hasta sistemas con
40 ó 45 ecuaciones e incógnitas. 17
Casos especiales:
Existen matrices A tales que el sistema Ax=b puede ser
resuelto fácilmente, son interesantes porque la mayoría de
los métodos directos modifican la matriz A para alcanzar
uno de estos.
b1
Matriz diagonal: x1 =
a1,1 x1 = b1 a1,1
a 2 , 2 x2 = b2 b2
x =
⇒ 2 a2 , 2
...
...
a n , n xn = bn bn
xn =
an , n

obteniendose la solución del sistema con n operaciones.

El determinante de A: A= a ∏
i
i ,i
18
Matriz triangular superior:
 bn 
 x n = 
an , n
 
b
 x = n −1 n −1,n n  − a x
a1,1 x1 + a1, 2 x2 + 1,3 3
a x ... + 1,n n = 1
a x b  n −1
an −1,n −1 
 = b  
a 2 , 2 x2 + a2,3 x3 ... + a2,n xn = b1
x
 a 
n ... n


 n ,n

 x = b − a x  
∑
n
+ = −
n −1 n −1, n n

a1,1 x1 a x ... a x
+ a3, 3x 3 + a x ... + a3, nx n = b b 3  ⇒  n −1
a bi  ai , k x k 
1, 2 2 1, 3 3 1, n n 1  x
 ... i =  k =i +1
n −1, n −1


a x
2, 2 2 + a x ... + a x = b
2,3 3
... 2,n n
... 
1 
 a
b −∑ i , i

n
a x
3, 3 3... + a x = b ⇒ 
3, n n 3 ia x  i ,k k
x = ...
k =i +1

an ,n...xn ... = bn  i

 a 
 i ,i 
a x
n,n n =b n 

 ... 


   
1 ∑k = 2 a1, k xk 
n
b − ∑ a xb −
n

 x =  a
 1 k =2 1, k k

 x1 = 
1

a1,1
1,1

 
obteniendose la solución del sistema con n2 operaciones.
El determinante de A:
A = ∏ ai ,i
i
19
Matriz triangular inferior:
 b1 
 x1 = 
a1,1 x1 = b1 a
 1,1

a2,1 x1 + a2, 2 x2 = b2  x = b2 − a2,1 x1 
 2 a2 , 2 
a3,1 x1 + a3, 2 x2 + a3,3 x3 = b3 ⇒x = b  
...
1
 
 
1
a x =b a
  1,1

a x...
∑
1,1 1 1
 x = b − a x  i −1
 b − ai ,k xk 
+a x =b
a
2 2 ,1 1
2 ,1 1 2, 2 2
 2 2
i
an ,1 x1 + an , 2 x2 + an ,3 x3 + ......+ an ,n xn = bn b − ∑...  axix =
a x +a x +a x =b ⇒   k
2, 2
=1

3,1 1 3, 2 2 3, 3 3 3

a
i −1

a 
i i ,k k
x =
k =1
a x + a x + a x + ... + a x = b
n ,1 1 n,2 2 n ,3 3 n,n n n i  2, 2
...  
  2, 2

  ...

obteniendose la solución del sistema con n 2 operaciones.

El determinante de A:
A = ∏ ai ,i
i
20
 Ejemplos:
1) Dado el sistema: 4x =4 Sustitución
2x-3y =8 progresiva
3x+2y-z =2
tendrá la solución: x=4/4=1 que produce el nuevo sistema:
-3y =8-2x= 6
2y-z=2-3x=-1 del que obtenemos y=6/(-3)=-2, y por
último –z=-1-2y=3, luego z=-3.
-----
2) Dado el sistema: 4x-2y-2z = 0 Sustitución
6y+ z = -1 regresiva
3z = 15
tendrá la solución: z=15/3=5 que produce el nuevo sistema:
4x-2y=0+2z = 10
6y=-1- z = -6 del que obtenemos y=-6/6 =-1, y por
último 4x=10+2y=8, luego x=-8/4=-2. 21
 Método de reducción de Gauss:
Dado el sistema Ax=b, con A no singular, se trata de
encontrar un sistema equivalente Rx=b’ con R triangular
superior.
Teóricamente consiste en encontrar una matriz S tal que
SA=R por lo que el sistema queda:
Ax=b Æ SAx=Sb Æ Rx=b’ (con b’=Sb).
En la práctica consiste en realizar operaciones por filas tal
que la matriz del sistema se haga triangular superior.
Se necesitan dos tipos de transformaciones:
- Las que hacen cero debajo de un término diagonal.
- Las que intercambian filas.
NOTA: El intercambio de filas puede ser controlado
mediante punteros y no necesita ser físico. 22
 Condiciones para evitar el intercambio de filas:
La condición necesaria y suficiente para que no haya que
intercambiar filas es que los sucesivos menores principales
excepto el último sean distintos de cero.
Son condiciones suficientes:
1) A sea definida positiva.
2) A sea de diagonal estrictamente dominante por filas.
3) A sea de diagonal estrictamente dominante por
columnas.
Cálculo del determinante:
El determinante de A podrá calcularse como el producto
de los término de la matriz triangular superior con signo
cambiado si el número de intercambios de fila es impar:
A = (−1) t ∏ ai ,i
23
i
 Método de Gauss (sin pivotaje):
Resolver el sistema: 3x+2y+z=5
2x+2y -z=7
4x+ y -z=4
Formamos la matriz ampliada y triangularizamos:

 3 2 1 5  f ← f − 2 f  3 2 1 5 
   2 2 1  
 2 2 −1 7  →  3
4  →  0 2 / 3 − 5 / 3 11 / 3 
 4 1 − 1 4   f 3 ← f 3 − f1   0 − 5 / 3 − 7 / 3 − 8 / 3 
   3   
 
3 2 1 5   z = −1 
 5     11 / 3 + (5 / 3)(−1) 2 
→  f 3 ← f 3 + f 2  →  0 2 / 3 − 5 / 3 11 / 3  ⇒  y = = = 3
 2     2/3 2/3 
 0 0 − 13 / 2 13 / 2   z= 5 − 2 (3 ) − ( −1) 0
= =0 
 3 3 

24
 Dado un sistema de ecuaciones lineales, para hacer cero
debajo de la diagonal k, se toma como pivote el término ak,k
si es diferente de cero en otro caso se intercambia la fila k
con otra j>k cuyo |aj,i| sea máximo.
A continuación se hace cero debajo de la diagonal en la
columna k, mediante:
k
k +1 k ai ,k k  i = k + 1,.., n 
ai , j = ai , j − ak , j  
 j = k + 1,.., n, n + 1
k
ak , k
Número de operaciones: Viene dado por:
4n 3 + 9n 2 − 7 n 2 3
Nop = ≈ n
6 3
NOTA: La fila que está en su lugar no debe ser multi-
plicada por ningún factor y solo se suma o resta una
combinación lineal de la fila pivote. 25
 Gauss con pivotaje:
Si en el paso k-ésimo ak,k=0 debemos intercambiar filas,
pero si es próximo a cero, cualquier error en los coeficien-
tes se incrementa innecesariamente para evitarlo se emplea
el pivotaje.
Pivotaje parcial:
Al elegir el pivote en la columna k-ésima tomamos el
término de mayor valor absoluto de los que están en dicha
columna en o por debajo de la diagonal. En caso de que esté
por debajo se intercambiarán filas.
Pivotaje total:
Al elegir el pivote en la columna k-ésima tomamos el
término de mayor valor absoluto de la submatriz formada
por las filas y columnas k, k+1,...n. Usualmente será
necesario intercambiar filas y columnas. 26
NOTA: En el pivotaje total se deberá controlar especial-
mente los intercambios de columna para asignar el valor
correcto a cada variable al efectuar la sustitución regresiva.

Ejemplo: Resolver el sistema: 3x+2y+z=5; 2x+2y -z=7;

4x+ y -z=4. a) Por pivotaje parcial. b) Por pivotaje
total.
a) El primer pivote será el 4 de la posición (3,1) por lo que
se intercambian las filas primera y tercera y hacemos
cero debajo del término (1,1):
 3 2 1 5  4 1 − 1 4   f ← f − 2 f   4 1 − 1 4 
     2 2 1
 2 2 − 1 7  → { f1 ↔ f 3 } →  2 2 − 1 7  →  4
 → 0 3
2
−1 5 
2 
 4 1 −1 4  3 2 1 5   f 3 ← f 3 − 3 f1   5 7 2
     4   0 4 4 
El nuevo pivote es 3/2 que se encuentra en la posición (2,2)
y por tanto no se intercambian filas. 27
 4 1 − 1 4 

 5/ 4 5  
→  f3 ← f3 − f 2 = f3 − f 2  → 0 3 −1 5 
 3/ 2 6   2 2 
0 0 13 − 13
 6 6 

que mediante sustitución regresiva produce la solución:

z=-1; y=(5+(1/2)(-1))/(3/2)=3; x=(4+z-y)/4=0.
b) El primer pivote será el máximo de la matriz A (sin
contar la columna cuarta) y coincide con el 4 de antes en la
posición (3,1), así:
 3 2 1 5  4 1 − 1 4   f ← f − 2 f   4 1 − 1 4 
     2 2 1
 2 2 − 1 7  → { f1 ↔ f 3 } →  2 2 − 1 7  →  4  → 0 3
2
−1 5 
2 
 4 1 −1 4  3 2 1 5   f 3 ← f 3 − 3 f1   5 7 2
     4   0 4 4 

pero ahora el pivote es 7/4 de la posición (3,3): 28

 4 1 − 1 4   4 −1 1 4 
  
c ↔ c 
→ { f 2 ↔ f3} →  0 5 7 2 →  2 3
 → 0 7 5 2 →
 4 4   ( x, z , y )   4 4 
0 3 −1 5   0 −1 3 5
 2 2   2 2 

Hay que considerar el intercambio en el orden de las

variables al permutar columnas.
 4 −1 1 4 
 
 1/ 2  
→  f3 ← f3 + f2  → 0 7 5 2 
 7/4   4 4 
 0 0 13 39 
 7 7 
que mediante sustitución regresiva obtiene la solución:
y=3; z=(2-(5/4)3)/(7/4)= -1; x=(4-3+(-1))/4=0
Aviso: La fila i que va a ser sustituida no puede
multiplicarse por ningún factor: f i ← f i + kf j 29
 Normalización (scaling):
Cuando la matriz contiene números de magnitudes muy
diferentes conviene normalizarla, es decir, operar con una
matriz equivalente con términos de similar magnitud.
Ejemplo:
El sistema: 0.000023x-0.000007y=0.000034
1230000x+4230000y= 9210000
queda mejor condicionado de la forma:
23x - 7y = 34
123x + 423y = 921
Usualmente cada ecuación se divide por la componente
de A mayor en valor absoluto:
x - (7/23)y = 34/23
123/423x + y = 921/423
30
Comentarios:

1) El objetivo de las técnicas de pivotaje y de normalización

es evitar la propagación innecesaria de errores.
2) Normalmente se emplea el pivotaje parcial, aunque el
pivotaje total con normalización es el que más reduce el
error. Esto aumenta la complejidad del algoritmo.
3) Para el cálculo del determinante de A, debemos llevar la
cuenta de los intercambios y considerar los factores de
normalización.

31
 Método de Gauss-Jordan:
Trata de reducir el sistema en diagonal mediante opera-
ciones por filas:
k
k +1 k ai ,k k  i = 1,2 ,...,k-1,k + 1,...,n 
ai , j = ai , j − ak , j  
 j = k + 1, k + 2,..., n, n + 1
k
ak , k

posteriormente basta resolver el sistema diagonal.

El número de operaciones es n 3+n 2 -n que es mayor que el
de Gauss, por lo que acumula más errores.
Sin embargo:
- Es más fácil de programar que Gauss (al evitar sustitución
regresiva)
- Permite paralelismo ya que con n-1 procesadores cada uno
puede hacer 0 en una fila.
- Es el método usual para calcular la inversa de una matriz.
32
Ejemplo: Calcular la inversa de A y resolver conjunta-
mente los sistemas por Gauss-Jordan:
2x –y +4z= -13 2x – y+4z= 14
x +y + z= 4 x + y+ z= 5
3x +2y- z= 19 3x +2y- z= 7
Para ello montamos la matriz y operamos por filas:

 2 − 1 4 − 1314 1 0 0   f ← f − 1 f   2 − 1 4 − 13 14 1 
0 0
   2 1
→ 0 3 1 0
2
 1 1 1 4 5 0 1 0 →  2  − 1 21 − 2 − 1
   f ← f − 3 f   2 2 2 
 3 2 − 1 19 7 0 0 1   3 3
2  0
1  7 − 7 77 − 14 − 3 0 1
2 2 2 

 2   2 0 10 − 6 38 2 2 0  
 10 
f
 1 ← f + f 2
3 3 3 3  f ← f + f3 
→
1
3 
→ 0 3 − 1 21 − 2 − 1 1 
0 → 
1 1
14 
7  2 2 2 3
 f3 ← f3 − f 2   − 14 14 − 28 −1 −7   f ← f − f3 
0 0 1 2 2
 3   3 3 3 3   14 
33
  
 3 −1 5   f ← f /2   3 −1 5 
2 0 0 4 6 7 7   1 1  1 0 0 2 3 14 2 14 
0   2  
3 0 15 0 −3 3 −3 →  f ← f  → 0 1 0 5 0 −2 1 −1 
 2 2 7 2 14  2
3  
2 7 7
 0 − 28 
0 − 14 14 3
−1 −7 1   − 3  0 0 1 − 32 1 1 −3 
 3 3 3   f3 ← f3   14 2 14 
 14 

Por tanto la solución del primer sistema es:

x=2, y=5, z=-3
la del segundo es:
x=3, y=0, z=2
y la matriz inversa es:
3 −1 5 
 14 2 14 
A−1 =  − 2 1 −1 
 7 7
1 1 −3 
 14 2 14 
34
Métodos de descomposición LU (Crout y Doolittle):
El método consiste en descomponer la matriz A en el
producto de dos matrices triangulares L (triangular
inferior) y U (triangular superior).
Esto no siempre es posible, pero:
Dada una matriz cuadrada A, siempre es posible
encontrar una matriz de permutación P tal que
PA=LU
Matriz de permutación:
Una matriz cuadrada es de permutación si y solo si en
cada fila y columna existe un único término con valor 1
siendo los restantes 0 .
Se llama así pues PA permuta las filas de A y AP
permuta las columnas. 35
Propiedades:
1) Si A es no singular entonces L y U son también no
singulares.
2) Si A es no singular, entonces P será la matriz identidad si
y solo si los sucesivos menores principales son distintos
de cero.
1) La descomposición no es única así cuando:
- Los términos diagonales de U valen 1 se llama
método de Crout.
- Los términos diagonales de L valen 1 se llama
método de Doolittle.
Resolución numérica: Una vez descompuesta
PA=LU el sistema Ax=b se transforma en: PAx=Pb=b’
LUx=b’ y llamando y=Ux se resuelve Ly=b’ y Ux=y lo
que resulta fácil por sustitución progresiva y regresiva. 36
 Ejemplo:
Resolver el sistema Ax=b por el método de Crout.
 4 2 1 7
   
A =  1 −1 1 b = - 2
 3 2 5 2
   

Planteamos la descomposición sin considerar la matriz de

permutación (observar el orden de cálculo):
 a 0 0  1 g h   4 2 1  a = 4 ag = 2 ⇒ g = 1 
        2 
A ⇒  b c 0  0 1 i  =  1 − 1 1  ⇒  b = 1 ⇒ 
1
 d e f  0 0 1   3 2 5  d = 3  ah = 1 ⇒ h = 
        4
 1 - 3  1 3i −1 
 bg + c = -1 ⇒ + c = −1 ⇒ c =   bh + ci = 1 ⇒ − = 1 ⇒ i = 
⇒ 2 2 ⇒ 4 2 2
3 1   3 1 9
 dg + e = 2 ⇒ + e = 2 ⇒ e =  dh + ei + f = 5 ⇒ − + f = 5 ⇒ f = 
 2 2   4 4 37 2 
  
 7 
 y1 = 
4 4
 0 0   y1   7   7

     −2− 5 

Resolvemos Ly = b 4 : 0 1 0  − 3  1 10  1 y = - 2 ⇒  y2 = 4 =
 2  2  4 2  − 3 
L = 1 − 3 0  U =  0 1 −1   2  2 
2 3  1  9   2y3  2
 
( ) ( )( )

3 1  9  2  0 02  1   7 − 1 5 
 2 2
y = 2 − 3

4 2 2 = −1
  3 7 9 
4 0  y


 y =
4


1
2 
0
   1    
7
−2−7 
4=5
Resolvemos Ly = b : 1 − 3
 2     

0   y2  =  - 2  ⇒  y2 =
−3 2
 
1  x  7   ( ) ( )( ) z = −1 
3 1 9   y3   2  2 
1 1
( )
 2 2  7 1 5 
 2 4     4    y = 2 − 3 4 − 2 2 = −1
5 + 1 (−1) 
  3 9 

Re solvemos Ux = y :  0 1 − 1   y  =  5  ⇒  y = 2  2 2 
=2 
 2  2  1


Re solvemos Ux = y :  − 1



 0 0 1 1 21 14 zx   74 - 1 
     
y = 5


⇒ y= x = =2 

 ( )
5 2 + (12 )(−41)− 2 2 − (− 1) 4
z = −
7 1
1 1 ( )
= 1



0 1
2   2  1 
0 0

1   z   - 1  
   

7 − 2(1 )− (− 1)(1 ) 
 1 
4 2 4 =1
x = 
 1 

La solución será x=1, y=2, z=-1.

38
 Determinante: Podrá calcularse mediante:
Si A = LU entonces : A = L U = ∏ li ,i
i

en caso de que P sea distinta de matriz unidad deberá con-

siderarse el signo de |P|.
Comentarios:
El método de Crout es similar al de Gauss y necesita de las
mismas operaciones. De hecho la diferencia es el calculo de
la matriz L tal que inv(L)A=U.
La matriz P es la unidad en las mismas condiciones que el
método de Gauss evita los intercambios de filas.
Este método permite empaquetamiento al poder guardarse
L y U en una matriz como (L\U).
Permite algo de paralelismo en el cálculo de coeficientes.
39
 Ventajas sobre el método de Gauss:

Aunque tiene la misma complejidad que Gauss, la des-

3
composición tiene complejidad del orden de 2n /3, y las
sustituciones progresivas y regresivas posteriores de 2n 2.

Esto permite resolver fácilmente nuevos sistemas con la

misma matriz A (aunque diferente b), pues la parte de
mayor complejidad computacional ya está calculada.

Esto resultará interesante para calcular la matriz inversa

mediante vectores columnas (por ejemplo).

40
 Método de Cholesky o de la raíz cuadrada:

Es una variante del método LU para matrices simétricas

definidas positivas. Así cuando A es simétrica definida
positiva puede ponerse como A=TT’ con T triangular
inferior (luego T’ es triangular superior).
3
El número de operaciones resulta Nop=n /6+n que es
menor que Gauss, LU y Gauss-Jordan por lo que controla
mejor los errores.

Resolución:
Dado el sistema Ax=b, se transforma en TT’x=b, y
llamando T’x=y, resultan dos sistemas triangulares Ty=b y
T’x=y que se resolverán mediante sustitución progresiva y
regresiva. 41
 Ejemplo:
Comprobar que la matriz A es simétrica definida positiva y
resolver el sistema Ax=b, siendo:  4 2 − 1  11 
   
A= 2 5 1 b =  11 
 −1 1 3   − 2
Descomponemos A:   

 2 
 a 0 0  a b c   4 2 − 1  a = 4 → a = 2 
     
 b d 0  0 d e  =  2 5 1  ⇒  ba = 2 → b = 1 
 c e f  0 0 f   − 1 1 3  ca = −1 → c = − 1
     
 2 
 
 b 2
+ d 2
= 5 → d 2
= 4 → d = 2 
 −1 
 3 
→ bc + de = 1 → + 2e = 1 → e = 
 2 4 
c 2 + e 2 + f 2 = 3 → 1 + 9 + f 2 = 3 → f 2 = 35 → f = 35 
 4 16 16 4  42
 La matriz A es simétrica definida positiva por ser simétrica
y de diagonal estrictamente positiva y términos positivos.
Hemos obtenido la descomposición:
 −1 
  2 1 
 2 0 0  2   4 2 − 1
  3   
 1 2 0  0 2 = 2 5 1 
 4 
 −1 3 35    −1 1 3 
  35  
 2 4 4  0 0 
 4 
que se resuelve mediante dos sustituciones:
 11 
   y1 = 
 2 0 
0  y1   11   2 
      11 11 
progresiva :  1 2 0  y 2  =  11  ⇒  + 2y 2 = 11 → y 2 = 
35  y   − 2   2 4 
 −1 3
  3     − 11 33 35 − 21 
 2 4 4  + + y 3 = −2 → y 3 =
 4 16 4 4 35 
43
  −1   11   35 z = − 21 → z = − 3 = −0.6 
2 1     
 2  x   2   4 4 35 5 
3     11  ⇒  2y + 3  − 3  = 11 → y = 8 = 1.6 
regresiva :  0 2 y
  =  
 4    4   4 5  4 5
  z   
35   − 21 8 1 − 3  11
 2x + −  
0 0  9
  = → x = = 1.8
 4   4 35   5 2 5  2 5 


Cálculo del determinante:

El determinante de la matriz A puede calcularse como el
producto de los términos diagonales de la matriz T eleva-
dos al cuadrado:
| A |= ∏ ti ,i
2

Así en el ejemplo anterior |A|=4*4*35/16=35.

44
Comentarios:

Reduce el cálculo para sucesivas b, como en el LU.

Reduce el espacio de almacenamiento y cálculo respec-

to al método LU.

Podría pensarse que multiplicando el sistema por A’:

A’Ax=A’b
tiene la matriz A’A simétrica definida positiva, pero sin
embargo está mal condicionada.

45
 Métodos de ortogonalización. Sistemas
sobredeterminados.
Matriz ortogonal: Una matriz Q se dice ortogonal
si y solo si QQ’=I. (Esto implica Q’Q=I)
-1
La condición anterior es equivalente a: Q =Q’.
En una matriz ortogonal los vectores fila que la forman
deben ser de norma 1 y ortogonales entre sí.

Método: Primero habrá que descomponer la matriz A

como A=QR, donde Q sea ortogonal y R triangular
superior.
Ax=b Æ QRx=b Æ Q’Qx=Q’b=b’ Æ Rx=b’
que se resuelve por sustitución regresiva. 46
 Comentarios:
Los métodos QR necesitan más operaciones que el de
Gauss.
Tienen a su favor que mantienen el condicionamiento de
la matriz y no lo empeoran como le pasa a los gaussianos.
Por tanto es el mejor método directo para problemas mal
condicionados como los sistemas resultantes en las apro-
ximaciones mínimo-cuadráticas.
Es el método que se emplea para los sistemas sobrede-
terminados.
Sistema sobredeterminado: Un sistema se dice
sobredeterminado si carece de solución (sistema incom-
patible o imposible). Un vector x se llama solución del
sistema sobredeterminado Ax=b si minimiza la norma 2
del vector residuo r=Ax-b. 47
Ejemplos:
1) Las instrucciones:
A=[4 4 5; 2 1 3;5 6 4]; [Q,R]=qr(A)
de MATLAB obtiene las matrices:

 - 0.5963 - 0.1988 - 0.7778   - 6.7082 - 7.1554 - 6.8573 

   
Q =  - 0.2981 - 0.8447 0.4444  R =  0 1.3416 - 1.5404 
 - 0.7454 0.4969 0.4444   0 0 - 0.7778 
  

Para resolver el sistema Ax=b con b=(2 3 5)’,

realizaremos b’=Q’b y mediante sustitución regresiva
Rx=b’, es decir:

bprima=Q’*b, x=sustreg( R, bprima )

48
 2) Buscamos una recta y=mx+b que pase por los
puntos (2,3), (-1,2), (-2,2), (0,2), (3,4).
Por lo que se debe verificar:
3=2m+b; 2=-m+b; 2=-2m+b; 2=b; 4=3m+b que no tiene
solución (sistema sobredeterminado). La mejor solución
es la que se obtiene por mínimos cuadrados y se obtiene
de forma eficiente por el método QR.
 2 1  3  - 0.4714 - 0.3558 0.2170 - 0.1729 - 0.7577   - 4.2426 - 0.4714 
       
 -1 1   2  0.2357 - 0.5083 - 0.7079 - 0.4298 - 0.0126   0 - 2.1858 
A= - 2 1 b= 2  ⇒ Q =  0.4714 - 0.5592 0.6410 - 0.1414 0.1851  R = 0 0 
       
 0 1  2  0 - 0.4575 - 0.1967 0.8663 - 0.0392   0 0 
 3 1   4  - 0.3050 0.0467 - 0.1222 0.6244   
   - 0.7071  0 0 

 - 2.8284 
 
 - 5.3374 
   0.3953 
b' = Q' b = 0.3104 ⇒ x =   La recta es : y = 0.3953x + 2.4419
   2.4419 
 - 0.4174 
 0.4910  49
 
 Métodos iterativos:
Para grandes sistemas los métodos directos son inservibles
a causa de la propagación de los errores de cálculo. Puede
considerarse que para sistemas con n>50 el método a em-
plear debe ser iterativo.
Los iterativos toman como vector inicial en cada paso el
obtenido en el paso anterior y por tanto no propagan errores.
Es fácil resolver un sistema mediante un método iterativo
sin guardar la matriz A en memoria.
En los sistemas grandes que surgen en la práctica tienen la
matriz A esparcida y aunque existen métodos directos espe-
ciales, son resueltos usualmente por métodos iterativos.
Los problemas de los iterativos serán: convergencia, velo-
cidad de la misma, vector inicial, condición de parada.
50
 Generalidades:
Dado un sistema Ax=b, buscamos una sucesión de
vectores x (m) que converja hacia la solución del sistema.
Usualmente el vector siguiente se obtiene como función
del anterior: ( m +1)
x = F (x
(m)
)
y para los sistemas lineales resulta:
( m +1)
x = Bx (m)
+C
donde B es la matriz del método y C es un vector columna.
Condición de convergencia:
La condición necesaria y suficiente para que un
método iterativo converja es que el radio espectral
de la matriz B sea menor que 1. ρ (B) < 1 51
 Los métodos iterativos usuales están preparados para que
el punto fijo sea la solución del sistema, es decir:
x=Bx+C cuando x verifica Ax=b
despejando podemos obtener: C=(I-B)x=(I-B)A-1 b, aunque
no es necesario comprobarlo.
NOTA: La convergencia no depende del vector b.
Métodos de partición regular:
Se basan en descomponer la matriz A como A=M-N de
forma que el sistema queda:
-1 -1
Ax=b Æ (M-N)x=b Æ Mx=Nx+b Æ x=M Nx+M b
que resultará útil si M es inversible y el radio espectral de
( )
-1
la matriz del método M N es menor que 1. ρ M −1 N < 1
El algoritmo iterativo resultante es:
( m +1) −1 −1
x = M Nx (m)
+M b 52
 Método de Jacobi:
Dado un sistema de a1,1 x1 + a1, 2 x2 + a1,3 x3 + ... + a1,n xn = b1 
ecuaciones lineales: a2,1 x1 + a2, 2 x2 + a2,3 x3 + ... + a2,n xn = b2 
 
a3,1 x1 + a3, 2 x2 + a3,3 x3 + ... + a3,n xn = b3 
... 

an ,1 x1 + an , 2 x2 + an ,3 x3 + ... + an ,n xn = bn 

x1 =
1
(b1 − a1,2 x2 − a1,3 x3 − ... − a1,n xn ) 
despejamos en la fila i a1,1

1
la incógnita i-ésima: x2 = (b2 − a2,1 x1 − a2,3 x3 − ... − a2,n xn ) 
a2 , 2

...
1 
xi = (bi − ai,1 x1 − ai,2 x2 − ... − ai,n xn ) 
ai ,i 
... 
1 
xn = (bn − an,2 x2 − an,3 x3 − ... − an,n−1 xn−153)
an , n 
 que conduce al método iterativo de Jacobi:
x1
( m +1)
=
1
(
a1,1
(m) (m)
b1 − a1, 2 x2 − a1,3 x3 − ... − a1,n xn
(m)
) 


x2
( m +1)
=
1
(
a2 , 2
(m) (m)
b2 − a2,1 x1 − a2,3 x3 − ... − a2,n xn
(m)
) 


...

xi
( m +1)
= (
1
ai ,i
(m) (m)
bi − ai ,1 x1 − ai , 2 x2 − ... − ai ,n xn
(m)
) 

... 

xn
( m +1)
= (
1
a2 , 2
(m) (m)
b2 − an ,1 x1 − an , 2 x2 − ... − an ,n −1 xn −1
(m)
)

Matricialmente: Descomponemos la matriz A=D+L+R
con D= matriz con sus términos 0 excepto la diagonal en
la que coincide con A. L con sus términos cero excepto
los que están por debajo de la diagonal y R lo mismo con
los que están por encima de la diagonal coincidentes con A.
54
 Expresión matricial del método de Jacobi:
Ax=b Æ (D+L+R)x=b Æ Dx=-(L+R)x+b Æ
-1 -1
x=-D (L+R)x+D b
-1 -1
Por tanto: BJ = =-D (L+R) y CJ= D b

Ejemplo: Hallar la matriz del método de Jacobi y estudiar

la convergencia del mismo para el sistema Ax=b. Dar 2 itera-
ciones iniciando en el vector nulo.  3 2 −1
 
3
 
A = 1 4 2  b= 9 
 2 2 − 5  − 1
   
Descomponemos: 3 0 0 
 
0 0

0

 0 2 − 1
 
D = 0 4 0  L = 1 0 0 R = 0 0 2 
 0 0 − 5 2 2 0  0 0 0 
    
1 0   0 2 − 1  0 −2 1 
 3 0   3 3
Calculo BJ: BJ = -D-1 (L + R) = -  0 1 0   1 0 2  =  −1 0 −1 
 4   4 2
 0 0 − 1   2 2 0   2 2 55 0 
 5  5 5 
  
1 0 0  3   1 
 3    
−1  1   9
CJ = D b = 0 0  9 =
 4    4 
 0 0 − 1  − 1  1 
 5  
5
Para estudiar la convergencia del método, calculo el radio espectral de BJ :
−λ −2 1
 0.5355
3 3
−1 1 1 
− λ − 1 = 0 ⇒ -λ 3 + λ + = 0 ⇒ λ = − 0.2677 + 0.3392i ⇒
4 2 10 10
2 2 − 0.2677 − 0.3392i
−λ 
5 5
ρ(BJ ) ≈ 0.5355 < 1 ⇒ El método de Jacobi es convergente

3− 2y + z  ( m +1) 3 − 2 y ( m) + z ( m) 
x= x = 
3x + 2 y − z = 3  3  3
(m) 
 9 − x − 2 z  ( m +1) 9 − x (m)
− 2 z 
x + 4 y + 2z = 9  ⇒ y =  ⇒ y = 
4 4
2 x + 2 y − 5 z = −1 −1− 2x − 2 y 

1 + 2 x (m)
+ 2 y (m) 
z= z ( m +1)
= 
−5  5  56
  3 − 2( 9 ) + ( 1 ) 
 3 − 2(0) + (0)     4 5   − 13 
     3   
 0  3  1  1 )   30 
  9 − 0 − 2( 0)   9  9 − 1 − 2(
x
(0)
=  0 ⇒ x
(1)
=  = ⇒x
( 2)
=  5  =  19 
 4   4  4   20 
 0 
   1 + 2(0) + 2(0)   1  1 + 2 (1) + 2 ( 9 )   3 
     4   
 5  5  5   2 
 

El vector residuo: r=Ax-b y su norma infinita tras las dos

iteraciones es:

  − 13   19  3 
 3  + 2  − − 3   − 68 
  30   20  2   15 
  − 13   19   3    − 39 
r = Ax − b =    + 4  + 2  − 9  = ⇒ r = 6.1
  ∞
  30   20   2    −10 
 2 − 13  + 2 19  − 5 3  + 1  61 
  30   20   2    10 
       
57
 Método de Gauss-Seidel:
Pretende acelerar el método anterior empleando la
versión más actualizada de la variable en cada cálculo.
Continuando con el ejemplo:
3− 2y + z  ( m +1) 3 − 2 y (m)
+ z (m)

x= x = 
3x + 2 y − z = 3  3  3
(m) 
 9 − x − 2 z  ( m +1) 9 − x ( m +1)
− 2 z 
x + 4 y + 2z = 9  ⇒ y =  ⇒ y = 
4 4
2 x + 2 y − 5 z = −1 
−1 − 2x − 2 y  1 + 2x ( m +1)
+ 2y ( m +1) 
z= z ( m +1)
= 
−5  5 

 3 − 2(2) + ( 7 ) 
 5 
 3 − 2(0) + (0)   3   2 

3
    2   15 
 0   1  9 − − 2( 7 5 )   91 
  9 − 1 − 2(0)   
 = 
(0) (1) ( 2)
x =  0 ⇒ x = = 2 ⇒x = 15
 4  7  4   60 
 0  1 + 2(1) + 2(2)   
   1 + 2( 2 ) + 2(91 )   43 
  5  15 60   
5 50
   5   
  58
 
   2   91  43 
 3  + 2  − − 3   − 32 
  15   60  50   75 
  2   91   43   − 27 
r = Ax − b =    + 4  + 2  − 9  =  ⇒ r = 27/25 = 1.08
  ∞
  15   60   50    25 
 2 2  + 2 91  − 5 43  + 1  0 
  15   60   50    
       

Las fórmulas para iterar en el caso general son:

( m +1) 1 (m) (m) (m) 
x1 = (b1 − a1, 2 x2 − a1,3 x3 − ... − a1,n xn ) 
a1,1

1
− a2,3 x3 − ... − a2,n xn ) 
( m +1) ( m +1) (m) (m)
x2 = (b2 − a2,1 x1
a2 , 2 

( m +1) 1 ( m +1) ( m +1) (m) 
x3 =+ (b3 − a3,1 x1 − a3, 2 x2 − ... − a3,n xn ) 
a3,3

... 
( m +1) 1 ( m +1) ( m +1) ( m +1) 
xn = (b2 − an ,1 x1 − a n , 2 x2 − ... − an ,n −1 xn −1 )
an , n 

59
Matricialmente:
Ax=b Æ (D+L+R)x=b Æ (D+L)x=-Rx+b Æ
-1 -1
x= -(D+L) Rx+(D+L) b
por lo que la matriz del método de Gauss-Seidel es:
-1 -1
BGS= -(D+L) R y CGS =(D+L) b
Como es de esrperar la condición necesaria y sufi-
ciente de convergencia es que ρ(BGS ) < 1

Comparaciones:
El método de Gauss-Seidel es en general más rápido que
el de Jacobi.
Sin embargo el método de Jacobi permite paralelismo y
el de Gauss-Seidel no.
También puede ocurrir que uno sea convergente y el
otro no lo sea. 60
 Método de sobre-relajación (SOR):
Consiste en añadir un parámetro w de relajación que ace-
lere la convergencia del método de Gauss-Seidel.
Se ha comprobado que usualmente los pasos que da el mé-
todo de Gauss-Seidel son demasiado cortos y lo que este
método hace es agrandar dichos pasos.
Esto se hace componente a componente. Así, si una com-
ponente en la iteración m vale ‘a’ y cuando aplicamos la
fórmula de Gauss-Seidel pasaría a valer ‘b’, el método SOR
le asigna en la iteración m+1 el valor a +w(b-a)
Este sería el valor que se emplea para calcular las
sucesivas componentes.
NOTA: El método SOR con w=1 es el de Gauss-Seidel.
Si w<1 (sub-relajación) y si w>1 (sobre-relajación).
61
 Fórmulas para iterar:
( m +1) (m) ( m +1) ( m) ( m +1) (m) (m) ( m +1)
∆ SOR = w∆ GS = w( xi ,GS − xi ) ⇒ xi , SOR = xi + w( xi ,GS − xi ) = (1 − w) xi + wxi ,GS ⇒
( m +1) w
(m) (m) (m) (m) 
x1 = (1 − w) x1 + (b1 − a1, 2 x2 − a1,3 x3 − ... − a1,n xn ) 
a1,1

w
− a2,3 x3 − ... − a2,n xn ) 
( m +1) (m) ( m +1) (m) (m)
x2 = (1 − w) x2 + (b2 − a2,1 x1
a2 , 2 

( m +1) (m) w ( m +1) ( m +1) (m) 
x3 = (1 − w) x3 + (b3 − a3,1 x1 − a3, 2 x2 − ... − a3,n xn ) 
a3,3

... 
( m +1) ( m) w ( m +1) ( m +1) ( m +1) 
xn = (1 − w) xn + (b2 − an ,1 x1 − a n , 2 x2 − ... − an ,n −1 xn −1 )
an , n 


Matricialmente: /
Ax = b ⇒ w( D + L + R ) x = wb ⇒ (1 − w) Dx + w( D + L + R) x = (1 − w) Dx + wb ⇒
( D + wL + wR ) x = (1 − w) Dx + wb ⇒ ( D + wL ) x = − wRx + (1 − w) Dx + wb ⇒
( D + wL) x = [(1 − w) D − wR ]x + wb ⇒ x = ( D + wL ) −1[(1 − w) D − wR ]x + ( D + wL) −1 wb

Luego : x ( m +1) = ( D + wL) −1[(1 − w) D − wR ]x ( m ) + ( D + wL) −1 wb

donde la matriz del método es : ( D + wL) −1[(1 − w) D − wR ] 62
Ejemplo:
Dar 2 iteraciones por SOR con w=1.2, al sistema:
3x+2y-z=3; x+4y+2z=9; 2x+2y-5z=-1, iniciando en (0,0,0).

3 − 2 y + z   ( m +1) 3 − 2 y (m) + z (m) 

x=  x = (1 − w) x + w
(m)

3x + 2 y − z = 3  3  3
  ( m +1) (m) 
 9 − x − 2 z   ( m +1) 9 − x − 2 z 
x + 4 y + 2z = 9  ⇒ y = ⇒ y = (1 − w) y ( m ) + w ⇒
 4   4
2 x + 2 y − 5 z = −1 ( m +1) 
1 + 2 x + 2 y   ( m +1) 1 + 2 x ( m + 1)
+ 2 y 
z= z = (1 − w) z (m)
+ w
5   5 


 3 − 2( 0) + 0 
 − 0. 2 ( 0 ) + 1 . 2 
 x   0  x  
(0) (1)
3   1.2 
 ( 0 )     (1)  9 − 1.2 − 2(0)   
 y  =  0 ⇒  y  =  − 0.2(0) + 1.2 =  2.34 
 ( 0 )     (1)   4 
z 0 z   1.9392 
      − 0.2(0) + 1.2 1 + 2 (1. 2 ) + 2( 2. 34 )  
 
 5 
63
  3 − 2(2.34) + 1.9392 
 − 0 . 2(1 .2 ) + 1 .2 
x  
( 2)
3   − 0.1363 
 ( 2)  9 − (−0.1363) − 2(1.9392)   

 y  = − 0.2(2.34) + 1.2 =  1.1094  ⇒
 z ( 2)   4 
   0.3192 
  − 0.2(1.9392) + 1.2 1 + 2 ( −0 .1363 ) + 2 (1 .1094 )  
 
 5 

 3(−0.1363) + 2(1.1094) − (0.3192) − 3   − 1.5094 

   
res =  (−0.1363) + 4(1.1094) + 2(0.3192) − 9  =  − 4.0604  ⇒ res ∞ = 4.0604
 2(−0.1363) + 2(1.1094) − 5(0.3192) + 1  1.35 
   

64
Condiciones de convergencia para los métodos
iterativos usuales:
1) La condición necesaria y suficiente de convergencia
será que el radio espectral de la matriz del método sea
estrictamente menor que la unidad. Asi:
Jacobi converge si y solo si : ρ(-D-1 (R + L)) < 1
Gauss - Seidel converge si y solo si : ρ(-(D + L)-1 R)) < 1
Para Jacobi y Gauss-Seidel obtenemos:
2) Jacobi converge si y solo si ningún término de la dia-
gonal principal de A es nulo y
Las soluciones de λD + L + R = 0 verifican λ i < 1
3) Gauss-Seidel converge si y solo si ningún término de
la diagonal principal de A es nulo y
Las soluciones de λD + λL + R = 0 verifican λ i < 1
65
 Demostración para Jacobi:
ρ(BJ ) < 1 ⇔ Todos los autovalores de BJ son menores que 1 ⇔ - D -1 (L + R) - λI = 0
tiene sus raices verificando λ i < 1 pero : - D -1 (L + R) - λI = - D -1 (L + R) - λD -1D =
= - D -1 (L + R + λD) = - D -1 | | L + R + λD luego como - D -1 ≠ 0, la ecuación
- D -1 (L + R) - λI = 0 tiene las mismas raices que : | L + R + λD |= 0 y Jacobi
converge si son menores que 1.

Demostración para Gauss-Seidel:

ρ(BGS ) < 1 ⇔ Todos los autovalores de BGS son menores que 1 ⇔ - (D + L)-1 R - λI = 0
tiene sus raices verificando λ i < 1 pero : - (D + L)-1 R - λI = - (D + L)-1 R - λ(D + L)-1 (D + L) =
= - (D + L)-1 (R + λ(L + D)) = - (D + L)-1 | | (R + λ(L + D)) luego como - (D + L)-1 ≠ 0, la ecuación
- (D + L)-1 R - λI = 0 tiene las mismas raices que : | λL + R + λD |= 0 y
Gauss - Seidel converge si son menores que 1.

66
 Ejemplos:
Estudiar la convergencia de los métodos de Jacobi y Gauss-
Seidel para los sistemas: Ax=b y Cx=d siendo:
 3 −1 2  − 1  3 −1 1 1
       
A =  4 2 5  b =  2  C =  4 9 5 d =  5 
 0 −1 3 3  1 −1 3 7
       

3λ − 1 2
Para el sistema Ax = b y Jacobi : λD + L + R = 0 ⇒ 4 2λ 5 = 0 ⇒
0 − 1 3λ
18λ 3 − 8 + 15λ + 12λ = 0 ⇒ 18λ 3 + 27λ - 8 = 0 ⇒ λ1 = -0.1407 + 1.2488i;
λ 2 = -0.1407 - 1.2488i; λ 3 = 0.2814 ⇒| λ1 |=| λ 2 |= 1.2567 > 1

Por tanto Jacobi no converge.

67
 3λ −1 2
Para el sistema Ax = b y Gauss - Seidel : λD + λL + R = 0 ⇒ 4λ 2λ 5 = 0 ⇒
0 − λ 3λ
3 −1 2 3 −1 2
λ4 2λ 5 = 0 ⇒ λ 2 4 2λ 5 = 0 ⇒ λ1 = λ 2 = 0 (doble)
0 − λ 3λ 0 −1 3
19
18λ − 8 + 15 + 12 = 0 ⇒ 18λ + 19 = 0 ⇒ λ 3 = - ⇒| λ 3 |> 1
18
Por tanto Gauss-Seidel no es convergente.
3λ − 1 1
Para el sistema Cx = d y Jacobi : λD + L + R = 0 ⇒ 4 9λ 5 = 0 ⇒
1 − 1 3λ
81λ 3 − 5 - 4 − 9λ + 15λ + 12λ = 0 ⇒ 81λ 3 + 18λ - 9 = 0 ⇒ λ1 = -0.1667 + 0.5528i;
λ 2 = - 0.1667 - 0.5528i; λ 3 = 0.3333 ⇒| λ1 |=| λ 2 |≈ 0.5774 < 1

Por tanto Jacobi es convergente.

68
 3λ − 1 1
Para el sistema Cx = d y Gauss - Seidel : λD + λL + R = 0 ⇒ 4λ 9λ 5 = 0 ⇒
λ − λ 3λ
3 −1 1
λ 4 9λ 5 = 0 ⇒ λ1 = 0, 81λ 2 − 5 - 4λ - 9λ + 15λ + 12λ = 0 ⇒ 81λ 2 + 18λ - 5 = 0 ⇒
1 − λ 3λ
λ 2 = -0.3833, λ 3 = 0.1611 ⇒| λ i | ≤ 0.3833 < 1

Por tanto Gauss-Seidel es convergente.

NOTA: Pronto veremos que el valor del radio espectral
de la matriz del método, no solo asegura la convergencia
sino que da una medida de la velocidad de la misma. Por
tanto para el sistema Cx=d, en el que resultan ambos
métodos convergentes, el método de Gauss-Seidel es más
rápido.
69
 4) Si el método SOR es convergente para un valor w,
entonces es convergente para cualquier otro valor w’
verificando: w’<w. Análogamente si diverge para un w,
divergerá para un w’ que verifique w’> w.

5) Si la matriz A es de diagonal estrictamente dominante

convergen los métodos de Jacobi, Gauss-Seidel y SOR
(w61).
Ejemplo: Para los sistemas Ax=b y Cx=d, los métodos
de Jacobi y Gauss-Seidel son convergentes.
 4 −1 2   − 1  3 −1 1  1
       
A =  4 10 5  b =  2  C =  4 − 9 − 4  d =  5 
 2 −1 − 4 3  1 −1 3  7
       
70
6) Si A es simétrica definida positiva entonces converge
el método SOR para 0<w<2. Para w=1 es Gauss-Seidel.

7) Si D-L-R es simétrica definida positiva entonces

converge Jacobi.

8) Si A es tridiagonal y simétrica definida positiva en-

tonces convergen los métodos de Jacobi y Gauss-Seidel,
verificandose además:
- El valor óptimo para w en el método SOR es:
w=
2
⇒ ρ(BSOR ) = w - 1 /
1 + 1 - ρ(BJ ) 2
- Los radios espectrales de las matrices de los métodos
de Gauss-Seidel y Jacobi verifican: ρ(BGS ) = ρ(BJ ) 2
71
Errores en la solución de un sistema lineal:
Al efectuar cualquier operación aritmética en la resolución
de un sistema aparecen errores que se propagarán a poste-
riores cálculos incrementando su efecto. Por otra parte los
métodos iterativos, aunque son insensibles a la propagación
del error, necesitarían de infinitas iteraciones.
Por tanto sea cual sea el método empleado obtendremos
una solución aproximada.
Acotación del error al resolver un sistema:
Al intentar resolver Ax=b, obtenemos x’ que verifica:
Ax’=b+r luego el residuo r=Ax’-Ax=A(x’-x).
∆x = x'-x ⇒ r = A ∆x ⇒ r = A ∆ x ≤ A ∆x
−1 −1
r = A ∆x ⇒ A r = ∆x ⇒ ∆ x = A r ≤ A -1
r
/
72
Luego obtenemos la acotación del error cometido:
r
≤ ∆x ≤ A −1 r
A
Análogamente: b
Ax = b ⇒ b = Ax ≤ A x ⇒ x ≥
A
−1 −1 −1 −1
/
x =A b⇒ x = A b ≤ A b ⇒ x ≤ A b

r
y por tanto: 1 r A ∆x A −1 r r
= −1 ≤ ≤ = c(A) /
c(A) b A b x b b
A

1 r r ∆x
≤ δx ≤ c(A) con δx = , c(A) = A A -1

c(A) b b x
73
 Número de condición de una matriz:
Al producto de la norma de una matriz por la norma de su
inversa se denomina número de condición
−1
c(A) = A A
y da una medida de los errores que se van a producir al
resolver un sistema con dicha matriz.
Propiedades:
1) El número de condición siempre es mayor que 1.
c(A) = A A −1 ≥ AA −1 = I = 1
2) Cuanto mayor es c(A), peor condicionada estará el
sistema. Puede considerarse que para una matriz con
c(A)=10 la última cifra de la solución no será fiable, para
c(A)=100, las dos últimas, etc,... 74
 Ejemplo:
Resolver el sistema Ax=b con 1
para n=10 y acotar el error A i, j = , bi = 1
cometido en cada variable.
i + j −1
Con MATLAB resulta: x=(-9.99629431, 989.67863835,
-23753.13525466, 240177.44656281, -1260961.06074487,
3782956.91108005, -6725367.75688302, 6999971.72699231,
-3937532.49575707, 923628.67899107) que produce un vector
residuo: r=1e-9*(0.01455191522837, -0.04365574568510
-0.05820766091347, -0.08731149137020, 0.02910383045673
-0.05093170329928, 0.12369127944112, -0.05820766091347
0.08003553375602, -0.14551915228367)
luego: r ∞ = 1.455191522836685 10-10
Calculando la inversa y su norma infinita resulta:
A −1 = 1.207036531817908 1013 ⇒ ∆x ∞
≤ A -1 r ∞
= 1756.469328855612
75
∞ ∞
 Esto significa que el valor obtenido en cada incógnita
oscila en más menos 1756.4, lo que deja algunas incóg-
nitas sin cifras significativas.
La matriz del problema propuesto es conocida como
matriz de Hilbert y está muy mal condicionada:
c(A)=3.535371683074594e+013 (norma infinita)

Ejemplo: Lo mismo para el sistema: 0.5x-1/3y=0

0.5x-0.333333y=0.1;
Resolviendo: x=(2.00000000005351, 3.00000000008027)
r=1.0e-010*(0.05820760540232, 0.14551915228367)
r ∞
≈ 1.455 10-11; A -1 ≈ 6000000 ⇒ ∆x ≤ 0.0000873
∞

es decir tenemos 3 decimales exactos.

NOTA: c(A)w5000000 (sistema mal condicionado) 76
 Acotación del error en los métodos
iterativos:
Dado el método iterativo: x ( m +1) = Bx ( m ) + C la solu-
ción exacta x verificará: x=Bx+C.
El error tras m+1 iteraciones verificará:
∆( m +1) = x ( m +1) − x = Bx ( m ) + C − Bx + C = B( x ( m ) − x) ≤ B x ( m ) − x = B ∆( m )

(Esto significa que en cada iteración el error se reduce, al

menos, en el factor ||B||).
Además:
( m +1) 2 ( m −1) 3 ( m−2) m +1
∆ ≤ B∆ (m)
≤ B ∆ ≤ B ∆ ≤ ... ≤ B ∆(0)

(Esto acota el error en función del error inicial).

77
 x ( m +1) − x ( m ) = Bx ( m ) + C − Bx ( m −1) − C = B ( x ( m ) − x ( m −1) ) ≤ B x ( m ) − x ( m −1)
2 3 m
≤ B x ( m −1) − x ( m − 2) ≤ B x ( m − 2 ) − x ( m −3) ≤ ... ≤ B x (1) − x ( 0 )

(Los vectores se aproximan entre sí con el factor ||B||). /

∆( m ) = x ( m ) − x ≤ x ( m +1) − x ( m ) + x ( m + 2 ) − x ( m +1) + x ( m +3) − x ( m + 2 ) + ... ≤
m m +1 m+ 2
≤ B x (1) − x ( 0 ) + B x (1) − x ( 0 ) + B x (1) − x ( 0 ) + ...
{Suma de infinitos términos de una progresión geométrica ilimitada con razón < 1}
Obtenemos la fórmula: /
m
B x (1) − x ( 0 )
∆( m ) ≤
1− B
que expresa el error tras m iteraciones en función de la
norma de B y la diferencia en las iteraciones iniciales.
78
 Ejemplos:
Acotar el error cometido al dar 100 iteraciones por los
métodos de Jacobi y de Gauss-Seidel al sistema (n=100)
Ax=b (iniciando en el vector 0) siendo:
 2 (i = j = 1) ∨ (i = j = 100)
4 (i = j) ∧ (1 < i < 100)

Ai , j =  bi = i
− 1 | i - j |= 1
 0 cualquier otro caso
El sistema es tridiagonal con matriz A tridiagonal y
 2 − 1 0 ... 0 0  con diagonal estrictamente domi-
 
 − 1 4 − 1 ... 0 0  minante por lo que convergen los
 0 − 1 4 ... 0 0 
A=  métodos de Jacobi y Gauss-Seidel.
 ... ... ... ... ... ... 
0 0 0 ... 4 − 1

0 0 0 ... − 1 2  79

 La matriz del método de Jacobi es:
 1 
0 0 ... 0 0
 2 
1 0
1
... 0 0
4 4 
 1 
BJ =  0 4
0 ... 0 0  ⇒ ρ(B ) ≤ B
J J = 0.5 ⇒ ρ(BGS ) = ρ(BJ ) 2 ≤ 0.25
  ∞

 ... ... ... ... ... ... 

0 1
0 0 ... 0
 4
 1 
0 0 0 ... 0
 2 
 1 
 
 2 
 2 
Jacobi:  4  0.5100 50
x (1) =  3  ⇒ x − x = 50 ⇒ x
(1) (0) (100)
−x ≤ ≈ 7.889 10-29
  1 - 0.5
 4 
 ... 
 100  80
 2 
 El cálculo para Gauss-Seidel es engorroso y da:
 0.5 
 
 0.625  100
0. 25 66.4444
x (1) =  0.9063  ⇒ x − x = 66.4444 ⇒ x
(1) (0) (100)
−x ≤ ≈ 5.5131 10-59
  1 - 0.25
 ... 
 66.4444 
 

NOTA: Donde se ha usado el valor 0.25 como norma

infinita de la matriz del método de Gauss-Seidel que es su
verdadero valor, aunque hasta ahora solo tenemos el radio
espectral que es una cota inferior.

81
 Acotación del error cuando existen
errores en los datos:
Suponemos ahora que los coeficientes del sistema están
afectados de errores y queremos estimar como se propaga-
rán al resultado.
Estamos resolviendo en realidad el sistema:
(A + ∆A)(x + ∆x) = b + ∆b
del que obtenemos la solución exacta x’, pero en realidad
queremos resolver Ax=b, entonces el error ∆x = x'-x
verifica:
(A + ∆A)(x + ∆x) = b + ∆b ⇒ Ax + ∆A x + A∆ x + ∆A∆x = b + ∆b ⇒
{despreciando el término ∆A∆x} ⇒ A∆ x + ∆A x = ∆b ⇒ /

∆x ≈ A -1 (∆ b - ∆A x)
82
 Para el error relativo tenemos el siguiente teorema:
TEOREMA:
∆x ∆b ∆A
Llamando : rel(x) = , rel(b) = , rel(A) = :
x b A
se obtiene la fórmula:
rel(A) + rel(b)
rel(x) ≤ c(A)
1 - c(A)rel(A)
Nota:
Las fórmulas obtenidas para los errores son válidas para
cualquier norma (siempre la misma).
La más usada es la norma infinito.
83
 Ejemplo: Dado el sistema Ax=b donde los coeficientes
vienen afectados de errores, encontrar intervalos para las
incógnitas.
(3 ± 0.(201 + ± − ± = 14 ± 0.005 
(3 ± 0.01) x + (4 ± 0.01) y − (1 ± 0.001) z = 14 ± 0.005

) x ( 4
± 0.001) x − (1 ± 0.001) y0 .01 ) y = 0 ± 0.005  0.001) z
( 1 

(1 ± 0.001) x + (2 ± 0.001) y + (2 ± 0.01) z = −1 ± 0.005 

(2 ± 0.001  3 )4 x −−

1 (1 ±

 140
 
 .001 ) y 0.01 0.001  0.005= 0 ± 0.005 
 0.01

 −1 ± 0.005
A =  2 − 1 0  ; b =  0 ; ∆A =  0.001 0.001 0 ; ∆b =  0.005  ⇒
(1 ± 0.001  1 )2x +
 2 ( 2 ±

 − 10
 
 .001 ) y
 0.001
 + (
0.0012 ±
0.01 0
 .01

)
 z
0.005= 
1
  ∆A ∞ 0.021
3 4 − 1 o el sistema
resolviend  14  Ax = b : x=0.01
 2  ⇒ x 0.01
= 3, 0.001
rel(A) =  = =0.005 
0.002625
      - 3 ∞
A ∞ 8  
A =  2 −1 0  ; b =  0 ; ∆A =  0.001 0.001
  0 ; ∆b =  0.005  ⇒
1 2 2  ∆b ∞  −0.005 1  0.001-1 0.001  0.03704  0.005
0.07407 0.37037
0.01 
 =
rel(b)   ≈ 0.00035714
=  ; A =  0.14815 - 0.25926 0.07407 ⇒ A−1 = 0.66667
b∞ 14  - 0.18519 0.07407 0.40741 
1  
−1   ∆ A 0.021
resolviendo el ⇒ c(A) = A ∞ A
sistema Ax = ∞b=:8(0x.66667=  2) ≈5.⇒33333
x ∞ = 3, rel(A) = ∞
= = 0.002625
rel(x) ≤ c(A)
rel(A) + rel(b)
= 5.33333 0.002625 + 0.00035714
=
A ∞
8
1 − c(A)rel(A)  - 3  1 − 5.33333(0.002625)
 0.07407 0.37037 0.03704 
∆b 0.005  
rel(b) = ∞
= ≈ 0.00035714; A =  0.14815 - 0.25926 0.07407  ⇒ A−1 = 0.66667
-1

b ∞
14  - 0.18519 0.07407 0.40741 
 
⇒ c(A) = A ∞
A −1 = 8(0.66667) ≈ 5.33333 84
∞
 rel(A) + rel(b) 0.002625 + 0.00035714
rel(x) ≤ c(A) = 5.33333 ≈ 0.01613
1 − c(A)rel(A) 1 − 5.33333(0.002625)

Luego cada incognita tiene, al menos, 2 cifras significativas.

∆x
rel(x) = ∞
≤ 0.01613 ⇒ ∆x ∞
≤ 3(0.01613) = 0.04239 < 0.05
x ∞

Luego cada incógnita tiene al menos 1 decimal exacto, es

decir:
0.957 6x6 1.043
1.957 6y6 2.043
-3.043 6z6 -2.957

85
Inversión por vectores columnas:
El cálculo de la matriz inversa puede realizarse columna
a columna resolviendo el sistema: Aci = e i
c i = (columna i - ésima de la inversa ), ei = (0,0,0,...,1,0,...,0,0)'

donde el único 1 del vector e estará en la posición i.

i

Así calcular la inversa de una matriz equivaldrá a resol-

ver n sistemas con la misma matriz A y resultará apropia-
do cualquier método como Crout, Cholesky o QR que
guarde la descomposición de la matriz A.

86
 Ejemplo:
 3 4 - 1
 
Hallar la inversa de la matriz A =  2 − 1 0  mediante vectores columnas :
1 2 2 
 
Aunque no es el método usual emplearemos Gauss, para
mostrar como puede realizarse una rutina que resuelva
varios sistemas simultáneamente (en eso consiste el cálculo
de la inversa).
 3 4 -1 1 0 0   2 
  f
 2 ← f 2 − f1 
Partimos de : (A | I) =  2 - 1 0 0 1 0  y triangularizamos :  3 →
1 
 1 2 2 0 0 1  f 3 ← f 3 − f1 
   3 
 3 4 -1 1 0 0   2   3 4 -1 1 0 0 
   3  
 0 - 11/3 2/3 - 2/3 1 0  →  f 3 ← f 3 − -11 f 2  →  0 - 11/3 2/3 - 2/3 1 0 
 0 2/3 7/3 - 1/3 0 1   3   0 0 27/11 - 5/11 2/11 1 
 
87
 3 4 -1 1 0 0 
 11     −3 2 
Sustitución hacia atrás :  f 3 ← f 3 *  →  0 - 11/3 2/3 - 2/3 1 0  →  f2 ←  f 2 − f 3 
 27    11  3 
0 0 1 - 5/27 2/27 11/27 
 3 4 -1 1 0 0   1 0 0 2/27 10/27 1/27   2/27 10/27 1/27 
   1     
→  0 1 0 4/27 - 7/27 2/27  →  f1 ← ( f1 − 4 f 2 + f 3 ) →  0 1 0 4/27 - 7/27 2/27  ⇒ A−1 =  4/27 - 7/27 2/27 
 0 0 1 - 5/27 2/27 11/27   3     - 5/27 2/27 11/27 
   0 0 1 - 5/27 2/27 11/27   

 0.07407407 407407 0.37037037 037037 0.03703703 703704 

 
A =  0.14814814 814815 - 0.25925925 925926 0.07407407 407407 
-1

 - 0.18518518 518519 0.07407407 407407 0.40740740 740741 

 

NOTA: Las rutinas numéricas para resolver sistemas linea-

les usan como parámetros de llamada dos matrices A y B, en
lugar de b vector, devolviendo como resultado una matriz X
tal que AX=B.
De esta forma X=rutina(A,B) devuelve la inversa de A
cuando la llamamos con B=I (matriz identidad).
88
 Acotación de los autovalores mediante los
círculos de Gerschgorin:
Para cada fila de la matriz A, encontraremos un autovalor
de A en el círculo de centro el término diagonal ai,i y radio
la suma de los restantes términos en valor absoluto.
Ejemplo:
Acotar los autovalores de la matriz A:
 4 − 1 0 1   centro = 4, radio = 2 
   
 1 − 5 2 0   centro = -5, radio = 3 
A=  ⇒  
0 −1 3 1 centro = 3, radio = 2 
  
 1 0 0 − 4  centro = −4, radio = 1
   

NOTA: Los autovalores de A son los mismos que los de A’.

89
 3

-1

-2

-3
-8 -6 -4 -2 0 2 4 6

En cada círculo habrá un autovalor, en efecto estos son:

-4, 4, -4.6056 y 2.6056 (marcados con x en la figura).
NOTA: El método usual para cálculo computacional de
autovalores se basa en la descomposición QR de A y con-
siste en obtener una sucesión de matrices equivalentes
mediante: [Q,R]=qr(Ai), Ai+1=R*Q 90