Ap Algnum

E.T.S.
DE INGENIERÍA INFORMÁTICA
Apuntes de
ÁLGEBRA NUMÉRICA
para la titulación de
INGENIERÍA INFORMÁTICA
Curso 2002-2003
por
Fco. Javier Cobos Gavala
DEPARTAMENTO DE
MATEMÁTICA APLICADA I
Contenido
1 Ecuaciones no lineales 1
1.1 Errores y condicionamiento en problemas numéricos . . . . . . 2
1.2 Método y algoritmo de la bisección: análisis de errores . . . . 4
1.2.1 Algoritmo . . . . . . . . . . . . . . . . . . . . . . . . . 5
1.3 Punto fijo e iteración funcional . . . . . . . . . . . . . . . . . . 8
1.3.1 Cota del error “a posteriori” . . . . . . . . . . . . . . . 9
1.4 Análisis del método de Newton-Raphson . . . . . . . . . . . . 12
1.4.1 Algoritmo . . . . . . . . . . . . . . . . . . . . . . . . . 14
1.4.2 Análisis de la convergencia: Regla de Fourier . . . . . . 15
1.4.3 Método de Newton para raı́ces múltiples . . . . . . . . 18
1.5 Un problema mal condicionado: ceros de un polinomio . . . . 20
1.5.1 Sucesiones de Sturm . . . . . . . . . . . . . . . . . . . 23
1.5.2 Algoritmo de Horner . . . . . . . . . . . . . . . . . . . 26
1.6 Sistemas de ecuaciones no lineales . . . . . . . . . . . . . . . . 28
1.6.1 Método de Newton . . . . . . . . . . . . . . . . . . . . 30
1.7 Ejercicios propuestos . . . . . . . . . . . . . . . . . . . . . . . 33
2 Sistemas de ecuaciones lineales 39

2.1 Normas vectoriales y matriciales . . . . . . . . . . . . . . . . . 39
2.1.1 Normas vectoriales . . . . . . . . . . . . . . . . . . . . 39
2.1.2 Distancia inducida por una norma . . . . . . . . . . . . 40
2.1.3 Convergencia en espacios normados . . . . . . . . . . . 41
i
ii Contenido
2.1.4 Normas matriciales . . . . . . . . . . . . . . . . . . . . 41

2.1.5 Transformaciones unitarias . . . . . . . . . . . . . . . . 43
2.1.6 Radio espectral . . . . . . . . . . . . . . . . . . . . . . 44
2.2 Sistemas de ecuaciones lineales . . . . . . . . . . . . . . . . . . 44
2.2.1 Número de condición . . . . . . . . . . . . . . . . . . . 46
2.3 Factorización LU . . . . . . . . . . . . . . . . . . . . . . . . . 53
2.4 Factorización de Cholesky . . . . . . . . . . . . . . . . . . . . 56
2.5 Métodos iterados . . . . . . . . . . . . . . . . . . . . . . . . . 59
2.5.1 Método de Jacobi . . . . . . . . . . . . . . . . . . . . . 63
2.5.2 Método de Gauss-Seidel . . . . . . . . . . . . . . . . . 64
2.5.3 Métodos de relajación (SOR) . . . . . . . . . . . . . . 64
2.6 Métodos del descenso más rápido y del gradiente conjugado . 65
2.6.1 Método del descenso más rápido . . . . . . . . . . . . . 67
2.6.2 Método del gradiente conjugado . . . . . . . . . . . . . 67
3 Sistemas inconsistentes y sistemas indeterminados 73

3.1 Factorizaciones ortogonales . . . . . . . . . . . . . . . . . . . . 73
3.2 Interpretación matricial del método de
Gram-Schmidt: factorización QR . . . . . . . . . . . . . . . . 73
3.3 Rotaciones y reflexiones . . . . . . . . . . . . . . . . . . . . . 75
3.4 Transformaciones de Householder . . . . . . . . . . . . . . . . 76
3.4.1 Interpretación geométrica en Rn . . . . . . . . . . . . . 77
3.4.2 Householder en Cn . . . . . . . . . . . . . . . . . . . . 78
3.4.3 Factorización QR de Householder . . . . . . . . . . . . 80
3.5 Sistemas superdeterminados. Problema de los mı́nimos cuadrados 85
3.5.1 Transformaciones en sistemas superdeterminados . . . 87
3.6 Descomposición en valores singulares y seudoinversa de Penrose 90
3.6.1 Seudoinversa de Penrose . . . . . . . . . . . . . . . . . 91
Contenido iii
4 Autovalores y autovectores 103

4.1 Conceptos básicos . . . . . . . . . . . . . . . . . . . . . . . . . 103
4.2 Método interpolatorio para la obtención del polinomio carac-
terı́stico . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105
4.3 Método de la potencia y variantes . . . . . . . . . . . . . . . . 106
4.4 Cociente de Rayleigh . . . . . . . . . . . . . . . . . . . . . . . 110
4.5 Teorema de Gerschgorin . . . . . . . . . . . . . . . . . . . . . 111
4.6 Sensibilidad de los autovalores a las transformaciones de seme-
janza . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112
4.7 Métodos iterados para la obtención de autovalores y autovectores118
4.7.1 Algoritmo QR de Francis . . . . . . . . . . . . . . . . . 118
4.7.2 Método de Jacobi para matrices simétricas reales . . . 121
4.8 Reducción del problema a matrices hermı́ticas . . . . . . . . . 126
4.9 Reducción del problema a matrices simétricas reales . . . . . . 128
4.10 Aplicación al cálculo de las raı́ces de un polinomio . . . . . . . 129
Índice 139
1. Ecuaciones no lineales
Dada una función no nula f : C → C, resolver la ecuación f (x) = 0 es hallar

los valores x que anulan a dicha función. A estos valores x se les denomina
raı́ces o soluciones de la ecuación, o también, ceros de la función f (x).
Los métodos de resolución de ecuaciones y sistemas de ecuaciones se clasifican
en directos e iterados . Los del primer grupo nos proporcionan la solución me-
diante un número finito de operaciones elementales, mientras que los iterados
producen una sucesión convergente a la solución del problema.
Un ejemplo de método directo es la conocida fórmula de resolución de las
ecuaciones de segundo grado ax2 + bx + c = 0, cuyas soluciones vienen dadas
por la fórmula √
−b ± b2 − 4ac
x=
2a
Sin embargo, el siglo pasado Abel probó que no existe ninguna fórmula equi-
valente (en término de raı́ces) para resolver ecuaciones de grado superior a
cuatro. Además, si la ecuación no es polinómica no podemos resolverla más
que mediante métodos iterados que, incluso en el caso de las polinómicas de
grado no superior a cuatro, son más eficientes.
Definición 1.1 Una solución x de la ecuación f (x) = 0 se dice que tiene

multiplicidad n si
f (x) = f 0 (x) = f 00 (x) = · · · = f (n−1 (x) = 0 y f (n (x) 6= 0
Si la multiplicidad de una raı́z es 1, diremos que es simple .
Todos los métodos numéricos de resolución de ecuaciones presentan dificul-

tades cuando la ecuación tiene raı́ces múltiples, ya que todos ellos se basan en
los cambios de signo de la función y éstos son difı́cilmente detectables en un
entorno de una raı́z múltiple.
Ese hecho produce que en estos casos el problema esté mal condicionado.
1
2 Ecuaciones no lineales
1.1 Errores y condicionamiento en problemas

numéricos
Cualquier problema numérico se resuelve a través de un algoritmo que nos
proporciona unos resultados a partir de unos datos iniciales. Es decir, se trata
de realizar un proceso del tipo
Datos =⇒ Algoritmo =⇒ Resultados
Dado que cualquier algoritmo puede cometer errores, no sólo por el algoritmo
en sı́, sino porque los datos pueden venir afectados de algún tipo de error
(redondeo, etc.) es muy importante el estudio de los distintos tipos de error
que puedan cometerse con vista a la fiabilidad de los resultados.
Se denomina error absoluto de un número x que aproxima a otro x a la
distancia entre ellos. Ası́, por ejemplo, en el caso real vendrá dado por |x − x|,
es decir, por el valor absoluto de la diferencia entre ambos.
Obsérvese que si sólo disponemos del dato de que el error es, por ejemplo,
de 1m. no sabemos nada acerca de la fiabilidad del resultado, ya que no es lo
mismo decir que se ha cometido un error de un metro al medir la altura de
una persona que al medir la distancia entre dos galaxias.
Debemos reflejar de alguna manera “lo que se está evaluando” en el dato del
error. Para ello se utiliza el denominado error relativo que es el cociente entre
el error absoluto y el objeto evaluado, es decir, en el caso real
x−x
.
x
En el caso x = 0 sólo se utiliza el error absoluto. En la mayorı́a de los procesos
numéricos utilizaremos como error el error absoluto ya que lo que nos interesa
conocer es el número de cifras decimales exactas que posee.
Evidentemente cualquier algoritmo que trabaje con unos datos afectados de
algún tipo de error nos proporcionará unos resultados que también vendrán
afectados de errores. Estos errores pueden depender sólo de los datos iniciales
o también del proceso que se ha realizado.
Supongamos que, queremos evaluar f (x) y damos un dato aproximado x. Es
evidente que, en general, si x 6= x será f (x) 6= f (x).
Dado que f (x) − f (x) ' (x − x)f 0 (x), se tiene que
|f (x) − f (x)| ' |x − x| · |f 0 (x)|

Errores y condicionamiento en problemas numéricos 3
por lo que aunque el error del dato sea muy pequeño, si la derivada f 0 (x) es
muy grande, el resultado obtenido f (x) puede diferir mucho del valor exacto
f (x). Además el problema no está en el algoritmo que se aplique para evaluar
f (x) sino en el propio problema a resolver.
Diremos que un problema está mal condicionado si pequeños errores en los
datos producen grandes errores en los resultados.
Se trata entonces de definir algún número que nos indique el condiciona-
miento del problema. A éste número lo llamaremos número de condición del
problema y lo denotaremos por κ. En el ejemplo anterior es evidente que
κ(x) = |f 0 (x)|
Para el problema inverso, es decir conocido f (x) buscar el valor de x (resolver

una ecuación) se tendrá que
1
|x − x| ' |f (x) − f (x)|
|f 0 (x)|
por lo que si |f 0 (x)| fuese muy pequeño el problema estarı́a mal condicionado.
Ası́ pues, un problema (en ambos sentidos) estará mejor condicionado mien-
tras más se acerque a 1 su número de condición.
Respecto al algoritmo que se utiliza para resolver un determinado problema,
diremos que es inestable cuando los errores que se cometen en los diferentes
pasos del algoritmo hacen que el error total que se genera sea muy grande. Si,
por el contrario los errores que se producen en los distintos pasos no alteran
de forma significativa el resultado del problema, diremos que el algoritmo es
estable.
Obsérvese que si el algoritmo es inestable no va a generar un resultado fiable,
por lo que deberemos utilizar otro algoritmo. Sin embargo, por muy estable
que sea el algoritmo, si el problema está mal condicionado lo único que pode-
mos hacer es plantear un problema equivalente al nuestro pero con la seguridad
de que se trata de un problema bien condicionado.
Ası́, por ejemplo, si queremos calcular la tangente de 89◦ 590 lo primero que
debemos hacer es expresar el ángulo en radianes, por lo que necesariamente
debemos redondear el dato (el número π hay que redondearlo), por lo que
el dato vendrá afectado de un error de redondeo. Utilicemos el método que
utilicemos, dado que | tg x − tg x| ' |1 + tg 2 x| · |x − x| y tg π = +∞ el
problema estará mal condicionado. Sin embargo si tenemos en cuenta que
tg a + tg b
tg (a + b) = podemos reducir nuestro problema al cálculo de
1 − tg a tg b
la tangente de 44◦ 590 resultando éste un proceso bien condicionado, ya que

tg 45 = 1.
1.2 Método y algoritmo de la bisección: análisis

de errores
Este método consiste en la aplicación directa del teorema de Bolzano.
Teorema 1.1 [Teorema de Bolzano] Si f es una función continua en el

intervalo cerrado [a, b] y f (a) · f (b) < 0, existe un punto α ∈ (a, b) en el cual
f (α) = 0.
Nuestro problema se reduce a localizarla. Para ello, supongamos que está

separada, es decir, que en el intervalo [a, b] es la única raı́z que existe. Esto
podemos garantizarlo, por ejemplo, viendo que f 0 (x) 6= 0 en todo el intervalo,
ya que entonces, el Teorema de Rolle (que se enuncia a continuación) nos
garantiza la unicidad de la raı́z.
Teorema 1.2 [Teorema de Rolle] Si f (x) es una función continua en el

intervalo [a, b], derivable en (a, b) y f (a) = f (b), existe un punto α ∈ (a, b)
para el que f 0 (α) = 0.
En efecto, si f (x) tuviese dos raı́ces α1 y α2 en el intervalo [a, b], verificarı́a

las hipótesis del teorema de Rolle en el intervalo [α1 , α2 ] ⊂ [a, b], por lo que
deberı́a existir un punto α ∈ (α1 , α2 ) =⇒ α ∈ (a, b) en el que se anulara la
derivada, por lo que si f 0 (x) 6= 0 en todo el intervalo [a, b], no pueden existir
dos raı́ces de la ecuación en dicho intervalo.
Supongamos, sin pérdida de generalidad, que f es creciente en [a, b].
a+b b−a
a) Tomamos α0 = y ε= .
2 2
b) Si f (α0 ) = 0 entonces FIN. α = α0 es la raı́z exacta.
Si f (α0 ) > 0 entonces hacemos b = α0 .
Si f (α0 ) < 0 entonces hacemos a = α0 .
a+b b−a
Se repite el paso 1, es decir, hacemos α0 = y ε= .
2 2
Método y algoritmo de la bisección: análisis de errores 5
c) Si ε < 10−k (error prefijado), entonces FIN. El valor de α0 es la raı́z

buscada con k cifras decimales exactas.
Si ε > 10−k , entonces repetimos el paso 2.
El error cometido, tomando como raı́z de la ecuación el punto medio del inter-
b−a
valo obtenido en la en la iteración n-ésima, viene dado por εn = n+1 , por lo
2
1 −3
que si b − a = 1 y n = 9 se tiene que ε9 < 10 < 10 , es decir, en 9 iteraciones
2
obtenemos tres cifras decimales exactas.
1.2.1 Algoritmo
ai + b i
Para i = 0, 1, 2, . . . , n, . . ., Ii = [ai , bi ] y mi = (punto medio del
2
intervalo Ii ) con

 [ai , mi ] si sig(f (ai )) 6= sig(f (mi ))
I0 = [a, b] y Ii+1 =
 [m , b ] si sig(f (b )) 6= sig(f (m ))
i i i i
El proceso debe repetirse hasta que
f (mi ) = 0 o bien bi − ai < ε con ε > 0 prefijado.
Se tiene, por tanto:
Input: a, b, ε, f (x)
Output: m
while (b − a)/2 > ε

m ← a + (b − a)/2
if f (m) = 0
a←m
b←m
end if
if signf (a) = signf (m)
a←m
end if
if signf (b) = signf (m)
b←m
end if
end
print m
El hecho de calcular el punto medio de [a, b] como m = a + b-a/2 es debido a

que para valores muy pequeños de a y b puede darse el caso de que a+b/2 se
encuentre fuera del intervalo [a, b].
Ejemplo 1.1 Supongamos que se quiere calcular la raı́z cuadrada de 3, para lo

que vamos a buscar la raı́z positiva de la ecuación f (x) = 0 con f (x) = x2 − 3.
Dado que f (1) = −2 < 0 y f (2) = 1 > 0, el teorema de Bolzano nos garantiza
la existencia de una raı́z (que además sabemos que es única ya que f 0 (x) = 2x
no se anula en el intervalo [1, 2]).
Para obtener la raı́z con 14 cifras decimales exactas, es decir, con un error
menor que 10−14 tendrı́amos que detener el proceso cuando
2−1
< 10−14 =⇒ 2n+1 > 1014 =⇒ n ≥ 46
2n+1
es decir, tendrı́amos que detenernos en m46 para poder garantizar la precisión
exigida.
Una variante del método de la bisección es el método de la regula falsi o de

la falsa posición consistente en dividir el intervalo [a, b] en dos subintervalos
[a, c] ∪ [c, b] donde el punto c, a diferencia del punto medio m del método de
la bisección, es el punto de corte de la recta secante que pasa por los puntos
(a, f (a)) y (b, f (b)) con el eje de abscisas OX.
Y
(b, fs(b))
f (b)

(c,s0)

O
X

s

f (a)
(a, f (a))
Figura 1.1: Método de la régula f alsi.
La pendiente de dicha secante viene determinada por

f (b) − f (a) 0 − f (b)
m= =
b−a c−b
Según se utilicen los puntos (a, f (a)) y (b, f (b)) ó (c, 0) y (b, f (b)) respectiva-
mente.
Método y algoritmo de la bisección: análisis de errores 7
Despejando el valor de c obtenemos que

b−a
c = b − f (b) ·
f (b) − f (a)
pudiéndose dar los mismos casos que en el método de la bisección, es decir:
• Si f (c) = 0 la raı́z buscada es c.
• Si f (a) y f (c) tienen signos contrarios, la raı́z se encuentra en el intervalo

[a, c].
• Si son f (c) y f (b) los que tienen signos contrarios, la raı́z está en el
intervalo [c, b].
El algoritmo quedarı́a de la siguiente forma:
Input: a, b, ε, f (x)
Output: c
c←a
while absf (c) > ε
b−a
c ← b − f (b) ·
f (b) − f (a)
if f (c) = 0
a←c
b←c
end if
if signf (a) = signf (c)
a←c
end if
if signf (b) = signf (c)
b←c
end if
end
print c
Aplicando el algoritmo al Ejemplo 1.1 obtenemos la raı́z de 3, con 14 cifras

decimales exactas, en sólo 13 iteraciones frente a las 46 necesarias mediante el
método de la bisección.
En el método de la regula falsi la convergencia suele ser más rápida que en

el de la bisección aunque ambos métodos son de convergencia lenta siendo
necesario, por tanto, buscar otros métodos cuya convergencia sea más rápida.
Estos métodos están basados en el denominado Teorema del punto fijo que
estudiamos en la siguiente sección.
1.3 Punto fijo e iteración funcional

Ya se comentó que los métodos iterados consisten en crear una sucesión con-
vergente a la solución del problema.
Una función f : R → R se dice contractiva si verifica que
|f (x1 ) − f (x2 )| < |x1 − x2 | ∀ x 1 , x2 ∈ R
Si la función es derivable, basta comprobar que |f 0 (x)| ≤ q < 1 cualquiera

que sea el valor de x ∈ R para poder garantizar que se trata de una función
contractiva.
Si se desea resolver la ecuación f (x) = 0, se escribe esta de la forma x = ϕ(x),
donde ϕ(x) es una función contractiva, y partiendo de un determinado valor
inicial x0 , se construye la sucesión xn+1 = ϕ(xn ). La convergencia de esta
sucesión la garantiza el siguiente teorema.
Teorema 1.3 [Teorema del punto fijo] Dada la ecuación x = ϕ(x) en

la que |ϕ0 (x)| ≤ q < 1 cualquiera que sea x ∈ [a, b] con ϕ([a, b]) ⊆ [a, b] y un
punto x0 ∈ [a, b], la sucesión x0 , x1 , . . . , xn , . . . en la que xn+1 = ϕ(xn ) converge
a un valor x que es la única solución de la ecuación en dicho intervalo.
Demostración. Dado que xn+1 = ϕ(xn ) y xn = ϕ(xn−1 ) se tiene que
xn+1 − xn = ϕ(xn ) − ϕ(xn−1 ) = (xn − xn−1 )ϕ0 (c)
con c ∈ (xn−1 , xn ) (xn , xn−1 ) y, por tanto, en [a, b], por lo que
|xn+1 − xn | ≤ q |xn − xn−1 |
es decir
|x2 − x1 | ≤ q |x1 − x0 |
|x3 − x2 | ≤ q |x2 − x1 | ≤ q 2 |x1 − x0 |
...................................
|xn+1 − xn | ≤ q n |x1 − x0 |
Punto fijo e iteración funcional 9
Si construimos la serie
x0 + (x1 − x0 ) + (x2 − x1 ) + · · · + (xn+1 − xn ) + · · ·
observamos que la suma parcial Sn+1 = xn y que la serie es absolutamente con-

vergente por estar acotados sus términos, en valor absoluto, por los términos
de una serie geométrica de razón q < 1.
Sea x la suma de la serie. Entonces x = lim Sn+1 = lim xn , es decir, la
sucesión x0 , x1 , . . . , xn , . . . converge a x.
El valor obtenido x es solución de la ecuación, ya que por la continuidad de
la función ϕ(x) se verifica que
x = lim xn = lim ϕ(xn−1 ) = ϕ(lim xn−1 ) = ϕ(x)
Además, es la única solución de la ecuación en el intervalo [a, b], ya que de

existir otra x0 se tendrı́a que
x − x0 = ϕ(x) − ϕ(x0 ) = (x − x0 )ϕ0 (c) con c ∈ (x, x0 )
por lo que
(x − x0 )(1 − ϕ0 (c)) = 0
y dado que el segundo paréntesis es no nulo, se tiene que x = x0 .
En la Figura 1.2 puede observarse que el método converge si |ϕ0 (x)| ≤ q < 1,
mientras que si |ϕ0 (x)| > 1 el método es divergente.
En los casos (a) y (b), en los que |ϕ0 (x)| ≤ q < 1 el método converge
monótonamente en (a) y de forma oscilatoria o en espiral en (b).
En los casos (c) y (d), en los que |ϕ0 (x)| > 1 el método diverge de forma
monótona en (a) y de forma oscilatoria en (b).
1.3.1 Cota del error “a posteriori”

Si f (x) es una función continua en el intervalo [a, b] y derivable en el abierto
(a, b), sabemos por el Teorema del Valor Medio que existe un punto c ∈ (a, b)
f (b) − f (a)
tal que = f 0 (c).
b−a
Sea x una solución de la ecuación f (x) = 0 y sea xn una aproximación de ella
obtenida por un método iterado cualquiera. Supongamos f (x) continua en el
(a) (b)

x0 x0
(c) (d) ?
x0 x0
Figura 1.2: Esquema de la convergencia para el teorema del punto fijo.
intervalo cerrado [xn , x] ó [x, xn ] (dependiendo de que x sea mayor o menor

que xn ) y derivable en el abierto.
Existe entonces un punto c ∈ (xn , x) ó c ∈ (x, xn ) tal que
f (x) − f (xn )
= f 0 (c).
x − xn
f (xn )
Como f (x) = 0 y (x − xn ) = εn , nos queda que εn = − , obteniéndose:
f 0 (c)

|f (xn )| |f (xn )| |f (xn )| (x, xn )
|εn | = 0 ≤ 0 ≤ con ∈ (a, b)
|f (c)| n mı́n |f (x)| mı́n |f 0 (x)| (xn , x)
x∈
(x, xn ) x∈(a,b)
(xn , x)
Lo único que debemos exigir es que la derivada de la función no se anule en

ningún punto del intervalo (a, b).
Ejemplo 1.2 El cálculo de la raı́z cuadrada de 3 equivale al cálculo de la raı́z

positiva de la ecuación x2 = 3. Aunque más adelante veremos métodos cuya
convergencia es más rápida, vamos a realizar los siguientes cambios:
3+x
x2 = 3 =⇒ x + x2 = x + 3 =⇒ x(1 + x) = 3 + x =⇒ x =
1+x
Punto fijo e iteración funcional 11
Es decir, hemos escrito la ecuación de la forma x = ϕ(x) con

3+x
ϕ(x) =
1+x
Dado que sabemos que la raı́z de 3 está comprendida entre 1 y 2 y que
2 2 1
|ϕ0 (x)| = 2
≤ 2 = <1 para cualquier x ∈ [1, 2]
(1 + x) 2 2
podemos garantizar que partiendo de x0 = 1 el método convergerá a la raı́z
cuadrada de 3.
3 + xn
Ası́ pues, partiendo de x0 = 1 y haciendo xn+1 = obtenemos:
1 + xn
x1 = 2 x14 = 1.73205079844084
x2 = 1.66666666666667 x15 = 1.73205081001473
x3 = 1.75000000000000 x16 = 1.73205080691351
x4 = 1.72727272727273 x17 = 1.73205080774448
x5 = 1.73333333333333 x18 = 1.73205080752182
x6 = 1.73170731707317 x19 = 1.73205080758148
x7 = 1.73214285714286 x20 = 1.73205080756550
x8 = 1.73202614379085 x21 = 1.73205080756978
x9 = 1.73205741626794 x22 = 1.73205080756863
x10 = 1.73204903677758 x23 = 1.73205080756894
x11 = 1.73205128205128 x24 = 1.73205080756886
x12 = 1.73205068043172 x25 = 1.73205080756888
x13 = 1.73205084163518 x26 = 1.73205080756888
|f (xn )|
El error vendrá dado por εn < donde f (x) = x2 − 3, por lo que
mı́n |f 0 (x)|
x∈[1,2]
|x226 − 3|
ε26 < = 4.884981308350688 · 10−15 < 10−14
2
√
es decir, 3 = 1.73205080756888 con todas sus cifras decimales exactas.
Obsérvese que en el Ejemplo 1.1 vimos cómo eran necesarias 46 iteraciones

para calcular la raı́z cuadrada de 3 (con 14 cifras decimales exactas) mediante
el método de la bisección y 13 mediante el de la regula falsi, mientras que ahora
hemos necesitado 26. El hecho de que la convergencia de éste último método
sea más lenta que cuando se utiliza el método de la regula falsi estriba en la
mala elección de la función ϕ(x), por lo que vamos a ver cómo el método de
Newton-Raphson, generalmente conocido como método de Newton, nos permite
escribir la ecuación f (x) = 0 de la forma x = ϕ(x) de forma que la
convergencia sea muy rápida.
1.4 Análisis del método de Newton-Raphson

Si tratamos de resolver la ecuación f (x) = 0 y lo que obtenemos no es la
solución exacta x sino sólo una buena aproximación xn tal que x = xn + h
tendremos que
f (xn )
f (x) ' f (xn ) + h · f 0 (xn ) ⇒ h ' −
f 0 (xn )
por lo que
f (xn )
x ' xn −
f 0 (xn )
obteniéndose la denominada fórmula de Newton-Raphson
f (xn )
xn+1 = xn − (1.1)
f 0 (xn )
Si construimos, utilizando la fórmula de Newton-Raphson, la sucesión {xn } y

ésta converge, se tendrá que lim xn = x, ya que nos quedarı́a, aplicando lı́mites
en (1.1)
f (lim xn )
lim xn+1 = lim xn − ⇒ f (lim xn ) = 0
f 0 (lim xn )
siempre que f 0 (lim xn ) 6= 0, lo cual se verifica si exigimos que la función posea

una única raı́z en [a, b]. Dado que la raı́z de la ecuación en el intervalo [a, b] es
única, necesariamente lim xn = x.
Este método es también conocido como método de la tangente, ya que si
trazamos la tangente a la curva y = f (x) en el punto (xn , f (xn )) obtenemos
la recta y = f (xn ) + f 0 (xn )(x − xn ), que corta al eje y = 0 en el punto de
f (xn )
abscisa x = xn − 0 que es precisamente el valor de xn+1 de la fórmula de
f (xn )
Newton-Raphson.
En la Figura1.3 puede observarse cómo actúa geométricamente el método de

Newton-Raphson.
Lo más dificultoso del método consiste en el cálculo de la derivada de la
función ası́ como la obtención del valor inicial x0 .
Busquemos, a continuación, alguna cota del error.
Análisis del método de Newton-Raphson 13

x!
!
x2 x1 x0
Figura 1.3: Interpretación geométrica del método de Newton.

f (xn ) f (xn ) f (xn )
εn+1 = x − xn+1 = x − xn − 0 = (x − xn ) + 0
= εn + 0
f (xn ) f (xn ) f (xn )
Desarrollando f (x) en un entorno de xn se obtiene

(
f 00
(t) (x, xn ) si x < xn
0 = f (x) = f (xn + εn ) = f (xn )+f 0 (xn )εn+ ε2n con t ∈
2! (xn , x) si x > xn
Supuesto que f 0 (xn ) 6= 0 podemos dividir por dicha derivada para obtener
f (xn ) f 00 (t) 2 f 00 (t) 2

0= + ε n + ε = ε n+1 + ε
f 0 (xn ) 2f 0 (xn ) n 2f 0 (xn ) n
por lo que
|f 00 (t)| 2
|εn+1 | = ε ≤ k · ε2n (1.2)
2 |f 0 (xn )| n
|f 00 (x)|
donde k ≥ máx siendo [a, b] cualquier intervalo, en caso de existir,
x∈[a,b] 2|f 0 (x)|
que contenga a la solución x y a todas las aproximaciones xn .
Esta última desigualdad podemos (no queriendo precisar tanto) modificarla
para escribir
máx |f 00 (x)|
k≥ con x ∈ [a, b] y f 0 (x) 6= 0
2 mı́n |f 0 (x)|
Supuesta la existencia de dicho intervalo [a, b], el valor de k es independiente

de la iteración que se realiza, por lo que
n+1
k · |εn+1 | ≤ |k · εn |2 ≤ |k · εn−1 |4 ≤ · · · ≤ |k · ε0 |2
o lo que es lo mismo:
1 n
|εn | ≤ · |k · ε0 |2
k
donde es necesario saber acotar el valor de ε0 = x − x0 .
Es decir, si existe un intervalo [a, b] que contenga a la solución y a todas
las aproximaciones xn se puede determinar a priori una cota del error, o lo
que es lo mismo, se puede determinar el número de iteraciones necesarias para
obtener la solución con un determinado error.
1
Evidentemente, el proceso convergerá si |k · ε0 | < 1, es decir, si |ε0 | < . En
k
caso de ser convergente, la convergencia es de segundo orden como puede verse
en la ecuación (1.2).
1.4.1 Algoritmo
Una vez realizado un estudio previo para ver que se cumplen las condiciones
que requiere el método, establecer el valor inicial x0 y calcular el valor de
m = mı́n |f 0 (x)|, el algoritmo es el siguiente
x∈[a,b]
Input: a, b, x0 , ε, f (x), m
Output: x
x ← x0
e ← abs( f(x)/m)
while e > ε
f (x)
x←x−
f 0 (x)
e ← abs( f(x)/m)
end
print x
Ejemplo 1.3 En el Ejemplo 1.2 calculamos la raı́z de 3 con 14 cifras decimales

exactas en 26 iteraciones. Vamos a ver cómo se disminuye considerablemente
el número de iteraciones cuando se utiliza la fórmula de Newton-Raphson.
Partimos de la ecuación f (x) = x2 − 3 = 0, por lo que la fórmula de Newton-
Raphson nos dice que
1 3
xn+1 = xn +
2 xn
Dado que la raı́z de 3 es un número comprendido entre 1 y 2 y la función

0
f (x) = 2x no se anula en dicho intervalo, podemos aplicar el método de
Newton tomando como valor inicial x0 = 2. (Más adelante veremos porqué
debemos tomar 2 como valor inicial), obteniéndose:
x0 =2
x1 = 1.75000000000000
x2 = 1.73214285714286
x3 = 1.73205081001473
x4 = 1.73205080756888
|f (xn )|
El error vendrá dado, al igual que en el Ejercicio 1.2 por εn < ,
mı́n |f 0 (x)|
x∈[1,2]
por lo que
|x24 − 3|
ε4 < = 4.884981308350688 · 10−15 < 10−14
2
√
es decir, 3 = 1.73205080756888 con todas sus cifras decimales exactas.

1 A
Nota: La fórmula xn+1 = xn + es conocida como fórmula de Heron
2 xn
ya que este matemático la utilizaba para aproximar la raı́z cuadrada de un
número real positivo A hacia el año 100 a.C. pues sabı́a que si tomaba un
valor inicial x0 y éste fuese mayor que la raı́z buscada, necesariamente A/x0
serı́a menor que su raı́z y viceversa, por lo que la media entre ambos debı́a ser
una mejor aproximación que x0 . En otras palabras, en cada iteración tomaba
como aproximación de la raı́z la media aritmética entre el valor xn anterior y
el cociente A/xn .
Puede observarse cómo la convergencia del método de Newton-Raphson es

mucho más rápida que la del método de la bisección y que el de la regula falsi,
ya que sólo hemos necesitado 5 iteraciones frente a las 46 que se necesitan en
el método de la bisección o las 13 del método de la regula falsi.
1.4.2 Análisis de la convergencia: Regla de Fourier

Hay que tener en cuenta que la naturaleza de la función puede originar difi-
cultades, llegando incluso a hacer que el método no converja.
Ejemplo 1.4 Tratemos de determinar, por el método de Newton-Raphson, la

raı́z positiva de la función f (x) = x10 −1, tomando como valor inicial x0 = 0.5.
La fórmula de Newton-Raphson es, en este caso,
x10
n −1
xn+1 = xn − .
10x9n
Aplicando el algoritmo se obtienen los valores
x1 = 51.65 x20 = 6.97714912329906
x2 = 46.485 x30 = 2.43280139954230
x3 = 41.8365 x40 = 1.00231602417741
x4 = 37.65285 x41 = 1.00002393429084
x5 = 33.887565 x42 = 1.00000000257760
x10 = 20.01026825685012 x43 =1
Puede observarse que la convergencia es muy lenta y sólo se acelera (a partir
de x40 ) cuando estamos muy cerca de la raı́z buscada.
a) Si en las proximidades de la raı́z existe un punto de inflexión, las itera-

ciones divergen progresivamente de la raı́z.

9

r
x
0

b) El método de Newton-Raphson oscila en los alrededores de un máximo

o un mı́nimo local, persistiendo o llegando a encontrarse con pendientes
cercanas a cero, en cuyo caso la solución se aleja del área de interés.
b
b
r
b
b
x0
c) Un valor inicial cercano a una raı́z puede converger a otra raı́z muy
distante de la anterior como consecuencia de encontrarse pendientes
cercanas a cero. Una pendiente nula provoca una división por cero
(geométricamente, una tangente horizontal que jamás corta al eje de
abscisas).
#
# Z -
@ #
I # r r Z
Z
@ x0 x0
Estos problemas pueden detectarse previamente a la aplicación del método.

Supongamos que tenemos acotada, en el intervalo [a, b], una única raı́z x de
la ecuación f (x) = 0 y que f 0 (x) y f 00 (x) no se anulan en ningún punto del
intervalo [a, b], es decir, que ambas derivadas tienen signo constante en dicho
intervalo.
Obsérvese que si f (a)f (b) < 0, dado que f 0 (x) no se anula en el intervalo (a, b)
sabemos, por los teoremas de Bolzano y Rolle, que existe una única raı́z en
dicho intervalo. Además, por las condiciones exigidas sabemos que no existe,
en (a, b) ningún punto crı́tico (ni extremo relativo ni punto de inflexión), con
lo que habremos evitado los problemas expuestos anteriormente.
En cualquiera de los cuatro casos posibles (véase la Figura 1.4), la función
cambia de signo en los extremos del intervalo (debido a que la primera derivada
no se anula en dicho intervalo), es decir, dado que la segunda derivada tiene
signo constante en [a, b], en uno de los dos extremos la función tiene el mismo
signo que su segunda derivada.
En estos casos, el método de Newton es convergente debiéndose tomar como
valor inicial 
 a si f (a) · f 00 (a) > 0
x0 =
 b si f (b) · f 00 (b) > 0
es decir, el extremo en el que la función tiene el mismo signo que su derivada

segunda.
Este resultado, que formalizamos a continuación en forma de teorema es
conocido como regla de Fourier.
Teorema 1.4 [Regla de Fourier] Sea f (x) una función continua y dos
veces derivable [a, b]. Si sig f (a) 6= sig f (b) y sus dos primeras derivadas f 0 (x)
y f 00 (x) no se anulan en [a, b] existe una única raı́z de la ecuación f (x) = 0 en
dicho intervalo y se puede garantizar la convergencia del método de Newton-
Raphson tomando como valor inicial x0 el extremo del intervalo en el que la
función y su segunda derivada tienen el mismo signo.
a a b b
b b a a
f 0 (x) > 0 f 0 (x) > 0 f 0 (x) < 0 f 0 (x) < 0

f 00 (x) < 0 f 00 (x) > 0 f 00 (x) > 0 f 00 (x) < 0
x0 = a x0 = b x0 = a x0 = b
Figura 1.4: Los cuatro casos posibles
Gracias a que la convergencia es de segundo orden, es posible modificar el

método de Newton para resolver ecuaciones que poseen raı́ces múltiples.
1.4.3 Método de Newton para raı́ces múltiples

Cuando el método de Newton converge lentamente nos encontramos con una
raı́z múltiple y, a diferencia de lo que ocurrı́a con otros métodos, podemos
modificar el método para acelerar la convergencia.
Sea x una raı́z de multiplicidad k de la ecuación f (x) = 0. En este caso,
el método de Newton converge muy lentamente y con grandes irregularidades
debido al mal condicionamiento del problema.
f (xn )
Si en vez de hacer xn+1 = xn − 0 hacemos
f (xn )
f (xn )
xn+1 = xn − k
f 0 (xn )
donde k representa el orden de la primera derivada que no se anula para x = x
(multiplicidad de la raı́z x), el método sigue siendo de segundo orden.
En la práctica, el problema es que no conocemos k pero a ello nos ayuda la
rapidez del método.
Ejemplo 1.5 Para resolver la ecuación x − sen x = 0 comenzamos por ex-

presarla de la forma x = sen x, por lo que las soluciones serán los puntos de
intersección de la curva y = sen x con y = x (Fig.1.5).
Aunque es conocido que la solución de la ecuación es x = 0, supondremos que
sólo conocemos que está comprendida entre -1 y 1 y vamos aplicar el método
6
y=x
y = sen x
-1
-
0 1
Figura 1.5: Las gráficas de y = x e y = sen x.
de Newton.
xn − sen xn sen xn − xn cos xn
xn+1 = xn − =
1 − cos xn 1 − cos xn
Comenzando con x0 = 1 se obtiene:
x0 = 1
.....................
 0
 f (x10 ) = 00 0001 . . .

x10 = 00 016822799 . . . f 00 (x10 ) = 00 016 . . .

 000
f (x10 ) = 00 9998 . . .
.....................
 0
 f (x20 ) = 00 00000001 . . .

x20 = 00 0000194 . . . f 00 (x20 ) = 00 0019 . . .

 000
f (x20 ) = 00 9999 . . .
Como la convergencia es muy lenta, hace pensar que se trata de una raı́z
múltiple. Además, como la primera y la segunda derivadas tienden a cero y
la tercera lo hace a 1, parece que nos encontramos ante una raı́z triple, por lo
que aplicamos el método generalizado de Newton.
xn − sen xn
xn+1 = xn − 3
1 − cos xn
que comenzando, al igual que antes, por x0 = 1 se obtiene:
x0 =1
x1 = −00 034 . . .
x2 = 00 000001376 . . .
x3 = 00 00000000000009 . . .
que se ve que converge rápidamente a x = 0.

Aplicamos ahora la cota del error a posteriori a este valor y obtenemos:
x = 0 =⇒ f (x) = x − sen x = 0 =⇒ la solución es exacta.
f 0 (x) = 1 − cos x =⇒ f 0 (x) = 0 




00 00
f (x) = sen x =⇒ f (x) = 0 =⇒ se trata de una raı́z triple.



f 000 (x) = cos x =⇒ f 000 (x) = 1


1.5 Un problema mal condicionado: ceros de

un polinomio
Si queremos resolver la ecuación P (x) = 0 el problema viene condicionado
por la función 1/P 0 (x) (su número de condición), por lo que si la derivada es
muy pequeña el problema estará mal condicionado, pero si la derivada es muy
grande cualquier método se hace muy lento, por lo que lo ideal serı́a que la
derivada estuviese muy próxima a 1, pero claro está, ese derivada ha de estar
muy próxima a 1 en todas las raı́ces del polinomio, cosa que es estadı́sticamente
casi imposible. Por tanto, el cálculo de los ceros de un polinomio es un ejemplo
de problema mal condicionado. Sin embargo, vamos a estudiar cómo podemos
aproximar un determinado cero del polinomio.
Hemos visto que uno de los mayores problemas que presenta la resolución de
una ecuación es encontrarnos que posee raı́ces múltiples ya que, en un entorno
de ellas, o bien la función no cambia de signo, o bien se aproxima demasiado
a cero y, por tanto, cualquier método puede dar soluciones erróneas.
En el caso de las ecuaciones algebraicas (Pn (x) = 0) este problema puede
solventarse buscando otra ecuación que posea las mismas raı́ces que la dada
pero todas ellas simples, es decir, eliminando las raı́ces múltiples.
Por el teorema fundamental del Álgebra sabemos que Pn (x) posee n raı́ces y,
por tanto, puede ser factorizado de la forma
Pn (x) = a0 (x − x1 )(x − x2 ) · · · (x − xn )
donde {x1 , x2 , . . . , xn } son los ceros del polinomio.

Si existen raı́ces múltiples, las podemos agrupar para obtener:
Pn (x) = a0 (x − x1 )m1 (x − x2 )m2 · · · (x − xk )mk

Un problema mal condicionado: ceros de un polinomio 21
donde mi (i = 1, 2, . . . k) representa la multiplicidad de la raı́z xi y verificándo-

se que m1 + m2 + · · · mk = n
Derivando esta expresión obtenemos:
P 0 (x) = na0 (x − x1 )m1 −1 · · · (x − xk )mk −1 Qk−1 (x)
con Qk−1 (xi ) 6= 0 i = 1, 2, . . . , k
Por tanto, si x es una raı́z de la ecuación P (x) = 0 con multiplicidad k, es
también una raı́z de P 0 (x) = 0 pero con multiplicidad k − 1.
D(x) = mcd [P (x), P 0 (x)] = (x − x1 )m1 −1 · · · (x − xk )mk −1
por lo que
P (x)
Q(x) = = a0 (x − x1 )(x − x2 ) · · · (x − xk )
D(x)
Es decir, hemos encontrado un polinomio cuyas raı́ces son las mismas que las
de P (x) pero todas ellas simples.
Si ya conocemos que una ecuación sólo tiene raı́ces simples y queremos encon-
trarlas, parece apropiado que un primer paso consista en detectar las posibles
situaciones en éstas. Ası́ por ejemplo, si son reales, determinar intervalos de
una amplitud reducida en los que se encuentren las raı́ces de la ecuación.
Definición 1.2 Dada una ecuación f (x) = 0

(en general compleja) se denomina acotar las
raı́ces a buscar dos números reales positivos r
y R tales que r ≤ |x| ≤ R para cualquier raı́z
x de la ecuación.
Geométricamente consiste en determinar una

corona circular de radios r y R dentro de la
cual se encuentran todas las raı́ces.
En el caso real se reduce a los intervalos
(−R, −r) y (r, R).
Veamos, a continuación, una cota para las raı́ces de una ecuación algebraica.
Proposicin 1.1 Si x es una raı́z de la ecuación

P (x) = a0 xn + a1 xn−1 + · · · + an−1 x + an = 0
se verifica que:
A
|x| < 1 + siendo A = máx |ai |
|a0 | i≥1
Demostración. Sea P (x) = a0 xn + a1 xn−1 + · · · + an−1 x + an . Tomando

módulos tenemos
|P (x)| = |a0 xn + a1 xn−1 + · · · + an−1 x + an | ≥
≥ |a0 xn | − |a1 xn−1 + · · · + an−1 x + an | ≥
h i
≥ |a0 xn | − |a1 xn−1 | + · · · + |an−1 x| + |an | =
h i
= |a0 xn | − |a1 | |x|n−1 + · · · + |an−1 | |x| + |an | ≥
≥ |a0 xn | − A |x|n−1 + · · · + |x| + 1

(Para considerar el último paréntesis como una progresión geométrica habrı́a

que añadir los términos que, probablemente, falten en P (x) y suponer que,
además, es |x| =
6 1).
n |x|n − 1
|P (x)| ≥ |a0 | |x| − A
|x| − 1
Dado que el teorema es trivial para |x| < 1, supondremos que |x| > 1 y
entonces:
|x|n

n n A
|P (x)| > |a0 | |x| − A = |x| |a0 | −
|x| − 1 |x| − 1
Como la expresión anterior es cierta para cualquier |x| > 1, sea |x| > 1 con
P (x) = 0. Entonces

n A A
0 > |x| |a0 | − =⇒ |a0 | − < 0 =⇒
|x| − 1 |x| − 1
A A
|a0 | < =⇒ |x| − 1 < =⇒
|x| − 1 |a0 |
A
|x| < 1 + con |x| > 1
|a0 |
Es evidente que esta cota también la verifican las raı́ces x con |x| < 1.
Proposicin 1.1 [Regla de Laguerre] Consideremos la ecuación
P (x) = a0 xn + a1 xn−1 + · · · + an−1 x + an = 0
Sean C(x) = b0 xn−1 + · · · + bn−2 x + bn−1 el cociente y r el resto de la división

de P (x) entre x − c. Si r ≥ 0 y bi ≥ 0 para 0 ≤ i ≤ n − 1, el número real c es
una cota superior para las raı́ces positivas de la ecuación. (Trivialmente lo es
también para las raı́ces negativas).
El procedimiento consiste en comenzar con la cota obtenida anteriormente

(que no suelen ser muy buena) e ir disminuyéndola hasta afinarla todo lo que
podamos.
Las cotas obtenidas anteriormente nos delimitan la zona en la que debemos
estudiar la existencia de soluciones de la ecuación pero, en realidad, lo que más
nos acerca a nuestro problema (resolver la ecuación) es separar cada raı́z en un
intervalo. A este proceso se le conoce como separación de raı́ces y estudiaremos
un método que se conoce como método de Sturm que nos permite separar las
raı́ces de una ecuación.
1.5.1 Sucesiones de Sturm

Una sucesión de Sturm en [a, b] es un conjunto de funciones continuas en dicho
intervalo f0 (x), f1 (x), . . . , fn (x) que cumplen las siguientes condiciones:
• fn (x) 6= 0 cualquiera que sea x ∈ [a, b]. Es decir, el signo de fn (x)

permanece constante en el intervalo [a, b].
• Las funciones fi (x) y fi+1 (x) no se anulan simultáneamente. En otras

palabras, si fi (c) = 0 entonces fi−1 (c) 6= 0 y fi+1 (c) 6= 0.
• Si fi (c) = 0 entonces fi−1 (c) y fi+1 (c) tienen signos opuestos, es decir,
fi−1 (c) · fi+1 (c) < 0. (Engloba al apartado anterior).
f0 (x)
• Si f0 (c) = 0 con c ∈ [a, b] entonces pasa de negativa a positiva
f1 (x)
en c. (Está bien definida en c por ser f1 (c) 6= 0 y es creciente en dicho
punto).
Teorema 1.5 [Teorema de Sturm]. Sea f0 (x), f1 (x), . . . , fn (x) una sucesión
de Sturm en el intervalo [a, b] y consideremos las sucesiones
sig[f0 (a)] sig[f1 (a)] · · · sig[fn (a)]
sig[f0 (b)] sig[f1 (b)] · · · sig[fn (b)]
teniendo en cuenta que si alguna de las funciones se anula en uno de los
extremos del intervalo [a, b] pondremos en su lugar, indistintamente, signo +
o - y denotemos por N1 al número de cambios de signo de la primera sucesión
y por N2 al de la segunda (siempre es N1 ≥ N2 ).
En estas condiciones, el número de raı́ces existentes en el intervalo [a, b] de
la ecuación f0 (x) = 0 viene dado por N1 − N2 .
La construcción de una sucesión de Sturm es, en general, complicada. Sin

embargo, cuando se trabaja con funciones polinómicas, el problema es mucho
más simple además de que siempre es posible construir una sucesión de Sturm
válida para cualquier intervalo.
Dada la ecuación algebraica Pn (x) = 0, partimos de
f0 (x) = Pn (x) y f1 (x) = Pn0 (x)
Para determinar las demás funciones de la sucesión vamos dividiendo fi−1 (x)
entre fi (x) para obtener
fi−1 (x) = ci (x) · fi (c) + ri (x)
donde ri (x) tiene grado inferior al de fi (x) y hacemos
fi+1 (x) = −ri (x)
Como el grado de fi (x) va decreciendo, el proceso es finito. Si se llega a un

resto nulo (el proceso que estamos realizando es precisamente el del algoritmo
de Euclides) la ecuación posee raı́ces múltiples y se obtiene el máximo común
divisor D(x) de Pn (x) y Pn0 (x). Dividiendo Pn (x) entre D(x) obtenemos una
nueva ecuación que sólo posee raı́ces simples. La sucesión fi (x)/D(x) es una
sucesión de Sturm para la ecuación P (x)/D(x) = Q(x) = 0 que posee las
mismas raı́ces que P (x) = 0 pero todas simples.
Si llegamos a un resto constante, no nulo, es éste quien nos determina la
finalización del proceso. Hemos obtenido, de esta manera, una sucesión de
Sturm válida para cualquier intervalo [a, b].
Nota: Obsérvese que, al igual que en el algoritmo de Euclides, podemos ir

multiplicando los resultados parciales de las divisiones por cualquier constante
positiva no nula, ya que sólo nos interesa el resto (salvo constantes positivas)
de la división.
Ejemplo 1.6 Vamos a construir una sucesión de Sturm que nos permita se-
parar las raı́ces de la ecuación x4 + 2x3 − 3x2 − 4x − 1 = 0.
f0 (x) = x4 + 2x3 − 3x2 − 4x − 1. f00 (x) = 4x3 + 6x2 − 6x − 4.
f1 (x) = 2x3 + 3x2 − 3x − 2.
2x4 + 4x3 − 6x2 − 8x − 2 |2x3 + 3x2 − 3x − 2

− 2x4 − 3x3 + 3x2 + 2x x+1
x3 − 3x2 − 6x − 2 multiplicando por 2
2x3 − 6x2 − 12x − 4
−2x3 − 3x2 + 3x + 2
−9x2 − 9x − 2
f2 (x) = 9x2 + 9x + 2.
18x3 + 27x2 − 27x − 18 |9x2 + 9x + 2

− 18x3 − 18x2 − 4x 2x + 1
9x2 − 31x − 18
− 9x2 − 9x − 2
−40x − 20
f3 (x) = 2x + 1.
18x2 + 18x + 4 |2x + 1

− 18x2 − 9x 9x + 9
9x + 4 multiplicando por 2
18x + 8
− 18x − 9
−1
f4 (x) = 1.
−∞ −3 −2 −1 0 1 2 +∞
f0 (x) = x4 + 2x3 − 3x2 − 4x − 1 + + − − − − + +
f1 (x) = 2x3 + 3x2 − 3x − 2 − − ± + − ± + +
f2 (x) = 9x2 + 9x + 2 + + + + + + + +
f3 (x) = 2x + 1 − − − − + + + +
f4 (x) = 1 + + + + + + + +
cambios de signo 4 4 3 3 1 1 0 0
Sabemos, por ello, que existe una raı́z en el intervalo (−3, −2), dos raı́ces en
el intervalo (−1, 0) y una cuarta raı́z en el intervalo (1, 2).
Como f0 (−1) = −1 < 0, f0 (−00 5) = 00 0625 > 0 y f0 (0) = −1 < 0 podemos

separar las raı́ces existentes en el intervalo (−1, 0) y decir que las cuatro raı́ces
de la ecuación dada se encuentran en los intervalos
(−3 − 2) (−1, −00 5) (−00 5, 0) y (1, 2)
Si, una vez eliminadas las raı́ces múltiples y separadas las raı́ces, queremos
resolver la ecuación, utilizaremos (excepto en casos muy determinados como el
del Ejemplo 1.4) el método de Newton-Raphson. Al aplicarlo nos encontramos
con que tenemos que calcular, en cada paso, los valores de P (xn ) y P 0 (xn ) por
lo que vamos a ver, a continuación, un algoritmo denominado algoritmo de
Horner que permite realizar dichos cálculos en tiempo lineal.
1.5.2 Algoritmo de Horner

Supongamos un polinomio P (x) y un número real (en general también puede
ser complejo) x0 ∈ R. Si dividimos P (x) entre x − x0 sabemos que el resto es
un polinomio de grado cero, es decir, un número real, por lo que

 r∈R
P (x) = (x − x0 )Q(x) + r con y
gr[Q(x)] = gr[P (x)] − 1

Haciendo x = x0 obtenemos que
P (x0 ) = 0 · Q(x0 ) + r =⇒ P (x0 ) = r
Este resultado es conocido como teorema del resto y lo enunciamos a conti-

nuación.
Teorema 1.6 [Teorema del resto] El valor numérico de un polinomio

P (x) para x = x0 viene dado por el resto de la división de P (x) entre x − x0 .
Sea
P (x) = a0 xn + a1 xn−1 + · · · + an−1 x + an .
Si llamamos bi (0 ≤ i ≤ n − 1) a los coeficientes del polinomio cociente
P (x) − P (x0 )
= Q(x) = b0 xn−1 + b1 xn−2 + · · · + bn−2 x + bn−1
x − x0
se tiene que
b 0 = a0
b 1 = a1 + x 0 b 0
b 2 = a2 + x 0 b 1
..
.
bn−1 = an−1 + x0 bn−2
r = P (x0 ) = an + x0 bn−1
Este procedimiento para calcular el polinomio cociente Q(x) y el valor numé-

rico de P (x0 ) es conocido como algoritmo de Horner .
Una regla útil para hacer los cálculos a mano es la conocida regla de Ruffini
que consiste en disponer las operaciones como se indica a continuación.
a0 a1 a2 ··· an−1 an
x0 x0 b0 x0 b1 · · · x0 bn−2 x0 bn−1
b0 b1 b2 ··· bn−1 P (x0 )
Además, dado que
P (x) = Q(x)(x − x0 ) + P (x0 ) =⇒ P 0 (x) = Q0 (x)(x − x0 ) + Q(x)
se tiene que
P 0 (x0 ) = Q(x0 )
y el cálculo de Q(x0 ) es análogo al que hemos realizado para P (x0 ), es decir,
aplicando el algoritmo de Horner a Q(x) obtenemos
Q(x) = C(x)(x − x0 ) + Q(x0 ) donde Q(x0 ) = P 0 (x0 ).
Si utilizamos la regla de Ruffini sólo tenemos que volver a dividir por x0 como
se muestra a continuación.
a0 a1 a2 ··· an−2 an−1 an

x0 x0 b0 x0 b1 · · · x0 bn−3 x0 bn−2 x0 bn−1
b0 b1 b2 ··· bn−2 bn−1 P (x0 )
x0 x0 c0 x0 c1 · · · x0 cn−3 x0 cn−2
c0 c1 c2 ··· cn−2 P 0 (x0 )
Ejemplo 1.7 Consideremos el polinomio P (x) = 2x4 + x3 − 3x2 + 4x − 5 y

vamos a calcular los valores de P(2) y P’(2). Aplicando reiteradamente al regla
de Ruffini obtenemos
2 1 −3 4 −5
2 4 10 14 36
2 5 7 18 31
2 4 18 50
2 9 25 68
por lo que
P (2) = 31 y P 0 (2) = 68
Evidentemente, la regla de Ruffini nos es útil para realizar cálculos a mano

con una cierta facilidad, pero cuando los coeficientes de P (x) y el valor de x0
no son enteros sino que estamos trabajando con varias cifras decimales, deja
de ser efectivo y debemos recurrir al algoritmo de Horner en una máquina.
1.6 Sistemas de ecuaciones no lineales

Dado un sistema de ecuaciones no lineales
f1 (x1 , x2 , . . . , xm ) = 0
f2 (x1 , x2 , . . . , xm ) = 0
..
.
fm (x1 , x2 , . . . , xm ) = 0
podemos expresarlo de la forma f (x) = 0 donde x es un vector de Rm y f una

función vectorial de variable vectorial, es decir:
f : D ⊂ Rm → Rm
o lo que es lo mismo, f = (f1 , f2 , . . . , fm ) con fi : Rm → R para 1 ≤ i ≤ m.

Ası́, por ejemplo, el sistema
)
x2 − 2x − y + 00 5 = 0
(1.3)
x2 + 4y 2 + 4 = 0
Sistemas de ecuaciones no lineales 29
puede considerarse como la ecuación f (x) = 0 (obsérvese que 0 representa

ahora al vector nulo, es decir, que 0 = (0, 0)T ) con x = (x, y)T y f = (f1 , f2 )
siendo 
 f1 (x) = x2 − 2x − y + 00 5
 f (x) = x2 + 4y 2 + 4
2
Como hemos transformado nuestro sistema en una ecuación del tipo f (x) = 0,
parece lógico que tratemos de resolverla por algún método paralelo a los estu-
diados para ecuaciones no lineales como puedan ser la utilización del teorema
del punto fijo o el método de Newton.
Si buscamos un método iterado basado en el teorema del punto fijo, debe-
mos escribir la ecuación f (x) = 0 de la forma x = F (x) (proceso que puede
realizarse de muchas formas, la más sencilla es hacer F (x) = x + f (x)) para,
partiendo de un vector inicial x0 construir la sucesión de vectores


 x1n+1 = F1 (x1n , x2n , . . . , xm
n)



 x2 = F2 (x1 , x2 , . . . , xm )

n+1 n n n
xn+1 = F (xn )
..
.






 m
xn+1 = Fm (x1n , x2n , . . . , xm n)
En el ejemplo (1.3) podemos hacer

x2 − y + 00 5
x2 − 2x − y + 00 5 = 0 =⇒ 2x = x2 − y + 00 5 =⇒ x =
2
x2 + 4y 2 − 4 = 0 =⇒ x2 + 4y 2 + y − 4 = y =⇒ y = x2 + 4y 2 + y − 4
es decir,


 x = (x, y)T



x = F (x) con y
2 0

 F (x) = ( x − y + 0 5 , x2 + 4y 2 + y − 4)T



2
Si x es una solución de la ecuación y xn+1 es una aproximación obtenida, se
tiene que
kx − xn+1 k = kF (x) − F (xn )k = kF 0 (α)(x − xn )k ≤ kF 0 (α)k · kx − xn k
por lo que si kF 0 (x)k ≤ q < 1 para cualquier punto de un determinado entorno
de la solución, se tiene que
kx − xn+1 k ≤ kx − xn k
y la sucesión
x0 , x1 , . . . , xn , . . .
converge a la única raı́z de la ecuación x = F (x) en la bola considerada
(intervalo de Rm ).
Es importante observar que:
a) Se ha utilizado el concepto de norma vectorial al hacer uso de kx − xn )k.
b) Se ha utilizado una versión del teorema del valor medio para varias va-
riables al decir que
kF (x) − F (xn )k ≤ kF 0 (α)kk(x − xn )k
c) Se ha utilizado el concepto de norma matricial al hacer uso de kF 0 (α)k

ya que F 0 (x) es la matriz jacobiana de la función F , es decir
∂F1 ∂F1
 
···
 ∂x1 ∂xm 
 
0
 .. . . .
. 
F (x) =  . . .  
 
 ∂F ∂F 
m m
···
∂x1 ∂xm
d) Se supone que kF 0 (α)(x − xn )k ≤ kF 0 (α)k · kx − xn k o de forma más

general, que para cualquier matriz A (cuadrada de orden n) y cualquier
vector de Rn se verifica que
kAxk ≤ kAk · kxk
e) Que el teorema del punto fijo es generalizable a funciones vectoriales de

variable vectorial.
1.6.1 Método de Newton

Consideremos la ecuación f (x) = 0 (donde f es una función vectorial de
variable vectorial) equivalente a un sistema de ecuaciones no lineales.
Sea x la solución exacta del sistema y xn una aproximación de ella. Si
denotamos por h = x − xn se tiene que
x = xn + h
Sistemas de ecuaciones no lineales 31
y haciendo uso del desarrollo de Taylor obtenemos que
0 = f (x) = f (xn + h) ' f (xn ) + f 0 (xn )h
de donde
h ' −f 0−1 (xn )f (xn )
y, por tanto
x ' xn − f 0−1 (xn )f (xn ).
Esta aproximación es que utilizaremos como valor de xn+1 , es decir
xn+1 = xn − f 0−1 (xn )f (xn )
Obsérvese entonces que cada iteración del método de Newton se reduce al

cálculo del vector h correspondiente y éste no es más que la solución del sistema
de ecuaciones lineales
f 0 (xn )h + f (xn ) = 0
En el ejemplo (1.3) se tiene que f (x) = 0 con x = (x, y)T y f = (f1 , f2 )T
donde
f1 (x) = x2 − 2x − y + 00 5 y f2 (x) = x2 + 4y 2 − 4
f1 (x)

t 1

t
-1 1
-2 2
-0’5
f2 (x)
-1
Tomando como valor inicial x0 = (−00 5, 1)T se tiene que
f (x0 ) = (00 75, 00 25)T

   
2x − 2 −1 −3 −1
J(x) =   =⇒ f 0 (x0 ) = J(x0 ) =  
2x 8y −1 8
por lo que debemos resolver el sistema

    
−3 −1 h11 −00 75
  = 
2 0
−1 8 h1 −0 25
   
h11 00 25
cuya solución es h1 =  =  y, por tanto
h21 0
 
−00 25
x1 = x0 + h1 =  
1
Aplicando de nuevo el método se obtiene

   
00 0625 −00 25 −1
f (x1 ) =   f 0 (x1 ) = J(x1 ) =  
0 0
0 0625 −0 5 8
obteniéndose el sistema
    
−00 25 −1 h12 −00 0625
  = 
0 2 0
−0 5 8 h2 −0 0625
   
h12 00 022561 . . .
cuya solución es h2 =  =  y, por tanto
0
h22 −0 006 . . .
 
−00 227439 . . .
x2 = x1 + h2 =  
00 994 . . .
En definitiva, podemos observar que la resolución de un sistema de ecuaciones

no lineales mediante el método de Newton se reduce, en cada iteración, a la
resolución de un sistema de ecuaciones lineales por lo que el tema siguiente lo
dedicaremos al estudio de dichos sistemas de ecuaciones.
Ejercicios propuestos 33
1.7 Ejercicios propuestos

Ejercicio 1.1 Dada la ecuación xex − 1 = 0, se pide:
a) Estudiar gráficamente sus raı́ces reales y acotarlas.
b) Aplicar el método de la bisección y acotar el error después de siete ite-

raciones.
c) Aplicar el método de Newton, hasta obtener tres cifras decimales exactas.
Ejercicio 1.2 Probar que la ecuación x2 + ln x = 0 sólo tiene una raı́z real y
hallarla, por el método de Newton, con 6 cifras decimales exactas.
Ejercicio 1.3 Eliminar las raı́ces múltiples en la ecuación x6 − 2x5 + 3x4 −

4x3 + 3x2 − 2x + 1 = 0. Resolver, exactamente, la ecuación resultante y
comprobar la multiplicidad de cada raı́z en la ecuación original.
Ejercicio 1.4 Dada la ecuación 8x3 − 4x2 − 18x + 9 = 0, acotar y separar sus
raı́ces reales.
Ejercicio 1.5 Dado el polinomio P (x) = x3 + 3x2 + 2 se pide:
a) Acotar sus raı́ces reales.
b) Probar, mediante una sucesión de Sturm, que P (x) sólo posee una raı́z
real y determinar un intervalo de amplitud 1 que la contenga.
c) ¿Se verifican, en dicho intervalo, las hipótesis del teorema de Fourier?

En caso afirmativo, determinar el extremo que debe tomarse como valor
inicial x0 para garantizar la convergencia del método de Newton.
d) Sabiendo que en un determinado momento del proceso de Newton se ha

obtenido xn = −3.1958, calcular el valor de xn+1 ası́ como una cota del
error en dicha iteración.
Ejercicio 1.6 Aplicar el método de Sturm para separar las raı́ces de la ecuación
2x6 − 6x5 + x4 + 8x3 − x2 − 4x − 1 = 0
y obtener la mayor de ellas con seis cifras decimales exactas por el método de
Newton.
Ejercicio 1.7 Se considera el polinomio P (x) = x3 − 6x2 − 3x + 7.
a) Probar, mediante una sucesión de Sturm, que posee una única raı́z en el
intervalo (6, 7).
b) Si expresamos la ecuación P (x) = 0 de la forma x = F (x) = 13 (x3 −

6x2 + 7), ¿podemos asegurar su convergencia?
c) Probar, aplicando el criterio de Fourier, que tomando como valor inicial

x0 = 7, el método de Newton es convergente.
d) Aplicando Newton con x0 = 7 se ha obtenido, en la segunda iteración,

x2 = 60 3039. ¿Qué error se comete al aproximar la raı́z buscada por el
valor x3 que se obtiene en la siguiente iteración?
√
Ejercicio 1.8 En este ejercicio se pretende calcular 10 1 por el método de
Newton. Consideramos, para ello, la función f (x) = x10 − 1 cuya gráfica se da
en la Figura 1.
Fig. 1 Fig. 2
a) Probar, analı́ticamente, que en el intervalo [00 5, 10 5] posee una única raı́z

real.
b) Si tomamos x0 = 00 5 obtenemos la raı́z x = 1 en la iteración número 43,

mientras que si tomamos x0 = 10 5 se consigue el mismo resultado en la
iteración número 9. ¿Cómo podrı́amos haber conocido a priori el valor
que se debe elegir para x0 ?
c) ¿Sabrı́as justificar el porqué de la extremada lentitud de la convergencia

cuando iniciamos el proceso en x0 = 00 5? y ¿por qué sigue siendo lento
el proceso si comenzamos en x0 = 10 5? Justifica las respuestas.
d) Dado que en el intervalo [00 5, 10 5] no se anula la función x5 , las raı́ces de

f (x) son las mismas que las de g(x) = f (x)/x5 = x5 − x−5 cuya gráfica
se da en la Figura 2. ¿Se puede aplicar a g(x) la regla de Fourier en
dicho intervalo?
e) Si resolvemos, por el método de Newton, la ecuación g(x) = 0, ¿ se

obtendrá la raı́z con mayor rapidez que cuando lo hicimos con f (x) = 0?
Justifica la respuesta sin calcular las iteraciones.
Ejercicio 1.9 Dada la ecuación x7 − 14x + 7 = 0 se pide:
a) Probar que sólo tiene una raı́z real negativa.
b) Encontrar un entero a de tal forma que el intervalo [a, a + 1] contenga a

la menor de las raı́ces positivas de la ecuación.
c) ¿Cuál de los extremos del intervalo [a, a + 1] debe tomarse como valor
inicial para asegurar la convergencia del método de Newton?
d) Aplicar el método de Newton para obtener la menor de las raı́ces positivas

de la ecuación con seis cifras decimales exactas.
Ejercicio 1.10 Sea el polinomio p(x) = x4 − x2 + 1/8.
a) Utilizar el método de Sturm para determinar el número de raı́ces reales

positivas del polinomio p, ası́ como para separarlas.
b) Hallar los 2 primeros intervalos de la sucesión ([a1 , b1 ], [a2 , b2 ], . . .) obte-

nida de aplicar el método de dicotomı́a para obtener la mayor raı́z, r, del
polinomio p. Elegir el intervalo [a1 , b1 ] de amplitud 1/2 y tal que uno de
sus extremos sea un número entero.
c) Sea la sucesión definida por la recurrencia x0 = 1, xn+1 = F (xn ), donde

la iteración es la determinada por el método de Newton. Estudiar si
la regla de Fourier aplicada al polinomio p en el intervalo [a1 , b1 ] del
apartado anterior garantiza la convergencia de la sucesión a la raı́z r.
¿Y en el intervalo [a2 , b2 ]?
d) Hallar la aproximación x1 del apartado anterior, determinando una cota

del error cometido.
e) ¿Cuántas iteraciones se deben realizar para garantizar una aproximación

de r con veinte cifras decimales exactas?
1 n máx |f 00 (x)|
Indicación: En+1 = (kE1 )2 , con k = en un intervalo
k 2 mı́n |f 0 (x)|
adecuado.
Ejercicio 1.11 Dado el polinomio P (x) = x4 + 4x3 − 2x2 + 4x − 3 se pide:
a) Acotar las raı́ces y construir una sucesión de Sturm para probar que sólo
posee dos raı́ces reales, una positiva y otra negativa, dando intervalos de
amplitud 1 que las contengan.
b) Partiendo de que la raı́z positiva se encuentra en el intervalo (0, 1) y

despejando la x del término independiente
1 1 3
x = − x4 − x3 + x2 + ⇐⇒ x = ϕ(x)
4 2 4
¿se puede asegurar la convergencia de la sucesión x1 , x2 , . . . , xn , . . . defi-
nida de la forma x1 = 0, xn+1 = ϕ(xn ) ?
c) Aplicar Fourier para determinar el valor inicial que debe tomarse para
garantizar la convergencia del método de Newton en el cálculo de la raı́z
negativa. ¿Tenemos las tres cifras exactas si tomamos como raı́z -4.646 ?
Ejercicio 1.12 Sea el polinomio p(x) = −3 − x + x3 .
a) Utilizar una sucesión de Sturm para probar que el polinomio p(x) sólo
tiene una raı́z α ∈ R y que ésta se encuentra en el intervalo I = [0, 3].
b) Comprobar que la gráfica adjunta se corresponde con la función y =

ϕ(x) cuya iteración, xn+1 = ϕ(xn ) = xn − p(xn )/p0 (xn ), es la obtenida
con el método de Newton para resolver p(x) = 0. Tomando x1 = 0,
estudiar geométricamente (sobre el dibujo) si se obtendrı́a una sucesión
(xn ) convergente a α. ¿Y empezando en x1 = 3?
c) Tomar un subintervalo de I en el que la regla de Fourier garantice la

convergencia del Método de Newton y, con un valor inicial apropiado,
obtener una aproximación de α con, al menos, tres cifras decimales exac-
tas.
Ejercicio 1.13 Dado el polinomio P (x) = x3 + 6x2 + 9x + k con k ∈ R se

pide:
a) ¿Puede carecer de raı́ces reales? ¿y tener dos y sólo dos raı́ces reales?
b) Utilizar el método de Sturm para probar que tiene una única raı́z real
si, y sólo si, k < 0 o k > 4, y que sólo tiene raı́ces múltiples si k = 0 o
k = 4 no existiendo, en ningún caso, una raı́z triple.
c) Para k = −4 admite una única raı́z real en el intervalo [0, 1]. Si tomamos
como valor aproximado de la raı́z x = 00 3553 ¿de cuántas cifras decimales
exactas disponemos?
d) Si, en el caso anterior en que k = −4, aplicamos el método de Newton
para hallar la raı́z del polinomio, ¿cuál de los extremos del intervalo [0, 1]
deberı́amos tomar como valor inicial x0 para garantizar la convergencia?
y ¿qué valor obtendrı́amos para x2 ?
2. Sistemas de ecuaciones linea-
les
2.1 Normas vectoriales y matriciales

Sea E un espacio vectorial definido sobre un cuerpo K. Se define una norma
como una aplicación, que denotaremos por k k, de E en R que verifica las
siguientes propiedades:
1.- kxk ≥ 0 ∀x ∈ E siendo kxk = 0 ⇐⇒ x = 0. (Definida positiva).
2.- kc xk = |c| kxk ∀c ∈ K, ∀x ∈ E. (Homogeneidad).
3.- kx + yk ≤ kxk + kyk ∀x, y ∈ E. (Desigualdad triangular).
Un espacio E en el que hemos definido una norma recibe el nombre de espacio

normado.
Es frecuente que en el espacio E se haya definido también el producto de dos
elementos. En este caso, si se verifica que
kx · yk ≤ kxk kyk
se dice que la norma es multiplicativa. Esta propiedad de las normas es fun-

damental cuando trabajamos en el conjunto Cn×n de las matrices cuadradas
de orden n. Sin embargo no tiene mucha importancia cuando se trabaja en el
espacio C[a, b] de las funciones continuas en [a, b].
2.1.1 Normas vectoriales

Sea E un espacio normado de dimensión n y sea B = {u1 , u2 , . . . , un } una base
suya. Cualquier vector x ∈ E puede ser expresado de forma única en función
39
40 Sistemas de ecuaciones lineales
de los vectores de la base B.

n
X
x= xi ui
i=1
donde los escalares (x1 , x2 , . . . , xn ) se conocen como coordenadas del vector x

respecto de la base B.
Utilizando esta notación, son ejemplos de normas los siguientes:
n
X
• Norma-1 kxk1 = |xi |
i=1
v
u n
uX
• Norma euclı́dea kxk2 = t |xi |2
i=1
• Norma infinito kxk∞ = máx |xi |

i
2.1.2 Distancia inducida por una norma

Dado un espacio vectorial E, se define una distancia como una aplicación
d : E × E → R cumpliendo que:
• d(x, y) ≥ 0 ∀x, y ∈ E siendo d(x, y) = 0 ⇐⇒ x = y.
• d(x, y) = d(y, x) ∀x, y ∈ E.
• d(x, y) ≤ d(x, z) + d(z, y) ∀x, y, z ∈ E.

Si E, k k es un espacio normado, la norma k k induce una distancia en E
que se conoce como distancia inducida por la norma k k y viene definida por:
d(x, y) = kx − yk
Veamos que, en efecto, se trata de una distancia:
• d(x, y) ≥ 0 por tratarse de una norma, y además:
d(x, y) = 0 ⇐⇒ kx − yk = 0 ⇐⇒ x − y = 0 ⇐⇒ x = y.
• d(x, y) = kx − yk = k−1(y − x)k = |−1| ky − xk = ky − xk = d(y, x).
• d(x, y) = kx − yk = kx − z + z − yk ≤ kx − zk + kz − yk =
= d(x, z) + d(z, y).
Normas vectoriales y matriciales 41
2.1.3 Convergencia en espacios normados

Una sucesión de vectores v1 , v2 , . . . de un espacio vectorial normado (V, k k) se
dice que es convergente a un vector v si
lim kvk − vk = 0
k→∞
Esta definición coincide con la idea intuitiva de que la distancia de los vectores
de la sucesión al vector lı́mite v tiende a cero a medida que se avanza en la
sucesión.
Teorema 2.1 Para un espacio vectorial normado de dimensión finita, el con-

cepto de convergencia es independiente de la norma utilizada.
2.1.4 Normas matriciales

Dada una matriz A y un vector x, consideremos el vector transformado Ax.
La relación existente entre la norma del vector transformado y la del vector
original va a depender de la matriz A. El mayor de los cocientes entre dichas
normas, para todos los vectores del espacio, es lo que vamos a definir como
norma de la matriz A, de tal forma que de la propia definición se deduce que
kAxk ≤ kAkkxk
cualquiera que sea el vector x del espacio. (Obsérvese que no es lo mismo que
la propiedad multiplicativa de una norma, ya que aquı́ se están utilizando dos
normas diferentes, una de matriz y otra de vector).
Ası́ pues, definiremos
kAxk
kAk = máx = máx{kAxk : kxk = 1}
x∈V −{0} kxk
de tal forma que a cada norma vectorial se le asociará, de forma natural, una
norma matricial.
• Norma-1
Si utilizamos la norma-1 de vector obtendremos
kAk1 = máx{kAxk1 : kxk1 = 1}.

n
X
Dado que Ax = y =⇒ yi = aik xk se tiene que
k=1
n X
X n
kAk1 = máx{ |aik xk | : kxk1 = 1}.
i=1 k=1
Por último, si descargamos todo el peso sobre una coordenada, es decir,

si tomamos un vector de la base canónica, obtenemos que
X n
kAk1 = máx |aij |.
j
i=1
• Norma euclı́dea
Utilizando la norma euclı́dea de vector se tendrá que
√ √
kAk2 = máx{ x∗ A∗ Ax : x∗ x = 1}
Descargando ahora el peso en los autovectores de la matriz A∗ A obtene-
mos que
p √ p
kAk2 = máx{ x∗ λi x : x∗ x = 1} = máx λi = máx σi
i i i
donde σi representa a los valores singulares de la matriz A.

• Norma infinito
Utilizando ahora la norma infinito de vector se tiene que
Xn
kAk∞ = máx{ |aij xj | : kxk∞ = 1}.
j=1
Como ahora se dará el máximo en un vector que tenga todas sus coor-
denadas iguales a 1, se tiene que
n
X
kAk∞ = máx |aij |.
i
j=1
Tenemos pues, que las normas asociadas (algunas veces llamadas subordina-
das) a las normas de vector estudiadas anteriormente son:
n
X
Norma-1 kAk1 = máx kAxk1 = máx |aij |.
kxk1 =1 j
i=1
Norma euclı́dea kAk2 = máx kAxk2 = máx σi .

kxk2 =1 i
n
X
Norma infinito kAk∞ = máx kAxk∞ = máx |aij |.
kxk∞ =1 i
j=1
Normas vectoriales y matriciales 43
Si consideramos la matriz como un vector de m × n coordenadas, podemos

definir (de manera análoga a la norma euclı́dea de vector) la denominada
sX √
Norma de Frobenius kAkF = |aij |2 = tr A∗ A
i,j
2.1.5 Transformaciones unitarias

Una matriz U ∈ Cn×n se dice unitaria si
U ∗U = U U ∗ = I
es decir, si U ∗ = U −1 .
Proposicin 2.1 La norma de Frobenius y la norma euclı́dea de vector son

invariantes mediante transformaciones unitarias.
Demostración.
• Para la norma de Frobenius de matrices.

p p p
kU AkF = tr [(U A)∗ (U A)] = tr (A∗ U ∗ U A) = tr (A∗ A) = kAkF .
p p p
kA U kF = tr [(A U )(A U )∗ ] = tr (A U U ∗ A∗ ) = tr (A A∗ ) = kAkF .
• Para la norma euclı́dea de vector.

p √ √
kU xk2 = (U x)∗ (U x) = x∗ U ∗ U x = x∗ x = kxk2 .
p p
xT U 2 = (xT U )T 2 = U T x 2 = (U T x)∗ (U T x) = x∗ (U T )∗ U T x =
p p √ √
= xT (U ∗ )T U T x = x∗ (U U ∗ )T x = x∗ I T x = x∗ x = kxk2 .
Lema 2.2 Las matrices A, A∗ , AU y U A, donde U es una matriz unitaria,

poseen los mismos valores singulares.
Proposicin 2.2 La norma euclı́dea es invariante mediante transformaciones

de semejanza unitarias.
kAxk2
Demostración. Dado que kAk2 = máx si U es unitaria, se tiene
x∈V −{0} kxk2
que:
kU xk2 kxk2
kU k2 = máx = máx = 1.
x∈V −{0} kxk2 x∈V −{0} kxk2
Es decir, si U es unitaria kU k2 = kU ∗ k2 = 1.
Si B = U ∗ AU tenemos: kBk2 ≤ kU ∗ k2 kAk2 kU k2 = kAk2
Como A = U BU ∗ , es: kAk2 ≤ kU k2 kBk2 kU ∗ k2 = kBk2
De ambas desigualdades se deduce que kBk2 = kAk2 .
2.1.6 Radio espectral

Se define el radio espectral de una matriz A, y se denota por ρ(A) como el
máximo de los módulos de los autovalores de la referida matriz.
ρ(A) = máx |λi |

i
Geométricamente representa el radio del cı́rculo mı́nimo que contiene a todos

los autovalores de la matriz.
Teorema 2.3 El radio espectral de una matriz es una cota inferior de todas
las normas multiplicativas de dicha matriz.
Demostración. Sean {λ1 , λ2 , . . . , λr } los autovalores de la matriz A y sean

{x1 , x2 , . . . , xr } autovectores asociados a dichos autovalores.
kAxi k = kλi xi k = |λi | kxi k .
Por otra parte sabemos que kAxi k ≤ kAk kxi k. Por tanto:
|λi | kxi k ≤ kAk kxi k siendo kxi k 6= 0 por tratarse de autovectores. Se obtiene
entonces que |λi | ≤ kAk para cualquiera que sea i = 1, 2, . . . , r, de donde
máx |λi | ≤ kAk, es decir: ρ(A) ≤ kAk.
i
2.2 Sistemas de ecuaciones lineales

En el capı́tulo anterior se estudiaron métodos iterados para la resolución de
ecuaciones no lineales. Dichos métodos se basaban en el teorema del punto
fijo y consistı́an en expresar la ecuación en la forma x = ϕ(x) exigiendo que
ϕ0 (x) ≤ q < 1 para cualquier x del intervalo en el cual se trata de buscar la
solución.
Sistemas de ecuaciones lineales 45
Para los sistemas de ecuaciones lineales, de la forma Ax = b, trataremos de

buscar métodos iterados de una forma análoga a como se hizo en el caso de las
ecuaciones, es decir, transformando el sistema en otro equivalente de la forma
x = F (x) donde F (x) = M x + N . Evidentemente habrá que exigir algunas
condiciones a la matriz M para que el método sea convergente (al igual que se
exigı́a que ϕ0 (x) ≤ q < 1 en el caso de las ecuaciones) y estas condiciones se
basan en los conceptos de normas vectoriales y matriciales.
Dada una aplicación f : Rm → Rn y un vector b ∈ Rn , resolver el sistema
de ecuaciones f (x) = b no es más que buscar el conjunto de vectores de Rm
cuya imagen mediante f es el vector b, es decir, buscar la imagen inversa de b
mediante f .
Un sistema de ecuaciones se dice lineal en su componente k-ésima si verifica
que
f (x1 , . . . , xk−1 , αx1k + βx2k , xk+1 , . . . , xm ) = αf (x1 , . . . , xk−1 , x1k , xk+1 , . . . , xm )+
+ βf (x1 , . . . , xk−1 , x2k , xk+1 , . . . , xm )
Diremos que un sistema es lineal si lo es en todas sus componentes, pudién-
dose, en este caso, escribir de la forma Ax = b.
Si la aplicación f se define de Cm en Cn resulta un sistema complejo que
puede ser transformado en otro sistema real. Ası́, por ejemplo, si el sistema
es lineal, es decir, de la forma M z = k con M ∈ Cm×n , x ∈ Cn×1 y k ∈
Cm×1 , podemos descomponer la matriz M en suma de otras dos de la forma
M = A + iB con A, B ∈ Rm×n y análogamente z = x + iy con x, y ∈ Rn×1
k = k1 + ik2 con k1 , k2 ∈ Rm×1 , por lo que (A + iB)(x + iy) = k1 + ik2 es decir
( ! ! !
Ax − By = k1 A −B x k1
=⇒ =
Bx + Ay = k2 B A y k2
sistema real de 2m ecuaciones con 2n incógnitas. Es por ello, que centraremos

nuestro estudio en los sistemas reales.
Podemos clasificar los sistemas de ecuaciones lineales atendiendo a
a) Su tamaño
a.1) Pequeños: n ≤ 300 donde n representa el número de ecuaciones.

a.2) Grandes: n > 300
(Esta clasificación corresponde al error de redondeo)

b) Su estructura
b.1) Si la matriz posee pocos elementos nulos diremos que se trata de

un sistema lleno.
b.2) Si, por el contrario, la matriz contiene muchos elementos nulos,
diremos que la matriz y, por tanto, que el sistema es disperso o
sparce. Matrices de este tipo son las denominadas
 
a11 a12 0 0
 a21 a22 a23 0 
• Tridiagonales: 
 0 a32 a33 a34 

0 0 a43 a44
 
a11 a12 a13 a14
 0 a22 a23 a24 
• Triangulares superiores:   0

0 a33 a34 
0 0 0 a44
 
a11 0 0 0
 a12 a22 0 0 
• Triangulares inferiores: 
 a31 a32 a33 0 

a41 a42 a43 a44
En cuanto a los métodos de resolución de sistemas de ecuaciones lineales,

podemos clasificarlos en
a) Métodos directos
b) Métodos iterados
Se denominan métodos directos a aquellos métodos que resuelven un sistema

de ecuaciones lineales en un número finito de pasos. Se utilizan para resolver
sistemas pequeños.
Los denominados métodos iterados crean una sucesión de vectores que con-
vergen a la solución del sistema. Estos métodos se utilizan para la resolución
de sistemas grandes, ya que al realizar un gran número de operaciones los
errores de redondeo pueden hacer inestable al proceso, es decir, pueden alterar
considerablemente la solución del sistema.
2.2.1 Número de condición

Un sistema de ecuaciones lineales Ax = b se dice bien condicionado cuando
los errores cometidos en los elementos de la matriz A y del vector b producen
en la solución un error del mismo orden, mientras que diremos que el sistema
está mal condicionado si el error que producen en la solución del sistema es de
orden superior al de los datos. Es decir:

A−A <ε  kx − xk ' ε sistema bien condicionado
=⇒
b−b <ε 
kx − xk ε sistema mal condicionado
Consideremos el sistema cuadrado Ax = b con A regular, es decir, un sistema

compatible determinado. En la práctica, los elementos de A y de b no suelen
ser exactos bien por que procedan de cálculos anteriores, o bien porque sean
irracionales, racionales periódicos, etc. Es decir, debemos resolver un sistema
aproximado cuya solución puede diferir poco o mucho de la verdadera solución
del sistema.
Ası́, por ejemplo, en un sistema de orden dos, la solución representa el punto
de intersección de dos rectas en el plano. Un pequeño error en la pendiente de
una de ellas puede hacer que dicho punto de corte se desplace sólo un poco o
una distancia considerable (véase la Figura 2.1), lo que nos dice que el sistema
está bien o mal condicionado, respectivamente.
Podemos ver que el sistema está mal condicionado cuando las pendientes
de las dos rectas son muy similares y que mientras más ortogonales sean las
rectas, mejor condicionado estará el sistema.
Se puede observar entonces que si, en un sistema mal condicionado, susti-
tuimos una de las ecuaciones por una combinación lineal de las dos, podemos
hacer que el sistema resultante esté bien condicionado.
r1 r10 r1 r10
%
%r
@ %
@ %
@r r

%
@ %
@
r r
%
@ 2 %
@
%
r2%
%

Sistema bien condicionado Sistema mal condicionado
Figura 2.1: Condicionamiento de un sistema.
Ejemplo 2.1 Si consideramos el sistema
3x + 4y = 7
x

1

de solución =
3x + 4.00001y = 7.00001 y 1
y cometemos un pequeño error en los datos, podemos obtener el sistema
3x + 4y = 7
x

7.6̄

de solución =
3x + 3.99999y = 7.00004 y −4
o bien este otro

3x + 4y = 7
x

9.6̄

de solución =
3x + 3.99999y = 7.000055 y −5.5
lo que nos dice que estamos ante un sistema mal condicionado.

Si sustituimos la segunda ecuación por la que resulta de sumarle la primera
multiplicada por −10 0000016 (la ecuación resultante se multiplica por 106 y se
divide por −10 2) nos queda el sistema
3x + 4y = 7
x

1

de solución =
4x − 3y = 1 y 1
siendo éste un sistema bien condicionado.
El estudio del condicionamiento de un sistema se realiza a través del deno-

minado número de condición que estudiamos a continuación.
Sea A una matriz cuadrada y regular. Se define el número de condición de
la matriz A y se denota por κ(A) como
κ(A) = kAk · A−1
donde la norma utilizada ha de ser una norma multiplicativa. Este número

nos permite conocer el condicionamiento del sistema Ax = b.
Dado que en la práctica el cálculo de la matriz inversa A−1 presenta grandes
dificultades lo que se hace es buscar una cota del número de condición.
κ(A) = kAk · A−1 < kAk · k
siendo k una cota de la norma de la matriz inversa.

kIk
Si kI − Ak < 1 entonces kA−1 k ≤ . En efecto:
1 − kI − Ak
A · A−1 = I =⇒ [I − (I − A)]A−1 = I =⇒
A−1 − (I − A)A−1 = I =⇒ A−1 = I + (I − A)A−1 =⇒
kA−1 k = kI + (I − A)A−1 k ≤ kIk+k(I − A)A−1 k ≤ kIk+kI − Ak kA−1 k ⇒
kA−1 k − kI − Ak kA−1 k ≤ kIk =⇒ (1 − kI − Ak) kA−1 k ≤ kIk =⇒
kIk
A−1 ≤
1 − kI − Ak
Es decir:
kIk
κ(A) ≤ kAk · k con k=
1 − kI − Ak
Debemos tener cuidado con esta acotación ya que si tenemos una matriz casi
regular, es decir, con det(A) ' 0, quiere decir que tiene un autovalor próximo
a cero, por lo que la matriz I − A tiene un autovalor próximo a 1 y será el
mayor de todos. En este caso kI − Ak ' 1, por lo que k → ∞ y darı́a lugar a
un falso condicionamiento, ya que A no tiene que estar, necesariamente, mal
condicionada.
Ejemplo 2.2 Para estudiar el condicionamiento del sistema
3x + 4y = 7
3x + 4.00001y = 7.00001
Se tiene que
3 4
A= =⇒ det(A) = 0.00003
3 4.00001

−1 1 4.00001 −4
A =
0.00003 −3 3
n
X
Utilizando la norma infinito kAk∞ = máx |aij | se tiene que
i
j=1

kAk∞ = 7.00001 
 56
8.00001  =⇒ κ∞ (A) ' · 105 > 1.8 · 106
kA−1 k∞ =  3
0.00003
Se trata pues, de un sistema mal condicionado.
n
X
Si utilizamos la norma-1 kAk1 = máx |aij | obtenemos:
j
i=1

kAk1 = 8.00001 
 56
7.00001 =⇒ κ1 (A) ' · 105 > 1.8 · 106
kA−1 k1 = 
 3
0.00003
obteniéndose, también, que se trata de un sistema mal condicionado.
Propiedades del número de condición κ(A).
a) Dado que kxk = kIxk ≤ kIkkxk, se verifica que kIk ≥ 1, cualquiera que
sea la norma utilizada. Como, por otra parte AA−1 = I se tiene que
1 ≤ kIk = kAA−1 k ≤ kAkkA−1 k = κ(A)
es decir κ(A) ≥ 1 cualquiera que sea la matriz cuadrada y regular A.
b) Si B = zA, con z ∈ C no nulo, se verifica que κ(B) = κ(A). En efecto:

1 −1 kA−1 k
κ(B)=kBk B −1 = kzAk A = |z| kAk = kAk A−1 = κ(A).
z |z|
Dado que det(B) = z n det(A), donde n representa el orden de la matriz
A, y κ(B) = κ(A) se ve que el condicionamiento de una matriz no
depende del valor de su determinante.
σn
c) Utilizando la norma euclı́dea k k2 se tiene que κ2 (A) = donde σ1
σ1
y σn representan, respectivamente, al menor y al mayor de los valores
singulares de la matriz A.
En efecto: sabemos que los valores singulares σi de la matriz A son las
raı́ces cuadradas positivas de los autovalores de la matriz A∗ A.
p
σi = λi (A∗ A)
Si suponemos σ1 ≤ σ2 ≤ · · · ≤ σn se tiene que
q
kAk2 = máx λi (A∗ A) = σn
i
q q
kA−1 k2 = máx λi (A−1 )∗ A−1 = máx λi (A∗ )−1 A−1 =

i i
r r
q
∗ −1 1 1
= máx λi (AA ) = máx ∗
= =
i i λi (AA ) mı́ni λi (AA∗ )
r r
1 1
= = =⇒
mı́ni λi (A∗ A) mı́ni σi2
1
A−1 2
=
σ1
Podemos concluir, por tanto, que


kAk2 = σn 
 σn
1 =⇒ κ2 (A) =
kA−1 k2 = 
 σ1
σ1
En cuanto a su relación con los números de condición obtenidos con otras

normas de matriz se tiene que:
√ √
kAk2 ≤ kAkF ≤ n kAk2 =⇒ kA−1 k2 ≤ kA−1 kF ≤ n kA−1 k2
κ2 (A) = kAk2 kA−1 k2 ≤ kAkF kA−1 kF = κ(A)F
√ √
κF (A) = kAkF kA−1 kF ≤ n n kAk2 kA−1 k2 = nκ2 (A) =⇒
κ2 (A) ≤ κF (A) ≤ nκ2 (A)

 p
 kAk2 ≤ kAk1 kAk∞ p
Además: p ⇒ κ2 (A) ≤ κ1 (A)κ∞ (A)
 kA−1 k ≤ kA−1 k1 kA−1 k∞
2
d) Una condición necesaria y suficiente para que κ2 (A) = 1 es que A = zU

siendo z ∈ C (no nulo) y U una matriz unitaria (U U ∗ = U ∗ U = I).
⇐) A = zU =⇒ κ2 (A) = 1. En efecto:
A = zU =⇒ A∗ A = zU ∗ zU = |z|2 U ∗ U = |z|2 I =⇒
λi (A∗ A) = |z|2 cualquiera que sea i = 1, 2, . . . , n y, por tanto,
σ1 = σ2 = · · · = σn = |z|
por lo que
σn
κ2 (A) = =1
σ1
⇒) κ2 (A) = 1 =⇒ A = zU .
En efecto: sabemos que si A es diagonalizable existe una matriz
regular R tal que R−1 AR = D con D = diag(λi ) (R es la matriz
de paso cuyas columnas son los autovectores correspondientes a los
autovalores λi ). Por otra parte sabemos que toda matriz hermı́tica

es diagonalizable mediante una matriz de paso unitaria.
Como la matriz A∗ A es hermı́tica y tiene los mismos autovalores
que AA∗ . Existe entonces una matriz R tal que
 
σ12
−1 ∗
 σ22 
R A AR = 
 
. . 
 . 
σn2
σn
Como κ2 (A) = = 1 =⇒ σ1 = σ2 = · · · = σn = σ, por lo que
σ1
R−1 A∗ AR = σ 2 I
Entonces
A∗ A = R(σ 2 I)R−1 = σ 2 (RIR−1 ) = σ 2 I
Ahora bien
1 ∗ 1
A A =I
σ σ
1 1
Llamando U = A se tiene que U ∗ = A∗ , ya que σ ∈ R ⇒ σ = σ.
σ σ
∗ 1 ∗ 1
Se tiene entonces que A = σU con U U = A A = I, es
σ σ
decir, con U unitaria.
Los sistemas mejor condicionados son aquellos que tienen sus filas o columnas
ortogonales y mientras mayor sea la dependencia lineal existente entres ellas
peor es el condicionamiento del sistema.
Al ser κ(AU ) = κ(U A) = κ(A) trataremos de buscar métodos de resolución
de sistemas de ecuaciones lineales que trabajen con matrices unitarias que
no empeoren el condicionamiento del sistema como lo hace, por ejemplo, el
método de Gauss basado en la factorización LU . Sin embargo, dado que ha
sido estudiado en la asignatura de Álgebra Lineal, comenzaremos estudiando
dicho método aunque pueda alterarnos el condicionamiento del problema.
Empezaremos estudiando pues, como métodos directos, los basados en la
factorización LU y la de Cholesky.
Factorización LU 53
2.3 Factorización LU
Al aplicar el método de Gauss al sistema Ax = b realizamos transformaciones
elementales para conseguir triangularizar la matriz del sistema. Si este proceso
puede realizarse sin intercambios de filas, la matriz triangular superior U obte-
nida viene determinada por el producto de un número finito de transformacio-
nes fila Fk Fk−1 · · · F1 aplicadas a la matriz A. Llamando L−1 = Fk Fk−1 · · · F1
(ya que el determinante de una transformación fila es ±1 y, por tanto, su pro-
ducto es inversible) se tiene que L−1 A = U , o lo que es lo mismo, A = LU .
Además, la matriz L es una triangular inferior con unos en la diagonal.
Esta factorización es única ya que de existir otra tal que A = L0 U 0 = LU
se tendrı́a que L−1 L0 = U U 0−1 . Como L−1 también es triangular inferior
con unos en la diagonal, el producto L−1 L0 también es una matriz del mismo
tipo. Análogamente, el producto U U 0−1 resulta ser una triangular superior. El
hecho de que L−1 L0 = U U 0−1 nos dice que necesariamente L−1 L0 = I, ya que
es simultáneamente triangular inferior y superior y su diagonal es de unos. Ası́
pues L−1 L0 = I, por lo que L = L0 y, por tanto U = U 0 es decir, la factorización
es única.
Debido a la unicidad de la factorización, ésta puede ser calculada por un
método directo, es decir, haciendo
  
1 0 0 ··· 0 u11 u12 u13 · · · u1n
 l21 1 0 · · · 0   0 u22 u23 · · · u2n 
  
 l31 l32 1 · · · 0   0 0 u33 · · · u3n
A=

 
 .. .. .. . . ..   .. .. .. .. .
. ..

 . . . . .  . . . 
ln1 ln2 ln3 · · · 1 0 0 0 · · · unn
y calculando los valores de los n2 elementos que aparecen entre las dos matrices.
 
3 1 2
Ası́, por ejemplo, para A =  6 3 2  tenemos
−3 0 −8
    
3 1 2 1 0 0 u11 u12 u13
 6 3 2  =  l21 1 0   0 u22 u23  =
−3 0 −8 l31 l32 1 0 0 u33
 
u11 u12 u13
=  l21 u11 l21 u12 + u22 l21 u13 + u23 
l31 u11 l31 u12 + l32 u22 l31 u13 + l32 u23 + u33
por lo que de la primera fila obtenemos que
u11 = 3 u12 = 1 u13 = 2
de la segunda (teniendo en cuenta los valores ya obtenidos) se tiene que


3l21 = 6  l21 = 2
l21 + u22 = 3 =⇒ u22 = 1
2l21 + u23 = 2 u23 = −2

y de la tercera (teniendo también en cuenta los resultados ya obtenidos)


3l31 = −3  l31 = −1
l31 + l32 = 0 =⇒ l32 = 1
2l31 − 2l32 + u33 = −8 u33 = −4

es decir:
    
3 1 2 1 0 0 3 1 2
 6 3 2  =  2 1 0   0 1 −2 
−3 0 −8 −1 1 1 0 0 −4
Se denominan matrices fundamentales de una matriz A, y se denotan por

Ak , a las submatrices constituidas por los elementos de A situados en las k
primeras filas y las k primeras columnas, es decir:
 
a11 a12 a13
a11 a12
A1 = (a11 ) A2 = A3 =  a21 a22 a23 
a21 a22
a31 a32 a33
Teorema 2.4 Una matriz regular A admite factorización LU si, y sólo si, sus
matrices fundamentales Ai (i = 1, 2, . . . , n) son todas regulares.
Demostración. Supongamos que A admite factorización LU . En ese caso

Ak Lk Uk
A= = =⇒
Ak = Lk Uk =⇒ det(Ak ) = det(Lk ) det(Uk ) = 1 · r11 r22 · · · rkk 6= 0

ya que, por sea A regular, todos los pivotes rii i = 1, 2, . . . , n son no nulos.
Recı́procamente, si todas las matrices fundamentales son regulares, A ad-
mite factorización LU , o lo que es equivalente, se puede aplicar Gauss sin
intercambio de filas. En efecto:
Dado que, por hipótesis es a11 6= 0, se puede utilizar dicho elemento como
pivote para anular al resto de los elementos de su columna quedándonos la
matriz  (2) (2) (2) 
a11 a12 · · · a1n
 0 a(2) · · · a(2) 
(2) 22 2n 
A = .

. .
. . . . 
 . . . .. 
(2) (2)
0 an2 · · · ann
Factorización LU 55
(2)
donde a1i = a1i para i = 1, 2, . . . , n.
a21
Si nos fijamos ahora en a222 = a22 − a12 podemos ver que es no nulo, ya
a11
que de ser nulo serı́a
a11 a22 − a12 a21 = det(A2 ) = 0
en contra de la hipótesis de que todas las matrices fundamentales son regulares.

(2)
Por tanto, podemos utilizar a22 6= 0 como nuevo pivote para anular a los
elementos de su columna situados bajo él.
Reiterando el procedimiento se puede ver que todos los elementos que vamos
obteniendo en la diagonal son no nulos y, por tanto, válidos como pivotes. Es
decir, puede aplicarse el método de Gauss sin intercambio de filas.
Comprobar si una matriz admite factorización LU estudiando si todas sus
matrices fundamentales son regulares es un método demasiado costoso debido
al número de determinantes que hay que calcular.
Definición 2.1 Una matriz cuadrada de orden n A = (aij )i = 1, 2, . . . , n se dice

j = 1, 2, . . . , n
que es una matriz de Hadamard o de diagonal dominante :
n
X
a) Por filas: si |aii | > |aik | i = 1, 2, . . . , n
k=1
k 6= i
n
X
b) Por columnas: si |aii | > |aki | i = 1, 2, . . . , n
k=1
k 6= i
 
3 1 1
Ası́, por ejemplo, la matriz A =  0 2 1  es de diagonal dominante
2 −1 5
por filas pero no por columnas.
Teorema 2.5 Toda matriz de Hadamard es regular.
Demostración. Supongamos que A es de diagonal dominante por filas

(de igual forma podrı́a probarse si lo fuese por columnas) y que su deter-
minante fuese nulo. En ese caso, el sistema Ax = 0 posee solución no trivial
(α1 , α2 , . . . , αn ) 6= (0, 0, . . . , 0).
Sea |αk | = máx |αi | > 0 y consideremos la k-ésima ecuación:
i
ak1 α1 + ak2 α2 + · · · + akk αk + · · · + akn αn = 0 =⇒
α1 α2 αn
ak1 + ak2 + · · · + akk + · · · + akn = 0 =⇒
αk αk αk
α1 αk − 1 αk + 1 αn
akk = −ak1 − · · · − ak k−1 − ak k+1 − · · · − akn =⇒
αk αk αk αk
n n
X αi X
|akk | ≤ |aki | ≤ |aki |
i=1
αk i=1
i 6= k i 6= k
en contra de la hipótesis de que A es de diagonal dominante por filas. Por

tanto, toda matriz de Hadamard, es regular.
Teorema 2.6 Las matrices fundamentales Ak de una matriz A de Hadamard,

son también de Hadamard.
Demostración. La demostración es trivial, ya que si A es de Hadamard se

verifica que
n
X X k
|aii | > |aij | ≥ |aij |
j =1 j =1
j 6= i j 6= i
luego Ak también lo es.

Como consecuencia de los Teoremas 2.4, 2.6 y 2.5, podemos deducir el si-
guiente corolario.
Corolario 2.7 Toda matriz de Hadamard admite factorización LU .
Otro tipo de matrices de las que se puede asegurar que admiten factorización
LU son las hermı́ticas definidas positivas, ya que las matrices fundamentales de
éstas tienen todas determinante positivo, por lo que el Teorema 2.4 garantiza
la existencia de las matrices L y U .
2.4 Factorización de Cholesky

Una vez visto el método de Gauss basado en la factorización LU vamos a
estudiar otros métodos que se basan en otros tipos de descomposiciones de la
matriz del sistema.
Factorización de Cholesky 57
Es conocido que toda matriz hermı́tica y definida positiva tiene sus autova-
lores reales y positivos y, además, en la factorización LU todos los pivotes son
reales y positivos.
Teorema 2.8 [Factorización de Cholesky] Toda matriz A hermı́tica y

definida positiva puede ser descompuesta de la forma A = BB ∗ siendo B una
matriz triangular inferior.
Demostración. Por tratarse de una matriz hermı́tica y definida positiva,

sabemos que admite factorización LU . Sea
  
1 0 ··· 0 u11 u12 · · · u1n
 l21 1 · · · 0   0 u22 · · · u2n 
A =  .. ..  =
  
.. . . ..   .. .. . .
 . . . .  . . . . 
ln1 ln2 · · · 1 0 0 · · · unn
   
1 0 ··· 0 u11 0 · · · 0 1 u12/u11 · · · u1n/u11
 l21 1 · · · 0   0 u22 · · · 0   0 1 · · · u2n/u22 
=  .. ..  =
   
.. . . ..   .. .. . . ..   .. .. . .
 . . . . . . . . . . . . 
ln1 ln2 · · · 1 0 0 · · · unn 0 0 ··· 1
 
u11 0 · · · 0
 0 u22 · · · 0 
= L  .. ..  R =
 
.. . .
 . . . . 
0 0 · · · unn
 √  √ 
u11 0 ··· 0 u11 0 ··· 0
√ √
 0 u22 · · · 0   0
 u22 · · · 0 
= L  .. ..  R ⇒
 
.. . . .
.   .
. .
. . .
 . . . .  . . . . 
√ √
0 0 ··· unn 0 0 ··· unn
A = BC
 √ 
u11 0 0 ··· 0
√
 0 u22 0 ··· 0 
 √ 
donde B = L 
 0 0 u33 · · · 0  es una matriz triangu-

 .. .. .. .. .. 
. . . . .
√
 
0 0 0 ··· unn
 √ 
u11 0 0 ··· 0
√
 0 u 22 0 ··· 0 
 √ 
lar inferior y C =
 0
 0 u33 ··· 0  R es una triangular

 .. .. .. ... .. 
 . . . .
√

0 0 0 ··· unn
superior.
Como A es hermı́tica, BC = A = A∗ = C ∗ B ∗ , por lo que (C ∗ )−1 B = B ∗ C −1 ,
y dado que (C ∗ )−1 B es triangular inferior y B ∗ C −1 es triangular superior,
ambas han de ser diagonales.
Por otra parte, B = LD y C = DR, por lo que C ∗ = R∗ D∗ = R∗ D y, por
tanto, (C ∗ )−1 = D−1 (R∗ )−1 . Ası́ pues, (C ∗ )−1 B = D−1 (R∗ )−1 LD.
Como las matrices diagonales conmutan,
(C ∗ )−1 B = D−1 D(R∗ )−1 L = (R∗ )−1 L.
Al ser (R∗ )−1 L triangular inferior con diagonal de unos y (C ∗ )−1 B diagonal,
podemos asegurar que (R∗ )−1 L = I o, lo que es lo mismo, R∗ = L. Además,
B ∗ C −1 = I, por lo que C = B ∗ , luego A = BB ∗ donde B = LD.
La unicidad de las matrices L y U implica la unicidad de la matriz B y, por
tanto, ésta puede ser calculada por un método directo.
Ejemplo 2.3 Consideremos el sistema

    
4 2i 4 + 2i x1 0
 −2i 2 2 − 2i   x2  =  0 
4 − 2i 2 + 2i 10 x3 −4
Realicemos la factorización BB ∗ directamente, es decir

    
4 2i 4 + 2i b11 0 0 b11 b21 b31
 −2i 2 2 − 2i  =  b21 b22 0   0 b22 b32 
4 − 2i 2 + 2i 10 b31 b32 b33 0 0 b33
Se obtiene multiplicando, que |b11 |2 = 4 por lo que b11 = 2. Utilizando este

resultado tenemos que 2b21 = 2i, por lo que b21 = −i y que 2b31 = 4 + 2i por
lo que b31 = 2 − i.
Por otro lado, |b21 |2 + |b22 |2 = 2, por lo que |b22 |2 = 1 y, por tanto, b22 = 1.
Como b21 b31 + b22 b32 = 2 − 2i tenemos que 1 − 2i + b32 = 2 − 2i, es decir,
b32 = 1 y, por tanto, b32 = 1.
Por último, |b31 |2 + |b32 |2 + |b33 |2 = 10, por lo que 5 + 1 + |b33 |2 = 10, es decir
|b33 |2 = 4 y, por tanto, b33 = 2. Ası́ pues, el sistema nos queda de la forma
     
2 0 0 2 i 2+i x1 0
 −i 1 0   0 1 1   x2  =  0 
2−i 1 2 0 0 2 x3 −4
Métodos iterados 59
    
2 0 0 y1 0
Haciendo ahora  −i 1 0   y2  =  0 , se obtiene
2−i 1 2 y3 −4
   
y1 0
 y2  =  0 
y3 −2
y de aquı́, que     
2 i 2+i x1 0
 0 1 1   x2  =  0 
0 0 2 x3 −2
de donde obtenemos que la solución del sistema es
x1 = 1 x2 = 1 x3 = −1
Hemos visto que toda matriz hermı́tica y definida positiva admite factori-
zación de Cholesky, pero podemos llegar más lejos y enunciar el siguiente
teorema (que no probaremos).
Teorema 2.9 Una matriz hermı́tica y regular A es definida positiva si, y sólo
si, admite factorización de Cholesky.
2.5 Métodos iterados

Un método iterado de resolución del sistema Ax = b es aquel que genera, a
partir de un vector inicial x0 , una sucesión de vectores x1 , x2 , . . .. El método se
dirá que es consistente con el sistema Ax = b, si el lı́mite de dicha sucesión, en
caso de existir, es solución del sistema. Se dirá que el método es convergente si
la sucesión generada por cualquier vector inicial x0 es convergente a la solución
del sistema.
Es evidente que si un método es convergente es consistente, sin embargo, el
recı́proco no es cierto como prueba el siguiente ejemplo.
Ejemplo 2.4 El método xn+1 = 2xn − A−1 b es consistente con al sistema

Ax = b pero no es convergente. En efecto:
xn+1 − x = 2xn − A−1 b − x = 2xn − 2x − A−1 b + x = 2(xn − x) − (A−1 b − x)

y como A−1 b = x, se tiene que
xn+1 − x = 2(xn − x)
Si existe lim xn = x∗ tendremos que

n→∞
x∗ − x = 2(x∗ − x) =⇒ x∗ − x = 0 =⇒ x∗ = x
es decir, el lı́mite es solución del sistema Ax = b, por lo que el método es

consistente.
Sin embargo, de xn+1 − x = 2(xn − x) obtenemos que
kxn+1 − xk = 2kxn − xk
es decir, el vector xn+1 dista de x el doble de lo que distaba xn , por lo que el

método no puede ser convergente.
Los métodos iterados que trataremos son de la forma
xn+1 = Kxn + c
en los que K será la que denominemos matriz del método y que dependerá de
A y de b y en el que c es un vector que vendrá dado en función de A, K y b.
Teorema 2.10 Un método iterado, de la forma xn+1 = Kxn +c, es consistente

con el sistema Ax = b si, y sólo si, c = (I − K)A−1 b y la matriz I − K es
invertible
Demostración.
a) Supongamos que el método es consistente con el sistema Ax = b.

Como x = Kx + (I − K)x = Kx + (I − K)A−1 b, se tiene que
xn+1 − x = K(xn − x) + c − (I − K)A−1 b (2.1)
Por ser consistente el método, de existir x∗ = lim xn ha de ser x∗ = x.

n→∞
Pasando al lı́mite en la Ecuación (2.1) obtenemos que
x∗ − x = K(x∗ − x) + c − (I − K)A−1 b
por lo que
(I − K)(x∗ − x) = c − (I − K)A−1 b (2.2)
y dado que x∗ = x nos queda que 0 = c − (I − K)A−1 b, es decir,
c = (I − K)A−1 b.
Además, dado que x = Kx + c, el sistema (I − K)x = c posee solución

única x y, por tanto, la matriz I − K es invertible.
b) Si c = (I − K)A−1 b y la matriz I − K es invertible, cuando exista

lim xn = x∗ se tendrá de (2.2) que
n→∞
(I − K)(x∗ − x) = 0
y como I − K es invertible, x∗ = x, por lo que el método es consistente.
Teorema 2.11 Un método iterado de la forma xn+1 = Kxn + c consistente

con el sistema Ax = b es convergente si, y sólo si, lim K n = 0.
n→∞
Demostración.
a) Por tratarse de un método consistente con el sistema Ax = b, se verifica

que c = (I − K)A−1 b, por lo que
xn+1 = Kxn + (I − K)A−1 b
restando el vector solución x a ambos miembros, podemos escribir
xn+1 − x = Kxn − (K + I − K)x + (I − K)A−1 b =
= K(xn − x) + (I − K)(A−1 b − x)
y dado que A−1 b − x = 0 obtenemos que xn+1 − x = K(xn − x).
Reiterando el proceso se obtiene:
xn − x = K(xn−1 − n) = K 2 (xn−2 − x) = · · · = K n (x0 − x)
Pasando al lı́mite
lim (xn − x) = (x0 − x) lim K n

n→∞ n→∞
Al suponer el método convergente, lim (xn − x) = x − x = 0, por lo que

n→∞
lim K n = 0
n→∞
b) Recı́procamente, si limn→∞ K n = 0, obtenemos que
lim (xn − x) = 0
n→∞
o lo que es lo mismo,
lim xn = x
n→∞
por lo que el método es convergente.
Teorema 2.12 Si para alguna norma matricial es kKk < 1, el proceso xn+1 =
Kxn + c, donde x0 ∈ Rn es un vector cualquiera, converge a la solución de la
ecuación x = Kx + c que existe y es única.
Demostración.
a) Veamos, en primer lugar, que la ecuación x = Kx + c posee solución

única.
En efecto: x = Kx + c =⇒ (I − K)x = c. Este sistema tiene solución
única si, y sólo si, el sistema homogéneo asociado (I − K)z = 0 admite
sólo la solución trivial z = 0, es decir, si I − K es invertible.
La solución z no puede ser distinta del vector nulo ya que de serlo, como
kKk < 1 se tiene que al ser (I − K)z = 0, o lo que es lo mismo, z = Kz
kzk = kKzk ≤ kKkkzk < kzk
lo cual es un absurdo, por lo que el sistema homogéneo sólo admite la

solución trivial y, por tanto, el sistema completo x = Kx + c posee
solución única.
b) Probaremos ahora que la sucesión {xn } converge a x.

Dado que xn+1 − x = (Kxn + c) − (Kx + c) = K(xn − x), podemos
reiterar el proceso para obtener que xn − x = K n (x0 − x) por lo que
kxn − xk = kK n kkx0 − xk ≤ kKkn kx0 − xk
y dado que kKk < 1, pasando al lı́mite se obtiene
lim kxn − xk = 0 =⇒ lim xn = x

n→∞ n→∞
Los métodos que vamos a estudiar consisten en descomponer la matriz inver-

tible A del sistema Ax = b de la forma A = M − N de manera que la matriz
M sea fácilmente invertible, por lo que reciben el nombre genérico de métodos
de descomposición . El sistema queda entonces de la forma
(M − N )x = b =⇒ M x = N x + b =⇒ x = M −1 N x + M −1 b
es decir, expresamos el sistema de la forma x = Kx + c con K = M −1 N y
c = M −1 b. Dado que
(I −K)A−1 b = (I −M −1 N )(M −N )−1 b = M −1 (M −N )(M −N )−1 b = M −1 b = c
y la matriz (I − K) = (I − M −1 N ) = M −1 (M − N ) = M −1 A es invertible,
estamos en las condiciones del Teorema 2.10 por lo que el método xn+1 =
Kxn + c es consistente con el sistema Ax = b. Es decir, si el proceso converge,
lo hace a la solución del sistema.
Sabemos también, por el Teorema 2.12, que el proceso será convergente si se
verifica que kM −1 N k < 1 para alguna norma.
Para el estudio de los métodos que trataremos a continuación, vamos a des-
componer la matriz A de la forma A = D − E − F siendo
   
a11 0 0 ··· 0 0 0 0 ··· 0
 0 a22
 0 ··· 0  
 a21 0
 0 ··· 0 

D= 0
 0 a33 · · · 0   − E =  a31 a32 0
 ··· 0 

 .. .. .. . . .   .. .. .. .. ..
. .. 

 . . .  . . . . . 
0 0 0 · · · ann an1 an2 an3 ··· 0
 
0 a12 a13 · · · a1n
 0 0 a23 · · · a2n 
 
−F =  0 0
 0 · · · a3n  
 .. .. .. . . .. 
 . . . . . 
0 0 0 ··· 0
2.5.1 Método de Jacobi

Consiste en realizar la descomposición A = M − N = D − (E + F ). El sistema
queda de la forma
Ax = b =⇒ Dx = (E + F )x + b =⇒ x = D−1 (E + F )x + D−1 b
La matriz J = D−1 (E + F ) = D−1 (D − A) = I − D−1 A se denomina matriz
de Jacobi.
Teorema 2.13 Si A es una matriz de diagonal dominante, el método de Ja-

cobi es convergente.
2.5.2 Método de Gauss-Seidel

Este método es el resultado de realizar la descomposición A = M − N =
(D − E) − F . El sistema nos queda
Ax = b =⇒ (D − E)x = F x + b =⇒ x = (D − E)−1 F x + (D − E)−1 b
La matriz
L1 = (D −E)−1 F = (A+F )−1 (A+F −A) = I −(A+F )−1 A = I −(D −E)−1 A
recibe el nombre de matriz de Gauss-Seidel.
Teorema 2.14 Si A es una matriz de diagonal estrictamente dominante, el

método de Gauss-Seidel es convergente.
2.5.3 Métodos de relajación (SOR)

Este método realiza la descomposición

1 1−ω 1 1−ω
A= D− D − E − F = (D − ωE) − D+F =M −N
ω ω ω ω
El sistema se transforma entonces en

1 1−ω
(D − ωE)x = D + F x + b =⇒
ω ω

(D − ωE)x = (1 − ω)D + ωF x + ωb =⇒

x = (D − ωE)−1 (1 − ω)D + ωF x + ω(D − ωE)−1 b
La matriz del método

−1
Lω = (D − ωE) (1 − ω)D + ωF
recibe el nombre de matriz de relajación.
• Si ω = 1 la matriz se reduce a L1 = (D − E)−1 F , es decir, se trata del

método de Gauss Seidel.
Métodos del descenso más rápido y del gradiente conjugado 65
• Si ω > 1 se dice que se trata de un método de sobre-relajación
• Si ω < 1 se dice que se trata de un método de sub-relajación
Teorema 2.15 Una condición necesaria para que converja el método de rela-
jación es que ω ∈ (0, 2).
Teorema 2.16 Si A es de diagonal dominante, el método de relajación es

convergente cualquiera que sea ω ∈ (0, 1].
Teorema 2.17 Si A es simétrica y definida positiva, el método de relajación

converge si, y sólo si, ω ∈ (0, 2)
2.6 Métodos del descenso más rápido y del

gradiente conjugado
Los métodos que vamos a tratar a continuación son válidos para sistemas
Ax = b cuya matriz A es simétrica y definida positiva, es decir, para matrices
tales que AT = A (simétrica) y xT Ax > 0 cualquiera que sea el vector x 6= 0
(definida positiva).
Lema 2.18 Si A es simétrica y definida positiva, el problema de resolver el

sistema Ax = b es equivalente al de minimizar la forma cuadrática
q(x) = h x, Ax i − 2h x, b i
donde h x, y i = xT y representa el producto escalar de los vectores x e y.
Demostración. Fijemos una dirección v (rayo unidimensional) y vamos a ver

cómo se comporta la forma cuadrática q para vectores de la forma x+tv donde
t es un escalar.
tv -
3

x

x + tv

q(x + tv) = h x + tv, A(x + tv) i − 2h x + tv, b i

= h x, Ax i + 2th x, Av i + t2 h v, Av i − 2h x, b i − 2th v, b i
(2.3)
= q(x) + 2th v, Ax i − 2th v, b i + t2 h v, Av i
= q(x) + 2th v, Ax − b i + t2 h v, Av i
ya que AT = A.
La ecuación (2.3) (ecuación de segundo grado en t con el coeficiente de t2
positivo, tiene un mı́nimo que se calcula igualando a cero la derivada
d
q(x + tv) = 2h v, Ax − b i + 2th v, Av i
dt
es decir, en el punto
t̂ = h v, b − Ax i/h v, Av i.
El valor mı́nimo que toma la forma cuadrática sobre dicho rayo unidimensional
viene dado por
q(x + t̂v) = q(x) + t̂ [2h v, Ax − b i + t̂h v, Av i]
= q(x) + t̂ [2h v, Ax − b i + h v, b − Ax i]
= q(x) − t̂h v, b − Ax i
= q(x) − h v, b − Ax i2 /h v, Av i
Esto nos indica que al pasar de x a x + t̂v siempre hay una reducción en el
valor de q excepto si v ⊥ (b − Ax), es decir, si h v, b − Ax i = 0. Ası́ pues, si
x no es una solución del sistema Ax = b existen muchos vectores v tales que
h v, b − Ax i =
6 0 y, por tanto, x no minimiza a la forma cuadrática q. Por el
contrario, si Ax = b, no existe ningún rayo que emane de x sobre el que q tome
un valor menor que q(x), es decir, x minimiza el valor de q.
El lema anterior nos sugiere un método iterado para resolver el sistema Ax =
b procediendo a minimizar la forma cuadrática q a través de una sucesión de
rayos.
En el paso k del algoritmo se dispondrá de los vectores
x(0) , x(1) , x(2) , . . . , x(k) .
Estos vectores nos permitirán buscar una dirección apropiada v (k) y el siguiente
punto de la sucesión vendrá dado por
x(k+1) = x(k) + tk v (k) ,
donde
h v (k) , b − Ax(k) i
tk =
h v (k) , Av (k) i
Gráficamente, si kv (k) k = 1, tk mide la distancia que nos movemos de x(k) para
obtener x(k+1)
Métodos del descenso más rápido y del gradiente conjugado 67
tk v (k) - (k)
v
-
"3
"
"
x(k)
"
"
"
" x(k+1)
"
"
"
"
"
2.6.1 Método del descenso más rápido

Si tomamos v (k) como el gradiente negativo de q en x(k) , es decir, como la
dirección del residuo r(k) = b − Ax(k) obtenemos el denominado método del
descenso más rápido.
Teniendo en cuenta que los diferentes vectores x(i) no es necesario conservar-
los, los podemos sobreescribir obteniéndose el siguiente algoritmo:
input x, A, b, n
for k = 1, 2, 3 . . . , n do
v ← b − Ax
t ← h v, v i/h v, Av i
x ← x + tv
output k, x
end
Este método resulta, en general, muy lento si las curvas de nivel de la forma
cuadrática están muy próximas, por lo que no suele utilizarse en la forma
descrita.
Sin embargo, utilizando condiciones de ortogonalidad en las denominadas
direcciones conjugadas, puede ser modificado de forma que se convierta en un
método de convergencia rápida que es conocido como método del gradiente
conjugado.
2.6.2 Método del gradiente conjugado

Por no profundizar en el concepto de direcciones conjugadas y en cómo se
determinan, nos limitaremos a dar el algoritmo correspondiente al método.
input x, A, b, n, ε, δ
r ← b − Ax
v←r
c ← h r, r i
for k = 1, 2, 3 . . . , n do
if h v, v i1/2 < δ then stop
z ← Av
t ← c/h v, z i
x ← x + tv
r ← r − tz
d ← h r, r i
if d2 < ε then stop
v ← r + (d/c)v
c←d
output k, x, r
end

Ejercicio
2.1 Estudiar
el número de condición de Frobenius de la matriz
a −b
A= .
a + ε −b
Ejercicio 2.2 Dado el sistema:

x + y = 2
2x + y = 3
a) Calcular su número de condición de Frobenius.
b) Calcular “a” para que el número de condición del sistema resultante de

sumarle a la segunda ecuación la primera multiplicada por dicha cons-
tante “a”, sea mı́nimo.

3x + 4y = 7
3x + 5y = 8
a) Calcular su número de condición euclı́deo.

b) Sustituir la segunda ecuación por una combinación lineal de ambas, de

forma que el número de condición sea mı́nimo.
Ejercicio 2.4 Comprobar que la matriz:

 
1 2 0 0 0
 1 4 3 0 0 
 
A=  0 4 9 4 0 

 0 0 9 16 5 
0 0 0 16 25
admite factorización LU y realizarla.
Ejercicio 2.5 Resolver, por el método de Cholesky, el sistema de ecuaciones:

    
6 −1 + 3i 1 − 2i x1 −1 − 2i
 −1 − 3i 3 −1 + i   x2  =  1 + i 
1 + 2i −1 − i 2 x3 1 − 2i
 
p −p 2p
Ejercicio 2.6 Dada la matriz A =  −p p + 2 −1  se pide:
2p −1 6p − 1
a) Determinar para qué valores de p es hermı́tica y definida positiva.
b) Para p = 1, efectuar la descomposición de Cholesky y utilizarla para

resolver el sistema Ax = b siendo b = (1 0 3)t
Ejercicio 2.7 Resolver por los métodos de Jacobi, Gauss-Seidel y SOR con
ω = 1.2, el sistema:
10x1 − x2 + 2x3 = 6
−x1 + 11x2 − x3 + 3x4 = 25
2x1 − x2 + 10x3 − x4 = −11
3x2 − x3 + 8x4 = 15
Ejercicio 2.8 Al resolver por el método de Gauss-Seidel, utilizando MAT-

LAB, el sistema 
 x − 3y + 5z = 5
8x − y − z = 8
−2x + 4y + z = 4

observamos que el programa se detiene en la iteración 138 dándonos el vector

(inf inf − inf )T .
a) El método de Gauss-Seidel realiza el proceso xn+1 = L1 xn +c. Determina

la matriz L1 .
b) Utilizar los cı́rculos de Gerschgorin para estimar el módulo de los auto-

valores de L1 .
c) Justificar el porqué de la divergencia del método. (Indicación: utilizar

el apartado anterior).
d) ¿Existe alguna condición suficiente que deba cumplir la matriz de un

sistema para garantizar la convergencia del método de Gauss-Seidel?
Hacer uso de ella para modificar el sistema de forma que el proceso sea
convergente?
Ejercicio 2.9 Sea α ∈ {0.5, 1.5, 2.5} y consideremos el sistema iterado

   
1 1
 α −1 1  x  1− α 
   
xn+1 n
 =  
 + 


yn+1  1  y
n
 1 
−1 +1 1−
α α
Se pide
a) Resolver el sistema resultante de tomar lı́mites para probar que, en caso

de que converja, el lı́mite de la sucesión

x0 x1 x2
, , ...
y0 y1 y2
no depende de α.
b) ¿Para qué valores de α converge la sucesión?
c) Para los valores anteriores que hacen que la sucesión sea convergente,
¿con cuál lo hace más rápidamente?

x0 0.5
d) Comenzando con el vector = , aproximar iteradamente
y0 0.5
el lı́mite de la sucesión utilizando el valor de α que acelere más la con-
vergencia.
Ejercicio 2.10 Sea el sistema Ax = b, donde

1000 999 x1 1999
A= , x= y b= .
999 998 x2 1997
a) Obtener la factorización LU de la matriz A. ¿Se puede conseguir la

factorización de Choleski?
b) Resolver el sistema Ax = b utilizando la factorización A = LU obtenida

en el apartado anterior.
c) Calcular ||A||∞ , ||A−1 ||∞ y el número de condición de la matriz κ∞ (A).

¿Se puede decir que está bien condicionada?
d) Comprueba que ||Ax||∞ = ||A||∞ para la solución x = (1, 1)T del sistema
Ax = b.
¿Cuál es el máximo valor que puede tomar ||Ax||∞ , cuando x es un
vector unitario para la norma || ||∞ ?
e) Si se perturba b en b+δb = (19980 99, 19970 01)T, calcular ||δb||∞ /||b||∞ .

Si x + δx es la solución obtenida para el nuevo sistema Ax = b + δb, ¿es
el error relativo ||δx||∞ /||x||∞ el máximo que se puede cometer?
||δx||∞ ||δb||∞
Indicación: ≤ κ∞ (A) .
||x||∞ ||b||∞
3. Sistemas inconsistentes y sis-
temas indeterminados
3.1 Factorizaciones ortogonales

Si tratamos de resolver un sistema Ax = b mediante la factorización LU (o la
de Cholesky), lo que hacemos es transformar el sistema en Ax = LU x = b para
hacer U x = L−1 b que es un sistema triangular que se resuelve por sustitución
regresiva. Sin embargo, la matriz del nuevo sistema es U = L−1 A y dado que
L−1 no es una matriz unitaria (ortogonal en el caso real) el número de condición
de la matriz del sistema ha cambiado pudiendo estar peor condicionada que
la matriz A del sistema original.
Vamos a estudiar, a continuación, otro tipo de factorización A = QR donde
R es, al igual que U , una matriz triangular superior, pero donde Q va a ser
una matriz unitaria, por lo que el sistema Ax = b lo transformaremos en
Rx = Q−1 b = Q∗ b y, a diferencia del caso anterior, R = Q∗ A tiene el mismo
número de condición que la matriz A del sistema original, ya que Q∗ es unitaria.
3.2 Interpretación matricial del método de

Gram-Schmidt: factorización QR
 
a11 a12 · · · a1n
 a21 a22 · · · a2n 
Consideremos la matriz regular A =  .. ..  = (a1 a2 · · · an )
 
.. . .
 . . . . 
an1 an2 · · · ann
donde ai representa la columna i-ésima de la matriz A.
Aplicando Gram-Schmidt existe un sistema ortonormal {y1 , y2 , . . . , yn } tal
73
74 Sistemas inconsistentes y sistemas indeterminados
que L{y1 , y2 , . . . , yk } = L{a1 , a2 , . . . , ak }, por lo que el vector yk+1 pertenece

a la variedad L⊥ {a1 , a2 , . . . , ak }.
Sea Q la matriz cuyas columnas son los vectores yi , Q = (y1 y2 · · · yn ).
Entonces,  
y1∗
 y∗ 
∗  2 
Q A =  ..  (a1 a2 · · · an )
 . 
yn∗
es decir:  
h a1 , y 1 i h a2 , y 1 i · · · h an , y 1 i
 h a1 , y 2 i h a2 , y 2 i · · · h an , y 2 i 
Q∗ A = 
 
.. .. ... .. 
 . . . 
h a1 , y n i h a2 , y n i · · · h an , y n i
Como yk+1 ∈ L⊥ {a1 , a2 , . . . , ak }, se tiene que h ai , yj i = 0 si, y sólo si, i < j,
por lo que la matriz Q∗ A es una triangular superior.
 
r11 r12 r13 · · · r1n
 0 r22 r23 · · · r2n 
 
Q∗ A = R =  0
 0 r33 · · · r3n 

 .. .. .. . . .
. ..

 . . . 
0 0 0 · · · rnn
Como las columnas de Q constituyen un sistema ortonormal de vectores, Q

es unitaria, es decir Q∗ Q = I, por lo que A = QR.
El problema que plantea la descomposición QR es que la matriz Q no es otra
que la constituida por una base ortonormal obtenida a partir de las columnas
de A por el método de Gram-Schmidt. Las transformaciones que se realizan
para ortonormalizar los vectores columna de la matriz A mediante el método
de Gram-Schmidt son transformaciones no unitarias, por lo que aunque el
resultado sea una factorización ortogonal, los pasos que se han dado para
llegar a ella han sido transformaciones no ortogonales. Ello nos lleva a tratar
de buscar un método por el que podamos realizar una factorización QR en la
que todos los pasos que se realicen sean transformaciones ortogonales.
Rotaciones y reflexiones 75
3.3 Rotaciones y reflexiones

Consideremos la matriz
 
a11 · · · a1i · · · a1j · · · a1n
 .. . . .. .. .. 
 . . . . . 
 ai1 · · · aii · · · aij · · · ain 
 
 . .. . . . .. 
 ..
A= . . .. . 
 j1 · · · aji · · · ajj · · · ajn 
 a

 . .. .. . . .. 
 .. . . . . 
an1 · · · ani · · · anj · · · ann
y tratemos de anular el elemento aji 6= 0. Para ello vamos a aplicar una
rotación  
1 ··· 0 ··· 0 ··· 0
 .. . . .. .. .. 
 . . . . . 
 0 · · · cos α · · · sen α · · · 0 
 
 . .. ... .. .. 
 ..
Qji =  . . . 

 0 · · · − sen α · · · cos α · · · 0 
 
 . .. .. . . .. 
 .. . . . . 
0 ··· 0 ··· 0 ··· 1
La matriz Qji A nos queda
 0
a11 · · · a01i · · · a01j · · · a01n

 .. . . . .. .. 
 . . .. . . 
 0
 ai1 · · · a0ii · · · aij · · · a0in 
0 
 . .. . . .. .. 
 ..
Qji A =  . . . . 

 a0 · · · a0 · · · a0jj · · · a0jn 
 j1 ji 
 . .. .. . . .. 
 .. . . . . 
an1 · · · a0ni · · ·
0
a0nj · · · a0nn
con a0ji = −aii sen α + aji cos α, por lo que
• Si aii = 0 se tiene que a0ji = aji cos α y si queremos anularlo aji cos α = 0.
Dado que suponemos que aji 6= 0, basta con hacer cos α = 0, es decir:
α = π/2.
• Si aii 6= 0 tendremos que hacer −aii sen α + aji cos α = 0, por lo que
aji
t = tg α = y, por tanto,
aii
t 1
sen α = √ cos α = √
1 + t2 1 + t2
Obsérvese además que los únicos elementos que se alteran en la matriz Qji A
son los correspondientes a la filas y columnas i y j.
Podemos, por tanto, mediante rotaciones, anular todos los elementos subdia-
gonales y llegar a una matriz R triangular superior
Qk · · · Q2 Q1 A = R ⇐⇒ Q∗ A = R con Q∗ = Qk · · · Q2 Q1
Dado que las matrices Qi de las rotaciones son ortogonales, su producto

también los es, por lo que Q∗ es una matriz ortogonal y, por tanto, A = QR.
En éste método de factorización QR, a diferencia del aplicado anteriormente
mediante el método de Gram-Schmidt todos los pasos que se dan están asocia-
dos a transformaciones ortogonales, sin embargo resulta costoso ya que cada
rotación hace un único cero en la matriz A, por lo que para una matriz de
n2 − n
orden n serı́an necesarias rotaciones.
2
Otra posibilidad es hacer reflexiones en vez de rotaciones, ya que éstas con-
siguen anular todos los elementos situados por debajo de uno prefijado de una
determinada columna. Este tipo de transformaciones vamos a estudiarlas a
continuación con las denominadas transformaciones de Householder.
3.4 Transformaciones de Householder

Consideremos un espacio vectorial de dimensión n definido sobre un cuerpo
K, que denotaremos por Kn (en general trabajaremos en Rn ó Cn ). Dado
un vector v ∈ Kn se define la transformación H de Householder asociada al
vector v a la que viene definida por la matriz:
n×n

 I∈K
 si v = 0
H=
 I − 2 vv ∗ si v 6= 0

v∗v
Proposicin 3.1 La transformación H de Householder asociada a un vector

v ∈ Kn posee las siguientes propiedades:
a) H es hermı́tica (H ∗ = H).
b) H es unitaria (H ∗ H = HH ∗ = I).
c) H 2 = I o lo que es lo mismo, H −1 = H.
Transformaciones de Householder 77
Demostración.
∗
∗ 2 ∗ ∗ 2 2
a) H = I − ∗ vv =I − ∗
(vv ∗ )∗ = I − ∗ vv ∗ = H
v v v v v v
Obsérvese que v ∗ v = h v, v i = kvk2 ∈ R, por lo que v ∗ v = v ∗ v
2
∗ 2 ∗ 2 ∗ 4 ∗ 2
b) HH = HH = I − ∗ vv I − ∗ vv = I− ∗ vv + ∗ vv ∗ vv ∗ =
v v v v v v v v
4 4 4 4
= I− ∗ vv ∗+ ∗ 2 v(v ∗ v)v ∗ = I− ∗ vv ∗+ ∗ 2 (v ∗ v)vv ∗ = I.
v v (v v) v v (v v)
c) Basta observar que, por los apartados anteriores, H 2 = HH = HH ∗ = I.
3.4.1 Interpretación geométrica en Rn

Sean v ∈ Rn un vector tal que kvk2 = 1 y H la transformación de Householder
asociada a él:
H = I − 2vv T
Dado un vector x ∈ Rn se tiene que
Hx = I − 2vv T x = x − 2vv T x = x − 2vh x, v i =
x
3

= x − 2v(kxk kvk cos α) = x − 2v(kxk cos α) =
= x − 2λv
α v
con λ = kxk cos α, donde α representa el ángulo -
λ -
que forman los vectores x y v.
Sea y el vector simétrico de x respecto del hiperplano perpendicular a v.

Podemos observar que y + λv = x − λv, o lo que es lo mismo, que y =
x − 2λv = Hx. Es decir, H transforma a un vector x en su simétrico respecto
del hiperplano perpendicular al vector v.
k
Q
6Q y x 36

y + λv Q
Q

x − λv
Q
−v Q
Q
v
Q-
-
λv −λv
Figura 3.1: x y su transformado.
Es fácil observar que si x es ortogonal a v se verifica que Hx = x, ası́ como

que Hv = −v.
Ası́ pues, si x e y son dos vectores de Rn tales que

x 6= y con kxk = kyk, la transformación de Hou- 6
x−y x+y
seholder asociada al vector v = transforma
kx − yk
el vector x en y, es decir, Hx = y.
En efecto: dado que ambos vectores tienen la y x
misma norma, h x + y, x − y i = kxk2 − h x, y i + KA
h y, x i − kyk2 = 0. Además, los vectores x e y son A

simétricos respecto de la dirección del vector x + y,
A
A
por lo que la transformación de Householder aso- A
x−y A x−y
ciada al vector v = transforma a x en y. A -
kx − yk

 x = (x1 , x2 , . . . , xn )

Consideremos los vectores q
 y = (x1 , x2 , . . . , xk , x2k+1 + · · · + x2n , 0, . . . , 0)

que poseen la misma norma.

0
 
..
.
 
 
0
 
 
x−y 1  q 
 xk+1 − x2 + · · · + x2  la transformación
Si v = = k+1 n
kx − yk kx − yk  
xk+2
 
 

 .
..


xn
H de Householder asociada a v transforma a x en y.
3.4.2 Householder en Cn
Sean v un vector de Cn y H la transformación de Householder asociada a él:
2
H=I− vv ∗
v∗v
• Si x = λv con λ ∈ C entonces

2 2λ
Hx = I − ∗ vv λv = λv − ∗ vv ∗ v = −λv = −x
∗
v v v v
• Si x ⊥ v se verifica que h x, v i = v ∗ x = 0 y, por tanto

2 2
Hx = I − ∗ vv x = x − ∗ vv ∗ x = x
∗
v v v v
Es decir, los vectores ortogonales a v son invariantes mediante H.
Cualquier vector x ∈ Cn puede ser descompuesto de forma única en la suma

de uno proporcional a v y otro w perteneciente a la variedad W ortogonal a
v, es decir x = λv + w con w ⊥ v. Ası́ pues,
Hx = H(λv + w) = H(λv) + Hw = −λv + w
por lo que Hx es el vector simétrico de x respecto del hiperplano ortogonal a

v.
Si x = (x1 , x2 , . . . , xn ) ∈ Cn y pretendemos encontrar un vector v tal que la

transformación de Householder Hv asociada a v transforme dicho vector x en
otro y = (y1 , 0, . . . , 0) es evidente que como
kyk2 = kHv xk2 = kxk2
por ser Hv unitaria, ambos vectores x e y han de tener igual norma, es decir,
ha de verificarse que |y1 | = kxk2 o lo que es lo mismo, y1 = kxk2 eiα con α ∈ R.
Tomemos un vector v unitario, es decir, tal que |x1 |2 + |x2 |2 + · · · + |xn |2 = 1.
Entonces
Hv x = (I − 2vv ∗ ) x = x − 2vv ∗ x = x − (2v ∗ x)v = y
Obligando a que 2v ∗ x = 1 se tiene que x − v = y, o lo que es lo mismo, que

v = x − y, es decir:    
v1 x1 − y1
 v2   x2 
v =  ..  = 
   
.. 
 .   . 
vn xn
De este vector son conocidas todas sus componentes excepto la primera v1 =
x1 − y1 .
Sabemos que 2v ∗ x = 1 y que v ∗ v = 1, por lo que
2(v1 x1 + |x2 |2 + · · · + |xn |2 ) = 1

) )
2(v1 x1 + v2 x2 + · · · + vn xn ) = 1
⇒ ⇒
v1 v1 + v2 v2 + · · · + vn vn = 1 |v1 |2 + |x2 |2 + · · · + |xn |2 = 1
2(v1 x1 + kxk2 − |x1 |2 ) = 1 (3.1)

|v1 |2 + kxk2 − |x1 |2 = 1 (3.2)
De la ecuación 3.1 se deduce que v1 x1 es un número real, por lo que el

argumento del producto v1 x1 ha de ser 0 ó π. Como, por otra parte, v1 x1 =
∧
|v1 | |x1 | ei(v1 x1 ) , los complejos v1 y x1 han de tener igual argumento, por lo que
v1 = λx1 .
Llevando este resultado a las ecuaciones 3.1 y 3.2 se obtiene que
2(λ |x1 |2 + kxk2 − |x1 |2 ) = 1

=⇒ |x1 |2 (λ2 − 2λ + 1) = kxk2 =⇒
λ2 |x1 |2 + kxk2 − |x1 |2 = 1
kxk
λ=1± supuesto x1 6= 0
|x1 |
(Si x1 = 0 basta con tomar v = (kxk , x2 , . . . , xn )).

kxk

1± x1

 |x1 | 

El vector que buscamos es, por tanto, v = 
 x2  obteniéndose

 .. 
 . 
xn
que  
y1
 0 
 
iα
 ..  con y1 = kxk e
Hv x =  
 . 
0
que resulta ser  x1 
∓ kxk
 |x1 | 
0
 
Hv x = y = x − v = 
 
.. 
.
 
 
0
3.4.3 Factorización QR de Householder

 
a11 a12 · · · a1n
 a21 a22 · · · a2n 
Supongamos que tenemos el sistema Ax = b con A =  .. .. .
 
.. ..
 . . . . 
   p   
a11 a211 + · · · + a2n1 r11
 a21   0   0 
Sean x1 =   e y1 =  = .
     
.. .. ..
 .   .   . 
an1 0 0
x 1 − y1
Sea H1 la transformación de Householder asociada al vector v1 = ,
kx1 − y1 k
por lo que H1 x1 = y1 . Tenemos entonces que
 
(1) (1)
r11 a12 · · · a1n
 
(1) (1) 
 0 a22 · · · a2n  (1) (1)

(1)
H1 A =  .. .. . . . 
= a
1 a 2 · · · a n
 . . . .. 
 
(1) (1)
0 an2 · · · ann
(1)
en la que ai representa la columna i-ésima de la matriz H1 A.
Busquemos ahora otro vector v2 tal que la transformación de Householder
(1) (1)
H2 asociada a él, deje invariante al vector a1 y transforme al vector a2 en
otro de la forma (r12 , r22 , 0, . . . , 0).
(1)
Como se quiere que deje invariante al vector a1 , v2 ha de ser ortogonal a él,
es decir, debe ser de la forma (0, u2 , . . . , un ).
   
(1) (1)
a21 a21
   q 
 (1)  (1) 2 (1) 2 
 a22   (a22 ) + · · · + (an2 ) 

   
(1)  (1) 
Tomando x2 = a2 =  a32  e y2 =  , la
 
0
 .. 
   
 .. 
 .   . 
   
(1)
an2 0
x2 − y2
transformación H2 asociada al vector v2 = deja invariante al vector
kx2 − y2 k
(1) (1)
a1 y transforma x2 = a2 en y2 .
Reiterando al procedimiento se puede triangularizar la matriz A llegar a una
matriz triangular superior R. Llegados a ese paso se tiene que
Hk Hk−1 · · · H1 A = R ⇐⇒ Q∗ A = R con Q∗ = Hk Hk−1 · · · H1
de donde
A = QR.
Si lo que nos interesa es resolver el sistema aplicamos las transformaciones

al sistema y no sólo a la matriz A.
Hk Hk−1 · · · H1 Ax = Hk Hk−1 · · · H1 b ⇐⇒ Rx = b0
sistema triangular de fácil resolución.

 
1 −1 −1
 
Consideremos, por ejemplo, la matriz A =  1 .
 
2 0
 
−2 7 1
   p   
1 12 + 22 + (−2)2 3
     
Como x1 =  2  e y1 =   =  0 , se tiene que
     
0
     
−2 0 0
   
1
−2 − 3
√
x1 − y1 1    
v1 = = √  2 = √  1
   
kx1 − y1 k 2 3   3 
−2 − √13
 
− √13
2 ∗
 
H1 = I − v v = I − 2 √1  − √13 √1 − √13 =
 
∗ 1 1
v1 v1
 3 3
 
1
− 3
√
     
1
1 0 0 3
− 13 1
3
1
3
2
3
− 23
     
−
 1 1
= 0 2  −3 − 3  =  23
1 1 2
    
1 0  3 3 3

     
1 1 1
0 0 1 3
−3 3
− 23 2
3
1
3
    
1 2
3 3
− 23 1 −1 −1 3 −5 − 31
    
H1 A =  2 1 2 1
1 = 0
    
3 3 3
 2 0 4 3

    
− 23 2
3
1
3
−2 7 1 0 3 5
3
     
−5 −5 −5
   √   
Ahora son x2 =  4  e y2 =  42 + 32 = 5 , por lo que
     
     
3 0 0
 
0
x2 − y2 1  
v2 = = √  −1 
 
kx2 − y2 k 10  
3
   
0 1 0 0
2 ∗ 2    
H2 = I − ∗ v2 v2 = I −  −1  0 −1 3 =  0 4 3
   
v2 v2 10  5 5

  
3
3 0 5
− 45
    
1 0 0 3 −5 − 13 3 −5 − 13
    
H2 H1 A =  0
 4 3 
  1 = 0
  19 
5 5 
0 4 3
5 15

   
3
0 5
− 45 0 3 5
3
0 0 − 17
15
 
a11 a12 · · · a1n
 a21 a22 · · · a2n 
Si partimos de una matriz A =  .  y tomamos
 
.. .. ..
 . . . .. 
T
x 1 = a1 = a11 a21 · · · an1
 
kxk
 1± |a11 |
a11 
 
definimos el vector v1 = 
 a21 , obteniéndose que

 .. 
 . 
an1
 
(1) (1)
α11 a12 ··· a1n
 
(1) (1)
· · · a2n
 
 0 a21 
H1 A =  .. .. .. .

. ..
 
 . . 
 
(1) (1)
0 an2 · · · ann
   
(1)
y a12
   1 
 (1)  
 a22   y2 

   
Buscamos ahora otro vector v2 tal que H2  a(1)   0  y de tal forma
=
   
 32  
 ..   .. 

 .   . 
   
(1)
an2 0
T
que mantenga invariante al vector α11 0 · · · 0 . Bastará coger, para
T
ello, v2 = 0 v2 · · · vn , que es ortogonal a él.
En este caso, la transformación de Householder viene dada por
   
0 0 0 ··· 0
   
· · · v2 vn 
   
 v2   0 v2 v2
H2 = I − 2  .  0 v 2 · · · vn = I − 2 
 
 .. .. ... ..  =

 ..   . . . 
   
vn 0 vn v2 · · · vn vn
   
1 0 ··· 0 1 0 ··· 0
   
0 1 − 2v2 v2 · · · −2v2 vn
   
   0 
= .. .. ... ..
=
..

H
   
 . . .   . 
   
0 −2vn v2 · · · 1 − 2vn vn 0
Aplicando a la matriz A ambas transformaciones, se tiene:

 
(1) (1)
α11 a12 ··· a1n
 
(1) (1)
···
 
 0 a21 a2n 
H2 H1 A = H2  .. .. .. ..
=
.
 
 . . . 
 
(1) (1)
0 an2 · · · ann
    
(1) (1) (1) (1)
1 0 ··· 0 1 a12 · · · a1n 1 a12 ··· a1n
    
    
0  0  0 
=
 ..   ..
 =
  ..

H A1 HA1

. .  . 
    
0 0 0
Sistemas superdeterminados. Problema de los mı́nimos cuadrados 85
Es decir, se trata de realizar un proceso análogo al anterior sólo que ahora

en Cn−1 . Posteriormente se realizará otro en Cn−2 y ası́ sucesivamente hasta
haber triangularizado la matriz A.
3.5 Sistemas superdeterminados. Problema de

los mı́nimos cuadrados
Dado un sistema de ecuaciones de la forma Ax = b en el A es una matriz
cuadrada de orden n, A ∈ Kn×n , y x y b son vectores de Kn sabemos,
por el teorema de Rouche-Fröbenius , que tiene solución si, y sólo si, exis-
ten x1 , x2 , . . . , xn ∈ Kn tales que
          
a11 · · · a1n x1 b a a1n b1
  1   11
.. ..  ..   ..   .. ..   .. 
      
..
. . = .  ⇔ x1  . + · · · +x . =
  
 . .   n . 
          
an1 · · · ann xn bn an1 ann bn
En otras palabras, el vector b puede expresarse como una combinación lineal

de las columnas de la matriz A, por lo que b ∈ C(A) (espacio columna de A).
Sin embargo, existen problemas en los que no ocurre ası́. Supongamos que
se tienen tres puntos en el plano y se desea calcular la recta que pasa por
ellos. Evidentemente, y dado que una recta la determinan sólo dos puntos, el
problema no tiene solución (salvo que los tres puntos estén alineados). Desde
el punto de vista algebraico este problema se expresa de la siguiente forma:
sean P = (a1 , b1 ), Q = (a2 , b2 ) y R = (a3 , b3 ). Si tratamos de hallar la ecuación
de la recta y = mx + n que pasa por ellos se obtiene
   
ma1 + n = b1 a1 1   b1
  m  
⇐⇒ =
   
ma2 + n = b2  a2 1    b2 
  n  
ma3 + n = b3 a3 1 b3
Decir
 queel sistema no posee solución equivale a decir que el vector

b a 1
 1   1 
b =  b2  no pertenece al espacio columna de la matriz A =  a2 1 .
   
   
b3 a3 1
Se define un sistema superdeterminado como aquel sistema de ecuaciones
lineales Ax = b en el que A ∈ Km×n , x ∈ Kn y b ∈ Km , donde K = R ó C.
Supongamos que se tiene un sistema superdeterminado

   
a · · · a1n b1
 11
 .. ..   .. 
   
..
 . . .  x1  . 
..
    
· · · ann    =  bn 
    
 an1 .
 .. ..   .. 
    
 . .  xn  . 
   
am1 · · · amn bm
con m > n, en el que rg A = n es decir, en el

que la matriz del sistema tiene rango máximo,
6 b
y denotemos por a1 , a2 , . . . , an las columnas de
c
A.

Si el sistema es incompatible es debido a que

el vector b no pertenece al espacio columna de -
C(A) b̃
A. Tomando cualquier vector b ∈ C(A) se sabe
que el sistema Ax = b posee solución única.
De entre todos los vectores del espacio columna de A se trata de buscar

aquel que minimiza su distancia al vector b, es decir, aquel vector b̃ ∈ C(A)
tal que b − b̃ es mı́nima (problema de los mı́nimos cuadrados). Dicho vector
sabemos que es la proyección ortogonal de b sobre el espacio C(A) y, que
respecto de la base formada por las columnas ai (1 ≤ i ≤ n) de la matriz A,
tiene por coordenadas
b̃ = (h b, a1 i, h b, a2 i, . . . , h b, an i}
Dado que b 6∈ C(A) y b̃ ∈ C(A) (b̃ proyección ortogonal de b sobre C(A)),

podemos expresar b como suma de b̃ más otro vector c de la variedad ortogonal
a C(A) y, además, de forma única. Entonces:
h b, ai i = h b̃ + c, ai i = h b̃, ai i + h c, ai i = h b̃, ai i 1≤i≤n
El sistema Ax = b̃ posee solución única es decir, existen (α1 , α2 , . . . , αn )

únicos, tales que
 
α1
 . 
 
α1 a1 + α2 a2 + · · · + αn an = b̃ ⇐⇒ A  ..  = b̃
 
αn
Multiplicando esta ecación por a1 , a2 , . . ., an , obtenemos
α1 h a1 , a1 i + · · · + αn h a1 , an i = h b̃, a1 i = h b, a1 i
.................................................
.................................................
α1 h an , a1 i + · · · + αn h an , an i = h b̃, an i = h b, an i
que equivale a
    
a∗1 a1 ··· a∗1 an α a∗1 b
 1  
.. ..   ..   .. 
 
...
 .  = 

 . . . 
    
a∗n a1 ∗
· · · an an αn ∗
an b
es decir      
a∗1 α a∗1
 1  
 ..   ..   .. 
  
 .  a1 · · · an  .  =  . b
     
∗ ∗
an αn an
o, lo que es lo mismo:  
α1
 . 
 
A∗ A  ..  = A∗ b
 
αn
Ası́ pues, la solución del sistema A∗ Ax = A∗ b nos proporciona las coordena-

das, respecto de la base formada por las columnas de la matriz A, del vector b̃
proyección ortogonal de b sobre el espacio columna C(A). Estas coordenadas
constituyen lo que denominaremos seudosolución del sistema superdetermi-
nado Ax = b.
3.5.1 Transformaciones en sistemas superdeterminados

Sabemos que dado un sistema compatible Ax = b y mediante transformaciones
elementales puede obtenerse otro sistema BAx = Bb equivalente al anterior,
es decir, obtenemos un sistema que posee la misma (o las mismas) soluciones
que el sistema dado.
Si partimos de un sistema superdeterminado y realizamos, al igual que antes,
transformaciones elementales, puede que el sistema obtenido no posea la misma
seudosolución que el sistema dado.
Obsérvese que para que los sistemas superdeterminados Ax = b y BAx = Bb

posean la misma seudosolución, han de tener igual solución (han de ser equi-
valentes) los sistemas A∗ Ax = A∗ b y (BA)∗ BAx = (BA)∗ Bb. Dado que este
último puede expresarse de la forma A∗ (B ∗ B)Ax = A∗ (B ∗ B)b, sólo podremos
garantizar que ambos sistemas son equivalentes si B ∗ B = I ya que, en dicho
caso, ambos sistemas son el mismo. Es decir, las únicas transformaciones que
podemos garantizar que no alterarán la solución de un sistema superdetermi-
nado son las unitarias.
Dado que las transformaciones de Householder son unitarias, podemos utili-
zarlas para resolver sistemas superdeterminados.
Consideremos el sistema superdeterminado Ax = b (en el que suponemos A
de rango máximo). Mediante transformaciones de Householder H1 , H2 , . . . , Hn
podemos transformar la matriz A en otra de la forma
 
t t · · · t1n
 11 12 
 0 t22 · · · t2n 
 
 .. .. . . ..  
 
 . . . .  
  T
HA = Hn · · · H1 A =  0 0 · · · tnn  = 
  
  Θ
 0 0 ··· 0 
 
 .. .. .. 
 
 . . . 
 
0 0 ··· 0
La seudosolución de este sistema superdeterminado es la solución del sistema

 
T
(HA)∗ (HA)x = (HA)∗ Hb ⇐⇒ T∗ Θ   x = T ∗ Θ Hb
Θ
 
b01
.. 
 
. 
  

  T
o llamando Hb = b0 =  0  , ∗ x = T Θ b0 .
∗

bn  T Θ  

 .. 
 Θ
 . 
 
0
bm
 
b01
.. 
   
.  b01


∗  .. 
   
Es fácil comprobar que T ∗ 0  = T  . , por lo que el

Θ  bn 
.. 
   
b0n

 . 
 
0
bm
cálculo de la seudosolución del sistema superdeterminado Ax = b se hace
resolviendo el sistema  
b01
∗  .. 
 
∗
T Tx = T  . 
 
b0n
y dado que estamos suponiendo que A tiene rango máximo, la matriz T posee
inversa y por tanto T ∗ , por lo que la solución es la misma que la del sistema
triangular  
b0
 1 
 . 
T x =  ..  ⇐⇒ T x = b̃
 
0
bn
Una vez calculada la seudosolución, la norma del error está representada por
la distancia b − b̃ que viene dada por
   
 b01
 b01
b0  ..  
  
.. 

 1  .   . 
 ..    
   .    
 0  
 b0n+1
0
T

   bn   bn  
.

 b0n  = 
x −  −  ..
 = 
  

   0 
Θ  .. 
   0   bn+1   
 .   ..  
  
.. 
 b0m
   .   . 
b0m    
0
0 bm
Por último, si la matriz A no tiene rango máximo, sus columnas no son

linealmente independientes, por lo que sólo constituyen un sistema generador
(no una base) del espacio columna C(A). Ello nos lleva a la existencia de
infinitas n-uplas (α1 , α2 , . . . , αn ) soluciones del sistema Ax = b̃ y, por tanto,
a infinitas pseudosoluciones del sistema superdeterminado, pero teniendo en
cuenta que al ser única la proyección ortogonal b̃ de b sobre el espacio columna
C(A), todas ellas representan diferentes coordenadas del vector b̃ respecto del
sistema generador de C(A) dado por las columnas de A. Sin embargo, el error
cometido b − b̃ es el mismo para todas las seudosoluciones del sistema.
3.6 Descomposición en valores singulares y seu-

doinversa de Penrose
La descomposición en valores singulares es otra factorización matricial que
tiene muchas aplicaciones.
Teorema 3.1 Toda matriz compleja A, de orden m × n puede ser factorizada

de la forma A = P DQ donde P es una matriz unitaria m × m, D una matriz
diagonal m × n y Q una unitaria de orden n × n.
Demostración. La matriz A∗ A es hermı́tica de orden n × n y semidefinida

positiva, ya que
x∗ (A∗ A)x = (Ax)∗ (Ax) ≥ 0
Resulta, de ello, que sus autovalores son reales no negativos σ12 , σ22 , . . . , σn2
(pudiendo estar repetidos, pero ordenados de forma que los r primeros son no
nulos y los n − r últimos son nulos). Los valores σ1 , σ2 , . . . , σn son los valores
singulares de la matriz A.
Sea {u1 , u2 , . . . , un } un conjunto ortonormal de vectores propios de A∗ A,
dispuestos de forma que
A∗ Aui = σi2 ui
Se verifica entonces que
kAui k22 = u∗i A∗ Aui = u∗i σi2 ui = σi2
Esto nos muestra que Aui = 0 si i ≥ r + 1. Obsérvese que
r = rg (A∗ A) ≤ mı́n{rg (A∗ ), rg (A)} ≤ mı́n{m, n}
Construyamos la matriz Q de orden n × n cuyas filas son u∗1 , u∗2 , . . . , u∗n y

definamos
vi = σi−1 Aui 1≤i≤r
Los vectores vi constituyen un sistema ortonormal, ya que para 1 ≤ i, j ≤ r
se tiene que
vi∗ vj = σi−1 (Aui )∗ σj−1 (Auj ) = (σi σj )−1 (u∗i A∗ Auj ) = (σi σj )−1 (u∗i σj2 uj ) = δij
Descomposición en valores singulares y seudoinversa de Penrose 91
Eligiendo vectores adicionales vn+1 , vn+2 , . . . , vm de tal forma que {v1 , . . . , vm }

constituya una base ortonormal de Cm y construyendo las matrices P de orden
m × m cuyas columnas son los vectores vi y la matriz D de orden m × n cuyos
elementos diagonales dii = σi y los restantes elementos nulos, se tiene que
A = P DQ
Para probarlo vamos a ver que P ∗ AQ∗ = D. En efecto:
(P ∗ AQ∗ )ij = vi∗ Auj = vi∗ σj vj = σj vi∗ vj = σj δij = Dij
3.6.1 Seudoinversa de Penrose

Para las matrices D de orden m × n tales que dij = 0 si i 6= j y dii > 0, se
define la seudoinversa como la matriz n × m D+ cuyos elementos diagonales
son los inversos de los elementos diagonales de D y el resto de los elementos
son nulos.
En el caso general de una matriz A de orden m × n se define la seudoinversa
A+ a través de la factorización en valores singulares A = P DQ de la forma
A+ = Q∗ D+ P ∗
La seudoinversa comparte algunas propiedades con la inversa, pero sólo al-
gunas ya que, por ejemplo, si A es de orden m × n con m > n A+ es de orden
n × m y la matriz AA+ no puede ser la matriz unidad, ya que AA+ es de orden
m × m, por lo que si fuese la matriz unidad serı́a rg (AA∗ ) = m cosa que no
es posible por ser n < m el máximo rango que pueden tener las matrices A y
A+ (recuérdese que el rango de la matriz producto nunca es superior al rango
de las matrices que se multiplican).
Teorema 3.2 [Propiedades de Penrose] Para cada matriz A existe, a lo

más, una matriz X que verifica las siguientes propiedades:
a) AXA = A
b) XAX = X
c) (AX)∗ = AX
d) (XA)∗ = XA
Demostración. Sean X e Y dos matrices que verifique las cuatro propiedades.

Se tiene entonces que
X = XAX (b)
= XAY AX (a)
= XAY AY AY AX (a)
= (XA)∗ (Y A)∗ Y (AY )∗ (AX)∗ (d) y (c)
= A ∗ X ∗ A ∗ Y ∗ Y Y ∗ A∗ X ∗ A ∗
= (AXA)∗ Y ∗ Y Y ∗ (AXA)∗
= A ∗ Y ∗ Y Y ∗ A∗ (a)
= (Y A)∗ Y (AY )∗
= Y AY AY (d) y (c)
= Y AY (b)
= Y (b)
Teorema 3.3 La seudoinversa de una matriz tiene las cuatro propiedades de

Penrose y, por tanto, es única.
Demostración. Sea A = P DQ la descomposición en valores singulares de

una matriz A. Sabemos que A+ = Q∗ D+ P ∗ .
Si la matriz A es de orden m × n y tiene rango r, la matriz D es también del
mismo orden y tiene la forma

 σ si i = j ≤ r
i
Dij =
 0 en caso contrario
Se tiene entonces que

DD∗ D = D
ya que
n
X m
X
(DD∗ D)ij = Dik +
Dkl Dlj .
k=1 l=1
Los términos Dik y Dlj hacen que el segundo miembro de la igualdad sea nulo
excepto en los casos en que i, j ≤ r en cuyo caso
r
X r
X r
X
∗
(DD D)ij = Dik +
Dkl Dlj = σi Dil+ Dlj = σi σi−1 Dij = Dij .
k=1 l=1 l=1
Razonamientos análogos nos permiten probar que D+ verifica las otras tres
propiedades de Penrose.
Si nos fijamos ahora en A+ , se tiene que
AA+ A = P DQQ∗ D+ P ∗ P DQ = P DD+ DQ = P DQ = A
y razonamientos similares nos permiten probar las otras tres propiedades.
Si la matriz A es de rango máximo, la matriz A∗ A es invertible. Las ecua-

ciones normales del sistema Ax = b vienen dadas por A∗ Ax = A∗ b y dado que
A∗ A es invertible se tiene que
x = (A∗ A)−1 A∗ b. (3.3)
Por otra parte, si hubiésemos resuelto el sistema a través de la seudoinversa

de Penrose habrı́amos obtenido
x = A+ b. (3.4)
por lo que comparando las ecuaciones (3.3) y (3.4) obtenemos, teniendo en

cuenta la unicidad de la seudoinversa, que “si A es de rango máximo”, la
seudoinversa viene dada por
A+ = (A∗ A)−1 A∗ .

4x + 5y = 13
3x + 5y = 11
a) Realizar la factorización QR de la matriz, y resolverlo basándose en ella
a.1) Mediante el método de Gram-Schmidt,

a.2) Mediante transformaciones de Householder.
b) Calcular el número de condición euclı́deo del sistema inicial y del trans-

formado, comprobando que son iguales.
Ejercicio 3.2 Resolver por el método de Householder el sistema:

    
1 −1 −1 x 0
    
1  y  =  4 
    
 2 0
    
−2 7 1 z −7
Ejercicio 3.3 Buscar la solución de mı́nimos cuadrados del sistema Ax = b,

siendo:    
3 −1 0
   
A= 4 y b= 2 
   
2 
   
0 1 1
a) A través de sus ecuaciones normales.
b) Por el método de Householder.
Ejercicio 3.4 Se considera el sistema de ecuaciones Ax = b con

   
1 2 3
   
   
 1 0   2 
A= 
 y b =  .
  
 1 1   0 
   
1 1 1
Se pide:
a) Multiplicando el sistema por la traspuesta de A, calcular la pseudoso-

lución utilizando el método de Choleski.
b) Sea v = (−1, 1, 1, 1)T . Demostrar que la transformación de Householder

asociada al vector v transforma la primera columna de la matriz A en el
vector (2, 0, 0, 0)T dejando invariante la segunda columna de A ası́ como
al vector b.
c) Calcular la pseudosolución del sistema utilizando transformaciones de

Householder, ası́ como la norma del error.
d) Si la matriz A del sistema fuese cuadrada y su número de condición fuese

mayor que 1, ¿qué ventajas e inconvenientes tendrı́a el resolver el sistema
multiplicando por la traspuesta de A y el resolverlo por transformaciones
de Householder?
Ejercicio 3.5 Hallar la recta de regresión de los puntos:
(10 1, 5), (1, 50 1), (2, 70 3), (10 8, 60 9), (10 5, 60 1), (3, 80 8), (30 1, 9) y (20 9, 90 1)
Ejercicio 3.6 Hallar la parábola de regresión de los puntos:
(1, 0), (0, 0), (−1, 0), (1, 2) y (2, 3)
Ejercicio 3.7 Dado el sistema superdeterminado:

   
1 1 0   1
  x  
    
 1 0 1    2 
  y  =  
   
 1 1 1    0 
  z  
1 2 1 −1
calcular, mediante transformaciones de Householder, la solución en mı́nimos

cuadrados (pseudosolución) ası́ como la norma del error.
Ejercicio 3.8 Resolver el sistema

   
2 1   1
  x  
0   =  1
   
 2 
  y  
−1 2 −5
y obtener la norma del error:
a) Mediante sus ecuaciones normales.
b) Mediante transformaciones de Householder.
c) Hallando la inversa generalizada de la matriz del sistema.
Ejercicio 3.9 Se considera el sistema superdeterminado Ax = b con

   
1 7 15 7
   
   
 1 4 8   7 
A= 
 y b=
  
 −5 

 1 0 1 
   
1 3 6 −9
a) Resolverlo mediante transformaciones de Householder, dando la norma

del vector error.
b) Hallar la inversa generalizada A+ de la matriz A.
c) Utilizar la inversa generalizada para resolver el sistema y hallar la norma

del vector error.
Ejercicio 3.10 Resolver el sistema superdeterminado

   
−3 1 1   8
  x  
 1 −3
    
1    4 
 
 y = 

1 −3 
   0 
 1  
  z  
1 1 1 4
calculando la inversa generalizada de la matriz A.
Ejercicio 3.11 Dado sistema superdeterminado Ax = b con

   
1 5 5 7
   
   
 1 2 3   16 
A=   y b= 
 −3 
 
 1 1 3 
   
1 2 1 10
se pide:
a) Resolverlo mediante transformaciones de Householder, dando la norma

del vector error.
b) Teniendo en cuenta el rango de la matriz A, hallar su inversa generali-

zada.
c) Utilizar la inversa generalizada obtenida en el apartado anterior para

calcular la pseudosolución del sistema y hallar la norma del vector error.
Ejercicio 3.12 Consideremos el sistema de ecuaciones Ax = b, con

   
2 −2   6
  x1  
A =  1 −1  , x =   y b= 3  ,
  
  x2  
−2 2 3
y un vector unitario u. Se pide:
a) Demostrar que si H = I − 2uuT es la matriz de Householder, asociada al

vector u, entonces: H es ortogonal, H 2 = I y kHak2 = kak2 cualquiera
que sea el vector a.
b) Obtener la matriz de Householder que transforma el vector (2, 1, −2)T
en otro de la forma (α, 0, 0)T , con α > 0.
c) Aplicando el método de Householder, probar que el sistema Ax = b
posee infinitas soluciones en cuadrados mı́nimos y que el error cometido,
al considerar cualquiera de ellas, es el mismo.
d) Obtener la pseudosolución del sistema Ax = b. Es decir, la solución en
cuadrados mı́nimos, de entre las obtenidas en el apartado anterior, que
tenga menor norma euclı́dea.

   
0 3   −10
  x  
A =  −3 5 , x =   y b=
   
6 
  y  
4 0 −8
a) Probar que la matriz AT · A es definida positiva, obteniendo la factori-

zación de Choleski.
b) Plantear la iteración Xn+1 = L1 · Xn + c que se obtiene de aplicar el
método de Gauss-Seidel a las ecuaciones normales del sistema Ax = b.
¿Será convergente el proceso iterativo a la pseudosolución?
c) Hallar la matriz Hu = I − βuuT de la reflexión que transforma el vector
a = (0, −3, 4)T en el vector r = (−5, 0, 0).
d) Obtener la solución en mı́nimos cuadrados del sistema Ax = b, utilizando
el método de Householder, y determinar la norma del error.
e) Sin haber resuelto el apartado anterior, ¿podrı́an predecirse Hu A y Hu b
de las relaciones geométricas entre L =< u >, L⊥ y los vectores columnas
implicados?
Ejercicio 3.14 Se considera el sistema superdeterminado Ax = b con

   
3 2 3
   
A= 4 5  y b= 1 
   
   
12 0 13
a) Calcular la pseudosolución (solución de mı́nimos cuadrados) ası́ como la

norma del error utilizando transformaciones de Householder.
 
1 0 0
 
b) Sea T =  0 1 0  la matriz asociada a la transformación elemen-
 
 
0 0 1/12
tal que divide por 12 la tercera de las ecuaciones del sistema:
   
3 2   3
  x  
T Ax = T b ⇐⇒  4 5  = 1 
    
  y  
1 0 13/12
Calcular su pseudosolución haciendo uso de las ecuaciones normales. De-

terminar la norma del error.
c) ¿A qué se debe que no coincidan las pseudosoluciones obtenidas en los

dos apartados anteriores? ¿Qué habrı́a ocurrido si la matriz T hubiese
sido unitaria?

   
3 −2   2
  x  
A= 0 3 , x =   y b =  0 .
   
  y  
4 4 1
a) Probar que la matriz B = AT · A es definida positiva, obteniendo la

factorización de Cholesky B = GT · G.
b) Hacer uso de la factorización obtenida en el apartado anterior para hallar

la pseudosolución mediante las ecuaciones normales del sistema. Calcular
el número de condición, κ∞ (B), de la matriz B para la norma || ||∞ .
¿Hasta que punto se podrı́a considerar fiable la pseudosolución obtenida
con aritmética de ordenador?
c) Hallar la matriz de la reflexión (matriz de Householder) Hu que trans-

forma el vector a = (3, 0, 4)T en el vector r = (−5, 0, 0)T . Una vez de-
terminado el vector u, justificar que se pueden conocer Hu A y Hu b sin
necesidad de efectuar los productos.
d) Obtener la solución en mı́nimos cuadrados del sistema Ax = b, utilizando

el método de Householder y determinar la norma del error. Operando
con el ordenador, ¿puede obtenerse una pseudosolución distinta de la
obtenida en el apartado b? Si ası́ ocurriera, ¿puede ser mayor el error?

     
1 −1 2 x 0
     
A= 0 3 −3  , x =  y  y b =  1  .
     
     
0 −4 4 z 2
a) Hallar kAk∞ . ¿Qué se puede decir sobre el número de condición de la

matriz A para la norma infinito? ¿Qué estimación darı́a MATLAB para
el número de condición espectral obtenido con el comando cond(A)?
b) Utilizar la descomposición LU de la matriz AT A para resolver el sistema

AT Ax = AT b. ¿Qué propiedad caracteriza a las soluciones en relación al
sistema Ax = b? Interpreta geométricamente el resultado.
c) Encontrar una matriz ortogonal Q que transforme el vector a = (0, 3, −4)T

en el vector r = (0, 5, 0)T . Obtener la norma del error para las soluciones
en mı́nimos cuadrados del sistema QAx = Qb.
d) ¿Qué relación hay entre las soluciones obtenidas en los apartados ante-
riores?
Si se obtienen las soluciones en mı́nimos cuadrados del sistema Ax = b,
escalonando previamente la matriz A, ¿se debe obtener mismo resultado
que en alguno de los apartados anteriores?
 
2 3 4
−
 3 25 25

 
e) Probar que la matriz P =  1 3 4  es la pseudoinversa de A,
 3 25 − 25 
 
1
3
0 0
verificando las propiedades de Penrose. (Hacer la comprobación sólo con
dos de ellas).
De entre todas las soluciones en mı́nimos cuadrados del sistema Ax = b,
hallar la de menor norma euclı́dea.
Ejercicio 3.17
a) En lo que sigue, Hv denota la transformación de Householder asociada al

vector v. Sean x, y, v, z vectores no nulos, con Hv x = y y z ⊥ v. Probar
que Hv v = −v y Hv z = z. Determinar razonadamente todos los vectores
w tales que Hw x = y.
b) Se considera el sistema de ecuaciones dado por

    
1
− 1 0 x 2
 2    
1   y  =  −1 
    
 1 2
    
1 0 −1 z −1
b.1) Estudiar el condicionamiento del sistema, utilizando la norma 1.

b.2) Resolver el sistema por medio de transformaciones de Householder.
b.3) Desde un punto de vista numérico, ¿serı́a razonable resolver el sis-
tema escalonando por Gauss? Razonar la respuesta.
c) Demostrar que el vector c = (− 43 , 12 , − 4a

3
− 1)T y la matriz
 
2
0 −3 0
 
G= 0 1
0 −2 
 
 
0 − 2a 3
0
son los propios del método de Gauss-Seidel asociado al sistema

    
3
1 0 x −2
 2    
=
    
 0 2 1   y   1 
    
a 0 −1 z 1
d) Estudiar, en función del parámetro a, el carácter diagonal dominante

por filas de la matriz de coeficientes del sistema dado, ası́ como el radio
espectral de G. ¿Para qué valores de a es convergente el método anterior?
e) Para a = 0 el método resulta convergente. Utilizando aritmética exacta,

y toamndo como vector inicial x0 = (0, 0, 0)T , realizar dos iteraciones,
acotando el error cometido. Razonar qué ocurre cuando se itera por
tercera vez. ¿Hubiera ocurrido otro tanto al trabajar con aritmética de
ordenador?

   
1 1   2
  x  
A =  α 0 , x =  , b =  β , con α > 0 y β, γ ∈ R.
   
  y  
−2 2 γ
a) Hallar α sabiendo que que existe una matriz de Householder, Hv , que

transforma la primera columna de la matriz A en el vector r = (3, 0, 0)T .
¿Quién es Hv ?
b) Determinar el conjunto de vectores b para los que se verifica Hv b = b,

siendo Hv la matriz del apartado anterior. Encontrar, entre ellos, el que
tiene menor norma euclı́dea.
c) Hallar la pseudosolución del sistema Ax = bm , para α = 2 y bm =

(2, 1, −1)T , utilizando transformaciones ortogonales para determinar el
error.
d) Probar que si una matriz real B tiene sus columnas linealmente inde-
pendientes, entonces B T B es definida positiva.
e) Sea el sistema AT A x = AT bm , con α y bm como en el apartado (c).
e.1) ¿Serı́a posible utilizar una descomposición AT A = GGT , con G

triangular inferior, para resolver el sistema?
e.2) Utilizando la norma k k∞ para medir el condicionamiento, ¿es un
sistema mal condicionado para utilizar aritmética de ordenador en
su resolución?
e.3) Sea (s0 , s1 , s2 , . . .) la sucesión que se obtiene al aplicar el método
de Gauss-Seidel al sistema, con s0 = (0, 0)T . Probar que, operando
en aritmética exacta, la sucesión (sn ) es convergente y obtener su
lı́mite s.
4. Autovalores y autovectores
4.1 Conceptos básicos

Una variedad lineal V de Kn se dice que es invariante mediante una aplicación
A si cualquier vector x ∈ V se transforma, mediante A, en otro vector de V ,
es decir, si Ax ∈ V para cualquier vector x de V .
Las variedades invariantes más simples son aquellas que tienen dimensión 1,
por lo que una base está constituida por un único vector v 6= 0. Podemos
observar que en ese caso cualquier vector x ∈ V puede expresarse de la forma
x = αv con α 6= 0 si x 6= 0 y que, por tratarse de una variedad invariante, ha
de ser Ax = βv, pero entonces:
β
Ax = Aαv = αAv = βv =⇒ Av = v = λv
α
Las variedades invariantes de dimensión 1 vienen determinadas por un vector

v 6= 0 tal que Av = λv. Estos vectores reciben en nombre de autovectores o
vectores propios de la matriz A y los correspondientes valores de λ reciben el
nombre de autovalores o valores propios de A.
Es obvio que si A posee un autovector es porque existe un vector v 6= 0
tal que Av = λv, o lo que es lo mismo, tal que (λI − A)v = 0. Por tanto,
el sistema (λI − A)x = 0 es compatible y además indeterminado, ya que si
v 6= 0 es solución del sistema, cualquier vector proporcional a él también lo es.
Se verifica entonces que rg (λI − A) < n (donde n representa el orden de la
matriz) y, por tanto det(λI − A) = 0.
Al polinomio p(λ) = det(λI − A) se le denomina polinomio caracterı́stico de
la matriz A y a la ecuación p(λ) = 0 ecuación caracterı́stica .
Nótese que si λ es una raı́z de la ecuación caracterı́stica de la matriz A, existe
un autovector v asociado al autovalor λ. Por tanto, desde el punto de vista
103
104 Autovalores y autovectores
teórico, el problema del cálculo de los autovectores de una matriz se reduce

a la resolución de los sistemas (λi I − A)x = 0 obtenidos para las diferentes
raı́ces λi de la ecuación caracterı́stica.
Sea A una matriz cuadrada de orden n y supongamos que existen V1 , V2 , . . . , Vk

subespacios invariantes tales que Kn = V1 ⊕ V2 ⊕ · · · ⊕ Vk .
Si {xi1 , xi2 , . . . , xini } es una base de Vi se tiene que
B = {x11 , . . . , x1n1 , x21 , . . . , x2n2 , . . . , xk1 , . . . , xknk }
constituye una base de Kn . Si P es la matriz del cambio de base de la base

canónica a la base B, se tiene que
 
J Θ ··· Θ
 1 
J2 · · · Θ
 
−1
 Θ 
P AP =   .. .. . . .

. ..

 . . 
 
Θ Θ · · · Jk
donde cada Ji es una caja cuadrada de dimensión ri = dim(Vi ).

La justificación es que Axij = (0, . . . , 0, αi1 , . . . , αini , 0, . . . , 0)TB con 1 ≤ i ≤ k,
1 ≤ j ≤ ni y, por tanto, puede verse fácilmente que
 
J Θ ··· Θ
 1 
 Θ J2 · · · Θ 
 
AP = P   .. .. . . .. 

 . . . . 
 
Θ Θ · · · Jk
Proposicin 4.1 Autovectores asociados a autovalores diferentes son lineal-

mente independientes
Demostración. Sean u y v dos autovectores de la matriz A asociados a los

autovalores λ y µ respectivamente con λ 6= µ.
Si u y v fuesen linealmente dependientes se verificarı́a que u = αv con α 6= 0.
Entonces:
λu = Au = A(αv) = αAv = αµv = µ(αv) = µu
y por tanto (λ − µ)u = 0, pero dado que u 6= 0 se tiene que λ = µ en contra
de la hipótesis de que λ 6= µ, lo que prueba el resultado.
Método interpolatorio para la obtención del polinomio caracterı́stico 105
4.2 Método interpolatorio para la obtención

del polinomio caracterı́stico
El problema del cálculo de los autovectores de una matriz, una vez calculados
los autovalores, se reduce a la resolución de un sistema de ecuaciones por cada
autovalor.
Trataremos, por tanto y, en primer lugar, calcular sus autovalores a partir
de la propiedad de ser las raı́ces del polinomio caracterı́stico de la matriz.
Es decir, dada una matriz cuadrada A de orden n pretendemos calcular su
polinomio caracterı́stico P (λ) = det(λI − A) para, posteriormente, hallar sus
raı́ces mediante alguno de los métodos de resolución de ecuaciones estudiados.
El método interpolatorio consiste en calcular el polinomio caracterı́stico de la
matriz dando n valores a λ, para calcular n determinantes y, posteriormente,
resolver el sistema de n ecuaciones con n incógnitas resultante. Veámoslo con
detenimiento:
Sea P (λ) = det(λI − A) = λn + a1 λn−1 + · · · + an−2 λ2 + an−1 λ + an . Dando
a λ n valores λ1 , λ2 , . . . , λn se tiene que:
λ = λ1 =⇒ λn1 + a1 λn−1
1 + · · · + an−2 λ21 + an−1 λ1 + an = det(λ1 I − A)
λ = λ2 =⇒ λn2 + a1 λn−1
2 + · · · + an−2 λ22 + an−1 λ2 + an = det(λ2 I − A)
.........................................................................
.........................................................................
λ = λn =⇒ λnn + a1 λn−1
n + · · · + an−2 λ2n + an−1 λn + an = det(λn I − A)
 
1 3 2 5
 
2 −1
 
 4 0 
Ası́, por ejemplo, dada la matriz A =  
, su polinomio

 0 1 0 1 
 
2 −2 −1 1
4 3 2
caracterı́stico es de grado cuatro P (λ) = λ + a1 λ + a2 λ + a3 λ + a4 , por lo
que vamos a dar a λ cuatro valores:
λ=0 =⇒ a4 = det (−A) = 41
λ=1 =⇒ 1 + a1 + a2 + a3 + a4 = det(I − A) = 74
λ = −1 =⇒ 1 − a1 + a2 − a3 + a4 = det(−I − A) = −20
λ=2 =⇒ 16 + 8a1 + 4a2 + 2a3 + a4 = det(2I − A) = 67
dado que a4 = 41 el sistema se reduce a


a1 + a2 + a3 = 32  



−a1 + a2 − a3 = −62 =⇒ a1 = −4, a2 = −15 y a3 = 51




8a1 + 4a2 + 2a3 = 10

por lo que su polinomio caracterı́stico es
P (λ) = λ4 − 4λ3 − 15λ2 + 51λ + 41
4.3 Método de la potencia y variantes

Sea A ∈ Rn×n una matriz diagonalizable y supongamos que sus autovalores
verifican que:
|λ1 | > |λ2 | ≥ · · · ≥ |λn |
Sea B = {x1 , x2 , · · · , xn } una base de Rn formada por autovectores asociados
a λ1 , λ2 , . . . , λn respectivamente. Se verifica entonces que
A2 xi = A(Axi ) = A(λi xi ) = λi Axi = λi (λi xi ) = λ2i xi
por lo que es fácil probar, por inducción, que
Ak xi = λki xi para cualquier i = 1, 2, . . . , n
Dado un vector z0 ∈ Rn se define, a partir de él, la sucesión (zn ) con
zn = Azn−1 = A2 zn−2 = · · · = An z0 .
Si las coordenadas del vector z0 respecto de la base B son (α1 , α2 , . . . , αn ) se

tiene que z0 = α1 x1 + α2 x2 + · · · + αn xn , por lo que
Método de la potencia y variantes 107
zk = Ak z0 = Ak (α1 x1 + · · · + αn xn ) = α1 Ak x1 + · · · + αn Ak xn =
" n k #
k k k
X λi
= λ1 α1 x1 + · · · + λn αk xk = λ1 α1 x1 + αi xi
i=2
λ1
k
λi
Dado que |λ1 | > |λi | se tiene que lim = 0 ∀i = 2, 3, . . . , n.
k→∞ λ1
zk
Se verifica entonces que lim k = α1 x1 que es un autovector de A asociado
k→∞ λ1
al autovalor λ1 .
Si k es suficientemente grande, se tiene
Azk = zk+1 ' λk+1 k

1 α1 x1 = λ1 (λ1 α1 x1 ) = λ1 zk
por lo que la sucesión (zn ) nos proporciona un método para calcular el auto-
valor λ1 .
Nota: Para evitar que las coordenadas de los vectores zk sean demasiado
grandes se considera la sucesión formada por los vectores zn = Awn−1 donde
zn
wn = .
kzn k∞
Ejemplo 4.1 Para calcular, por el método

 de la potencia,
 el autovalor de
1 2 3
 
mayor valor absoluto de la matriz A =  2 2 3 , partiendo del vector
 
 
3 3 3
z0 = (10, 11, 1) obtenemos:
       
00 9091 30 1818 00 5303 40 8939
       
w0 =  10 0000 z1 =  40 0909 w1 =  00 6818 z2 =  50 4242
       
   
       
00 0909 60 0000 10 0000 60 6364
       
00 7374 50 3721 00 7009 50 2952
       
w2 =  00 8174 z3 =  60 1096 w3 =  00 7971 z4 =  50 9961
       
   
       
10 0000 70 6644 10 0000 70 4942
       
00 7066 50 3068 00 7057 50 3090
       
w4 =  00 8001 z4 =  60 0134 w5 =  00 7996 z6 =  60 0107
       
   
       
10 0000 70 5200 10 0000 70 5160
       
00 7058 50 3053 00 7058 50 3052
       
w6 =  00 7997 z7 =  60 0111 w7 =  00 7997 z8 =  60 0110
       
   
       
10 0000 70 5166 10 0000 70 5165
   
00 7058 50 3052
   
w8 =  00 7997 z9 =  60 0110
   
 
   
10 0000 70 5165
Tenemos, por tanto, que
 
00 7058
 
z9 = Aw8 = z8 = 70 5165  00 7997  = 70 5165 w8 = λw8
 
 
10 0000
es decir, el autovalor de mayor valor absoluto es 70 5165 y su autovector asociado

es w8 .
Este método sólo nos permite calcular, en caso de existir, el autovalor do-
minante (mayor valor absoluto) de una matriz. Existen sin embargo otras
variantes del método de la potencia que nos permiten calcular cualquiera de
sus autovalores a partir de una aproximación de ellos. Debido a la existencia de
dichas variantes el método estudiado anteriormente es conocido como método
de la potencia simple para distinguirlo de los que estudiaremos a continuación.
Conviene aquı́ recordar algunas propiedades de los autovalores de una matriz.
Teorema 4.1 Si A una matriz regular y λ un autovalor suyo asociado al

autovector v, se verifica que 1/λ es un autovalor de A−1 asociado al mismo
autovector v.
Demostración. Si λ es un autovalor de A asociado a v sabemos que Av = λv.

Al tratarse de una matriz invertible (det(A) 6= 0) sus autovalores son todos
no nulos, por lo que podemos dividir por λ y multiplicar por A−1 la última
1
igualdad para obtener que A−1 v = v es decir, 1/λ es un autovalor de A−1
λ
asociado a v.
Teorema 4.2 Si λ es un autovalor de una matriz A asociado a un autovector

v y α una constante cualquiera se verifica que λ − α es un autovalor de la
matriz A − αI asociado al mismo autovalor v.
Método de la potencia y variantes 109
Demostración. Sabemos, por hipótesis, que Av = λv, por lo que Av − αv =

λv − αv y, por tanto,
(A − αI)v = (λ − α)v
es decir, λ − α es un autovector de A − αI asociado a v.
Teorema 4.3 Si λ es un autovalor de una matriz A asociado a un autovector

v y α (una constante cualquiera) no es autovalor de A entonces 1/(λ − α) es
un autovalor de la matriz (A − αI)−1 asociado al autovector v.
La demostración se basa en los teoremas 4.1 y 4.2 y se deja al lector.
Sea A ∈ Rn×n una matriz diagonalizable regular para la que sus autovalores
verifican que:
0 < |λ1 | < |λ2 | ≤ · · · ≤ |λn | (4.1)
Los autovalores µ1 , µ2 , . . . , µn de la matriz A−1 son los inversos de los auto-

valores de A,
1
µi = para i = 1, 2, . . . , n.
λi
Invirtiendo la expresión (4.1) nos queda
1 1 1
> ≥ ··· ≥
|λ1 | |λ2 | |λn |
o lo que es lo mismo:
|µ1 | > |µ2 | ≥ · · · ≥ |µn |
es decir, el autovalor de menor valor absoluto de la matriz A se corresponde
con el de mayor valor absoluto de A−1 , por lo que aplicando el método de
1
la potencia a la matriz A−1 obtenemos el inverso (µ1 = ) del autovalor de
λ1
menor valor absoluto de la matriz A y, por tanto, dicho autovalor.
Obsérvese que debemos ir calculando, en cada paso, los vectores zn = A−1ωn−1
zn
y ωn = . Pues bien, al cálculo de zn se realiza resolviendo, por alguno
kzn k∞
de los métodos estudiados, el sistema Azn = ωn−1 lo que nos evita calcular
A−1 y arrastrar los errores que se cometan a lo largo de todo el proceso.
Este método es conocido como método de la potencia inversa y nos permite
calcular, en caso de existir, el autovalor de menor valor absoluto de una matriz
invertible A.
Veamos cómo una nueva variante de este último método nos permite calcular
un autovalor cualquiera de una matriz regular A a partir de una aproximación
suya. Sin embargo, hay que tener en cuenta que si el autovalor que vamos a
aproximar es múltiple o existen otros con el mismo módulo aparecen dificul-
tades que sólo podremos solventar utilizando otros métodos.
Consideremos una matriz A regular tal que todos sus autovalores tengan dife-
rente módulo y supongamos conocida una aproximación α de un determinado
autovalor λk . Si la aproximación es buena se verificará que |λk − α| < |λi − α|
para cualquier i = 1, 2, . . . n con i 6= k.
Dado que 1/(λi − α) son los autovalores de A − αI y ésta posee un autovalor
(λk − α) de menor valor absoluto que los demás, el método de la potencia
1
inversa nos proporcionará el valor µk = pudiéndose, a partir de éste
λk − α
1
último, hallar el valor de λk = + α.
µk
Está variante es conocida como método de la potencia inversa con traslación.
4.4 Cociente de Rayleigh

Hemos visto en el método de la potencia que el cálculo del autovalor y el
autovector asociado se realiza de forma paralela.
Si nos limitamos a calcular la sucesión (zn ) hasta obtener una aproximación x
adecuada del autovector podemos obtener el autovalor resolviendo el sistema
vectorial λ x = Ax. Este sistema resulta, en general, incompatible por no
ser x exactamente un autovector sino sólo una aproximación, por lo que la
solución que mejor se ajusta es la seudosolución del sistema, que nos dará una
aproximación del autovalor λ.
x∗ Ax
λ x∗ x = x∗ Ax =⇒ λ =
x∗ x
x∗ Ax
Al cociente ∗ se le denomina cociente de Rayleigh del vector x respecto
xx
de la matriz A.
Podemos, por tanto, obtener el autovector por el método de la potencia para
más tarde calcular el autovalor mediante el cociente de Rayleigh.
Para aplicar el método de la potencia inversa con traslación es necesario

disponer, previamente, de una aproximación del autovalor que se pretende
Teorema de Gerschgorin 111
calcular y para ello es necesario separar los autovalores de la matriz. El estudio

de los cı́rculos de Gerschgorin nos permite dar un primer paso en este sentido.
4.5 Teorema de Gerschgorin

Definición 4.1 Sea A ∈ Kn×n . Se definen los cı́rculos de Gerschgorin como
los conjuntos

C1 = {z : |z − a11 | ≤ r1 }  





C2 = {z : |z − a22 | ≤ r2 }  


 X n
. . . . . . . . . . . . . . . . . . . . . . . . .  en donde ri = |aij | para 1 ≤ i ≤ n


 j =1
......................... 


 j 6= i




C = {z : |z − a | ≤ r } 
n nn n
Teorema 4.4 Sean C1 , C2 , . . . , Cn los cı́rculos de Gerschgorin de una matriz

[n
A. Si λ es un autovalor de A, se verifica que λ ∈ Ci .
i=1
n
[
Demostración. Si λ 6∈ Ci es porque λ 6∈ Ci para ningún i = 1, 2, . . . , n,
i=1
por lo que
|λ − a11 | > r1
|λ − a22 | > r2
.............
.............
|λ − ann | > rn
y, por tanto, la matriz
 
λ − a11 −a12 ··· −a1n
 
 −a21 λ − a22 ··· −a2n
 

 = λI − A
.. .. ..

..
.
 
 . . . 
 
−an1 −an2 · · · λ − ann
es una matriz de Hadamard, por lo que det(λI − A) 6= 0 y, por tanto, λ no

puede ser un autovalor de la matriz A en contra de la hipótesis.
Teorema 4.5 Si un conjunto de k discos de Gerschgorin de una matriz A

constituyen un dominio conexo aislado de los otros discos, existen exactamente
k autovalores de la matriz A en dicho dominio.
4.6 Sensibilidad de los autovalores a las trans-

formaciones de semejanza
 
0 1 0 ··· 0 0
 
0 0 1 ··· 0 0
 
 
 
0 ···
 
 0 0 0 0 
Consideremos la matriz A =  .. .. .. . . .. ..
, que es una caja de
.
 
 . . . . . 
 
0 0 0 ··· 0 1
 
 
 
0 0 0 ··· 0 0
Jordan y, por tanto, no diagonalizable. Su polinomio caracterı́stico es λn . Si
introducimos en la matriz una perturbación
 y en vez de 
la matriz A tomamos
0 0 0 ··· 0 0
 
0 0 0 ··· 0 0 
 

 
0 0 0 ··· 0 0 
 

la matriz B = A + E con E =  .. .. .. . . .. ..  siendo ε > 0 muy

. . . 

 . . .
 
0 0 0 ··· 0 0 
 

 
ε 0 0 ··· 0 0
pequeño (incluso más pequeño que la resolución del ordenador) se obtiene
como polinomio caracterı́stico λn + (−1)n ε que posee n raı́ces distintas (las
raı́ces n-ésimas de (−1)n−1 ε), por lo que la matriz resulta ser diagonalizable.
Sea A una matriz diagonalizable y B una perturbación de dicha matriz (B =
A + E). Sean λi los autovalores de la matriz A y sea µ uno cualquiera de los
autovalores de la matriz B con λi 6= µ.
Teorema 4.6 En las condiciones anteriores se tiene que

mı́n |µ − λi | ≤ kP k P −1 kEk
i
Sensibilidad de los autovalores a las transformaciones de semejanza 113
Demostración. Por ser A diagonalizable existe P tal que P −1 AP = D. Sea

x un autovector de B asociado a µ. Entonces Bx = µx o lo que es lo mismo:
(A + E)x = µx ⇐⇒ (µI − A)x = Ex ⇐⇒ (µI − P DP −1 )x = Ex ⇐⇒
P (µI − D)P −1 x = Ex ⇐⇒ (µI − D)(P −1 x) = P −1 Ex = P −1 EP P −1 x
Supuesto que µ 6= λi cualquiera que sea i = 1, 2, . . . , n, la matriz µI − D es
regular, por lo que
P −1 x = (µI − D)−1 (P −1 EP )P −1 x, por lo que para cualquier norma multi-
plicativa se tiene
kP −1 xk ≤ k(µI − D)−1 k kP −1 EP k kP −1 xk ⇐⇒
1 ≤ (µI − D)−1 P −1 EP

−1 1 1
Dado que (µI − D) = diag ,..., tenemos que
µ − λ1 µ − λn

1
1 ≤ máx P −1 kEk kP k
i µ − λi
por lo que
mı́n |µ − λi | ≤ P −1 kP k kEk
i
Obsérvese que la perturbación cometida en los autovalores depende de la

matriz de paso P y dado que esta no es única, se trata de elegir, entre todas
las posibles matrices de paso, aquella que verifique que kP −1 k kP k sea mı́nima.
Es decir, aquella cuyo número de condición sea lo menor posible.
Corolario 4.7 Si A es unitariamente diagonalizable (diagonalizable mediante

una matriz de paso unitaria), la perturbación producida en los autovalores es
menor o igual que la perturbación producida en la matriz.
mı́n |µ − λi | ≤ kEk
i
Por tanto, las mejores matrices para el cálculo efectivo de sus autovalores y
autovectores son las diagonalizables unitariamente y éstas reciben el nombre
de matrices normales.
Es necesario aclarar que si una matriz A no es normal y le calculamos sus
autovalores, estos pueden reflejar, con la exactitud que se desee, sus verdaderos
valores, pero estos no van a representar la solución del problema que generó
dicha matriz, pues los elementos de A pueden arrastrar errores (redondeo,

medición, etc.) E y estos hacen los autovalores de las matrices A y A + E
puedan diferir bastante.
Por otra parte, es evidente que no podemos estudiar si una matriz es normal
calculando sus autovalores y autovectores para comprobar que se puede dia-
gonalizar por semejanza mediante una matriz de paso unitaria (matriz cons-
tituida por una base ortonormal de autovectores). Es necesario, por tanto,
encontrar otros métodos que detecten si una matriz es, o no es, normal.
Teorema 4.8 Sea T una matriz triangular. Si T ∗ T = T T ∗ entonces es dia-

gonal.
Demostración. Probaremos que se trata de una matriz diagonal por in-

ducción en el orden de la matriz.
   
a b a 0
Para n = 1 es obvio. Para n = 2 es T =   y T∗ =  .
0 c b c
    
∗
a 0 a b |a|2 ab
T T =  = 
b c 0 c ab |b|2 + |c|2
    
2 2
a b a 0 |a| + |b| bc
TT∗ =   = 
2
0 c b c bc |c|
Dado que se trata de una matriz normal es T ∗ T = T T ∗ por lo que igualando
ambas matrices se obtiene que |b|2 = 0 es decir b = 0 y, por tanto, T es
diagonal.
Supongamos ahora que el teorema es cierto para cualquier matriz triangular
y normal de orden n y vamos a probarlo para otra de orden n + 1.
   
a a · · · an+1 a 0 ··· 0
 1 2   1 
   
 0  ∗
 a2 
Sean T =  ..
yT =
 ..

Tn∗
 
 . Tn   . 
   
0 an
 n 
X  
2 2
|a i | |a1 |
TT∗ =  T ∗T = 
 

 i=1

∗
T T

Tn Tn∗ n
De la igualación de ambas obtenemos que |a2 |2 +|a3 |2 +· · ·+|an+1 |2 = 0, por lo

que a2 = a3 = · · · = an+1 = 0. Como, además, es Tn Tn∗ = Tn∗ Tn , por hipótesis
de inducción sabemos que Tn es diagonal y, por tanto, T es diagonal.
Teorema 4.9 Teorema de Schur

Cualquier matriz cuadrada A es unitariamente semejante a una triangular
superior T . Es decir, existe una unitaria U (U ∗ U = U U ∗ = I) tal que U ∗ AU =
T.
Demostración. Por inducción en el orden de A. Si n = 1 e obvio. Supuesto

cierto para n vamos a probarlo para una matriz de orden n + 1.
Sea A ∈ K(n+1)×(n+1) . Sea λ un autovalor de A y x un autovector unitario
(kxk = 1) asociado a λ.
Consideremos la matriz P = (x e2 · · · en ) en la que sus columnas x, e2 , . . . , en
constituyen una base ortonormal de Kn+1
 
λ αij
P ∗ AP =  
Θ An
 
1 Θ
Sea Q =   en donde Un es la matriz unitaria que, por hipótesis
Θ Un
de inducción, verifica que Un∗ An Un = Triangular superior.
Si consideremos la matriz U = P Q es fácil comprobar que U ∗ AU = Q∗ P ∗ AP Q
es una triangular superior.
Teorema 4.10 Una matriz A es normal si, y sólo si, AA∗ = A∗ A.
Demostración. Supongamos que AA∗ = A∗ A. Por el Teorema 4.9 sabemos

que existe una matriz unitaria U tal que U ∗ AU = T con T triangular superior.
Entonces A = U T U ∗ , por lo que
AA∗ = (U T U ∗ )(U T U ∗ )∗ = U T U ∗ U T ∗ U ∗ = U T T ∗ U ∗

=⇒
A∗ A = (U T U ∗ )∗ (U T U ∗ ) = U T ∗ U ∗ U T U ∗ = U T ∗ T U ∗
U T T ∗U ∗ = U T ∗T U ∗
es decir, T T ∗ = T ∗ T por lo que T es una matriz normal y triangular, por lo que
el Teorema 4.8 nos asegura que es diagonal y, por tanto, A es diagonalizable
unitariamente, es decir, es normal.
Recı́procamente, si A es unitariamente diagonalizable (normal) existe una

matriz unitaria U tal que U ∗ AU = D o lo que es lo mismo, A = U DU ∗ .
AA∗ = U DU ∗ (U DU ∗ )∗ = U DD∗ U ∗ = U diag(|d1 |2 , . . . , |dn |2 )U ∗
A∗ A = (U DU ∗ )∗ U DU ∗ = U D∗ DU ∗ = U diag(|d1 |2 , . . . , |dn |2 )U ∗
es decir, AA∗ = A∗ A.
Hemos visto que para que los autovalores de una matriz A reflejen la solución
del problema que la generó, ésta ha de ser normal. Si no podemos plantear
nuestro problema mediante una matriz normal lo que debemos hacer es tratar
de plantearlo mediante una matriz “lo más normal posible”. Vamos a estudiar
entonces el condicionamiento del problema.
Por el Teorema 4.9 cualquier matriz cuadrada A es unitariamente semejante
a una triangular superior T donde
     
t11 t12 · · · t1n t11 0 ··· 0 0 t12 · · · t1n
 0 t22 · · · t2n   0 t22 · · · 0   0 0 · · · t2n 
T =  .. ..  =  .. +
     
.. . . .. . . . .. .. .. .. 
 . . . .   . . . ..   . . . . 
0 0 · · · tnn 0 0 · · · tnn 0 0 ··· 0
Es decir, T = D + M y, obviamente, A es normal si, y sólo si, M = Θ.
En cierta manera (a kM k se le llama desviación de la normalidad de A)
la matriz M y más concretamente su norma, mide el condicionamiento del
problema.
Se tiene, además, que kM k2 es constante (no depende de U y de T sino sólo
de A), ya que
Xn
2 2 2 2 2 2
∗
kM k2 = kT k2 − kDk2 = kU AU k2 − kDk2 = kAk2 − |λi |2 =⇒
i=1
v
u n
X
2
|λi |2
u
kM k2 = kAk2 −
t
i=1
Definición 4.2 Se denomina matriz conmutatriz de una matriz cuadrada A

a la matriz C(A) = A∗ A − AA∗ .
Otra forma de medir el condicionamiento de una matriz es considerar la

norma euclı́dea de la matriz conmutatriz
kC(A)k2 = kAA∗ − A∗ Ak2
obteniéndose el siguiente resultado:
Teorema 4.11 Dada una matriz A ∈ Kn×n se verifica:
kC(A)k22 2
a) 2 ≤ kM k2 (Desigualdad de Eberlein).
6 kAk2
r
2 n3 − 3
b) kM k2 ≤ kC(A)k22 (Desigualdad de Heurici).
12
Los autovalores de una matriz son, en general, números complejos. Sin em-
bargo, las matrices hermı́ticas (en el caso real, las simétricas) tienen todos sus
autovalores reales como prueba el siguiente teorema.
Teorema 4.12 Los autovalores de una matriz hermı́tica son todos reales y
autovectores correspondientes a dos autovalores diferentes son ortogonales.
Demostración. Sea A una matriz hermı́tica, es decir, una matriz tal que
A∗ = A y sea λ un autovalor de A asociado al autovector x. Se verifica entonces
que Ax = λx.
Multiplicando la expresión anterior, por la izquierda por x∗ obtenemos que
x∗ Ax = x∗ λx = λx∗ x (4.2)
Trasponiendo y conjugando la expresión 4.2 obtenemos
(x∗ Ax)∗ = (λx∗ x)∗ =⇒ x∗ A∗ x = λ̄x∗ x =⇒ x∗ Ax = λ̄x∗ x (4.3)
Si comparamos las expresiones (4.2) y (4.3) obtenemos que
λx∗ x = λ̄x∗ x
y dado que x es un autovector (un vector no nulo) sabemos que x∗ x = kxk2 6= 0,
por lo que podemos dividir por x∗ x para obtener que λ = λ̄, es decir, λ ∈ R.
Por otra parte, si x e y son autovectores asociados a dos autovalores λ 6= µ
de una matriz hermı́tica A se verifica que
Ax = λx =⇒ (Ax)∗ = (λx)∗ =⇒ x∗ A = λx∗

x∗ Ay = λx∗ y =⇒ x∗ µy = λx∗ y =⇒ µx∗ y = λx∗ y =⇒ (λ − µ)x∗ y = 0
y dado que λ 6= µ =⇒ λ − µ 6= 0 obtenemos que
x∗ y = 0 ⇐⇒ x ⊥ y
Teorema 4.13 [Teorema espectral para matrices hermı́ticas] Sea

A una matriz hermı́tica de orden n. Existe una base de Cn constituida por au-
tovectores de A. Dicho de otra forma, toda matriz hermı́tica es diagonalizable
por semejanza (es normal).
4.7 Métodos iterados para la obtención de au-

tovalores y autovectores
4.7.1 Algoritmo QR de Francis

Dada una matriz A ∈ Kn×n , se pretende encontrar una sucesión de matrices
(An )n∈N convergente a la matriz triangular de Schur T (en cuya diagonal
aparecen los autovalores de la matriz A).
La construcción de dicha sucesión se hace de la siguiente manera: A0 = A.
Supuesta encontrada An , se realiza la factorización An = Qn Rn mediante
matrices de Householder y An+1 = Rn Qn . En otras palabras:
A0 = A = Q0 R0




 A1 = R0 Q0 = Q1 R1



 A2 = R1 Q1 = Q2 R2
 ...................



...................
Teniendo en cuenta que
A = Q0 R0 ⇒ A1 = R0 Q0 = Q∗0 AQ0
R0 Q0 = Q1 R1 ⇒ A2 = R1 Q1 = Q∗1 R0 Q0 Q1 = Q∗1 Q∗0 AQ0 Q1 =

= (Q0 Q1 )∗ A(Q0 Q1 )
.......................................................................
.......................................................................
Rn−1 Qn−1 = Qn Rn ⇒ An = Rn Qn = (Q0 Q1 · · · Qn−1 )∗ A(Q0 Q1 · · · Qn−1 )
.......................................................................
.......................................................................
En resumen, a partir de A se calcula Q0 y con ella A1 = Q∗0 AQ0 . Con

A1 calculamos Q1 y, a continuación, Q0 Q1 que nos permite calcular A2 =
(Q0 Q1 )∗ A(Q0 Q1 ). Con A2 calculamos Q2 y, a continuación Q0 Q1 Q2 (obsérvese
Métodos iterados para la obtención de autovalores y autovectores 119
que Q0 Q1 está almacenada en la memoria) que nos permite calcular ahora

A3 = (Q0 Q1 Q2 )∗ A(Q0 Q1 Q2 ) y ası́, sucesivamente.
 
1 2 3
Ejemplo 4.2 Para el cálculo de los autovalores de la matriz A =  2 2 3 
3 3 3
se tiene:
−10 9415 00 6671
   
1 2 3 7
A0 = A =  2 2 3  A1 =  −10 9415 −00 5385 00 1850 
3 3 3 00 6671 00 1850 −00 4615
 0
00 3365 00 0344
 0
7 5162 −00 0530 00 0016
 
7 5034
A2 =  00 3365 −10 1474 −00 1174  A3 = −00 0530 −10 1758 00 0348 
0 0 0 0
0 0344 −0 1174 −0 3554 0 0016 0 0348 −00 3404
0
 0
00 0083 00 0001
 0
7 5165 −00 0013 00 0000
 
7 5165
A4 =  00 0083 −10 1775 −00 0100  A5 = −00 0013 −10 1776 00 0029 
00 0001 −00 0100 −00 3390 00 0000 00 0029 −00 3389
 0
00 0002 70 5165
  
7 5165 0 0 0
A6 =  00 0002 −10 1776 −00 0008  A7 =  0 −10 1776 0 
0 0 0
0 −0 0008 −0 3389 0 0 −0 3389
Por lo que los autovalores de la matriz A son:
−10 1776 − 00 3389 y 70 5165
Nota: Si λ1 , λ2 , . . . , λn son los autovalores de una matriz A tales que |λ1 | >
|λ2 | > · · · > |λn | > 0, el algoritmo converge. En general se tiene que, si existen
autovalores de igual módulo, la sucesión converge a una matriz triangular por
cajas en la que cada caja contiene a los autovalores del mismo módulo.
En otras palabras, la sucesión converge, en general, a una matriz triangular
por cajas de tal manera que
a) Si todos los autovalores tienen distinto módulo el lı́mite de la sucesión

es una matriz triangular superior en la que los elementos de su diagonal
son los autovalores de la matriz.
b) Si existen autovalores de igual módulo la matriz lı́mite es una matriz

triangular superior por cajas en la que cada caja de orden k de su
diagonal es una matriz cuyos autovalores son todos los de la matriz A
de igual módulo.
Ası́, por ejemplo, si una matriz A de orden 3 tiene un autovalor real y dos auto-
valores complejos conjugados (por tanto de igual módulo) de distinto módulo
que el autovalor real se llegarı́a a una matriz del tipo
 
a11 a12 a13
 0 a22 a23 
 
0 a32 a33
!
a22 a23
donde a11 es el autovalor real y los autovalores de la matriz son
a32 a33
los dos autovalores complejos de la matriz A.
De forma más general, si A es una matriz cualquiera (normal o no), al apli-
carle el algoritmo, la sucesión converge (en la mayorı́a de los casos) a su forma
de Schur.
Definición 4.3 Se dice que una matriz A ∈ Kn×n es una matriz (superior)
de Hessenberg si aij = 0 para i > j + 1, es decir:
 
a11 a12 a13 · · · a1 n−1 a1n
 a21 a22 a23 · · · a2 n−1 a2n 
 
 0 a32 a33 · · · a3 n−1 a3n 
A =  ..
 
.. .. . . .. .. 
 . . . . . . 
 
 0 0 0 · · · an−1 n−1 an−1 n 
0 0 0 · · · an n−1 ann
Teorema 4.14 Dada una matriz A ∈ Kn×n existen n − 2 transformaciones

de Householder H1 , H2 , . . . , Hn−2 tales que la matriz
Hn−2 Hn−3 · · · H1 AH1 · · · Hn−3 Hn−2
es una matriz de Hessenberg.
Demostración. La demostración se hará por inducción en el orden n de la

matriz.
Si A es de orden 1 ó 2, el resultado es obvio. Supuesto cierto para matrices
de orden n vamos a probarlo para matrices de orden n + 1.
   
a11 a12 · · · a1 n+1 0
 a21 a22 · · · a2 n+1   v2 
Dada la matriz A =  .. sea v =  ..  el
   
.. .. .. 
 . . . .   . 
an+1 1 an+1 2 · · · an+1 n+1 vn+1
vector cuya transformación de Householder asociada H1 es tal que

   
a11 a11
 a21   ã21 
   
 a31   0 
H1  = .
 ..   .. 
 .   . 
an+1 1 0
 
1 0 0 ··· 0

 0 

Dicha transformación es de la forma H1 = 
 0 . Entonces:

..
H10
 
 . 
0
    
a11 ∗ ∗ ∗ ··· 1 0 0 ··· 0 a11 ∗ ∗ · · · ∗
 ã21
 ∗ ∗ ∗···
 0
   ã21
 


H1 AH1 =  0
 ∗ ∗ ···
∗  0

=
  0 

 .. .. .. ..  ..
.. 0   ..
H1   . An 

 . . . .
.  .
0 ∗ ∗ ··· ∗ 0 0
Es fácil comprobar que si para k = 2, 3, . . . , n − 1 hacemos
 
1 0 0 ··· 0
 0 
 
 0
Hk = 


 .. 0
Hk 

 .
0
0 0
siendo Hn−1 , Hn−2 , . . . , H20 las n−2 transformaciones de Householder tales que
0 0
Hn−1 Hn−2 0
· · · H20 An H20 · · · Hn−2 0
Hn−1
es una matriz de Hessenberg, la matriz
Hn−1 Hn−2 · · · H2 H1 An+1 H1 H2 · · · Hn−2 Hn−1
también es de Hessenberg.
4.7.2 Método de Jacobi para matrices simétricas reales

Comenzaremos viendo el método para dimensión dos y, más tarde, generaliza-
remos para una dimensión cualquiera.

a b
Dada una matriz simétrica y real A = el método trata de buscar
b c

cos α sen α
un valor para α de tal forma que la rotación P = haga
− sen α cos α
que la matriz P −1 AP sea diagonal.

−1 cos α − sen α a b cos α sen α ∗ 0
P AP = =
sen α cos α b c − sen α cos α 0 ∗
Se debe verifica, por tanto, que
(a − c) cos α sen α + b(cos2 α − sen2 α) = 0
es decir:
b(cos2 α − sen2 α) = (c − a) cos α sen α
π
Si c = a basta tomar α = .
4
Si c 6= a podemos dividir por c − a y obtenemos que
2b 2 cos α sen α sen 2α
= 2 2
= = tg 2α
c−a cos α − sen α cos 2α
2b
Llamando m = se tiene que tg 2α = m, por lo que
c−a
√
−1 ± 1 + m2
t = tg α =
m
y, a partir del valor de t obtenido, se pueden calcular
1 t
cos α = √ y sen α = √
1 + t2 1 + t2
valores, estos últimos, que nos permiten determinar la matriz P .
Algoritmo de cálculo

a b
a) Dada A =
b c
a.1) Si b = 0 FIN.
a.2) Si b 6= 0 y a = c
√ √ ! !
2/2 2/2 ∗ 0
P = √ √ y P −1 AP = . FIN.
− 2/2 2/2 0 ∗
√
2b −1 ± 1 + m2
b) m = t= .
c−a m
1 t
c) cos α = √ sen α = √ .
1 + t2 1 + t2
! !
cos α sen α ∗ 0
P = y P −1 AP = . F IN.
− sen α cos α 0 ∗
Para dimensiones superiores sea P la matriz diagonal por bloques

 
Ip−1
P = Tq−p+1 
In−q
en la que  
cos α sen α
Tq−p+1 =  Iq−p−1 
− sen α cos α
y todos los demás elementos nulos.
• Si apq = 0, hacemos cos α = 1 y sen α = 0.

√
2
• Si app = aqq , hacemos cos α = sen α = .
2
√
2apq −1 ± 1 + m2
• Si app 6= aqq llamando m = y t= hacemos
aqq − app m
1 t
cos α = √ y sen α = √
1 + t2 1 + t2
Entonces, los elementos que ocupan los lugares pq y qp de la matriz P −1 AP

son nulos.
El método de Jacobi consiste en buscar el elemento apq con p 6= q de mayor
módulo, anularlo mediante una matriz P1 , buscar el siguiente elemento de
mayor módulo y anularlo mediante otra matriz P2 y ası́ sucesivamente hasta
diagonalizar la matriz A.
Dado que las matrices Pk , o matrices de Jacobi, son ortogonales, se verifica
que Pk−1 = PkT por lo que se trata, en definitiva, de buscar P1 , P2 , . . . , Pk de la
forma descrita anteriormente, para obtener
PkT · · · P1T AP1 · · · PK = D

Nota: Lo más importante de este método es que el cálculo de las matrices

Pi puede hacerse simultáneamente ya que no se requiere el conocimiento de
ninguna de ellas en el cálculo de cualquier otra. Es decir, el método es para-
lelizable.
 
1 −2 3 1
 −2 0 4 2 
Ejemplo 4.3 Consideremos la matriz A =   3
 que es
4 1 1 
1 2 1 0
simétrica y real.
El elemento extradiagonal de mayor valor absoluto es a23 = 4, por lo que
comenzamos haciendo
√
2a23 −1 + 1 + m2 1 t
m= t= cos = √ sen = √
a33 − a22 m 1 + t2 1 + t2
P23 = I4 p22 = p33 = cos p23 = −p32 = sen
obteniéndose  
1 0 0 0
 0 0.7497 0.6618 0 
P23 =
 0 −0.6618

0.7497 0 
0 0 0 1
por lo que
 
1 −3.4848 0.9254 1
T
 −3.4848 −3.5311 0 0.8376 
P23 AP23 =
 0.9254

0 4.5311 2.0733 
1 0.8376 2.0733 0
Obsérvese que la transformación de semejanza sólo ha afectado a las filas y

columnas 2 y 3 dejando invariantes a los demás elementos. Si queremos aplicar
T
el método a la matriz P23 AP23 obtenida, utilizando como referencia el elemento
a14 = 1, podemos observar que todos los cálculos necesarios para obtener
la nueva matriz P14 podı́amos haberlos realizados al mismo tiempo que los
realizados para P23 ya que sólo necesitamos los valores de los elementos (1, 1),
T
(1, 4), (4, 1) y (4, 4) de la matriz P23 AP23 , que por haber quedado invariantes
son los mismos que los de la matriz A.
Realizando los cálculos necesarios de forma análoga al caso anterior obtene-
mos  
0.8507 0 0 −0.5257
 0 1 0 0 
P14 =  
 0 0 1 0 
0.5257 0 0 0.8507
T T T
que aplicarı́amos a P23 AP23 para obtener P14 P23 AP23 P14 = P T AP donde la
matriz P viene dada por P = P23 P14 .
Ahora bien, la matriz
 
0.8507 0 0 −0.5257
 0 0.7497 0.6618 0 
P = P23 P14 = 
 0 −0.6618 0.7497 0 
0.5272 0 0 0.8507
puede ser escrita directamente sin necesidad de construir previamente P23 y

P14 y multiplicarlas.
En resumen, se pueden enviar a un ordenador los valores de a22 , a23 = a32
y a33 para que calcule el ángulo de giro correspondiente y nos devuelva los
valores de p22 = p33 y p23 = −p32 , mientras que de forma simultánea enviamos
a otro ordenador los valores de a11 , a14 y a44 y nos devolverá los de p11 = p44
y p14 = −p41 .
En el ordenador central construimos la matriz P con los datos recibidos y
calculamos
 
1.6180 −2.5240 1.8772 0
 −2.5240 −3.5311 0 2.5445 
P T AP = 
 1.8772

0 4.5311 1.2771 
0 2.5445 1.2771 −0.6180
que volveremos a renombrar con A.

A la vista de la matriz enviarı́amos al Ordenador 1 los datos de los elementos
a22 , a24 = a42 y a44 mientras que enviarı́amos al Ordenador 2 los de a11 ,
a13 = a31 y a33 que nos devolverı́an los elementos necesarios para construir
una nueva matriz
 
0.8981 0 0.4399 0
 0 0.8661 0 0.5016 
P =  −0.4399

0 0.8981 0 
0 −0.5016 0 0.8651
y a partir de ella
 
0.6986 −1.6791 0 −1.6230
 −1.6791 −5.0065 −1.5358 0 
P T AP =  
 0 −1.5358 5.4506 0.4353 
−1.6230 0 0.4353 0.8573
Reiterando el proceso llegamos a la matriz

 
2.5892 0 0 0

 0 −5.6823 0 0 

 0 0 5.7118 0 
0 0 0 −0.6188
cuyos elementos diagonales son los autovalores de la matriz original.
Los autovalores de una matriz normal cualquiera pueden ser reales o com-
plejos mientras que los de una matriz hermı́tica son todos reales, por lo que
si pudiéramos transformar el cálculo de los autovalores de una matriz nor-
mal en los de otra hermı́tica habrı́amos simplificado el problema. Es más,
si pudiéramos transformarlo en el cálculo de los autovalores de una matriz
simétrica real, podrı́amos trabajar en el campo real en vez de hacerlo en el
campo complejo y, entre otras cosas, podrı́amos utilizar el método de Jacobi a
partir de una matriz normal cualquiera previa transformación en un problema
de autovalores de una matriz simétrica real.
En las siguientes secciones estudiaremos cómo pueden realizarse dichas trans-
formaciones.
4.8 Reducción del problema a matrices hermı́-

ticas
Proposicin 4.2 Dada cualquier matriz cuadrada A, existen dos matrices hermı́ticas
H1 y H2 tales que A = H1 + iH2 .
Demostración. Si se quiere que A = H1 + iH2 con H1 y H2 hermı́ticas,

se tendrá que A∗ = H1∗ − iH2∗ = H1 − iH2 , por lo que resolviendo el sistema
resultante se tiene que
A + A∗ A − A∗
H1 = H2 =
2 2i
además, dado que el sistema es compatible determinado, la solución es única.
Teorema 4.15 Sea A = H1 +iH2 una matriz normal con H1 y H2 hermı́ticas.

Si x 6= 0 es un autovector de A asociado al autovalor α + iβ (α, β ∈ R) se
verifica que
H1 x = αx y H2 x = βx
Reducción del problema a matrices hermı́ticas 127
Demostración. Por ser A normal, existe una matriz unitaria U tal que
U ∗ AU = D
 
α + iβ 0 ··· 0
 0 α2 + iβ2 ··· 0 
D = U ∗ AU =   =⇒
 
.. .. .. ..
 . . . . 
0 0 · · · αn + iβn
   
α 0 ··· 0 β 0 ··· 0
 0 α2 ··· 0   ∗
 0 β2 ··· 0 
 ∗
A=U U + iU  U =⇒
 
.. .. .. ..   .. .. .. ..
 . . . .   . . . . 
0 0 ··· αn 0 0 · · · βn
   
α 0 ··· 0 β 0 ··· 0
 0 α2 ··· 0 
 ∗
 0 β2 ··· 0 
 ∗
H1 = U  U H2 = U  U
 
.. .. .. .. .. .. .. ..
 . . . .   . . . . 
0 0 · · · αn 0 0 · · · βn
donde H1 y H2 son hermı́ticas. Se tiene entonces que α es un autovalor de
H1 asociado a x (primera columna de U ) y β un autovalor de H2 asociado
también a x.
Las componentes hermı́ticas H1 y H2 de una matriz cuadrada A nos propor-
cionan otro método para estudiar si la matriz A es, o no es, normal.
Teorema 4.16 Una matriz cuadrada A es normal si, y sólo si, sus compo-
nentes hermı́ticas conmutan.
A es normal ⇐⇒ H1 H2 = H2 H1
Demostración.
A + A∗ A − A∗ A2 − AA∗ + A∗ A + (A∗ )2
H1 H2 = · =
2 2i 4i
A − A∗ A + A∗ A2 + AA∗ − A∗ A + (A∗ )2
H 2 H1 = · =
2i 2 4i
a) Si A es normal se verifica que AA∗ = A∗ A por lo que
A + A∗ A − A∗ A2 + (A∗ )2 

H1 H2 = · = 
2 2i 4i


=⇒ H1 H2 = H2 H1
A − A∗ A + A∗ A2 + (A∗ )2 

H2 H1 = · =


2i 2 4i
b) Si H1 H2 = H2 H1 se verifica que
A2 − AA∗ + A∗ A + (A∗ )2 = A2 + AA∗ − A∗ A + (A∗ )2
por lo que
−AA∗ + A∗ A = A2 + AA∗ − A∗ A ⇐⇒ 2A∗ A = 2AA∗ ⇐⇒ A∗ A = AA∗
y, por tanto, A es normal.
4.9 Reducción del problema a matrices simé-

tricas reales
Sea H una matriz hermı́tica y sean hij = aij + ibij con 1 ≤ i, j ≤ n sus
elementos. Podemos descomponer la matriz H de la forma H = A + iB donde
A y B son matrices reales con A = (aij ) y B = (bij ). Por ser H hermı́tica se
verifica que

 A = AT
∗ ∗ ∗ T T
A+iB = H = H = A −iB = A −iB (por ser A y B reales) =⇒
 B = −B T
Sea x = a + ib un autovector de H asociado al autovalor α (recuérdese que

α ∈ R por ser H una matriz hermı́tica). Se verifica entonces que

Hx = (A + iB)(a + ib) = (Aa − Bb) + i(Ab + Ba) 
αx = α(a + ib) = αa + iαb 


 Aa − Bb = αa
Hx = αx ⇒
 Ab + Ba = αb
por lo que
A −B a a
=α
B A b b

A −B
Obsérvese además que la matriz real es simétrica, ya que
B A
T
A −B AT B T A −B
= =
B A −B T AT B A
Aplicación al cálculo de las raı́ces de un polinomio 129

A −B
Los autovalores de la matriz son, por tanto, los mismos que
B A
A −B
los de H y si un autovector de es (x1 , . . . , xn , xn+1 , . . . , x2n ), el
B A
correspondiente autovector de la matriz H es (x1 + ixn+1 ), . . . , (xn + ix2n ) .
Dada una matriz normal calcularemos sus autovalores de la siguiente forma:
Algoritmo de cálculo de los autovalores de una matriz normal
Paso 1 Calcular las componentes hermı́ticas

A + A∗ A − A∗
H1 = y H2 =
2 2i
de la matriz A.
Paso 2 Calcular los autovalores λi (1 ≤ i ≤ n) de la primera componente

hermı́tica H1 y sus autovectores asociados vi reduciendo previamente el
problema al caso de matrices reales y simétricas.
vi∗ H2 vi
Paso 3 Calcular los cocientes de Rayleigh µi = ∗ .
vi vi
Como los vectores vi (1 ≤ i ≤ n) son también autovectores de H2 , los µi
obtenidos son sus autovalores asociados.
Paso 4 Los autovectores de A son los mismos vectores vi y sus autovalores

asociados vienen dados por λi + iµi .
4.10 Aplicación al cálculo de las raı́ces de un

polinomio
Consideremos el polinomio
P (x) = a0 xn + a1 xn−1 + a2 xn−2 + · · · + an−1 x + an
Dado que el polinomio caracterı́stico de la matriz
 
−a1 /a0 −a2 /a0 · · · −an−1 /a0 −an /a0
 
 0 
A=
 
.. 

 In−1 . 

0
(donde In−1 representa la matriz unidad de orden n − 1) es precisamente P (x),

la raı́ces de dicho polinomio coinciden con los autovalores de la matriz A, por
lo que dichas raı́ces pueden ser obtenidas, en bloque, mediante un método
iterado de cálculo de autovalores.
Ası́, por ejemplo, MATLAB en su comando roots(P) para calcular las raı́ces
de un polinomio P construye, en primer lugar, la matriz A definida anterior-
mente y se limita luego a calcular sus autovalores aplicando internamente el
comando eig(A) de cálculo de los autovalores de una matriz. Este comando
lo que hace, en primer lugar, es aplicar un algoritmo para llevar la matriz A a
una matriz de Hessenberg y posteriormente aplicarle el algoritmo QR mediante
transformaciones de Householder.

Ejercicio 4.1 Realizar la descomposición de Schur de la matriz
 
−6 9 3
A =  −2 3 1  .
−4 6 2
00 6 00 8

Ejercicio 4.2 Comprobar que la matriz U = es unitaria (or-
−00 8 00 6
togonal) y obtener, basándose en ella, una matriz normal A que tenga por
autovalores 2 y 3i. Calcular la conmutatriz de A y comprobar que sus compo-
nentes hermı́ticas conmutan.
Ejercicio 4.3 Probar, basándose en el teorema de Gerchgorin, que la matriz:

 
9 1 −2 1
 0 8 1 1 
A=  −1

0 7 0 
1 0 0 1
tiene, al menos, dos autovalores reales.

2 + 3i 1 + 2i
Ejercicio 4.4 Dada la matriz A = se pide:
1 + 2i 2 + 3i
a) Comprobar que es normal sin calcular la matriz conmutatriz.

b) Calcular sus autovalores a partir de los de sus componentes hermı́ticas.
c) Comprobar que estos autovalores están en el dominio de Gerchgorin.
 
a 1 1
Ejercicio 4.5 Dada la matriz A =  1 a 1  donde a es un número com-
1 1 a
plejo cualquiera, se pide:
a) Obtener su polinomio caracterı́stico.
b) Probar que tiene por autovalores: λ = a − 1 doble y λ = a + 2 simple.
c) Calcular los autovectores y comprobar que no dependen de a.

 
2−i 0 −2 + 4i
Ejercicio 4.6 Dada la matriz A =  0 4 − 5i 2 − 4i , se pide:
−2 + 4i 2 − 4i 3 − 3i
a) Probar que es normal.
b) Obtener su primera componente hermı́tica H1 y calcular el polinomio

caracterı́stico de dicha componente.
c) Calcular los autovalores de H1 y, por el mismo método anterior, sus

autovectores.
d) Teniendo en cuenta que estos autovectores también lo son de la matriz

H2 (segunda componente hermı́tica), calcular sus autovalores.
e) Obtener, a partir de los resultados anteriores, los autovalores y autovec-

tores de A, ası́ como la matriz de paso unitaria U tal que U ∗ AU = D.

2 1
1 0
a) Calcular su polinomio caracterı́stico por el método interpolatorio.
b) Tomar una aproximación, con dos cifras decimales exactas, del mayor de
los autovalores y afinarla con el cociente de Rayleigh.
 
1 2 3
Ejercicio 4.8 Dada la matriz A =  2 2 3 
3 3 3
a) Hallar sus autovalores mediante el algoritmo QR.

b) Hallar el autovalor de mayor valor absoluto, por el método de la potencia,

partiendo del vector (10, 11, 1)

1 + i −2 + i
2−i 1+i

−i
a) Comprobar que es normal y que v1 = es un autovector de su
1
primera componente hermı́tica H1 asociado al autovalor λ1 = 0.
b) Calcular el otro autovalor (y un autovector asociado) de la matriz H1
aplicando el método de la potencia.

x y
c) Considérese la matriz real y simétrica S = . Probar que la
y x
t cos α sen α
transformación de Jacobi Q SQ con Q = y α = π/4
− sen α cos α
nos anula los elementos extradiagonales.
d) Transformar el problema del cálculo de los autovalores de la matriz H2
(segunda componente hermı́tica de la matriz A) al del cálculo de los au-
tovalores de una matriz C simétrica real y comprobar que son suficientes
dos transformaciones de Jacobi Q1 y Q2 , del tipo de las del apartado
anterior, para diagonalizar dicha matriz C y obtener los autovalores de
H2 .
e) Obtener, a partir de las columnas de la matriz Q = Q1 Q2 , los autovec-
tores de la matriz H2 . ¿Cuáles son los autovalores de la matriz A?
Ejercicio
4.10 Sea λ = α + iβ, α, β ∈ R, autovalor de la matriz A =

2i −2
.
2 2i
a) Utilizar el teorema de Gerschgorin para probar que el único autovalor

real de la matriz A sólo puede ser λ = 0.
b) Probar que A? = −A y deducir, a partir de ello, que A es una matriz nor-

mal. ¿Puede ser NO diagonalizable una matriz compleja que verifique
esa relación?
c) Utilizar la descomposición hermı́tica de la matriz, A = H1 + iH2 , para

deducir que la parte real de los autovalores de A tiene que ser α = 0.
d) Hallar el autovalor dominante (de máximo módulo) de la componente

hermı́tica H2 aplicando el método de la potencia. ¿Quién es el autovalor
dominante de A?
Sugerencia: Iniciar el método con el vector v1 = (1, 0)T .
e) Si se perturba la matriz A en la matriz
(2 − 10−3 ) i

−2
A + δA = ,
2 (2 + 10−2 ) i
hallar la norma euclı́dea de la matriz δA. ¿Puedes encontrar una cota

del error E = |µ − λ|, transmitido al autovalor dominante?
Indicación: |µ − λ| ≤ ||P || ||P −1 || ||δA||, siendo P −1 AP diagonal.
Ejercicio 4.11
a) Probar que las raı́ces del polinomio bλ + cλ2 + λ3 son los

 p(λ) = a +
0 1 0
autovalores de la matriz A(p) =  0 0 1 .
−a −b −c
b) Si el método de la potencia aplicado a la matriz A(p) converge a un

vector v, ¿qué relación tiene v con las raı́ces del polinomio p(λ)?
c) Si el algoritmo QR aplicado a la matriz A(p) converge a una matriz

triangular T , ¿qué relación tiene T con las raı́ces del polinomio p(λ)?
d) Si se puede obtener la factorización LU de la matriz P A(p), siendo P

una matriz de permutación, ¿quiénes tienen que ser P, L y U ?
Ejercicio 4.12 Sean el polinomio p(x) = x3 +2x2 −2x−4 y su correspondiente

matriz A = A(p), definido en el ejercicio 4.11
a) Utilizar una sucesión de Sturm para probar que el polinomio p(x) tiene
sus raı́ces reales y que sólo una de ellas, que denotaremos α, es positiva.
b) Utilizar el método de Newton para obtener una aproximación de la raı́z

α, garantizando 5 cifras decimales exactas.
c) Obtener la pseudosolución, β, del sistema (A2 v)x = A3 v, determinando

la norma del error, para v = (1, 1, 1)T . ¿Deberı́a ser β una aproximación
de α?
d) Obtener la matriz de Householder que transforma el vector a = (0, 0, 4)T

en el vector b = (4, 0, 0)T . ¿Se podı́a haber predicho el resultado?
e) Obtener la factorización QR de la matriz A, utilizando el método de

Householder. (Sugerencia: ¡el apartado anterior!)
f) Dar el primer paso del algoritmo QR aplicado a la matriz A. Indicar

cómo podrı́a el método de Gram-Schmidt utilizarse para los sucesivos
pasos del algoritmo y si esto serı́a una buena decisión para obtener las
raı́ces del polinomio p(x).
Ejercicio 4.13
a) ¿Qué pasos se dan para calcular los autovalores de una matriz cuadrada
A mediante el algoritmo QR? y ¿que forma tiene la matriz a la que
converge el algoritmo en los siguientes casos?
a.1) Si todos sus autovalores tienen distinto módulo.

a.2) Si existen autovalores de igual módulo.
 
−4 2 −4 3
 1 0 0 0
b) El polinomio caracterı́stico de la matriz A = 
 0
 es
1 0 0
0 0 1 0
precisamente el polinomio P (x) del Ejercicio 1. Calculando sus autova-
lores mediante el algoritmo QR el proceso converge a la matriz
 
−4.64575131106459 4.07664693269566 1.32820441231845 −2.21143157264058
 
−0.24888977635522 −0.86635866374600
 
 0 0.58988079050108 
 
 
 0 1.22575806673700 0.24888977635522 0.03848978825890 
 
0 0 0 0.64575131106459
Calcular, a partir de dicha matriz, las raı́ces del polinomio (autovalores

de A).
c) Al aplicar
 el método
 de la potencia y comenzando el proceso con el vec-

1 1
 1   −0.2152 
tor x = 
 1  se obtiene en la cuarta iteración el vector  0.0463 .
  
1 −0.0100
Determinar una aproximación de la raı́z de mayor valor absoluto del po-
linomio P (x) (autovalor correspondiente) utilizando el cociente de Ray-
leigh.
Ejercicio 4.14 Sean las matrices A, An y B definidas como:

 0
1 671 00 242 20 164
    
0 1 0 0 1 0
A= 0 0 1 , An =  00 00 −00 50 10 47  y B =  0 0 1 .
0 0 0 0
3 1 0 0 00 −0 81 −1 16 3 1 01
a) Aplicando el algoritmo QR real a la matriz A se obtiene (iterando sufi-

cientemente), como aproximación “aceptable” del método, la matriz An .
¿Por qué las matrices A y An deben tener los mismos autovalores?
Hallar las aproximaciones de los autovalores de la matriz A que se ob-
tienen de An .
b) Tomando v0 aleatoriamente, ¿se debe esperar convergencia o divergencia

en el método de la potencia aplicado a la matriz A?
Empezar en v0 = (1, 1, 1)T y determinar los tres primeros vectores v1 ,
v2 y v3 que proporciona el método. Hallar la mejor aproximación del
autovalor dominante de A, en norma || ||2 , que se obtiene con v3 .
c) Estudiar si A es una matriz normal. Si se perturba A en la matriz

B = A + δA, hallar la medida de la perturbación ||δA||2 .
¿Se podrı́a asegurar que los autovalores dominantes de las matrices A y
B difieren, a lo más, en 0’1?
 
1 1 2
Ejercicio 4.15 Sean A =  0 1 1  con 0 < ε ≤ 1,
   1 −1 ε
0 −1 −2 0
J = 0 0 −1  y b =  −1 .
−1 1 0 2
a) Obtener la factorización A = LU . Utilizar la factorización obtenida para

resolver el sistema Ax = b.
b) Hallar el número de condición κ∞ (A) de la matriz A para la norma || ||∞ .

Razonar si el resultado del apartado anterior, obtenido con aritmética
de ordenador, podrı́a ser considerado fiable para ε próximo a cero.
c) Para ε = 1, comprobar que J es la matriz de la iteración xn+1 = J ·

xn + c que se obtiene al aplicar el método de Jacobi al sistema Ax = b.
Determinar c y, empezando en x1 = (1, 0, 0)T , hallar el vector x3 .
d) Hallar la aproximación λ3 del autovalor dominante λ de la matriz J

utilizando el método de la potencia, con v0 = (1, 0, 0)T , y el cociente de
Rayleigh para determinar λ3 con el valor obtenido para v3 .
Sabiendo que λ3 tiene una cota de error estimada en e < 00 5. ¿Es
suficiente dicha aproximación para analizar la convergencia de la sucesión
(xn ) del método de Jacobi?
e) Para ε = 0, hallar la solución en mı́nimos cuadrados del sistema A0 x = b
que se obtiene al suprimir la primera columna de A, utilizando las ecua-
ciones normales. Determinar el error y justificar el resultado obtenido.
f) Analizar si es posible encontrar la matriz H de Householder que trans-
forma la segunda columna de A0 en el vector b. En caso afirmativo, ¿es
normal la matriz H resultante?
 
3 0 −1
Ejercicio 4.16 Considérese la matriz A =  1 −2 2 .
−1 −1 8
a) Hacer uso de los cı́rculos de Gerschgorin para estudiar el número de
autovalores reales que posee.
Obtener su polinomio caracterı́stico P (λ) y un intervalo de amplitud 1
que contenga a su autovalor dominante.
b) Comprobar que la fórmula de Newton-Raphson asociada a dicho polino-
mio es
2λ3 − 9λ2n − 39
λn+1 = ϕ(λn ) = n2 .
3λn − 18λn + 3
Sabiendo que la gráfica de la función y = ϕ0 (λ)
en el intervalo [8, 9] viene dada por la figura
adjunta, ¿podemos garantizar la convergencia
del método de Newton partiendo de cualquier
punto λ0 ∈ [8, 9]?
Nota: ϕ(λ) es la función que aparece en la

fórmula de Newton-Raphson.
c) Si tomamos λ0 = 8 ¿con qué error se obtiene la aproximación λ1 ?

d) ¿Existe algún vector v0 para el que podamos garantizar la convergencia
del método de la potencia simple aplicado a la matriz A? ¿Qué aproxi-
mación se obtiene para el autovalor dominante aplicando el cociente de
Rayleigh al vector v1 si partimos de v0 = (1 0 0)T ?
e) Aplicando el método QR para el cálculo de los autovalores de A, con

una aritmética de ordenador con una precisión de cuatro decimales, el
método se estabiliza en la matriz An en la que hemos omitido dos de sus
elementos x e y
 0
8 0195 −00 5134 20 7121

An =  0 20 7493 10 5431 
0 x y
¿Puede ser nulo el elemento x?, ¿se pueden determinar los elementos
que faltan sin necesidad de volver a aplicar el algoritmo QR? ¿Sabrı́as
decir cuál es la aproximación obtenida para los otros dos autovalores de
la matriz A?
 
0 1 0
Ejercicio 4.17 Se considera la matriz A =  0 0 1 .
−0.25 −0.125 1
a) Demostrar que las raı́ces del polinomio P (x) = 2+x−8x2 +8x3 coinciden
con los autovalores de A. Acotar y separar las raı́ces de P (x), indicando
cuántas raı́ces reales y complejas tiene. Comparar los resultados con la
información que se desprende del estudio de los cı́rculos de Gerschgorin.
b) Determinar un intervalo de amplitud 0.5 con un extremo entero que

contenga a la raı́z negativa de P (x). Razonar si se verifican en dicho in-
tervalo las condiciones de Fourier. Aproximar por el método de Newton-
Raphson dicha raı́z con 2 cifras decimales exactas.
c) Tomando como vector inicial z0 = (0, 1, 0)T , realizar dos iteraciones del
método de la potencia inversa. Por medio del cociente de Rayleigh aso-
ciado al vector hallado, determinar una aproximación del autovalor de A
correspondiente. ¿Qué relación existe entre éste valor y la aproximación
hallada en el apartado anterior? ¿Puede haber autovalores de la matriz
A en el cı́rculo de centro 0 y radio 14 ? Razonar las respuestas.
d) Al aplicar el algoritmo QR a la matriz A se obtiene como salida la matriz

T = Q∗ AQ, para cierta matriz unitaria Q. ¿Puede ser T una matriz
triangular superior? Justificar la respuesta.
Ejercicio 4.18
a) Utilizar el método interpolatorio para determinar el polinomio carac-

terı́stico p(λ) de la matriz
 
2 −1 0
A =  0 −2 1 
1 0 5
b) A la vista de los cı́rculos de Gerschgorin, ¿se puede garantizar que el

algoritmo QR aplicado a la matriz A convergerá a una matriz triangular
con sus autovalores en la diagonal? ¿Se puede garantizar la convergencia
del método de la potencia simple si comenzamos a iterar con el vector
v = (1 1 1)T ?
c) Haciendo uso de los cı́rculos de Gerschgorin, determinar cuántos auto-

valores reales posee y calcular un intervalo de amplitud 1 y extremos
enteros que contenga al autovalor dominante.
d) Comprobar que, en dicho intervalo, se verifican TODAS las hipótesis de

Fourier para la convergencia del método de Newton. ¿En qué extremo
deberı́amos comenzar a iterar?
e) Tomando x0 = 5 y aplicando el método de Newton, ¿con cuántas cifras

exactas se obtiene x1 ?
Índice
Algoritmo Desigualdad
de Horner, 26, 27 de Eberlein, 117
de la bisección, 5 de Heurici, 117
de Newton, 14 Desviación de la normalidad, 116
QR de Francis, 118 Distancia, 40
Autovalor, 103 inducida, 40
Autovector, 103
Eberlein
Bisección desigualdad de, 117
algoritmo de la, 5 Ecuaciones
método de la, 4 algebraicas, 20
Bolzano Ecuación
teorema de, 4 caracterı́stica, 103
Error
Ceros
absoluto, 2
de polinomios, 20
relativo, 2
de una función, 1
Espacio normado, 39
Cholesky
Estabilidad, 3
factorización de, 56
Cociente de Rayleigh, 110 Factorización
Condicionamiento, 3 de Cholesky, 56
Condición LU, 53
número de, 3, 46, 48 ortogonal, 73
Convergencia, 41 Fórmula
de Heron, 15
Descenso más rápido
de Newton-Raphson, 12
método de, 65
Fourier
Descomposición
regla de, 15, 17
en valores singulares, 90
Función
Descomposición
ceros de una, 1
de Schur, 115
contractiva, 8
método de, 63
139
140 Índice
Gauss-Seidel Método
método de, 64 consistente, 59
Gerschgorin convergente, 59
cı́rculos de, 111 de descomposición, 63
teorema de, 111 de Gauss-Seidel, 64
Gradiente conjugado de Jacobi, 63
método de, 65 para matrices simétricas rea-
les, 121
Hadamard de la bisección, 4
martiz de, 55 de la potencia, 106
Heron de la regula falsi, 6
fórmula de, 15 de Newton, 12
Hessenberg para raı́ces múltiples, 18
matriz de, 120 de relajación, 64
Heurici de Sturm, 23
desigualdad de, 117 del descenso más rápido, 65
Horner del gradiente conjugado, 65
algoritmo de, 26, 27 directo, 1, 46, 52
Householder interpolatorio, 105
transformación de, 76 iterado, 1, 46, 59
transformación en el campo com-
plejo, 78 Newton
algoritmo de, 14
Jacobi método de, 12
método de, 63, 121 para raı́ces múltiples, 18
Laguerre Newton-Raphson
regla de, 22 fórmula de, 12
Norma, 39
Matriz euclı́dea, 40
conmutatriz, 116 infinito, 40
de diagonal dominante, 55 matricial, 41
de Hadamard, 55 euclı́dea, 42
de Hessenberg, 120 infinito, 42
fundamental, 54 uno, 41
sparce, 46 multiplicativa, 39
triangular uno, 40
inferior, 46 vectorial, 39
superior, 46 Número de condición, 3, 46, 48
tridiagonal, 46
unitaria, 43 Penrose
Índice 141
seudoinversa de, 91 compatible

Pivote, 54 determinado, 47
Polinomio mal condicionado, 47
caracterı́stico, 103 superdeterminado, 85
Potencia Sturm
método de la, 106 método de, 23
Punto fijo sucesión de, 23
teorema del, 8 Sucesión
de Sturm, 23
Radio espectral, 44
Rayleigh Teorema
cociente de, 110 de Bolzano, 4
Raı́ces de Gerschgorin, 111
acotación de, 21, 22 de Rolle, 4
de una ecuación, 1 de Rouche-Fröbenius, 85
múltiples, 1 del punto fijo, 8
separación de, 23 Fundamental del Álgebra, 20
simples, 1 Transformación
Regla de Householder, 76
de Fourier, 15, 17 en el campo complejo, 78
de Laguerre, 22 unitaria, 43
de Ruffini, 27
Relajación Valor propio, 103
método de, 64 Variedad
Rolle invariante, 103
teorema de, 4 Vector
Rouche-Fröbenius propio, 103
teorema de, 85
Ruffini
regla de, 27
Régula falsi
método de la, 6
Schur
descomposición de, 115
Seudoinversa
de Penrose, 90, 91
Seudosolución, 87
Sistema
bien condicionado, 46

Ap Algnum

Cargado por

Copyright:

Formatos disponibles

Ap Algnum

Cargado por

Información del documento

Título original

Derechos de autor

Formatos disponibles

Compartir este documento

Compartir o incrustar documentos

Opciones para compartir

¿Le pareció útil este documento?

¿Este contenido es inapropiado?

Copyright:

Formatos disponibles

Ap Algnum

Cargado por

Copyright:

Formatos disponibles

E.T.S.

Fco. Javier Cobos Gavala

2 Sistemas de ecuaciones lineales 39

2.1.4 Normas matriciales . . . . . . . . . . . . . . . . . . . . 41

3 Sistemas inconsistentes y sistemas indeterminados 73

4 Autovalores y autovectores 103

Dada una función no nula f : C → C, resolver la ecuación f (x) = 0 es hallar

Definición 1.1 Una solución x de la ecuación f (x) = 0 se dice que tiene

Todos los métodos numéricos de resolución de ecuaciones presentan dificul-

1.1 Errores y condicionamiento en problemas

Datos =⇒ Algoritmo =⇒ Resultados

|f (x) − f (x)| ' |x − x| · |f 0 (x)|

Para el problema inverso, es decir conocido f (x) buscar el valor de x (resolver

la tangente de 44◦ 590 resultando éste un proceso bien condicionado, ya que

1.2 Método y algoritmo de la bisección: análisis

Teorema 1.1 [Teorema de Bolzano] Si f es una función continua en el

Nuestro problema se reduce a localizarla. Para ello, supongamos que está

Teorema 1.2 [Teorema de Rolle] Si f (x) es una función continua en el

En efecto, si f (x) tuviese dos raı́ces α1 y α2 en el intervalo [a, b], verificarı́a

c) Si ε < 10−k (error prefijado), entonces FIN. El valor de α0 es la raı́z

El proceso debe repetirse hasta que

f (mi ) = 0 o bien bi − ai < ε con ε > 0 prefijado.

Se tiene, por tanto:

while (b − a)/2 > ε

El hecho de calcular el punto medio de [a, b] como m = a + b-a/2 es debido a

Ejemplo 1.1 Supongamos que se quiere calcular la raı́z cuadrada de 3, para lo

Una variante del método de la bisección es el método de la regula falsi o de

Figura 1.1: Método de la régula f alsi.

La pendiente de dicha secante viene determinada por

Despejando el valor de c obtenemos que

pudiéndose dar los mismos casos que en el método de la bisección, es decir:

• Si f (c) = 0 la raı́z buscada es c.

• Si f (a) y f (c) tienen signos contrarios, la raı́z se encuentra en el intervalo

El algoritmo quedarı́a de la siguiente forma:

Aplicando el algoritmo al Ejemplo 1.1 obtenemos la raı́z de 3, con 14 cifras

En el método de la regula falsi la convergencia suele ser más rápida que en

1.3 Punto fijo e iteración funcional

|f (x1 ) − f (x2 )| < |x1 − x2 | ∀ x 1 , x2 ∈ R

Si la función es derivable, basta comprobar que |f 0 (x)| ≤ q < 1 cualquiera

Teorema 1.3 [Teorema del punto fijo] Dada la ecuación x = ϕ(x) en

Demostración. Dado que xn+1 = ϕ(xn ) y xn = ϕ(xn−1 ) se tiene que

xn+1 − xn = ϕ(xn ) − ϕ(xn−1 ) = (xn − xn−1 )ϕ0 (c)

|xn+1 − xn | ≤ q |xn − xn−1 |

x0 + (x1 − x0 ) + (x2 − x1 ) + · · · + (xn+1 − xn ) + · · ·

observamos que la suma parcial Sn+1 = xn y que la serie es absolutamente con-

x = lim xn = lim ϕ(xn−1 ) = ϕ(lim xn−1 ) = ϕ(x)

Además, es la única solución de la ecuación en el intervalo [a, b], ya que de

x − x0 = ϕ(x) − ϕ(x0 ) = (x − x0 )ϕ0 (c) con c ∈ (x, x0 )

1.3.1 Cota del error “a posteriori”

intervalo cerrado [xn , x] ó [x, xn ] (dependiendo de que x sea mayor o menor

Lo único que debemos exigir es que la derivada de la función no se anule en

Ejemplo 1.2 El cálculo de la raı́z cuadrada de 3 equivale al cálculo de la raı́z

Es decir, hemos escrito la ecuación de la forma x = ϕ(x) con

Obsérvese que en el Ejemplo 1.1 vimos cómo eran necesarias 46 iteraciones

1.4 Análisis del método de Newton-Raphson

Si construimos, utilizando la fórmula de Newton-Raphson, la sucesión {xn } y

siempre que f 0 (lim xn ) 6= 0, lo cual se verifica si exigimos que la función posea

En la Figura1.3 puede observarse cómo actúa geométricamente el método de

Figura 1.3: Interpretación geométrica del método de Newton.