Mns
Mns
Mns
[email protected] [email protected]
www.hectormora.info
1 Preliminares 1
1.1 Notación . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.2 Repaso de algunos conceptos de cálculo . . . . . . . . . . . . 2
1.3 Sucesiones . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
1.4 Polinomio de Taylor . . . . . . . . . . . . . . . . . . . . . . . 7
1.5 Notación O grande . . . . . . . . . . . . . . . . . . . . . . . . 11
1.6 Orden de convergencia . . . . . . . . . . . . . . . . . . . . . . 12
1.7 Números en un computador . . . . . . . . . . . . . . . . . . . 15
1.8 Truncamiento y redondeo . . . . . . . . . . . . . . . . . . . . 17
1.9 Error absoluto y relativo . . . . . . . . . . . . . . . . . . . . . 18
1.10 Errores lineal y exponencial . . . . . . . . . . . . . . . . . . . 19
1.11 Condicionamiento de un problema . . . . . . . . . . . . . . . 21
i
ii ÍNDICE GENERAL
3 Métodos iterativos 76
3.1 Método de Gauss-Seidel . . . . . . . . . . . . . . . . . . . . . 76
3.2 Normas vectoriales . . . . . . . . . . . . . . . . . . . . . . . . 82
3.2.1 En Scilab . . . . . . . . . . . . . . . . . . . . . . . . . 83
3.3 Normas matriciales . . . . . . . . . . . . . . . . . . . . . . . . 83
3.3.1 En Scilab . . . . . . . . . . . . . . . . . . . . . . . . . 96
3.4 Condicionamiento de una matriz . . . . . . . . . . . . . . . . 96
3.5 Método de Jacobi . . . . . . . . . . . . . . . . . . . . . . . . . 101
3.6 Método iterativo general . . . . . . . . . . . . . . . . . . . . . 102
3.7 Método de sobrerrelajación . . . . . . . . . . . . . . . . . . . 103
ÍNDICE GENERAL iii
Preliminares
1.1 Notación
[a, b] = {x ∈ R : a ≤ x ≤ b},
]a, b[ = {x ∈ R : a < x < b}.
También es usual denotar el intervalo abierto por (a, b) pero puede con-
fundirse con la pareja ordenada (a, b).
C[a, b] es el conjunto de funciones continuas en el intervalo [a, b]; C n [a, b]
es el conjunto de funciones con n derivadas continuas sobre [a, b] (con esta
notación C[a, b] = C 0 [a, b]). Algunas veces, por brevedad, se dice que f es
de clase C n .
De manera análoga, C ∞ [a, b] es el conjunto de funciones que se pueden
derivar tantas veces como se desee. Ejemplos de estas funciones son: f (x) =
ex , g(x) = 3x5 − 8x2 + 12, h(x) = sen(2x).
Pn es el conjunto de todos los polinomios de grado menor o igual a n.
1
2 1. PRELIMINARES
f (x̃)
f (x̄)
a x̃ x̄ b
M
t
a c c̃ b
Teorema 1.2. Teorema del valor intermedio. Sea f continua en [a, b],
m = min{f (a), f (b)}, M = max{f (a), f (b)}. Si t es un valor intermedio,
m ≤ t ≤ M , entonces existe por lo menos un c ∈ [a, b] tal que
f (c) = t.
f (b) − f (a)
f ′ (c) =
b−a
.
4 1. PRELIMINARES
a c b
a c b
1.3 Sucesiones
Una sucesión es simplemente una función que va del conjunto de los números
naturales en los reales:
u :N → R
n 7→ un
Algunas veces las sucesiones están definidas para los naturales positivos (no
para el 0). Como se observa, habitualmente se escribe un en lugar de u(n).
Es frecuente denotar la sucesión {un }n∈N o {un }∞
n=0 o, si no hay confusión,
de manera aún más simple, {un } o un .
(−1)n m3 − 2m
Ejemplos de sucesiones: un = 1/n2 , vn = 5+ , w m = .
n3 100m2 + 20m
Una sucesión se puede definir de manera recurrente a partir del primer
término o de los primeros términos. Por ejemplo, la sucesión de números de
Fibonacci (Leonardo de Pisa) se define por:
u0 = 0
u1 = 1
un = un−2 + un−1 , para n ≥ 2.
Es usual escribir
lim xn = L
n→∞
xn −→ L
n→∞
6 1. PRELIMINARES
xn −→ L
Ejemplo 1.1. Sea
1
xn = 5 + .
n2
Veamos que el lı́mite es 5. Si ε = 0.01, se requiere que
5 + 1 − 5 ≤ 0.01
n2
1
≤ 0.01
n2
1
≤ n2
0.01
100 ≤ n2
10 ≤ n.
Es decir para ε = 0.01 basta conrtomar N ≥ 10. En general para un ε
1
cualquiera, basta con tomar N ≥ . 3
ε
lim xn = +∞
n→∞
o simplemente
xn −→ +∞
p(x̄) = f (x̄) y
p′ (x̄) = f ′ (x̄).
p(x̄) = f (x̄),
p′ (x̄) = f ′ (x̄),
p′′ (x̄) = f ′′ (x̄).
Entonces
f ′′ (x̄)
p(x) = f (x̄) + f ′ (x̄)(x − x̄) + (x − x̄)2 .
2
8 1. PRELIMINARES
p(x̄) = f (x̄),
p′ (x̄) = f ′ (x̄),
p′′ (x̄) = f ′′ (x̄),
..
.
p(n) (x̄) = f (n) (x̄).
Este polinomio es
f (n+1) (ξ(x))
Rn (x) = (x − x̄)n+1 (1.2)
(n + 1)!
∞
X f (k) (x̄)
f (x) = (x − x̄)k
k!
k=0
La anterior expresión es el desarrollo en serie de Taylor de f alrededor de
x̄.
|x − x̄|n+1
|Rn (x)| ≤ max f (n+1) (t)
(n + 1)! t∈I(x,x̄)
f ′ (x) = ex
f ′′ (x) = ex
f (n) (x) = ex
f (0) = 1
f ′ (0) = 1
f ′′ (x) = 1
f (n) (x) = 1
x2 x3 x4
ex = 1 + x + + + + ···
2 3! 4!
∞
X
x xn
e =
n!
n=0
f ′ (x) = cos(x)
f ′′ (x) = − sen(x)
f ′′′ (x) = − cos(x)
f (4) (x) = sen(x)
f (5) (x) = cos(x)
10 1. PRELIMINARES
f (0) = 0
f ′ (0) = 1
f ′′ (x) = 0
f ′′′ (x) = −1
f (4) (0) = 0
f (5) (0) = 1
x3 x5 x7
sen(x) = x − + − + ···
3! 5! 7!
Ejemplo 1.4. Obtener la serie de Taylor de f (x) = cos(x) alrededor de
x̄ = 0.
f ′ (x) = − sen(x)
f ′′ (x) = − cos(x)
f ′′′ (x) = sen(x)
f (4) (x) = cos(x)
f (5) (x) = − sen(x)
f (0) = 1
f ′ (0) = 0
f ′′ (x) = −1
f ′′′ (x) = 0
f (4) (0) = 1
f (5) (0) = 0
x2 x4 x6
cos(x) = 1 − + − + ···
2! 4! 6!
Ejemplo 1.5. Obtener el polinomio de Taylor de orden 2 de cos(x) alrede-
dor de π, acotar el error para x = 3 y calcular el error.
cos(π)
p2 (x) = cos(π) − sen(π)(x − π) − (x − π)2
2
1
p2 (x) = −1 + (x − π)2
2
1.5. NOTACIÓN O GRANDE 11
|3 − π|3
|error| ≤ max | sen(t)|
6 t∈[3,π]
n
X f (k) (x̄)
pn (x̄ + h) = hk (1.3)
k!
k=0
f (n+1) (ξ(h))
Rn (x̄ + h) = hn+1 , ξ(h) ∈ I(0, h), (1.4)
(n + 1)!
∞
X f (k) (x̄)
f (x̄ + h) = hk . (1.5)
k!
k=0
f (x) = O(g(x))
Ejemplo 1.6. Sea f (x) = 4x3 +5x6 . Recordemos que si 0 < y < 1, entonces
y > y 2 > y 3 > · · · . Entonces, si |x| < 1,
|x3 | ≤ |x|
|4x3 | ≤ 4|x|
|x6 | ≤ |x|
|5x6 | ≤ 5|x|
|4x3 + 5x6 | ≤ 9|x|
4x3 + 5x6 = O(x).
Sea {xk } una sucesión de números reales con lı́mite L. Se dice que la con-
vergencia tiene orden de convergencia p ≥ 1, si p es el mayor valor tal que
el siguiente lı́mite existe.
|xk+1 − L|
lim = β < ∞.
k→∞ |xk − L|p
|xk+1 − L|
lim = 0.
k→∞ |xk − L|
1
|xk+1 − L| |π + − π|
lim = lim k+1
k→∞ |xk − L| k→∞ 1
|π + − π|
k
1
= lim k + 1
k→∞ 1
k
k
= lim
k→∞ k + 1
= 1.
Luego la sucesión tiene orden de convergencia por lo menos igual a 1. Veamos
que pasa con p > 1. Se puede suponer que p = 1 + ε, con ε > 0.
1
|xk+1 − L|
lim = lim k + 1
k→∞ |xk − L|1+ε k→∞ 1
k 1+ε
k 1+ε
= lim
k→∞ 1 + k
k kε
= lim
k→∞ 1 + k
k
= lim lim k ε
k→∞ 1 + k k→∞
= (1)(+∞).
1
|xk+1 − L| k+1
lim = lim 2
k→∞ |xk − L|1+ε k→∞ 1
(2k )1+ε
(2k )1+ε
= lim
k→∞ 2k+1
2k 2kε
= lim k+1
k→∞ 2
2k kε
= lim k+1 lim 2
k→∞ 2 k→∞
1
= lim 2kε
2 k→∞
Ejemplo 1.9.
69
x1 =
10
xn = 6 + (xn−1 − 6)2 , n ≥ 2.
1 6.900000000000000
2 6.810000000000000
3 6.656100000000000
4 6.430467210000001
5 6.185302018885185
6 6.034336838202925
7 6.001179018457774
8 6.000001390084524
9 6.000000000001933
10 6.000000000000000
1.7. NÚMEROS EN UN COMPUTADOR 15
xn = 6 + yn , n = 1, 2, ...
9
y1 =
10
2
yn = yn−1 , n = 2, 3, ...
2n−1
9
yn = , n = 1, 2, ...
10
Como yn → 0, entonces xn → 6.
|xk+1 − L| yk+1
lim = lim p
k→∞ |xk − L|p k→∞ yk
yk2
= lim p
k→∞ yk
= lim yk2−p
k→∞
|xk+1 − L| 2−(2+ε)
lim p
= lim yk
k→∞ |xk − L| k→∞
1
= lim ε
k→∞ yk
= +∞
format(30)
x = 1/3
El resultado es
0.3333333333333333148296
Únicamente hay 16 dı́gitos correctos, los demás son “basura” producida por
Scilab para satisfacer el formato deseado. Esto nos indica que en Scilab, en la
representación interna de un número, no hay más de 16 cifras significativas.
Relacionado con el concepto anterior, está el épsilon de la máquina, que se
define ası́:
x = 1;
while x/10 > 0.0
x0 = x;
x = x/10;
end
x_final = x0
truncamiento(1234.56789, 2) = 1200
truncamiento(1234.56789, 6) = 1234.56
redondeo(1234.56789, 2) = 1200
redondeo(1234.56789, 6) = 1234.57
x y x̃ ỹ z z̃ ea er
1/7 2/3 0.14286 0.66667 −11/21 −0.52381 4.8e-7 9.1e-7
1/7 0.14284 0.14286 0.14284 0.00001714... 0.00002 2.9e-6 1.7e-1
Los principales casos en los que los errores pueden ser grandes, son:
Estos casos, en cuanto sea posible, deben ser evitados y, si no es posible, los
resultados deben ser interpretados de manera muy cuidadosa.
x = [1 1]T ,
x′ = [−1998 2002]T ,
Uno de los problemas numéricos más frecuentes, o tal vez el más frecuente,
consiste en resolver un sistema de ecuaciones de la forma
Ax = b (2.1)
2.1 En Scilab
x = a\b
Por ejemplo
22
2.1. EN SCILAB 23
da como resultado
x =
2.
- 3.
n = 500;
a = rand(n,n);
x = rand(n,1);
b = a*x;
tic()
x1 = a\b;
t_sol = toc();
tic()
x2 = inv(a)*b;
t_inv = toc();
error1 = norm(x1-x);
error2 = norm(x2-x);
2.2 Notación
x = A−1 b.
Para resolver un sistema de ecuaciones por la regla de Cramer, hay que calcu-
lar n + 1 determinantes, luego el número total de multiplicaciones necesarias
para resolver un sistema de ecuaciones por la regla de Cramer, calculando
los determinantes por cofactores, es superior a (n + 1)!.
21! = 5.1091E19.
Siendo muy optimistas (sin tener en cuenta las sumas y otras operaciones
concomitantes), supongamos que un computador del año 2000 hace 1000
millones de multiplicaciones por segundo. Entonces, el tiempo necesario
para resolver un sistema de ecuaciones de orden 20 por la regla de Cramer
y el método de cofactores es francamente inmanejable:
El caso más sencillo de (2.1) corresponde a una matriz diagonal. Para ma-
trices triangulares, en particular para las diagonales, el determinante es el
26 2. SOLUCIÓN DE SISTEMAS LINEALES
xn = bn /ann
para i = n − 1, ..., 1
xi = (bi − A(i, i + 1 : n) x(i + 1 : n))/aii
fin-para
x(n) = b(n)/a(n,n)
for i=n-1:-1:1
x(i) = ( b(i) - a(i,i+1:n)*x(i+1:n) )/a(i,i)
end
res = 0
if min(abs(diag(a))) <= eps, return, end
res = 1
n = size(a,1)
x = zeros(n,1)
x(n) = b(n)/a(n,n)
for k = n-1:-1:1
x(k) = (b(k) - a(k,k+1:n)*x(k+1:n) )/a(k,k)
end
endfunction
for k = n:-1:1
x(k) = (b(k) - a(k,k+1:n)*x(k+1:n) )/a(k,k)
end
aunque hay más operaciones fuera de las anteriores, por ejemplo las com-
paraciones y las operaciones entre enteros. Las cuatro operaciones se cono-
cen con el nombre genérico de operaciones de punto flotante flops (floating
point operations). Algunas veces se hacen dos grupos: por un lado sumas y
restas, y por otro multiplicaciones y divisiones. Si se supone que el tiempo
necesario para efectuar una multiplicación es bastante mayor que el tiempo
de una suma, entonces se acostumbra a dar el número de multiplicaciones
(o divisiones). El diseño de los procesadores actuales muestra tendencia al
hecho de que los dos tiempos sean comparables. Entonces se acostumbra a
evaluar el número de flops.
i−1
X
bi − aij xj
j=1
xi = · (2.3)
aii
para i = 1, ..., n
xi = (bi − A(i, 1 : i − 1) x(1 : i − 1))/aii
fin-para
para k = 1, ..., n − 1
buscar ceros en la columna k, por debajo de la diagonal.
fin-para k
para k = 1, ..., n − 1
para i = k + 1, ..., n
buscar ceros en la posición de aik .
fin-para i
fin-para k
Ejemplo 2.3. Consideremos el siguiente sistema de ecuaciones:
4x1 + 3x2 − 2x3 + x4 = 4
3x1 + 2x2 + x3 + 5x4 = −8
−2x1 + 3x2 + x3 + 2x4 = −7
−5x1 + x3 + x4 = −8
2.7. MÉTODO DE GAUSS 31
Inicialmente hay que buscar ceros en la primera columna. Para buscar cero
en la posición (2,1), fila 2 y columna 1, se hace la siguiente operación:
4 3 -2 1 4
0 -0.25 2.5 4.25 -11
0 0 45 79 -203
0 3.75 -1.5 2.25 -3
Para buscar cero en la posición (4,2) se hace siguiente operación:
para k = 1, ..., n − 1
para i = k + 1, ..., n
lik = aik /akk , aik = 0
A(i, k + 1 : n) = A(i, k + 1 : n)−lik∗A(k, k + 1 : n)
bi = bi −lik∗bk
fin-para i
fin-para k
para k = 1, ..., n − 1
para i = k + 1, ..., n
si |akk | ≤ ε ent
buscar m, k + 1 ≤ m ≤ n, tal que |amk | > ε
si no fue posible ent salir
intercambiar(A(k, k : n), A(m, k : n))
intercambiar(bk , bm )
fin-si
lik = aik /akk , aik = 0
A(i, k + 1 : n) = A(i, k + 1 : n)−lik∗A(k, k + 1 : n)
bi = bi −lik∗bk
fin-para i
fin-para k
si |ann | ≤ ε ent salir
para k = 1, ..., n
para i = k + 1, ..., n
si |akk | ≤ ε ent
buscar m, k + 1 ≤ m ≤ n, tal que |amk | > ε
si no fue posible ent salir
intercambiar(A(k, k : n), A(m, k : n))
intercambiar(bk , bm )
fin-si
lik = aik /akk , aik = 0
A(i, k + 1 : n) = A(i, k + 1 : n)−lik∗A(k, k + 1 : n)
bi = bi −lik∗bk
fin-para i
fin-para k
//
// indic valdra 1 si todo funciono bien,
// 0 si la matriz es singular o casi.
//
n = size(a,1)
if argn(2) < 3, eps = 1.0e-10, end
for k=1:n
if abs(a(k,k)) <= eps
m = posNoNulo(a, k)
if m == 0
indic = 0
return
end
t = a(k,k:n)
a(k,k:n) = a(m,k:n)
a(m,k:n) = t
t = b(k)
b(k) = b(m)
b(m) = t
end
for i=k+1:n
lik = a(i,k)/a(k,k)
a(i,k) = 0
a(i,k+1:n) = a(i,k+1:n) - lik*a(k,k+1:n)
b(i) = b(i) - lik*b(k)
end
end
indic = 1
endfunction
//----------------------------------------------------------
function m = posNoNulo(a, k, eps)
// Busca la posicion del primer elemento no nulo en la
// columna k, debajo de la diagonal.
//
// Si no es posible encontrarlo, m valdra 0.
//
if argn(2) < 3, eps = 1.0e-10, end
n = size(a,1)
for i = k+1:n
if abs(a(i,k)) >= eps
36 2. SOLUCIÓN DE SISTEMAS LINEALES
m = i
return
end
end
m = 0
endfunction
//----------------------------------------------------------
function [x, indic] = Gauss(a, b, eps)
// Solucion de un sistema de ecuaciones
// por el metodode Gauss.
//
// indic valdra 1 si todo funciono bien,
// en este caso el vector columna x
// sera la solucion.
// 0 si la matriz es singular o casi
// -1 los tamanos son incompatibles.
//
indic = -1
x = []
n = verifTamanoAb(a, b)
if n == 0, return, end
indic = 0
x = []
[a, b, res] = triangulariza(a, b, eps)
if res == 0, return, end
indic = 1
x = solTriSup(a, b, eps)
endfunction
//----------------------------------------------------------
function n = verifTamanoAb(a, b)
// Esta funcion verifica si los tamanos de a, b
// corresponden a un sistema cuadrado a x = b.
// Devuelve n (num. de filas) si todo esta bien,
// devuelve 0 si hay errores.
Triangularización
Consideremos inicialmente la búsqueda de cero en la posición (2, 1). Para
efectuar A(2, 2 : n) = A(2, 2 : n) − lik ∗ A(1, 2 : n) es necesario hacer n − 1
sumas y restas. Para b2 = b2 −lik∗b1 es necesario una resta. En resumen n
sumas (o restas). Multiplicaciones y divisiones: una división para calcular
lik; n − 1 multiplicaciones para lik ∗ A(1, 2 : n) y una para lik∗b1 . En
resumen, n + 1 multiplicaciones (o divisiones).
Para obtener un cero en la posición (3, 1) se necesita exactamente el mismo
número de operaciones. Entonces para la obtener ceros en la primera columna:
n−1
X n−1
X n3 n n3
i(i + 1) = (i2 + i) = − ≈ ·
3 3 3
i=1 i=1
n−1
X n−1
X n3 n2 5n n3
i(i + 2) = (i2 + 2i) = + − ≈ ·
3 2 6 3
i=1 i=1
Número de operaciones:
n3 n n3 n2 5n 2n3 n2 7n 2n3
− + + − = + − ≈ ·
3 3 3 2 6 3 2 6 3
Proceso completo
El número de operaciones para las dos partes, triangularización y solución
del sistema triangular, es
2.8 Factorización LU
Ax = b̃,
LU x = b̃,
Ly = b̃,
donde U x = y.
En resumen:
40 2. SOLUCIÓN DE SISTEMAS LINEALES
Al resolver
1 0 0 0 y1 8
0.75 1 0 0 y2 30
=
-0.5 -18 1 0 y3 15
-1.25 -15 0.8 1 y4 2
T
se obtiene y = 8 24 451 11.2 . Al resolver
4 3 -2 1 x1 8.0
0 -0.25 2.5 4.25 x2 24.0
=
0 0 45 79 x3 451.0
0 0 0 2.8 x4 11.2
T
se obtiene la solución final x = 1 2 3 4 . 3
..
.
lik = aik /akk
aik = lik
A(i, k + 1 : n − 1) = A(i, k + 1 : n − 1)−lik∗A(k, k + 1 : n − 1)
bi = bi −lik∗bk
..
.
Buscar ceros en las posiciones de a32 , a42 se hace usando los valores de
lik= 2/3 = 0.6666 y 1. Se obtiene
-5 0 1 1 -8
0 3 0.6 1.6 -3.8
0 0 1.2 4.5333 -10.2667
0 0 -1.8 0.2 1.4
kx1 − x∗ k , kx2 − x∗ k.
Si kx1 −x∗ k < kx2 −x∗ k, entonces x1 es, entre x1 y x2 , la mejor aproximación
de x∗ . Cuando no se conoce x∗ , entonces se utiliza la norma del vector
2.9. MÉTODO DE GAUSS CON PIVOTEO PARCIAL 45
• Resolver P T z = b, o sea, z = P b.
• Resolver Ly = z.
• Resolver U x = y.
2.11. MÉTODO DE CHOLESKY 49
-8
-7
z = Pb =
4
-8
-8
-3.8
Ly = z , entonces y =
1.4
-9.3333
1
0
U x = y , entonces x =
-1 3
-2
[L, U, P] = lu(A)
xT Ax > 0, ∀ x ∈ Rn , x 6= 0. (2.5)
50 2. SOLUCIÓN DE SISTEMAS LINEALES
Si A es simétrica,
n
X n−1
X n
X
xT Ax = aii x2i + 2 aij xi xj .
i=1 i=1 j=i+1
• A es definida positiva.
• λi > 0, ∀i.
• δi > 0, ∀i.
• Existe U matriz triangular superior e invertible tal que A = U T U .
a = [ 4 -6; -6 25]
u = chol(a)
Sea
1 2
A= .
2 5
Entonces
u11 0 u11 u12 1 2
=
u12 u22 0 u22 2 5
u211 = 1
u11 u12 = 2,
u212 + u222 = 5
2.11. MÉTODO DE CHOLESKY 53
Se deduce que
u11 = 1
u12 = 2,
u22 = 1,
1 2
U = .
0 1
Sea
1 2
A= .
2 4
Entonces
u11 0 u11 u12 1 2
=
u12 u22 0 u22 2 4
u211 = 1
u11 u12 = 2,
u212 + u222 = 4
Se deduce que
u11 = 1
u12 = 2,
u22 = 0,
1 2
U = .
0 0
Sea
1 2
A= .
2 3
54 2. SOLUCIÓN DE SISTEMAS LINEALES
Entonces
u11 0 u11 u12 1 2
=
u12 u22 0 u22 2 3
u211 = 1
u11 u12 = 2,
u212 + u222 = 3
Se deduce que
u11 = 1
u12 = 2,
u222 = −1.
En el caso general,
u11 u11 · · · u1k · · · u1j · · · u1n
.. ..
. .
u1k · · · ukk ukk · · · ukj · · · ukn
.. ..
. .
u1j · · · ukj · · · ujj ujj · · · ujn
.. ..
. .
u1n · · · ukn · · · ujn · · · unn unn
u211 = a11 .
Luego
√
u11 = a11 . (2.6)
El producto de la fila 1 de U T por la columna j de U da:
Luego
a1j
u1j = , j = 2, ..., n. (2.7)
u11
Al hacer el producto de la fila 2 de U T por la columna 2 de U , se puede
calcular u22 . Al hacer el producto de la fila 2 de U T por la columna j de
2.11. MÉTODO DE CHOLESKY 55
k
X
u2ik = akk
i=1
k−1
X
u2ik + u2kk = akk .
i=1
Luego
v
u k−1
u X
ukk t
= akk − u2ik , k = 2, ..., n. (2.8)
i=1
k
X
uik uij = akj .
i=1
Luego
k−1
X
akj − uik uij
i=1
ukj = , k = 2, ..., n, j = k + 1, ..., n. (2.9)
ukk
√
u11 = 16 = 4
56 2. SOLUCIÓN DE SISTEMAS LINEALES
−12
u12 = = −3
4
8
u13 = =2
4
−16
u14 = = −4
4
p
u22 = 18 − (−3)2 = 3
−6 − (−3)(2)
u23 = =0
3
9 − (−3)(−4)
u24 = = −1
3
p
u33 = 5 − (22 + 02 ) = 1
Entonces,
4 −3 2 −4
0 3 0 −1
U =
0
. 3
0 1 −2
0 0 0 5
Para ahorrar espacio de memoria, los valores ukk y ukj se pueden almacenar
sobre los antiguos valores de akk y akj . O sea, al empezar el algoritmo se
tiene la matriz A. Al finalizar, en la parte triangular superior del espacio
2.11. MÉTODO DE CHOLESKY 57
datos: A, ε
para k = 1, ..., n
cálculo de t según (2.10)
si t ≤ √ ε ent salir
akk = t
para j = k + 1, ..., n
cálculo de akj según (2.12)
fin-para j
fin-para k
n = size(A,1)
U = zeros(n,n)
58 2. SOLUCIÓN DE SISTEMAS LINEALES
for k = 1:n
t = A(k,k) - U(1:k-1,k)’*U(1:k-1,k)
if t <= eps
printf(’Matriz no definida positiva.\n’)
ind = 0
return
end
U(k,k)= sqrt(t)
for j = k+1:n
U(k,j) = ( A(k,j) - U(1:k-1,k)’*U(1:k-1,j) )/U(k,k)
end
end
ind = 1
endfunction
n3 n2 n n3
+ + ≈ .
3 2 6 3
resolver U T y = b, (2.13)
resolver U x = y. (2.14)
Resolver cada uno de los dos sistemas es muy fácil. El primero es triangular
inferior, el segundo triangular superior. El número total de operaciones
para resolver el sistema está dado por la factorización más la solución de
dos sistemas triangulares.
n3 n3
Número de operaciones ≈ + 2 n2 ≈ ·
3 3
60 2. SOLUCIÓN DE SISTEMAS LINEALES
Esto quiere decir que para valores grandes de n, resolver un sistema, con
A definida positiva, por el método de Cholesky, gasta la mitad del tiempo
requerido por el método de Gauss.
El método de Cholesky se utiliza para matrices definidas positivas. Pero
no es necesario tratar de averiguar por otro criterio si la matriz es definida
positiva. Simplemente se trata de obtener la factorización de Cholesky de
A simétrica. Si fue posible, entonces A es definida positiva y se continúa
con la solución de los dos sistemas triangulares. Si no fue posible obtener la
factorización de Cholesky, entonces A no es definida positiva y no se puede
aplicar el método de Cholesky para resolver Ax = b.
Al resolver U T y = b se obtiene
//
// info valdra 1 si a es definida positiva,
// asi x sera un vector columna
// con la solucion,
// 0 si a no es definida positiva.
//
[a, info] = Cholesky(a)
if info == 0, return, end
y = sol_UT_y_b(a, b)
x = solTriSup(a, y)
endfunction
2.12.1 En Scilab
La orden para hallar por la solución por mı́nimos cuadrados es la misma que
para resolver sistemas de ecuaciones cuadrados, a saber, a\b . Por ejemplo,
para resolver el sistema
2 3 43
x
4 5 1 “=” 77
x2
6 7 109
basta con
a = [ 2 3; 4 5; 7 6 ], b = [ 43 77 109 ]’
x = a\b
El resultado obtenido es
x =
7.6019417
9.3009709
∂ϕ
(x̄) ·
∂xi
∂ϕ
= 9(4x31 + 6x4 )8 (12x21 ) + 5x2 ,
∂x1
∂ϕ
= 5x1 ,
∂x2
∂ϕ
= 0,
∂x3
∂ϕ
= 54(4x31 + 6x4 )8 + 8.
∂x4
∂f
=2(a11 x1 + a12 x2 + a13 x3 − b1 )a11
∂x1
+ 2(a21 x1 + a22 x2 + a23 x3 − b2 )a21
+ 2(a31 x1 + a32 x2 + a33 x3 − b3 )a31
+ 2(a41 x1 + a42 x2 + a43 x3 − b4 )a41 .
∂f
=2(A1·x − b1 )a11 + 2(A2·x − b2 )a21 + 2(A3·x − b3 )a31
∂x1
+ 2(A4·x − b4 )a41 .
64 2. SOLUCIÓN DE SISTEMAS LINEALES
De manera semejante
∂f
= 2 AT (Ax − b) 2 ,
∂x2
∂f
= 2 AT (Ax − b) 3
∂x3
Igualando a cero las tres derivadas parciales y quitando el 2 se tiene
AT (Ax − b) 1 = 0,
AT (Ax − b) 2 = 0,
AT (Ax − b) 3 = 0
Es decir,
AT (Ax − b) = 0,
AT A x = AT b . (2.16)
Por ser AT A invertible, hay una única solución de (2.16), o sea, hay un
solo vector x que hace que las derivadas parciales sean nulas. En general,
las derivadas parciales nulas son simplemente una condición necesaria para
obtener el mı́nimo de una función (también lo es para máximos o para puntos
de silla), pero en este caso, como AT A es definida positiva, f es convexa, y
entonces anular las derivadas parciales se convierte en condición necesaria y
suficiente para el mı́nimo.
En resumen, si las columnas de A son linealmente independientes, entonces
la solución por mı́nimos cuadrados existe y es única. Para obtener la solución
por mı́nimos cuadrados se resuelven las ecuaciones normales.
Como AT A es definida positiva, (2.16) se puede resolver por el método de
Cholesky. Si m ≥ n y al hacer la factorización de Cholesky resulta que
AT A no es definida positiva, entonces las columnas de A son linealmente
dependientes.
Si el sistema Ax = b tiene solución exacta, ésta coincide con la solución por
mı́nimos cuadrados.
El error, Ax − b, es:
−0.0628
0.0196
−0.0039 . 3
0.0275
66 2. SOLUCIÓN DE SISTEMAS LINEALES
x = (2, −1) .
El error, Ax − b, es:
0
0
.
0
0
En este caso, el sistema inicial tenı́a solución exacta y la solución por
mı́nimos cuadrados coincide con ella. 3
2.13. SISTEMAS TRIDIAGONALES 67
F1 C1 : d1 = a11
F1 C2 : u1 = a12
F2 C1 : f1 d1 = a21
F2 C2 : f1 u1 + d2 = a22
F2 C3 : u2 = a23
F3 C2 : f2 d2 = a32
F3 C3 : f2 u2 + d3 = a33
F3 C4 : u3 = a34
..
.
Fi Ci−1 : fi−1 di−1 = ai,i−1
Fi Ci : fi−1 ui−1 + di = aii
Fi Ci+1 : ui = ai,i+1
2.13. SISTEMAS TRIDIAGONALES 69
d1 = a11 ,
y 1 = b1 ,
fi−1 yi−1 + yi = bi ,
dn xn = yn ,
di xi + ui xi+1 = yi .
y 1 = b1 ,
yi = bi − fi−1 yi−1 , i = 2, ..., n,
yn (2.18)
xn = ,
dn
yi − ui xi+1
xi = , i = n − 1, n − 2, ..., 1.
di
2 4 0 0 −8
3 5 6 0 1
A=
0 −4 −5
, b=
−2 .
1
0 0 −1 −2 −10
70 2. SOLUCIÓN DE SISTEMAS LINEALES
Entonces
d1 = 2 ,
u1 = 4 ,
3
f1 = = 1.5 ,
2
d2 = 5 − 1.5 × 4 = −1 ,
u2 = 6 ,
−4
f2 = = 4,
−1
d3 = −5 − 4 × 6 = −29 ,
u3 = 1,
−1
f3 = = 0.034483 ,
−29
d4 = −2 − 0.034483 × 1 = −2.034483 ,
Las fórmulas (2.17) y (2.18) se pueden utilizar sin ningún problema si todos
los di son no nulos. Algún elemento diagonal de U resulta nulo si la matriz
A no es invertible o si simplemente A no tiene factorización LU .
Ejemplo 2.20. Consideremos las dos matrices siguientes:
2 −3 ′ 0 2
A= , A = .
−8 12 3 4
2.13. SISTEMAS TRIDIAGONALES 71
datos: d, f, u, b, ε
datos: A, ε
resultados: la inversa almacenada en A, indic
Factorización:
p = (1, 2, ..., n − 1)
para k = 1 : n − 1
determinar m tal que |amk | = max{ |aik | : i = k, ..., n}
si |amk | ≤ ε
indic = 0, parar
fin-si
pk = m
si m > k
A(k, : ) ↔ A(m, : )
fin-si
A(k + 1 : n, k) = A(k + 1 : n, k)/akk
A(k + 1 : n, k + 1 : n) = A(k + 1 : n, k + 1 : n) − A(k + 1 : n, k)A(k, k + 1 : n)
fin-para
si |ann | ≤ ε
indic = 0, parar
fin-si
indic = 1
2.14. CÁLCULO DE LA INVERSA 73
Cálculo de U −1 :
para k = 1 : n
akk = 1/akk
para i = 1 : k − 1
aik = −akk A(i, i : k − 1)A(i : k − 1, k)
fin-para
fin-para
Cálculo de U −1 L−1 :
para k = n − 1 : −1 : 1
t = A(k + 1 : n, k)
A(k + 1 : n, k) = 0
A( : , k) = A( : , k) − A( : , k + 1 : n) t
fin-para
Reordenamiento de columnas :
para k = n − 1 : −1 : 1
si pk 6= k
A( : , k) ↔ A( : , pk )
fin-si
fin-para
Factorisacion
k = 1
m = 2
p :
2 2 3
intercambio de filas : 1 2
A despues de intercambio
-5.0000 1.0000 2.0000 1.0000
-2.0000 -4.0000 4.0000 -2.0000
74 2. SOLUCIÓN DE SISTEMAS LINEALES
k = 2
m = 2
p :
2 2 3
A despues de operaciones
-5.0000 1.0000 2.0000 1.0000
0.4000 -4.4000 3.2000 -2.4000
-0.8000 0.5000 0.0000 -2.0000
0.4000 0.7727 -2.2727 0.4545
k = 3
m = 4
p :
2 2 4
intercambio de filas : 3 4
A despues de intercambio
-5.0000 1.0000 2.0000 1.0000
0.4000 -4.4000 3.2000 -2.4000
0.4000 0.7727 -2.2727 0.4545
-0.8000 0.5000 0.0000 -2.0000
A despues de operaciones
-5.0000 1.0000 2.0000 1.0000
0.4000 -4.4000 3.2000 -2.4000
0.4000 0.7727 -2.2727 0.4545
-0.8000 0.5000 -0.0000 -2.0000
Expresiones explicitas de L, U, P
L
1.0000 0.0000 0.0000 0.0000
0.4000 1.0000 0.0000 0.0000
0.4000 0.7727 1.0000 0.0000
-0.8000 0.5000 -0.0000 1.0000
U
-5.0000 1.0000 2.0000 1.0000
0.0000 -4.4000 3.2000 -2.4000
0.0000 0.0000 -2.2727 0.4545
0.0000 0.0000 0.0000 -2.0000
P :
0 1 0 0
1 0 0 0
0 0 0 1
0 0 1 0
3
Métodos iterativos
lim xk = x∗ .
k→∞
Generalmente los métodos indirectos son una buena opción cuando la matriz
es muy grande y dispersa o rala (sparse), es decir, cuando el número de
elementos no nulos es pequeño comparado con n2 , número total de elementos
de A. En estos casos se debe utilizar una estructura de datos adecuada que
permita almacenar únicamente los elementos no nulos.
76
3.1. MÉTODO DE GAUSS-SEIDEL 77
n
X
|aii | > |aij | , ∀i.
j=1,j6=i
n
X
bi − aij xj
j=1,j6=i
xi ← ,
aii
n
X
bi − aij xj + aii xi
j=1
xi ← ,
aii
bi − Ai· x
xi ← xi + .
aii
Sean
ri = bi − Ai· x,
ri
δi = ·
aii
ε
kδk ≤ ·
max |aii |
datos: A, b, x0 , ε, maxit
x = x0
para k = 1, ...,maxit
nrmD← 0
para i = 1, ..., n
δi = (bi − Ai· x)/aii
xi ← xi + δi
nrmd←nrmD+|δi |
fin-para i
si nrmD ≤ ε ent x∗ ≈ x, salir
fin-para k
A continuación hay una versión, no muy eficiente, que permite mostrar los
resultados intermedios
for k = 1:maxit
//printf(’\n k = %d\n’, k)
D = 0
for i = 1:n
di = ( b(i) - A(i,:)*x )/A(i,i)
x(i) = x(i) + di
D = max(D, abs(di))
end
disp(x’)
if D < eps, return, end
end
ind = 0
endfunction
p √
22 + 32 + (−4)2 = 29
µ:V →R
µ(x) ≥ 0, ∀x ∈ V,
µ(x) = 0 sssi x = 0,
µ(αx) = |α| µ(x), ∀α ∈ R, ∀x ∈ V,
µ(x + y) ≤ µ(x) + µ(y), ∀x, y ∈ V. (desigualdad triangular)
n
!1/2
X
||x||2 = ||x|| = x2i norma euclidiana,
i=1
n
!1/p
X
p
||x||p = |xi | norma de Holder de orden p ≥ 1,
i=1
||x||∞ = ||x||max = max |xi |,
1≤i≤n
α||x|| con α > 0 y || || una norma,
√
||x||A = xT Ax con A definida positiva.
||x||1 = 7,
||x||2 = 5,
||x||∞ = 4.
3.2.1 En Scilab
está definido el producto, es interesante contar con normas que tengan carac-
terı́sticas especiales relativas al producto entre matrices y al producto entre
una matriz y un vector. En particular en algunos casos es conveniente que
se tengan estas dos propiedades:
entonces
19 22 17
AB = , Ax = ,
43 50 39
pero
Sean || ||m una norma matricial sobre Rn×n y || ||v una norma sobre Rn .
Estas dos normas se llaman compatibles o consistentes si, para toda matriz
A ∈ Rn×n y para todo x ∈ Rn
||Ax||
|||A||| = sup (3.1)
x6=0 ||x||
||Ax||
|||A||| = max (3.2)
x6=0 ||x||
Proposición 3.1. La definición anterior está bien hecha, es decir, ||| ||||
es una norma, es matricial y es compatible con || ||.
Demostración. Sea
||Ax||
µ(A) = sup
x6=0 ||x||
Ante todo es necesario mostrar que la función µ está bien definida, o sea,
para toda matriz A,
||Ax||
µ(A) = sup < ∞.
x6=0 ||x||
Veamos
86 3. MÉTODOS ITERATIVOS
x
A ||x||
||x||
µ(A) = sup
x6=0 ||x||
||x||A x
||x||
= sup
x6=0 ||x||
x
||x|| A
||x||
= sup
x6=0 ||x||
x
= sup A
x6=0 ||x||
= sup ||Aξ||
||ξ||=1
ej
µ(A) ≥ ||Av|| = A j
||e ||
A·j
= j
||e ||
||A·j ||
= > 0.
||ej ||
µ(λA) = max ||λAx|| = max |λ| ||Ax|| = |λ| max ||Ax|| = |λ|µ(A).
||x||=1 ||x||=1 ||x||=1
µ(A + B) = sup ||(A + B)x|| = sup ||Ax + Bx|| ≤ sup (||Ax|| + ||Bx||)
||x||=1 ||x||=1 ||x||=1
Para las 3 normas vectoriales más usadas, las normas matriciales generadas
son:
n
X
||A||1 = max |aij |, (3.5)
1≤j≤n
i=1
p
||A||2 = ρ(AT A) (norma espectral), (3.6)
Xn
||A||∞ = max |aij |. (3.7)
1≤i≤n
j=1
88 3. MÉTODOS ITERATIVOS
1/2
X
||A||F = (aij )2 . (3.8)
i,j
√
Para cualquier norma generada ||I|| = 1. Como ||I||F = n, entonces esta
norma no puede ser generada por ninguna norma vectorial
1 2
A=
3 −4
Entonces
T 10 −10
A A=
−10 20
||A||1 = 6,
||A||2 = 5.1166727,
||A||∞ = 7.
n
X
Proposición 3.2. . ||A||1 = max |aij |
j
i=1
3.3. NORMAS MATRICIALES 89
Demostración.
Pn
donde sj = i=1 |aij |. Si αj , βj ≥ 0 para todo j, entonces
n
X n
X
αj βj ≤ max βj αj .
j
j=1 j=1
Luego
90 3. MÉTODOS ITERATIVOS
Xn
||A||1 ≤ max max sj |xj |
||x||1 =1 j
j=1
= max max sj
||x||1 =1 j
= max sj
j
n
X
= max |aij |
j
i=1
En resumen
n
X
||A||1 ≤ max |aij |.
j
i=1
n
X n
X
|aik | = max |aij |
j
i=1 i=1
es decir,
n
X
||A||1 ≥ max |aij |.
j
i=1
p
Proposición 3.3. ||A||2 = ρ(AT A).
3.3. NORMAS MATRICIALES 91
Demostración.
||A||22 = max xT AT Ax
||x||2 =1
λ1 ≥ λ2 ≥ ... ≥ λn ≥ 0.
(AT A)v i = λi v i ,
T
v i v j = δij .
n
X
x= αi v i .
i=1
Entonces
92 3. MÉTODOS ITERATIVOS
n
X
AT Ax = AT A αi v i
i=1
n
X
= αi AT Av i
i=1
n
X
= αi λi v i
i=1
T !
n
X n
X
xT AT Ax = αj v j αi λi v i
j=1 i=1
n
X
= αi2 λi
i=1
n
X
≤ λ1 αi2
i=1
= λ1
En resumen,
||A||22 ≤ λ1
p
||A||2 ≤ λ1
√
||A||2 ≥ xT AT Ax pra todo x con ||x||2 = 1
p
||A||2 ≥ v 1 T AT Av 1
p
||A||2 ≥ v 1 T λ1 v 1
p
||A||2 ≥ λ1 v 1 T v 1
p
||A||2 ≥ λ1 .
n
X
Proposición 3.4. ||A||∞ = max |aij |
i
j=1
3.3. NORMAS MATRICIALES 93
Demostración.
n
X
||A||∞ ≤ max max |aij | ||x||∞
||x||∞ =1 i
j=1
n
X
= max ||x||∞ max |aij |
||x||∞ =1 i
j=1
n
X
= max |aij |
i
j=1
n
X n
X
|akj | = max |aij |
i
j=1 j=1
0 si akj = 0
x̄j = |akj |
signo(akj ) = si akj 6= 0.
akj
94 3. MÉTODOS ITERATIVOS
||x||∞ = 1,
||A||∞ ≥ ||Ax||∞ si ||x||∞ = 1,
||A||∞ ≥ ||Ax̄||∞
= max |Ai·x̄|
i
= |Ai·x̄| para todo i,
||A||∞ ≥ |Ak·x̄|
X
n |akj |
= akj
j=1 akj
X
n
= |akj |
j=1
n
X
= |akj |
j=1
n
X
= max |aij |.
i
j=1
RESUMEN DE RESULTADOS
• ||I|| = 1.
3.3. NORMAS MATRICIALES 95
n
X
• ||A||1 = max |aij |
j
i=1
n
X
• ||A||∞ = max |aij |
i
j=1
p
• ||A||2 = ρ(AT A)
• ||A||2 = max{σ1 , σ2 , ..., σn } = max{valores singulares de A} (ver [AlK02]).
• ||A||2 = ρ(A) si A 0.
v
u n
p uX
• ||A||F = tr(AT A) = t σ2 i
i=1
3.3.1 En Scilab
b′ = b + ∆b,
x̄′ = x̄ + ∆x.
∆x = x̄′ − x̄
= A−1 b′ − A−1 b
= A−1 (b + ∆b) − A−1 b
= A−1 ∆b.
b = Ax
||b|| ≤ ||A|| ||x̄||
1 ||A||
≤
||x̄|| ||b||
||∆x|| ||∆b||
≤ ||A|| ||A−1 || . (3.9)
||x̄|| ||b||
Entonces
||∆x|| ||∆b||
≤ κ(A) . (3.10)
||x̄|| ||b||
−10 −7
A= .
6 4
98 3. MÉTODOS ITERATIVOS
Entonces
−1 2 7/2
A =
−3 −5
136 94
AT A =
94 65
−1 T −1 13 22
A A =
22 149/4
esp(AT A) = {0.0199025, 200.9801}
T
esp(A−1 A−1 ) = {0.0049756, 50.245024}
||A||2 = 14.176745
||A−1 ||2 = 7.0883725
κ2 (A) = 100.49005
||A||1 = 16
||A−1 ||1 = 17/2
κ1 (A) = 136
||A||∞ = 17
||A−1 ||∞ = 8
κ∞ (A) = 136. 3
• κ(A) ≥ 1.
• κ(αA) = κ(A) si α 6= 0.
• κ2 (A) = 1 si y solamente si A es un múltiplo de una matriz ortogonal
(o unitaria).
Entonces
∆b = [0.01 − 0.01]T ,
||∆b||
= 0.0005,
||b||
κ(A) = 1.0752269.
x = [1.9999474 0.0010526]T ,
x′ = [1.9998947 0.0021053]T ,
∆x = [−0.0000526 .0010526]T ,
||∆x||
= 0.0005270,
||x||
||∆b||
κ(A) = 0.0005376.
||b||
Entonces
∆b = [0.01 − 0.01]T ,
||∆b||
= 0.0005,
||b||
−1 −99900 100000
A = ,
100000 −100100
κ(A) = 4000002.
100 3. MÉTODOS ITERATIVOS
x = [1 1]T ,
x′ = [−1998 2002]T ,
∆x = [−1999 2001]T ,
||∆x||
= 2000.0002,
||x||
||∆b||
κ(A) = 2000.0008.
||b||
Entonces
∆b = [0.01 0.01]T ,
||∆b||
= 0.0005,
||b||
−99900 100000
A−1 = ,
100000 −100100
κ(A) = 4000002.
x = [1 1]T ,
x′′ = [2 0]T ,
∆x = [1 − 1]T ,
||∆x||
= 1,
||x||
||∆b||
κ(A) = 2000.0008.
||b||
3.5. MÉTODO DE JACOBI 101
Ejemplo 3.9.
4 1 −1 7 1.2
A= 2 5 0 , b = 19 , x0 = 1.5 .
−2 3 10 45 1.6
Gauss-Seidel Jacobi
x1 x2 x3 x1 x2 x3
1.2 1.5 1.6 1.2 1.5 1.6
1.775 1.5 1.6 1.775 1.5 1.6
1.775 3.09 1.6 1.775 3.32 1.6
1.775 3.09 3.928 1.775 3.32 4.29
1.9595 3.09 3.928 1.9925 3.32 4.29
1.9595 3.0162 3.928 1.9925 3.09 4.29
1.9595 3.0162 3.98704 1.9925 3.09 3.859
19 − 2 × 1.775 − 0 × 1.6
x2 = = 3.09 .
5
102 3. MÉTODOS ITERATIVOS
En el método de GS:
45 + 2 × 1.775 − 3 × 3.09
x3 = = 3.928 .
10
En el método de Jacobi:
45 + 2 × 1.2 − 3 × 1.5
x3 = = 4.29 .
10
xk+1 = M xk + p. (3.11)
||M || < 1.
Para el método GS
MGS = −(D + L)−1 U,
pGS = (D + L)−1 b.
ri = bi − Ai·x ,
ri
δi = , (3.12)
aii
xi = xi + ωδi .
Una condición necesaria para que el método SOR converja, ver [Dem97], es
que
0 < ω < 2.
5 −1 2 −2 25
0 4 2 3 −10
A=
3
, b=
35 .
3 8 −2
−1 4 −1 6 −33
3.7. MÉTODO DE SOBRERRELAJACIÓN 105
Entonces
r1 = b1 − A1·x = 25 − 4 = 21
21
δ1 = = 4.2
5
ωδ1 = 5.88
x1 = 1 + 5.88 = 6.88
r2 = −10 − 9 = −19
−19
δ2 = = −4.75
4
ωδ2 = −6.65
x2 = 1 − 6.65 = −5.65
r3 = 35 − 9.69 = 25.31
25.31
δ3 = = 3.163750
8
ωδ3 = 4.429250
x3 = 1 + 4.429250 = 5.429250
r1 = 25 − 50.817517 = −25.817517
−25.817517
δ1 = = −5.163503
5
ωδ1 = −7.228905
x1 = 6.880000 + −7.228905 = −0.348905
Sobrerrelajación, ω = 1.4.
k x1 x2 x3 x4
0 1.000000 1.000000 1.000000 1.000000
1 6.880000 -5.650000 5.429250 0.045492
2 -0.348905 -5.088241 6.823724 -1.458380
3 1.076876 -4.710011 4.792473 -1.351123
4 1.810033 -3.552048 4.649676 -2.337041
5 1.368852 -2.880061 4.240550 -2.768266
6 1.721105 -2.409681 3.821389 -3.050409
7 1.788640 -2.008170 3.644054 -3.337915
8 1.812353 -1.742759 3.462571 -3.507443
9 1.883878 -1.543881 3.333868 -3.638593
10 1.909584 -1.395632 3.248121 -3.738508
11 1.932877 -1.289998 3.179762 -3.807650
12 1.952699 -1.211802 3.131447 -3.859624
13 1.964616 -1.154687 3.096340 -3.897553
14 1.974261 -1.113133 3.070228 -3.925007
15 1.981287 -1.082649 3.051371 -3.945238
Gauss-Seidel
k x1 x2 x3 x4
0 1.000000 1.000000 1.000000 1.000000
1 5.200000 -3.750000 4.081250 -1.453125
2 2.036250 -3.450781 4.542168 -2.103076
3 1.651746 -3.193777 4.427492 -2.357609
4 1.647204 -2.945539 4.272474 -2.549694
5 1.682025 -2.723966 4.128304 -2.715634
6 1.717631 -2.527427 3.999765 -2.862150
3
7 1.749749 -2.353270 3.885783 -2.991898
8 1.778274 -2.198968 3.784786 -3.106845
9 1.803554 -2.062259 3.695303 -3.208684
10 1.825953 -1.941139 3.616023 -3.298912
11 1.845798 -1.833828 3.545783 -3.378851
12 1.863381 -1.738753 3.483552 -3.449676
13 1.878958 -1.654519 3.428416 -3.512425
14 1.892760 -1.579890 3.379568 -3.568019
15 1.904987 -1.513770 3.336289 -3.617274
2
ωopt = π
1 + sin
m+1
Se puede mostrar que el método SOR se puede expresar como una iteracion
108 3. MÉTODOS ITERATIVOS
Ax = b
D + ωL (ω − 1)D + ωU
+ x=b
ω ω
D + ωL + (ω − 1)D + ωU x = ωb
(D + ωL)x = − (ω − 1)D + ωU x + ωb
(D + ωL)x = (1 − ω)D − ωU x + ωb
x = (D + ωL)−1 (1 − ω)D − ωU x + ω(D + ωL)−1 b
La siguiente tabla nos muestra los valores del número de iteraciones y del
radio espectral para diferentes valores de ω. El criterio de parada utilizado
fue max{|δi | : i = 1, ..., n} ≤ 0.000001.
3.7. MÉTODO DE SOBRERRELAJACIÓN 109
ρ(M )
1
ω
1 2
ω k ρ(M )
0.10 999 0.994
0.20 641 0.987
0.30 415 0.979
0.40 301 0.970
0.50 232 0.961
0.60 185 0.950
0.70 151 0.937
0.80 125 0.923
0.90 105 0.906
1.00 88 0.886
1.10 74 0.862
1.20 61 0.831
1.30 50 0.790
1.40 40 0.731
1.50 29 0.620
1.60 33 0.662
1.70 50 0.765
1.80 92 0.867
1.90 408 0.969
SOR: SOBRERRELAJACIÓN
datos: A, b, ω, x0 , ε, maxit
x = x0
para k = 1, ...,maxit
difX = 0
para i = 1, ..., n
ri = bi − Ai· x
ri
δi =
aii
xi = xi + ωδi
difX = max{difX, |ωδi |}
fin-para i
si difX ≤ ε ent x∗ ≈ x, salir
fin-para k
Ax = b (3.13)
∇f (x) = f ′ (x) = Ax − b = 0.
Ax = b
es equivalente a resolver
AT Ax = AT b
y es equivalente a minimizar
1
f (x) = xT AT Ax − (AT b)T x.
2
La matriz AT A es definida positiva, luego siempre se puede pensar en re-
solver un sistema de ecuaciones donde la matriz es definida positiva, prob-
lema equivalente a minimizar una función cuadrática estrictamente convexa
(3.14).
f ′ (xk )T dk < 0.
Esto garantiza que la dirección sea de descenso, es decir, que para t suficien-
temente pequeño
f (xk + tdk ) < f (xk ).
El segundo paso consiste en encontrar el mejor t posible, o sea, encontrar
xk+1 = xk + tk dk .
1
ϕ(t) = (xk + tdk )T A(xk + tdk ) − bT (xk + tdk )
2
t2 T T
ϕ(t) = dk Adk + tdk (Axk − b) + f (xk )
2
T T
ϕ′ (t) = tdk Adk + dk (Axk − b)
entonces T
dk (Axk − b)
tk = tc = − T (3.16)
dk Adk
dk = −f ′ (xk )
= −(Axk − b)
Entonces
dk = b − Axk (3.17)
T
dk dk
tk = T (3.18)
dk Adk
xk+1 = xk + tk dk . (3.19)
Ejemplo 3.11. Aplicar el método del descenso más pendiente para resolver
Ax = b, sabiendo que A es definida positiva, donde
4 1 2 13 1
A = 1 5 −2 , b = −21 , x0 = 1 .
2 −2 10 50 1
k = 0
d : 6.000000 -25.000000 40.000000
3.9. MÉTODO DEL DESCENSO MÁS PENDIENTE 113
t = 0.094488
x1 : 1.566927 -1.362196 4.779514
k = 1
d : -1.464541 -6.196916 -3.653391
t = 0.190401
x2 : 1.288078 -2.542093 4.083907
k = 2
d : 2.221969 -1.409801 1.500593
t = 0.135469
x3 : 1.589087 -2.733078 4.287191
k = 3
d : 0.802349 -0.349316 -1.516240
t = 0.164510
x4 : 1.721081 -2.790544 4.037754
k = 4
d : 0.830711 -0.692854 0.599209
t = 0.135907
x5 : 1.833980 -2.884707 4.119191
k = 5
d : 0.310405 -0.172063 -0.629281
t = 0.164543
x6 : 1.885055 -2.913019 4.015647
Ejemplo 3.12. Aplicar el método del descenso más pendiente para resolver
Ax = b, sabiendo que A es definida positiva, donde
19 6 8 55 1
A = 6 5 2 , b = 22 , x0 = 1 .
8 2 4 24 1
k = 0
d : 22.000000 9.000000 10.000000
t = 0.040905
x1 : 1.899920 1.368149 1.409055
k = 1
d : -0.579812 0.941625 0.428123
t = 0.531990
x2 : 1.591466 1.869085 1.636812
k = 2
d : 0.453147 -0.167842 0.982857
t = 0.089118
x3 : 1.631849 1.854127 1.724402
k = 3
d : -0.925117 -0.510535 0.339342
t = 0.068514
x4 : 1.568466 1.819148 1.747652
k = 4
d : 0.303036 -0.001843 0.823366
t = 0.091249
x5 : 1.596118 1.818980 1.822783
k = 5
d : -0.822384 -0.317174 0.301965
t = 0.069496
x6 : 1.538966 1.796938 1.843768
primer ejemplo
v =
2.3714059
5.5646277
11.063966
coc = 4.6655726
Segundo ejemplo:
valores propios
0.4250900
3.0722446
24.502665
coc = 57.641129
de la dirección anterior,
Para el caso de la solución de un sistema lineal por medio del método GC,
es corriente denominar el vector residuo
rk = Axk − b. (3.21)
α1 = 0, (3.22)
||rk ||22
αk = , k = 2, ..., n, (3.23)
||rk−1 ||22
||rk ||22
tk = T , k = 1, ..., n. (3.24)
dk Adk
• dk es dirección de descenso.
19 6 8 55
A = 6 5 2 , b = 22 .
8 2 4 24
118 3. MÉTODOS ITERATIVOS
x1 = x4 = (1, 2, 3),
r1 = (0, 0, 0).
para calcular tk . Las otras operaciones necesarias son producto escalar en-
tre vectores, sumas o restas de vectores y multiplicación de un escalar por
un vector. Todo esto hace que sea un método muy útil para matrices muy
grandes pero muy poco densas.
4
Solución de ecuaciones no
lineales
f (x) = 0,
f : R → R.
x5 − 3x4 + 10x − 8 = 0,
ex − x3 + 8 = 0,
x2 + x
− x = 0.
cos(x − 1) + 2
120
121
|xk − x∗ | ≤ ε,
|f (xk )| ≤ ε.
4.1 En Scilab
Para resolver
f (x) = 0,
x − ex
− cos(x) + 0.1 = 0,
1 + x2
es necesario definir una función de Scilab donde esté f y después utilizar
fsolve.
function fx = func156(x)
fx = ( x - exp(x) )/( 1 + x*x ) - cos(x) + 0.1
endfunction
function y = f123(x)
y = x*x*x - 4*x*x + 10*x - 20
endfunction
//------------------------------------------------
function d1 = der123(x)
d1 = 3*x*x - 8*x +10
endfunction
function y = f13(x)
y = exp(x) - 2.7*x
endfunction
x = fsolve(1, f13)
info =
4.
fx =
0.0182285
x =
0.9957334
y = f (x)
(x0 , f (x0 )) b
(x1 , f (x1 )) b
x2 x1 x0
Entonces
f (xk )
xk+1 = xk − (4.1)
f ′ (xk )
k xk f (xk ) f ′ (xk )
0 3.000000 2.200000E+01 91.000000
1 2.758242 5.589425E+00 47.587479
2 2.640786 9.381331E-01 32.171792
3 2.611626 4.892142E-02 28.848275
4 2.609930 1.590178E-04 28.660840
5 2.609924 1.698318E-09 28.660228
6 2.609924 -2.838008E-15 28.660227
Las raı́ces reales del polinomio x5 − 3x4 + 10x − 8 son: 2.6099, 1.3566, 1.
Tomando otros valores iniciales el método converge a estas raı́ces. Si se toma
x0 = 2.1, se esperarı́a que el método vaya hacia una de las raı́ces cercanas,
2.6099 o 1.3566 . Sin embargo, hay convergencia hacia 1.
k xk f (xk ) f ′ (xk )
0 2.100000 -4.503290e+00 -3.891500
1 0.942788 -1.974259e-01 3.894306
2 0.993484 -1.988663e-02 3.103997
3 0.999891 -3.272854e-04 3.001745
4 1.000000 -9.509814e-08 3.000001
5 1.000000 -7.993606e-15 3.000000
• Sencillez.
• Generalmente converge.
• Puede no converger.
número máximo de iteraciones, llamado por ejemplo maxit. Para una pre-
cisión εf , la detención deseada para el proceso iterativo se tiene cuando
|f (xk )| ≤ εf . Otra detención posible se da cuando dos valores de x son casi
iguales, es decir, cuando |xk −xk−1 | ≤ εx . Se acostumbra a utilizar el cambio
relativo, o sea, |xk − xk−1 |/|xk | ≤ εx . Para evitar las divisiones por cero,
se usa |xk − xk−1 |/(1 + |xk |) ≤ εx . Finalmente, siempre hay que evitar las
divisiones por cero o por valores casi nulos. Entonces, otra posible parada,
no deseada, corresponde a |f ′ (xk )| ≤ ε0 . El algoritmo para el método de
Newton puede tener el siguiente esquema:
x = x0
for k=0:maxit
[fx, der] = func(x)
//printf(’%3d %10.6f %10.6f %10.6f\n’, k, x, fx, der)
x = x - fx/der
end
ind = 2
endfunction
lim xk = x∗ , (4.2)
k→∞
|xk+1 − x∗ | |f ′′ (x∗ )|
lim = (4.3)
k→∞ |xk − x∗ |2 2|f ′ (x∗ )|
(x − xn )2
f (x) = f (xn ) + f ′ (xn )(x − xn ) + f ′′ (ξ) , ξ ∈ I(x, xn )
2
tomando x = x∗
(x∗ − xn )2
f (x∗ ) = 0 = f (xn ) + f ′ (xn )(x∗ − xn ) + f ′′ (ξ) , ξ ∈ I(x∗ , xn )
2
f (xxn ) ′′
∗ ∗ 2 f (ξ)
0= + (x − xn ) + (x − x n ) ,
f ′ (xxn ) 2f ′ (xn )
f (xxn ) f ′′ (ξ)
0 = x∗ − xn − ′ + (x∗ − xn )2 ′ ,
f (xxn ) 2f (xn )
f ′′ (ξ)
0 = x∗ − xn+1 + (x∗ − xn )2 ′ ,
2f (xn )
f ′′ (ξ)
x∗ − xn+1 = −(x∗ − xn )2 ′ . (4.4)
2f (xn )
Sea
I = [x∗ − ε, x + ε]
max |f ′′ (x)|
x∈I
M=
2 min |f ′ (x)|
x∈I
4.2. MÉTODO DE NEWTON 129
En particular,
|x∗ − x1 | = M |x∗ − x0 |2 ,
M |x∗ − x1 | = (M |x∗ − x0 |)2 .
Entonces
M |x∗ − x1 | < 1 ,
M |x∗ − x1 | < (M |x∗ − x0 |)2 ,
M |x∗ − x1 | < M |x∗ − x0 |, ya que 0 < t < 1 ⇒ t2 < t,
|x∗ − x1 | < ε,
..
.
M |x∗ − xn | < 1 ,
|x∗ − xn | < ε .
Luego
|x∗ − xn | ≤ M |x∗ − xn−1 |,
M |x∗ − xn | ≤ (M |x∗ − xn−1 |)2 ,
n
M |x∗ − xn | ≤ (M |x∗ − x0 |)2 ,
1 n
|x∗ − xn | ≤ (M |x∗ − x0 |)2 ,
M
lim |x∗ − xn | = 0,
n→∞
130 4. SOLUCIÓN DE ECUACIONES NO LINEALES
es decir
lim xn = x∗ .
n→∞
Reescribiendo (4.4),
x∗ − xn+1 f ′′ (ξ)
= − , ξ ∈ I(x∗ , xn )
(x∗ − xn )2 2f ′ (xn )
x∗ − xn+1 f ′′ (x∗ )
lim = − .
n→∞ (x∗ − xn )2 2f ′ (x∗ )
y = f (x)
(x0 , f (x0 )) b
(x1 , f (x1 )) b
x2 x1 x0
k xk f (xk )
0 3.000000 2.200000e+01
1 3.010000 2.292085e+01
2 2.761091 5.725624e+00
3 2.678210 2.226281e+00
4 2.625482 4.593602e-01
5 2.611773 5.317368e-02
6 2.609979 1.552812e-03
7 2.609925 5.512240e-06
8 2.609924 5.747927e-10
9 2.609924 -2.838008e-15
//*************
eps0 = 1.0e-12
//*************
4.4. MÉTODO DE LA BISECCIÓN 133
x = x0
h = 0.1
x1 = x0 + h
f0 = f(x0)
f1 = f(x1)
for k=1:maxit
den = f1-f0
if abs(den) <= eps0
ind = 0
return
end
d2 = f1*(x1-x0)/den
x2 = x1 - d2
f2 = f(x2)
disp(k,x2,f2)
if abs(f2) <= epsf | abs(d2) <= epsx
x = x2
ind = 1
return
end
x0 = x1, f0 = f1
x1 = x2, f1 = f2
end
x = x2
ind = 2
endfunction
ek = |xk − x∗ |.
bk − ak ≤ ε,
k
1
(b0 − a0 ) ≤ ε,
2
k
1 ε
≤ ,
2 b0 − a0
b0 − a0
2k ≥ ,
ε
b0 − a0
k log 2 ≥ log ,
ε
log b0 −ε
a0
k ≥ ·
log 2
4.5. MÉTODO DE REGULA FALSI 135
y = f (x)
(bk , f (bk )) b
ak ck
b
bk
f (ak )(bk − ak )
ck = ak − (4.7)
f (bk ) − f (ak )
• f (ak )f (ck ) < 0; en este caso hay una raı́z en el intervalo [ak , ck ] =
[ak+1 , bk+1 ];
• f (ak )f (ck ) > 0; en este caso hay una raı́z en el intervalo [ck , bk ] =
[ak+1 , bk+1 ].
lim (bk − ak ) = 0.
k→∞
Los dos métodos, bisección y Regula Falsi, se pueden combinar en uno solo
de la siguiente manera. En cada iteración se calcula mk y ck . Esto define
tres subintervalos en [ak , bk ]. En por lo menos uno de ellos se tiene una raı́z.
Si los tres subintervalos sirven, se puede escoger cualquiera, o mejor aún, el
de menor tamaño. En un caso muy especial, cuando mk y ck coinciden, se
tiene simplemente una iteración del método de bisección.
nmodif = 2 + 2 kmodif ·
g(x) = x. (4.8)
x0 = -1
x1 = -0.833333
x2 = -0.852623
x3 = -0.851190
x4 = -0.851303
x5 = -0.851294
x6 = -0.851295
x7 = -0.851295
x8 = -0.851295
x0 = -0.8510
x1 = -0.8488
x2 = -0.8294
x3 = -0.6599
x4 = 1.1415
x5 = -11.6832
x6 = -142.0691
x7 = -2.0190e+4
En este caso se observa que, aun partiendo de una muy buena aproximación
de la solución, no hay convergencia. 3
Antes de ver un resultado sobre convergencia del método de punto fijo, ob-
servemos su interpretación gráfica. La solución de g(x) = x está determinada
por el punto de corte, si lo hay, entre las gráficas y = g(x) y y = x.
Después de dibujar las dos funciones, la construcción de los puntos x1 , x2 ,
x3 ... se hace de la siguiente manera. Después de situar el valor x0 sobre el
eje x, para obtener el valor x1 , se busca verticalmente la curva y = g(x). El
punto obtenido tiene coordenadas (x0 , g(x0 )), o sea, (x0 , x1 ). Para obtener
x2 = g(x1 ) es necesario inicialmente resituar x1 sobre el eje x, para lo cual
basta con buscar horizontalmente la recta y = x para obtener el punto
140 4. SOLUCIÓN DE ECUACIONES NO LINEALES
y=x
y = g(x)
x∗
y=x
y = g(x)
x0 x1 x2 x3 x∗
y=x
y = g(x)
x0 x2 x∗ x3 x1
y = g(x) y=x
x∗ x0 x2 x3 x4
y = g(x) y=x
x3 x1 x∗ x0 x2 x4
x2 = 3,
x2 + x2 = x2 + 3,
x2 + 3
x = ,
2x
x + 3/x
x = ·
2
x0 = 3
x1 = 2
x2 = 1.75000000000000
x3 = 1.73214285714286
x4 = 1.73205081001473
x5 = 1.73205080756888
x6 = 1.73205080756888
g(1) = 2,
g(4) = 2.375,
1 3
g ′ (x) = − 2,
√ 2 2x
g ′ ( 3) = 0,
g ′ (1) = −1,
13
g ′ (4) = ,
32
3
g ′′ (x) = ·
x3
Entonces g ′′ (x) > 0 para todo x positivo. Luego g ′ (x) es creciente para
x > 0. Como g ′ (1) = −1, entonces −1 < g ′ (1 + ε). De nuevo por ser g ′ (x)
creciente, entonces −1 < g ′ (x) ≤ 13/32 para todo x ∈ I. En resumen,
|g ′ (x)| < 1 cuando x ∈ I.
√ √
Entre 1 + ε y 3 el valor de g ′ (x) es negativo.√ Entre 3 y√ 4 el
valor de g ′ (x) es positivo.
√ Luego g decrece
√ en√[1 + ε, √ 3] y crece √en [ 3, 4].
Entonces g([1 + ε, 3]) = [g(1 √ + ε), 3] ⊆ [2, 3] y g([ 3, 4]) = [ 3, 2.375].
En consecuencia g(I)√ = [ 3, 2.375] ⊆ I. Entonces el método de punto
fijo converge a x∗ = 3 para cualquier x√ 0 ∈]1, 4]. Este resultado se puede
generalizar al intervalo [1 + ε, b] con b > 3.
Si se empieza con x0 = 1/2, no se cumplen las condiciones del teorema; sin
embargo, el método converge. 3
y = mx + b
g(xj ) − g(xi )
m= (4.10)
xj − xi
g(xi ) = mxi + b
b = g(xi ) − mxi (4.11)
xk = mxk + b
b
xk = . (4.12)
1−m
Ahora se usan los puntos (xj , g(xj )), (xk , g(xk )), para obtener un nuevo xm ,
y ası́ sucesivamente. Usualmente, j = i + 1 y k = j + 1.
0 = cf (x),
x = x + cf (x) = g(x). (4.13)
1 + cf ′ (x∗ ) = 0,
1
c = − .
f ′ (x∗ )
f (xk )
xk+1 = xk − . (4.14)
f ′ (x∗ )
x21 + x1 x2 + x3 − 3 = 0
2x1 + 3x2 x3 − 5 = 0 (4.15)
2
(x1 + x2 + x3 ) − 10x3 + 1 = 0.
F1 (x1 , x2 , ..., xn ) = 0
F2 (x1 , x2 , ..., xn ) = 0
..
.
Fn (x1 , x2 , ..., xn ) = 0,
F (x) = 0. (4.16)
f (xk )
xk+1 = xk − ,
f ′ (xk )
f (xk )
xk+1 = xk − .
f ′ (xk )
148 4. SOLUCIÓN DE ECUACIONES NO LINEALES
El método empieza con un vector x0 = (x01 , x02 , ..., x0n ), aproximación inicial
de la solución x∗ . Mediante (4.17) se construye una sucesión de vectores
{xk = (xk1 , xk2 , ..., xkn )} con el deseo de que xk → x∗ . En palabras, el vector
xk+1 es igual al vector xk menos el producto de la inversa de la matriz
jacobiana F ′ (xk ) y el vector F (xk ). Para evitar el cálculo de una inversa, la
fórmula se puede reescribir
x21 + x1 x2 + x3 − 3 = 0
2x1 + 3x2 x3 − 5 = 0
(x1 + x2 + x3 )2 − 10x3 + 1 = 0
k=3
−0.1213 2.9718 0.7833 1.0000 0.1824 0.9658
1.4751 , 2.0000 3.4931 4.2156 , −0.3454 , 1.0598
0.5981 6.7057 6.7057 −3.2943 −0.1502 1.0141
150 4. SOLUCIÓN DE ECUACIONES NO LINEALES
k=4
−0.0297 2.9913 0.9658 1.0000 0.0335 0.9993
0.1557 , 2.0000 3.0424 3.1793 , −0.0587 , 1.0011
0.0981 6.0793 6.0793 −3.9207 −0.0139 1.0002
k=5
−0.0008 2.9997 0.9993 1.0000 0.0007 1.0000
0.0025 , 2.0000 3.0006 3.0033 , −0.0011 , 1.0000
0.0015 6.0012 6.0012 −3.9988 −0.0002 1.0000
0 1
F (x6 ) ≈ 0 , luego x∗ ≈ 1 . 3
0 1
Este método sirve para hallar raı́ces reales o complejas de polinomios. Sea
p(x) un polinomio real (con coeficientes reales), de grado n, es decir,
p(x) = a0 + a1 x + a2 x2 + ... + an xn , ai ∈ R, i = 0, 1, ..., n, an 6= 0.
En general no se puede garantizar que p(x) tenga raı́ces reales. Sin embargo
(teorema fundamental del Álgebra) se puede garantizar que tiene n raı́ces
complejas (algunas de ellas pueden ser reales). De manera más precisa,
existen r1 , r2 , ..., rn ∈ C tales que
p(ri ) = 0, i = 1, 2, ..., n.
El polinomio p se puede expresar en función de sus raı́ces:
p(x) = an (x − r1 )(x − r2 ) · · · (x − rn ).
Las raı́ces complejas, no reales, siempre vienen por parejas, es decir si r =
a + ib, b 6= 0, es una raı́z entonces r̄ = a − ib, el conjugado de r, también es
raı́z. Esto garantiza que los polinomios de grado impar tienen por lo menos
una raı́z real. Para los polinomios de grado par, el número de raı́ces reales
es par y el número de raı́ces estrictamente complejas también es par. Ası́
un polinomio de grado par puede tener cero raı́ces reales.
D = b2 − 4ac,
√
R= D
−b ± R −b ∓ R
x3 − x2 =
2a −b ∓ R
b − R2
2 b2 − b2 + 4ac 2c
x3 − x2 = = =
2a − b ∓ R 2a − b ∓ R −b ∓ R
2c
x3 − x2 = −
b±R
2c
x3 = x2 − (4.20)
b + signo(b)R
Si en una iteración
D = b2 − 4ac < 0
es necesario utilizar, a partir de ahı́, aritmética compleja (Scilab lo hace
automáticamente). Eso hace que los siguientes valores a, b y c no sean
necesariamente reales. Muy posiblemente b2 − 4ac tampoco es real. Para
utilizar (4.20) es necesario obtener la raı́z cuadradada de un complejo.
Sean z un complejo, θ el ángulo (en radianes) formado con el eje real (“eje
x”), llamado con frecuencia argumento de z, y ρ la norma o valor absoluto
de z. La dos raı́ces cuadradas de z son:
√ √
z = ζ1 = ρ cos(θ/2) + i sen(θ/2) ,
ζ2 = −ζ1 .
ρ = 20,
θ = tan−1 (16/12) = 0.927295,
√
ζ1 = 20 cos(0.927295/2) + i sen(0.927295/2) = 4 + 2i,
ζ2 = −4 − 2i. 3
4.9. MÉTODO DE MULLER 153
D = b2 − 4ac
√
R= D
(
b+R si |b + R| ≥ |b − R| (4.21)
δ=
b−R si |b + R| < |b − R|
2c
x3 = x2 − ·
δ
Ejemplo 4.9. Hallar las raı́ces de p(x) = 2x5 + x4 + 4x3 + 19x2 − 18x + 40
partiendo de x0 = 0, x1 = 0.5, x2 = 1.
f (x0 ) = 40
f (x1 ) = 36.375
f (x2 ) = 48
d = −0.25
a = 30.5
b = 38.5
c = 48
D = −4373.75
R = 66.134333i
δ = 38.5 + 66.134333i
x3 = 0.368852 + 1.084169i
f (x3 ) = 12.981325 − 9.579946i
154 4. SOLUCIÓN DE ECUACIONES NO LINEALES
Ahora utilizamos x1 , x2 y x3
d = 0.546325 + 0.413228i
a = 27.161207 + 11.293018i
b = −21.941945 + 50.286087i
c = 12.981325 − 9.579946i
D = −3890.341507 − 1752.330850i
R = 13.719321 − 63.863615i
δ = −35.661266 + 114.149702i
x4 = 0.586513 + 1.243614i
f (x4 ) = 3.760763 − 6.548104i
..
.
x5 = 0.758640 + 1.246582i
f (x5 ) = −2.013839 − 1.490220i
x6 = 0.748694 + 1.196892i
f (x6 ) = 0.123017 + 0.025843i
x7 = 0.750002 + 1.198942i
f (x7 ) = 0.000535 + 0.000636i
x8 = 0.750000 + 1.198958i
f (x8 ) = 0
Ahora se trabaja con p(x) = 2x3 + 4x2 + 6x2 + 20. Sean x0 = −3, x1 = −2.5
y x2 = −2. También se hubiera podido volver a utilizar x0 = 0, x1 = 0.5 y
4.9. MÉTODO DE MULLER 155
x2 = 1.
f (x0 ) = −16
f (x1 ) = −1.25
f (x2 ) = 8
d = −0.25
a = −11
b = 13
c=8
D = 521
R = 22.825424
δ = 35.825424
x3 = −2.446610
f (x3 ) = −0.026391
Ahora utilizamos x1 , x2 y x3
d = 0.011922
a = −9.893220
b = 22.390216
c = −0.026391
D = 500.277428
R = 22.366882
δ = 44.757098
x4 = −2.445431
f (x4 ) = −0.000057
..
.
x5 = −2.445428
f (x5 ) = 0
Ahora se trabaja con p(x) = 2x2 − 0.890857x + 8.178526. Sus raı́ces son
0.2227142 + 2.009891i y 0.2227142 − 2.009891i. En resumen, las 5 raı́ces de
p(x) son:
0.75 + 1.1989579i
0.75 − 1.1989579i
− 2.445428
0.222714 + 2.009891i
0.222714 − 2.009891i. 3
MÉTODO DE MULLER
datos: p, x0 , εf , ε0 , maxit
r = x0 , h = 0.5
mientras grado(p) ≥ 3
x0 = r, x1 = x0 + h, x2 = x1 + h
(r, inf o) = M uller1(p, xo , x1 , x2 , εf , ε0 , maxit)
si inf o = 0, ent parar
si |imag(r)| ≤ ε0 ent q(x) = (x − r)
sino q(x) = (x − r)(x − r̄)
p(x) = p(x)/q(x)
fin-mientras
calcular raı́ces de p (de grado no superior a 2)
reescrito como
R(d, e) = 0
S(d, e) = 0
Sea
q(x) = βn−2 xn−2 + βn−3 xn−3 + ... + β1 x + β0
reescrito como
el cociente. Entonces
p(x) = q(x)(x2 − dx − e) + Rx + S.
Es decir,
u1 xn + u2 xn−1 + ... + un x + un+1 = (v1 xn−2 + v2 xn−3 + ... + vn−2 x + vn−1 )(x2 − dx − e) + R
4.10. MÉTODO DE BAIRSTOW 159
u 1 = v1
u2 = v2 − dv1
u3 = v3 − dv2 − ev1
u4 = v4 − dv3 − ev2
ui = vi − dvi−1 − evi−2
vn = R
vn+1 = S + dvn
Entonces:
un = vn − dvn−1 − evn−2
un+1 = dvn − dvn − evn−1 + S
o sea, un+1 = vn+1 − dvn − evn−1
u 1 = v1
u2 = v2 − dv1
ui = vi − dvi−1 − evi−2 , i = 3, ..., n + 1.
v1 = u 1
v2 = u2 + dv1 (4.23)
vi = ui + dvi−1 + evi−2 , i = 3, ..., n + 1.
R = vn
S = vn+1 − dvn
u1 u2 u3 u4 ··· un+1
d dv1 dv2 dv3 ··· dvn
e ev1 ev2 ··· evn−1
v1 = u 1 v2 = Σ v3 = Σ v4 = Σ vn+1 = Σ
R = vn , S = vn+1 − dvn
4 5 1 0 -1 2
2 8 26 30 -18 -128
-3 -12 -39 -45 27
4 13 15 -9 -64 -99
vn (d, e) = 0
vn+1 (d, e) = 0
∆dk vn (dk , ek )
resolver el sistema J =− (4.24)
∆ek vn+1 (dk , ek )
k+1 k k
d d ∆d
k+1 = k + (4.25)
e e ∆ek
∂vn (dk , ek ) ∂vn k k
(d , e )
∂d ∂e
J =
.
∂vn+1 ∂vn+1 k k
(dk , ek ) (d , e )
∂d ∂e
∂v1
=0
∂d
∂v2
= v1
∂d
∂vi ∂vi−1 ∂vi−2
= vi−1 + d +e
∂d ∂d ∂d
∂v1
=0
∂e
∂v2
=0
∂e
∂vi ∂vi−1 ∂vi−2
=d + vi−2 + e
∂e ∂e ∂e
∂vi ∂vi−1 ∂vi−2
= vi−2 + d +e
∂e ∂e ∂e
∂v1
=0
∂d
∂v2
= v1
∂d
∂v3 ∂v2 ∂v1
= v2 + d +e
∂d ∂d ∂d
∂v3 ∂v2
= v2 + d
∂d ∂d
∂v4 ∂v3 ∂v2
= v3 + d +e
∂d ∂d ∂d
∂vi ∂vi−1 ∂vi−2
= vi−1 + d +e
∂d ∂d ∂d
Sea
w1 = v 1
w2 = v2 + dw1 (4.26)
wi = vi + dwi−1 + ewi−2 , i = 3, ..., n.
∂v1
=0
∂d
∂v2
= w1
∂d
∂v3
= w2
∂d
∂vi
= wi−1
∂d
∂v1
=0
∂e
∂v2
=0
∂e
∂v3
= v1
∂e
∂v4
= v2 + dv1
∂e
∂v5 ∂v4 ∂v3
= v3 + d +e
∂e ∂e ∂e
∂v1
=0
∂e
∂v2
=0
∂e
∂v3
= w1
∂e
∂v4
= w2
∂e
∂v5
= w3
∂e
∂vi
= wi−2
∂e
Entonces
∂vn
= wn−1
∂d
∂vn
= wn−2
∂e
∂vn+1
= wn
∂d
∂vn+1
= wn−1
∂e
164 4. SOLUCIÓN DE ECUACIONES NO LINEALES
MÉTODO DE BAIRSTOW
k = 0
4.0000 5.00001.0000 0.0000 -1.0000 2.0000
2.0000 8.0000
26.0000 30.0000 -18.0000 -128.0000
-3.0000 -12.0000 -39.0000 -45.0000 27.0000
------------------------------------------------------------
4.0000 13.0000 15.0000 -9.0000 -64.0000 -99.0000
2.0000 8.0000 42.0000 90.0000 36.0000
-3.0000 -12.0000 -63.0000 -135.0000
--------------------------------------------------
4.10. MÉTODO DE BAIRSTOW 165
Entonces
d = 1.0127362
e = −0.4344745
x2 − 1.0127362 x + 0.4344745 divide a p ,
r1 = 0.5063681 + 0.4219784 i es raı́z de p,
r2 = 0.5063681 − 0.4219784 i es raı́z de p,
q(x) = 4 x3 + 9.0509449 x2 + 8.4283219 x + 4.6032625 .
d = −0.9339455
e = −0.8660624
x2 + 0.9339455 x + 0.8660624 divide a p ,
r3 = −0.4669728 + 0.8049837 i es raı́z de p,
r4 = −0.4669728 − 0.8049837 i es raı́z de p,
q̃(x) = 4 x + 5.3151629 .
Ejercicios
4.1 x3 + 2x2 + 3x + 4 = 0.
4.2 x3 + 2x2 − 3x − 4 = 0.
4.6
3x − 6 x−2
− 2
cos(x) + 2 x + 1
+ x3 − 8 = 0.
x2 + x + 10
ex + x2
4.7
(1 + i)12
1000000 i = 945560.
(1 + i)12 − 1
4.9 x1 + x2 + 2x1 x2 − 31 = 0,
6x1 + 5x2 + 3x1 x2 − 74 = 0.
5
Interpolación y aproximación
Año Población
1930 3425
1940 5243
1950 10538
1960 19123
1970 38765
1980 82468
1985 91963
1990 103646
1995 123425
x1 f (x1 )
x2 f (x2 )
.. ..
. .
xn f (xn )
169
170 5. INTERPOLACIÓN Y APROXIMACIÓN
y
b
b
b
b
b
b
b
b
b
Cuando se desea que la función f˜ pase exactamente por los puntos conocidos,
f˜(xi ) = f (xi ) ∀i,
se habla de interpolación o de métodos de colocación, figura 5.2.
En los demás casos se habla de aproximación, figura 5.3. En este capı́tulo
se verá aproximación por mı́nimos cuadrados.
5.1. INTERPOLACIÓN 171
b
b
b
b
5.1 Interpolación
5.1.1 En Scilab
clear, clf
x = [ 0.5 1 1.5 2.1 3 3.6]’
y = [ 1 2 1.5 2.5 2.1 2.4]’
t = 0.8
ft = interpln( [x’; y’], t)
n = length(x);
xx = ( x(1):0.1:x(n) )’;
y1 = interpln( [x’; y’], xx);
plot2d(xx, y1)
3
b
b
b
2 b
1 b
1 2 3 4
clear, clf
x = [ 0.5 1 1.5 2.1 3 3.6]’
y = [ 1 2 1.5 2.5 2.1 2.4]’
n = length(x);
xx = ( x(1):0.1:x(n) )’;
d = splin(x, y);
ys = interp(xx, x, y, d);
plot2d(xx, ys)
3
b
b
b
2 b
1 b
1 2 3 4
Como las funciones de la base son conocidas, para conocer f˜ basta conocer
los escalares a1 , a2 , ..., an .
Las funciones de la base deben ser linealmente independientes. Si
n ≥ 2, la independencia lineal significa que no es posible que una de las
funciones sea combinación lineal de las otras. Por ejemplo, las funciones
ϕ1 (x) = 4, ϕ2 (x) = 6x2 − 20 y ϕ3 (x) = 2x2 no son linealmente independi-
entes.
Los escalares a1 , a2 , ..., an se escogen de tal manera que f˜(xi ) = yi , para
i = 1, 2, ..., n. Entonces
Φ a = y. (5.1)
174 5. INTERPOLACIÓN Y APROXIMACIÓN
Al plantear Φa = y, se tiene
1 −1 1 a1 1
1 2 4 a2 = −2
1 3 9 a3 5
Entonces
−4
a = −3 , f˜(x) = −4 − 3x + 2x2 ,
2
que efectivamente pasa por los puntos dados. 3
x = x(:);
y = y(:);
n = size(x,1);
n1 = n - 1;
F = ones(n,n);
for i=1:n1
F(:,i+1) = x.^i;
end
a = F\y
p = poly(a, ’x’, ’c’)
5.2. INTERPOLACIÓN POLINOMIAL DE LAGRANGE 175
xx = (x(1):0.05:x(n))’;
yp = horner(p, xx);
Hay ejemplos clásicos de los problemas que se pueden presentar con valores
relativamente pequeños, n = 20.
Ejemplo 5.2. Dados los puntos mismos (−1, 1), (2, −2), (3, 5) y la base
formada por las funciones ϕ1 (x) = 1, ϕ2 (x) = ex , ϕ3 (x) = e2x , encontrar la
función de interpolación.
Al plantear Φa = y, se tiene
1 0.3679 0.1353 a1 1
1 7.3891 54.5982 a2 = −2
1 20.0855 403.4288 a3 5
Entonces
−1.2921
a = −0.8123 , f˜(x) = 1.2921 − 0.8123ex + 0.0496e2x ,
0.0496
Los valores xi deben ser todos diferentes entre sı́. Sin perder generalidad, se
puede suponer que x1 < x2 < · · · < xn .
El problema 5.2 se puede resolver planteando n ecuaciones con n incógnitas
(los coeficientes del polinomio). Este sistema lineal se puede resolver y se
tendrı́a la solución. Una manera más adecuada de encontrar p es por medio
de los polinomios de Lagrange.
f (b) − f (a)
= f ′ (c).
b−a
f ′ (c) = 0.
(x − 2)(x − 3) x2 − 5x + 6
L1 (x) = = ,
(−1 − 2)(−1 − 3) 12
(x − −1)(x − 3) x2 − 2x − 3
L2 (x) = = ,
(2 − −1)(2 − 3) −3
(x − −1)(x − 2) x2 − x − 2
L3 (x) = = .
(3 − −1)(3 − 2) 4
n
X
p(x) = yk Lk (x). (5.5)
k=1
Para el ejemplo,
Obviamente p(1) = 3.8 y p(2) = 3.95. Sin embargo p(1.35) = 6.946. Ver
figura (5.6). 3
x = [-1 1 2 3]’;
n = length(x)
k = 2
4 b
b b b b
0
0 1 2 3 4
q(xi ) = yi , i = 1, 2, ..., n.
para algún ξ ∈ It .
180 5. INTERPOLACIÓN Y APROXIMACIÓN
Φ(x) = (x − x1 )(x − x2 ) · · · (x − xn ),
Φ(x)
G(x) = E(x) − E(t).
Φ(t)
Entonces
G ∈ CInt ,
Φ(xi )
G(xi ) = E(xi ) − E(t) = 0, i = 1, ..., n
Φ(t)
Φ(t)
G(t) = E(t) − E(t) = 0.
Φ(t)
Como G tiene por lo menos n + 1 ceros en It , aplicando el corolario del
teorema del valor medio, se deduce que G′ tiene por lo menos n + 1 − 1 ceros
en It . Ası́ sucesivamente se concluye que G(n) tiene por lo menos un cero en
It . Sea ξ tal que
G(n) (ξ) = 0.
De acuerdo con las definiciones
Ejemplo 5.5. Considere los valores de la función seno en los puntos 5, 5.2,
5.5 y 6. Sea p el polinomio de interpolación. Obtenga una cota para error
cometido al aproximar sen(5.8) por p(5.8). Compare con el valor real del
error.
se puede construir pn−2 ∈ Pn−2 . Sea c(x) la corrección que permite pasar
de pn−2 a pn−1 ,
pn−1 (x) = pn−2 (x) + c(x), es decir, c(x) = pn−1 (x) − pn−2 (x).
Esta fórmula es la que se utiliza para calcular pn−1 (x), una vez que se sepa
calcular, de manera sencilla, f [x1 , x2 , ..., xn ].
• A partir de pn−2 (x), con el valor f [x1 , x2 , ..., xn ], se calcula pn−1 (x).
5.3. DIFERENCIAS DIVIDIDAS DE NEWTON 183
Obviamente
p0 (x) = f (x1 ). (5.9)
f [x1 ] := f (x1 ),
que se generaliza a
f [xi ] := f (xi ), ∀i. (5.10)
Para x = x2 ,
p1 (x2 ) − po (x2 )
f [x1 , x2 ] = ,
x2 − x1
f (x2 ) − f (x2 )
f [x1 , x2 ] = ,
x2 − x1
f [x2 ] − f [x1 ]
f [x1 , x2 ] = .
x2 − x1
f [xi+1 ] − f [xi ]
f [xi , xi+1 ] = · (5.11)
xi+1 − xi
Deducción de f [x1 , x2 , x3 ] :
Entonces
(0, 0), (0.5, 0.7071), (1, 1), (2, 1.4142), (3, 1.7321), (4, 2).
186 5. INTERPOLACIÓN Y APROXIMACIÓN
para i = 1, ..., n
D0 f [xi ] = f (xi )
fin-para i
para j = 1, ..., m
para i = 1, ..., n − j
calcular Dj f [xi ] según (5.18)
fin-para i
fin-para j
x = x(:)
y = y(:)
n = size(x,1)
DD = zeros(n,m+1);
5.3. DIFERENCIAS DIVIDIDAS DE NEWTON 187
DD(:,1) = y;
for j=1:m
for i=1:n-j
Djfi = ( DD(i+1,j) - DD(i,j) )/( x(i+j) - x(i) );
DD(i,j+1) = Djfi;
end
end
disp(DD)
• Si para algún m todos los valores f [xk , xk+1 , ..., xk+m ] son iguales (o
aproximadamente iguales), entonces f es (aproximadamente) un poli-
nomio de grado m.
• Si para algún r todos los valores f [xk , xk+1 , ..., xk+r ] son nulos (o
aproximadamente nulos), entonces f es (aproximadamente) un poli-
nomio de grado r − 1.
xi
0
.5
1
2
3
4
xi + xi+m
ui = , (5.21)
2
xi + xi+1 + xi+2 + · · · + xi+m
vi = , (5.22)
m+1
|x̄ − uk | = min{|x̄ − ui | : x̄ ∈ [xi , xi+m ]}, (5.23)
i
|x̄ − vk | = min{|x̄ − vi | : x̄ ∈ [xi , xi+m ]}. (5.24)
i
Los valores ui y vi son, de alguna forma, indicadores del centro de masa del
intervalo [xi , xi+m ]. Con frecuencia, los dos criterios, (5.23) y (5.24), definen
el mismo xk , pero en algunos casos no es ası́. De todas formas son criterios
razonables y para trabajar se escoge un solo criterio, lo cual da buenos
resultados. Se puede preferir la utilización de vi que, aunque requiere más
operaciones, tiene en cuenta todos los xj pertenecientes a [xi , xi+m ].
xi ui |x̄ − ui | vi |x̄ − vi |
0
.5 1.25 0.44 1.1667 0.5233
√ √
1 2.00 0.31 2.0000 0.3100
2
3
4
xi ui |x̄ − ui | vi |x̄ − vi |
γ0 = 1, p0 (x) = 0.7071,
γ1 = 1(1.69−0.5) = 1.19, p1 (x) = 0.7071+0.5858(1.19)
p1 (x) = 1.404202
γ2 = 1.19(1.69−1) = 0.8211, p2 (x) = 1.404202−0.1144(0.8211)
p2 (x) = 1.310268
γ3 = 0.8211(1.69−2) =−0.254541, p3 (x) = 1.310268+0.0265(−0.254541)
p3 (x) = 1.303523. 3
determinar k
px = f (xk )
gi = 1.0
para j = 1, ..., m
gi = gi ∗ (x̄ − xk+j−1 )
px = px + gi ∗ Dj f [xk ]
fin-para j
n = length(x);
if t <= x(1)
k = 1
else if t >= x(n)
k = n-m;
else
5.3. DIFERENCIAS DIVIDIDAS DE NEWTON 191
distmin = 1.0e10;
k = -1;
for i=1:n-m
if ( x(i) <= t & t <= x(i+m) ) | m == 0
vi = mean(x(i:i+m));
di = abs( t - vi );
if di < distmin
distmin = di;
k = i;
end
end // if
end // for i
end // else
end // else
pt = DD(k,1)
gi = 1
for j=1:m
gi = gi*(t-x(k+j-1))
pt = pt + gi*DD(k,j+1)
end
Cuando los puntos (x1 , f (x1 )), (x2 , f (x2 )), (x3 , f (x3 )), ..., (xn , f (xn )), están
igualmente espaciados en x, es decir, existe un h > 0 tal que
xi = xi−1 + h, i = 2, ..., n
xi = x1 + (i − 1)h, i = 1, ..., n
∆0 fi = fi (5.25)
∆fi = fi+1 − fi (5.26)
k+1 k k k
∆ fi = ∆ (∆fi ) = ∆ fi+1 − ∆ fi (5.27)
k
X
k j k
∆ fi = (−1) fi+k−j , (5.28)
j
j=0
k
X k
fi+k = ∆j fi . (5.29)
j
j=0
D0 f [xi ] = f [xi ] = fi = ∆0 fi
fi+1 − fi ∆1 fi
D1 f [xi ] = f [xi , xi+1 ] = =
xi+1 − xi h
f [xi+1 , xi+2 ] − f [xi , xi+1 ] ∆2 fi
D2 f [xi ] = f [xi , xi+1 , xi+2 ] = = ··· =
xi+2 − xi 2h2
∆m fi
Dm f [xi ] = f [xi , ..., xi+m ] = (5.30)
m! hm
5.4. DIFERENCIAS FINITAS 193
xi fi ∆fi ∆2 fi ∆3 fi
x1 f1
∆f1
x2 f2 ∆2 f1
∆f2 ∆3 f1
x3 f3 ∆2 f2
∆f3 ∆3 f2
x4 f4 ∆2 f3
∆f4
x5 f5
xi fi ∆fi ∆2 fi ∆3 fi
0 0.0000
0.7071
.5 0.7071 −0.4142
0.2929 0.3460
1 1.0000 −0.0682
0.2247 0.0330
1.5 1.2247 −0.0352
0.1895 0.0126
2 1.4142 −0.0226
0.1669
2.5 1.5811
194 5. INTERPOLACIÓN Y APROXIMACIÓN
para i = 1, ..., n
∆0 fi = f (xi )
fin-para i
para j = 1, ..., m
para i = 1, ..., n − j
∆j fi = ∆j−1 fi+1 − ∆j−1 fi
fin-para i
fin-para j
Q
El valor i! se puede escribir i−1
j=0 (j + 1). Además, sea s = (x − xk )/h, es
decir, x = xk + sh. Entonces, x − xk+j = xk + sh − xk − jh = (s − j)h.
Xm i i−1
Y
pm (x) = ∆ fk (s − j)h
i! hi
i=0 j=0
Xm if Y i−1
∆ k
= (s − j)
i!
i=0 j=0
m
X i−1
Y s−j
= ∆i fk
j+1
i=0 j=0
xi + xi+m
ui = , i = 1, ..., n − m,
2
|x − uk | = min{|x − ui | : i = 1, ..., n − m}.
x − xk
s= .
h
γ0 = 1, p0 (x) = fk
γ1 = γ0 s, p1 (x) = p0 (x) + ∆1 fk γ1
γ2 = γ1 (s − 1)/2, p2 (x) = p1 (x) + ∆2 fk γ2
γ3 = γ2 (s − 2)/3, p3 (x) = p2 (x) + ∆3 fk γ3
γ4 = γ3 (s − 3)/4, p4 (x) = p3 (x) + ∆4 fk γ4
..
.
Ejemplo 5.9. Calcular p3 (1.96) y p2 (1.96) a partir de los puntos (0, 0),
(0.5, 0.7071), (1, 1), (1.5, 1.2247), (2, 1.4142), (2.5, 1.5811).
La tabla de diferencias finitas es la misma del ejemplo anterior. Para calcular
p3 (1.96) se tiene xk = x2 = 1. Entonces s = (1.96 − 1)/0.5 = 1.92 .
γ0 = 1, p0 (x) = f2 = 1
γ1 = 1(1.92) = 1.92, p1 (x) = 1 + .2247(1.92) = 1.431424
γ2 = 1.92(1.92 − 1)/2 = .8832, p2 (x) = 1.431424 − .0352(.8832)
p2 (x) = 1.400335
γ3 = γ2 (1.92 − 2)/3 = −.023552, p3 (x) = 1.400335 + .0126(−.023552)
p3 (x) = 1.400039
γ0 = 1, p0 (x) = f3 = 1.2247
γ1 = 1(0.92) = 0.92, p1 (x) = 1.2247 + .1895(.92) = 1.39904
γ2 = 0.92(0.92 − 1)/2 = −.0368, p2 (x) = 1.39904 − .0226(−0.0368)
p2 (x) = 1.399872
5.5. TRAZADORES CÚBICOS, INTERPOLACIÓN POLINOMIAL POR TROZOS, SPLINES197
3 •
2 • • • •
0
0 1 2 3 4 5
Para interpolar por polinomios de orden 2, si t < 2.5 se utilizan los puntos
(1, 2), (2, 2) y (3, 2). Entonces, por ejemplo, p2 (2.49) = 2. Si 2.5 < t < 3.5,
se utilizan los puntos (2, 2), (3, 2) y (4, 3). Después de algunos cálculos se ob-
tiene p2 (2.51) = 1.87505. Para t = 2.501 se obtiene p2 (2.501) = 1.8750005.
El lı́mite de p2 (t), cuando t → 2.5+ , es 1.875. Esto nos muestra una discon-
tinuidad. En t = 3.5 también se presenta una discontinuidad.
Estas discontinuidades se pueden evitar utilizando en el intervalo [1, 3] un
polinomio p2 (t) y en el intervalo [3, 5] otro polinomio p2 (t).
198 5. INTERPOLACIÓN Y APROXIMACIÓN
3 •
2 • • • •
0
0 1 2 3 4 5
3 •
2 • • • •
0
0 1 2 3 4 5
Si (x) = ai (x − xi )3 + bi (x − xi )2 + ci (x − xi ) + di , i = 1, 2, ..., n − 1.
(5.33)
Se requiere que S(x) pase por los puntos (xi , yi ), y que sea doblemente
diferenciable. Además, es necesario tratar algunos detalles adicionales en
los extremos de los intervalos. Entonces,
S(xi ) = yi , i = 1, ..., n,
Si (xi+1 ) = Si+1 (xi+1 ), i = 1, ..., n − 2,
Si′ (xi+1 ) = Si+1
′
(xi+1 ), i = 1, ..., n − 2,
Si′′ (xi+1 ) = Si+1
′′
(xi+1 ), i = 1, ..., n − 2.
di = yi i = 1, ..., n, (5.34)
ai h3i 2
+ bi hi + ci hi + di = di+1 i = 1, ..., n − 1, (5.35)
3ai h2i + 2bi hi + ci = ci+1 i = 1, ..., n − 2, (5.36)
3ai hi + bi = bi+1 i = 1, ..., n − 2. (5.37)
De (5.37):
bi+1 − bi
ai = (5.38)
3hi
Reemplazando (5.38) en (5.35):
h2i
(bi+1 − bi ) + bi h2i + ci hi + di = di+1
3
h2i
(bi+1 + 2bi ) + ci hi + di = di+1 (5.39)
3
Reemplazando (5.38) en (5.36):
(bi+1 − bi )hi + 2bi hi + ci = ci+1
(bi+1 + bi )hi + ci = ci+1 (5.40)
Despejando ci de (5.39):
1 hi
ci = (di+1 − di ) − (2bi + bi+1 ) (5.41)
hi 3
Cambiando i por i − 1:
1 hi−1
ci−1 = (di − di−1 ) − (2bi−1 + bi ) (5.42)
hi−1 3
5.5. TRAZADORES CÚBICOS, INTERPOLACIÓN POLINOMIAL POR TROZOS, SPLINES201
Multiplicando por 3:
3 3
hi−1 bi−1 + 2(hi−1 + hi )bi + hi bi+1 = (di−1 − di ) + (−di + di+1 ).
hi−1 hi
(5.44)
S1′′ (x1 ) = 0,
′′ (5.45)
Sn−1 (xn ) = 0.
Además del resultado anterior (b1 = 0), se puede introducir una variable
adicional bn = 0. Esto permite que la ecuación (5.44) se pueda aplicar para
i = n − 1. Recuérdese que ya se introdujo dn = yn y que para todo i se
tiene di = yi . Entonces se tiene un sistema de n ecuaciones con n incógnitas,
escrito de la forma
Ab = ζ, (5.49)
donde
1 0 0 0
h1 2(h1 + h2 ) h2 0
0 h2 2(h2 + h3 ) h3
0 0 h3 2(h3 + h4 ) h4
A=
0 0 h 2(h + h ) h
n−2 n−2 n−1 n−1
0 0 0 1
b1 0
3 3
(y − y ) + (−y + y )
b2 1 2 2 3
h1 h2
3 3
b3 (y2 − y3 ) + (−y3 + y4 )
h 2 h3
b= , ζ= .
.. .
.
. .
3 3
bn−1 (yn−2 − yn−1 ) + (−yn−1 + yn )
h hn−1
n−2
bn 0
Una vez conocidos los valores b1 , b2 , ..., bn−1 , bn , se puede aplicar (5.41)
para calcular los ci :
1 hi
ci = (yi+1 − yi ) − (2bi + bi+1 ), i = 1, ..., n − 1. (5.50)
hi 3
5.5. TRAZADORES CÚBICOS, INTERPOLACIÓN POLINOMIAL POR TROZOS, SPLINES203
• di = yi , i = 1, ..., n − 1.
• Obtener b1 , b1 , ..., bn resolviendo (5.49).
En particular b1 = 0 y, de acá en adelante, no se considera bn = 0.
• Para i = 1, ..., n − 1 calcular ci según (5.50).
• Para i = 1, ..., n − 1 calcular ai según (5.51).
Ejemplo 5.10. Construir el trazador cúbico para los puntos (1, 2), (2, 2),
(3, 2), (4, 3) y (5, 2).
3 •
2 • • • •
0
0 1 2 3 4 5
Figura 5.10: Interpolación con trazadores cúbicos o splines
Entonces
S1 (x) = −0.107143(x − 1)3 + 0(x − 1)2 + 0.107143(x − 1) + 2
S2 (x) = 0.535714(x − 2)3 − 0.321429(x − 2)2 − 0.214286(x − 2) + 2
S3 (x) = −1.035714(x − 3)3 + 1.285714(x − 3)2 + 0.75(x − 3) + 2
S4 (x) = 0.607143(x − 4)3 − 1.821429(x − 4)2 + 0.214286(x − 4) + 3 .
es decir,
m
X 2
min f˜(xi ) − yi .
i=1
Esto significa que se está buscando una función f˜, combinación lineal de
las funciones de la base, tal que minimiza la suma de los cuadrados de las
distancias entre los puntos (xi , f˜(xi )) y (xi , yi ).
Ejemplo 5.11. Dadas las funciones ϕ1 (x) = 1, ϕ2 (x) = x, ϕ3 (x) = x2 ,
encontrar la función f˜ que aproxima por mı́nimos cuadrados la función dada
por los puntos (0, 0.55), (1, 0.65), (1.5, 0.725), (2, 0.85), (3, 1.35).
Como las funciones de la base son 1, x, x2 , en realidad se está buscando
aproximar por mı́nimos cuadrados por medio de un parábola. El sistema
inicial es
1 0 0 0.55
1 1 1 a1 0.65
1 1.5 2.25 a2 ≈ 0.725
1 2 4 a3 0.85
1 3 9 1.35
Las ecuaciones normales dan:
5 7.5 16.25 a1 4.1250
7.5 16.25 39.375 a2 = 7.4875
16.25 39.375 103.0625 a3 17.8313
La solución es:
0.56
a = −0.04 , f˜(x) = 0.56 − 0.04x + 0.1x2 .
0.10
f˜(x1 ) 0.56 0.55
f˜(x2 )
0.62 0.65
˜
f (x3 ) = Φ a =
0.725 ,
y=
0.725 3
f˜(x4 ) 0.88 0.85
˜
f (x5 ) 1.34 1.35
5.6. APROXIMACIÓN POR MÍNIMOS CUADRADOS 207
Ejercicios
(1, −5),
(2, −4),
(4, 4).
(−1, −5),
(1, −5),
(2, −2),
(4, 40).
(−1, 10),
(1, 8),
(2, 4),
(4, −10).
5.7 Considere los mismos puntos del ejercicio anterior. Construya la tabla
de diferencias finitas hasta el orden 3. Halle p2 (0.11), p2 (0.08),
p2 (0.25), p2 (0.12), p2 (0.33), p2 (0.6), p3 (0.25), p3 (0.33), p3 (0.6).
5.9 Considere los mismos puntos del ejercicio anterior. Obtenga la parábola
de aproximación por mı́nimos cuadrados.
5.10 Considere los mismos puntos de los dos ejercicios anteriores. Use
otra base y obtenga la correspondiente función de aproximación por
mı́nimos cuadrados.
6
Integración y diferenciación
Esta técnica sirve para calcular el valor numérico de una integral definida,
es decir, parar obtener el valor
Z b
I= f (x)dx.
a
Z 0.5
2
ex dx.
0.1
209
210 6. INTEGRACIÓN Y DIFERENCIACIÓN
y = f (x)
a b
6.2 En Scilab
Para obtener una aproximación del valor de una integral definida, por ejem-
plo,
Z 0.5
2
e−x dx
0.1
function fx = f57(x)
fx = exp(-x*x)
endfunction
Z 2π
sen(x) dx .
0
b−a
xi = a + ih , i = 0, 1, 2, ..., n , h = ,
n
y supongamos conocidos yi = f (xi ).R Supongamos además que n es un
x
múltiplo de m, n = km. La integral x0n f (x)dx se puede separar en inter-
212 6. INTEGRACIÓN Y DIFERENCIACIÓN
x0 x1 x2 xm x2m xn−m xn
a b
Figura 6.2: División en subintervalos
Z xn Z xm Z x2m Z xn
f (x)dx = f (x)dx + f (x)dx + · · · + f (x)dx.
x0 x0 xm xn−m
Z xm Z xm
f (x)dx ≈ pm (x)dx.
x0 x0
x − x1 x − x0
p1 (x) = y0 + y1 ,
x0 − x1 x1 − x0
x−h x
p1 (x) = y0 + y1 ,
−h h
y1 − y 0
p1 (x) = y0 + x .
h
6.3. FÓRMULA DEL TRAPECIO 213
y0 y1
x0 x1
h
Figura 6.3: Fórmula del trapecio
Entonces
Z x1 Z h
y1 − y0
p1 (x)dx = (y0 + x )dx
x0 0 h
h2 y1 − y0
= y0 h + ,
2 h
y0 y1
= h( + ),
Z x1 2 2
y0 + y1
f (x)dx ≈ h · (6.1)
x0 2
De la fórmula (6.1) o de la gráfica se deduce naturalmente el nombre de
fórmula del trapecio.
Ejemplo 6.1.
Z 0.2
1 1
ex dx ≈ 0.2( e0 + e0.2 ) = 0.22214028 . 3
0 2 2
Ejemplo 6.2.
Z 0.8
1 1
ex dx ≈ 0.2( e0 + e0.2 + e0.4 + e0.6 + e0.8 ) = 1.22962334 . 3
0 2 2
= (f (x) − p1 (x))dx .
x0
h3
|eloc | ≤ M, M = max{|f ′′ (z)| : z ∈ [x0 , x1 ]}. (6.4)
12
En el ejemplo 6.1, f ′′ (x) = ex , max{|f ′′ (z)| : z ∈ [0, 0.2]} = 1.22140276,
luego el máximo error que se puede cometer, en valor absoluto, es (0.2)3 ×
1.22140276/12 = 8.1427 · 10−4 . En este ejemplo, se conoce el valor exacto
I = e0.2 − 1 = 0.22140276, luego |e| = 7.3752 · 10−4 .
En algunos casos, la fórmula del error permite afinar un poco más. Si
f ′′ (x) > 0 (f estrictamente convexa) en [x0 , x1 ] y como I = I˜+eloc , entonces
la fórmula del trapecio da un valor aproximado pero superior al exacto.
En el mismo ejemplo, f ′′ (x) varı́a en el intervalo [1, 1.22140276] cuando
x ∈ [0, 0.2]. Luego
entonces
I ∈ [0.22132601, 0.22147361].
El error global es el error correspondiente al hacer la aproximación de la
integral sobre todo el intervalo [x0 , xn ], o sea, el error en la fórmula 6.2,
Z xn
y0 yn
eglob = f (x)dx − h( + y1 + y2 + · · · + yn−2 + yn−1 + )
x0 2 2
Xn
f ′′ (zi ) h3
= (− ) , zi ∈ [xi−1 , xi ]
12
i=1
n
h3 X ′′
=− f (zi ) , zi ∈ [xi−1 , xi ]
12
i=1
216 6. INTEGRACIÓN Y DIFERENCIACIÓN
Sean
M1 = min{f ′′ (x) : x ∈ [a, b]} , M2 = max{f ′′ (x) : x ∈ [a, b]}.
Entonces
M1 ≤ f ′′ (zi ) ≤ M2 , ∀i
X n
nM1 ≤ f ′′ (zi ) ≤ nM2 ,
i=1
n
1 X ′′
M1 ≤ f (zi ) ≤ M2 .
n
i=1
2 , entonces, aplicando el teorema del valor intermedio a f ′′ , existe
Si f ∈ C[a,b]
ξ ∈ [a, b] tal que
n
′′ 1 X ′′
f (ξ) = f (zi ) .
n
i=1
Entonces
h3 ′′
eglob = −
nf (ξ) , ξ ∈ [a, b].
12
Como h = (b − a)/n, entonces n = (b − a)/h.
(b − a)f ′′ (ξ)
eglob = −h2 , ξ ∈ [a, b]. (6.5)
12
Z 2h
1 8h3 4h2
p2 (x)dx = (y 0 − 2y 1 + y 2 ) + h (−3y0 + 4y1 − y2 )
0 2h2 3 2
+ 2h2 (2h)y0 ,
Z 2h
1 4 1
p2 (x)dx = h( y0 + y1 + y2 ).
0 3 3 3
Entonces
Z x2
h
f (x)dx ≈ (y0 + 4y1 + y2 ) (6.6)
x0 3
Ejemplo 6.3.
Z 0.8
0.2 0
ex dx ≈ (e + 4(e0.2 + e0.6 ) + 2 e0.4 + e0.8 ) = 1.22555177 .
0 3
Z h Z h
e(h) = eloc (h) = f (x) dx − p2 (x) dx,
−h −h
Z h
h
= f (x) dx − f (−h) + 4f (0) + f (h) .
−h 3
218 6. INTEGRACIÓN Y DIFERENCIACIÓN
Rh
Sea F tal que F ′ (x) = f (x), entonces −h f (x) dx = F (h) − F (−h). Al
derivar con respecto a h se tiene f (h) + f (−h).
1
e′ (h) = f (h) + f (−h) − f (−h) + 4f (0) + f (h)
3
h
− − f (−h) + f ′ (h) ,
′
3
3e (h) = 2f (h) + 2f (−h) − 4f (0) − h(f ′ (h) − f ′ (−h)).
′
3e′′′ (h) = f ′′ (h) + f ′′ (−h) − (f ′′ (h) + f ′′ (−h)) − h(f ′′′ (h) − f ′′′ (−h)),
= −h(f ′′′ (h) − f ′′′ (−h)),
h
e′′′ (h) = − ( f ′′′ (h) − f ′′′ (−h) ),
3
2h 2 f ′′′ (h) − f ′′′ (−h)
e′′′ (h) = − .
3 2h
De los resultados anteriores se ve claramente que e(0) = e′ (0) = e′′ (0) =
e′′′ (0) = 0. Además, como f ∈ C 4 , entonces f ′′′ ∈ C 1 . Por el teorema del
valor medio, existe β ∈ [−h, h], β = αh, α ∈ [−1, 1], tal que
Z h
′′
e (h) = e′′′ (t) dt + e′′ (0),
0
Z
′′ 2 h 2
e (h) = − t g4 (t) dt.
3 0
6.4. FÓRMULA DE SIMPSON 219
M1 ≤ f (4) (zj ) ≤ M2 , ∀j
X k
kM1 ≤ f (4) (zj ) ≤ kM2 ,
j=1
Xk
1
M1 ≤ f (4) (zj ) ≤ M2 ,
k
j=1
Entonces
(b − a)f (4) (ξ)
eglob = −h4 , ξ ∈ [a, b]. (6.9)
180
6.5. OTRAS FÓRMULAS DE NEWTON-COTES 221
Z x3
h 3
f (x)dx = (3y0 + 9y1 + 9y2 + 3y3 ) − h5 f (4) (z) , z ∈ [x0 , x3 ],
x0 8 80
Z xm Z xm
f (x)dx ≈ pm (x)dx.
x0 x0
222 6. INTEGRACIÓN Y DIFERENCIACIÓN
m error
h f ′′ (z) 3
1 (y0 + y1 ) − h
2 12
h f (4) (z) 5
2 (y0 + 4y1 + y2 ) − h
3 90
3h 3 f (4) (z) 5
3 (y0 + 3y1 + 3y2 + y3 ) − h
8 80
2h 8 f (6) (z) 7
4 (7y0 + 32y1 + 12y2 + 32y3 + 7y4 ) − h
45 945
En todos los casos, z ∈ [x0 , xm ].
En general, las fórmulas cerradas son más precisas que las abiertas, en-
tonces, siempre que se pueda, es preferible utilizar las fórmulas cerradas.
Las fórmulas abiertas se usan cuando no se conoce el valor de la función
f en los extremos del intervalo de integración; por ejemplo, en la solución
numérica de algunas ecuaciones diferenciales ordinarias.
b − a p (q)
I = In + F (b − a)( ) f (ξ),
n
(b − a)p+1 (q)
= In + F f (ξ).
np
Sea m = 2n,
(b − a)p+1 (q)
I = Im + F f (ζ),
np 2p
Supongamos que
Entonces
224 6. INTEGRACIÓN Y DIFERENCIACIÓN
I ≈ In + 2p G ≈ In + en ,
I ≈ Im + G ≈ In + em ,
p+1
donde G = F (b−a)
np 2p f
(q) (ζ), e y e son los errores. Se puede despejar G:
n m
Im − In
em ≈ G = (6.10)
2p − 1
Im − In
= trapecio
3
Im − In
= Simpson
15
I ≈ Im + G. (6.11)
Ejemplo 6.5. Z π
I= sen(x)dx,
0
1 0.0000000000000002
2 1.5707963267948966 0.5235987755982988
4 1.8961188979370398 0.1084408570473811
8 1.9742316019455508 0.0260375680028370
16 1.9935703437723395 0.0064462472755962
32 1.9983933609701441 0.0016076723992682
64 1.9995983886400375 0.0004016758899645
128 1.9998996001842038 0.0001004038480554
256 1.9999749002350531 0.0000251000169498
512 1.9999937250705768 0.0000062749451746
1024 1.9999984312683834 0.0000015687326022
2048 1.9999996078171378 0.0000003921829181
4096 1.9999999019542845 0.0000000980457155
8192 1.9999999754885744 0.0000000245114300
16384 1.9999999938721373 0.0000000061278543
Método de Simpson:
n In G
2 2.0943951023931953
4 2.0045597549844207 -0.0059890231605850
8 2.0002691699483881 -0.0002860390024022
16 2.0000165910479355 -0.0000168385933635
32 2.0000010333694127 -0.0000010371785682
64 2.0000000645300013 -0.0000000645892941
128 2.0000000040322572 -0.0000000040331829
Los valores wi se llaman los pesos o ponderaciones y los xi son las abscisas.
Si se desea integrar en otro intervalo,
Z b
ϕ(ξ) dξ
a
Z b Z 1
b−a b−a
ϕ(ξ) dξ = ϕ( (t + 1) + a) dt,
a 2 −1 2
Z b n
X
b−a b−a
ϕ(ξ) dξ ≈ wi ϕ( (xi + 1) + a), (6.13)
a 2 2
i=1
Z b Xn
b−a
ϕ(ξ) dξ ≈ wi ϕ(ξi ), (6.14)
a 2
i=1
b−a
ξi = (xi + 1) + a. (6.15)
2
En la cuadratura de Gauss se desea que la fórmula (6.12) sea exacta para los
polinomios de grado menor o igual que m = mn , y se desea que este valor
mn sea lo más grande posible. En particular,
Z 1 n
X
f (x) dx = wi f (xi ) , si f (x) = 1, x, x2 , ..., xmn .
−1 i=1
Recordemos que
Z 0 si k es impar,
1
xk dx =
−1
2
si k es par.
k+1
Para n = 1, se debe cumplir
Z 1
w1 = 1 dx = 2,
−1
Z 1
w1 x1 = x dx = 0.
−1
w1 = 2 , x1 = 0. (6.16)
Para n = 2,
Z 1
w1 + w2 = 1 dx = 2,
−1
Z 1
w1 x1 + w2 x2 = x dx = 0,
−1
Z 1
2
w1 x21 + w2 x22 = x2 dx = ,
−1 3
Z 1
w1 x31 + w2 x32 = x3 dx = 0.
−1
x1 < 0 < x2 ,
x1 = −x2 ,
w1 = w2 .
228 6. INTEGRACIÓN Y DIFERENCIACIÓN
Entonces
2w1 = 2,
2
2w1 x21 = .
3
Finalmente,
r
1
w1 = 1,x1 = − ,
3
r
1
w2 = 1,x2 = .
3
Para n = 3,
w1 + w2 + w3 = 2,
w1 x1 + w2 x2 + w3 x3 = 0,
2
w1 x21 + w2 x22 + w3 x23 = ,
3
3 3 3
w1 x1 + w2 x2 + w3 x3 = 0,
2
w1 x41 + w2 x42 + w3 x43 = ,
5
5 5 5
w1 x1 + w2 x2 + w3 x3 = 0.
x1 < 0 = x2 < x3 ,
x1 = −x3 ,
w1 = w3 .
Entonces
2w1 + w2 = 2,
2
2w1 x21 = ,
3
2
2w1 x41 = .
5
6.7. CUADRATURA DE GAUSS-LEGENDRE 229
Finalmente,
r
5 3
w1 = ,x1 = − ,
9 5
8
w2 = ,x2 = 0,
9 r
5 3
w3 = ,x3 = .
9 5
La siguiente tabla contiene los valores wi , xi , para valores de n menores o
iguales a 6.
n wi xi
1 2 0
2 1 ±0.577350269189626
3 0.888888888888889 0
0.555555555555556 ±0.774596669241483
4 0.339981043584856 ±0.652145154862546
0.861136311594053 ±0.347854845137454
5 0.568888888888889 0
0.478628670499366 ±0.538469310105683
0.236926885056189 ±0.906179845938664
6 0.467913934572691 ±0.238619186083197
0.360761573048139 ±0.661209386466265
0.171324492379170 ±0.932469514203152
22n+1 (n!)4
en = f (2n) (ξ) , −1 < ξ < 1 . (6.17)
(2n + 1)((2n)!)3
Para 6.14 el error está dado por:
Hay varias maneras de definir los polinomios de Legendre; una de ellas es:
P0 (x) = 1, (6.19)
n
1 d
Pn (x) = (x2 − 1)n . (6.20)
2n n! dxn
Por ejemplo,
P0 (x) = 1,
P1 (x) = x,
1
P2 (x) = (3x2 − 1),
2
1
P3 (x) = (5x3 − x),
2
1
P4 (x) = (35x4 − 30x2 + 3).
8
P0 (x) = 1, (6.21)
P1 (x) = x, (6.22)
2n + 1 n
Pn+1 (x) = x Pn (x) − Pn−1 (x) . (6.23)
n+1 n+1
Algunas de las propiedades de los polinomios de Legendre son:
Z 1
• xk Pn (x) dx = 0 , k = 0, 1, 2, ..., n − 1, (6.24)
−1
Z 1
• Pm (x)Pn (x) dx = 0 , m 6= n, (6.25)
−1
Z 1
2
• (Pn (x))2 dx = · (6.26)
−1 2n + 1
232 6. INTEGRACIÓN Y DIFERENCIACIÓN
En general,
(yi+1 − yi ) h
f ′ (xi ) = − f ′′ (ξ), ξ ∈ [xi , xi+1 ] (6.29)
h 2
(y i − y i−1 ) h
f ′ (xi ) = + f ′′ (ζ), ζ ∈ [xi−1 , xi ] (6.30)
h 2
El primer término después del signo igual corresponde al valor aproximado.
El segundo término es el error. Se acostumbra decir que el error es del orden
de h. Esto se escribe
(yi+1 − yi )
f ′ (xi ) = + O(h),
h
(yi − yi−1 )
f ′ (xi ) = + O(h).
h
Para n = 2, sea s = (x − x0 )/h,
s(s − 1) ∆2 f0
p2 (x) = y0 + s∆f0 + ,
2 2
x − x0 x − x0 x − x0 − h ∆2 f0
p2 (x) = y0 + ∆f0 + ,
h h h 2
∆f0 2x − 2x0 − h ∆2 f0
p′2 (x) = + ,
h h2 2
∆f0 ∆2 f0
p′2 (x1 ) = + = ···
h 2h
y2 − y0
p′2 (x1 ) = ·
2h
Entonces
y2 − y0 h2 ′′′
f ′ (x1 ) = − f (ξ) , ξ ∈ [x0 , x2 ].
2h 6
De manera general,
yi+1 − yi−1 h2 ′′′
f ′ (xi ) = − f (ξ) , ξ ∈ [xi−1 , xi+1 ], (6.31)
2h 6
′ yi+1 − yi−1
f (xi ) = + O(h2 ).
2h
234 6. INTEGRACIÓN Y DIFERENCIACIÓN
En [YoG72], página 357, hay una tabla con varias fórmulas para diferen-
ciación numérica. Para la segunda derivada, una fórmula muy empleada
es:
yi+1 − 2yi + yi−1 h2 (4)
f ′′ (xi ) = − f (ξ) , ξ ∈ [xi−1 , xi+1 ], (6.32)
h2 12
yi+1 − 2yi + yi−1
f ′′ (xi ) = + O(h2 ).
h2
La deducción de las fórmulas de derivación numérica se hizo a partir de
una tabla de valores (xi , yi ), pero para el uso de éstas solamente se requiere
conocer o poder evaluar f en los puntos necesarios. Por esta razón, algunas
veces las fórmulas aparecen directamente en función de h:
f (x + h) − f (x)
f ′ (x) = + O(h), (6.33)
h
f (x) − f (x − h)
f ′ (x) = + O(h), (6.34)
h
f (x + h) − f (x − h)
f ′ (x) = + O(h2 ), (6.35)
2h
f (x + h) − 2f (x) + f (x − h)
f ′′ (x) = + O(h2 ). (6.36)
h2
√
Ejemplo 6.7. Dada f (x) = x, evaluar aproximadamente f ′ (4) y f ′′ (4),
utilizando h = 0.2.
2.0494 − 2
f ′ (4) ≈ = 0.2470
0.2
2 − 1.9494
f ′ (4) ≈ = 0.2532
0.2
2.0494 − 1.9494
f ′ (4) ≈ = 0.2501
2 × 0.2
2.0494 − 2 × 2 + 1.9494
f ′′ (4) ≈ = −0.0313 . 3
0.22
El error de las dos primeras aproximaciones no es el mismo, pero es del
mismo orden de magnitud O(h). La tercera aproximación es mejor que
las anteriores; su error es del orden de O(h2 ). Los valores exactos son
f ′ (4) = 0.25, f ′′ (4) = −0.03125.
∂f 1
(x̄) = f (x̄1 , ..., x̄i−1 , x̄i + h, x̄i+1 , ..., x̄n )
∂xi 2h
− f (x̄1 , ..., x̄i−1 , x̄i − h, x̄i+1 , ..., x̄n ) + O(h2 ) (6.37)
6.9.2 En Scilab
function y = func246(x)
y = sqrt(x)
endfunction
236 6. INTEGRACIÓN Y DIFERENCIACIÓN
der = derivative(func246, 4)
function y = func245( x )
y = exp(x(1)) * sin(x(2))
endfunction
x = [2 3]’
g = derivative(func245, x)
x = [2 3]’
[g, A] = derivative(func245, x, H_form =’blockmat’)
function fx = func247( x )
fx = zeros(3,1)
fx(1) = exp(x(1)) * sin(x(2))
fx(2) = 3*x(1) + 4*x(2)
fx(3) = x(1)*x(1) + 5*x(1)*x(2) + 3*x(2)*x(2)
endfunction
x = [2 3]’
J = derivative(func247, x)
Ejercicios
6.1 Calcule Z 1
ex dx
0.2
6.2 Calcule Z 1
2
e−x dx
0
6.7 Considere los mismos puntos del ejercicio anterior. Calcule una aprox-
imación de f ′ (0.25), f ′ (0.225), f ′′ (0.30).
Ecuaciones diferenciales
239
240 7. ECUACIONES DIFERENCIALES
7.0.3 En Scilab
x+y
y′ = + 4 + cos(x) ,
x2 + y 2
y(2) = 3.
function Dy = func158(x, y)
Dy = ( x + y )/( x*x + y*y ) + 4 + cos(x)
endfunction
x0 = 2
y0 = 3
t = 2:0.05:3;
yt = ode(y0, x0, t, func158)
plot2d(t, yt)
7.1. MÉTODO DE EULER 241
y ′ (x0 ) = f (x0 , y0 ).
Entonces
y(x0 + h) ≈ y0 + hf (x0 , y0 ).
y1 = y0 + hf (x0 , y0 ).
ϕ(t, h) = O(hp )
b
(x1 , y1 )
y(x0 + h)
y0 = y(x0 ) b
x0 x0 + h
ϕ(t, h) ≈ chp .
El error local tiene que ver con el error cometido para calcular y(xi+1 )
suponiendo que yi es un valor exacto, es decir, yi = y(xi ). El error global
es el error que hay al considerar yn como aproximación de y(xn ) (n indica
el número de intervalos).
Los resultados sobre el error en el método de Euler son:
y ′ = 2x2 − 4x + y
y(1) = 0.7182818
1
b b
0 b
b
b
b
b
b
−1
1 2 3
xi ỹ(xi ) y(xi )
1.00 0.7182818 0.7182818
1.25 0.3978523 0.3653430
1.50 0.0285654 -0.0183109
1.75 -0.3392933 -0.3703973
2.00 -0.6428666 -0.6109439
2.25 -0.8035833 -0.6372642
2.50 -0.7232291 -0.3175060
2.75 -0.2790364 0.5176319
3.00 0.6824545 2.0855369
h = (xf-x0)/n
X = zeros(1,n+1)
Y = X
X(1) = x0
Y(1) = y0
xi = x0
yi = y0
for i=1:n
yi = yi + h*f(xi,yi)
xi = xi+h
Y(i+1) = yi
X(i+1) = xi
end
endfunction
b
(x1 , y1 )
y(x0 + h)
y0 = y(x0 ) b
x0 x0 + h
y ′ (x) + y ′ (x + h)
y(x + h) ≈ y(x) + h
2
o sea,
K1 = hf (xi , yi )
K2 = hf (xi + h, yi + K1 ) (7.5)
1
yi+1 = yi + (K1 + K2 ).
2
246 7. ECUACIONES DIFERENCIALES
y ′ = 2x2 − 4x + y
y(1) = 0.7182818
K1 = hf (x0 , y0 )
= 0.25f (1, 0.7182818)
= −0.320430
K2 = hf (x0 + h, y0 + K1 )
= 0.25f (1.25, 0.397852)
= −0.369287
y1 = y0 + (K1 + K2 )/2
= 0.3734236
K1 = hf (x1 , y1 )
= 0.25f (1.25, 0.3734236)
= −0.375394
K2 = hf (x1 + h, y1 + K1 )
= 0.25f (1.500000, −0.001971)
= −0.375493
y2 = y1 + (K1 + K2 )/2
= −0.0020198
K1 = ...
7.3. MÉTODO DEL PUNTO MEDIO 247
b
2
1
b
b
b
0 b
b
b
b b
−1
1 2 3
xi ỹ(xi ) y(xi )
1.00 0.7182818 0.7182818
1.25 0.3734236 0.3653430
1.50 -0.0020198 -0.0183109
1.75 -0.3463378 -0.3703973
2.00 -0.5804641 -0.6109439
2.25 -0.6030946 -0.6372642
2.50 -0.2844337 -0.3175060
2.75 0.5418193 0.5176319
3.00 2.0887372 2.0855369
En este ejemplo los resultados son mucho mejores. Por un lado, el método
es mejor, pero, por otro, es natural tener mejores resultados pues hubo
que evaluar 16 veces la función f (x, y), 2 veces en cada iteración. En el
ejemplo del método de Euler hubo simplemente 8 evaluaciones de la función
f (x, y). Al aplicar el método de Heun con h = 0.5 (es necesario evaluar 8
veces la función) se obtiene ỹ(3) = 2.1488885, resultado no tan bueno como
2.0887372, pero netamente mejor que el obtenido por el método de Euler. Si
se trabaja con h = 0.1 se obtiene ỹ(3) = 2.0841331; con h = 0.01 se obtiene
ỹ(3) = 2.0855081; con h = 0.001 se obtiene ỹ(3) = 2.0855366. 3
y(x0 + h) b
(x1 , y1 )
y0 = y(x0 ) b
x0 x0 + h/2 x0 + h
la aproximación
y(x + h) = y(x) + hy ′ (x).
o sea,
h
y(x + h/2) ≈ y(x) + f (x, y)
2
h
y(x + h) ≈ y(x) + h f (x + h/2, y(x) + f (x, y))·
2
La anterior aproximación suele escribirse de la siguiente manera:
K1 = hf (xi , yi )
K2 = hf (xi + h/2, yi + K1 /2) (7.6)
yi+1 = yi + K2 .
7.3. MÉTODO DEL PUNTO MEDIO 249
Ejemplo 7.3. Resolver, por el método del punto medio, la ecuación difer-
encial
y ′ = 2x2 − 4x + y
y(1) = 0.7182818
K1 = hf (x0 , y0 )
= 0.25f (1, 0.7182818)
= −0.320430
K2 = hf (x0 + h/2, y0 + K1 /2)
= 0.25f (1.125, 0.558067)
= −0.352671
y1 = y0 + K2
= 0.3656111
K1 = hf (x1 , y1 )
= 0.25f (1.25, 0.3656111)
= −0.377347
K2 = hf (x1 + h/2, y1 + K1 /2)
= 0.25f (1.375, 0.176937)
= −0.385453
y2 = y1 + K2
= −0.0198420
K1 = ...
250 7. ECUACIONES DIFERENCIALES
2 b
1
b
b b
0 b
b b
b b
−1
1 2 3
Figura 7.6: Ejemplo del método del punto medio
xi ỹ(xi ) y(xi )
1.00 0.7182818 0.7182818
1.25 0.3656111 0.3653430
1.50 -0.0198420 -0.0183109
1.75 -0.3769851 -0.3703973
2.00 -0.6275434 -0.6109439
2.25 -0.6712275 -0.6372642
2.50 -0.3795415 -0.3175060
2.75 0.4121500 0.5176319
3.00 1.9147859 2.0855369
Método de Euler:
K1 = hf (xi , yi ) (7.8)
yi+1 = yi + K1 .
Método de Heun:
K1 = hf (xi , yi )
K2 = hf (xi + h, yi + K1 ) (7.9)
1 1
yi+1 = yi + K1 + K2 .
2 2
Método del punto medio:
K1 = hf (xi , yi )
1 1
K2 = hf (xi + h, yi + K1 ) (7.10)
2 2
yi+1 = yi + 0K1 + K2 .
y ′ = 2x2 − 4x + y
y(1) = 0.7182818
K1 = hf (x0 , y0 )
= 0.25f (1, 0.7182818)
= −0.320430
K2 = hf (x0 + h/2, y0 + K1 /2)
= 0.25f (1.125, 0.558067)
= −0.352671
K3 = hf (x0 + h/2, y0 + K2 /2)
= 0.25f (1.125, 0.541946)
= −0.356701
K4 = hf (x0 + h, y0 + K3 )
= 0.25f (1.25, 0.361581)
= −0.378355
y1 = y0 + (K1 + 2K2 + 2K3 + K4 )/6
= 0.3653606
K1 = hf (x1 , y1 )
= 0.25f (1.25, 0.3653606)
= −0.377410
K2 = hf (x1 + h/2, y1 + K1 /2)
= 0.25f (1.375, 0.176656)
= −0.385524
7.4. MÉTODO DE RUNGE-KUTTA 253
b
2
1
b
b
b
0 b
b
b
b b
−1
1 2 3
Figura 7.7: Ejemplo del método Runge-Kutta 4
xi ỹ(xi ) y(xi )
1.00 0.7182818 0.7182818
1.25 0.3653606 0.3653430
1.50 -0.0182773 -0.0183109
1.75 -0.3703514 -0.3703973
2.00 -0.6108932 -0.6109439
2.25 -0.6372210 -0.6372642
2.50 -0.3174905 -0.3175060
2.75 0.5175891 0.5176319
3.00 2.0853898 2.0855369
En este ejemplo, los resultados son aún mejores. Hubo que evaluar 32 veces
la función f (x, y), 4 veces en cada iteración. Si se trabaja con h = 0.1 se
obtiene ỹ(3) = 2.0855314; con h = 0.01 se obtiene ỹ(3) = 2.0855369; con
h = 0.001 se obtiene ỹ(3) = 2.0855369. 3
254 7. ECUACIONES DIFERENCIALES
h = (xf-x0)/n
X = zeros(1,n+1)
Y = X
X(1) = x0
Y(1) = y0
xi = x0
yi = y0
for i=1:n
K1 = h*f(xi, yi)
K2 = h*f(xi+h/2, yi+K1/2);
K3 = h*f(xi+h/2, yi+K2/2);
K4 = h*f(xi+h, yi+K3);
xi = xi+h
yi = yi + (K1 + 2*K2 + 2*K3 + K4)/6
Y(i+1) = yi
X(i+1) = xi
end
endfunction
7.5. DEDUCCIÓN DE RK2 255
Entonces
y ′′ = fx + f fy
h2 h2
y(xi + h) ≈ y + hf + fx + f fy . (7.15)
2 2
Por otro lado, el método RK2 se puede reescribir:
yi+1 = yi + R1 hf (xi , yi ) + R2 hf (xi + αh, yi + βK1 ).
Utilizando (7.12):
yi+1 = yi + R1 hf (xi , yi )
∂f ∂f
+ R2 h f (xi , yi ) + αh (xi , yi ) + βK1 (xi , yi ) .
∂x ∂y
Utilizando la notación se obtiene:
yi+1 = y + R1 hf + R2 h (f + αhfx + βK1 fy )
yi+1 = y + (R1 + R2 )hf + R2 h2 αfx + R2 hβK1 fy .
Como K1 = hf , entonces
yi+1 = y + (R1 + R2 )hf + R2 αh2 fx + R2 βh2 f fy . (7.16)
Al hacer la igualdad y(xi + h) = yi+1 , en las ecuaciones (7.15) y (7.16) se
comparan los coeficientes de hf , de h2 fx y de h2 f fy y se deduce:
R1 + R2 = 1,
1
R2 α = ,
2
1
R2 β = .
2
7.6. CONTROL DEL PASO 257
Entonces
β = α, (7.17)
1
R2 = · (7.18)
2α
R1 = 1 − R 2 . (7.19)
C0 = Φ2 (yA , yB , h, p)
C = ϕ(C0 , ...)
h′ = Ch.
1 1
2
4 4
3 3 9
3
8 32 32
12 1932 7200 7296
4 −
13 2197 2197 2197
(7.20)
439 3680 845
5 1 −8 −
216 513 4104
1 8 3544 1859 11
6 − 2 − −
2 27 2565 4104 40
25 1408 2197 1
yA = y i + K1 + 0K2 + K3 + K4 − K5
216 2565 4104 5
16 6656 28561 9 2
yB = y i + K1 + 0K2 + K3 + K4 − K5 + K6
135 12825 56430 50 55
Los errores locales son respectivamente O(h5 ) y O(h6 ). Realmente hay varias
fórmulas RK5 y RK6; las anteriores están en [BuF85] y [EnU96]. Hay otras
fórmulas diferentes en [ChC99].
La aproximación del error está dada por
|yA − yB |
|error| ≈ e = . (7.21)
h
El coeficiente para la modificación del valor de h está dado por:
ε 1/4
C0 = 0.84 ,
e
C = min{C0 , 4}, (7.22)
C = max{C, 0.1}.
Las fórmulas anteriores buscan que C no sea muy grande ni muy pequeño.
Más especı́ficamente, C debe estar en el intervalo [0.1, 4].
En la descripción del algoritmo usaremos la siguiente notación de Matlab y
de Scilab. La orden
u = [u; t]
260 7. ECUACIONES DIFERENCIALES
MÉTODO RUNGE-KUTTA-FEHLBERG
datos: x0 , y0 , b, h0 , ε, hmin
x = x0 , y = y0 , h = h0
X = [x0 ], Y = [y0 ]
mientras x < b
h = min{h, b − x}
hbien = 0
mientras hbien = 0
calcular ya , yB según (7.20)
e = |yA − yB |/h
si e ≤ ε
x = x + h, y = yB
bienh = 1
X = [X; x], Y = [Y ; y]
fin-si
C0 = 0.84(ε/e)1/4
C = max{C0 , 0.1}, C = min{C, 4}
h = Ch
si h < hmin ent parar
fin-mientras
fin-mientras
Ejemplo 7.5. Resolver, por el método RKF con control de paso, la ecuación
7.6. CONTROL DEL PASO 261
diferencial
y ′ = 2x2 − 4x + y
y(1) = 0.7182818
yA = −0.01834063
yB = −0.01830704
e = 0.00006717
h = 0.5 no sirve.
C0 = 0.29341805
C = 0.29341805
h = 0.14670902
yA = 0.51793321
yB = 0.51793329
e = 0.00000057
h = 0.14670902 sirve.
x = 1.14670902
y = 0.51793329
C0 = 0.96535578
C = 0.96535578
h = 0.14162640
yA = 0.30712817
yB = 0.30712821
e = 0.00000029
h = 0.14162640 sirve.
x = 1.28833543
y = 0.30712821
..
.
262 7. ECUACIONES DIFERENCIALES
b
2 b
1 b
b
b b
b
b
b
0 b
b
b b
b b
b b
b b
−1
1 2 3
Figura 7.8: Ejemplo del método Runge-Kutta-Fehlberg
x h ỹ(x) y(x)
1.0000000 0.1467090 0.7182818 0.7182818
1.1467090 0.1416264 0.5179333 0.5179333
1.2883354 0.1622270 0.3071282 0.3071282
1.4505624 0.1686867 0.0572501 0.0572501
1.6192491 0.1333497 -0.1946380 -0.1946380
1.7525988 0.1329359 -0.3736279 -0.3736279
1.8855347 0.1191306 -0.5206051 -0.5206051
2.0046653 0.1092950 -0.6137572 -0.6137571
2.1139603 0.1024064 -0.6566848 -0.6566847
2.2163666 0.0971218 -0.6506243 -0.6506241
2.3134884 0.0928111 -0.5948276 -0.5948275
2.4062996 0.0891591 -0.4877186 -0.4877184
2.4954587 0.0859853 -0.3273334 -0.3273332
2.5814440 0.0831757 -0.1114979 -0.1114977
2.6646196 0.0806534 0.1620898 0.1620900
2.7452730 0.0783639 0.4958158 0.4958160
2.8236369 0.0762674 0.8921268 0.8921270
2.8999043 0.0743333 1.3535162 1.3535164
2.9742376 0.0257624 1.8825153 1.8825156
3.0000000 2.0855366 2.0855369
7.7. ORDEN DEL MÉTODO Y ORDEN DEL ERROR 263
Para algunos de los métodos hasta ahora vistos, todos son métodos de RK,
se ha hablado del orden del método, del orden del error local y del orden
del error global.
El orden del método se refiere al número de evaluaciones de la función f en
cada iteración. Ası́ por ejemplo, el método de Euler es un método de orden
1 y el método de Heun es un método de orden 2.
El orden del error local se refiere al exponente de h en el error teórico
cometido en cada iteración. Si la fórmula es
El orden del error global es generalmente igual al orden del error local menos
una unidad. Por ejemplo, el error global en el método de Euler es O(h).
A medida que aumenta el orden del método, aumenta el orden del error, es
decir, el error disminuye. Pero al pasar de RK4 a RK5 el orden del error
no mejora. Por eso es más interesante usar el RK4 que el RK5 ya que se
hacen solamente 4 evaluaciones y se tiene un error semejante. Ya con RK6
se obtiene un error más pequeño, pero a costa de dos evaluaciones más.
264 7. ECUACIONES DIFERENCIALES
−5
b
−6 b
−7 b
−8
−2.5 −2.0 −1.5
y1 = 0.3653606
y2 = −0.0182773.
Entonces
f0 = f (x0 , y0 ) = −1.2817182
f1 = f (x1 , y1 ) = −1.5096394
f2 = −1.5182773
y3 = y2 + h(5f0 − 16f1 + 23f2 )/12
= −0.3760843
f3 = f (x3 , y3 ) = −1.2510843
y4 = −0.6267238
..
.
xi ỹ(xi ) y(xi )
1.00 0.7182818 0.7182818
1.25 0.3653606 0.3653430
1.50 -0.0182773 -0.0183109
1.75 -0.3760843 -0.3703973
2.00 -0.6267238 -0.6109439
2.25 -0.6681548 -0.6372642
2.50 -0.3706632 -0.3175060
2.75 0.4320786 0.5176319
3.00 1.9534879 2.0855369
En este caso hubo que evaluar 8 veces la función para los dos valores de RK4
y en seguida 6 evaluaciones para un total de 14 evaluaciones de la función
f. 3
268 7. ECUACIONES DIFERENCIALES
2 b
1
b
b
b
0 b
b b
b b
−1
1 2 3
Figura 7.10: Ejemplo del método de Adams-Bashforth 2
n error
1 ′′ 2
0 yi+1 = yi + hfi 2 y (ξ)h
h 5 ′′′ 3
1 yi+1 = yi + (−fi−1 + 3fi ) 12 y (ξ)h
2
h 3 (4) 4
2 yi+1 = yi + (5fi−2 − 16fi−1 + 23fi ) 8 y (ξ)h
12
h 251 (5) 5
3 yi+1 = yi + (−9fi−3 + 37fi−2 − 59fi−1 + 55fi ) 720 y (ξ)h
24
h 95 (6) 6
4 yi+1 = yi + (251fi−4 − 1274fi−3 + 2616fi−2 288 y (ξ)h
720
−2774fi−1 + 1901fi )
1 2 2
p2 (x) = (f0 − 2f1 + f2 )x + (−3f0 + 4f1 − f2 )hx + 2h f0 .
2h2
Z x2
y ′ (x)dx = y(x2 ) − y(x1 )
x1
Z x2
y(x2 ) = y(x1 ) + y ′ (x)dx
x1
Z 2h
y(x2 ) = y(x1 ) + f (x, y)dx.
h
Z 2h
y(x2 ) ≈ y(x1 ) + p2 (x)dx
h
Z 2h
1
y(x2 ) ≈ y(x1 ) + (f0 − 2f1 + f2 )x2 +
h 2h2
2
(−3f0 + 4f1 − f2 )hx + 2h f0 dx
270 7. ECUACIONES DIFERENCIALES
1 7
y2 = y1 + 2 (f0 − 2f1 + f2 ) h3 +
2h 3
3 3 3
(−3f0 + 4f1 − f2 ) h + 2h f0
2
h
y2 = y1 + (−f0 + 8f1 + 5f2 ). (7.27)
12
La anterior igualdad se conoce con el nombre de fórmula de Adams-Moul-
ton de orden 2 (se utiliza un polinomio de orden 2). También recibe el
nombre de método multipaso implı́cito o método multipaso cerrado de orden
2.
Si los valores y0 , y1 y y2 son exactos, o sea, si y0 = y(x0 ), y1 = y(x1 ) y
y2 = y(x2 ), entonces los valores fi son exactos, o sea, f (xi , yi ) = f (xi , y(xi ))
y el error está dado por
h 1 (3)
y(x2 ) = y(x1 ) + (−f0 + 8f1 + 5f2 ) − y (z)h4 , z ∈ [x0 , x2 ].
12 24
(7.28)
La fórmula (7.27) se escribe en el caso general
h
yi+1 = yi + (−fi−1 + 8fi + 5fi+1 ). (7.29)
12
Para empezar a aplicar esta fórmula es indispensable conocer los valores fj
anteriores. Entonces se requiere utilizar un método RK el número de veces
necesario. El método RK escogido debe ser de mejor calidad que el método
de Adams-Bashforth que estamos utilizando. Para nuestro caso podemos
utilizar RK4.
Una dificultad más grande, y especı́fica de los métodos implı́citos, está dada
por el siguiente hecho: para calcular yi+1 se utiliza fi+1 , pero este valor es
justamente f (xi+1 , yi+1 ). ¿Cómo salir de este cı́rculo vicioso? Inicialmente
0 , una primera aproximación, por el método de Euler. Con este
se calcula yi+1
valor se puede calcular fi+10 0 ) y en seguida y 1 . De nuevo se
= f (xi+1 , yi+1 i+1
1 = f (x
calcula fi+1 1 2
i+1 , yi+1 ) y en seguida yi+1 . Este proceso iterativo acaba
k k+1
cuando dos valores consecutivos, yi+1 y yi+1 , son muy parecidos. Este
método recibe también el nombre de método predictor-corrector. La
fórmula queda entonces ası́:
k+1 h k
yi+1 = yi + (−fi−1 + 8fi + 5fi+1 ). (7.30)
12
7.9. MÉTODOS MULTIPASO IMPLÍCITOS 271
|yik+1 − yik |
≤ ε.
max{1, |yik+1 |}
y ′ = 2x2 − 4x + y
y(1) = 0.7182818
y1 = 0.3653606
Entonces
f0 = f (x0 , y0 ) = −1.2817182
f1 = f (x1 , y1 ) = −1.5096394
y20 = −0.0120493
f20 = −1.5120493
y21 = −0.0170487
f21 = −1.5170487
y22 = −0.0175694
f22 = −1.5175694
y23 = −0.0176237 = y2
f1 = −1.5096394
f2 = −1.5176237.
272 7. ECUACIONES DIFERENCIALES
y30 = −0.3970296
f30 = −1.2720296
y31 = −0.3716132
f31 = −1.2466132
y32 = −0.3689657
f32 = −1.2439657
y33 = −0.3686899
f33 = −1.2436899
y34 = −0.3686612 = y3
..
.
xi ỹ(xi ) y(xi )
1.00 0.7182818 0.7182818
1.25 0.3653606 0.3653430
1.50 -0.0176237 -0.0183109
1.75 -0.3686612 -0.3703973
2.00 -0.6076225 -0.6109439
2.25 -0.6315876 -0.6372642
2.50 -0.3084043 -0.3175060
2.75 0.5316463 0.5176319
3.00 2.1065205 2.0855369
En este caso hubo que evaluar 4 veces la función para el valor de RK4 y en
seguida, en cada uno de los otros 7 intervalos, una evaluación fija más las
requeridas al iterar. En este ejemplo hubo, en promedio, 4 por intervalo,
para un total de 32 evaluaciones de f . El valor final y8 es más exacto que
el obtenido por Adams-Bashforth, pero a costa de más evaluaciones. 3
Teóricamente, los dos métodos multipaso de orden 2 tienen un error local del
mismo orden, O(h4 ), pero el coeficiente en el método multipaso explı́cito,
3/8, es nueve veces el coeficiente en el error del método implı́cito, 1/24.
7.9. MÉTODOS MULTIPASO IMPLÍCITOS 273
1
b
b
b
0 b
b
b
b b
−1
1 2 3
Figura 7.11: Ejemplo del método de Adams-Moulton 2
n error
h 1 ′′
1 yi+1 = yi + (fi + fi+1 ) − 12 y (ξ)h3
2
h 1 (3)
2 yi+1 = yi + (−fi−1 + 8fi + 5fi+1 ) − 24 y (ξ)h4
12
h 19 (4)
3 yi+1 = yi + (fi−2 − 5fi−1 + 19fi + 9fi+1 ) − 720 y (ξ)h5
24
h 27 (5)
4 yi+1 = yi + (−19fi−3 + 106fi−2 − 264fi−1 − 1440 y (ξ)h6
720
+646fi + 251fi+1 )
dy1
= f1 (x, y1 , y2 , ..., ym )
dx
dy2
= f2 (x, y1 , y2 , ..., ym )
dx
..
.
dym
= fm (x, y1 , y2 , ..., ym )
dx
y1 (x0 ) = y10
y2 (x0 ) = y20
..
.
0
ym (x0 ) = ym .
y = (y1 , y2 , ..., ym )
y 0 = (y10 , y20 , ..., ym
0
)
f (x, y) = f (x, y1 , y2 , ..., ym )
= ( f1 (x, y1 , ..., ym ), f2 (x, y1 , ..., ym ), ..., fm (x, y1 , ..., ym ) ).
y ′ = f (x, y), x0 ≤ x ≤ b
y(x0 ) = y 0 .
yji ≈ yj (xk ).
K 1 = hf (xi , y i )
K 2 = hf (xi + h/2, y i + K 1 /2)
K 3 = hf (xi + h/2, y i + K 2 /2) (7.31)
4 i 3
K = hf (xi + h, y + K )
y i+1 = y i + (K 1 + 2K 2 + 2K 3 + K 4 )/6.
2y1
y1′ = + x3 y2 , 1 ≤ x ≤ 2
x
3
y2′ = − y2
x
y1 (1) = −1
y2 (1) = 1
con h = 0.2.
La solución (exacta) de este sencillo sistema de ecuaciones es:
y1 (x) = −x
y2 (x) = x−3 .
276 7. ECUACIONES DIFERENCIALES
K 1 = (−0.2001062, −0.2895317)
K 2 = (−0.2093988, −0.2004450)
K 3 = (−0.1912561, −0.2210035)
K 4 = (−0.2011961, −0.1534542)
y 2 = (−1.4011269, 0.3647495)
..
.
7.10.1 En Scilab
2y1
y1′ = + x3 y2 ,
x
3
y2′ = − y2
x
y1 (1) = −1
y2 (1) = 1.
x0 = 1
y0 = [-1 1]’
t = (1:0.2:2)’
yt = ode(y0, x0, t, func43)
// n = numero de subintervalos.
//
// Y sera una matriz con p filas, n+1 columnas.
// X sera un vector fila de n+1 elementos.
// Cada columna de Y contendra las aproximaciones de
278 7. ECUACIONES DIFERENCIALES
h = (xf-x0)/n
p = size(y0,1)
disp(p, ’p’)
X = zeros(1,n+1)
Y = zeros(p,n+1)
X(1) = x0
Y(:,1) = y0
xi = x0
yi = y0
for i=1:n
K1 = h*f(xi, yi)
K2 = h*f(xi+h/2, yi+K1/2);
K3 = h*f(xi+h/2, yi+K2/2);
K4 = h*f(xi+h, yi+K3);
xi = xi+h
yi = yi + (K1 + 2*K2 + 2*K3 + K4)/6
Y(:,i+1) = yi
X(i+1) = xi
end
endfunction
(m−1)
donde κ0 = [y0 y0′ y0′′ ... y0 ]T . Este sistema se puede resolver por
los métodos para sistemas de ecuaciones diferenciales de primer orden.
4y − xy ′
y ′′ = , 1 ≤ x ≤ 2,
x2
y(1) = 3
y ′ (1) = 10,
u′1 = u2
4u1 − xu2
u′2 = , 1 ≤ x ≤ 2,
x2
u1 (1) = 3
u2 (1) = 10.
K 1 = (2, 0.4)
K 2 = (2.04, 0.7900826)
K 3 = (2.0790083, 0.7678437)
K 4 = (2.1535687, 1.0270306)
u1 = (5.0652642, 10.7571472)
..
.
y ′′ = f (x, y, y ′ ), a ≤ x ≤ b,
y(a) = ya (7.32)
y(b) = yb .
yb − ya
ya′ ≈ .
b−a
donde:
v = aproximación de ya′ ,
n = número de intervalos para la solución numérica,
ỹ = (ỹ0 , ỹ1 , ..., ỹn ) = solución numérica del siguiente problema:
y ′ = f (x, y), a ≤ x ≤ b
y(a) = ya P(v)
′
y (a) = v,
2 cos(2x) − y ′ − 4x2 y
y ′′ = , 0.2 ≤ x ≤ 0.7
x2
y(0.2) = 0.3894183
y(0.7) = 0.9854497,
ỹ5 = 0.94935663.
ϕ0 = 0.03609310
v1 = 1.19206278 + 0.03609310/(0.7 − 0.5) = 1.26424897
ỹ5 = 0.95337713
ϕ1 = 0.03207260
v2 = 1.84009748
ỹ5 = 0.98544973
Este disparo fue preciso (no siempre se obtiene la solución con una sola
iteración de la secante). El último vector ỹ es la solución. La solución
exacta es y = sen(2x).
284 7. ECUACIONES DIFERENCIALES
xi ỹ(xi ) y(xi )
0.2 0.3894183 0.3894183
0.3 0.5647741 0.5646425
0.4 0.7174439 0.7173561 3
0.5 0.8415217 0.8414710
0.6 0.9320614 0.9320391
0.7 0.9854497 0.9854497
Entonces:
yi−1 − 2yi + yi+1 −yi−1 + yi+1
pi + qi + ri yi = si , i = 1, ..., n − 1.
h2 2h
7.13. ECUACIONES LINEALES CON CONDICIONES DE FRONTERA285
d1 u1 y1 β1
l1 d2 u2 y2 β2
l2 d3 u3 y3 β3
= , (7.35)
ln−3 dn−2 un−2 yn−2 βn−2
ln−2 dn−1 yn−1 βn−1
donde
d1 = −4p1 + 2h2 r1
d1 = −4(0.3)2 + 2(0.1)2 4(0.3)2 = −0.3528
u1 = 2p1 + hq1
u1 = 2(0.3)2 + 0.1(1) = 0.28
l1 = 2p2 − hq2
l1 = 2(0.4)2 − 0.1(1) = 0.22
d = (−0.3528, −0.6272, −0.98, −1.4112),
u = (0.28, 0.42, 0.6),
l = (0.22, 0.4, 0.62),
β = (0.00186, 0.0278683, 0.0216121, −0.7935745).
Su solución es
xi ỹ(xi ) y(xi )
0.2 0.3894183 0.3894183
0.3 0.5628333 0.5646425
0.4 0.7158127 0.7173561 3
0.5 0.8404825 0.8414710
0.6 0.9315998 0.9320391
0.7 0.9854497 0.9854497
Ejercicios
7.1
y
y ′ = ex −
x
y(1) = 0.
ex
Su solución es y = ex − ·
x
7.2
y1′ = 2y1 + y2 + 3
y2′ = 4y1 − y2 + 9
y1 (0) = −3
y2 (0) = 5.
7.3
2
y ′′ = y′
x(2 − x)
y(1) = −2
y ′ (1) = 1.
7.4
y ′′′ + y ′′ + y ′ + y = 4ex
y(0) = 1
y ′ (0) = 2
y ′′ (0) = 1
y ′′′ (0) = 0.
Su solución es y = ex + sen(x).
288 7. ECUACIONES DIFERENCIALES
7.5
Su solución es y = ex + sen(x).
7.6
Su solución es y = ex + sen(x).
8
Ecuaciones diferenciales
parciales
8.1 Generalidades
Sea u = u(x, y) una función de dos variables con derivadas parciales de orden
dos. Una ecuación diferencial se llama cuasi-lineal si es de la forma
uxx + uyy = 0.
289
290 8. ECUACIONES DIFERENCIALES PARCIALES
utt = c2 uxx .
△u(x, y) = f (x, y) en Ω,
(8.1)
u(x, y) = g(x, y) en ∂Ω.
nx ∈ Z, nx ≥ 1,
ny ∈ Z, ny ≥ 1,
b−a
hx = ,
nx + 1
d−c
hy = ,
ny + 1
xi = a + ihx , i = 1, ..., nx ,
yj = c + jhy , j = 1, ..., ny ,
uij ≈ u(xi , yj ), i = 1, ...nx , j = 1, ..., ny .
8.2. ELÍPTICAS: ECUACIÓN DE POISSON 291
d
ui,j+1
yny b
hy
b
ui,j−1
y1
c
hx
a x1 x2 xi xnx b
Usando la aproximación
ϕ(t + h) − 2ϕ(t) + ϕ(t − h)
ϕ′′ (t) ≈
h2
se obtiene
ui+1,j − 2uij + ui−1,j ui,j+1 − 2uij + ui,j−1
△u(xi , yj ) ≈ 2
+ . (8.2)
hx h2y
Sea η = hx /hy .
ui+1,j − 2uij + ui−1,j ui,j+1 − 2uij + ui,j−1
△u(xi , yj ) ≈ + η2
h2x h2x
ui+1,j + ui−1,j + η 2 ui,j+1 + η 2 ui,j−1 − (2 + 2η 2 )uij
△u(xi , yj ) ≈ . (8.3)
h2x
En el caso particular cuando h = hx = hy
ui+1,j + ui−1,j + ui,j+1 ui,j−1 − 4uij
△u(xi , yj ) ≈ . (8.4)
h2
Al aplicar la aproximación (8.3) en (8.1), y cambiando el signo aproximación
por el signo de igualdad, se obtiene
−ui+1,j − ui−1,j − η 2 ui,j+1 − η 2 ui,j−1 + (2 + 2η 2 )uij = −h2x fij , (8.5)
292 8. ECUACIONES DIFERENCIALES PARCIALES
n = nx
m = ny
N = nm
h = hx
h
η=
hy
ρ = η2
σ = 2 + 2η 2
αj = g(a, yj )
βj = g(b, yj )
γi = g(xi , c)
δi = g(xi , d)
Entonces
..., en (xn , y1 ), en (x1 , y2 ), ... Para las variables utilizaremos el mismo orden
ξ1 = u11
ξ2 = u21
..
.
ξn = un1
ξn+1 = u12
ξn+2 = u22
..
.
ξ2n = un2
..
.
ξN = unm
Con el anterior orden para las variables la igualdad (8.6) se reescribe ası́:
Aξ = v. (8.7)
Es necesario cambiar u10 por el valor conocido γ1 y cambiar u01 por el valor
conocido α1 . Utilizando la notación ξk se obtiene:
σ −1 0 −ρ 0 0 0 0 0 0 0 0
−1 σ −1 0 −ρ 0 0 0 0 0 0 0
0 −1 σ 0 0 −ρ 0 0 0 0 0 0
−ρ 0 0 σ −1 0 −ρ 0 0 0 0 0
0 −ρ 0 −1 σ −1 0 −ρ 0 0 0 0
0 0 −ρ 0 −1 σ 0 0 −ρ 0 0 0
A=
0 0 0 −ρ 0 0 σ −1 0 −ρ 0 0
0 0 0 0 −ρ 0 −1 σ −1 0 −ρ 0
0 0 0 0 0 −ρ 0 −1 σ 0 0 −ρ
0 0 0 0 0 0 −ρ 0 0 σ −1 0
0 0 0 0 0 0 0 −ρ 0 −1 σ −1
0 0 0 0 0 0 0 0 −ρ 0 −1 σ
En algunas filas
n
X
|aij | < aii .
j=1
j6=i
con nx = 3 y ny = 4.
296 8. ECUACIONES DIFERENCIALES PARCIALES
Entonces hx = 3, hy = 1, ρ = 9, σ = 20,
v = [235 2529 10531 −519 −810 1353 −505 −918 1367 6319 8235 16615]T .
u = [118 397 1054 192 471 1128 314 593 1250 496 775 1432]T .
xi = i hx , i = 0, 1, 2, ..., m
tj = j ht , j = 0, 1, 2, ...
donde
L
hx = ·
m
8.3. PARABÓLICAS: ECUACIÓN DEL CALOR 297
b b b b
t3 b b b b
u13
t2 b b b b
u12
t1 b b b b
0
0 x1 x2 xm−1 L
0≤t≤T (8.12)
T
ht = (8.13)
n
tj = j ht , j = 0, 1, 2, ..., n. (8.14)
b b b
Para calcular los valores u11 , u21 , ..., um−1,1 se necesitan valores uk0 , pero
estos son conocidos por corresponder a la condición (8.11). Entonces los
valores u11 , u21 , ..., um−1,1 , se calculan sin ningún problema.
Para calcular los valores u12 , u22 , ..., um−1,2 se necesitan los valores u01 , u11 ,
u21 , ..., um−1,1 , um1 . El primero y el último están dados por las condiciones
(8.9) y (8.10); los otros se acaban de calcular. Después, de manera semejante,
se calculan los valores ui3 y ası́ sucesivamente.
Ejemplo 8.2. Aplicar las fórmulas anteriores a la ecuación diferencial
∂u 2 ∂2u π
(x, t) = (x, t), 0<x< , 0<t≤2
∂t 9 ∂x2 3
con las condiciones
u(0, t) = 5, t ≥ 0
π
u( , t) = 5, t ≥ 0
3
π
u(x, 0) = sen(3x) + 5, 0≤x≤ .
3
con m = 10 y n = 50.
La solución exacta de la ecuación diferencial es
u(x, t) = e−2t sen(3x) + 5 .
Al aplicar las fórmulas se obtiene:
hx = π/30 = 0.1047198
ht = 0.04
α = 0.8105695
β = −0.6211389
Para empezar, los valores u00 , u10 , u20 , ..., u10,0 son: 5, 5.309017, 5.5877853,
5.809017, 5.9510565, 6, 5.9510565, 5.809017, 5.5877853, 5.309017, 5.
Para t = t1 = 0.04 son datos: u01 = 5, u10,0 = 5.
u11 = αu00 + βu10 + αu20
u11 = 5.2844983
u21 = αu10 + βu20 + αu30
u21 = 5.5411479
u91 = αu80 + βu90 + αu10,0
u91 = 5.2844983
300 8. ECUACIONES DIFERENCIALES PARCIALES
Para t = t2 = 0.08
5.0000000 5.0000000
5.0056598 9.2699374
5.0107657 - 3.1030622
5.0148177 16.178842
5.0174192 - 8.1110307
5.0183156 18.817809
5.0174192 - 8.1110307
5.0148177 16.178842
5.0107657 - 3.1030622
5.0056598 9.2699374
5.0000000 5.0000000
Se observa que los resultados obtenidos por las fórmulas no son buenos. Pero
si se utiliza n = 100 sı́ lo son. En la tabla siguiente están los valores teóricos,
los obtenidos con n = 100 y los errores:
c2 ht
= 0.8105695 si n = 50
h2x
c2 ht
= 0.4052847 si n = 100
h2x
∂u ui,j − ui,j−1
(xi , tj ) ≈ (8.21)
∂t ht
2
∂ u ui+1,j − 2uij + ui−1,j
2
(xi , tj ) ≈ (8.22)
∂x h2x
ui,j−1
γ −α 0 0 ... 0 u1j u1,j−1 + αu0j
−α γ −α 0 ... 0 u2j u2,j−1
0 −α γ −α ... 0
u3j = u3,j−1 (8.26)
0 0 −α γ um−1,j um−1,j−1 + αumj
Este sistema tridiagonal se puede resolver por cualquier método, pero es más
eficiente resolverlo por un método especı́fico para sistemas tridiagonales, ver
sección (2.13). Además, como la matriz del sistema es la misma para todas
las iteraciones, entonces la factorización LU tridiagonal es la misma para
todas las iteraciones y se calcula únicamente una vez. Ası́, en cada iteración
se resuelve el sistema conociendo ya la factorización LU .
Los valores u0j y umj están dados por los valores v(tj ) y w(tj ) provenientes
de las condiciones de frontera.
Ejemplo 8.3. Aplicar este método a la misma ecuación diferencial
∂u 2 ∂2u π
(x, t) = (x, t), 0<x< , 0<t≤2
∂t 9 ∂x2 3
con las condiciones
u(0, t) = 5, t ≥ 0
π
u( , t) = 5, t ≥ 0
3
π
u(x, 0) = sen(3x) + 5, 0≤x≤ .
3
8.3. PARABÓLICAS: ECUACIÓN DEL CALOR 303
con m = 10 y n = 50.
Al aplicar las fórmulas se obtiene:
hx = π/30 = 0.1047198
ht = 0.04
α = 0.8105695
γ = 2.6211389
Para empezar, los valores u00 , u10 , u20 , ..., u10,0 son: 5, 5.309017, 5.5877853,
5.809017, 5.9510565, 6, 5.9510565, 5.809017, 5.5877853, 5.309017, 5.
Los valores α y γ definen la matriz del sistema (8.26) para todas las ite-
raciones. Para t = t1 = 0.04, los términos independientes son: 9.3618643,
5.5877853, 5.809017, 5.9510565, 6, 5.9510565, 5.809017, 5.5877853, 9.3618643 .
La solución del sistema es : 5.2863007, 5.5445763, 5.749545, 5.8811429,
5.9264885, 5.8811429, 5.749545, 5.5445763, 5.2863007 . Estos valores cor-
responden a u11 , u21 , ..., u91 .
La siguiente tabla muestra, para t = 2, los valores teóricos, los valores
obtenidos por el método y las diferencias:
Esta fórmula tiene un error de orden O(h2t +h2x ). Ahora agruparemos dejando
a izquierda los valores buscados, t = tj , y a derecha los que se suponen
conocidos, t = tj−1 :
b b b
hx = 0.1047198
ht = 0.04
α = 0.8105695
ρ = 0.405284
µ = 1.8105695
ϕ = 0.1894305
Para empezar, los valores u00 , u10 , u20 , ..., u10,0 son: 5, 5.309017, 5.5877853,
5.809017, 5.9510565, 6, 5.9510565, 5.809017, 5.5877853, 5.309017, 5.
Los valores µ y ρ definen la matriz del sistema (8.32) para todas las ite-
raciones. Para t = t1 = 0.04, los términos independientes son: 7.3231813,
5.5644666, 5.7769216, 5.9133261, 5.9603279, 5.9133261, 5.7769216,
5.5644666, 7.3231813 .
La solución del sistema es : 5.2854339, 5.5429275, 5.7472756, 5.8784752,
5.9236835, 5.8784752, 5.7472756, 5.5429275, 5.2854339.
Estos valores corresponden a u11 , u21 , ..., u91 .
La siguiente tabla muestra, para t = 2, los valores teóricos, los valores
obtenidos por el método y las diferencias:
∂2u 2
2∂ u
(x, t) = c (x, t), 0 < x < L, 0 < t , (8.33)
∂t2 ∂x2
u(0, t) = a, t ≥ 0, (8.34)
u(L, t) = b, t ≥ 0, (8.35)
u(x, 0) = f (x), 0 ≤ x ≤ L, (8.36)
∂u
(x, 0) = g(x), 0 ≤ x ≤ L. (8.37)
∂t
∂2u ∂2u
2
(xi , tj ) = c2 2 (xi , tj ) (8.38)
∂t ∂x
8.4. HIPERBÓLICAS: ECUACIÓN DE ONDA 309
La molécula es:
bc ui,j+1
ui−1,j ui+1,j
b b b
uij
b
ui,j−1
La fórmula (8.39), con error de orden O(h2x +h2t ), se puede aplicar fácilmente,
salvo para la obtención de los valores ui1 ya que serı́a necesario conocer los
valores ui,−1 . Aproximaciones de estos valores se obtienen utilizando las
condiciones (8.37) sobre la velocidad inicial, mediante la siguiente aproxi-
mación cuyo error es del orden de O(h2t ),
ui1 − ui,−1
gi = g(xi ) ≈ ,
2ht
ui,−1 = ui1 − 2ht gi , (8.42)
∂2u ∂2u
(x, t) = 3 (x, t), 0 < x < 2, 0 < t ,
∂t2 ∂x2
con las condiciones
u(0, t) = 0, t ≥ 0,
u(2, t) = 0, t ≥ 0,
1
u(x, 0) = x(2 − x), 0 ≤ x ≤ 2,
4
∂u
(x, 0) = 0, 0 ≤ x ≤ 2,
∂t
hx = 0.2
ht = 0.04
α = 0.12
β = 1.76
Los valores inciales ui0 son: 0, 0.09, 0.16, 0.21, 0.24, 0.25, 0.24, 0.21, 0.16,
0.09, 0.
Para calcular los valores ui1 se utiliza (8.43) y se obtiene: 0, 0.0888, 0.1588,
0.2088, 0.2388, 0.2488, 0.2388, 0.2088, 0.1588, 0.0888, 0.
Los demás valores uik se calculan usando (8.39). Ası́ los valores ui2 son: 0,
0.085344, 0.1552, 0.2052, 0.2352, 0.2452, 0.2352, 0.2052, 0.1552, 0.085344,
0.
8.4. HIPERBÓLICAS: ECUACIÓN DE ONDA 311
A continuación aparecen: los valores xi , los valores u(xi , 12) exactos y los
valores aproximados obtenidos por el método:
0.00 0.000000 0.
0.20 0.021539 -3366402376293274.
0.40 0.043078 6403277832890416.
0.60 0.064617 -8813355840944206.
0.80 0.086097 10360721173025462.
1.00 0.096097 -10893906929301864.
1.20 0.086097 10360721173025462.
1.40 0.064617 -8813355840944206.
1.60 0.043078 6403277832890416.
1.80 0.021539 -3366402376293274.
2.00 0.000000 0.
γ −α u1,j+1
−α γ −α
u2,j+1
0 −α γ −α ..
=
.
um−2,j+1
0 0 0 −α γ um−1,j+1
4u1j − γu1,j−1 + α(a + u2,j−1 )
4u2j − γu2,j−1 + α(u1,j−1 + u3,j−1 )
..
. (8.48)
4um−2,j − γum−2,j−1 + α(um−3,j−1 + um−1,j−1 )
4um−1,j − γum−1,j−1 + α(um−2,j−1 + b)
8.4. HIPERBÓLICAS: ECUACIÓN DE ONDA 313
Ejercicios
∂2u ∂2u
8.1 Considere la ecuación diferencial (x, t) = (x, t), para 0 < x < 1
∂t2 ∂x2
y 0 < t con las condiciones u(0, t) = u(1, t) = 0, u(x, 0) = sen(πx),
∂u
(x, 0) = 0. La solución exacta es u(x, t) = sen(πx) cos(πt). Obtenga
∂t
la solución aproximada para t = 3 con m = 10, n = 60 y con n = 20.
Compare con la solución exacta.
9
Valores propios
9.1 Preliminares
314
9.1. PRELIMINARES 315
Av = λv
Av − λv = 0
Av − λIv = 0
(A − λI)v = 0
Como v 6= 0 es solución de un sistema homogéneo, entonces A no puede ser
invertible, es decir
det(A − λI) = 0. (9.1)
Como A es real, se puede demostrar que p(λ) = det(A − λI) es un polinomio
real de grado n:
p(λ) = pA (λ) = α0 + α1 λ + α2 λ2 + · · · + αn λn . (9.2)
Este polinomio se llama el polinomio caracterı́stico de A (algunas veces se
considera que el polinomio caracterı́stico es det(λI − A)). Entonces, para
matrices pequeñas, se pueden calcular los valores propios, obteniendo las
raı́ces del polinomio caracterı́stico de A.
Ejemplo 9.2.
1 −2
A=
3 4
1 − λ −2
det(A − λI) = det
3 4−λ
= (1 − λ)(4 − λ) + 8
= λ2 − 5λ + 10
√
λ1 = 2.5 + i 15/2
√
λ2 = 2.5 − i 15/2
A = C −1 BC.
b. A es ortogonal si A−1 = AT .
A∗ = (Ā)T = AT .
d. A es hermitiana (o hermı́tica) si A = A∗ .
e. A es unitaria si A−1 = A∗ .
D = B −1 AB.
Resultados:
1.
αn = (−1)n
αn−1 = (−1)n−1 traza(A)
α0 = det(A).
n
X
2. λi = traza(A)
i=1
n
Y
3. λi = det(A)
i=1
15. Sea A simétrica. La matriz es definida positiva sssi los valores propios
son positivos.
|λ| ≤ ||A||.
9.1.1 En Scilab
Los valores propios se calculan por medio de la función spec (en Matlab
se usa eig). Si se ha definido una matriz cuadrada a , entonces la orden
spec(a)
[V, L] = spec(a)
produce una matriz L diagonal, cuyos elementos diagonales son los valores
propios y una matriz V cuyas columnas son vectores propios normalizados
asociados correspondientes.
Sea {v 1 , v 2 , ..., v n } una base formada por vectores propios asociados a los
valores propios λ1 , λ2 , ..., λn respectivamente. Entonces x0 6= 0 se puede
9.2. MÉTODO DE LA POTENCIA 319
x0 = α1 v 1 + α2 v 2 + ... + αn v n
x1 = Ax0
x1 = A(α1 v 1 + α2 v 2 + ... + αn v n )
x1 = α1 Av 1 + α2 Av 2 + ... + αn Av n
x1 = α1 λ1 v 1 + α2 λ2 v 2 + ... + αn λn v n
x2 = Ax1
= A(α1 λ1 v 1 + α2 λ2 v 2 + ... + αn λn v n )
x2 = α1 λ1 Av 1 + α2 λ2 Av 2 + ... + αn λn Av n
x2 = α1 λ21 v 1 + α2 λ22 v 2 + ... + αn λ2n v n
..
.
xk = α1 λk1 v 1 + α2 λk2 v 2 + ... + αn λkn v n
X n k !
α i λi
xk = α1 λk1 v 1 + vi
α1 λ1
i=2
xk ≈ α1 λk1 v 1 .
De manera análoga
xk+1 ≈ α1 λk+1 1
1 v .
Entonces
xk+1 ≈ λ1 xk .
xk+1
j
≈ λ1 .
xkj
320 9. VALORES PROPIOS
Ejemplo 9.4. Usar las fórmulas anteriores, partiendo de x0 = (1, 1, 1), para
hallar el valor propio dominante de
−1 −2 −3
A = −4 −5 −6 .
−7 −8 −8
Algoritmo de la potencia
para k = 1, ...,maxit
z k = Axk−1
zk
xk = k
||z ||2
T
λk1 = xk z k
si |λk1 − λ1k−1 | ≤ ε, parar
fin-para
Este método se puede aplicar para hallar λn , el valor propio menos domi-
nante de una matriz diagonalizable e invertible A, cuando éste existe, o sea,
si
|λ1 | ≥ |λ2 | ≥ |λ3 | ≥ · · · > |λn | > 0.
Si A es invertible y tiene valores propios λ1 , λ2 ,..., λn , entonces los valores
propios de A−1 son
1 1 1
, , ..., .
λ1 λ2 λn
El valor propio dominante de A−1 es justamente 1/λn . Entonces se puede
aplicar el método de la potencia a A−1 . En lugar de escribir explı́citamente
9.3. MÉTODO DE LA POTENCIA INVERSA 323
Potencia inversa
para k = 1, ...,maxit
resolver Az k = xk−1
zk
xk = k
||z ||2
T
σ1k = xk z k
si |σ1k − σ1k−1 | ≤ ε, parar
fin-para
9.4 Factorización QR
• A = QR.
• Q ∈ Rm×m es ortogonal.
• R ∈ Rm×n es triangular superior ( rij = 0 si i > j ).
H = H(v, β) = In − βv v T .
Hv x ∈ < e1 > .
Entonces,
Hv x = αe1 .
v T p = 0 , o sea, p ∈ U .
Si
z =x−p
pT z = 0.
(x − y)T ξ = 0 .
Ejemplo 9.7.
−2 3 −5 −2/3 −1/3 2/3
x = −1 , y = 0 , v = −1 , H = −1/3 14/15 2/15
2 0 2 2/3 2/15 11/15
O también,
−2 −3 1 2/3 1/3 −2/3
x = −1 , y = 0 , v = −1 , H = 1/3 2/3 2/3
2 0 2 −2/3 2/3 −1/3
9.4. FACTORIZACIÓN QR 327
Es usual escoger v tal que v1 = 1, ası́ sólo se requiere almacenar los valores v2 ,
v3 , ..., vn . Generalmente estos valores se pueden almacenar donde estaban
x2 , x3 , xn . Además, no es necesario construir explı́citamente H, basta con
conocer v y β.
Denotaremos por H(x) la matriz que proyecta x sobre el subespacio < e1 >.
La siguiente función, ver [Par80] y [GoVa96], presenta una manera eficiente
de calcular v y β a partir de un vector columna x. Está escrita en seu-
docódigo utilizando parcialmente notación de Scilab.
[v, β] = vHouse(x)
n = dim (x)
t = x(2 : n)T x(2 : n)
v = [ 1 ; x(2 : n) ]
si t = 0
β=0
sino p
ν = x21 + t
si x1 ≤ 0
v1 = x1 − ν
sino
v1 = −t/(x1 + ν)
fin-si
β = 2v12 /(t + v12 )
v = v/v1
fin-si
fin vHouse
En resumen, dado x =∈ Rn ,
[v, β] = vHouse(x)
H(x) = H(v, β) = I − βvv T .
Ejemplo 9.8.
−2 1
5
x = −1 , v = 1/5 , β= .
3
2 −2/5
328 9. VALORES PROPIOS
c = cos(θ)
s = sen(θ),
cxi − sxk
si j = i,
yj = sxi + cxk si j = k,
xj en los demás casos.
xi
c= q ,
2 2
xi + xk
−xk
s= q .
x2i + x2k
[c, s] = csGivens(a, b)
si b = 0
c=1
s=0
sino
si |b| > |a|
t = −a/b√
s = 1/ 1 + t2
c = st
sino
t = −b/a√
c = 1/ 1 + t2
s = ct
fin-si
fin-si
fin csGivens
a 2 c = 0.5547002
= por medio de la función se obtiene
b −3 s = 0.8320503
y ası́
T
c s 2 3.6055513
= .
−s c −3 0
H
" # si p = n
b = H(n,
H b H) = In−p 0
si p < n
0 H
330 9. VALORES PROPIOS
En general,
Como se supone que en la iteración k, las columnas 1, ..., k−1 de A son nulas
debajo de la diagonal, entonces no es necesario recalcularlas. La presentación
9.4. FACTORIZACIÓN QR 331
k = 1
beta = 0.717157
v : 1 -0.98598563 -0.39439425 0.19719713 0.78878851
H =
0.2828427 0.7071068 0.2828427 -0.1414214 -0.5656854
0.7071068 0.3028029 -0.2788789 0.1394394 0.5577577
0.2828427 -0.2788789 0.8884485 0.0557758 0.2231031
-0.1414214 0.1394394 0.0557758 0.9721121 -0.1115515
-0.5656854 0.5577577 0.2231031 -0.1115515 0.5537938
A =
7.0710678 7.0710678 5.939697
0 -0.0140144 1.0874867
0 -0.6056057 -0.7650053
0 -1.1971971 -2.6174973
0 -1.7887885 -2.4699893
Q =
0.2828427 0.7071068 0.2828427 -0.1414214 -0.5656854
0.7071068 0.3028029 -0.2788789 0.1394394 0.5577577
0.2828427 -0.2788789 0.8884485 0.0557758 0.2231031
-0.1414214 0.1394394 0.0557758 0.9721121 -0.1115515
-0.5656854 0.5577577 0.2231031 -0.1115515 0.5537938
-----------------------------------------------------------------
k = 2
beta = 1.006267
v : 1 0.26914826 0.53206814 0.79498802
H =
9.4. FACTORIZACIÓN QR 333
A =
7.0710678 7.0710678 5.939697
0 2.236068 3.5777088
0 0 -0.0947664
0 0 -1.2925295
0 0 -0.4902926
Q=
0.2828427 0.4472136 0.2128929 -0.2797022 -0.7722974
0.7071068 -0.4472136 -0.4807445 -0.2596204 -0.0384964
0.2828427 -0.4472136 0.8431415 -0.0337898 0.0892790
-0.1414214 -0.4472136 -0.1021209 0.6599727 -0.5779337
-0.5656854 -0.4472136 -0.0473832 -0.6462647 -0.2451463
-----------------------------------------------------------------
k=3
beta = 1.068392
v : 1 0.87309062 0.33118770
H =
-0.0683918 -0.9328028 -0.3538382
-0.9328028 0.1855786 -0.3089328
-0.3538382 -0.3089328 0.8828131
A =
7.0710678 7.0710678 5.939697
0 2.236068 3.5777088
0 0 1.3856406
0 0 0
0 0 0
Q =
0.2828427 0.4472136 0.5196152 -0.0119059 -0.6707147
0.7071068 -0.4472136 0.2886751 0.4121526 0.2163259
0.2828427 -0.4472136 -0.0577350 -0.8203366 -0.2090802
-0.1414214 -0.4472136 -0.4041452 0.3962781 -0.6779604
334 9. VALORES PROPIOS
Observaciones:
[Q, R] = QR_Givens(A)
[m, n] = tamaño(A)
Q = Im
para k = 1 : min(m, n)
para i = m : −1 : k + 1
[c, s] = csGivens(ai−1,k , aik )
G = G(i − 1, i, c, s, m)
A = GT A
Q = QG
fin-para
fin-para
R=A
fin QR_Givens
9.4. FACTORIZACIÓN QR 335
k = 1
i = 5
c = -0.242536 s = 0.970143
G =
1 0 0 0 0
0 1 0 0 0
0 0 1 0 0
0 0 0 -0.2425356 0.9701425
0 0 0 -0.9701425 -0.2425356
A =
2 3 4
5 4 3
2 1 0
4.1231056 5.3357838 4.6081769
0 -0.7276069 -1.940285
Q =
1 0 0 0 0
0 1 0 0 0
0 0 1 0 0
0 0 0 -0.2425356 0.9701425
0 0 0 -0.9701425 -0.2425356
i = 4
c = -0.436436 s = 0.899735
G =
1 0 0 0 0
0 1 0 0 0
0 0 -0.4364358 0.8997354 0
0 0 -0.8997354 -0.4364358 0
0 0 0 0 1
336 9. VALORES PROPIOS
A =
2 3 4
5 4 3
-4.5825757 -5.2372294 -4.1461399
0 -1.4289915 -2.0111733
0 -0.7276069 -1.940285
Q =
1 0 0 0 0
0 1 0 0 0
0 0 -0.4364358 0.8997354 0
0 0 0.2182179 0.1058512 0.9701425
0 0 0.8728716 0.4234049 -0.2425356
...
k = 3
...
i = 4
c = -0.612372 s = 0.790569
A =
-7.0710678 -7.0710678 -5.939697
0 -2.236068 -3.5777088
0 0 -1.3856406
0 0 0
0 0 0
Q =
-0.2828427 -0.4472136 -0.5196152 0.6708204 0
-0.7071068 0.4472136 -0.2886751 -0.2236068 0.4082483
-0.2828427 0.4472136 0.0577350 0.2236068 -0.8164966
0.1414214 0.4472136 0.4041452 0.6708204 0.4082483
0.5656854 0.4472136 -0.6928203 0 0
x = B(i, : )
B(i, : ) = c B(i, : ) − s B(j, : )
B(j, : ) = s x + c B(j, : )
x = B( : , i)
B( : , i) = c B( : , i) + s B( : , j)
B( : , j) = −s x + c B( : , j)
||Qx|| = ||x||.
Sea A ∈ Rm×n , c = QT b,
U d
R= , c= ,
0qn r
Q =
-0.2828427 -0.4472136 -0.5196152 0.6708204 0
-0.7071068 0.4472136 -0.2886751 -0.2236068 0.4082483
-0.2828427 0.4472136 0.0577350 0.2236068 -0.8164966
0.1414214 0.4472136 0.4041452 0.6708204 0.4082483
0.5656854 0.4472136 -0.6928203 0 0
R =
-7.0710678 -7.0710678 -5.939697
0 -2.236068 -3.5777088
0 0 -1.3856406
0 0 0
0 0 0
U =
-7.0710678 -7.0710678 -5.939697
0 -2.236068 -3.5777088
0 0 -1.3856406
r : 0.0223607 -0.0816497
El método más popular para obtener los valores propios de una matriz
simétrica (todos reales) es el método QR. Es posiblemente el más eficente
para casos generales. El proceso tiene dos pasos:
Ejemplo 9.13.
A =
2 3 4 5
3 -1 0 1
4 0 -2 8
5 1 8 10
H =
1 0 0 0
0 0.4242641 0.5656854 0.7071068
0 0.5656854 0.4441896 -0.6947630
0 0.7071068 -0.6947630 0.1315463
H A =
2 3 4 5
7.0710678 0.2828427 4.5254834 12.020815
0 -1.2604484 -6.4464829 -2.8284271
0 -0.5755605 2.4418963 -3.5355339
H A H =
2 7.0710678 0 0
7.0710678 11.18 -6.1814444 -1.3628445
0 -6.1814444 -1.6113918 3.2154369
0 -1.3628445 3.2154369 -2.5686082
9.5. MÉTODO QR PARA VALORES PROPIOS DE MATRICES SIMÉTRICAS341
A = triHouse(A)
n = dim (A)
para k = 1 : n − 2
x = A(k + 1 : n, k)
H = H(x)
b = H(n,
H b H)
b
A = HAH b
fin-para
fin triHouse
Ejemplo 9.14.
A =
2. 3. 4. 5.
3. - 1. 0. 1.
4. 0. - 2. 8.
5. 1. 8. 10.
k = 1
H =
0.4242641 0.5656854 0.7071068
0.5656854 0.4441896 - 0.6947630
0.7071068 - 0.6947630 0.1315463
A =
2. 7.0710678 0. 0.
7.0710678 11.18 - 6.1814444 - 1.3628445
0. - 6.1814444 - 1.6113918 3.2154369
0. - 1.3628445 3.2154369 - 2.5686082
k = 2
H =
- 0.9765473 - 0.2153028
- 0.2153028 0.9765473
342 9. VALORES PROPIOS
A =
2. 7.0710678 0. 0.
7.0710678 11.18 6.3298973 0.
0. 6.3298973 - 0.3036510 - 2.7160739
0. 0. - 2.7160739 - 3.876349
A = triHouse(A)
n = dim (A)
para k = 1 : n − 2
x = A(k + 1 : n, k)
[v, β] = vHouse(x)
p = β A(k + 1 : n, k + 1 : n) v
w = p − (β/2) (pT v) v
ak+1,k = ak,k+1 = ||x||
A(k + 2 : n, k) = 0
A(k, k + 2 : n) = 0
A(k + 1 : n, k + 1 : n) = A(k + 1 : n, k + 1 : n) − v wT − w v T
fin-para
fin triHouse
A = triGivens(A)
n = dim (A)
para k = 1 : n − 2
para i = n : −1 : k + 2
[c, s] = csGivens(ai−1,k , aik )
G = G(i − 1, i, c, s, n)
A = GT AG
fin-para
fin-para
fin triHouse
Ejemplo 9.15.
A =
2 3 4 5
3 -1 0 1
4 0 -2 8
5 1 8 10
k = 1
i = 4
c = -0.624695 s = 0.780869
A =
2 3 -6.4031242 0
3 -1 -0.7808688 -0.6246950
-6.4031242 -0.7808688 13.121951 4.097561
0 -0.6246950 4.097561 -5.1219512
i = 3
c = 0.424264 s = 0.905539
A =
2 7.0710678 0 0
7.0710678 11.18 -4.9257204 -3.9755349
0 -4.9257204 0.9419512 1.1727625
0 -3.9755349 1.1727625 -5.1219512
k = 2
i = 4
344 9. VALORES PROPIOS
c = 0.778168 s = -0.628057
A =
2 7.0710678 0 0
7.0710678 11.18 -6.3298973 0
0 -6.3298973 -0.3036510 -2.7160739
0 0 -2. -3.876349
T siempre se puede expresar como una matriz diagonal por bloques, donde
cada bloque es de tamaño es 1 × 1 o de mayor tamaño pero tridiagonal no
reducido.
En el primer caso del ejemplo anterior hay un solo bloque, en el segundo
hay dos. En el tercer caso hay tres bloques, uno de ellos es 1 × 1.
Para encontrar los valores propios de T basta con encontrar los de cada
bloque tridiagonal simétrico no reducido, agregando los bloques 1 × 1 que
son valores propios de T .
9.5. MÉTODO QR PARA VALORES PROPIOS DE MATRICES SIMÉTRICAS345
QT QR = QT T
R = QT T
T + = RQ = QT T Q.
Ejemplo 9.17.
T =
2 3 0 0
3 4 5 0
0 5 6 7
0 0 7 8
R =
-3.6055513 -4.9923018 -4.1602515 0
0 -5.0076864 -5.8371805 -6.9892556
0 0 -7.6563459 -7.4712474
0 0 0 2.8863072
Q =
-0.5547002 -0.0460830 0.3365427 0.7595545
-0.8320503 0.0307220 -0.2243618 -0.5063697
0 -0.9984651 -0.0224362 -0.0506370
0 0 -0.9142743 0.4050957
T+ =
6.1538462 4.1666469 0 0
4.1666469 5.6743747 7.644594 0
0 7.644594 7.0025484 -2.6388764
0 0 -2.6388764 1.1692308
repetir
QR = T factorización QR de T
T = RQ
fin-repetir
T =
2 3 0 0
3 4 5 0
0 5 6 7
0 0 7 8
k = 1
T+ =
9.8718663 -4.486006 0 0
-4.486006 10.134151 -4.5625729 0
0 -4.5625729 -1.1770851 -0.7764250
0 0 -0.7764250 1.1710681
k = 2
T+ =
13.296028 -3.5861468 0 0
-3.5861468 8.2428763 1.7266634 0
0 1.7266634 -2.7961816 0.3062809
0 0 0.3062809 1.2572771
k = 10
T+ =
15.191934 -0.0059687 0 0
-0.0059687 6.6303783 0.0035727 0
0 0.0035727 -3.100073 0.0002528
0 0 0.0002528 1.2777606
9.5. MÉTODO QR PARA VALORES PROPIOS DE MATRICES SIMÉTRICAS347
k = 20
T+ =
15.191938 -0.0000015 0 0
-0.0000015 6.6303755 0.0000018 0
0 0.0000018 -3.1000743 3.577E-08
0 0 3.577E-08 1.2777606
T+ =
15.191938 -4.514E-09 0 0
-4.514E-09 6.6303755 -8.713E-09 0
0 -8.713E-09 -3.1000743 -7.230E-11
0 0 -7.229E-11 1.2777606
T =
2 3 0 0
3 4 5 0
0 5 6 7
0 0 7 8
T - s I =
1 3 0 0
3 3 5 0
0 5 5 7
0 0 7 7
348 9. VALORES PROPIOS
k = 9, matriz reducida:
T+ =
14.191935 -0.0052796 0 0
-0.0052796 5.5882374 0.6389663 0
0 0.6389663 -4.057933 -8.844E-12
0 0 -8.844E-12 0.2777606
T + s I
15.191935 -0.0052796 0 0
-0.0052796 6.5882374 0.6389663 0
0 0.6389663 -3.057933 -8.844E-12
0 0 -8.844E-12 1.2777606
0 0 −6 0
k = 1
mu = 2.8102497
T -mu I
5.1897503 3 0 0
3 3.1897503 -4 0
0 -4 -12.81025 -6
0 0 -6 -2.8102497
T+ = RQ
7.2885019 2.0988427 0 0
2.0988427 -9.5701241 8.9042431 0
0 8.9042431 -4.1976395 -0.6390185
0 0 -0.6390185 -0.7617370
T + mu I
10.098752 2.0988427 0 0
2.0988427 -6.7598744 8.9042431 0
0 8.9042431 -1.3873898 -0.6390185
0 0 -0.6390185 2.0485127
k = 2
mu = 2.1635102
350 9. VALORES PROPIOS
T -mu I
7.9352413 2.0988427 0 0
2.0988427 -8.9233846 8.9042431 0
0 8.9042431 -3.5509 -0.6390185
0 0 -0.6390185 -0.1149975
T+ = RQ
7.8706324 -3.26714 0 0
-3.26714 -14.885642 -2.4468061 0
0 -2.4468061 2.5541744 0.0357613
0 0 0.0357613 -0.1932052
T + mu I
10.034143 -3.26714 0 0
-3.26714 -12.722132 -2.4468061 0
0 -2.4468061 4.7176845 0.0357613
0 0 0.0357613 1.970305
k = 3
mu = 1.9698396
T -mu I
8.064303 -3.26714 0 0
-3.26714 -14.691972 -2.4468061 0
0 -2.4468061 2.7478449 0.0357613
0 0 0.0357613 0.0004654
T+ = RQ
7.1298463 5.6488809 0 0
5.6488809 -14.048752 0.5009906 0
0 0.5009906 3.0394919 0.0000006
0 0 0.0000006 0.0000557
T + mu I
9.0996859 5.6488809 0 0
5.6488809 -12.078913 0.5009906 0
0 0.5009906 5.0093315 0.0000006
0 0 0.0000006 1.9698954
9.5. MÉTODO QR PARA VALORES PROPIOS DE MATRICES SIMÉTRICAS351
k = 4
mu = 1.9698954
T -mu I
7.1297905 5.6488809 0 0
5.6488809 -14.048808 0.5009906 0
0 0.5009906 3.0394362 0.0000006
0 0 0.0000006 1.379E-13
T+ = RQ
4.4614948 -9.0220625 0 0
-9.0220625 -11.390431 -0.1052167 0
0 -0.1052167 3.049355 -2.585E-17
0 0 1.656E-22 7.811E-16
T + mu I
6.4313901 -9.0220625 0 0
-9.0220625 -9.4205358 -0.1052167 0
0 -0.1052167 5.0192503 -2.585E-17
0 0 1.656E-22 1.9698954
T reducida
6.4313901 -9.0220625 0 0
-9.0220625 -9.4205358 -0.1052167 0
0 0.1052167 5.0192503 0
0 0 0 1.9698954
−2 8 0 0 0 0
8 −2 0 0 0 0
0 0 8 3 0 0
A=
0 0 3 6 −4 0
0 0 0 −4 −10 −6
0 0 0 0 −6 0
i1 i2 : 3 6
T inicial
8 3 0 0
3 6 -4 0
0 -4 -10 -6
0 0 -6 0
mu = 2.810250
T final
10.098752 2.0988427 0 0
2.0988427 -6.7598744 8.9042431 0
0 8.9042431 -1.3873898 -0.6390185
0 0 -0.6390185 2.0485127
i1 i2 : 3 6
T inicial
10.098752 2.0988427 0 0
2.0988427 -6.7598744 8.9042431 0
0 8.9042431 -1.3873898 -0.6390185
0 0 -0.6390185 2.0485127
mu = 2.163510
T final
10.034143 -3.26714 0 0
-3.26714 -12.722132 -2.4468061 0
0 -2.4468061 4.7176845 0.0357613
0 0 0.0357613 1.970305
i1 i2 : 3 6
9.5. MÉTODO QR PARA VALORES PROPIOS DE MATRICES SIMÉTRICAS353
T inicial
10.034143 -3.26714 0 0
-3.26714 -12.722132 -2.4468061 0
0 -2.4468061 4.7176845 0.0357613
0 0 0.0357613 1.970305
mu = 1.969840
T final
9.0996859 5.6488809 0 0
5.6488809 -12.078913 0.5009906 0
0 0.5009906 5.0093315 0.0000006
0 0 0.0000006 1.9698954
i1 i2 : 3 6
T inicial
9.0996859 5.6488809 0 0
5.6488809 -12.078913 0.5009906 0
0 0.5009906 5.0093315 0.0000006
0 0 0.0000006 1.9698954
mu = 1.969895
T final
6.4313901 -9.0220625 0 0
-9.0220625 -9.4205358 -0.1052167 0
0 -0.1052167 5.0192503 8.383E-17
0 0 -1.058E-22 1.9698954
A =
-2 8 0 0 0 0
8 -2 0 0 0 0
0 0 6.4313901 -9.0220625 0 0
0 0 -9.0220625 -9.4205358 -0.1052167 0
0 0 0 -0.1052167 5.0192503 0
0 0 0 0 0 1.9698954
i1 i2 : 3 5
T inicial
6.4313901 -9.0220625 0
354 9. VALORES PROPIOS
T final
-6.2865541 11.012094 0
11.012094 3.2972548 -0.0000058
0 -0.0000058 5.0194039
i1 i2 : 3 5
T inicial
-6.2865541 11.012094 0
11.012094 3.2972548 -0.0000058
0 -0.0000058 5.0194039
mu = 5.019404
T final
-12.629095 -4.5002992 0
-4.5002992 9.6397959 2.575E-17
0 2.079E-17 5.0194039
A =
-2 8 0 0 0 0
8 -2 0 0 0 0
0 0 -12.629095 -4.5002992 0 0
0 0 -4.5002992 9.6397959 0 0
0 0 0 0 5.0194039 0
0 0 0 0 0 1.9698954
i1 i2 : 3 4
T inicial
-12.629095 -4.5002992
-4.5002992 9.6397959
mu = 10.514870
T final
-13.50417 -2.914E-16
3.384E-16 10.51487
9.5. MÉTODO QR PARA VALORES PROPIOS DE MATRICES SIMÉTRICAS355
A =
-2 8 0 0 0 0
8 -2 0 0 0 0
0 0 -13.50417 0 0 0
0 0 0 10.51487 0 0
0 0 0 0 5.0194039 0
0 0 0 0 0 1.9698954
i1 i2 : 1 2
T inicial
-2 8
8 -2
mu = -10.000000
T final
6 -8.782E-17
-1.735E-18 -10
A =
6 0 0 0 0 0
0 -10 0 0 0 0
0 0 -13.50417 0 0 0
0 0 0 10.51487 0 0
0 0 0 0 5.0194039 0
0 0 0 0 0 1.9698954
Bibliografı́a
357
358 ÍNDICE ANALÍTICO
LU, 39 matriz
PA=LU, 46 de diagonal estrictamente domi-
fórmula nante por filas, 79
de Adams-Bashforth, 266 de Givens, 67
de Simpson, 216 de Householder, 67
deAdams-Moulton, 270 definida positiva, 49, 51, 79
del trapecio, 211 jacobiana, 147
formulas positivamente definida, see ma-
de Newton-Cotes, 211, 221 triz definida positiva
fórmulas método
de Newton-Cotes, 216 de Cholesky, 49, 59
abiertas, 222 de colocación, 170
cerradas, 221 de Euler, 241, 251
funciones de la base, 172, 205 de Gauss, 30
de Gauss con pivoteo parcial, 41,
Gauss, see método de Gauss 46
Gauss-Seidel, see método de Gauss- de Gauss-Seidel, 76
Seidel de Heun, 244, 251
de la bisección, 133
Heun de la secante, 130
método de, 244, 251 de Newton, 124, 145
orden del método de, 263 de Newton en Rn , 146, 147
integración numérica, 209 de punto fijo, 138, 145
interpl, 171 de Regula Falsi, 135
de Regula Falsi modificado, 137
interpolacón polinomial, 174
de Runge-Kutta (RK), 250
interpolación, 169–171
de Runge-Kutta-Fehlberg, 258, 260
de Lagrange, 175
del disparo (shooting), 281, 282
linea, 171
del punto medio, 247, 251
por diferencias divididas, 187
del trapecio, 244
por diferencias finitas, 194
multipaso abierto, 266
interpolación polinomial por trozos,
multipaso cerrado, 270
197
multipaso explı́cito, 266
intg, 210
multipaso implı́cito, 270
Lagrange, see interpolación de La- orden del, 263
grange predictor-corrector, 270
RK, 250
Matlab, 24 RK2, 255
matrices deducción del, 255
ortogonales, 67 RK4, 251
ÍNDICE ANALÍTICO 359