Apuntes Calculo Numerico

E.T.S.
DE INGENIERÍA INFORMÁTICA
INGENIERÍAS TÉCNICAS EN INFORMÁTICA

DE
SISTEMAS Y GESTIÓN
Apuntes de
CÁLCULO NUMÉRICO
por
Fco. Javier Cobos Gavala
DEPARTAMENTO DE
MATEMÁTICA APLICADA I
Contenido
1 Ecuaciones no lineales 1
1.1 Introducción . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.2 Método y algoritmo de la bisección: análisis de errores . . . . 4
1.2.1 Algoritmo . . . . . . . . . . . . . . . . . . . . . . . . . 6
1.3 Punto fijo e iteración funcional . . . . . . . . . . . . . . . . . . 7
1.3.1 Cota del error “a posteriori” . . . . . . . . . . . . . . . 9
1.4 Método de Newton: análisis de errores . . . . . . . . . . . . . 11
1.4.1 Algoritmo . . . . . . . . . . . . . . . . . . . . . . . . . 13
1.4.2 Regla de Fourier . . . . . . . . . . . . . . . . . . . . . 14
1.4.3 Método de Newton para raı́ces múltiples . . . . . . . . 17
1.5 Cálculo de ceros de polinomios . . . . . . . . . . . . . . . . . . 19
1.5.1 Sucesiones de Sturm . . . . . . . . . . . . . . . . . . . 19
1.5.2 Algoritmo de Horner . . . . . . . . . . . . . . . . . . . 22
1.6 Sistemas de ecuaciones no lineales . . . . . . . . . . . . . . . . 24
1.6.1 Método de Newton . . . . . . . . . . . . . . . . . . . . 26
1.7 Ejercicios propuestos . . . . . . . . . . . . . . . . . . . . . . . 29
2 Sistemas de ecuaciones lineales 35

2.1 Normas vectoriales y matriciales . . . . . . . . . . . . . . . . . 35
2.1.1 Normas vectoriales . . . . . . . . . . . . . . . . . . . . 36
2.1.2 Distancia inducida por una norma . . . . . . . . . . . . 36
2.1.3 Convergencia en espacios normados . . . . . . . . . . . 37
i
ii Contenido
2.1.4 Normas matriciales . . . . . . . . . . . . . . . . . . . . 37

2.1.5 Transformaciones unitarias . . . . . . . . . . . . . . . . 39
2.1.6 Radio espectral . . . . . . . . . . . . . . . . . . . . . . 40
2.2 Sistemas de ecuaciones lineales . . . . . . . . . . . . . . . . . . 41
2.3 Número de condición . . . . . . . . . . . . . . . . . . . . . . . 43
2.4 Factorización LU . . . . . . . . . . . . . . . . . . . . . . . . . 49
2.5 Factorización de Cholesky . . . . . . . . . . . . . . . . . . . . 52
2.6 Métodos iterados . . . . . . . . . . . . . . . . . . . . . . . . . 55
2.6.1 Método de Jacobi . . . . . . . . . . . . . . . . . . . . . 59
2.6.2 Método de Gauss-Seidel . . . . . . . . . . . . . . . . . 60
2.6.3 Métodos de relajación (SOR) . . . . . . . . . . . . . . 60
2.7 Métodos del descenso más rápido y del gradiente conjugado . 61
2.7.1 Método del descenso más rápido . . . . . . . . . . . . . 63
2.7.2 Método del gradiente conjugado . . . . . . . . . . . . . 63
2.8 Factorizaciones ortogonales . . . . . . . . . . . . . . . . . . . . 64
2.9 Transformaciones de Householder . . . . . . . . . . . . . . . . 65
2.9.1 Interpretación geométrica en Rn . . . . . . . . . . . . . 66
2.9.2 Householder en Cn . . . . . . . . . . . . . . . . . . . . 68
2.10 Factorización QR de Householder . . . . . . . . . . . . . . . . 70
2.11 Sistemas superdeterminados. Problema de los mı́nimos cuadrados 74
2.11.1 Transformaciones en sistemas superdeterminados . . . 77
2.12 Descomposición en valores singulares y seudoinversa de Penrose 79
2.12.1 Seudoinversa de Penrose . . . . . . . . . . . . . . . . . 81
2.13 Ejercicios propuestos . . . . . . . . . . . . . . . . . . . . . . . 83
3 Interpolación 93
3.1 Introducción . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93
3.2 Interpolación polinomial . . . . . . . . . . . . . . . . . . . . . 94
3.2.1 Interpolación de Lagrange . . . . . . . . . . . . . . . . 94
3.2.2 Interpolación de Newton . . . . . . . . . . . . . . . . . 98
Contenido iii
• Diferencias divididas . . . . . . . . . . . . . . . . . . 98
• Diferencias finitas . . . . . . . . . . . . . . . . . . . 102
3.2.3 Interpolación de Hermite . . . . . . . . . . . . . . . . . 106
3.3 Interpolación por splines . . . . . . . . . . . . . . . . . . . . . 109
3.3.1 Splines cúbicos . . . . . . . . . . . . . . . . . . . . . . 109
3.3.2 Cálculo de los splines cúbicos de interpolación . . . . . 111
3.4 Ejercicios . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113
4 Integración numérica 115

4.1 Introducción . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115
4.2 Fórmulas de cuadratura . . . . . . . . . . . . . . . . . . . . . 116
4.3 Fórmulas de Newton-Cötes . . . . . . . . . . . . . . . . . . . . 119
4.3.1 Fórmula del trapecio . . . . . . . . . . . . . . . . . . . 120
4.3.2 Fórmula de Simpson . . . . . . . . . . . . . . . . . . . 121
4.4 Fórmulas compuestas . . . . . . . . . . . . . . . . . . . . . . . 122
4.4.1 Simpson para n par . . . . . . . . . . . . . . . . . . . . 122
4.4.2 Trapecios para n impar . . . . . . . . . . . . . . . . . . 123
4.5 Ejercicios . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123
Índice 127
1. Ecuaciones no lineales
1.1 Introducción
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 .
1
2 Ecuaciones no lineales
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.
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
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.
Introducción 3
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.
Proposición 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).
|x|n − 1
|P (x)| ≥ |a0 | |x|n − 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.
Proposición 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, aunque en la práctica sólo se utiliza en el caso de las
ecuaciones algebraicas, por lo que lo veremos en la Sección 1.5.
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
Método y algoritmo de la bisección: análisis de errores 5
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
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

b−a
intervalo obtenido en la en la iteración n-ésima, viene dado por εn = n+1 ,
2
1 −3
por lo que si b − a = 1 y n = 9 se tiene que ε9 < 10 < 10 , es decir, en 9
2
iteraciones 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

 f (mi ) = 0

El proceso debe repetirse hasta que o bien

 b − a < ε con ε > 0 prefijado.
i i
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 sign(f (a)) = sign(f (m))
a←m
end if
if sign(f (b)) = sign(f (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.
Punto fijo e iteración funcional 7
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
n+1
< 10−14 =⇒ 2n+1 > 1014 =⇒ n ≥ 46
2
es decir, tendrı́amos que detenernos en m46 para poder garantizar la precisión
exigida.
Vamos a ver a continuación otros métodos que reducen, de forma considera-

ble, el número de operaciones.
1.3 Punto fijo e iteración funcional

Ya se comentó que los métodos iterados consisten en crear una sucesión
convergente 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] 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 |
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.1 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).
Punto fijo e iteración funcional 9
(a) (b)

x0 x0
(c) (d) ?
x0 x0
Figura 1.1: Esquema de la convergencia para el teorema del punto fijo.
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 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)
f (x) − f (xn )
ó c ∈ (x, xn ) tal que = f 0 (c).
x − xn
Como f (x) = 0 y (x − xn ) = εn , nos queda que εn = − ff(x n)
0 (c) , obteniéndose que

|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
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 < 0 donde f (x) = x2 − 3, por lo que
mı́n |f (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.
Método de Newton: análisis de errores 11
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, mientras que ahora sólo hemos necesitado 26. Sin
embargo vamos a ver a continuación cómo se puede reducir aún más el numero
de iteraciones aplicando el método conocido como método de Newton.
1.4 Método de Newton: análisis de errores

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 − 0 ⇒ f (lim xn ) = 0
f (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

f (xn )
punto de abscisa x = xn − 0 que es pre-
f (xn )
cisamente el valor de xn+1 de la fórmula de

Newton-Raphson.
En la figura adjunta puede observarse cómo

actúa geométricamente el método de Newton- x x2 x1 x0
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.

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 = ε n+1 + ε
f 0 (xn ) 2f 0 (xn ) 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
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.
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, ya que sólo hemos necesi-
tado 5 iteraciones frente a las 46 que se necesitan en el método de la bisección.
De hecho, existen métodos para determinar el valor inicial x0 que debe to-
marse para que en la segunda iteración se disponga ya de 8 cifras decimales
exactas.
1.4.2 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.2), 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.2: 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 escribién-
dola de la forma sen x = x, por lo que las soluciones serán los puntos de
intersección de la curva y = sen x con y = x.
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
de Newton. y = sen x
-1
-
0 1
xn − sen xn
xn+1 = xn − =
1 − cos xn
sen xn − xn 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.
Cálculo de ceros de polinomios 19
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 Cálculo de ceros de polinomios

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.
Si la función es un polinomio P (x) hemos visto que el polinomio
P (x)
Q(x) =
mcd (P (x), P 0 (x))
posee los mismos ceros que P (x) pero todos simples. Con lo que el primer paso
a la hora de calcular los ceros de un polinomio es eliminar sus raı́ces múltiples.
Vamos a estudiar, por tanto, un método que nos permite separar las raı́ces
de una ecuación algebraica.
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 25
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 el teorema del valor medio para varias variables al decir

que
F (x) − F (xn ) = F 0 (α)(x − xn )
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 ∂Fm 
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 27
y haciendo uso del desarrollo de Taylor obtenemos que
0 = f (x) = f (xn + h) ≈ f (xn ) + hf 0 (xn )
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 29
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 Se considera la ecuación real 2 cos(2x) + 4x − k = 0.
a) Determinar el valor de k para que tenga una única raı́z triple en el

intervalo [0, 1].
b) Para k = 3, probar que posee una única raı́z simple en el intervalo [0, 1], y
calcularla con 6 cifras decimales exactas utilizando el método de Newton.
Ejercicio 1.3 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.4 Resolver, por los métodos de la bisección y Newton, la ecuación

ln x − sen x = 0, acotando previamente sus raı́ces.
Ejercicio 1.5 Separar las raı́ces reales de la ecuación xe−x − x2 + 1 = 0, y ob-

tenerlas con ocho cifras decimales exactas por el método de Newton, aplicando
previamente la Regla de Fourier.
Ejercicio 1.6 Dada la ecuación ex − (x − 2)2 = 0, probar que sólo posee

una raı́z real y obtenerla, por el método de Newton, con seis cifras decimales
exactas.
Ejercicio 1.7 Dada la ecuación ex − (x + 1)2 = 0, se pide:
a) Estudiar gráficamente sus raı́ces reales y acotarlas.
b) Obtener la mayor de ellas con dos cifras decimales exactas por el método
de la bisección.
c) Obtenerla con seis cifras decimales exactas por el método de Newton.

Ejercicio 1.8 La ecuación 00 81(x−1)−ln x = 0, tiene dos raı́ces reales, una de

las cuales es la unidad. Calcular la otra por el método de Newton, estudiando
previamente el campo de convergencia.
Ejercicio 1.9 Se considera la ecuación (x−1) ln x2 −2x2 +7x−7 = 0. Separar

sus raı́ces y obtener la mayor de ellas con seis cifras decimales exactas por el
método de Newton aplicando, previamente, la regla de Fourier.
2 x2 − 7x + 7
Ejercicio 1.10 Dada la ecuación e−x − = 0 se pide:
10 · (x − 1)2
a) Determinar el número de raı́ces reales que posee y separarlas.
b) Demostrar que para cualquier x > 10 6 es f 0 (x) < 0 y f 00 (x) > 0.
c) Calcular la mayor de las raı́ces, con dos cifras decimales exactas, por el
método de Newton.
Ejercicio 1.11 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.12 Dada la ecuación 8x3 − 4x2 − 18x + 9 = 0, acotar y separar

sus raı́ces reales.
Ejercicio 1.13 Dada la ecuación x3 − 6x2 + 3x + 9 = 0, acotar y separar sus

raı́ces reales.
Ejercicio 1.14 Dada la ecuación x3 − 3ax − 2b = 0 y basándose en el método

de Sturm, discutir para qué valores de a y b, existe una única raı́z real.
Ejercicio 1.15 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.16 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.17 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).
1
b) Si expresamos la ecuación P (x) = 0 de la forma x = F (x) = (x3 −6x2 +
3
7) y, partiendo de un x0 ∈ (6, 7), realizamos el proceso xn+1 = F (xn ),
¿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.18 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.19 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.20 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?
n 00
Indicación: En+1 = 1 (kE )2 , con k = max |f (x)| en un intervalo ade-
k 1 2 min |f 0 (x)|
cuado.
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].
35
36 Sistemas de ecuaciones lineales
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 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:
Xn
• 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.
Normas vectoriales y matriciales 37
• 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).
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
n
X
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
Normas vectoriales y matriciales 39
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
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 .
Proposición 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 .
T p p
x U = (xT U )T = U T x = (U T x)∗ (U T x) = x∗ (U T )∗ U T x =
2 2 2
p p √ √
= x (U ) U x = x∗ (U U ∗ )T x = x∗ I T x = x∗ x = kxk2.
T ∗ T T
Lema 2.2 Las matrices A, A∗ , AU y U A, donde U es una matriz unitaria,

poseen los mismos valores singulares.
Proposición 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
Sistemas de ecuaciones lineales 41
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.
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
Número de condición 43
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.3 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:

 kx − xk ≈ ε sistema bien condicionado

A − A < ε
=⇒
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.
Ejemplo 2.1 Si consideramos el sistema
3x + 4y = 7
x

1

de solución =
3x + 4.00001y = 7.00001 y 1
r1 r10 r1 r10

r

%

%
@ %
@ %

@rr

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

Sistema bien condicionado Sistema mal condicionado
Figura 2.1: Condicionamiento de un sistema.
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 =⇒
−1
A ≤ kIk
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 4.00001 −4
A−1 =
0.00003 −3 3
Utilizando la norma infinito kAk = n · máx |aij | se tiene que
i,j

kAk = 2 · 4.00001 
 64
4.00001  =⇒ κ(A) ' · 105 > 2 · 106
−1
kA k = 2 ·  3
0.00003
Se trata pues, de un sistema mal condicionado.

Xn
Si utilizamos la norma fila kAk = máx |aij | obtenemos:
i
j=1

kAk = 7.00001 
 56
8.00001 =⇒ κ(A) ' · 105 > 1.9 · 106
kA−1 k = 
 3
0.00003
obteniéndose, también, que se trata de un sistema mal condicionado.
Propiedades del número de condición κ(A).
a) Como ya se ha visto anteriormente κ(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:

kA−1 k

−1 1 −1 −1
κ(B)=kBk B = kzAk A
= |z| kAk = kAk A = κ(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
A = 1
−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 49
2.4 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 51
(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.5 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 53
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 u22 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 55
    
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.6 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 subordinada 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 subordinada.
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.6.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.6.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.6.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 =⇒

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

Lω = (D − ωE)−1 (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 61
• 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.7 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 63
tk v (k) - (k)
v
-

"3
"

"
x(k)
"
"
"
" (k+1)

"" x

""
"

2.7.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.7.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
2.8 Factorizaciones ortogonales

 
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
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 , y1 i h a2 , y1 i · · · h an , y 1 i
 h a1 , y2 i h a2 , y2 i · · · h an , y 2 i 
Q∗ A = 
 
.. .. .. .. 
 . . . . 
h a1 , yn i h a2 , yn i · · · h an , y n i
Transformaciones de Householder 65
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.
Al resolver el sistema Ax = b, si realizamos la descomposición QR obtenemos
QRx = b o lo que es lo mismo, Rx = Q∗ b que se trata de un sistema triangular.
Podemos observar que con la descomposición LU o la de Cholesky, lo que

conseguimos es transformar la resolución de nuestro sistema en la de dos sis-
temas triangulares. Ası́, por ejemplo, si descomponemos A como un producto
LU , lo que nos queda es LU x = b. Haciendo U x = y podemos resolver el sis-
tema triangular Ly = b para obtener su solución y. Nos queda ahora U x = y
que se trata de otro sistema triangular y, por tanto, de fácil resolución.
En la descomposición QR el proceso es otro. Dado que Q es unitaria se tiene,
inmediatamente la matriz Q∗ , por lo que se puede conocer de forma directa
el valor de Q∗ b, con lo que el sistema lo transformamos en Rx = Q∗ b que
es un sistema triangular. Es decir, la ventaja de esta descomposición es que
sólo debemos resolver unsistema triangular, frente a los dos que es necesario
resolver en los casos LU ó Cholesky.
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. Por tanto, se nos plantean las mismas
dificultades que en el referido método de ortonormalización.
Ello nos lleva a tratar de buscar un método por el que podamos realizar una
descomposición QR sin necesidad de aplicar Gram-Schmidt. Este proceso lo
lograremos mediante las denominadas Transformaciones de Householder.
2.9 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
Proposición 2.3 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.
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.
2.9.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
Q
6
Q y x
36

y + λv Q
Q

x − λv
Q
Q
−v Q
Q
v
Q
-
-
λv −λv
Figura 2.2: 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.
2.9.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 (2.4)
|v1 |2 + kxk2 − |x1 |2 = 1 (2.5)
De la ecuación 2.4 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 2.4 y 2.5 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
 
 
Hv x =  ..  con y1 = kxk eiα
.
 
 
0
que resulta ser  x1 
∓ kxk
 |x1 | 
0
 
Hv x = y = x − v = 
 
.. 
.
 
 
0
2.10 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
x1 − 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)
r a · · · a1n
 11 12 
(1) (1) 
 0 a22 · · · a2n  (1) (1)

(1)
H1 A =  .

.. . . .  = a1 a2 · · · an

 .. . . .. 
 
(1) (1)
0 an2 · · · ann
Factorización QR de Householder 71
(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 =  e y =  =  0 , se tiene que
     
2  1  0
     
−2 0 0
   
−2 − √13
x1 − y1 1    
v1 = = √  2  =  √13
   
kx1 − y1 k

2 3   
−2 − √13
 
− √13
2 ∗
 
H1 = I − ∗ v 1 v1 = I − 2  √1  − √13 √1 − √13 =
 
v1 v1  3  3
1
− √3
     
1
1 0 0 3
− 13 1
3
1
3
2
3
− 23
     
= 0 0  − 2  − 13 1
− 1 
 =  23 1 2
    
1 3 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
Factorización QR de Householder 73
    
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
 
.. . .
 . . . . 
 
a11
 a21 
x1 = a1 =  .. 
 
 . 
an1
 
kxk
 1 ± |a | a11 
 11 
definimos el vector v1 = 
 a 21 , obteniéndose que

 .. 
 . 
an1
 
(1) (1)
α a · · · a1n
 11 12 
 0 a(1) · · ·
(1) 

21 a 2n 
H1 A = 
 .. .. . . .. 

 . . . . 
 
(1) (1)
0 an2 · · · ann
   
(1)
a12 y1
   
 (1) 
 a22   y2 
 
   
Buscamos ahora otro vector v2 tal que H2  a(1) =  y de tal
   
32   0
 ..   .. 
   
 .   . 
   
(1)
a 0
 n2
α
 11 
 
 0 
forma que mantenga invariante al vector 
 .. . Bastará coger, para ello,

 . 
 
0
T
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
   
 0 v2 v2 · · · v2 vn
   
 v2  
H2 = I − 2 
 .. 
 0 v ··· v
2 n = 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
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.
2.11 Sistemas superdeterminados. Problema

de los mı́nimos cuadrados
Dado un sistema de ecuaciones de la forma Ax = b en el A es una ma-
triz cuadrada de orden n, A ∈ Kn×n , y x y b son vectores de Kn sabemos,
Sistemas superdeterminados. Problema de los mı́nimos cuadrados 75
por el teorema de Rouche-Fröbenius , que tiene solución si, y sólo si, existen
x1 , x2 , . . . , xn ∈ Kn tales que
          
a · · · a1n x1 b a11 a b1
 11   1  1n
 .. ..  ..   ..  .   .   .. 
      
..
.  . = .  ⇔ x1  .. + · · · +xn  .. = . 

 . .
          
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 a 1   b
 1  m  1 
⇐⇒  =  b2 
  
ma2 + n = b2  a2 1  
  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
Sistemas superdeterminados. Problema de los mı́nimos cuadrados 77
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
     
a∗n αn a∗n
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.
2.11.1 Transformaciones en sistemas superdetermina-

dos
Sabemos que dado un sistema compatible Ax = b y mediante transforma-
ciones 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 transformaciones unitarias.
Dado que las transformaciones de Householder son unitarias, podemos resol-
ver un sistema superdeterminado mediante transformaciones de Householder.
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
Θ
 
0
b
 1 
 .. 
 .   
  T
o llamando Hb = b0 =  b0n , T ∗ Θ   x = T ∗ Θ b0 .
 
 .. 
  Θ
 . 
 
0
bm
 
0
b
 1 
 .. 
 
 .  0
b
 1 
. 
 
Es fácil comprobar que T ∗ Θ  b0n  = T ∗  .. , por lo que el
  
 .. 
   
 .  0
bn
 
b0m
cálculo de la seudosolución del sistema superdeterminado Ax = b se hace
Descomposición en valores singulares y seudoinversa de Penrose 79
resolviendo el sistema  
b01
 . 
 
T ∗ T x = T ∗  .. 
 
0
bn
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̃
 
b0n
Una vez calculada

la seudosolución, la norma del error está representada por
la distancia b − b̃ que viene dada por

   
0 0
 b b

0  1   1 

b

 ..   ..  

 1 

 ..   .   .  
   .      b 0

T   0   0   n+1 
  bn   bn   . 

  0 
x −  bn  =  −  =  .. 
 0   b0n+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 peudosoluciones del sistema.

2.12 Descomposición en valores singulares y

seudoinversa de Penrose
La descomposición en valores singulares es otra factorización matricial que
tiene muchas aplicaciones.

Teorema 2.19 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
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

Descomposición en valores singulares y seudoinversa de Penrose 81
2.12.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 2.20 [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 2.21 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. (2.6)
Por otra parte, si hubiésemos resuelto el sistema a través de la seudoinversa

de Penrose habrı́amos obtenido
x = A+ b. (2.7)
por lo que comparando las ecuaciones (2.6) y (2.7) 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∗ .
2.13 Ejercicios propuestos

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 =
−1 + i   x2   1 + i 
    
3
    
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   1− α 
   
xn+1   xn
  
= +
     
     
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
      
x x x
  0 , 1 , 2 ... 
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||∞
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 2.12 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 2.13 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 2.14 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 2.15 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 2.16 Hallar la parábola de regresión de los puntos:
(1, 0), (0, 0), (−1, 0), (1, 2) y (2, 3)
Ejercicio 2.17 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 2.18 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 2.19 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 2.20 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 2.21 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 2.22 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.
Ejercicio 2.23 Sea el sistema Ax = b, donde

   
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 2.24 Se considera el sistema superdeterminado Ax = b con

   
3 2 3
   
A= 4 y b =
   
5   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 /121
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. Interpolación
3.1 Introducción
Supongamos que se conocen los n + 1 valores que toma una función f (x), en
los puntos del conjunto {x0 , x1 , . . . , xn } denominado soporte , es decir
f (x0 ) = y0 f (x1 ) = y1 ··· f (xn ) = yn
El problema de la interpolación consiste en encontrar una función g(x) con

determinadas caracterı́sticas y tal que g(xi ) = yi para i = 0, 1, . . . , n. En caso
de existir, se dice que g(x) interpola a f (x) en x0 , x1 , . . . , xn .
Al decir con determinadas caracterı́sticas nos referimos a que se exige que
g(x) sea, por ejemplo, un polinomio, un cociente de polinomios, una función
trigonométrica, etc.
La finalidad de encontrar una función g(x) que interpola a otra f (x) en los
puntos x0 , x1 , . . . , xn es la de aproximar la función f (x) en un punto x de tal
forma que se pueda decir que f (x) ≈ g(x) una vez encontrada g(x). (Otra
cosa es la evaluación de f (x) − g(x)).
Si los valores de x se encuentran en el intervalo [x0 , xn ] se dice que estamos
interpolando. Si se encuentran fuera de dicho intervalo, se dice que estamos
extrapolando.
Como aplicaciones más directas tenemos:
• Evaluación: (una aproximación) de una función complicada f , en un
cierto punto x.
• Si g(x) es cómoda de derivar o integrar, la sustitución, en cierta medida,

Z b Z b
0 0
de f por g o f por g.
a a
93
94 Interpolación
En este capı́tulo sólo trataremos la interpolación polinomial y la interpolación

polinomial a trozos (splines).
x 0 1 2
Ejemplo 3.1 Dada la tabla de valores
y 1 3 7
a) Dado que los tres puntos no están alineados, no existe ninguna recta que
interpole a dichos valores.
b) Si queremos calcular la parábola y = ax2 + bx + c que interpola a dichos

valores, planteando el correspondiente sistema se obtiene, como solución
única, y = x2 + x + 1.
c) Si nuestra intención es buscar una parábola cúbica y = ax3 + bx2 + cx + d

nos encontramos con que existen infinitas soluciones que son de la forma
y = x2 + x + 1 + α x(x − 1)(x − 2) ∀α ∈ R.
d) Por último, para calcular la función polinómica de grado n que interpola

a dichos valores obtenemos y = x2 + x + 1 + α xn1 (x − 1)n2 (x − 2)n3 para
cualesquiera n1 + n2 + n3 = n y cualquier α ∈ R.
3.2 Interpolación polinomial

Trataremos en esta sección los tres tipos más generalizados de interpolación
polinomial, a saber: Lagrange, Newton y Hermite.
3.2.1 Interpolación de Lagrange
Como en cualquier problema de interpolación, consideremos la tabla
x x0 x1 · · · xn
y y0 y1 · · · yn
y construyamos el polinomio de grado n que interpola a dichos valores. Para

ello, consideremos los denominados polinomios de Lagrange
Interpolación polinomial 95
(x − x1 )(x − x2 ) · · · (x − xn )
L0 (x) =
(x0 − x1 )(x0 − x2 ) · · · (x0 − xn )
(x − x0 )(x − x2 ) · · · (x − xn )
L1 (x) =
(x1 − x0 )(x1 − x2 ) · · · (x1 − xn )
..
.
n
Y x − xj
Li (x) =
xi − xj
j =0
j 6= i
..
.
(x − x0 )(x − x1 ) · · · (x − xn−1 )
Ln (x) =
(xn − x0 )(xn − x1 ) · · · (xn − xn−1 )
Teorema 3.1 Los polinomios de Lagrange, definidos anteriormente, verifi-

can:

 0 si i 6= j
a) Li (xj ) =
1 si i = j

b) gr(Li (x)) = n cualquiera que sea 0 ≤ i ≤ n.
c) El polinomio Pn (x) = y0 L0 (x) + y1 L1 (x) + · · · + yn Ln (x) interpola los

x x0 x1 · · · xn
valores de la tabla con x0 < x1 < · · · < xn , siendo
y y0 y1 · · · yn
gr(P (x)) ≤ n.
Demostración. Los dos primeros apartados son triviales. Para el tercero es

fácil observar que (ver el apartado 1) Pn (xi ) = yi y que Pn (x), al ser una
combinación lineal de polinomios de grado n (apartado 2), no puede ser de
grado superior a n (aunque sı́ inferior).
Ejemplo 3.2 Para interpolar los valores de la tabla
x 1 2 3 4
y 0 −1 2 −5
96 Interpolación
los polinomios de Lagrange son

(x − 2)(x − 3)(x − 4) 1
L0 (x) = = − (x3 − 9x2 + 26x − 24)
(1 − 2)(1 − 3)(1 − 4) 6
(x − 1)(x − 3)(x − 4) 1
L1 (x) = = (x3 − 8x2 + 19x − 12)
(2 − 1)(2 − 3)(2 − 4) 2
(x − 1)(x − 2)(x − 4) 1
L2 (x) = = − (x3 − 7x2 + 14x − 8)
(3 − 1)(3 − 2)(3 − 4) 2
(x − 1)(x − 2)(x − 3) 1
L3 (x) = = (x3 − 6x2 + 11x − 6)
(4 − 1)(4 − 2)(4 − 3) 6
y como P3 (x) = y0 · L0 (x) + y1 · L1 (x) + y2 · L2 (x) + y3 · L3 (x) obtenemos que

7 98
P3 (x) = − x3 + 16x2 − x + 19
3 3
El cálculo de los polinomios de Lagrange, puede verse con el Ejemplo 3.2, no

es un proceso dinámico, en el sentido de que si ahora añadiéramos un nuevo
punto al soporte, habrı́a que comenzar de nuevo todo el proceso.
Teorema 3.2 Dados los números reales x0 < x1 < · · · < xn y los n + 1
números reales cualesquiera y0 , y1 , . . . , yn , existe un único polinomio Pn (x) de
grado no superior a n tal que Pn (xi ) = yi para i = 0, 1, . . . , n.
Demostración. La existencia del polinomio es obvia por el Teorema 3.1. Para

probar la unicidad supongamos que existiesen dos polinomios Pn (x) y Qn (x)
de grados no superiores a n y tales que Pn (xi ) = Qn (xi ) = yi para 0 ≤ i ≤ n.
Consideremos el polinomio Dn (x) = Pn (x) − Qn (x).



 Dn (x0 ) = 0 =⇒ Dn (x) = (x − x0 )Dn−1 (x)

 Dn (x1 ) = 0 =⇒ Dn−1 (x) = (x − x1 )Dn−2 (x)




..
 .

Dn (xn−1 ) = 0 =⇒ D1 (x) = (x − xn−1 )D0 (x) = (x − xn−1 ) · k






 D (x ) = 0 =⇒ k = 0
n n
Se obtiene, por tanto, que

Dn (x) = (x − x0 )(x − x1 ) · · · (x − xn−1 ) · 0 = 0 =⇒ Pn (x) = Qn (x).
Dada una función f (x) de la que se conocen los transformados de n + 1

puntos x0 , x1 , . . . , xn y su polinomio de interpolación Pn (x), sólo nos falta dar
una medida del error que se comete al sustituir la función f (x) por el polinomio
Pn (x).
Teorema 3.3 Sean x0 < x1 < · · · < xn y sea f una función n + 1 veces
derivable tal que la derivada f (n + 1 (x) es continua.
Sean y0 = f (x0 ), y1 = f (x1 ), . . . , yn = f (xn ) y Pn (x) el polinomio de interpo-
x x0 x1 · · · xn
lación de los valores de la tabla y sea x un número real
y y0 y1 · · · yn
cualquiera. Se verifica que
f (n + 1 (c)
f (x) − Pn (x) = (x − x0 ) · · · (x − xn )
(n + 1)!
donde el punto c se encuentra en el intervalo determinado por los puntos

x, x0 , x1 , . . . , xn .
Demostración. Sean ε(x) = f (x) − Pn (x), z(x) = (x − x0 ) · · · (x − xn ) y x

z(t)
un punto cualquiera. Consideremos la función g(t) = ε(t) − ε(x)
z(x)
• g(xi ) = 0 cualquiera que sea 0 ≤ i ≤ n ya que
z(xi )
g(xi ) = ε(xi ) − ε(x) =0
z(x)
por serlo ε(xi ) y z(xi ).

z(x)
• g(x) = 0 ya que g(x) = ε(x) − ε(x) =0
z(x)
Por tanto, la función g(x) se anula en x, x0 , . . . , xn+1 , es decir, en n + 2

puntos. Si ordenamos estos puntos y los denotamos por t1 , t2 , . . . , tn+2 , en
cada intervalo (ti , ti+1 ) se tiene, por Rolle, que existe un punto ci tal que
g 0 (xi ) = 0 (ver Figura 3.1)
t1 t2 t3 tn+2
u u u u
c1 c2 c3 cn+1
Figura 3.1: Los n + 1 puntos en los que se anula g 0 (x).
98 Interpolación
Razonando de igual forma obtenemos n puntos donde se anula la derivada

segunda g 00 (x), n − 1 donde se anula la derivada tercera, y ası́ sucesivamente
hasta obtener un punto c en el que se anula la derivada de orden n + 1, siendo
c un punto del intervalo que determinan los puntos x, x0 , x1 , . . . , xn .
ε(x) (n + 1
Como g (n + 1 (t) = ε(n + 1 (t) − z (t), podemos particularizar en t = c y
z(x)
obtenemos
ε(x) (n + 1 ε(x)
0 = ε(n + 1 (c) − z (c) = f (n + 1 (c) − Pn(n + 1 (c) − (n + 1)!
z(x) z(x)
ε(x)
por lo que (n + 1)! = f (n + 1 (c), o lo que es lo mismo,
z(x)
f (n + 1 (c) f (n + 1 (c)
ε(x) = z(x) = (x − x0 ) · · · (x − xn )
(n + 1)! (n + 1)!
3.2.2 Interpolación de Newton
1.- Diferencias divididas

Consideremos una función f (x) y un soporte {x0 , x1 , . . . , xn } de n+1 puntos.
Denotemos por fi = f (xi ) y consideremos la tabla
x x0 x1 · · · xn
y f0 f1 · · · fn
Vamos a probar que el polinomio de grado no superior a n que interpola a

estos valores es de la forma
P (x) = c0 + c1 (x − x0 ) + c2 (x − x0 )(x − x1 ) + c3 (x − x0 )(x − x1 )(x − x2 ) + · · · +
+ · · · + cn (x − x0 )(x − x1 )(x − x2 ) · · · (x − xn−1 )
para después, calcular los valores de los coeficientes c0 , c1 , . . . , cn .
Proposición 3.1 Los coeficientes c0 , c1 , . . . , cn , descritos más arriba, depen-
den de los valores x0 , x1 , . . . , xn y f0 , f1 , . . . , fn .
Demostración. Como el polinomio P (x) interpola a la tabla de valores, se

tiene que
P (x0 ) = f0 =⇒ c0 = f0
f1 − f0
P (x1 ) = f1 =⇒ f1 = f0 + c1 (x1 − x0 ) =⇒ c1 =
x1 − x0
f1 − f0
P (x2 ) = f2 =⇒ f2 = f0 + (x2 − x0 ) + c2 (x2 − x0 )(x2 − x1 ) =⇒
x1 − x0
c2 depende de x0 , x1 , x2 , f0 , f1 y f2
Supuesto que ck−1 depende de x0 , . . . , xk−1 , f0 , . . . , fk−1 , se tiene que
P (xk ) = fk = c0 + c1 (xk − x0 ) + · · · + ck (xk − x0 ) · · · (xk − xk−1 ) por lo que
1 h i
ck = fk − c0 − · · · − ck−1 (xk − x0 ) · · · (xk − xk−2 )
(xk − x0 ) · · · (xk − xk−1 )
por lo que ck depende de x0 , x1 , . . . , xk , f0 , f1 , . . . , xk .
En lo que sigue utilizaremos la notación ck = f [x0 , x1 , . . . , xk ], con lo que el

polinomio quedará de la forma
P (x) = f [x0 ]+f [x0 , x1 ](x−x0 )+· · ·+f [x0 , x1 , . . . , xn ](x−x0 )(x−x1 ) · · · (x−xn )
y quedará determinado una vez que se determinen los valores de los coeficientes
f [x0 , x1 , . . . , xk ] con k = 0, 1, . . . , n.
Proposición 3.2 Sea P (x) el polinomio de interpolación correspondiente a la

x x0 x1 · · · xn
tabla . Si Q(x) y R(x) son los polinomios que interpolan
y f0 f1 · · · fn
x x0 x1 · · · xn−1 x x1 x2 · · · xn
las tablas y respectivamente, se
y f0 f1 · · · fn−1 y f1 f2 · · · fn
verifica que:
x − x0
P (x) = Q(x) + R(x) − Q(x)
xn − x0
Demostración. Consideremos el polinomio
x − x0
T (x) = Q(x) + R(x) − Q(x) .
xn − x0
Si probamos que T (x) es de grado no superior a n y que interpola la tabla
x x0 x1 · · · xn
, dado que dicho polinomio es único, se habrá probado
y f0 f1 · · · fn
que T (x) ≡ P (x) y, por tanto, nuestra proposición.
100 Interpolación
a) Es obvio que como los grados de Q(x) y R(x) no son superiores a n − 1,

T (x) no puede tener grado superior a n.
b) • T (x0 ) = Q(x0 ) = f0
• Si xi ∈ {x1 , x2 , . . . , xn−1 } se tiene que
xi − x0 xi − x0
T (xi ) = Q(xi ) + R(xi ) − Q(xi ) = fi + (fi − fi ) =
xn − x0 xn − x0
fi .
xn − x0
• T (xn ) = Q(xn ) + R(xn ) − Q(xn ) = R(xn ) = fn .
xn − x0
Proposición 3.3 Para cualquiera que sea k = 0, 1, . . . , n se verifica que

f [x1 , . . . , xk ] − f [x0 , . . . , xk−1 ]
f [x0 , x1 , . . . , xk ] =
xk − x0
siendo f [xi ] = fi para 0 ≤ i ≤ n.
Demostración. f [x0 , x1 , . . . , xk ] es el coeficiente ck del término de mayor

x x0 x1 · · · xk
grado del polinomio de interpolación de la tabla .
y f0 f1 · · · fk
Basta ir haciendo tomar a k los valores de 0 a n para que la Proposición 3.2
nos justifique el resultado.
Ejemplo 3.3 Calculemos, por el método de las diferencias divididas de New-

ton, el polinomio de interpolación de la tabla
x 1 3 4 5 7
y 0 1 −1 2 3
xi fi f [xi , xi+1 ] f [xi , xi+1 , xi+2 ] f [xi , xi+1 , xi+2 , xi+3 ] f [xi , xi+1 , xi+2 , xi+3 , xi+4 ]
1 0
1/2
3 1 − 5/6
− 2/1 5/6
4 −1 5/2 − 5/18
3/1 − 5/6
5 2 − 5/6
1/2
7 3
Por lo que el polinomio de interpolación es

1 5 5
P (x) = (x − 1) − (x − 1)(x − 3) + (x − 1)(x − 3)(x − 4)−
2 6 6
5
− (x − 1)(x − 3)(x − 4)(x − 5) =
18
1
= − (5x4 − 80x3 + 430x2 − 889x + 534).
18
La ventaja de este método es que si ahora introducimos un nuevo dato, por
x 9
ejemplo que f (9) = 5, es decir, el polinomio que se obtiene es
y 5
Q(x) = P (x) + f [x0 , x1 , x2 , x3 , x4 , x5 ](x − 1)(x − 3)(x − 4)(x − 5)(x − 7)
y tan sólo habrı́a que calcular el coeficiente f [x0 , x1 , x2 , x3 , x4 , x5 ] añadiendo
una nueva lı́nea a la tabla anterior.
x x0 · · · xn−1 xn
Puede observarse que dada la tabla , el polinomio de
y y0 · · · yn−1 yn
interpolación es de la forma
Pn (x) = Pn−1 (x) + f [x0 , . . . , xn ](x − x0 ) · · · (x − xn−1 )
Se tenı́a, también, que para dicho polinomio era

f (n (ξ)
f (x) − Pn−1 (x) = (x − x0 ) · · · (x − xn−1 )
n!
Sustituyendo x por xn tenemos:
f (n (ξ)
f (xn ) − Pn−1 (xn ) = (xn − x0 ) · · · (xn − xn−1 )
n!
y dado que f (xn ) = fn = Pn (xn ), se tiene que:
f (n (ξ)
Pn (xn ) − Pn−1 (xn ) = (xn − x0 ) · · · (xn − xn−1 ).
n!
Podemos, por tanto, enunciar la siguiente proposición.
x x0 · · · xn−1 xn
Proposición 3.1 Dada la tabla , con x0 < · · · < xn ,
y y0 · · · yn−1 yn
existe un punto c en el intervalo [x0 , xn ] para el que
f (n (c)
f [x0 , . . . , xn ] =
n!
102 Interpolación
2.- Diferencias finitas

x x0 x1 · · · xn
Consideremos la tabla en donde el soporte x0 , . . . , xn
f (x) f0 f1 · · · fn
es regular, es decir, en donde xi+1 − xi = cte.
Definición 3.1 Dados y0 , y1 , . . . , yn , se definen las diferencias finitas ∆k yi
como
∆yi = yi+1 − yi ∆k yi = ∆(∆k−1 yi )
Ası́, por ejemplo, para y0 , y1 , y2 , y3 se tendrı́an:
∆y0 = y1 − y0
∆2 y0 = ∆y1 − ∆y0
∆y1 = y2 − y1 ∆3 y0 = ∆2 y1 − ∆2 y0
∆2 y1 = ∆y2 − ∆y1
∆y2 = y3 − y2
x x0 x1 · · · xn
Proposición 3.4 Dada la tabla en la que x0 , . . . , xn
f (x) f0 f1 · · · fn
es un soporte regular con xi+1 − xi = h, se verifica que, para cualquier valor
de k = 1, 2, . . . , n, es:
∆k f0
f [x0 , . . . , xk ] = k
h · k!
Demostración. Haremos inducción en k. Para k = 1 sabemos que f [x0 , x1 ] =

f1 − f0 ∆f0
= . Supuesto cierto para k vamos a probarlo para k + 1.
x1 − x0 h · 1!
f [x1 , . . . , xk+1 ] − f [x0 , . . . , xk ] ∆k f1 − ∆k f0
f [x0 , . . . , xk+1 ] = = k =
xk+1 − x0 h · k! · k · h
∆k+1 f0
= .
hk+1 · (k + 1)!
El polinomio de interpolación del soporte x0 , x1 , . . . , xn es, por tanto:

∆2 f0 x − x0

x − x0 x − x1
Pn (x) = f0 + ∆f0 + + ···
h 2! h h
∆n f0 x − x0

x − xn−1
··· + ···
n! h h
Teniendo en cuenta que x − xk = x − (x0 + k · h) = (x − x0 ) − k · h, podemos

poner
∆2 f0 x − x0

x − x0 x − x0
Pn (x) = f0 + ∆f0 + − 1 + ···
h 2! h h
∆n f0 x − x0

x − x0
··· + ··· −k−1
n! h h
x − x0
por lo que, si denotamos por t = , se tiene que
h
∆f0 ∆2 f0 ∆n f0
Pn (x) = f0 + t+ t(t − 1) + · · · + t(t − 1) · · · (t − (n − 1))
1! 2! n!
! ! !
t t t
Es decir: Pn (x) = f0 + ∆f0 + ∆2 f0 + · · · + ∆n f0 donde
1 2 n
!
t t(t − 1) · · · (t − (k − 1))
= .
k k!
Fenómeno de Runge
Dada una función continua en [a, b], podrı́a pensarse que la sucesión Pn (x)
con n ∈ N de polinomios de interpolación, obtenidos al aumentar el número
de puntos del soporte, converge a la función f (x) es decir, podrı́amos pensar
que
lim |f (x) − Pn (x)| = 0 ∀ x ∈ [a, b]
n→∞
cosa que, sin embargo, no es cierta. En realidad, al aumentar el número de

puntos del soporte se mejora la aproximación en la parte central del intervalo,
pero la diferencia entre la función y el polinomio interpolador puede aumentar
rápidamente en los extremos. Ello nos dice que no es bueno hacer demasiado
extenso el soporte, ya que además de aumentar el número de operaciones con la
consecuente acumulación de errores, podemos aumentar la pérdida de precisión
en los extremos. Este fenómeno es conocido como fenómeno de Runge.
1
Ejemplo 3.4 Si aproximamos la función f (x) = por un polinomio de
1 + x2
segundo grado, en el soporte {−4, 0, 4}, obtenemos que P2 (x) = 1 − x2 /17. En
la Figura 3.2 (Izda.) podemos ver ambas gráficas.
104 Interpolación
Figura 3.2: Las funciones {f (x), P2 (x)} y {f (x), P4 (x)}
Si la aproximación la hacemos mediante un polinomio de grado 4 en el soporte

{−4, −2, 0, 2, 4} obtenemos
1
P4 (x) = (85 − 21x2 + x4 )
85
que podemos ver representada junto a la función f (x) en la Figura 3.2 (Dcha.).
Si afinamos aún más y aproximamos mediante un polinomio de grado 8 en
el soporte {−4, −3, −2, −1, 0, 1, 2, 3, 4} obtenemos
1
P8 (x) = (1700 − 1124x2 + 304x4 − 31x6 + x8 )
1700
cuya gráfica podemos observar en la Figura 3.3.
Figura 3.3: Las funciones f (x) y P8 (x)
Puede verse el hecho comentado anteriormente del fenómeno de Runge .

Vamos mejorando la aproximación en la parte central del intervalo, pero vamos
empeorándola en los extremos.
Ejemplo 3.5 Si construimos con MATHEMATICA el polinomio de interpo-

lación de la función Log(1 + x) en el soporte {0, 1, 2, 3, 4} y representamos el
resultado en el intervalo [0.5,1.5] obtenemos la gráfica de la Figura 3.4 (Izda.).
Si utilizamos los nueve puntos del soporte {0, 0.5, 1, 1.5, 2, 2.5, 3, 3.5, 4} y repre-
sentamos el polinomio obtenido en el intervalo [0.5, 1.5] obtenemos la gráfica
de la Figura 3.4 (Dcha.).
Figura 3.4: Las funciones {log(1 + x), P4 (x)} y {log(1 + x), P8 (x)}
Intercalando los puntos medios para obtener ahora un soporte de 17 pun-

tos y realizando la gráfica correspondiente obtenemos como resultado el que
se muestra en la Figura 3.5 (Izda.). Sin embargo, si volvemos a utilizar los
puntos medios para obtener un soporte de 33 puntos, podemos ver en la Fi-
gura 3.5 (Dcha.) que el fenómeno de Runge junto a la acumulación de errores
hace que el polinomio obtenido deje de ser una buena aproximación de la
función.
Figura 3.5: Las funciones {log(1 + x), P16 (x)} y {log(1 + x), P32 (x)}

106 Interpolación
3.2.3 Interpolación de Hermite

Este método consiste en buscar un polinomio que interpole a una función
f (x) en el soporte x0 , x1 , . . . , xn pero haciendo que coincidan, en los puntos
del soporte, no sólo los valores de la función con los del polinomio, sino que
también coincidan los valores de sus respectivas derivadas.
Consideremos, por tanto, la tabla
x x0 x1 · · · xn
f (x) f0 f1 · · · fn
f 0 (x) f00 f10 · · · fn0
donde fi = f (xi ) y fi0 = f 0 (xi ) para 0 ≤ i ≤ n.
Se tienen, en este caso, 2n + 2 condiciones, por lo que debemos buscar un
polinomio de grado 2n + 1
P (x) = a2n+1 x2n+1 + a2n x2n + · · · + a1 x + a0
que verifique las condiciones:
P (x0 ) = f0 P 0 (x0 ) = f00

P (x1 ) = f1 P 0 (x1 ) = f10
.. ..
. .
P (xn ) = fn P 0 (xn ) = fn0
Teorema 3.4 Dada la tabla

x x0 x1 · · · xn
f (x) f0 f1 · · · fn
f 0 (x) f00 f10 · · · fn0
sean Lk (x) (k = 0, 1, . . . , n los polinomios de Lagrange para el soporte dado.
X n h i
El polinomio P2n+1 (x) = ak + bk (x − xk ) L2k (x) en el que
k=0
ak = fk
bk = fk0 − 2fk L0k (xk )
verifica que )
P2n+1 (xk ) = fk
0
k = 0, 1, . . . , n
P2n+1 (xk ) = fk0
siendo, además, el único polinomio de grado 2n+1 que verifica las condiciones
anteriores.
La demostración del Teorema 3.4 debe realizarla el alumno a modo de ejercicio.
Los polinomios de Lagrange para el soporte x0 , x1 , . . . , xn son
(x − x1 )(x − x2 ) · · · (x − xn )
L0 (x) =
(x0 − x1 )(x0 − x2 ) · · · (x0 − xn )
(x − x0 )(x − x2 ) · · · (x − xn )
L1 (x) =
(x1 − x0 )(x1 − x2 ) · · · (x1 − xn )
..
.
n
Y x − xj
Li (x) =
xi − xj
j =0
j 6= i
..
.
(x − x0 )(x − x1 ) · · · (x − xn−1 )
Ln (x) =
(xn − x0 )(xn − x1 ) · · · (xn − xn−1 )
Llamando z(x) = (x − x0 )(x − x1 ) · · · (x − xn ) se tiene que

h i0
z 0 (x) = 1 · (x − x1 ) · · · (x − xn ) + (x − x0 ) (x − x1 ) · · · (x − xn )
por lo que
z 0 (x0 ) = (x0 − x1 ) · · · (x0 − xn )
y de manera análoga se obtiene que
z 0 (xk ) = (xk − x1 ) · · · (xk − xk−1 )(xk − xk+1 ) · · · (xk − xn )
por lo que los polinomios de Lagrange pueden escribirse de la forma
z(x)
Lk (x) =
(x − xk )z 0 (xk )
Teorema 3.5 Sea f (x) una función 2n + 2 veces derivable con derivada de
orden 2n+2 continua y sea P2n+1 el polinomio de Hermite que interpola a f (x)
en el soporte x0 , x1 , . . . , xn . Existe un punto c del intervalo que determinan
los puntos x, x0 , . . . , xn en el que se verifica que
f (2n + 2 (c)
f (x) − P2n+1 (x) = (x − x0 )2 · · · (x − xn )2
(2n + 2)!
108 Interpolación
Demostración. Sean ε(x) = f (x) − P2n+1 (x) y z(x) = (x − x0 ) · · · (x − xn ) y

considérese la función
ε(x)
g(t) = ε(t) − 2 z 2 (t)
z (x)
La demostración es similar a la del Teorema 3.3.
Ejemplo 3.6 Si aplicamos este método a la función del Ejercicio 3.4, en el

soporte {−4, −2, 0, 2, 4} obtenemos el polinomio de grado 8 (en realidad se
busca de grado 9 pero al ser una función par, el término de grado 9 se anula)
Figura 3.6: Las funciones {f (x), P8 (x)} y {f (x), P16 (x)}
1
P8 (x) = (7225 − 3129x2 + 569x4 − 41x6 + x8 )
7225
cuya gráfica puede verse en la Figura 3.6 (Izda.).
Si lo hacemos en el soporte {−4, −3, −2, −1, 0, 1, 2, 3, 4} obtenemos
1
P16 (x) = ( 2890000 − 2558224x2 + 1613584x4 − 626688x6 + 144408x8
2890000
−19527x10 + 1507x12 − 61x14 + x16 )
que podemos ver en la Figura 3.6 (Dcha.).

Si comparamos con los resultados obtenidos en el Ejercicio 3.4, podemos ob-
servar la mejora que produce la imposición de que coincidan no sólo los valores
de la función, sino que también lo hagan los de su derivada, en los puntos del
soporte. Sin embargo, sigue manifestándose el fenómeno de Runge, es decir,
se mejora el resultado en la parte central del intervalo, pero en los extremos,
la diferencia entre el polinomio interpolador y la función es considerable.
La manera de evitar el fenómeno de Runge es hacer una interpolación poli-

nomial a trozos, es decir, lo que se conoce como una interpolación por splines.
Interpolación por splines 109
3.3 Interpolación por splines

Consideremos una partición del intervalo [a, b]
∆ = {x0 = a < x1 < · · · < xn−1 < xn = b}
en la que los puntos xi reciben el nombre de nodos. Una interpolación por spli-
nes no es más que tomar un soporte en cada subintervalo [xi−1 , xi ] y construir
un polinomio de interpolación, de grado no superior a k (para un k prefijado)
sobre dicho soporte, por lo que el método se conoce también como interpo-
lación polinomial a trozos. Damos a continuación una definición formal de lo
que denominaremos función spline.
Definición 3.2 Una función spline de grado k con nodos en x0 , x1 , . . . , xn es
una función S(x) formada por varios polinomios, cada uno de ellos definido
sobre un subintervalo y que se unen entre sı́ bajo ciertas condiciones de conti-
nuidad. Las condiciones que debe cumplir S(x) son las siguientes:
a) En cada intervalo [xi−1 , xi ), S(x) es un polinomio de grado gr[S(x)] ≤ k,
b) S(x) admite derivada continua de orden k − 1 en [x0 , xn ].
En general, pueden crearse funciones spline de grado k, pero la interpolación

más frecuente es a través de funciones spline de grado 3, es decir, de splines
cúbicos.
3.3.1 Splines cúbicos

Dado que a partir de ahora vamos a trabajar con splines cúbicos, vamos a
concretar la Definición 3.2 al caso de k = 3.
Definición 3.3 Dado el conjunto de puntos ∆ = {x0 , x1 , . . . , xn }, diremos
que la función Sδ es un spline cúbico asociado a ∆ si cumple las siguientes
condiciones:
a) La restricción de S∆ a cada intervalo [xi−1 , xi ) para i = 1, 2, . . . n es

un polinomio de grado no superior a tres. Es decir, S∆|[xi−1 ,xi ] ∈ P3 [x],
donde P3 [x] representa al conjunto de los polinomios de grado menor o
igual a tres.
b) S∆ ∈ C 2 [a, b], es decir, S∆ es una función continua, dos veces derivable

y con derivadas continuas en el intervalo [a, b].
110 Interpolación
Definición 3.4 Diremos que S∆ (x) es un spline de interpolación en x según

la partición ∆, si
a) S∆ (x) es un spline cúbico asociado a ∆.
b) S∆ (xi ) = f (xi ) = yi para i = 0, 1, . . . , n, es decir, cumple las condiciones

de interpolación.
Antes de construir un spline cúbico vamos a ver cuántas condiciones ha de

cumplir y cuántas incógnitas van a hacernos falta. Si en cada intervalo de
la partición intentamos construir un polinomio de grado tres que aproxime a
la función, deberemos calcular cuatro incógnitas (los cuatro coeficientes del
polinomio de grado tres) por intervalo, es decir, 4n incógnitas. Por otro lado,
estos polinomios deben cumplir, en cada uno de los nodos, las condiciones:
S∆|[ xi−1 ,xi ] (xi ) = S∆|[ xi ,xi+1 ] (xi )

0 0
S∆| [ xi−1 ,xi ]
(xi ) = S∆| [ xi ,xi+1 ]
(xi ) i = 1, 2, . . . , n − 1 (3.1)
00 00
S∆| [ xi−1 ,xi ]
(xi ) = S∆| [ xi ,xi+1 ]
(xi )
Es decir, se deben cumplir un total de 3(n − 1) condiciones además de las n + 1

condiciones de interpolación
S∆ (xi ) = f (xi ) i = 0, 1, . . . , n
Dado que tenemos un total de 4n incógnitas para 4n − 2 condiciones, debemos

imponer dos nuevas condiciones para poder determinar los coeficientes de la
función spline. Dependiendo de las condiciones que impongamos, obtendremos
un tipo de spline u otro.
• Si exigimos que las derivadas segundas se anulen en los extremos, es

decir, si
00 00
S∆ (a) = S∆ (b) = 0
diremos que S∆ (x) es el spline natural asociado a la partición ∆.
• Si exigimos que
0 0 00 00
S∆ (a) = S∆ (b), S∆ (a) = S∆ (b)
diremos que se trata de un spline periódico .

Interpolación por splines 111
3.3.2 Cálculo de los splines cúbicos de interpolación

Nos centraremos en el cálculo de los splines naturales y con al fin de simpli-
ficar la notación, llamaremos
hi = xi − xi−1 i = 1, 2, . . . , n
00
M i = S∆ (xi ) i = 0, 1, . . . , n
Los valores Mi se denominan momentos y determinarán completamente los

splines cúbicos.
Obsérvese, en primer lugar, que como en cada intervalo [xi , xi+1 ] el spline S∆
es un polinomio de grado tres, su segunda derivada es una recta (un polinomio
de grado uno). En consecuencia, al imponer las condiciones (3.1) sobre la
igualdad de las derivadas segundas en los nodos, obligamos a que la segunda
00
derivada de la función spline S∆ (x) constituya un conjunto de rectas que se
intersecan en los nodos de la partición elegida. Ahora bien, dado que cada
recta queda determinado por dos puntos, podemos escribir el valor de las
00
restricciones (3.1) sobre S∆| [x ,x ]
como
i i+1
00 xi+1 − x x − xi
S∆| (x) = Mi + Mi+1 .
[x ,x
i i+1 ] hi+1 hi+1
Integrando respecto a x obtenemos el valor de la primera derivada del spline
en este intervalo
0 Mi (xi+1 − x)2 Mi+1 (x − xi )2

S∆| (x) = − + + Ai .
[x ,x
i i+1 ]
2 hi+1 2 hi+1
Volviendo a integrar respecto a x obtenemos
Mi (xi+1 − x)3 Mi+1 (x − xi )3

S∆|[xi ,xi+1 ] (x) = + + Ai (x − xi ) + Bi .
6 hi+1 6 hi+1
Si imponemos ahora las condiciones de interpolación
S∆ (xi ) = yi , S∆ (xi+1 ) = yi+1
obtenemos
Mi 2 Mi 2
hi+i + Bi = yi =⇒ Bi = yi − h
6 6 i+1
Mi+1 2 yi+1 − yi hi+1
hi+i + Ai hi+1 + Bi = yi+1 =⇒ Ai = − (Mi+1 − Mi ).
6 hi+1 6
112 Interpolación
Podemos, por tanto, determinar los valores de las constantes Ai y Bi , que

determinan el valor de S∆ (x) en el intervalo [xi , xi+1 ], en función de los mo-
mentos.
El problema se reduce, por tanto, a calcular los momentos para cada uno
de los intervalos, para lo que utilizaremos la única condición de (3.1) que no
hemos utilizado:
0 0
S∆| [x ,x ]
(xi ) = S∆| [x ,x ]
(xi ).
i−1 i i i+1
Esta condición nos da, para cada i = 1, 2, . . . , n − 1, una ecuación:

hi hi+1 6 yi+1 − yi yi − yi−1
Mi−1 + 2Mi + Mi+1 = −
hi + hi+1 hi + hi+1 hi + hi+1 hi+1 hi
En el caso del spline natural tenemos que M0 = M − n = 0, quedándonos el
sistema tridiagonal de n − 1 ecuaciones con n − 1 incógnitas
 
h2
2
h1 + h2  
M1

h2 h3
 
2
 
  M 
 h2 + h3 h2 + h3  2 
... ... ... .. 
  



 . =
 . 

 ... ... h n−1   .. 
 

 hn−2 − hn−1 
 Mn−1
 
 hn−1
2
hn−1 + hn
 
6 y2 − y1 y1 − y0
−

 h1 + h2 h2 h1 

 6 y 3 − y 2 y 2 − y 1

 − 

 h2 + h3 h3 h2 

= .
.. 

..
 
 
 . 

yn − yn−1 yn−1 − yn−2 

 6
−
hn−1 + hn hn hn−1
Este sistema puede resolverse por cualquiera de los métodos iterados estudia-
dos en el Capı́tulo 2 ya que, al ser la matriz del sistema de diagonal dominante,
todos ellos son convergentes.
Ejemplo 3.7 Si aplicamos le interpolación por splines cúbicos a la función
del Ejemplo 3.4
1
f (x) = en la partición ∆ = {−4, −3, −2, −1, 0, 1, 2, 3, 4}
1 + x2
Ejercicios 113
Figura 3.7: Las funciones f (x) y S∆ (x)
obtenemos, utilizando MATHEMATICA, el resultado de la Figura 3.7, que

puede verse que, independientemente de ser mejor que el que se obtuvo en
la Figura 3.6 (Dcha.) con el método de Hermite, no aparece el fenómeno de
Runge.
3.4 Ejercicios
Ejercicio 3.1 Calcular los polinomios de Lagrange para el soporte canónico
con 1 ≤ n ≤ 3.
Ejercicio 3.2 Hallar el polinomio de interpolación de la función f (x) = 2x4

en el soporte canónico S = {0, 1, 2, 3}. Obtener una expresión del error.
Ejercicio 3.3 Hallar el polinomio de interpolación de la función f (x) = ex en

√
el soporte {0, 1k} y con él, aproximar e estimando el error cometido.
Ejercicio 3.4 Obtener el polinomio de interpolación de los puntos:
(0, −5), (1, −3), (2, 1), (3, 13)
a) Mediante resolución de un sistema de ecuaciones.
b) Mediante la fórmula de Lagrange
c) Mediante la fórmula de Newton para diferencias divididas.
d) Mediante la fórmula de Newton para diferencias finitas.

114 Interpolación
Ejercicio 3.5 Probar que F (n) = 12 + 22 + 32 + · · · + n2 es un polinomio en

n y obtenerlo por interpolación.
Ejercicio 3.6 Obtener el polinomio de interpolación de Hermite de la función

f (x) = ln x en el soporte S = {1, 2} y, supuesto conocido ln 2, aproximar el
valor de ln 10 5 acotando el error cometido.
Ejercicio 3.7 Dada la función f (x) = ex , se pide:
a) Hallar el polinomio de interpolación en el soporte {−1, 0, 1} y una cota

del error en el intervalo [−1, 1].
Calcular P (00 01) y compararlo con el valor dado por la calculadora para
0
e0 01 .
b) Hallar el polinomio de segundo grado, de mejor aproximación, en el

intervalo [−1, 1], respecto a la norma inducida por el producto escalar
Z 1
< f, g >= f (x)g(x)dx.
−1
0
Calcular el error cometido en la aproximación de e0 01 y comparar su
valor con el que nos da la calculadora.
x 0 1 2 3
Ejercicio 3.8 Dada la tabla se pide:
y −1 6 31 98
a) Hallar su polinomio de interpolación por el método de los polinomios de

Lagrange.
b) Determinar la forma general de todos los polinomios de cuarto grado que

satisfacen dicha tabla, determinando aquel que verifica que para x = 4
es y = 255.
c) Determinar los polinomios anteriores (para los soportes {0, 1, 2, 3} y

{0, 1, 2, 3, 4}) por el método de las diferencias divididas de Newton.
4. Integración numérica
4.1 Introducción
Se pretende,
Z b en este tema, dar una aproximación numérica del valor de una
integral f (x)dx en los distintos problemas que se presentan en la práctica,
a
como son:
a) Conocida una primitiva F (x) de la función f (x) sabemos que
Z b
f (x)dx = F (b) − F (a)
a
pero necesitamos aproximar el valor de F (b) − F (a).

Ası́, por ejemplo,
Z 2
1 h i2
dx = log x = log 2 − log 1 = log 2
1 x 1
pero hay que aproximar el valor de log 2.
b) Si se conoce la función f (x), pero no se conoce ninguna primitiva suya,

se busca otra función g(x) que aproxime a la función f (x) y de la cual
sı́ se conozcan primitivas.
Z 2 x
e
Ası́, por ejemplo, para calcular dx, se desarrolla en serie de poten-
1 x
cias
xn
e x 1 + x + · · · + n−1
f (x) = = n! + ε(x) = 1 + 1 + · · · + x + ε(x) =
x x x n!
= g(x) + ε(x)
Z 2 Z 2 Z 2
f (x)dx = g(x)dx + ε(x)dx
1 1 1
115
116 Integración numérica
Z 2
en donde habrá que evaluar ε(x)dx.
1
c) Sólo se conocen los valores de f (x) en un soporte {x0 , x1 , . . . , xn }.

En éste caso, se interpola la función (por ejemplo mediante la interpo-
lación polinómica).
Z b Z b Z b
f (n + 1 (c)
f (x)dx = Pn (x)dx + (x − x0 )(x − x1 ) · · · (x − xn )dx
a a a (n + 1)!
Z b Z b Z b
f (x)dx = Pn (x)dx + ε(x)dx
a a a
4.2 Fórmulas de cuadratura

Si realizamos la interpolación de Lagrange, y llamamos
z(x) = (x − x0 )(x − x1 ) · · · (x − xn ),
el polinomio de interpolación es
Pn (x) = y0 L0 (x) + y1 L1 (x) + · · · + yn Ln (x)
en donde los polinomios de Lagrange Li (x) pueden expresarse de la forma
z(x)
Li (x) =
(x − xi ) z 0 (xi )
Además,
Z b n Z
X b n
X Z b n
X
Pn (x) = yi Li (x)dx = yi Li (x)dx = ai y i
a i=0 a i=0 a i=0
Rb
donde los coeficientes ai = a
Li (x)dx no dependen de la función, sino sólo del
soporte.
Además, si f es un polinomio de grado no superior a n, ε(x) = 0, por lo que
para polinomios es
Z b n
X
P (x)dx = ai P (xi )
a i=1
Fórmulas de cuadratura 117
Por tanto:

P (x) = 1 =⇒ b − a = a0 + a1 + · · · + an 


b 2 − a2



P (x) = x =⇒ = a0 x 0 + a1 x 1 + · · · + an x n


2



....................................................... (4.1)


.......................................................  


n+1 n+1
b − a


n n n n
= a0 x 0 + a1 x 1 + · · · + an x n 

P (x) = x =⇒ 
n+1
sistema que, en forma matricial es
    
1 1 ··· 1 a0 b−a
    
 x x ··· x  b2 −a2
0 1 n a1  
2

=
    
 . .. . . ..  .. ..
 ..
 
 . . . 
 .  
  . 

n n n an+1 −bn+1
x0 x1 · · · xn an n+1
cuyo determinante es un Vandermonde.
Una vez calculados los coeficientes ai se obtiene una fórmula de aproximación

que sólo dependerá del soporte. Para cada soporte, las fórmulas reciben el
nombre de fórmulas de cuadratura .
Ejemplo 4.1 Vamos a integrar una función f (x) en [0, 1] considerando los
soportes:
1 1 1 1 3
S1 = {0, , } y S2 = { , , }.
3 2 4 2 4
1 1
a) En el soporte S1 = {0, , }
3 2
Z 1
1 1
f (x)dx ' a0 f (0) + a1 f ( ) + a2 f ( )
0 3 2
El sistema a resolver es, en este caso:
P (x) = 1 =⇒ a0 + a1 + a2 = 1
1 1 1
P (x) = x =⇒ 0 · a0 + a1 + a2 =
3 2 2
1 1 1
P (x) = x2 =⇒ 0 · a0 + a1 + a2 =
9 4 3
cuya solución es
1 3
a0 = a1 = − a2 = 2
2 2
luego Z 1
1 3 1 1
f (x)dx ' f (0) − f ( ) + 2f ( )
0 2 2 3 2
1 1 3
b) En el soporte S2 = { , , }
4 2 4
Z 1
1 1 3
f (x)dx ' b0 f ( ) + b1 f ( ) + b2 f ( )
0 4 2 4
El sistema a resolver es, en este caso:
P (x) = 1 =⇒ b0 + b1 + b2 = 1
1 1 3 1
P (x) = x =⇒ b0 + b1 + b2 =
4 2 4 2
1 1 9 1
P (x) = x2 =⇒ b0 + b1 + b2 =
16 4 16 3
cuya solución es
2 1 2
b0 = b1 = − b2 =
3 3 3
luego Z 1
2 1 1 1 2 3
f (x)dx ' f ( ) − f ( ) + f ( )
0 3 4 3 2 3 4
Las formas más generalizadas de aproximación de la integral de una función

f (x) se realizan mediante uno de los dos procesos siguientes:
a) Dando un soporte (generalmente regular) y los valores de la función en

los puntos del soporte. Fórmulas de Newton-Cötes .
b) Dando diferentes soportes y buscando el polinomio P (x) que hace más

Z b
pequeña la integral (f (x)−P (x))dx. Fórmulas de Gauss que no se
a
verán en este curso.
Fórmulas de Newton-Cötes 119
4.3 Fórmulas de Newton-Cötes

Partamos del soporte regular {x0 , x1 , . . . , xn } con
x0 = a x1 = a + h ··· xi = a + ih ··· xn = a + nh = b
Si llamamos z(x) = (x − x0 )(x − x1 ) · · · (x − xn ) se tiene que los polinomios

de Lagrange son
(x − x0 )(x − x1 ) · · · (x − xn )
Li (x) = =
(x − xi )(xi − x0 ) · · · (xi − xi−1 )(xi − xin1 ) · · · (xi − xn )
(x − a)(x − a − h) · · · (x − a − nh)
= =
(x − a − ih) ih (i − 1)h · · · h (−h) (−2h) · · · (−(n − i)h)
(x − a)(x − a − h) · · · (x − a − nh)
= =
(x − a − ih) · i! · (n − i)! · hn−1 · (−1)n−i
x − a x − a x − a
− 1 ··· −n
=h h h h =
(x − a − ih) · i! · (n − 1)! · (−1)n−i
x − a x − a x − a
− 1 ··· −n
= h h h
x−a
h
− i · i! · (n − 1)! · (−1)n−i
x−a
Por lo que haciendo t = se tiene que
h
t(t − 1) · · · (t − n)
Li (x) =
(t − i) · i! · (n − i)! · (−1)n−i
Por tanto:
Z n n
t(t − 1) · · · (t − n)
Z
ai = Li (x)dx = hdt =
0 0 (t − i) · i! · (n − i)! · (−1)n−i
n−i n
t(t − 1) · · · (t − n)
Z
(−1) h
= dt =⇒
i! · (n − i)! 0 t−i

n
Z n
n−i
i z(t)
ai = h(−1) dt
n! 0 t−i
que son los coeficientes de Newton-Cötes .

Proposición 4.1 Los coeficientes de Cötes verifican que ak = an−k .
Demostración.

n
n
n−k
Z
k z(t)
an−k = h(−1) dt
n! 0 t − (n − k)

 t − n = −u
Haciendo se tiene que
 dt = −du

n
Z n
k+1
n−k (−u + n)(−u + n − 1) · · · (−u)
an−k = h(−1) du =
n! 0 −u + k

n
Z n
k+1+n+1
n−k u(u − 1) · · · (u − n)
= h(−1) du =
n! 0 u−k

n
Z n
n+k
k u(u − 1) · · · (u − n)
= h(−1) du =
(−1)2k n! 0 u−k

n
Z n
n−k
k z(u)
= h(−1) du = ak
n! 0 u−k
Teniendo en cuanta la Proposición 4.1, sólo hay que calcular la mitad de los
coeficientes.
Las Fórmulas de Newton-Cötes en los casos n = 1 y n = 2 son conocidas

como Fórmula del trapecio y Fórmula de Simpson respectivamente.
4.3.1 Fórmula del trapecio

La fórmula de Newton-Cötes en el caso n = 1 sólo tiene dos coeficientes.
Como por la Proposición 4.1 es a0 = a1 y por las ecuaciones (4.1) es ao + a1 =
b − a, se tiene que
b−a
a0 = a1 =
2
por lo que
Z b
b−a b−a f (a) + f (b)
f (x)dx = f (a) + f (b) = (b − a)
a 2 2 2
Es decir, el método del trapecio nos aproxima la integral por el área de la

Fórmulas de Newton-Cötes 121
f (x)
f (b)
f (a)
a b
Figura 4.1: Método del trapecio
región plana limitada por las rectas x = a, x = b, y = 0 y la recta que pasa por
los puntos (a, f (a)) y (b, f (b)), es decir, el área de un trapecio (ver Figura 4.1).
4.3.2 Fórmula de Simpson

a+b
Para el caso n = 2 tenemos que x0 = a, x1 = y x2 = b.
2
Z
n 2
n−0 0 t(t − 1)(t − 2) b − a h t3 t2 i2 b − a
a0 = h(−1) dt = − 3 + 2t =
n! 0 t−0 4 3 2 0 6
b−a
a2 = a0 =
6
b−a 2(b − a)
Dado que a0 + a1 + a2 = b − a, se tiene que a1 = b − a − 2 = .
6 3
Se tiene, por tanto, que
Z b
b−a 2(b − a) a + b b−a
f (x)dx = f (a) + f( )+ f (b)
a 6 3 2 6
o, lo que es lo mismo:
Z b
b − ah a+b i
f (x)dx = f (a) + 4f ( ) + f (b)
a 6 2
Teorema 4.1 Al aplicar la fórmula de Newton-Cötes para un entero n, el

error que se comete viene dado por:
a) Si n es par:
n
hn+3 f (n + 2 (c)
Z
εn = t · t(t − 1) · · · (t − n)dt
(n + 2)! 0
b) Si n es impar:
n
hn+2 f (n + 1 (c)
Z
εn = (t − 1) · · · (t − n)dt
(n + 1)! 0
Corolario 4.2 El error cometido en la aproximación numérica de una integral

es:
a) Para la fórmula del trapecio:
h3 · f 00 (c)
ε=−
12
b) Para el método de Simpson:
h5 · f (iv (c)
ε=−
90
4.4 Fórmulas compuestas
4.4.1 Simpson para n par

Descomponiendo el soporte en {x0 , x1 , x2 }∪{x2 , x3 , x4 }∪· · ·∪{xn−2 , xn−1 , xn }
se obtiene que
Z b
x2 − x0 h i
f (x)dx ' f (x0 ) + 4f (x1 ) + f (x2 ) +
a 6
x4 − x2 h i
+ f (x2 ) + 4f (x3 ) + f (x4 ) + · · · +
6
xn − xn−2 h i
+ ··· + f (xn−2 ) + 4f (xn−1 ) + f (xn ) =
6
hh
= f (x0 ) + 4f (x1 ) + 2f (x2 ) + 4f (x3 ) + 2f (4 ) + · · · +
3
i
+ · · · + 2f (xn−2 ) + 4f (xn−1 ) + f (xn ) =

h
= f (x0 )+f (xn ) +2 f (x2 )+· · ·+f (xn−2 ) +4 f (x3 )+· · ·+f (xn−1 )
3
El error viene dado por

Ejercicios 123
h5 h (iv (iv (iv

i h5 n
|ε| = f (c0 ) + f (c2 ) + · · · + f (c 2 ) ≤
n máx |f (iv (x)| · =⇒
90 90 x∈[x0 ,xn ] 2
(b − a)5
|ε| ≤ máx |f (iv (x)|
180 n4 x∈[x0 ,xn ]
4.4.2 Trapecios para n impar

Con un proceso análogo al anterior obtenemos que
Z b
h
f (x)dx ' f (x0 ) + f (x1 ) + f (x1 ) + f (x2 ) + · · · +
a 2

+ · · · + f (xn−1 + f (xn ) =⇒
b
b−a
Z
f (x)dx ' f (x0 ) + f (xn ) + 2 f (x1 ) + f (x2 ) + · · · + f (xn−1 )
a 2n
El error que se comete viene dado por

h3 h 00 i h3
|ε| ≤ f (c1 ) + f 00 (c2 ) + · · · + f 00 (cn ) ≤ · n · máx |f 00 (x)| =⇒
12 12 x∈[x0 ,xn ]
(b − a)3
|ε| = máx |f 00 (x)|
12 n2 x∈[x0 ,xn ]
4.5 Ejercicios
Ejercicio 4.1 Probar que los coeficientes ak de las fórmulas de Newton-Cötes
verifican: n
X (−1)k ak

n
= 0.
k=0 k
1
1 − x2
Z
Ejercicio 4.2 Dada la integral dx se pide:
0 1 + x2
a) Calcularla exactamente.
b) Calcularla, aproximadamente, por la fórmula básica de Simpson.
c) Calcularla por la fórmula compuesta de Simpson de 11 sumandos.

d) Aplicar la siguiente fórmula:

Z 1
1h p p i
f (x) dx ' 5f (− /5) + 8f (0) + 5f ( 3/5)
3
−1 9
comprobando que integra, exactamente, polinomios de grado menor o

igual que 5.
Ejercicio 4.3 Se considera el soporte {−1, c, 1} donde c ∈ (−1, 1) es fijo. Sea

f (x) ∈ C[−1, 1]
a) Obtener el polinomio de interpolación de f (x) y una expresión del error.
b) Determinar los coeficientes a0 , a1 , a2 en la fórmula de cuadratura

Z 1
f (x) dx ' a0 f (−1) + a1 f (c) + a2 f (1)
−1
para que integre, exactamente, polinomios del mayor grado posible.
c) Dar una condición, necesaria y suficiente, para que dicha fórmula sea
exacta para polinomios de tercer grado.
r
5x + 13
d) Aplicar la fórmula a f (x) = con c = 00 1 y comparar con el
2
valor exacto.
Z 1
Ejercicio 4.4 Calcular f (x) ln x dx interpolando f (x), por un polinomio
0
de tercer grado, en el soporte {0, 1/3, 2/3, 1} y aplicar el resultado al cálculo de
Z 1
sen x ln x dx.
0 Z 1
−1
Ayuda: xm ln x dx = (m ≥ 0)
0 (m + 1)2
Ejercicio 4.5 Se considera una fórmula de cuadratura del tipo:

Z 1
f (x) dx = a0 f (− 1/2) + a1 f (0) + a2 f ( 1/2) + E
−1
a) Determinar los coeficientes a0 , a1 y a2 , para que sea exacta para polino-

mios del mayor grado posible.
Ejercicios 125
b) Sea P3 (x) el polinomio de interpolación “mixta” tal que:
P3 (− 1/2) = f (− 1/2), P3 (0) = f (0), P30 (0) = f 0 (0) y P3 ( 1/2) = f ( 1/2).
Razonar que la fórmula de cuadratura obtenida en el apartado a) es la

misma que se obtendrı́a integrando el polinomio P3 (x).
c) Como consecuencia, deducir que el error de la fórmula puede expresarse

de la forma:
Z 1 (IV
f (cx ) 2 1 1
E= x (x − )(x + ) dx.
−1 4! 2 2
Ejercicio 4.6 Probar que la fórmula compuesta de los trapecios para el in-
tervalo [0, 2π]:
Z 2π
h
f (x) dx = [f (0) + 2f (h) + 2f (2h) + · · · + 2f ((n − 1)h) + f (2π)] + E
0 2
(h = 2π/n) integra, exactamente, las funciones:
1, sin x, cos x, sin 2x, cos 2x, . . . , sin(n − 1)x, cos(n − 1)x.
Ejercicio 4.7 Determinar el número de sumandos necesarios, en las fórmulas

compuestas de los trapecios y Simpson, para calcular, con seis cifras decimales
exactas, las siguientes integrales:
Z 2 Z 3 x
e
a) I = ln x dx. b) I = dx.
1 2 x
Z 1
Ejercicio 4.8 Se considera la integral ex (4 − x) dx :
0
a) Calcularla exactamente (se supone conocido el número e).
b) Determinar el número mı́nimo de sumandos necesarios, en la fórmula

compuesta de Simpson, para que el error de discretización sea menor
que 10−m con m = 2, 3, 4, 5 y 6.
c) Calcular la integral, por la fórmula compuesta de Simpson, con cuatro

cifras decimales exactas.
Ejercicio 4.9 El recinto de la figura adjunta, que se encuentra inmerso en

una cuadrı́cula, está limitado por una recta y una curva de la que se conoce
que se trata de un polinomio de cuarto grado.
a) Calcular el área exacta del recinto sin determinar el polinomio que la

delimita.
b) Determinar, por el método de las diferencias divididas, el polinomio que

la delimita y comprobar que el área calculada en el apartado anterior
coincide con la que se obtiene por integración directa del polinomio.
Índice
Algoritmo Ecuaciones
de Horner, 22, 23 algebraicas, 2
de la bisección, 6 Espacio normado, 35
de Newton, 13
Factorización
Bisección de Cholesky, 52
algoritmo de la, 6 LU, 49
método de la, 4 ortogonal, 64
Bolzano Fórmula
teorema de, 4 de cuadratura, 117
de Gauss, 118
Ceros de Heron, 14
de polinomios, 19 de Newton-Cötes, 118
de una función, 1 de Newton-Raphson, 11
Cholesky de Simpson, 121
factorización de, 52 para n par, 122
Coeficientes del trapecio, 120
de Newton-Cötes, 119 para n impar, 123
Condición Fourier
número de, 43, 44 regla de, 14, 16
Convergencia, 37 Función
Descenso más rápido ceros de una, 1
método de, 61 contractiva, 7
Descomposición spline, 109
en valores singulares, 79 Gauss
Descomposición fórmulas de, 118
método de, 59 Gauss-Seidel
Diferencias finitas, 102 método de, 60
Distancia, 36 Gradiente conjugado
inducida, 36 método de, 61
127
128 Índice
Hadamard consistente, 55
martiz de, 51 convergente, 55
Heron de descomposición, 59
fórmula de, 14 de Gauss-Seidel, 60
Horner de Jacobi, 59
algoritmo de, 22, 23 de la bisección, 4
Householder de Newton, 11
transformación de, 65 para raı́ces múltiples, 17
transformación en el campo com- de relajación, 60
plejo, 68 de Sturm, 4
del descenso más rápido, 61
Interpolación, 93 del gradiente conjugado, 61
de Lagrange, 94 directo, 1, 42, 48
de Newton iterado, 1, 42, 55
diferencias divididas, 98
diferencias finitas, 102 Newton
polinomial, 94 algoritmo de, 13
por splines, 94, 109 interpolación de
spline de, 110 diferencias divididas, 98
diferencias finitas, 102
Jacobi método de, 11
método de, 59 para raı́ces múltiples, 17
Lagrange Newton-Cötes
interpolación de, 94 coeficientes de, 119
polinomios de, 94 fórmulas de, 118
Laguerre Newton-Raphson
regla de, 4 fórmula de, 11
Nodos, 109
Matriz Norma, 35
de diagonal dominante, 51 euclı́dea, 36
de Hadamard, 51 infinito, 36
fundamental, 50 matricial, 37
sparce, 42 euclı́dea, 38
triangular infinito, 38
inferior, 42 uno, 38
superior, 42 multiplicativa, 35
tridiagonal, 42 uno, 36
unitaria, 39 vectorial, 36
Momentos, 111 Número de condición, 43, 44
Método
Índice 129
Penrose mal condicionado, 43

seudoinversa de, 81 superdeterminado, 75
Pivote, 50 Soporte, 93
Polinomios regular, 102
de Lagrange, 94 Spline
Punto fijo cúbico, 109
teorema del, 7 de interpolación, 110
función, 109
Radio espectral, 40 natural, 110
Raı́ces periódico, 110
acotación de, 3, 4 Sturm
de una ecuación, 1 método de, 4
múltiples, 1 sucesión de, 19
separación de, 4 Sucesión
simples, 1 de Sturm, 19
Regla
de Fourier, 14, 16 Teorema
de Laguerre, 4 de Bolzano, 4
de Ruffini, 23 de Rolle, 5
Relajación de Rouche-Fröbenius, 75
método de, 60 del punto fijo, 7
Rolle Fundamental del Álgebra, 2
teorema de, 5 Transformación
Rouche-Fröbenius de Householder, 65
teorema de, 75 en el campo complejo, 68
Ruffini unitaria, 39
regla de, 23 Trapecio
Runge fórmula del, 120
fenómeno de, 103, 104 para n impar, 123
Seudoinversa
de Penrose, 79, 81
Seudosolución, 77
Simpson
fórmula de, 121
para n par, 122
Sistema
bien condicionado, 43
compatible
determinado, 43

Apuntes Calculo Numerico

Cargado por

Copyright:

Formatos disponibles

Apuntes Calculo Numerico

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

Apuntes Calculo Numerico

Cargado por

Copyright:

Formatos disponibles

E.T.S.

INGENIERÍAS TÉCNICAS EN INFORMÁTICA

Fco. Javier Cobos Gavala

2 Sistemas de ecuaciones lineales 35

2.1.4 Normas matriciales . . . . . . . . . . . . . . . . . . . . 37

4 Integración numérica 115

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

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-

donde {x1 , x2 , . . . , xn } son los ceros del polinomio.

Pn (x) = a0 (x − x1 )m1 (x − x2 )m2 · · · (x − xk )mk

donde mi (i = 1, 2, . . . k) representa la multiplicidad de la raı́z xi y verificándo-

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

D(x) = mcd [P (x), P 0 (x)] = (x − x1 )m1 −1 · · · (x − xk )mk −1

Definición 1.2 Dada una ecuación f (x) = 0

Geométricamente consiste en determinar una

Proposición 1.1 Si x es una raı́z de la ecuación

P (x) = a0 xn + a1 xn−1 + · · · + an−1 x + an = 0

Demostración. Sea P (x) = a0 xn + a1 xn−1 + · · · + an−1 x + an . Tomando

≥ |a0 xn | − A |x|n−1 + · · · + |x| + 1

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

Proposición 1.1 [Regla de Laguerre] Consideremos la ecuación

El procedimiento consiste en comenzar con la cota obtenida anteriormente

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

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

El error cometido, tomando como raı́z de la ecuación el punto medio del

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

Vamos a ver a continuación otros métodos que reducen, de forma considera-

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

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”

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).

podemos garantizar que partiendo de x0 = 1 el método convergerá a la raı́z

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

1.4 Método de Newton: análisis de errores

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

Este método es también conocido como mé-

Supuesta la existencia de dicho intervalo [a, b], el valor de k es independiente

Ejemplo 1.3 En el Ejemplo 1.2 calculamos la raı́z de 3 con 14 cifras decimales

Newton tomando como valor inicial x0 = 2. (Más adelante veremos porqué

Puede observarse cómo la convergencia del método de Newton-Raphson es

1.4.2 Regla de Fourier